"""Geometry extruded from polygonal shapes."""
from __future__ import annotations
import math
from copy import copy
from functools import lru_cache
from typing import TYPE_CHECKING, Any, ClassVar
import autograd.numpy as np
import shapely
from autograd.tracer import getval
from numpy.polynomial.legendre import leggauss as _leggauss
from pydantic import Field, field_validator, model_validator
from tidy3d.components.autograd import TracedArrayFloat2D, get_static
from tidy3d.components.autograd.path_utils import traced_paths
from tidy3d.components.autograd.types import PathType, TracedFloat
from tidy3d.components.autograd.utils import hasbox
from tidy3d.components.base import cached_property
from tidy3d.components.transformation import ReflectionFromPlane, RotationAroundAxis
from tidy3d.components.types import TYPE_TAG_STR, ArrayFloat1D
from tidy3d.config import config
from tidy3d.constants import LARGE_NUMBER, MICROMETER, fp_eps
from tidy3d.exceptions import AdjointError, SetupError, Tidy3dImportError, ValidationError
from tidy3d.log import log
from tidy3d.packaging import verify_packages_import
from . import base, triangulation
from .vertex_utils import remove_adjacent_duplicate_vertices
if TYPE_CHECKING:
from gdstk import Cell
from numpy.typing import NDArray
from pydantic import PositiveFloat
from tidy3d.compat import Self
from tidy3d.components.autograd import AutogradFieldMap
from tidy3d.components.autograd.derivative_utils import DerivativeInfo
from tidy3d.components.autograd.path_utils import AutogradRoute
from tidy3d.components.types import (
ArrayFloat2D,
ArrayLike,
Axis,
Bound,
Coordinate,
MatrixReal4x4,
PlanePosition,
Shapely,
)
# sampling polygon along dilation for validating polygon to be
# non self-intersecting during the entire dilation process
_N_SAMPLE_POLYGON_INTERSECT = 5
_IS_CLOSE_RTOL = np.finfo(float).eps
# Warn for too many divided polyslabs
_COMPLEX_POLYSLAB_DIVISIONS_WARN = 100
# Warn before triangulating large polyslabs due to inefficiency
_MAX_POLYSLAB_VERTICES_FOR_TRIANGULATION = 500
_MIN_POLYGON_AREA = fp_eps
# Angular resolution for arc discretization (approximately 1 degrees)
_ARC_DISCRETIZATION_ANGULAR_STEP = np.pi / 180
# Minimum number of intermediate points per arc segment for discretization
_MIN_POINTS_PER_ARC = 5
@lru_cache(maxsize=128)
def leggauss(n: int) -> tuple[NDArray, NDArray]:
"""Cached version of leggauss with dtype conversions."""
g, w = _leggauss(n)
return g.astype(config.adjoint.gradient_dtype_float, copy=False), w.astype(
config.adjoint.gradient_dtype_float, copy=False
)
class _PolyBulgeUtil:
"""Pure-geometry utility methods for arc/bulge computations on polygons for internal use.
All methods in this class are static methods.
"""
@staticmethod
def _arcs_from_bulges(
arc_starts: ArrayFloat2D,
arc_ends: ArrayFloat2D,
bulges: ArrayFloat1D,
) -> dict:
"""Compute arc geometry for multiple arcs at once (vectorized).
Parameters
----------
arc_starts : ArrayFloat2D
Shape (N, 2) start points of each arc.
arc_ends : ArrayFloat2D
Shape (N, 2) end points of each arc.
bulges : ArrayFloat1D
Shape (N,) bulge values. Must all be non-zero; an error is raised if any are ≈ 0.
Returns
-------
dict
Dictionary with arrays of shape (N,) or (N, 2):
'centers', 'radii', 'start_angles', 'end_angles',
'included_angles'.
"""
arc_starts = np.asarray(arc_starts, dtype=float)
arc_ends = np.asarray(arc_ends, dtype=float)
bulges = np.asarray(bulges, dtype=float)
# Validate that no bulges are close to 0
if np.any(np.isclose(bulges, 0)):
raise ValueError(
"All bulges must be non-zero. Edges with bulge ≈ 0 should be filtered out "
"before calling this method, as they represent straight line segments."
)
chord_vecs = arc_ends - arc_starts # (N, 2)
chord_lengths = np.linalg.norm(chord_vecs, axis=1) # (N,)
abs_bulges = np.abs(bulges)
included_angles = 4.0 * np.arctan(bulges) # (N,)
# Sagitta: s = (chord/2) * |bulge|
sagittas = (chord_lengths / 2.0) * abs_bulges # (N,)
# Radius: r = ((chord/2)^2 + s^2) / (2*s)
half_chords = chord_lengths / 2.0
radii = (half_chords**2 + sagittas**2) / (2.0 * sagittas) # (N,)
# Midpoints of each chord
chord_midpoints = (arc_starts + arc_ends) / 2.0 # (N, 2)
# Perpendicular unit vectors (rotate chord by 90 degrees CCW)
chord_perpendiculars = (
np.column_stack([-chord_vecs[:, 1], chord_vecs[:, 0]]) / chord_lengths[:, np.newaxis]
) # (N, 2)
# Distance from midpoint to center
midpoint_to_center_dist = radii - sagittas # (N,)
# Center: positive bulge -> center at mid + d*perp, negative -> mid - d*perp
bulge_signs = np.sign(bulges) # (N,)
centers = (
chord_midpoints
+ (bulge_signs * midpoint_to_center_dist)[:, np.newaxis] * chord_perpendiculars
) # (N, 2)
# Start and end angles
start_angles = np.arctan2(
arc_starts[:, 1] - centers[:, 1], arc_starts[:, 0] - centers[:, 0]
) # (N,)
end_angles = np.arctan2(
arc_ends[:, 1] - centers[:, 1], arc_ends[:, 0] - centers[:, 0]
) # (N,)
return {
"centers": centers,
"radii": radii,
"start_angles": start_angles,
"end_angles": end_angles,
"included_angles": included_angles,
}
@staticmethod
def _compute_bulge_data(vertices: ArrayFloat2D, bulges: ArrayFloat1D) -> dict:
"""Precompute arc geometry from vertices and bulges, e.g. encoding which edges are arcs.
Parameters
----------
vertices : ArrayFloat2D
Shape (N, 2) polygon vertices.
bulges : ArrayFloat1D
Shape (N,) bulge values for each edge.
Returns
-------
dict
Always contains 'arc_mask' that encodes which edges are arcs. When arcs exist, also contains
'centers', 'radii', 'start_angles', 'end_angles', 'included_angles'.
"""
bulges = np.asarray(bulges)
arc_mask = ~np.isclose(bulges, 0)
if not np.any(arc_mask):
return {"arc_mask": arc_mask}
edge_starts = np.asarray(vertices)
edge_ends = np.roll(edge_starts, -1, axis=0)
arc_data = _PolyBulgeUtil._arcs_from_bulges(
edge_starts[arc_mask], edge_ends[arc_mask], bulges[arc_mask]
)
return {"arc_mask": arc_mask, **arc_data}
@staticmethod
def _polygon_arcs_bounds(
vertices: ArrayFloat2D,
bulge_data: dict,
) -> tuple[ArrayFloat2D, ArrayFloat2D]:
"""Compute analytical bounding boxes for arc segments in a polygon.
An arc's bounding box is determined by its endpoints plus any cardinal
extrema (rightmost, topmost, leftmost, bottommost points of the underlying
circle) that the arc actually passes through. The method:
1. Initializes bounds from the two endpoints of each arc.
2. Checks which of the 4 cardinal angles (0, pi/2, pi, 3pi/2) fall within
each arc's angular sweep, accounting for sweep direction (CCW vs CW)
and wraparound at 0/2pi.
3. For each cardinal angle the arc passes through, expands the bounding box
to include the corresponding extremal point (center +/- radius).
Parameters
----------
vertices : ArrayFloat2D
Shape (N, 2) polygon vertices.
bulge_data : dict
Precomputed arc data from ``_compute_bulge_data``.
Returns
-------
tuple[ArrayFloat2D, ArrayFloat2D]
(min_coords, max_coords) each of shape (N, 2) where N is the number of arc segments.
Raises
------
KeyError
If ``bulge_data`` contains no arcs (all bulges are zero).
"""
arc_mask = bulge_data["arc_mask"]
edge_starts = np.asarray(vertices)
edge_ends = np.roll(edge_starts, -1, axis=0)
arc_starts = edge_starts[arc_mask]
arc_ends = edge_ends[arc_mask]
# Initialize bounds from endpoints
min_coords = np.minimum(arc_starts, arc_ends) # (N, 2)
max_coords = np.maximum(arc_starts, arc_ends) # (N, 2)
# The 4 cardinal angles correspond to the extremal points of the circle:
# pi/2 (top: max y)
# |
# pi (left: min x) ---+--- 0 (right: max x)
# |
# 3pi/2 (bottom: min y)
cardinal_angles = np.array([0.0, np.pi / 2, np.pi, 3 * np.pi / 2]) # (4,)
# Normalize from arctan2's [-pi, pi] to [0, 2*pi) so the wraparound
# logic below can use a single consistent range.
start_angles_norm = bulge_data["start_angles"] % (2 * np.pi) # (N,)
end_angles_norm = bulge_data["end_angles"] % (2 * np.pi) # (N,)
# For each arc, check which cardinal angles lie within its angular sweep.
# Broadcast: (N, 1) vs (1, 4) -> (N, 4) boolean matrix
start_angles_norm_col = start_angles_norm[:, np.newaxis] # (N, 1)
end_angles_norm_col = end_angles_norm[:, np.newaxis] # (N, 1)
cardinal_row = cardinal_angles[np.newaxis, :] # (1, 4)
sweeps_ccw_col = (bulge_data["included_angles"] > 0)[:, np.newaxis] # (N, 1)
# CCW: angles increase from start to end.
# No wrap (end > start): cardinal in [start, end]
# Wrap (end <= start): cardinal in [start, 2*pi) or [0, end]
ccw_simple = (cardinal_row >= start_angles_norm_col) & (cardinal_row <= end_angles_norm_col)
ccw_wrapped = (end_angles_norm_col <= start_angles_norm_col) & (
(cardinal_row >= start_angles_norm_col) | (cardinal_row <= end_angles_norm_col)
)
ccw_hits_cardinal = ccw_simple | ccw_wrapped
# CW: angles decrease from start to end.
# No wrap (end < start): cardinal in [end, start]
# Wrap (end >= start): cardinal in [0, start] or [end, 2*pi)
cw_simple = (cardinal_row <= start_angles_norm_col) & (cardinal_row >= end_angles_norm_col)
cw_wrapped = (end_angles_norm_col >= start_angles_norm_col) & (
(cardinal_row <= start_angles_norm_col) | (cardinal_row >= end_angles_norm_col)
)
cw_hits_cardinal = cw_simple | cw_wrapped
arc_hits_cardinal = np.where(sweeps_ccw_col, ccw_hits_cardinal, cw_hits_cardinal) # (N, 4)
# Expand bounds for each cardinal the arc passes through.
# Cardinal directions: 0=right(0°), 1=top(90°), 2=left(180°), 3=bottom(270°)
for n, _ in enumerate(cardinal_angles):
coord_idx = n % 2 # Alternates x(0), y(1), x(0), y(1)
sign = (-1) ** (n // 2) # Pattern: +1, +1, -1, -1
bounds = max_coords if n < 2 else min_coords # max for 0,1; min for 2,3
op_func = np.maximum if n < 2 else np.minimum
hit = arc_hits_cardinal[:, n]
extremal_values = bulge_data["centers"][:, coord_idx] + sign * bulge_data["radii"]
bounds[:, coord_idx] = np.where(
hit, op_func(bounds[:, coord_idx], extremal_values), bounds[:, coord_idx]
)
return min_coords, max_coords
@staticmethod
def _polygon_discretize(
vertices: ArrayFloat2D, bulge_data: dict, num_points_per_arc: int | None = None
) -> ArrayFloat2D:
"""Discretize a polygon with arc segments into a pure polyline.
This method converts arc edges (defined by bulge factors) into straight line segments,
producing a polyline approximation. It is used for:
- Validation: checking for self-intersection
- Geometric queries: inside point testing, plane intersections
- Visualization: rendering polygons with arcs as polylines
Parameters
----------
vertices : ArrayFloat2D
Shape (N, 2) defining polygon vertices.
bulge_data : dict
Precomputed arc data from ``_compute_bulge_data``.
num_points_per_arc : int, optional
Number of intermediate points per arc segment. If None, the number of points
is automatically determined based on the arc's included angle, using approximately
one point per degree and a minimum of 5 points per arc.
Returns
-------
ArrayFloat2D
Discretized vertices with arcs converted to polylines.
"""
vertices = np.asarray(vertices)
arc_mask = bulge_data["arc_mask"]
num_vertices = len(vertices)
if not np.any(arc_mask):
return vertices.copy()
arc_indices = np.where(arc_mask)[0]
# Compute number of points per arc (vectorized)
if num_points_per_arc is not None:
points_per_arc = np.full(len(arc_indices), num_points_per_arc, dtype=int)
else:
abs_angles = np.abs(bulge_data["included_angles"])
points_per_arc = np.maximum(
_MIN_POINTS_PER_ARC, np.ceil(abs_angles / _ARC_DISCRETIZATION_ANGULAR_STEP)
).astype(int)
# Build a lookup: edge index -> position in arc arrays
edge_to_arc_idx = dict(zip(arc_indices, range(len(arc_indices))))
# Loop to assemble variable-length per-arc points
all_points = []
for vertex_idx in range(num_vertices):
all_points.append(vertices[vertex_idx])
if vertex_idx in edge_to_arc_idx:
arc_pos = edge_to_arc_idx[vertex_idx]
center = bulge_data["centers"][arc_pos]
radius = bulge_data["radii"][arc_pos]
start_angle = bulge_data["start_angles"][arc_pos]
end_angle = bulge_data["end_angles"][arc_pos]
# Generate angles
if bulge_data["included_angles"][arc_pos] > 0:
if end_angle < start_angle:
end_angle += 2 * np.pi
else:
if end_angle > start_angle:
end_angle -= 2 * np.pi
angles = np.linspace(start_angle, end_angle, points_per_arc[arc_pos] + 2)
# Skip first (edge_start) and last (edge_end) — last will be added as next vertex
angles = angles[1:-1]
if len(angles) > 0:
arc_points = np.column_stack(
[center[0] + radius * np.cos(angles), center[1] + radius * np.sin(angles)]
)
all_points.extend(arc_points)
return np.array(all_points)
@staticmethod
def _arc_segment_area(bulge_data: dict) -> float:
"""Compute the signed total area of circular segments from arc edges.
Each arc edge contributes a circular segment (the region between the
arc and its chord). This correction is meant to be added to the
shoelace (chord-based) polygon area.
Parameters
----------
bulge_data : dict
Precomputed arc data from ``_compute_bulge_data``.
Returns
-------
float
Signed sum of circular segment areas.
Raises
------
KeyError
If ``bulge_data`` contains no arcs (all bulges are zero).
"""
radii = bulge_data["radii"]
theta = np.abs(bulge_data["included_angles"])
# Circular segment area: A = r² * (θ - sin(θ)) / 2
segment_areas = (radii**2) * (theta - np.sin(theta)) / 2.0
# Sign: positive included_angle → outward bulge → adds area (for CCW polygon)
bulge_signs = np.sign(bulge_data["included_angles"])
return float(np.sum(bulge_signs * segment_areas))
@staticmethod
def _polygon_perimeter(vertices: ArrayFloat2D, bulge_data: dict) -> float:
"""Compute the polygon perimeter using arc lengths for curved segments.
Parameters
----------
vertices : ArrayFloat2D
Shape (N, 2) defining the polygon vertices in the xy-plane.
bulge_data : dict
Precomputed arc data from ``_compute_bulge_data``.
Returns
-------
float
Total perimeter including arc lengths.
"""
next_vertices = np.roll(vertices, shift=-1, axis=0)
edge_lengths = np.linalg.norm(next_vertices - vertices, axis=1) # (N,)
arc_mask = bulge_data["arc_mask"]
if not np.any(arc_mask):
return float(np.sum(edge_lengths))
# For arc edges: replace chord lengths with arc lengths
radii = bulge_data["radii"]
edge_lengths[arc_mask] = radii * np.abs(bulge_data["included_angles"])
return float(np.sum(edge_lengths))
[docs]
class PolySlab(base.Planar):
"""Polygon extruded with optional sidewall angle along axis direction.
Polygon edges may optionally be circular arcs via the ``bulges`` parameter, which
follows the convention: ``bulge = tan(θ/4)`` where ``θ`` is the
included angle of the arc. For CCW polygon, a positive bulge produces an arc that bulges outward
; a negative bulge bulges inward; and zero means a straight edge.
Example
-------
>>> vertices = np.array([(0,0), (1,0), (1,1)])
>>> p = PolySlab(vertices=vertices, axis=2, slab_bounds=(-1, 1))
>>> # With arc segments (a rounded edge from vertex 1 to vertex 2):
>>> p_arc = PolySlab(
... vertices=vertices, axis=2, slab_bounds=(-1, 1), bulges=[0, 0.5, 0]
... )
"""
_traced_supported_paths: ClassVar[tuple[PathType, ...]] = traced_paths(
"vertices",
"sidewall_angle",
("slab_bounds", 0),
("slab_bounds", 1),
)
slab_bounds: tuple[TracedFloat, TracedFloat] = Field(
title="Slab Bounds",
description="Minimum and maximum positions of the slab along axis dimension.",
json_schema_extra={"units": MICROMETER},
)
dilation: float = Field(
0.0,
title="Dilation",
description="Dilation of the supplied polygon by shifting each edge along its "
"normal outwards direction by a distance; a negative value corresponds to erosion.",
json_schema_extra={"units": MICROMETER},
)
vertices: TracedArrayFloat2D = Field(
title="Vertices",
description="List of (d1, d2) defining the 2 dimensional positions of the polygon "
"face vertices at the ``reference_plane``. "
"The index of dimension should be in the ascending order: e.g. if "
"the slab normal axis is ``axis=y``, the coordinate of the vertices will be in (x, z)",
json_schema_extra={"units": MICROMETER},
)
bulges: ArrayFloat1D | None = Field(
None,
title="Bulges In Segments",
description="List of values describing the bulge of each edge, following the typical convention: "
"``bulge = tan(θ/4)`` where θ is the included angle of the arc. "
"When ``bulges=None``, all segments are straight lines; otherwise, the number of bulge values "
"must equal the number of vertices. ``bulges[i]`` defines the arc on the edge from vertex ``i`` "
"to vertex ``(i+1) mod N``. "
"Sign convention: ``bulge=0`` means straight line; ``bulge>0`` bulges outward "
"(for CCW polygons); ``bulge<0`` bulges inward. "
"All values must be finite. "
"Canonicalization: bulges are kept synchronized with vertices during internal canonicalization. "
"When winding order is reversed to CCW (i.e., when the input polygon is CW), bulges are permuted to match the reversed edges and "
"their signs are flipped. When duplicate adjacent vertices are removed, the corresponding "
"zero-bulge entry is dropped; a zero-length edge with non-zero bulge raises ``SetupError``.",
)
[docs]
@staticmethod
def make_shapely_polygon(vertices: ArrayLike) -> shapely.Polygon:
"""Make a shapely polygon from some vertices, first ensures they are untraced."""
vertices = get_static(vertices)
return shapely.Polygon(vertices)
[docs]
@field_validator("slab_bounds")
@classmethod
def slab_bounds_order(cls, val: tuple[float, float]) -> tuple[float, float]:
"""Maximum position of the slab should be no smaller than its minimal position."""
if val[1] < val[0]:
raise SetupError(
"Polyslab.slab_bounds must be specified in the order of "
"minimum and maximum positions of the slab along the axis. "
f"But now the maximum {val[1]} is smaller than the minimum {val[0]}."
)
return val
[docs]
@field_validator("vertices")
@classmethod
def correct_shape(cls, val: ArrayFloat2D) -> ArrayFloat2D:
"""Makes sure vertices size is correct."""
# overall shape of vertices
if val.shape[1] != 2:
raise SetupError(
"PolySlab.vertices must be a 2 dimensional array shaped (N, 2). "
f"Given array with shape of {val.shape}."
)
return val
@model_validator(mode="after")
def _validate_bulge_array_length(self: Self) -> Self:
"""When ``bulges`` is not None, its length must match ``vertices``."""
if self.bulges is None:
return self
if len(self.bulges) != self.vertices.shape[0]:
self._raise_validation_error_at_loc(
SetupError(
"The number of bulge values in 'bulges' must match the number of 'vertices'. "
f"However, this polyslab contains {self.vertices.shape[0]} vertices, but {len(self.bulges)} "
"bulge values."
),
"bulges",
)
return self
@model_validator(mode="after")
def _validate_bulges_are_finite(self: Self) -> Self:
"""All bulge values must be finite (no inf or nan)."""
if self.bulges is None:
return self
if not np.all(np.isfinite(self.bulges)):
self._raise_validation_error_at_loc(
SetupError(
"All 'bulges' values must be finite. Infinite or NaN bulge values are not allowed."
),
"bulges",
)
return self
@model_validator(mode="after")
def _validate_no_bulge_on_zero_length_edges(self: Self) -> Self:
"""Zero-length edges cannot have non-zero bulge values."""
if self.bulges is None:
return self
# Detect zero-length edges
vertices = np.array(self.vertices)
bulges = np.array(self.bulges)
vertices_next = np.roll(vertices, shift=-1, axis=0)
edge_lengths = np.linalg.norm(vertices - vertices_next, axis=1)
is_zero_length = np.isclose(edge_lengths, 0, rtol=_IS_CLOSE_RTOL)
# Check for zero-length edges with non-zero bulge
has_nonzero_bulge = ~np.isclose(bulges, 0)
invalid_edges = is_zero_length & has_nonzero_bulge
if np.any(invalid_edges):
invalid_idx = np.where(invalid_edges)[0][0] # Report first invalid edge
self._raise_validation_error_at_loc(
SetupError(
f"Zero-length edge at vertex index {invalid_idx} has a non-zero bulge value "
f"({bulges[invalid_idx]}). A zero-length edge cannot define an arc segment."
),
"bulges",
)
return self
@cached_property
def _polygon_not_degenerate_or_self_intersecting(self) -> bool:
"""Raise if the polygon encloses no area or is self-intersecting."""
vertices = get_static(self.vertices)
if self.bulges is None or np.allclose(self.bulges, 0) or vertices.shape[0] < 3:
# make sure no polygon splitting, islands, 0 area
poly_heal = shapely.make_valid(self.make_shapely_polygon(vertices))
if poly_heal.area < _MIN_POLYGON_AREA:
self._raise_validation_error_at_loc(
SetupError("The polygon almost collapses to a 1D curve."),
"vertices",
)
if not poly_heal.geom_type == "Polygon" or len(poly_heal.interiors) > 0:
self._raise_validation_error_at_loc(
SetupError(
"Polygon is self-intersecting, resulting in "
"polygon splitting or generation of holes/islands. "
"A general treatment to self-intersecting polygon will be available "
"in future releases."
),
"vertices",
)
return True
# with arc segments, the chord polygon can be degenerate (e.g. collinear
# vertices) while the arc contour encloses a finite area; gate on the
# bulge-aware signed area instead. Self-intersection of the arc contour
# is checked in '_validate_no_self_intersection_with_arcs'.
if abs(PolySlab._signed_area_with_bulges(vertices, self.bulges)) < _MIN_POLYGON_AREA:
self._raise_validation_error_at_loc(
SetupError("The polygon almost collapses to a 1D curve."),
"vertices",
)
return True
@model_validator(mode="after")
def _validate_polygon_not_degenerate_or_self_intersecting(self: Self) -> Self:
"""Polygon must enclose a finite area and not be self-intersecting.
The check is cached on the instance so it runs once: embedding into a parent
model (e.g. 'GeometryGroup') revalidates the same instance and would otherwise
repeat the expensive shapely work for every embedding.
"""
_ = self._polygon_not_degenerate_or_self_intersecting
return self
@model_validator(mode="after")
def _vertical_sidewall_no_dilation_with_arc(self: Self) -> Self:
"""In the presence of arc segments, for now sidwall is limited to vertical, and no dilation allowed."""
if self.bulges is None:
return self
if not np.allclose(self.bulges, 0):
if not math.isclose(self.sidewall_angle, 0):
self._raise_validation_error_at_loc(
SetupError(
"Arc segments are present in the polygon. 'sidewall_angle' must be 0. "
"A general treatment for slanted polyslab containing arc segments "
"will be available in future releases."
),
"sidewall_angle",
)
if not math.isclose(self.dilation, 0):
self._raise_validation_error_at_loc(
SetupError(
"Arc segments are present in the polygon. 'dilation' value must be 0. "
"A general treatment will be available in future releases."
),
"dilation",
)
return self
@model_validator(mode="after")
def _validate_no_self_intersection_with_arcs(self: Self) -> Self:
"""Ensure polygon with arc segments is not self-intersecting."""
if self.bulges is None:
return self
if np.allclose(self.bulges, 0):
return self # No arcs, defer to existing vertex-based validation
# Use canonical (deduplicated, CCW) vertices and bulges
canon_verts, canon_bulges = PolySlab._canonicalize_vertices_and_bulges(
self.vertices, self.bulges
)
# Discretize the polygon with arc segments
bulge_data = _PolyBulgeUtil._compute_bulge_data(canon_verts, canon_bulges)
discretized = _PolyBulgeUtil._polygon_discretize(canon_verts, bulge_data)
# Check if the resulting polygon is valid
poly = shapely.Polygon(discretized)
if not poly.is_valid or not poly.is_simple:
# when the straight-edge polygon is itself self-intersecting, the
# vertex order -- not the bulge values -- is the actionable input
chord_heal = shapely.make_valid(self.make_shapely_polygon(self.vertices))
if chord_heal.area >= _MIN_POLYGON_AREA and (
chord_heal.geom_type != "Polygon" or len(chord_heal.interiors) > 0
):
self._raise_validation_error_at_loc(
SetupError(
"Polygon is self-intersecting, resulting in "
"polygon splitting or generation of holes/islands. "
"A general treatment to self-intersecting polygon will be available "
"in future releases."
),
"vertices",
)
self._raise_validation_error_at_loc(
SetupError(
"Polygon with arc segments is self-intersecting. "
"Please adjust the bulge values to avoid self-intersection."
),
"bulges",
)
return self
[docs]
@model_validator(mode="after")
def no_complex_self_intersecting_polygon_at_reference_plane(self: Self) -> Self:
"""At the reference plane, check if the polygon is self-intersecting.
There are two types of self-intersection that can occur during dilation:
1) the one that creates holes/islands, or splits polygons, or removes everything;
2) the one that does not.
For 1), we issue an error since it is yet to be supported;
For 2), we heal the polygon, and warn that the polygon has been cleaned up.
"""
val = self.vertices
# no need to validate anything here
if math.isclose(self.dilation, 0):
return self
val_np = PolySlab._proper_vertices(val)
dist = self.dilation
# 0) fully eroded
if dist < 0 and dist < -PolySlab._maximal_erosion(val_np):
self._raise_validation_error_at_loc(
SetupError("Erosion value is too large. The polygon is fully eroded."),
"dilation",
)
# no edge events
if not PolySlab._edge_events_detection(val_np, dist, ignore_at_dist=False):
return self
poly_offset = PolySlab._shift_vertices(val_np, dist)[0]
if PolySlab._area(poly_offset) < fp_eps**2:
self._raise_validation_error_at_loc(
SetupError("Erosion value is too large. The polygon is fully eroded."),
"dilation",
)
# edge events
poly_offset = shapely.make_valid(self.make_shapely_polygon(poly_offset))
# 1) polygon split or create holes/islands
if not poly_offset.geom_type == "Polygon" or len(poly_offset.interiors) > 0:
self._raise_validation_error_at_loc(
SetupError(
"Dilation/Erosion value is too large, resulting in "
"polygon splitting or generation of holes/islands. "
"A general treatment to self-intersecting polygon will be available "
"in future releases."
),
"dilation",
)
# case 2
log.warning(
"The dilation/erosion value is too large. resulting in a "
"self-intersecting polygon. "
"The vertices have been modified to make a valid polygon."
)
return self
[docs]
@model_validator(mode="after")
def no_self_intersecting_polygon_during_extrusion(self: Self) -> Self:
"""In this simple polyslab, we don't support self-intersecting polygons yet, meaning that
any normal cross section of the PolySlab cannot be self-intersecting. This part checks
if any self-interction will occur during extrusion with non-zero sidewall angle.
There are two types of self-intersection, known as edge events,
that can occur during dilation:
1) neighboring vertex-vertex crossing. This type of edge event can be treated with
``ComplexPolySlab`` which divides the polyslab into a list of simple polyslabs.
2) other types of edge events that can create holes/islands or split polygons.
To detect this, we sample _N_SAMPLE_POLYGON_INTERSECT cross sections to see if any creation
of polygons/holes, and changes in vertices number.
"""
val = self.vertices
# no need to validate anything here
# sidewall_angle may be autograd-traced; use static value for this check only
if math.isclose(getval(self.sidewall_angle), 0):
return self
# apply dilation
poly_ref = PolySlab._proper_vertices(val)
if not math.isclose(self.dilation, 0):
poly_ref = PolySlab._shift_vertices(poly_ref, self.dilation)[0]
poly_ref = PolySlab._heal_polygon(poly_ref)
slab_bounds = get_static(self.slab_bounds)
_slab_min, _slab_max = slab_bounds
# first, check vertex-vertex crossing at any point during extrusion
length = slab_bounds[1] - slab_bounds[0]
dist = [-length * np.tan(self.sidewall_angle)]
# reverse the dilation value if it's defined on the top
if self.reference_plane == "top":
dist = [-dist[0]]
# for middle, both direction needs to be examined
elif self.reference_plane == "middle":
dist = [dist[0] / 2, -dist[0] / 2]
# capture vertex crossing events
max_thick = []
for dist_val in dist:
max_dist = PolySlab._neighbor_vertices_crossing_detection(poly_ref, dist_val)
if max_dist is not None:
max_thick.append(max_dist / abs(dist_val) * length)
if len(max_thick) > 0:
max_thick = min(max_thick)
self._raise_validation_error_at_loc(
SetupError(
"Sidewall angle or structure thickness is so large that the polygon "
"is self-intersecting during extrusion. "
f"Please either reduce structure thickness to be < {max_thick:.3e}, "
"or use our plugin 'ComplexPolySlab' to divide the complex polyslab "
"into a list of simple polyslabs."
),
"sidewall_angle",
)
# vertex-edge crossing event.
for dist_val in dist:
if PolySlab._edge_events_detection(poly_ref, dist_val):
self._raise_validation_error_at_loc(
SetupError(
"Sidewall angle or structure thickness is too large, "
"resulting in polygon splitting or generation of holes/islands. "
"A general treatment to self-intersecting polygon will be available "
"in future releases."
),
"sidewall_angle",
)
return self
[docs]
@classmethod
def from_gds(
cls,
gds_cell: Cell,
axis: Axis,
slab_bounds: tuple[float, float],
gds_layer: int,
gds_dtype: int | None = None,
gds_scale: PositiveFloat = 1.0,
dilation: float = 0.0,
sidewall_angle: float = 0,
reference_plane: PlanePosition = "middle",
) -> list[PolySlab]:
"""Import :class:`PolySlab` from a ``gdstk.Cell``.
Parameters
----------
gds_cell : gdstk.Cell
``gdstk.Cell`` containing 2D geometric data.
axis : int
Integer index into the polygon's slab axis. (0,1,2) -> (x,y,z).
slab_bounds: tuple[float, float]
Minimum and maximum positions of the slab along ``axis``.
gds_layer : int
Layer index in the ``gds_cell``.
gds_dtype : int = None
Data-type index in the ``gds_cell``.
If ``None``, imports all data for this layer into the returned list.
gds_scale : float = 1.0
Length scale used in GDS file in units of MICROMETER.
For example, if gds file uses nanometers, set ``gds_scale=1e-3``.
Must be positive.
dilation : float = 0.0
Dilation of the polygon in the base by shifting each edge along its
normal outwards direction by a distance;
a negative value corresponds to erosion.
sidewall_angle : float = 0
Angle of the sidewall.
``sidewall_angle=0`` (default) specifies vertical wall,
while ``0<sidewall_angle<np.pi/2`` for the base to be larger than the top.
reference_plane : PlanePosition = "middle"
The position of the GDS layer. It can be at the ``bottom``, ``middle``,
or ``top`` of the PolySlab. E.g. if ``axis=1``, ``bottom`` refers to the
negative side of y-axis, and ``top`` refers to the positive side of y-axis.
Returns
-------
list[:class:`PolySlab`]
List of :class:`PolySlab` objects sharing ``axis`` and slab bound properties.
"""
all_vertices = PolySlab._load_gds_vertices(gds_cell, gds_layer, gds_dtype, gds_scale)
return [
cls(
vertices=verts,
axis=axis,
slab_bounds=slab_bounds,
dilation=dilation,
sidewall_angle=sidewall_angle,
reference_plane=reference_plane,
)
for verts in all_vertices
]
@staticmethod
def _load_gds_vertices(
gds_cell: Cell,
gds_layer: int,
gds_dtype: int | None = None,
gds_scale: PositiveFloat = 1.0,
) -> list[ArrayFloat2D]:
"""Import :class:`PolySlab` from a ``gdstk.Cell``.
Parameters
----------
gds_cell : gdstk.Cell
``gdstk.Cell`` containing 2D geometric data.
gds_layer : int
Layer index in the ``gds_cell``.
gds_dtype : int = None
Data-type index in the ``gds_cell``.
If ``None``, imports all data for this layer into the returned list.
gds_scale : float = 1.0
Length scale used in GDS file in units of MICROMETER.
For example, if gds file uses nanometers, set ``gds_scale=1e-3``.
Must be positive.
Returns
-------
list[ArrayFloat2D]
List of :class:`.ArrayFloat2D`
"""
import gdstk
gds_cell_class_name = str(gds_cell.__class__)
if not isinstance(gds_cell, gdstk.Cell):
if (
"gdstk" in gds_cell_class_name
): # Check if it might be a gdstk cell but gdstk is not found
raise Tidy3dImportError(
"Module 'gdstk' not found. It is required to import gdstk cells."
)
raise ValueError(
f"validate 'gds_cell' of type '{gds_cell_class_name}' "
"does not seem to be associated with 'gdstk' package "
"and therefore can't be loaded by Tidy3D."
)
all_vertices = base.Geometry.load_gds_vertices_gdstk(
gds_cell=gds_cell,
gds_layer=gds_layer,
gds_dtype=gds_dtype,
gds_scale=gds_scale,
)
# convert vertices into polyslabs
polygons = [PolySlab.make_shapely_polygon(vertices).buffer(0) for vertices in all_vertices]
polys_union = shapely.unary_union(polygons, grid_size=base.POLY_GRID_SIZE)
if polys_union.geom_type == "Polygon":
all_vertices = [np.array(polys_union.exterior.coords)]
elif polys_union.geom_type == "MultiPolygon":
all_vertices = [np.array(polygon.exterior.coords) for polygon in polys_union.geoms]
return all_vertices
@property
def center_axis(self) -> float:
"""Gets the position of the center of the geometry in the out of plane dimension."""
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
def length_axis(self) -> float:
"""Gets the length of the geometry along the out of plane dimension."""
zmin, zmax = self.slab_bounds
return zmax - zmin
@property
def finite_length_axis(self) -> float:
"""Gets the length of the PolySlab along the out of plane dimension.
First clips the slab bounds to LARGE_NUMBER and then returns difference.
"""
zmin, zmax = self.slab_bounds
zmin = max(zmin, -LARGE_NUMBER)
zmax = min(zmax, LARGE_NUMBER)
return zmax - zmin
@cached_property
def reference_polygon(self) -> NDArray:
"""The polygon at the reference plane.
Returns
-------
ArrayLike[float, float]
The vertices of the polygon at the reference plane.
"""
vertices = self._canonical_vertices_and_bulges[0]
if math.isclose(self.dilation, 0):
return vertices
offset_vertices = self._shift_vertices(vertices, self.dilation)[0]
return self._heal_polygon(offset_vertices)
@cached_property
def middle_polygon(self) -> NDArray:
"""The polygon at the middle.
Returns
-------
ArrayLike[float, float]
The vertices of the polygon at the middle.
"""
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":
return self._shift_vertices(self.reference_polygon, -dist)[0]
# middle case
return self.reference_polygon
@cached_property
def base_polygon(self) -> NDArray:
"""The polygon at the base, derived from the ``middle_polygon``.
Returns
-------
ArrayLike[float, float]
The vertices of the polygon at the base.
"""
if self.reference_plane == "bottom":
return self.reference_polygon
dist = self._extrusion_length_to_offset_distance(-self.finite_length_axis / 2)
return self._shift_vertices(self.middle_polygon, dist)[0]
@cached_property
def top_polygon(self) -> NDArray:
"""The polygon at the top, derived from the ``middle_polygon``.
Returns
-------
ArrayLike[float, float]
The vertices of the polygon at the top.
"""
if self.reference_plane == "top":
return self.reference_polygon
dist = self._extrusion_length_to_offset_distance(self.finite_length_axis / 2)
return self._shift_vertices(self.middle_polygon, dist)[0]
@cached_property
def _normal_2dmaterial(self) -> Axis:
"""Get the normal to the given geometry, checking that it is a 2D geometry."""
if self.slab_bounds[0] != self.slab_bounds[1]:
raise ValidationError("'Medium2D' requires the 'PolySlab' bounds to be equal.")
return self.axis
def _update_from_bounds(self, bounds: tuple[float, float], axis: Axis) -> PolySlab:
"""Returns an updated geometry which has been transformed to fit within ``bounds``
along the ``axis`` direction."""
if axis != self.axis:
raise ValueError(
f"'_update_from_bounds' may only be applied along axis '{self.axis}', "
f"but was given axis '{axis}'."
)
return self.updated_copy(slab_bounds=tuple(bounds))
@staticmethod
def _canonicalize_vertices_and_bulges(
vertices: ArrayFloat2D, bulges: ArrayFloat1D
) -> tuple[ArrayFloat2D, ArrayFloat1D]:
"""Canonicalize vertices and bulges together: remove zero-length edges,
enforce CCW winding, and synchronize bulge array.
Parameters
----------
vertices : ArrayFloat2D
Shape (N, 2) polygon vertices.
bulges : ArrayFloat1D
Shape (N,) bulge values for each edge.
Returns
-------
tuple[ArrayFloat2D, ArrayFloat1D]
(canonical_vertices, canonical_bulges) both in CCW order
with zero-length edges removed.
"""
vertices = np.array(vertices)
bulges = np.array(bulges)
# Detect and remove zero-length edges
next_vertices = np.roll(vertices, shift=-1, axis=0)
edge_lengths = np.linalg.norm(vertices - next_vertices, axis=1)
is_zero_length = np.isclose(edge_lengths, 0, rtol=_IS_CLOSE_RTOL)
if np.any(is_zero_length):
# Drop zero-length edges (both vertex and bulge)
# Note: validation for non-zero bulges on zero-length edges
# is handled by _validate_no_bulge_on_zero_length_edges validator
has_nonzero_length = ~is_zero_length
vertices = vertices[has_nonzero_length]
bulges = bulges[has_nonzero_length]
# Check winding and reverse if CW; arc segments contribute to the signed
# area so winding stays well-defined even for collinear chord vertices
if PolySlab._signed_area_with_bulges(vertices, bulges) < 0:
vertices = vertices[::-1]
bulges = -np.roll(bulges[::-1], -1)
return vertices, bulges
@cached_property
def _canonical_vertices_and_bulges(self) -> tuple[ArrayFloat2D, ArrayFloat1D]:
"""Canonical (CCW, deduplicated) vertices and bulges as a synchronized pair."""
if self.bulges is None:
proper = self._proper_vertices(self.vertices)
return proper, np.zeros(proper.shape[0])
return self._canonicalize_vertices_and_bulges(self.vertices, self.bulges)
@property
def _bulges(self) -> ArrayFloat1D:
"""Postprocessed bulge values."""
return self._canonical_vertices_and_bulges[1]
@cached_property
def _has_arc_segments(self) -> bool:
"""Whether this PolySlab has any non-zero bulge (arc) segments."""
return not np.allclose(self._bulges, 0)
@cached_property
def _bulge_data(self) -> dict:
"""Precomputed arc geometry data from bulges and reference polygon.
Note: currently assumes all polygon cross-sections are identical
(sidewall_angle must be 0 when bulges are present). This will need
per-polygon bulge_data when nonzero sidewall angle is supported with bulges.
"""
return _PolyBulgeUtil._compute_bulge_data(self.reference_polygon, self._bulges)
@cached_property
def _discretized_reference_polygon(self) -> ArrayFloat2D:
"""Reference polygon with arc segments discretized to polylines.
For polygons without arc segments, returns the regular reference_polygon.
For polygons with arcs, returns a higher-resolution polyline approximation.
"""
if not self._has_arc_segments:
return self.reference_polygon
return _PolyBulgeUtil._polygon_discretize(self.reference_polygon, self._bulge_data)
@cached_property
def is_ccw(self) -> bool:
"""Is this ``PolySlab`` CCW-oriented?"""
return PolySlab._area(self.vertices) > 0
[docs]
def inside(self, x: NDArray[float], y: NDArray[float], z: NDArray[float]) -> NDArray[bool]:
"""For input arrays ``x``, ``y``, ``z`` of arbitrary but identical shape, return an array
with the same shape which is ``True`` for every point in zip(x, y, z) that is inside the
volume of the :class:`~tidy3d.Geometry`, and ``False`` otherwise.
Note
----
For slanted sidewalls, this function only works if x, y, and z are arrays produced by a
``meshgrid call``, i.e. 3D arrays and each is constant along one axis.
Parameters
----------
x : np.ndarray[float]
Array of point positions in x direction.
y : np.ndarray[float]
Array of point positions in y direction.
z : np.ndarray[float]
Array of point positions in z direction.
Returns
-------
np.ndarray[bool]
``True`` for every point that is inside the geometry.
"""
self._ensure_equal_shape(x, y, z)
z, (x, y) = self.pop_axis((x, y, z), axis=self.axis)
z0 = self.center_axis
dist_z = np.abs(z - z0)
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):
return inside_height
# check what points are inside polygon cross section (face)
z_local = z - z0 # distance to the middle
dist = -z_local * self._tanq
if isinstance(x, np.ndarray):
inside_polygon = np.zeros_like(inside_height)
xs_slab = x[inside_height]
ys_slab = y[inside_height]
# vertical sidewall
if math.isclose(self.sidewall_angle, 0):
# Use discretized polygon for arc segments
poly_vertices = self._discretized_reference_polygon
face_polygon = shapely.Polygon(poly_vertices).buffer(fp_eps)
shapely.prepare(face_polygon)
inside_polygon_slab = shapely.contains_xy(face_polygon, x=xs_slab, y=ys_slab)
inside_polygon[inside_height] = inside_polygon_slab
# slanted sidewall, offsetting vertices at each z
else:
# a helper function for moving axis
def _move_axis(arr: NDArray) -> NDArray:
return np.moveaxis(arr, source=self.axis, destination=-1)
def _move_axis_reverse(arr: NDArray) -> NDArray:
return np.moveaxis(arr, source=-1, destination=self.axis)
inside_polygon_axis = _move_axis(inside_polygon)
x_axis = _move_axis(x)
y_axis = _move_axis(y)
for z_i in range(z.shape[self.axis]):
if not _move_axis(inside_height)[0, 0, z_i]:
continue
vertices_z = self._shift_vertices(
self.middle_polygon, _move_axis(dist)[0, 0, z_i]
)[0]
face_polygon = shapely.Polygon(vertices_z).buffer(fp_eps)
shapely.prepare(face_polygon)
xs = x_axis[:, :, 0].flatten()
ys = y_axis[:, :, 0].flatten()
inside_polygon_slab = shapely.contains_xy(face_polygon, x=xs, y=ys)
inside_polygon_axis[:, :, z_i] = inside_polygon_slab.reshape(x_axis.shape[:2])
inside_polygon = _move_axis_reverse(inside_polygon_axis)
else:
if math.isclose(self.sidewall_angle, 0):
poly_vertices = self._discretized_reference_polygon
else:
poly_vertices = self._shift_vertices(self.middle_polygon, dist)[0]
face_polygon = self.make_shapely_polygon(poly_vertices).buffer(fp_eps)
point = shapely.Point(x, y)
inside_polygon = face_polygon.covers(point)
return inside_height * inside_polygon
@verify_packages_import(["trimesh"])
def _do_intersections_tilted_plane(
self,
normal: Coordinate,
origin: Coordinate,
to_2D: MatrixReal4x4,
quad_segs: int | None = None,
) -> list[Shapely]:
"""Return a list of shapely geometries at the plane specified by normal and origin.
Parameters
----------
normal : Coordinate
Vector defining the normal direction to the plane.
origin : Coordinate
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
-------
list[shapely.geometry.base.BaseGeometry]
List of 2D shapes that intersect plane.
For more details refer to
`Shapely's Documentation <https://shapely.readthedocs.io/en/stable/project.html>`_.
"""
import trimesh
if math.isclose(self.sidewall_angle, 0):
base_poly = self._discretized_reference_polygon
top_poly = base_poly
else:
base_poly = self.base_polygon
top_poly = self.top_polygon
if len(base_poly) > _MAX_POLYSLAB_VERTICES_FOR_TRIANGULATION:
log.warning(
f"Processing PolySlabs with over {_MAX_POLYSLAB_VERTICES_FOR_TRIANGULATION} vertices can be slow.",
log_once=True,
)
base_triangles = triangulation.triangulate(base_poly)
top_triangles = (
base_triangles
if math.isclose(self.sidewall_angle, 0)
else triangulation.triangulate(top_poly)
)
n = len(base_poly)
faces = (
[[a, b, c] for c, b, a in base_triangles]
+ [[n + a, n + b, n + c] for a, b, c in top_triangles]
+ [(i, (i + 1) % n, n + i) for i in range(n)]
+ [((i + 1) % n, n + ((i + 1) % n), n + i) for i in range(n)]
)
x = np.hstack((base_poly[:, 0], top_poly[:, 0]))
y = np.hstack((base_poly[:, 1], top_poly[:, 1]))
z = np.hstack((np.full(n, self.slab_bounds[0]), np.full(n, self.slab_bounds[1])))
vertices = np.vstack(self.unpop_axis(z, (x, y), self.axis)).T
mesh = trimesh.Trimesh(vertices, faces)
section = mesh.section(plane_origin=origin, plane_normal=normal)
if section is None:
return []
path, _ = section.to_2D(to_2D=to_2D)
return path.polygons_full
def _intersections_normal(self, z: float, quad_segs: int | None = None) -> list[Shapely]:
"""Find shapely geometries intersecting planar geometry with axis normal to slab.
Parameters
----------
z : float
Position along the axis normal to slab.
Returns
-------
list[shapely.geometry.base.BaseGeometry]
List of 2D shapes that intersect plane.
For more details refer to
`Shapely's Documentation <https://shapely.readthedocs.io/en/stable/project.html>`_.
"""
if math.isclose(self.sidewall_angle, 0):
# Use discretized polygon for arc segments
poly_vertices = self._discretized_reference_polygon
return [self.make_shapely_polygon(poly_vertices)]
z0 = self.center_axis
z_local = z - z0 # distance to the middle
dist = -z_local * self._tanq
vertices_z = self._shift_vertices(self.middle_polygon, dist)[0]
return [self.make_shapely_polygon(vertices_z)]
def _intersections_side(self, position: float, axis: int) -> list[Shapely]:
"""Find shapely geometries intersecting planar geometry with axis orthogonal to slab.
For slanted polyslab, the procedure is as follows,
1) Find out all z-coordinates where the plane will intersect directly with a vertex.
Denote the coordinates as (z_0, z_1, z_2, ... )
2) Find out all polygons that can be formed between z_i and z_{i+1}. There are two
types of polygons:
a) formed by the plane intersecting the edges
b) formed by the plane intersecting the vertices.
For either type, one needs to compute:
i) intersecting position
ii) angle between the plane and the intersecting edge
For a), both are straightforward to compute; while for b), one needs to compute
which edge the plane will slide into.
3) Looping through z_i, and merge all polygons. The partition by z_i is because once
the plane intersects the vertex, it can intersect with other edges during
the extrusion.
Parameters
----------
position : float
Position along ``axis``.
axis : int
Integer index into 'xyz' (0,1,2).
Returns
-------
list[shapely.geometry.base.BaseGeometry]
List of 2D shapes that intersect plane.
For more details refer to
`Shapely's Documentation <https://shapely.readthedocs.io/en/stable/project.html>`_.
"""
# find out all z_i where the plane will intersect the vertex
z0 = self.center_axis
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.finite_length_axis)
h_base = 0.0
for h_top in height_list:
# length within between top and bottom
h_length = h_top - h_base
# coordinate of each subsection
z_min = z_base + h_base
z_max = np.inf if np.isposinf(h_top) else z_base + h_top
# for vertical sidewall, no need for complications
if math.isclose(self.sidewall_angle, 0):
ints_y, ints_angle = self._find_intersecting_ys_angle_vertical(
self._discretized_reference_polygon, position, axis_ordered
)
else:
if self._has_arc_segments:
raise NotImplementedError(
"Arc segments are not supported with slanted sidewalls in "
"_intersections_side. Update _shift_vertices and this method to handle "
"arc geometry before lifting the sidewall_angle restriction."
)
# for slanted sidewall, move up by `fp_eps` in case vertices are degenerate at the base.
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
)
# make polygon with intersections and z axis information
for y_index in range(len(ints_y) // 2):
y_min = ints_y[2 * y_index]
y_max = ints_y[2 * y_index + 1]
minx, miny = self._order_by_axis(plane_val=y_min, axis_val=z_min, axis=axis)
maxx, maxy = self._order_by_axis(plane_val=y_max, axis_val=z_max, axis=axis)
if math.isclose(self.sidewall_angle, 0):
polys.append(self.make_shapely_box(minx, miny, maxx, maxy))
else:
angle_min = ints_angle[2 * y_index]
angle_max = ints_angle[2 * y_index + 1]
angle_min = np.arctan(np.tan(self.sidewall_angle) / np.sin(angle_min))
angle_max = np.arctan(np.tan(self.sidewall_angle) / np.sin(angle_max))
dy_min = h_length * np.tan(angle_min)
dy_max = h_length * np.tan(angle_max)
x1, y1 = self._order_by_axis(plane_val=y_min, axis_val=z_min, axis=axis)
x2, y2 = self._order_by_axis(plane_val=y_max, axis_val=z_min, axis=axis)
x3, y3 = self._order_by_axis(
plane_val=y_max - dy_max, axis_val=z_max, axis=axis
)
x4, y4 = self._order_by_axis(
plane_val=y_min + dy_min, axis_val=z_max, axis=axis
)
vertices = ((x1, y1), (x2, y2), (x3, y3), (x4, y4))
polys.append(self.make_shapely_polygon(vertices).buffer(0))
# update the base coordinate for the next subsection
h_base = h_top
# merge touching polygons
polys_union = shapely.unary_union(polys, grid_size=base.POLY_GRID_SIZE)
if polys_union.geom_type == "Polygon":
return [polys_union]
if polys_union.geom_type == "MultiPolygon":
return polys_union.geoms
# in other cases, just return the original unmerged polygons
return polys
def _find_intersecting_height(self, position: float, axis: int) -> NDArray:
"""Found a list of height where the plane will intersect with the vertices;
For vertical sidewall, just return np.array([]).
Assumes axis is handles so this function works on xy plane.
Parameters
----------
position : float
position along axis.
axis : int
Integer index into 'xyz' (0,1,2).
Returns
np.ndarray
Height (relative to the base) where the plane will intersect with vertices.
"""
if math.isclose(self.sidewall_angle, 0):
return np.array([])
if self._has_arc_segments:
raise NotImplementedError(
"Arc segments are not supported with slanted sidewalls in "
"_find_intersecting_height. Update this method to use discretized "
"arc geometry before lifting the sidewall_angle restriction."
)
# shift rate
dist = 1.0
shift_x, shift_y = PolySlab._shift_vertices(self.middle_polygon, dist)[2]
shift_val = shift_x if axis == 0 else shift_y
shift_val[np.isclose(shift_val, 0, rtol=_IS_CLOSE_RTOL)] = np.inf # for static vertices
# distance to the plane in the direction of vertex shifting
distance = self.middle_polygon[:, axis] - position
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.finite_length_axis - fp_eps]
return height
def _find_intersecting_ys_angle_vertical(
self,
vertices: NDArray,
position: float,
axis: int,
exclude_on_vertices: bool = False,
) -> tuple[NDArray, NDArray]:
"""Finds pairs of forward and backwards vertices where polygon intersects position at axis,
Find intersection point (in y) assuming straight line,and intersecting angle between plane
and edges. (For unslanted polyslab).
Assumes axis is handles so this function works on xy plane.
Parameters
----------
vertices : np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
position : float
position along axis.
axis : int
Integer index into 'xyz' (0,1,2).
exclude_on_vertices : bool = False
Whether to exclude those intersecting directly with the vertices.
Returns
-------
Union[np.ndarray, np.ndarray]
List of intersection points along y direction.
List of angles between plane and edges.
"""
vertices_axis = vertices
# flip vertices x,y for axis = y
if axis == 1:
vertices_axis = np.roll(vertices_axis, shift=1, axis=1)
# get the forward vertices
vertices_f = np.roll(vertices_axis, shift=-1, axis=0)
# x coordinate of the two sets of vertices
x_vertices_f, _ = vertices_f.T
x_vertices_axis, _ = vertices_axis.T
# 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)
if exclude_on_vertices:
# 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]
# intersecting positions and angles
ints_y = []
ints_angle = []
for vertices_f_local, vertices_b_local in zip(iverts_b, iverts_f):
x1, y1 = vertices_f_local
x2, y2 = vertices_b_local
slope = (y2 - y1) / (x2 - x1)
y = y1 + slope * (position - x1)
ints_y.append(y)
ints_angle.append(np.pi / 2 - np.arctan(np.abs(slope)))
ints_y = np.array(ints_y)
ints_angle = np.array(ints_angle)
# 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[sort_index]
return ints_y_sort, ints_angle_sort
def _find_intersecting_ys_angle_slant(
self, vertices: NDArray, position: float, axis: int
) -> tuple[NDArray, NDArray]:
"""Finds pairs of forward and backwards vertices where polygon intersects position at axis,
Find intersection point (in y) assuming straight line,and intersecting angle between plane
and edges. (For slanted polyslab)
Assumes axis is handles so this function works on xy plane.
Parameters
----------
vertices : np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
position : float
position along axis.
axis : int
Integer index into 'xyz' (0,1,2).
Returns
-------
Union[np.ndarray, np.ndarray]
List of intersection points along y direction.
List of angles between plane and edges.
"""
if self._has_arc_segments:
raise NotImplementedError(
"Arc segments are not supported with slanted sidewalls in "
"_find_intersecting_ys_angle_slant. Update this method to use "
"discretized arc geometry before lifting the sidewall_angle restriction."
)
vertices_axis = vertices.copy()
# flip vertices x,y for axis = y
if axis == 1:
vertices_axis = np.roll(vertices_axis, shift=1, axis=1)
# get the forward vertices
vertices_f = np.roll(vertices_axis, shift=-1, axis=0)
# get the backward vertices
vertices_b = np.roll(vertices_axis, shift=1, axis=0)
## First part, plane intersects with edges, same as vertical
ints_y, ints_angle = self._find_intersecting_ys_angle_vertical(
vertices, position, axis, exclude_on_vertices=True
)
ints_y = ints_y.tolist()
ints_angle = ints_angle.tolist()
## Second part, plane intersects directly with vertices
# vertices on the intersection
intersects_on = np.isclose(vertices_axis[:, 0], position, rtol=_IS_CLOSE_RTOL)
iverts_on = vertices_axis[intersects_on]
# position of the neighbouring vertices
iverts_b = vertices_b[intersects_on]
iverts_f = vertices_f[intersects_on]
# shift rate
dist = -np.sign(self.sidewall_angle)
shift_x, shift_y = self._shift_vertices(self.middle_polygon, dist)[2]
shift_val = shift_x if axis == 0 else shift_y
shift_val = shift_val[intersects_on]
for vertices_f_local, vertices_b_local, vertices_on_local, shift_local in zip(
iverts_f, iverts_b, iverts_on, shift_val
):
x_on, y_on = vertices_on_local
x_f, y_f = vertices_f_local
x_b, y_b = vertices_b_local
num_added = 0 # keep track the number of added vertices
slope = [] # list of slopes for added vertices
# case 1, shifting velocity is 0
if np.isclose(shift_local, 0, rtol=_IS_CLOSE_RTOL):
ints_y.append(y_on)
# Slope w.r.t. forward and backward should equal,
# just pick one of them.
slope.append((y_on - y_b) / (x_on - x_b))
ints_angle.append(np.pi / 2 - np.arctan(np.abs(slope[0])))
continue
# case 2, shifting towards backward direction
if (x_b - position) * shift_local < 0:
ints_y.append(y_on)
slope.append((y_on - y_b) / (x_on - x_b))
num_added += 1
# case 3, shifting towards forward direction
if (x_f - position) * shift_local < 0:
ints_y.append(y_on)
slope.append((y_on - y_f) / (x_on - x_f))
num_added += 1
# in case 2, and case 3, if just num_added = 1
if num_added == 1:
ints_angle.append(np.pi / 2 - np.arctan(np.abs(slope[0])))
# if num_added = 2, the order of the two new vertices needs to handled correctly;
# it should be sorted according to the -slope * moving direction
elif num_added == 2:
dressed_slope = [-s_i * shift_local for s_i in slope]
sort_index = np.argsort(np.array(dressed_slope))
sorted_slope = np.array(slope)[sort_index]
ints_angle.append(np.pi / 2 - np.arctan(np.abs(sorted_slope[0])))
ints_angle.append(np.pi / 2 - np.arctan(np.abs(sorted_slope[1])))
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]
ints_angle_sort = ints_angle[sort_index]
return ints_y_sort, ints_angle_sort
@cached_property
def bounds(self) -> Bound:
"""Returns bounding box min and max coordinates. The dilation and slant angle are not
taken into account exactly for speed. Instead, the polygon may be slightly smaller than
the returned bounds, but it should always be fully contained.
Returns
-------
tuple[float, float, float], tuple[float, float float]
Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
"""
# check for the maximum possible contribution from dilation/slant on each side
max_offset = self.dilation
# sidewall_angle may be autograd-traced; unbox for this check
if not math.isclose(getval(self.sidewall_angle), 0):
if self.reference_plane == "bottom":
max_offset += max(0, -self._tanq * self.finite_length_axis)
elif self.reference_plane == "top":
max_offset += max(0, self._tanq * self.finite_length_axis)
elif self.reference_plane == "middle":
max_offset += max(0, abs(self._tanq) * self.finite_length_axis / 2)
# special care when dilated
if max_offset > 0:
dilated_vertices = self._shift_vertices(
self._proper_vertices(self.vertices), max_offset
)[0]
xmin, ymin = np.amin(dilated_vertices, axis=0)
xmax, ymax = np.amax(dilated_vertices, axis=0)
else:
# otherwise, bounds are directly based on the supplied vertices
xmin, ymin = np.amin(self.vertices, axis=0)
xmax, ymax = np.amax(self.vertices, axis=0)
# Account for arc segments extending beyond chord endpoints
if self._has_arc_segments:
vertices = self._canonical_vertices_and_bulges[0]
arc_mins, arc_maxs = _PolyBulgeUtil._polygon_arcs_bounds(vertices, self._bulge_data)
xmin = min(xmin, np.min(arc_mins[:, 0]))
ymin = min(ymin, np.min(arc_mins[:, 1]))
xmax = max(xmax, np.max(arc_maxs[:, 0]))
ymax = max(ymax, np.max(arc_maxs[:, 1]))
# get bounds in (local) z
zmin, zmax = self.slab_bounds
# rearrange axes
coords_min = self.unpop_axis(zmin, (xmin, ymin), axis=self.axis)
coords_max = self.unpop_axis(zmax, (xmax, ymax), axis=self.axis)
return (tuple(coords_min), tuple(coords_max))
def _extrusion_length_to_offset_distance(self, extrusion: float) -> float:
"""Convert extrusion length to offset distance."""
if math.isclose(self.sidewall_angle, 0):
return 0
return -extrusion * self._tanq
@staticmethod
def _area(vertices: NDArray) -> float:
"""Compute the signed polygon area (positive for CCW orientation).
Parameters
----------
vertices : np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
Returns
-------
float
Signed polygon area (positive for CCW orientation).
"""
vert_shift = np.roll(vertices, axis=0, shift=-1)
xs, ys = vertices.T
xs_shift, ys_shift = vert_shift.T
term1 = xs * ys_shift
term2 = ys * xs_shift
return np.sum(term1 - term2) * 0.5
@staticmethod
def _signed_area_with_bulges(vertices: NDArray, bulges: ArrayFloat1D) -> float:
"""Compute the signed polygon area including arc segment contributions.
Parameters
----------
vertices : np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
bulges : ArrayFloat1D
Shape (N,) bulge values for each edge.
Returns
-------
float
Signed polygon area (positive for CCW orientation).
"""
area = PolySlab._area(vertices)
bulge_data = _PolyBulgeUtil._compute_bulge_data(vertices, bulges)
if np.any(bulge_data["arc_mask"]):
area += _PolyBulgeUtil._arc_segment_area(bulge_data)
return area
@staticmethod
def _perimeter(vertices: NDArray) -> float:
"""Compute the polygon perimeter.
Parameters
----------
vertices : np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
Returns
-------
float
Polygon perimeter.
"""
vert_shift = np.roll(vertices, axis=0, shift=-1)
squared_diffs = (vertices - vert_shift) ** 2
# distance along each edge
dists = np.sqrt(squared_diffs.sum(axis=-1))
# total distance along all edges
return np.sum(dists)
@staticmethod
def _orient(vertices: NDArray) -> NDArray:
"""Return a CCW-oriented polygon.
Parameters
----------
vertices : np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
Returns
-------
np.ndarray
Vertices of a CCW-oriented polygon.
"""
return vertices if PolySlab._area(vertices) > 0 else vertices[::-1, :]
@staticmethod
def _remove_duplicate_vertices(vertices: NDArray) -> NDArray:
"""Remove redundant/identical nearest neighbour vertices.
Parameters
----------
vertices : np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
Returns
-------
np.ndarray
Vertices of polygon.
"""
return remove_adjacent_duplicate_vertices(vertices, rtol=_IS_CLOSE_RTOL)
@staticmethod
def _proper_vertices(vertices: ArrayFloat2D) -> NDArray:
"""convert vertices to np.array format,
removing duplicate neighbouring vertices,
and oriented in CCW direction.
Returns
-------
ArrayLike[float, float]
The vertices of the polygon for internal use.
"""
vertices_np = np.array(vertices)
return PolySlab._orient(PolySlab._remove_duplicate_vertices(vertices_np))
@staticmethod
def _edge_events_detection(
proper_vertices: NDArray, dilation: float, ignore_at_dist: bool = True
) -> bool:
"""Detect any edge events within the offset distance ``dilation``.
If ``ignore_at_dist=True``, the edge event at ``dist`` is ignored.
"""
# ignore the event that occurs right at the offset distance
if ignore_at_dist:
dilation -= fp_eps * dilation / abs(dilation)
# number of vertices before offsetting
num_vertices = proper_vertices.shape[0]
# 0) fully eroded?
if dilation < 0 and dilation < -PolySlab._maximal_erosion(proper_vertices):
return True
# sample at a few dilation values
dist_list = (
dilation
* np.linspace(
0, 1, 1 + _N_SAMPLE_POLYGON_INTERSECT, dtype=config.adjoint.gradient_dtype_float
)[1:]
)
for dist in dist_list:
# offset: we offset the vertices first, and then use shapely to make it proper
# in principle, one can offset with shapely.buffer directly, but shapely somehow
# automatically removes some vertices even though no change of topology.
poly_offset = PolySlab._shift_vertices(proper_vertices, dist)[0]
# flipped winding number
if PolySlab._area(poly_offset) < fp_eps**2:
return True
poly_offset = shapely.make_valid(PolySlab.make_shapely_polygon(poly_offset))
# 1) polygon split or create holes/islands
if not poly_offset.geom_type == "Polygon" or len(poly_offset.interiors) > 0:
return True
# 2) reduction in vertex number
offset_vertices = PolySlab._proper_vertices(poly_offset.exterior.coords)
if offset_vertices.shape[0] != num_vertices:
return True
# 3) some split polygon might fully disappear after the offset, but they
# can be detected if we offset back.
poly_offset_back = shapely.make_valid(
PolySlab.make_shapely_polygon(PolySlab._shift_vertices(offset_vertices, -dist)[0])
)
if poly_offset_back.geom_type == "MultiPolygon" or len(poly_offset_back.interiors) > 0:
return True
offset_back_vertices = poly_offset_back.exterior.coords
if PolySlab._proper_vertices(offset_back_vertices).shape[0] != num_vertices:
return True
return False
@staticmethod
def _neighbor_vertices_crossing_detection(
vertices: NDArray, dist: float, ignore_at_dist: bool = True
) -> float:
"""Detect if neighboring vertices will cross after a dilation distance dist.
Parameters
----------
vertices : np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
dist : float
Distance to offset.
ignore_at_dist : bool, optional
whether to ignore the event right at ``dist`.
Returns
-------
float
the absolute value of the maximal allowed dilation
if there are any crossing, otherwise return ``None``.
"""
# ignore the event that occurs right at the offset distance
if ignore_at_dist:
dist -= fp_eps * dist / abs(dist)
edge_length, edge_reduction = PolySlab._edge_length_and_reduction_rate(vertices)
length_remaining = edge_length - edge_reduction * dist
if np.any(length_remaining < 0):
index_oversized = length_remaining < 0
max_dist = np.min(
np.abs(edge_length[index_oversized] / edge_reduction[index_oversized])
)
return max_dist
return None
[docs]
@staticmethod
def array_to_vertices(arr_vertices: NDArray) -> ArrayFloat2D:
"""Converts a numpy array of vertices to a list of tuples."""
return list(arr_vertices)
[docs]
@staticmethod
def vertices_to_array(vertices_tuple: ArrayFloat2D) -> NDArray:
"""Converts a list of tuples (vertices) to a numpy array."""
return np.array(vertices_tuple)
@cached_property
def interior_angle(self) -> ArrayFloat1D:
"""Angle formed inside polygon by two adjacent edges."""
if self._has_arc_segments:
raise NotImplementedError(
"interior_angle does not account for arc segment geometry. "
"The angle at each vertex depends on the arc tangent directions, "
"not the chord directions."
)
def normalize(v: NDArray) -> NDArray:
return v / np.linalg.norm(v, axis=0)
vs_orig = self.reference_polygon.T
vs_next = np.roll(vs_orig, axis=-1, shift=-1)
vs_previous = np.roll(vs_orig, axis=-1, shift=+1)
asp = normalize(vs_next - vs_orig)
asm = normalize(vs_previous - vs_orig)
cos_angle = asp[0] * asm[0] + asp[1] * asm[1]
sin_angle = asp[0] * asm[1] - asp[1] * asm[0]
angle = np.arccos(cos_angle)
# concave angles
angle[sin_angle < 0] = 2 * np.pi - angle[sin_angle < 0]
return angle
@staticmethod
def _shift_vertices(
vertices: NDArray, dist: float
) -> tuple[NDArray, NDArray, tuple[NDArray, NDArray]]:
"""Shifts the vertices of a polygon outward uniformly by distances
`dists`.
Parameters
----------
np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
dist : float
Distance to offset.
Returns
-------
tuple[np.ndarray, np.narray,tuple[np.ndarray,np.ndarray]]
New polygon vertices;
and the shift of vertices in direction parallel to the edges.
Shift along x and y direction.
"""
# 'dist' may be autograd-traced; unbox for the zero-check only
if math.isclose(getval(dist), 0):
return vertices, np.zeros(vertices.shape[0], dtype=float), None
def rot90(v: tuple[NDArray, NDArray]) -> NDArray:
"""90 degree rotation of 2d vector
vx -> vy
vy -> -vx
"""
vxs, vys = v
return np.stack((-vys, vxs), axis=0)
def cross(u: NDArray, v: NDArray) -> Any:
return u[0] * v[1] - u[1] * v[0]
def normalize(v: NDArray) -> NDArray:
return v / np.linalg.norm(v, axis=0)
vs_orig = copy(vertices.T)
vs_next = np.roll(copy(vs_orig), axis=-1, shift=-1)
vs_previous = np.roll(copy(vs_orig), axis=-1, shift=+1)
asp = normalize(vs_next - vs_orig)
asm = normalize(vs_orig - vs_previous)
# the vertex shift is decomposed into parallel and perpendicular directions
perpendicular_shift = -dist
det = cross(asm, asp)
tan_half_angle = np.where(
np.isclose(det, 0, rtol=_IS_CLOSE_RTOL),
0.0,
cross(asm, rot90(asm - asp)) / (det + np.isclose(det, 0, rtol=_IS_CLOSE_RTOL)),
)
parallel_shift = dist * tan_half_angle
shift_total = perpendicular_shift * rot90(asm) + parallel_shift * asm
shift_x = shift_total[0, :]
shift_y = shift_total[1, :]
return (
np.swapaxes(vs_orig + shift_total, -2, -1),
parallel_shift,
(shift_x, shift_y),
)
@staticmethod
def _edge_length_and_reduction_rate(
vertices: NDArray,
) -> tuple[NDArray, NDArray]:
"""Edge length of reduction rate of each edge with unit offset length.
Parameters
----------
vertices : np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
Returns
-------
tuple[np.ndarray, np.narray]
edge length, and reduction rate
"""
# edge length
vs_orig = copy(vertices.T)
vs_next = np.roll(copy(vs_orig), axis=-1, shift=-1)
edge_length = np.linalg.norm(vs_next - vs_orig, axis=0)
# edge length remaining
dist = 1
parallel_shift = PolySlab._shift_vertices(vertices, dist)[1]
parallel_shift_p = np.roll(copy(parallel_shift), shift=-1)
edge_reduction = -(parallel_shift + parallel_shift_p)
return edge_length, edge_reduction
@staticmethod
def _maximal_erosion(vertices: NDArray) -> float:
"""The erosion value that reduces the length of
all edges to be non-positive.
"""
edge_length, edge_reduction = PolySlab._edge_length_and_reduction_rate(vertices)
ind_nonzero = abs(edge_reduction) > fp_eps
return -np.min(edge_length[ind_nonzero] / edge_reduction[ind_nonzero])
@staticmethod
def _heal_polygon(vertices: NDArray) -> NDArray:
"""heal a self-intersecting polygon."""
shapely_poly = PolySlab.make_shapely_polygon(vertices)
if shapely_poly.is_valid:
return vertices
elif hasbox(vertices):
raise NotImplementedError(
"The dilation caused damage to the polygon. "
"Automatically healing this is currently not supported when "
"differentiating w.r.t. the vertices. Try increasing the spacing "
"between vertices or reduce the amount of dilation."
)
# perform healing
poly_heal = shapely.make_valid(shapely_poly)
return PolySlab._proper_vertices(list(poly_heal.exterior.coords))
def _volume(self, bounds: Bound) -> float:
"""Returns object's volume within given bounds."""
z_min, z_max = self.slab_bounds
z_min = max(z_min, bounds[0][self.axis])
z_max = min(z_max, bounds[1][self.axis])
length = z_max - z_min
# Use arc-aware area calculation when arc segments present
if self._has_arc_segments:
arc_correction = _PolyBulgeUtil._arc_segment_area(self._bulge_data)
top_area = abs(self._area(self.top_polygon) + arc_correction)
base_area = abs(self._area(self.base_polygon) + arc_correction)
else:
top_area = abs(self._area(self.top_polygon))
base_area = abs(self._area(self.base_polygon))
# https://mathworld.wolfram.com/PyramidalFrustum.html
return 1.0 / 3.0 * length * (top_area + base_area + np.sqrt(top_area * base_area))
def _surface_area(self, bounds: Bound) -> float:
"""Returns object's surface area within given bounds."""
area = 0
top_polygon = self.top_polygon
base_polygon = self.base_polygon
# Use arc-aware calculations when arc segments present
if self._has_arc_segments:
arc_correction = _PolyBulgeUtil._arc_segment_area(self._bulge_data)
top_area = abs(self._area(top_polygon) + arc_correction)
base_area = abs(self._area(base_polygon) + arc_correction)
top_perim = _PolyBulgeUtil._polygon_perimeter(top_polygon, self._bulge_data)
base_perim = _PolyBulgeUtil._polygon_perimeter(base_polygon, self._bulge_data)
else:
top_area = abs(self._area(top_polygon))
base_area = abs(self._area(base_polygon))
top_perim = self._perimeter(top_polygon)
base_perim = self._perimeter(base_polygon)
z_min, z_max = self.slab_bounds
if z_min < bounds[0][self.axis]:
z_min = bounds[0][self.axis]
else:
area += base_area
if z_max > bounds[1][self.axis]:
z_max = bounds[1][self.axis]
else:
area += top_area
length = z_max - z_min
area += 0.5 * (top_perim + base_perim) * length
return area
""" Autograd code """
def _resolve_autograd_route(self, field_path: tuple[Any, ...]) -> AutogradRoute:
"""Resolve and validate one traced PolySlab path for adjoint routing."""
if self._has_arc_segments:
raise AdjointError(
"Automatic differentiation with respect to PolySlab geometry parameters is not "
"supported when vertices include non-zero bulge values. Use straight-edged "
"polyslabs for inverse design optimization."
)
return super()._resolve_autograd_route(field_path)
def _compute_derivatives(self, derivative_info: DerivativeInfo) -> AutogradFieldMap:
"""
Return VJPs while handling several edge-cases:
- If the slab volume does not overlap the simulation, all grads are zero
(one warning is issued).
- Faces that lie completely outside the simulation give zero ``slab_bounds``
gradients; this includes the +/- inf cases.
- A 2d simulation collapses the surface integral to a line integral
"""
vjps: AutogradFieldMap = {}
intersect_min, intersect_max = map(np.asarray, derivative_info.bounds_intersect)
sim_min, sim_max = map(np.asarray, derivative_info.simulation_bounds)
extents = intersect_max - intersect_min
is_2d = np.isclose(extents[self.axis], 0.0)
# early return if polyslab is not in simulation domain
slab_min, slab_max = self.slab_bounds
if (slab_max < sim_min[self.axis]) or (slab_min > sim_max[self.axis]):
log.warning(
"'PolySlab' lies completely outside the simulation domain.",
log_once=True,
)
for p in derivative_info.paths:
vjps[p] = np.zeros_like(self.vertices) if p == ("vertices",) else 0.0
return vjps
# create interpolators once for ALL derivative computations
# use provided interpolators if available to avoid redundant field data conversions
interpolators = derivative_info.interpolators or derivative_info.create_interpolators(
dtype=config.adjoint.gradient_dtype_float
)
for path in derivative_info.paths:
if path == ("vertices",):
vjps[path] = self._compute_derivative_vertices(
derivative_info, sim_min, sim_max, is_2d, interpolators
)
elif path == ("sidewall_angle",):
vjps[path] = self._compute_derivative_sidewall_angle(
derivative_info, sim_min, sim_max, is_2d, interpolators
)
elif path[0] == "slab_bounds":
idx = path[1]
face_coord = self.slab_bounds[idx]
# face entirely outside -> gradient 0
if (
np.isinf(face_coord)
or face_coord < sim_min[self.axis]
or face_coord > sim_max[self.axis]
or is_2d
):
vjps[path] = 0.0
continue
v = self._compute_derivative_slab_bounds(derivative_info, idx, interpolators)
# outward-normal convention
if idx == 0:
v *= -1
vjps[path] = v
return vjps
# ---- Shared helpers for VJP surface integrations ----
def _z_slices(
self, sim_min: NDArray, sim_max: NDArray, is_2d: bool, dx: float
) -> tuple[NDArray, float, float, float]:
"""Compute z-slice centers and spacing within bounds.
Returns (z_centers, dz, z0, z1). For 2D, returns single center and dz=1.
"""
if is_2d:
midpoint_z = np.maximum(
np.minimum(self.center_axis, sim_max[self.axis]),
sim_min[self.axis],
)
zc = np.array([midpoint_z], dtype=config.adjoint.gradient_dtype_float)
return zc, 1.0, self.center_axis, self.center_axis
z0 = max(self.slab_bounds[0], sim_min[self.axis])
z1 = min(self.slab_bounds[1], sim_max[self.axis])
if z1 <= z0:
return np.array([], dtype=config.adjoint.gradient_dtype_float), 0.0, z0, z1
n_z = max(1, int(np.ceil((z1 - z0) / dx)))
dz = (z1 - z0) / n_z
z_centers = np.linspace(
z0 + dz / 2, z1 - dz / 2, n_z, dtype=config.adjoint.gradient_dtype_float
)
return z_centers, dz, z0, z1
@staticmethod
def _clip_edges_to_bounds_batch(
segment_starts: NDArray,
segment_ends: NDArray,
sim_min: NDArray,
sim_max: NDArray,
*,
_edge_clip_tol: float | None = None,
_dtype: type | None = None,
) -> tuple[NDArray, NDArray, NDArray]:
"""
Compute parametric bounds for multiple segments clipped to simulation bounds.
Parameters
----------
segment_starts : NDArray
(N, 3) array of segment start coordinates.
segment_ends : NDArray
(N, 3) array of segment end coordinates.
sim_min : NDArray
(3,) array of simulation minimum bounds.
sim_max : NDArray
(3,) array of simulation maximum bounds.
Returns
-------
is_within_bounds : NDArray
(N,) boolean array indicating if the segment intersects the bounds.
t_starts : NDArray
(N,) array of parametric start values (0.0 to 1.0).
t_ends : NDArray
(N,) array of parametric end values (0.0 to 1.0).
"""
n = segment_starts.shape[0]
if _edge_clip_tol is None:
_edge_clip_tol = config.adjoint.edge_clip_tolerance
if _dtype is None:
_dtype = config.adjoint.gradient_dtype_float
t_starts = np.zeros(n, dtype=_dtype)
t_ends = np.ones(n, dtype=_dtype)
is_within_bounds = np.ones(n, dtype=bool)
for dim in range(3):
start_coords = segment_starts[:, dim]
end_coords = segment_ends[:, dim]
bound_min = sim_min[dim]
bound_max = sim_max[dim]
# check for parallel edges (faster than isclose)
parallel = np.abs(start_coords - end_coords) < 1e-12
# parallel edges: check if outside bounds
outside = parallel & (
(start_coords < (bound_min - _edge_clip_tol))
| (start_coords > (bound_max + _edge_clip_tol))
)
is_within_bounds &= ~outside
# non-parallel edges: compute t_min, t_max
not_parallel = ~parallel & is_within_bounds
if np.any(not_parallel):
denom = np.where(not_parallel, end_coords - start_coords, 1.0) # avoid div by zero
t_min = (bound_min - start_coords) / denom
t_max = (bound_max - start_coords) / denom
# swap if needed
swap = t_min > t_max
t_min_new = np.where(swap, t_max, t_min)
t_max_new = np.where(swap, t_min, t_max)
# update t_starts and t_ends for valid non-parallel edges
t_starts = np.where(not_parallel, np.maximum(t_starts, t_min_new), t_starts)
t_ends = np.where(not_parallel, np.minimum(t_ends, t_max_new), t_ends)
# still valid?
is_within_bounds &= ~not_parallel | (t_starts < t_ends)
is_within_bounds &= t_ends > t_starts + _edge_clip_tol
return is_within_bounds, t_starts, t_ends
@staticmethod
def _adaptive_edge_samples(
L: float,
dx: float,
t_start: float = 0.0,
t_end: float = 1.0,
*,
_sample_fraction: float | None = None,
_gauss_order: int | None = None,
_dtype: type | None = None,
) -> tuple[NDArray, NDArray]:
"""
Compute Gauss samples and weights along [t_start, t_end] with adaptive count.
Parameters
----------
L : float
Physical length of the full edge.
dx : float
Target discretization step size.
t_start : float, optional
Start parameter, by default 0.0.
t_end : float, optional
End parameter, by default 1.0.
Returns
-------
tuple[NDArray, NDArray]
Tuple of (samples, weights) for the integration.
"""
if _sample_fraction is None:
_sample_fraction = config.adjoint.quadrature_sample_fraction
if _gauss_order is None:
_gauss_order = config.adjoint.gauss_quadrature_order
if _dtype is None:
_dtype = config.adjoint.gradient_dtype_float
L_eff = L * max(0.0, t_end - t_start)
n_uniform = max(1, int(np.ceil(L_eff / dx)))
n_gauss = max(2, n_uniform) if n_uniform <= 3 else max(2, int(n_uniform * _sample_fraction))
if n_gauss <= _gauss_order:
g, w = leggauss(n_gauss)
half_range = 0.5 * (t_end - t_start)
s = (half_range * g + 0.5 * (t_end + t_start)).astype(_dtype, copy=False)
wt = (w * half_range).astype(_dtype, copy=False)
return s, wt
# composite Gauss with fixed local order
g_loc, w_loc = leggauss(_gauss_order)
segs = max(1, int(np.ceil(n_gauss / _gauss_order)))
edges_t = np.linspace(t_start, t_end, segs + 1, dtype=_dtype)
# compute all segments at once
a = edges_t[:-1] # (segs,)
b = edges_t[1:] # (segs,)
half_width = 0.5 * (b - a) # (segs,)
mid = 0.5 * (b + a) # (segs,)
# (segs, 1) * (order,) + (segs, 1) -> (segs, order)
S = (half_width[:, None] * g_loc + mid[:, None]).astype(_dtype, copy=False)
W = (half_width[:, None] * w_loc).astype(_dtype, copy=False)
return S.ravel(), W.ravel()
def _collect_sidewall_patches(
self,
vertices: NDArray,
next_v: NDArray,
edges: NDArray,
basis: dict,
sim_min: NDArray,
sim_max: NDArray,
is_2d: bool,
dx: float,
) -> dict:
"""
Collect sidewall patch geometry for batched VJP evaluation.
Parameters
----------
vertices : NDArray
Array of polygon vertices.
next_v : NDArray
Array of next vertices (forming edges).
edges : NDArray
Edge vectors.
basis : dict
Basis vectors dictionary.
sim_min : NDArray
Simulation minimum bounds.
sim_max : NDArray
Simulation maximum bounds.
is_2d : bool
Whether the simulation is 2D.
dx : float
Discretization step.
Returns
-------
dict
Dictionary containing:
- centers: (N, 3) array of patch centers.
- normals: (N, 3) array of patch normals.
- perps1: (N, 3) array of first tangent vectors.
- perps2: (N, 3) array of second tangent vectors.
- Ls: (N,) array of edge lengths.
- s_vals: (N,) array of parametric coordinates along the edge.
- s_weights: (N,) array of quadrature weights.
- zc_vals: (N,) array of z-coordinates.
- dz: float, slice thickness.
- edge_indices: (N,) array of original edge indices.
"""
# cache config values to avoid repeated lookups (overhead not insignificant here)
_dtype = config.adjoint.gradient_dtype_float
_edge_clip_tol = config.adjoint.edge_clip_tolerance
_sample_fraction = config.adjoint.quadrature_sample_fraction
_gauss_order = config.adjoint.gauss_quadrature_order
theta = get_static(self.sidewall_angle)
z_ref = self.reference_axis_pos
cos_th = np.cos(theta)
cos_th = np.clip(cos_th, 1e-12, 1.0)
tan_th = np.tan(theta)
dprime = -tan_th # dd/dz
# axis unit vector in 3D
axis_vec = np.zeros(3, dtype=_dtype)
axis_vec[self.axis] = 1.0
# densify along axis as |theta| grows, dz scales with cos(theta)
z_centers, dz, _z0, _z1 = self._z_slices(sim_min, sim_max, is_2d=is_2d, dx=dx * cos_th)
# early exit: no slices
if (not is_2d) and len(z_centers) == 0:
return {
"centers": np.empty((0, 3), dtype=_dtype),
"normals": np.empty((0, 3), dtype=_dtype),
"perps1": np.empty((0, 3), dtype=_dtype),
"perps2": np.empty((0, 3), dtype=_dtype),
"Ls": np.empty((0,), dtype=_dtype),
"s_vals": np.empty((0,), dtype=_dtype),
"s_weights": np.empty((0,), dtype=_dtype),
"zc_vals": np.empty((0,), dtype=_dtype),
"dz": dz,
"edge_indices": np.empty((0,), dtype=int),
}
# estimate patches for pre-allocation
n_edges = len(vertices)
estimated_patches = 0
denom_edge = max(dx * cos_th, 1e-12)
for ei in range(n_edges):
v0, v1 = vertices[ei], next_v[ei]
L = np.linalg.norm(v1 - v0)
if not np.isclose(L, 0.0):
# prealloc guided by actual step; ds_phys scales with cos(theta)
n_samples = max(1, int(np.ceil(L / denom_edge) * 0.6))
estimated_patches += n_samples * max(1, len(z_centers))
estimated_patches = int(max(1, estimated_patches) * 1.2)
# pre-allocate arrays
centers = np.empty((estimated_patches, 3), dtype=_dtype)
normals = np.empty((estimated_patches, 3), dtype=_dtype)
perps1 = np.empty((estimated_patches, 3), dtype=_dtype)
perps2 = np.empty((estimated_patches, 3), dtype=_dtype)
Ls = np.empty((estimated_patches,), dtype=_dtype)
s_vals = np.empty((estimated_patches,), dtype=_dtype)
s_weights = np.empty((estimated_patches,), dtype=_dtype)
zc_vals = np.empty((estimated_patches,), dtype=_dtype)
edge_indices = np.empty((estimated_patches,), dtype=int)
patch_idx = 0
# if the simulation is effectively 2D (one tangential dimension collapsed),
# slightly expand degenerate bounds to enable finite-length clipping of edges.
sim_min_eff = np.array(sim_min, dtype=_dtype)
sim_max_eff = np.array(sim_max, dtype=_dtype)
for dim in range(3):
if dim == self.axis:
continue
if np.isclose(sim_max_eff[dim] - sim_min_eff[dim], 0.0):
sim_min_eff[dim] -= 0.5 * dx
sim_max_eff[dim] += 0.5 * dx
# pre-compute values that are constant across z slices
n_z = len(z_centers)
z_centers_arr = np.asarray(z_centers, dtype=_dtype)
d_all = -(z_centers_arr - z_ref) * tan_th # (n_z,)
# _shift_vertices() treats positive distances as outward offsets for CCW polygons.
# Preserve the caller's vertex order, but flip the offset for CW inputs so the
# physical sidewall slices match the outward normals used by the VJP integrand.
offset_sign = 1.0 if PolySlab._area(get_static(vertices)) > 0 else -1.0
slice_vertices = np.asarray(
[self._shift_vertices(vertices, offset_sign * dist)[0] for dist in d_all],
dtype=_dtype,
)
slice_next_vertices = np.roll(slice_vertices, -1, axis=1)
# slanted local basis (constant across z for non-slanted case)
# for slanted: rz = axis_vec + dprime * n2d, but dprime is constant
for ei, (v0, v1) in enumerate(zip(vertices, next_v)):
edge_vec = v1 - v0
L = np.sqrt(np.dot(edge_vec, edge_vec))
if L < 1e-12:
continue
# constant along edge: unit tangent in 3D (no axis component)
t_edge = basis["perp1"][ei]
# outward in-plane normal from canonical basis normal
n2d = basis["norm"][ei].copy()
n2d[self.axis] = 0.0
nrm = np.linalg.norm(n2d)
if not np.isclose(nrm, 0.0):
n2d = n2d / nrm
else:
# fallback to right-handed construction if degenerate
tmp = np.cross(axis_vec, t_edge)
n2d = tmp / (np.linalg.norm(tmp) + 1e-20)
# compute basis vectors once per edge
rz = axis_vec + dprime * n2d
T1_vec = t_edge
N_vec = np.cross(T1_vec, rz)
N_norm = np.linalg.norm(N_vec)
if not np.isclose(N_norm, 0.0):
N_vec = N_vec / N_norm
# align N with outward edge normal
if float(np.dot(N_vec, basis["norm"][ei])) < 0.0:
N_vec = -N_vec
T2_vec = np.cross(N_vec, T1_vec)
T2_norm = np.linalg.norm(T2_vec)
if not np.isclose(T2_norm, 0.0):
T2_vec = T2_vec / T2_norm
# batch compute true offset-polygon segment starts and ends for all z slices
# rather than translating the reference edge, which over-samples near vertices
# when the cross-section shrinks.
segment_vertices = slice_vertices[:, ei, :]
segment_next_vertices = slice_next_vertices[:, ei, :]
segment_edges = segment_next_vertices - segment_vertices
segment_lengths = np.linalg.norm(segment_edges, axis=1)
segment_starts = np.empty((n_z, 3), dtype=_dtype)
segment_ends = np.empty((n_z, 3), dtype=_dtype)
plane_axes = [i for i in range(3) if i != self.axis]
segment_starts[:, self.axis] = z_centers_arr
segment_starts[:, plane_axes] = segment_vertices
segment_ends[:, self.axis] = z_centers_arr
segment_ends[:, plane_axes] = segment_next_vertices
# batch clip all z slices at once
is_within_bounds, t_starts, t_ends = self._clip_edges_to_bounds_batch(
segment_starts,
segment_ends,
sim_min_eff,
sim_max_eff,
_edge_clip_tol=_edge_clip_tol,
_dtype=_dtype,
)
# process only valid z slices (sampling has variable output sizes)
valid_indices = np.nonzero(is_within_bounds)[0]
if len(valid_indices) == 0:
continue
# group z slices by unique (t0, t1) pairs to avoid redundant quadrature calculations.
# since most z-slices will have identical clipping bounds (0.0, 1.0),
# we can compute the Gauss samples once and reuse them for almost all slices.
# rounding ensures we get cache hits despite tiny floating point differences.
t0_valid = np.round(t_starts[valid_indices], 10)
t1_valid = np.round(t_ends[valid_indices], 10)
# simple cache for sampling results: (t0, t1) -> (s_list, w_list)
sample_cache = {}
# process each z slice
for zi, t0, t1 in zip(valid_indices, t0_valid, t1_valid):
L_slice = segment_lengths[zi]
if L_slice < 1e-12:
continue
sample_cache_key = (t0, t1, np.round(L_slice, 10))
if sample_cache_key not in sample_cache:
sample_cache[sample_cache_key] = self._adaptive_edge_samples(
L_slice,
denom_edge,
t0,
t1,
_sample_fraction=_sample_fraction,
_gauss_order=_gauss_order,
_dtype=_dtype,
)
s_list, w_list = sample_cache[sample_cache_key]
if len(s_list) == 0:
continue
zc = z_centers_arr[zi]
segment_v0 = segment_vertices[zi]
segment_edge = segment_edges[zi]
pts2d = segment_v0 + s_list[:, None] * segment_edge # faster than np.outer
# inline unpop_axis_vect for xyz computation
n_pts = len(s_list)
xyz = np.empty((n_pts, 3), dtype=_dtype)
xyz[:, self.axis] = zc
xyz[:, plane_axes] = pts2d
n_patches = n_pts
new_size_needed = patch_idx + n_patches
if new_size_needed > centers.shape[0]:
# grow arrays by 1.5x to avoid frequent reallocations
new_size = int(new_size_needed * 1.5)
centers.resize((new_size, 3), refcheck=False)
normals.resize((new_size, 3), refcheck=False)
perps1.resize((new_size, 3), refcheck=False)
perps2.resize((new_size, 3), refcheck=False)
Ls.resize((new_size,), refcheck=False)
s_vals.resize((new_size,), refcheck=False)
s_weights.resize((new_size,), refcheck=False)
zc_vals.resize((new_size,), refcheck=False)
edge_indices.resize((new_size,), refcheck=False)
sl = slice(patch_idx, patch_idx + n_patches)
centers[sl] = xyz
normals[sl] = N_vec
perps1[sl] = T1_vec
perps2[sl] = T2_vec
Ls[sl] = L_slice
s_vals[sl] = s_list
s_weights[sl] = w_list
zc_vals[sl] = zc
edge_indices[sl] = ei
patch_idx += n_patches
# trim arrays to final size
centers = centers[:patch_idx]
normals = normals[:patch_idx]
perps1 = perps1[:patch_idx]
perps2 = perps2[:patch_idx]
Ls = Ls[:patch_idx]
s_vals = s_vals[:patch_idx]
s_weights = s_weights[:patch_idx]
zc_vals = zc_vals[:patch_idx]
edge_indices = edge_indices[:patch_idx]
return {
"centers": centers,
"normals": normals,
"perps1": perps1,
"perps2": perps2,
"Ls": Ls,
"s_vals": s_vals,
"s_weights": s_weights,
"zc_vals": zc_vals,
"dz": dz,
"edge_indices": edge_indices,
}
def _compute_derivative_sidewall_angle(
self,
derivative_info: DerivativeInfo,
sim_min: NDArray,
sim_max: NDArray,
is_2d: bool = False,
interpolators: dict | None = None,
) -> float:
"""VJP for dJ/dtheta where theta = sidewall_angle.
Use dJ/dtheta = integral_S g(x) * V_n(x; theta) * dA, with g(x) from
`evaluate_gradient_at_points`. For a ruled sidewall built by
offsetting the mid-plane polygon by d(z) = -(z - z_ref) * tan(theta),
the normal velocity is V_n = (dd/dtheta) * cos(theta) = -(z - z_ref)/cos(theta)
and the area element is dA = (dz/cos(theta)) * d_ell.
Therefore each patch weight is w = L * dz * (-(z - z_ref)) / cos(theta)^2.
"""
if interpolators is None:
interpolators = derivative_info.create_interpolators(
dtype=config.adjoint.gradient_dtype_float
)
# 2D sim => no dependence on theta (z_local=0)
if is_2d:
return 0.0
# The scalar sidewall-angle VJP should not depend on the user's input winding.
# Use the canonical reference polygon; vertex VJPs keep the raw ordering.
vertices, next_v, edges, basis = self._edge_geometry_arrays(
vertices=self.reference_polygon,
is_ccw=True,
)
dx = derivative_info.adaptive_vjp_spacing()
# collect patches once
patch = self._collect_sidewall_patches(
vertices=vertices,
next_v=next_v,
edges=edges,
basis=basis,
sim_min=sim_min,
sim_max=sim_max,
is_2d=False,
dx=dx,
)
if patch["centers"].shape[0] == 0:
return 0.0
# Shape-derivative factors:
# - Offset: d(z) = -(z - z_ref) * tan(theta)
# - Tangential rate: dd/dtheta = -(z - z_ref) * sec(theta)^2
# - Normal velocity (project to surface normal): V_n = (dd/dtheta) * cos(theta) = -(z - z_ref)/cos(theta)
# - Area element of slanted strip: dA = (dz/cos(theta)) * d_ell
# => Patch weight scales as: V_n * dA = -(z - z_ref) * dz * d_ell / cos(theta)^2
cos_theta = np.cos(get_static(self.sidewall_angle))
inv_cos2 = 1.0 / (cos_theta * cos_theta)
z_ref = self.reference_axis_pos
g = derivative_info.evaluate_gradient_at_points(
patch["centers"], patch["normals"], patch["perps1"], patch["perps2"], interpolators
)
z_local = patch["zc_vals"] - z_ref
weights = patch["Ls"] * patch["s_weights"] * patch["dz"] * (-z_local) * inv_cos2
return float(np.real(np.sum(g * weights)))
def _compute_derivative_slab_bounds(
self, derivative_info: DerivativeInfo, min_max_index: int, interpolators: dict
) -> TracedArrayFloat2D:
"""VJP for one of the two horizontal faces of a ``PolySlab``.
The face is discretized into a Cartesian grid of small planar patches.
The adjoint surface integral is evaluated on every retained patch; the
resulting derivative is split equally between the two vertices that bound
the edge segment.
"""
# rmin/rmax over the geometry and simulation box
if np.isclose(self.slab_bounds[1] - self.slab_bounds[0], 0.0):
log.warning(
"Computing slab face derivatives for flat structures is not fully supported and "
"may give zero for the derivative. Try using a structure with a small, but nonzero "
"thickness for slab bound derivatives."
)
rmin, rmax = derivative_info.bounds_intersect
_, (r1_min, r2_min) = self.pop_axis(rmin, axis=self.axis)
_, (r1_max, r2_max) = self.pop_axis(rmax, axis=self.axis)
ax_val = self.slab_bounds[min_max_index]
# planar grid resolution, clipped to polygon bounding box
face_verts = self.base_polygon if min_max_index == 0 else self.top_polygon
face_poly = shapely.Polygon(face_verts).buffer(fp_eps)
# limit the patch grid to the face that lives inside the simulation box
poly_min_r1, poly_min_r2, poly_max_r1, poly_max_r2 = face_poly.bounds
r1_min = max(r1_min, poly_min_r1)
r1_max = min(r1_max, poly_max_r1)
r2_min = max(r2_min, poly_min_r2)
r2_max = min(r2_max, poly_max_r2)
# intersect the polygon with the simulation bounds
face_poly = face_poly.intersection(shapely.box(r1_min, r2_min, r1_max, r2_max))
if (r1_max <= r1_min) and (r2_max <= r2_min):
# the polygon does not intersect the current simulation slice
return 0.0
# re-compute the extents after clipping to the polygon bounds
extents = np.array([r1_max - r1_min, r2_max - r2_min])
# choose surface or line integral
integral_fun = (
self.compute_derivative_slab_bounds_line
if np.isclose(extents, 0).any()
else self.compute_derivative_slab_bounds_surface
)
return integral_fun(
derivative_info,
extents,
r1_min,
r1_max,
r2_min,
r2_max,
ax_val,
face_poly,
min_max_index,
interpolators,
)
[docs]
def compute_derivative_slab_bounds_line(
self,
derivative_info: DerivativeInfo,
extents: NDArray,
r1_min: float,
r1_max: float,
r2_min: float,
r2_max: float,
ax_val: float,
face_poly: shapely.Polygon,
min_max_index: int,
interpolators: dict,
) -> float:
"""Handle degenerate line cross-section case"""
line_dim = 1 if np.isclose(extents[0], 0) else 0
poly_min_r1, poly_min_r2, poly_max_r1, poly_max_r2 = face_poly.bounds
if line_dim == 0: # x varies, y is fixed
l_min = max(r1_min, poly_min_r1)
l_max = min(r1_max, poly_max_r1)
else: # y varies, x is fixed
l_min = max(r2_min, poly_min_r2)
l_max = min(r2_max, poly_max_r2)
length = l_max - l_min
if np.isclose(length, 0):
return 0.0
dx = derivative_info.adaptive_vjp_spacing()
n_seg = max(1, int(np.ceil(length / dx)))
coords = np.linspace(
l_min, l_max, 2 * n_seg + 1, dtype=config.adjoint.gradient_dtype_float
)[1::2]
# build XY coordinates and in-plane direction vectors
if line_dim == 0:
xy = np.column_stack((coords, np.full_like(coords, r2_min)))
dir_vec_plane = np.column_stack((np.ones_like(coords), np.zeros_like(coords)))
else:
xy = np.column_stack((np.full_like(coords, r1_min), coords))
dir_vec_plane = np.column_stack((np.zeros_like(coords), np.ones_like(coords)))
inside = shapely.contains_xy(face_poly, xy[:, 0], xy[:, 1])
if not inside.any():
return 0.0
xy = xy[inside]
dir_vec_plane = dir_vec_plane[inside]
n_pts = len(xy)
centers_xyz = self.unpop_axis_vect(np.full(n_pts, ax_val), xy)
areas = np.full(n_pts, length / n_seg) # patch length
normals_xyz = self.unpop_axis_vect(
np.full(
n_pts, -1 if min_max_index == 0 else 1, dtype=config.adjoint.gradient_dtype_float
),
np.zeros_like(xy, dtype=config.adjoint.gradient_dtype_float),
)
perps1_xyz = self.unpop_axis_vect(np.zeros(n_pts), dir_vec_plane)
perps2_xyz = self.unpop_axis_vect(np.zeros(n_pts), np.zeros_like(dir_vec_plane))
vjps = derivative_info.evaluate_gradient_at_points(
centers_xyz, normals_xyz, perps1_xyz, perps2_xyz, interpolators
)
return np.real(np.sum(vjps * areas)).item()
[docs]
def compute_derivative_slab_bounds_surface(
self,
derivative_info: DerivativeInfo,
extents: NDArray,
r1_min: float,
r1_max: float,
r2_min: float,
r2_max: float,
ax_val: float,
face_poly: shapely.Polygon,
min_max_index: int,
interpolators: dict,
) -> float:
"""2d surface integral on a Gauss quadrature grid"""
dx = derivative_info.adaptive_vjp_spacing()
# uniform grid would use n1 x n2 points
n1_uniform, n2_uniform = np.maximum(1, np.ceil(extents / dx).astype(int))
# use ~1/2 Gauss points in each direction for similar accuracy
n1 = max(2, n1_uniform // 2)
n2 = max(2, n2_uniform // 2)
g1, w1 = leggauss(n1)
g2, w2 = leggauss(n2)
coords1 = (0.5 * (r1_max - r1_min) * g1 + 0.5 * (r1_max + r1_min)).astype(
config.adjoint.gradient_dtype_float, copy=False
)
coords2 = (0.5 * (r2_max - r2_min) * g2 + 0.5 * (r2_max + r2_min)).astype(
config.adjoint.gradient_dtype_float, copy=False
)
r1_grid, r2_grid = np.meshgrid(coords1, coords2, indexing="ij")
r1_flat = r1_grid.flatten()
r2_flat = r2_grid.flatten()
pts = np.column_stack((r1_flat, r2_flat))
in_face = shapely.contains_xy(face_poly, pts[:, 0], pts[:, 1])
if not in_face.any():
return 0.0
xyz = self.unpop_axis_vect(
np.full(in_face.sum(), ax_val, dtype=config.adjoint.gradient_dtype_float), pts[in_face]
)
n_patches = xyz.shape[0]
normals_xyz = self.unpop_axis_vect(
np.full(
n_patches,
-1 if min_max_index == 0 else 1,
dtype=config.adjoint.gradient_dtype_float,
),
np.zeros((n_patches, 2), dtype=config.adjoint.gradient_dtype_float),
)
perps1_xyz = self.unpop_axis_vect(
np.zeros(n_patches, dtype=config.adjoint.gradient_dtype_float),
np.column_stack(
(
np.ones(n_patches, dtype=config.adjoint.gradient_dtype_float),
np.zeros(n_patches, dtype=config.adjoint.gradient_dtype_float),
)
),
)
perps2_xyz = self.unpop_axis_vect(
np.zeros(n_patches, dtype=config.adjoint.gradient_dtype_float),
np.column_stack(
(
np.zeros(n_patches, dtype=config.adjoint.gradient_dtype_float),
np.ones(n_patches, dtype=config.adjoint.gradient_dtype_float),
)
),
)
w1_grid, w2_grid = np.meshgrid(w1, w2, indexing="ij")
weights_flat = (w1_grid * w2_grid).flatten()[in_face]
jacobian = 0.25 * (r1_max - r1_min) * (r2_max - r2_min)
# area-based correction for non-rectangular domains (e.g. concave polygon)
# for constant integrand, integral should equal polygon area
sum_weights = np.sum(weights_flat)
if sum_weights > 0:
area_correction = face_poly.area / (sum_weights * jacobian)
weights_flat = weights_flat * area_correction
vjps = derivative_info.evaluate_gradient_at_points(
xyz, normals_xyz, perps1_xyz, perps2_xyz, interpolators
)
return np.real(np.sum(vjps * weights_flat * jacobian)).item()
def _compute_derivative_vertices(
self,
derivative_info: DerivativeInfo,
sim_min: NDArray,
sim_max: NDArray,
is_2d: bool = False,
interpolators: dict | None = None,
) -> NDArray:
"""VJP for the vertices of a ``PolySlab``.
Uses shared sidewall patch collection and batched field evaluation.
"""
vertices, next_v, edges, basis = self._edge_geometry_arrays()
dx = derivative_info.adaptive_vjp_spacing()
# collect patches once
patch = self._collect_sidewall_patches(
vertices=vertices,
next_v=next_v,
edges=edges,
basis=basis,
sim_min=sim_min,
sim_max=sim_max,
is_2d=is_2d,
dx=dx,
)
# early return if no patches
if patch["centers"].shape[0] == 0:
return np.zeros_like(vertices)
dz = patch["dz"]
dz_surf = 1.0 if is_2d else dz / np.cos(self.sidewall_angle)
# use provided interpolators or create them if not provided
if interpolators is None:
interpolators = derivative_info.create_interpolators(
dtype=config.adjoint.gradient_dtype_float
)
# evaluate integrand
g = derivative_info.evaluate_gradient_at_points(
patch["centers"], patch["normals"], patch["perps1"], patch["perps2"], interpolators
)
# compute area-based weights and weighted vjps
areas = patch["Ls"] * patch["s_weights"] * dz_surf
patch_vjps = (g * areas).real
# distribute to vertices using vectorized accumulation
normals_2d = np.delete(basis["norm"], self.axis, axis=1)
edge_idx = patch["edge_indices"]
s = patch["s_vals"]
w0 = (1.0 - s) * patch_vjps
w1 = s * patch_vjps
edge_norms = normals_2d[edge_idx]
# Accumulate per-vertex contributions using bincount (O(N_patches))
num_vertices = vertices.shape[0]
contrib0 = w0[:, None] * edge_norms # (n_patches, 2)
contrib1 = w1[:, None] * edge_norms # (n_patches, 2)
idx0 = edge_idx
idx1 = (edge_idx + 1) % num_vertices
v0x = np.bincount(idx0, weights=contrib0[:, 0], minlength=num_vertices)
v0y = np.bincount(idx0, weights=contrib0[:, 1], minlength=num_vertices)
v1x = np.bincount(idx1, weights=contrib1[:, 0], minlength=num_vertices)
v1y = np.bincount(idx1, weights=contrib1[:, 1], minlength=num_vertices)
vjp_per_vertex = np.stack((v0x + v1x, v0y + v1y), axis=1)
return vjp_per_vertex
def _edge_geometry_arrays(
self,
dtype: np.dtype = config.adjoint.gradient_dtype_float,
vertices: NDArray | None = None,
is_ccw: bool | None = None,
) -> tuple[NDArray, NDArray, NDArray, dict[str, NDArray]]:
"""Return (vertices, next_v, edges, basis) arrays for sidewall edge geometry."""
if vertices is None:
vertices = self.vertices
vertices = np.asarray(vertices, dtype=dtype)
next_v = np.roll(vertices, -1, axis=0)
edges = next_v - vertices
basis = self.edge_basis_vectors(edges, is_ccw=is_ccw)
return vertices, next_v, edges, basis
[docs]
def edge_basis_vectors(
self,
edges: NDArray, # (N, 2)
is_ccw: bool | None = None,
) -> dict[str, NDArray]: # (N, 3)
"""Normalized basis vectors for ``normal`` direction, ``slab`` tangent direction and ``edge``."""
# ensure edges have consistent dtype
edges = edges.astype(config.adjoint.gradient_dtype_float, copy=False)
num_vertices, _ = edges.shape
zeros = np.zeros(num_vertices, dtype=config.adjoint.gradient_dtype_float)
ones = np.ones(num_vertices, dtype=config.adjoint.gradient_dtype_float)
# normalized vectors along edges
edges_norm_in_plane = self.normalize_vect(edges)
edges_norm_xyz = self.unpop_axis_vect(zeros, edges_norm_in_plane)
# normalized vectors from base of edges to tops of edges
cos_angle = np.cos(self.sidewall_angle)
sin_angle = np.sin(self.sidewall_angle)
slabs_axis_components = cos_angle * ones
# create axis_norm as array directly to avoid tuple->array conversion in np.cross
axis_norm = np.zeros(3, dtype=config.adjoint.gradient_dtype_float)
axis_norm[self.axis] = 1.0
slab_normal_xyz = -sin_angle * np.cross(edges_norm_xyz, axis_norm)
_, slab_normal_in_plane = self.pop_axis_vect(slab_normal_xyz)
slabs_norm_xyz = self.unpop_axis_vect(slabs_axis_components, slab_normal_in_plane)
# normalized vectors pointing in normal direction of edge
# cross yields inward normal when the extrusion axis is y, so negate once for axis==1
if is_ccw is None:
is_ccw = self.is_ccw
sign = (-1 if self.axis == 1 else 1) * (-1 if not is_ccw else 1)
normals_norm_xyz = sign * np.cross(edges_norm_xyz, slabs_norm_xyz)
return {
"norm": normals_norm_xyz,
"perp1": edges_norm_xyz,
"perp2": slabs_norm_xyz,
}
[docs]
def unpop_axis_vect(self, ax_coords: NDArray, plane_coords: NDArray) -> NDArray:
"""Combine coordinate along axis with coordinates on the plane tangent to the axis.
ax_coords.shape == [N]
plane_coords.shape == [N, 2]
return shape == [N, 3]
"""
n_pts = ax_coords.shape[0]
arr_xyz = np.zeros((n_pts, 3), dtype=ax_coords.dtype)
plane_axes = [i for i in range(3) if i != self.axis]
arr_xyz[:, self.axis] = ax_coords
arr_xyz[:, plane_axes] = plane_coords
return arr_xyz
[docs]
def pop_axis_vect(self, coord: NDArray) -> tuple[NDArray, tuple[NDArray, NDArray]]:
"""Combine coordinate along axis with coordinates on the plane tangent to the axis.
coord.shape == [N, 3]
return shape == ([N], [N, 2]
"""
arr_axis, arrs_plane = self.pop_axis(coord.T, axis=self.axis)
arrs_plane = np.array(arrs_plane).T
return arr_axis, arrs_plane
[docs]
@staticmethod
def normalize_vect(arr: NDArray) -> NDArray:
"""normalize an array shaped (N, d) along the `d` axis and return (N, 1)."""
norm = np.linalg.norm(arr, axis=-1, keepdims=True)
norm = np.where(norm == 0, 1, norm)
return arr / norm
[docs]
def translated(self, x: float, y: float, z: float) -> PolySlab:
"""Return a translated copy of this geometry.
Parameters
----------
x : float
Translation along x.
y : float
Translation along y.
z : float
Translation along z.
Returns
-------
:class:`PolySlab`
Translated copy of this ``PolySlab``.
"""
t_normal, t_plane = self.pop_axis((x, y, z), axis=self.axis)
translated_vertices = np.array(self.vertices) + np.array(t_plane)[None, :]
translated_slab_bounds = (self.slab_bounds[0] + t_normal, self.slab_bounds[1] + t_normal)
return self.updated_copy(vertices=translated_vertices, slab_bounds=translated_slab_bounds)
[docs]
def scaled(self, x: float = 1.0, y: float = 1.0, z: float = 1.0) -> PolySlab:
"""Return a scaled copy of this geometry.
Parameters
----------
x : float = 1.0
Scaling factor along x.
y : float = 1.0
Scaling factor along y.
z : float = 1.0
Scaling factor along z.
Returns
-------
:class:`~tidy3d.Geometry`
Scaled copy of this geometry.
"""
scale_normal, scale_in_plane = self.pop_axis((x, y, z), axis=self.axis)
scaled_vertices = self.vertices * np.array(scale_in_plane)
scaled_slab_bounds = tuple(scale_normal * bound for bound in self.slab_bounds)
return self.updated_copy(vertices=scaled_vertices, slab_bounds=scaled_slab_bounds)
[docs]
def rotated(self, angle: float, axis: Axis | Coordinate) -> PolySlab:
"""Return a rotated copy of this geometry.
Parameters
----------
angle : float
Rotation angle (in radians).
axis : Union[int, tuple[float, float, float]]
Axis of rotation: 0, 1, or 2 for x, y, and z, respectively, or a 3D vector.
Returns
-------
:class:`PolySlab`
Rotated copy of this ``PolySlab``.
"""
_, plane_axs = self.pop_axis([0, 1, 2], self.axis)
if (isinstance(axis, int) and axis == self.axis) or (
isinstance(axis, tuple) and all(axis[ax] == 0 for ax in plane_axs)
):
verts_3d = np.zeros((3, self.vertices.shape[0]))
verts_3d[plane_axs[0], :] = self.vertices[:, 0]
verts_3d[plane_axs[1], :] = self.vertices[:, 1]
rotation = RotationAroundAxis(angle=angle, axis=axis)
rotated_vertices = rotation.rotate_vector(verts_3d)
rotated_vertices = rotated_vertices[plane_axs, :].T
return self.updated_copy(vertices=rotated_vertices)
return super().rotated(angle=angle, axis=axis)
[docs]
def reflected(self, normal: Coordinate) -> PolySlab:
"""Return a reflected copy of this geometry.
Parameters
----------
normal : tuple[float, float, float]
The 3D normal vector of the plane of reflection. The plane is assumed
to pass through the origin (0,0,0).
Returns
-------
:class:`PolySlab`
Reflected copy of this ``PolySlab``.
"""
if math.isclose(normal[self.axis], 0):
_, plane_axs = self.pop_axis((0, 1, 2), self.axis)
verts_3d = np.zeros((3, self.vertices.shape[0]))
verts_3d[plane_axs[0], :] = self.vertices[:, 0]
verts_3d[plane_axs[1], :] = self.vertices[:, 1]
reflection = ReflectionFromPlane(normal=normal)
reflected_vertices = reflection.reflect_vector(verts_3d)
reflected_vertices = reflected_vertices[plane_axs, :].T
return self.updated_copy(vertices=reflected_vertices)
return super().reflected(normal=normal)
[docs]
class ComplexPolySlabBase(PolySlab):
"""Interface for dividing a complex polyslab where self-intersecting polygon can
occur during extrusion. This class should not be used directly. Use instead
:class:`plugins.polyslab.ComplexPolySlab`."""
[docs]
@model_validator(mode="after")
def no_self_intersecting_polygon_during_extrusion(self: Self) -> Self:
"""Turn off the validation for this class."""
return self
[docs]
@classmethod
def from_gds(
cls,
gds_cell: Cell,
axis: Axis,
slab_bounds: tuple[float, float],
gds_layer: int,
gds_dtype: int | None = None,
gds_scale: PositiveFloat = 1.0,
dilation: float = 0.0,
sidewall_angle: float = 0,
reference_plane: PlanePosition = "middle",
) -> list[PolySlab]:
"""Import :class:`.PolySlab` from a ``gdstk.Cell``.
Parameters
----------
gds_cell : gdstk.Cell
``gdstk.Cell`` containing 2D geometric data.
axis : int
Integer index into the polygon's slab axis. (0,1,2) -> (x,y,z).
slab_bounds: tuple[float, float]
Minimum and maximum positions of the slab along ``axis``.
gds_layer : int
Layer index in the ``gds_cell``.
gds_dtype : int = None
Data-type index in the ``gds_cell``.
If ``None``, imports all data for this layer into the returned list.
gds_scale : float = 1.0
Length scale used in GDS file in units of MICROMETER.
For example, if gds file uses nanometers, set ``gds_scale=1e-3``.
Must be positive.
dilation : float = 0.0
Dilation of the polygon in the base by shifting each edge along its
normal outwards direction by a distance;
a negative value corresponds to erosion.
sidewall_angle : float = 0
Angle of the sidewall.
``sidewall_angle=0`` (default) specifies vertical wall,
while ``0<sidewall_angle<np.pi/2`` for the base to be larger than the top.
reference_plane : PlanePosition = "middle"
The position of the GDS layer. It can be at the ``bottom``, ``middle``,
or ``top`` of the PolySlab. E.g. if ``axis=1``, ``bottom`` refers to the
negative side of y-axis, and ``top`` refers to the positive side of y-axis.
Returns
-------
list[:class:`.PolySlab`]
List of :class:`.PolySlab` objects sharing ``axis`` and slab bound properties.
"""
all_vertices = PolySlab._load_gds_vertices(gds_cell, gds_layer, gds_dtype, gds_scale)
polyslabs = [
cls(
vertices=verts,
axis=axis,
slab_bounds=slab_bounds,
dilation=dilation,
sidewall_angle=sidewall_angle,
reference_plane=reference_plane,
)
for verts in all_vertices
]
return [sub_poly for sub_polys in polyslabs for sub_poly in sub_polys.sub_polyslabs]
@property
def geometry_group(self) -> base.GeometryGroup:
"""Divide a complex polyslab into a list of simple polyslabs, which
are assembled into a :class:`.GeometryGroup`.
Returns
-------
:class:`.GeometryGroup`
GeometryGroup for a list of simple polyslabs divided from the complex
polyslab.
"""
return base.GeometryGroup(geometries=self.sub_polyslabs)
@property
def sub_polyslabs(self) -> list[PolySlab]:
"""Divide a complex polyslab into a list of simple polyslabs.
Only neighboring vertex-vertex crossing events are treated in this
version.
Returns
-------
list[PolySlab]
A list of simple polyslabs.
"""
sub_polyslab_list = []
num_division_count = 0
# initialize sub-polyslab parameters
sub_polyslab_dict = self.model_dump(exclude={TYPE_TAG_STR}).copy()
if math.isclose(self.sidewall_angle, 0):
return [PolySlab.model_validate(sub_polyslab_dict)]
sub_polyslab_dict.update({"dilation": 0}) # dilation accounted in setup
# initialize offset distance
offset_distance = 0
for dist_val in self._dilation_length:
dist_now = 0.0
vertices_now = self.reference_polygon
# constructing sub-polyslabs until reaching the base/top
while not math.isclose(dist_now, dist_val):
# bounds for sub-polyslabs assuming no self-intersection
slab_bounds = [
self._dilation_value_at_reference_to_coord(dist_now),
self._dilation_value_at_reference_to_coord(dist_val),
]
# 1) find out any vertices touching events between the current
# position to the base/top
max_dist = PolySlab._neighbor_vertices_crossing_detection(
vertices_now, dist_val - dist_now
)
# vertices touching events captured, update bounds for sub-polyslab
if max_dist is not None:
# max_dist doesn't have sign, so construct signed offset distance
offset_distance = max_dist * dist_val / abs(dist_val)
slab_bounds[1] = self._dilation_value_at_reference_to_coord(
dist_now + offset_distance
)
# 2) construct sub-polyslab
slab_bounds.sort() # for reference_plane=top/bottom, bounds need to be ordered
# direction of marching
reference_plane = "bottom" if dist_val / self._tanq < 0 else "top"
sub_polyslab_dict.update(
{
"slab_bounds": tuple(slab_bounds),
"vertices": vertices_now,
"reference_plane": reference_plane,
}
)
sub_polyslab_list.append(PolySlab.model_validate(sub_polyslab_dict))
# Now Step 3
if max_dist is None:
break
dist_now += offset_distance
# new polygon vertices where collapsing vertices are removed but keep one
vertices_now = PolySlab._shift_vertices(vertices_now, offset_distance)[0]
vertices_now = PolySlab._remove_duplicate_vertices(vertices_now)
# all vertices collapse
if len(vertices_now) < 3:
break
# polygon collapse into 1D
if self.make_shapely_polygon(vertices_now).buffer(0).area < fp_eps:
break
vertices_now = PolySlab._orient(vertices_now)
num_division_count += 1
if num_division_count > _COMPLEX_POLYSLAB_DIVISIONS_WARN:
log.warning(
f"Too many self-intersecting events: the polyslab has been divided into "
f"{num_division_count} polyslabs; more than {_COMPLEX_POLYSLAB_DIVISIONS_WARN} may "
f"slow down the simulation."
)
return sub_polyslab_list
@property
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.finite_length_axis)]
# reverse the dilation value if the reference plane is on the top
if self.reference_plane == "top":
dist = [-dist[0]]
# for middle, both directions
elif self.reference_plane == "middle":
dist = [dist[0] / 2, -dist[0] / 2]
return dist
def _dilation_value_at_reference_to_coord(self, dilation: float) -> float:
"""Compute the coordinate based on the dilation value to the reference plane."""
z_coord = -dilation / self._tanq + self.slab_bounds[0]
if self.reference_plane == "middle":
return z_coord + self.finite_length_axis / 2
if self.reference_plane == "top":
return z_coord + self.finite_length_axis
# bottom case
return z_coord
[docs]
def intersections_tilted_plane(
self,
normal: Coordinate,
origin: Coordinate,
to_2D: MatrixReal4x4,
cleanup: bool = True,
quad_segs: int | None = None,
section_tolerance_2d: bool = False,
) -> list[Shapely]:
"""Return a list of shapely geometries at the plane specified by normal and origin.
Parameters
----------
normal : Coordinate
Vector defining the normal direction to the plane.
origin : Coordinate
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.
section_tolerance_2d : bool = False
See :meth:`tidy3d.components.geometry.base.Geometry.intersections_tilted_plane`.
Returns
-------
list[shapely.geometry.base.BaseGeometry]
List of 2D shapes that intersect plane.
For more details refer to
`Shapely's Documentation <https://shapely.readthedocs.io/en/stable/project.html>`_.
"""
return [
shapely.unary_union(
[
base.Geometry.evaluate_inf_shape(shape)
for polyslab in self.sub_polyslabs
for shape in polyslab.intersections_tilted_plane(
normal,
origin,
to_2D,
cleanup=cleanup,
quad_segs=quad_segs,
section_tolerance_2d=section_tolerance_2d,
)
]
)
]
