"""Abstract base classes for geometry."""
from __future__ import annotations
import functools
import pathlib
from abc import ABC, abstractmethod
from typing import TYPE_CHECKING, Any, ClassVar
import autograd.numpy as np
import shapely
from pydantic import Field, NonNegativeFloat, field_validator, model_validator
from tidy3d.compat import _package_is_older_than
from tidy3d.components.autograd import TracedCoordinate, TracedFloat, TracedSize, get_static
from tidy3d.components.autograd.path_utils import (
AutogradRoute,
format_traced_paths,
indexed_traced_paths,
raise_unsupported_traced_path,
raise_with_traced_path_context,
resolve_delegated_autograd_route,
traced_paths,
validate_traced_path,
)
from tidy3d.components.autograd.types import PathType
from tidy3d.components.base import Tidy3dBaseModel, cached_property
from tidy3d.components.geometry.bound_ops import bounds_intersection, bounds_union
from tidy3d.components.geometry.float_utils import increment_float
from tidy3d.components.transformation import ReflectionFromPlane, RotationAroundAxis
from tidy3d.components.types import (
Axis,
ClipOperationType,
Coordinate,
MatrixReal4x4,
PlanePosition,
)
from tidy3d.components.types.base import discriminated_union
from tidy3d.components.viz import (
ARROW_LENGTH,
PLOT_BUFFER,
add_ax_if_none,
equal_aspect,
plot_params_geometry,
polygon_patch,
set_default_labels_and_title,
)
from tidy3d.constants import LARGE_NUMBER, MICROMETER, RADIAN, fp_eps, inf
from tidy3d.exceptions import (
AdjointError,
SetupError,
Tidy3dError,
Tidy3dImportError,
Tidy3dKeyError,
ValidationError,
format_chained_exception_message,
)
from tidy3d.log import log
from tidy3d.packaging import verify_packages_import
if TYPE_CHECKING:
from collections.abc import Callable, Iterable
from os import PathLike
import pydantic
from gdstk import Cell
from matplotlib.backend_bases import Event
from matplotlib.patches import FancyArrowPatch
from numpy.typing import ArrayLike, NDArray
from pydantic import NonNegativeInt, PositiveFloat
from typing_extensions import Self
from tidy3d.components.autograd import AutogradFieldMap
from tidy3d.components.autograd.derivative_utils import DerivativeInfo
from tidy3d.components.types import (
ArrayFloat2D,
ArrayFloat3D,
Ax,
Bound,
Coordinate2D,
LengthUnit,
Shapely,
Size,
)
from tidy3d.components.viz import PlotParams, VisualizationSpec
POLY_GRID_SIZE = 1e-12
POLY_TOLERANCE_RATIO = 1e-12
POLY_DISTANCE_TOLERANCE = 8e-12
# Tolerance for validating linear-only transforms (no translation)
LINEAR_TRANSFORM_TOL = 1e-12
GDS_MAX_COORDINATE_INDEX = 2**31 - 1
def _raise_unsupported_traced_geometry_path(
geometry_name: str,
field_path: tuple[Any, ...],
*,
supported_parameters: tuple[str, ...] = (),
) -> None:
"""Raise a user-facing validation error for an unsupported geometry trace."""
raise_unsupported_traced_path(
parameter_kind="geometry",
owner_kind="geometry type",
owner_name=geometry_name,
field_path=field_path,
supported_parameters=supported_parameters,
)
_shapely_operations = {
"union": shapely.union,
"intersection": shapely.intersection,
"difference": shapely.difference,
"symmetric_difference": shapely.symmetric_difference,
}
_bit_operations = {
"union": lambda a, b: a | b,
"intersection": lambda a, b: a & b,
"difference": lambda a, b: a & ~b,
"symmetric_difference": lambda a, b: a != b,
}
# Validators for geometry classes (defined here instead of validators.py to avoid circular imports)
def assert_geometry_finite(field_name: str = "geometry") -> Callable[[type, Geometry], Geometry]:
"""Validator that ensures a geometry field has finite bounds."""
@field_validator(field_name)
@classmethod
def geometry_has_finite_bounds(cls: type, val: Geometry) -> Geometry:
"""Raise validation error if geometry has non-finite bounds."""
if not np.isfinite(val.bounds).all():
raise ValidationError(
f"'{cls.__name__}' requires a geometry with finite dimensions. "
"Try using a large value instead of 'inf' when creating geometries."
)
return val
return geometry_has_finite_bounds
def check_transform_invertible(transform: MatrixReal4x4, index: int | None = None) -> None:
"""Check if a transform matrix is invertible.
Parameters
----------
transform : MatrixReal4x4
The 4x4 transformation matrix to check.
index : Optional[int]
If provided, includes the index in the error message (for array of transforms).
Raises
------
ValidationError
If the transform matrix is not invertible.
"""
try:
_ = np.linalg.inv(transform)
except np.linalg.LinAlgError as err:
if index is not None:
raise ValidationError(
format_chained_exception_message(
f"Transform at index {index} is not invertible", err
)
) from err
raise ValidationError(
format_chained_exception_message("Transform matrix is not invertible", err)
) from err
[docs]
class Geometry(Tidy3dBaseModel, ABC):
"""Abstract base class, defines where something exists in space."""
_traced_supported_paths: ClassVar[tuple[PathType, ...]] = ()
@classmethod
def _traced_autograd_supported_parameters(cls) -> tuple[str, ...]:
"""Return user-facing supported parameter names for setup validation."""
return format_traced_paths(cls._traced_supported_paths)
@cached_property
def plot_params(self) -> PlotParams:
"""Default parameters for plotting a Geometry object."""
return plot_params_geometry
[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.
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.
"""
def point_inside(x: float, y: float, z: float) -> bool:
"""Returns ``True`` if a single point ``(x, y, z)`` is inside."""
shapes_intersect = self.intersections_plane(z=z)
loc = self.make_shapely_point(x, y)
return any(shape.contains(loc) for shape in shapes_intersect)
arrays = tuple(map(np.array, (x, y, z)))
self._ensure_equal_shape(*arrays)
inside = np.zeros((arrays[0].size,), dtype=bool)
arrays_flat = map(np.ravel, arrays)
for ipt, args in enumerate(zip(*arrays_flat)):
inside[ipt] = point_inside(*args)
return inside.reshape(arrays[0].shape)
@staticmethod
def _ensure_equal_shape(*arrays: Any) -> None:
"""Ensure all input arrays have the same shape."""
shapes = {np.array(arr).shape for arr in arrays}
if len(shapes) > 1:
raise ValueError("All coordinate inputs (x, y, z) must have the same shape.")
[docs]
@staticmethod
def make_shapely_box(minx: float, miny: float, maxx: float, maxy: float) -> shapely.box:
"""Make a shapely box ensuring everything untraced."""
minx = get_static(minx)
miny = get_static(miny)
maxx = get_static(maxx)
maxy = get_static(maxy)
return shapely.box(minx, miny, maxx, maxy)
[docs]
@staticmethod
def make_shapely_point(minx: float, miny: float) -> shapely.Point:
"""Make a shapely Point ensuring everything untraced."""
minx = get_static(minx)
miny = get_static(miny)
return shapely.Point(minx, miny)
def _inds_inside_bounds(
self, x: NDArray[float], y: NDArray[float], z: NDArray[float]
) -> tuple[slice, slice, slice]:
"""Return slices into the sorted input arrays that are inside the geometry bounds.
Parameters
----------
x : np.ndarray[float]
1D array of point positions in x direction.
y : np.ndarray[float]
1D array of point positions in y direction.
z : np.ndarray[float]
1D array of point positions in z direction.
Returns
-------
tuple[slice, slice, slice]
Slices into each of the three arrays that are inside the geometry bounds.
"""
bounds = self.bounds
inds_in = []
for dim, coords in enumerate([x, y, z]):
inds = np.nonzero((bounds[0][dim] <= coords) * (coords <= bounds[1][dim]))[0]
inds_in.append(slice(0, 0) if inds.size == 0 else slice(inds[0], inds[-1] + 1))
return tuple(inds_in)
[docs]
def inside_meshgrid(
self, x: NDArray[float], y: NDArray[float], z: NDArray[float]
) -> NDArray[bool]:
"""Perform ``self.inside`` on a set of sorted 1D coordinates. Applies meshgrid to the
supplied coordinates before checking inside.
Parameters
----------
x : np.ndarray[float]
1D array of point positions in x direction.
y : np.ndarray[float]
1D array of point positions in y direction.
z : np.ndarray[float]
1D array of point positions in z direction.
Returns
-------
np.ndarray[bool]
Array with shape ``(x.size, y.size, z.size)``, which is ``True`` for every
point that is inside the geometry.
"""
arrays = tuple(map(np.array, (x, y, z)))
if any(arr.ndim != 1 for arr in arrays):
raise ValueError("Each of the supplied coordinates (x, y, z) must be 1D.")
shape = tuple(arr.size for arr in arrays)
is_inside = np.zeros(shape, dtype=bool)
inds_inside = self._inds_inside_bounds(*arrays)
coords_inside = tuple(arr[ind] for ind, arr in zip(inds_inside, arrays))
coords_3d = np.meshgrid(*coords_inside, indexing="ij")
is_inside[inds_inside] = self.inside(*coords_3d)
return is_inside
[docs]
@abstractmethod
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. If ``None``, uses
high-quality visualization settings.
section_tolerance_2d : bool = False
If ``True``, allow effectively zero-thickness geometries to contribute a section
when the requested plane is within ``fp_eps`` of the geometry bounds along that
axis. Intended for plotting paths where small transform or snap offsets should not
hide 2D structures; does not affect strictly 3D geometries.
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>`_.
"""
[docs]
def intersections_plane(
self,
x: float | None = None,
y: float | None = None,
z: float | None = None,
cleanup: bool = True,
quad_segs: int | None = None,
section_tolerance_2d: bool = False,
) -> list[Shapely]:
"""Returns list of shapely geometries at plane specified by one non-None value of x,y,z.
Parameters
----------
x : float = None
Position of plane in x direction, only one of x,y,z can be specified to define plane.
y : float = None
Position of plane in y direction, only one of x,y,z can be specified to define plane.
z : float = None
Position of plane in z direction, only one of x,y,z can be specified to define plane.
cleanup : bool = True
If True, removes extremely small features from each polygon's boundary.
quad_segs : Optional[int] = None
Number of segments used to discretize circular shapes. If ``None``, uses
high-quality visualization settings.
section_tolerance_2d : bool = False
If ``True``, allow effectively zero-thickness geometries to contribute a section
when the requested plane is within ``fp_eps`` of the geometry bounds along that
axis. Intended for plotting paths where small transform or snap offsets should not
hide 2D structures; does not affect strictly 3D geometries.
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>`_.
"""
axis, position = self.parse_xyz_kwargs(x=x, y=y, z=z)
origin = self.unpop_axis(position, (0, 0), axis=axis)
normal = self.unpop_axis(1, (0, 0), axis=axis)
to_2D = np.eye(4)
if axis != 2:
last, indices = self.pop_axis((0, 1, 2), axis)
to_2D = to_2D[[*list(indices), last, 3]]
return self.intersections_tilted_plane(
normal,
origin,
to_2D,
cleanup=cleanup,
quad_segs=quad_segs,
section_tolerance_2d=section_tolerance_2d,
)
[docs]
def intersections_2dbox(self, plane: Box) -> list[Shapely]:
"""Returns list of shapely geometries representing the intersections of the geometry with
a 2D box.
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>`_.
"""
log.warning(
"'intersections_2dbox()' is deprecated and will be removed in the future. "
"Use 'plane.intersections_with(...)' for the same functionality."
)
return plane.intersections_with(self)
[docs]
def intersects(
self, other: Geometry, strict_inequality: tuple[bool, bool, bool] = [False, False, False]
) -> bool:
"""Returns ``True`` if two :class:`~tidy3d.Geometry` have intersecting `.bounds`.
Parameters
----------
other : :class:`~tidy3d.Geometry`
Geometry to check intersection with.
strict_inequality : tuple[bool, bool, bool] = [False, False, False]
For each dimension, defines whether to include equality in the boundaries comparison.
If ``False``, equality is included, and two geometries that only intersect at their
boundaries will evaluate as ``True``. If ``True``, such geometries will evaluate as
``False``.
Returns
-------
bool
Whether the rectangular bounding boxes of the two geometries intersect.
"""
self_bmin, self_bmax = self.bounds
other_bmin, other_bmax = other.bounds
for smin, omin, smax, omax, strict in zip(
self_bmin, other_bmin, self_bmax, other_bmax, strict_inequality
):
# are all of other's minimum coordinates less than self's maximum coordinate?
in_minus = omin < smax if strict else omin <= smax
# are all of other's maximum coordinates greater than self's minimum coordinate?
in_plus = omax > smin if strict else omax >= smin
# if either failed, return False
if not all((in_minus, in_plus)):
return False
return True
[docs]
def contains(
self, other: Geometry, strict_inequality: tuple[bool, bool, bool] = [False, False, False]
) -> bool:
"""Returns ``True`` if the `.bounds` of ``other`` are contained within the
`.bounds` of ``self``.
Parameters
----------
other : :class:`~tidy3d.Geometry`
Geometry to check containment with.
strict_inequality : tuple[bool, bool, bool] = [False, False, False]
For each dimension, defines whether to include equality in the boundaries comparison.
If ``False``, equality will be considered as contained. If ``True``, ``other``'s
bounds must be strictly within the bounds of ``self``.
Returns
-------
bool
Whether the rectangular bounding box of ``other`` is contained within the bounding
box of ``self``.
"""
self_bmin, self_bmax = self.bounds
other_bmin, other_bmax = other.bounds
for smin, omin, smax, omax, strict in zip(
self_bmin, other_bmin, self_bmax, other_bmax, strict_inequality
):
# are all of other's minimum coordinates greater than self's minimim coordinate?
in_minus = omin > smin if strict else omin >= smin
# are all of other's maximum coordinates less than self's maximum coordinate?
in_plus = omax < smax if strict else omax <= smax
# if either failed, return False
if not all((in_minus, in_plus)):
return False
return True
[docs]
def intersects_plane(
self, x: float | None = None, y: float | None = None, z: float | None = None
) -> bool:
"""Whether self intersects plane specified by one non-None value of x,y,z.
Parameters
----------
x : float = None
Position of plane in x direction, only one of x,y,z can be specified to define plane.
y : float = None
Position of plane in y direction, only one of x,y,z can be specified to define plane.
z : float = None
Position of plane in z direction, only one of x,y,z can be specified to define plane.
Returns
-------
bool
Whether this geometry intersects the plane.
"""
axis, position = self.parse_xyz_kwargs(x=x, y=y, z=z)
return self.intersects_axis_position(axis, position)
[docs]
def intersects_axis_position(
self, axis: int, position: float, section_tolerance_2d: bool = False
) -> bool:
"""Whether self intersects plane specified by a given position along a normal axis.
Parameters
----------
axis : int = None
Axis normal to the plane.
position : float = None
Position of plane along the normal axis.
section_tolerance_2d : bool = False
If ``True``, allow effectively zero-thickness geometries to intersect a plane when
the plane is within ``fp_eps`` of the geometry bounds along that axis.
Returns
-------
bool
Whether this geometry intersects the plane.
"""
min_bound = self.bounds[0][axis]
max_bound = self.bounds[1][axis]
if min_bound <= position <= max_bound:
return True
if section_tolerance_2d and np.isclose(min_bound, max_bound, rtol=fp_eps, atol=fp_eps):
return bool(
np.isclose(position, min_bound, rtol=fp_eps, atol=fp_eps)
or np.isclose(position, max_bound, rtol=fp_eps, atol=fp_eps)
)
return False
@cached_property
@abstractmethod
def bounds(self) -> Bound:
"""Returns bounding box min and max coordinates.
Returns
-------
tuple[float, float, float], tuple[float, float float]
Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
"""
[docs]
@staticmethod
def bounds_intersection(bounds1: Bound, bounds2: Bound) -> Bound:
"""Return the bounds that are the intersection of two bounds."""
return bounds_intersection(bounds1, bounds2)
[docs]
@staticmethod
def bounds_union(bounds1: Bound, bounds2: Bound) -> Bound:
"""Return the bounds that are the union of two bounds."""
return bounds_union(bounds1, bounds2)
@cached_property
def bounding_box(self) -> Box:
"""Returns :class:`~tidy3d.Box` representation of the bounding box of a :class:`~tidy3d.Geometry`.
Returns
-------
:class:`~tidy3d.Box`
Geometric object representing bounding box.
"""
return Box.from_bounds(*self.bounds)
@cached_property
def zero_dims(self) -> list[Axis]:
"""A list of axes along which the :class:`~tidy3d.Geometry` is zero-sized based on its bounds."""
zero_dims = []
for dim in range(3):
if self.bounds[1][dim] == self.bounds[0][dim]:
zero_dims.append(dim)
return zero_dims
def _pop_bounds(self, axis: Axis) -> tuple[Coordinate2D, tuple[Coordinate2D, Coordinate2D]]:
"""Returns min and max bounds in plane normal to and tangential to ``axis``.
Parameters
----------
axis : int
Integer index into 'xyz' (0,1,2).
Returns
-------
tuple[float, float], tuple[tuple[float, float], tuple[float, float]]
Bounds along axis and a tuple of bounds in the ordered planar coordinates.
Packed as ``(zmin, zmax), ((xmin, ymin), (xmax, ymax))``.
"""
b_min, b_max = self.bounds
zmin, (xmin, ymin) = self.pop_axis(b_min, axis=axis)
zmax, (xmax, ymax) = self.pop_axis(b_max, axis=axis)
return (zmin, zmax), ((xmin, ymin), (xmax, ymax))
@staticmethod
def _get_center(pt_min: float, pt_max: float) -> float:
"""Returns center point based on bounds along dimension."""
if np.isneginf(pt_min) and np.isposinf(pt_max):
return 0.0
if np.isneginf(pt_min) or np.isposinf(pt_max):
raise SetupError(
f"Bounds of ({pt_min}, {pt_max}) supplied along one dimension. "
"We currently don't support a single ``inf`` value in bounds for ``Box``. "
"To construct a semi-infinite ``Box``, "
"please supply a large enough number instead of ``inf``. "
"For example, a location extending outside of the "
"Simulation domain (including PML)."
)
return (pt_min + pt_max) / 2.0
@cached_property
def _normal_2dmaterial(self) -> Axis:
"""Get the normal to the given geometry, checking that it is a 2D geometry."""
raise ValidationError("'Medium2D' is not compatible with this geometry class.")
def _update_from_bounds(self, bounds: tuple[float, float], axis: Axis) -> Geometry:
"""Returns an updated geometry which has been transformed to fit within ``bounds``
along the ``axis`` direction."""
raise NotImplementedError(
"'_update_from_bounds' is not compatible with this geometry class."
)
[docs]
@equal_aspect
@add_ax_if_none
def plot(
self,
x: float | None = None,
y: float | None = None,
z: float | None = None,
ax: Ax = None,
plot_length_units: LengthUnit = None,
viz_spec: VisualizationSpec = None,
**patch_kwargs: Any,
) -> Ax:
"""Plot geometry cross section at single (x,y,z) coordinate.
Parameters
----------
x : float = None
Position of plane in x direction, only one of x,y,z can be specified to define plane.
y : float = None
Position of plane in y direction, only one of x,y,z can be specified to define plane.
z : float = None
Position of plane in z direction, only one of x,y,z can be specified to define plane.
ax : matplotlib.axes._subplots.Axes = None
Matplotlib axes to plot on, if not specified, one is created.
plot_length_units : LengthUnit = None
Specify units to use for axis labels, tick labels, and the title.
viz_spec : VisualizationSpec = None
Plotting parameters associated with a medium to use instead of defaults.
**patch_kwargs
Optional keyword arguments passed to the matplotlib patch plotting of structure.
For details on accepted values, refer to
`Matplotlib's documentation <https://tinyurl.com/2nf5c2fk>`_.
Returns
-------
matplotlib.axes._subplots.Axes
The supplied or created matplotlib axes.
"""
# find shapes that intersect self at plane
axis, _position = self.parse_xyz_kwargs(x=x, y=y, z=z)
shapes_intersect = self.intersections_plane(x=x, y=y, z=z, section_tolerance_2d=True)
plot_params = self.plot_params
if viz_spec is not None:
plot_params = plot_params.override_with_viz_spec(viz_spec)
plot_params = plot_params.include_kwargs(**patch_kwargs)
# for each intersection, plot the shape
for shape in shapes_intersect:
ax = self.plot_shape(shape, plot_params=plot_params, ax=ax)
# clean up the axis display
ax = self.add_ax_lims(axis=axis, ax=ax)
ax.set_aspect("equal")
# Add the default axis labels, tick labels, and title
ax = Box.add_ax_labels_and_title(ax=ax, x=x, y=y, z=z, plot_length_units=plot_length_units)
return ax
[docs]
def plot_shape(self, shape: Shapely, plot_params: PlotParams, ax: Ax) -> Ax:
"""Defines how a shape is plotted on a matplotlib axes."""
if shape.geom_type in (
"MultiPoint",
"MultiLineString",
"MultiPolygon",
"GeometryCollection",
):
for sub_shape in shape.geoms:
ax = self.plot_shape(shape=sub_shape, plot_params=plot_params, ax=ax)
return ax
_shape = Geometry.evaluate_inf_shape(shape)
if _shape.geom_type == "LineString":
xs, ys = zip(*_shape.coords)
ax.plot(xs, ys, color=plot_params.facecolor, linewidth=plot_params.linewidth)
elif _shape.geom_type == "Point":
ax.scatter(shape.x, shape.y, color=plot_params.facecolor)
else:
patch = polygon_patch(_shape, **plot_params.to_kwargs())
ax.add_artist(patch)
return ax
@staticmethod
def _do_not_intersect(
bounds_a: float, bounds_b: float, shape_a: Shapely, shape_b: Shapely
) -> bool:
"""Check whether two shapes intersect."""
# do a bounding box check to see if any intersection to do anything about
if (
bounds_a[0] > bounds_b[2]
or bounds_b[0] > bounds_a[2]
or bounds_a[1] > bounds_b[3]
or bounds_b[1] > bounds_a[3]
):
return True
# look more closely to see if intersected.
if shape_b.is_empty or not shape_a.intersects(shape_b):
return True
return False
@staticmethod
def _get_plot_labels(axis: Axis) -> tuple[str, str]:
"""Returns planar coordinate x and y axis labels for cross section plots.
Parameters
----------
axis : int
Integer index into 'xyz' (0,1,2).
Returns
-------
str, str
Labels of plot, packaged as ``(xlabel, ylabel)``.
"""
_, (xlabel, ylabel) = Geometry.pop_axis("xyz", axis=axis)
return xlabel, ylabel
def _get_plot_limits(
self, axis: Axis, buffer: float = PLOT_BUFFER
) -> tuple[Coordinate2D, Coordinate2D]:
"""Gets planar coordinate limits for cross section plots.
Parameters
----------
axis : int
Integer index into 'xyz' (0,1,2).
buffer : float = 0.3
Amount of space to add around the limits on the + and - sides.
Returns
-------
tuple[float, float], tuple[float, float]
The x and y plot limits, packed as ``(xmin, xmax), (ymin, ymax)``.
"""
_, ((xmin, ymin), (xmax, ymax)) = self._pop_bounds(axis=axis)
return (xmin - buffer, xmax + buffer), (ymin - buffer, ymax + buffer)
[docs]
def add_ax_lims(self, axis: Axis, ax: Ax, buffer: float = PLOT_BUFFER) -> Ax:
"""Sets the x,y limits based on ``self.bounds``.
Parameters
----------
axis : int
Integer index into 'xyz' (0,1,2).
ax : matplotlib.axes._subplots.Axes
Matplotlib axes to add labels and limits on.
buffer : float = 0.3
Amount of space to place around the limits on the + and - sides.
Returns
-------
matplotlib.axes._subplots.Axes
The supplied or created matplotlib axes.
"""
(xmin, xmax), (ymin, ymax) = self._get_plot_limits(axis=axis, buffer=buffer)
# note: axes limits dont like inf values, so we need to evaluate them first if present
xmin, xmax, ymin, ymax = self._evaluate_inf((xmin, xmax, ymin, ymax))
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
return ax
[docs]
@staticmethod
def add_ax_labels_and_title(
ax: Ax,
x: float | None = None,
y: float | None = None,
z: float | None = None,
plot_length_units: LengthUnit = None,
) -> Ax:
"""Sets the axis labels, tick labels, and title based on ``axis``
and an optional ``plot_length_units`` argument.
Parameters
----------
ax : matplotlib.axes._subplots.Axes
Matplotlib axes to add labels and limits on.
x : float = None
Position of plane in x direction, only one of x,y,z can be specified to define plane.
y : float = None
Position of plane in y direction, only one of x,y,z can be specified to define plane.
z : float = None
Position of plane in z direction, only one of x,y,z can be specified to define plane.
plot_length_units : LengthUnit = None
When set to a supported ``LengthUnit``, plots will be produced with annotated axes
and title with the proper units.
Returns
-------
matplotlib.axes._subplots.Axes
The supplied matplotlib axes.
"""
axis, position = Box.parse_xyz_kwargs(x=x, y=y, z=z)
axis_labels = Box._get_plot_labels(axis)
ax = set_default_labels_and_title(
axis_labels=axis_labels,
axis=axis,
position=position,
ax=ax,
plot_length_units=plot_length_units,
)
return ax
@staticmethod
def _evaluate_inf(array: ArrayLike) -> NDArray[np.floating]:
"""Processes values and evaluates any infs into large (signed) numbers."""
array = get_static(np.array(array))
return np.where(np.isinf(array), np.sign(array) * LARGE_NUMBER, array)
[docs]
@staticmethod
def evaluate_inf_shape(shape: Shapely) -> Shapely:
"""Returns a copy of shape with inf vertices replaced by large numbers if polygon."""
if not any(np.isinf(b) for b in shape.bounds):
return shape
return shapely.transform(shape, Geometry._evaluate_inf, include_z=None)
[docs]
@staticmethod
def pop_axis(coord: tuple[Any, Any, Any], axis: int) -> tuple[Any, tuple[Any, Any]]:
"""Separates coordinate at ``axis`` index from coordinates on the plane tangent to ``axis``.
Parameters
----------
coord : tuple[Any, Any, Any]
Tuple of three values in original coordinate system.
axis : int
Integer index into 'xyz' (0,1,2).
Returns
-------
Any, tuple[Any, Any]
The input coordinates are separated into the one along the axis provided
and the two on the planar coordinates,
like ``axis_coord, (planar_coord1, planar_coord2)``.
"""
plane_vals = list(coord)
axis_val = plane_vals.pop(axis)
return axis_val, tuple(plane_vals)
[docs]
@staticmethod
def unpop_axis(ax_coord: Any, plane_coords: tuple[Any, Any], axis: int) -> tuple[Any, Any, Any]:
"""Combine coordinate along axis with coordinates on the plane tangent to the axis.
Parameters
----------
ax_coord : Any
Value along axis direction.
plane_coords : tuple[Any, Any]
Values along ordered planar directions.
axis : int
Integer index into 'xyz' (0,1,2).
Returns
-------
tuple[Any, Any, Any]
The three values in the xyz coordinate system.
"""
coords = list(plane_coords)
coords.insert(axis, ax_coord)
return tuple(coords)
[docs]
@staticmethod
def parse_xyz_kwargs(**xyz: Any) -> tuple[Axis, float]:
"""Turns x,y,z kwargs into index of the normal axis and position along that axis.
Parameters
----------
x : float = None
Position of plane in x direction, only one of x,y,z can be specified to define plane.
y : float = None
Position of plane in y direction, only one of x,y,z can be specified to define plane.
z : float = None
Position of plane in z direction, only one of x,y,z can be specified to define plane.
Returns
-------
int, float
Index into xyz axis (0,1,2) and position along that axis.
"""
xyz_filtered = {k: v for k, v in xyz.items() if v is not None}
if len(xyz_filtered) != 1:
raise ValueError("exactly one kwarg in [x,y,z] must be specified.")
axis_label, position = list(xyz_filtered.items())[0]
axis = "xyz".index(axis_label)
return axis, position
@staticmethod
def _validate_gds_precision(
*,
polygons: list[Any],
gds_precision: float,
context: str,
) -> float:
"""Validate that the requested GDS precision is safe for the written polygons."""
if not np.isfinite(gds_precision) or gds_precision <= 0:
raise SetupError(
f"Requested 'gds_precision={gds_precision:.6g} um' in {context} must be "
"positive and finite."
)
if not polygons:
return gds_precision
max_abs_coord = 0.0
for polygon in polygons:
bbox = polygon.bounding_box()
if bbox is None:
continue
for point in bbox:
for value in point:
coordinate = float(value)
if not np.isfinite(coordinate):
raise SetupError(
f"Cannot export non-finite GDS coordinate '{coordinate}' in "
f"{context}. Use finite geometry bounds before exporting to GDS."
)
max_abs_coord = max(max_abs_coord, abs(coordinate))
if max_abs_coord <= 0:
return gds_precision
min_safe_precision = float(np.nextafter(max_abs_coord / GDS_MAX_COORDINATE_INDEX, np.inf))
if gds_precision >= min_safe_precision:
return gds_precision
raise SetupError(
f"Requested 'gds_precision={gds_precision:.6g} um' in {context} is too fine for "
f"the export bounds (+/-{max_abs_coord:.6g} um). The minimum safe precision is "
f"'{min_safe_precision:.6g} um' to stay within the signed 32-bit GDS coordinate "
"range. Use a larger 'gds_precision'."
)
[docs]
@staticmethod
def parse_two_xyz_kwargs(**xyz: Any) -> list[tuple[Axis, float]]:
"""Turns x,y,z kwargs into indices of axes and the position along each axis.
Parameters
----------
x : float = None
Position in x direction, only two of x,y,z can be specified to define line.
y : float = None
Position in y direction, only two of x,y,z can be specified to define line.
z : float = None
Position in z direction, only two of x,y,z can be specified to define line.
Returns
-------
[(int, float), (int, float)]
Index into xyz axis (0,1,2) and position along that axis.
"""
xyz_filtered = {k: v for k, v in xyz.items() if v is not None}
assert len(xyz_filtered) == 2, "exactly two kwarg in [x,y,z] must be specified."
xyz_list = list(xyz_filtered.items())
return [("xyz".index(axis_label), position) for axis_label, position in xyz_list]
[docs]
@staticmethod
def rotate_points(points: ArrayFloat3D, axis: Coordinate, angle: float) -> ArrayFloat3D:
"""Rotate a set of points in 3D.
Parameters
----------
points : ArrayLike[float]
Array of shape ``(3, ...)``.
axis : Coordinate
Axis of rotation
angle : float
Angle of rotation counter-clockwise around the axis (rad).
"""
rotation = RotationAroundAxis(axis=axis, angle=angle)
return rotation.rotate_vector(points)
[docs]
def reflect_points(
self,
points: ArrayFloat3D,
polar_axis: Axis,
angle_theta: float,
angle_phi: float,
) -> ArrayFloat3D:
"""Reflect a set of points in 3D at a plane passing through the coordinate origin defined
and normal to a given axis defined in polar coordinates (theta, phi) w.r.t. the
``polar_axis`` which can be 0, 1, or 2.
Parameters
----------
points : ArrayLike[float]
Array of shape ``(3, ...)``.
polar_axis : Axis
Cartesian axis w.r.t. which the normal axis angles are defined.
angle_theta : float
Polar angle w.r.t. the polar axis.
angle_phi : float
Azimuth angle around the polar axis.
"""
# Rotate such that the plane normal is along the polar_axis
axis_theta, axis_phi = [0, 0, 0], [0, 0, 0]
axis_phi[polar_axis] = 1
plane_axes = [0, 1, 2]
plane_axes.pop(polar_axis)
axis_theta[plane_axes[1]] = 1
points_new = self.rotate_points(points, axis_phi, -angle_phi)
points_new = self.rotate_points(points_new, axis_theta, -angle_theta)
# Flip the ``polar_axis`` coordinate of the points, which is now normal to the plane
points_new[polar_axis, :] *= -1
# Rotate back
points_new = self.rotate_points(points_new, axis_theta, angle_theta)
points_new = self.rotate_points(points_new, axis_phi, angle_phi)
return points_new
[docs]
def volume(self, bounds: Bound = None) -> float:
"""Returns object's volume with optional bounds.
Parameters
----------
bounds : tuple[tuple[float, float, float], tuple[float, float, float]] = None
Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
Returns
-------
float
Volume in um^3.
"""
if not bounds:
bounds = self.bounds
return self._volume(bounds)
@abstractmethod
def _volume(self, bounds: Bound) -> float:
"""Returns object's volume within given bounds."""
[docs]
def surface_area(self, bounds: Bound = None) -> float:
"""Returns object's surface area with optional bounds.
Parameters
----------
bounds : tuple[tuple[float, float, float], tuple[float, float, float]] = None
Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
Returns
-------
float
Surface area in um^2.
"""
if not bounds:
bounds = self.bounds
return self._surface_area(bounds)
@abstractmethod
def _surface_area(self, bounds: Bound) -> float:
"""Returns object's surface area within given bounds."""
[docs]
def translated(self, x: float, y: float, z: float) -> Geometry:
"""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:`~tidy3d.Geometry`
Translated copy of this geometry.
"""
return Transformed(geometry=self, transform=Transformed.translation(x, y, z))
[docs]
def scaled(self, x: float = 1.0, y: float = 1.0, z: float = 1.0) -> Geometry:
"""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.
"""
return Transformed(geometry=self, transform=Transformed.scaling(x, y, z))
[docs]
def rotated(self, angle: float, axis: Axis | Coordinate) -> Geometry:
"""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:`~tidy3d.Geometry`
Rotated copy of this geometry.
"""
return Transformed(geometry=self, transform=Transformed.rotation(angle, axis))
[docs]
def reflected(self, normal: Coordinate) -> Geometry:
"""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:`~tidy3d.Geometry`
Reflected copy of this geometry.
"""
return Transformed(geometry=self, transform=Transformed.reflection(normal))
[docs]
def array(
self,
offsets: ArrayLike | None = None,
transforms: ArrayLike | None = None,
) -> GeometryArray:
"""Return an array of copies of this geometry with optional offsets and/or linear transforms.
This method creates a :class:`GeometryArray` containing multiple copies of this
geometry. When both ``offsets`` and ``transforms`` are provided, transforms are
applied to each copy before the translation given by offsets is applied.
Parameters
----------
offsets : Optional[ArrayLike] = None
Optional array of offset vectors with shape (N, 3) where N is the number of
geometries. Each row specifies the (x, y, z) translation for one geometry
(after any transform is applied). If not provided, no additional translation
is applied beyond any transforms.
transforms : Optional[ArrayLike] = None
Optional array of 4x4 linear-only transform matrices with shape (N, 4, 4).
Each transform must be a valid homogeneous linear transform
(rotation/reflection/scale/shear) with no translation component.
Typical transforms can be created using ``Transformed.rotation``,
``Transformed.reflection``, or ``Transformed.scaling``.
Returns
-------
:class:`GeometryArray`
Array containing N copies of this geometry.
Notes
-----
- ``offsets`` represent all per-instance translation.
- ``transforms`` represent linear transforms only and must not contain translation.
- If both ``offsets`` and ``transforms`` are ``None``, the array contains a single
instance of the base geometry.
- If both are provided, they must have the same length and transforms are applied
before the translation given by offsets.
- Adjoint/autodiff is not currently supported for ``GeometryArray``.
Example
-------
>>> import tidy3d as td
>>> import numpy as np
>>> box = td.Box(size=(1, 1, 1))
>>> # Create a 2x2 grid of boxes using offsets
>>> offsets = [[0, 0, 0], [2, 0, 0], [0, 2, 0], [2, 2, 0]]
>>> array = box.array(offsets=offsets)
>>> # Create array using linear transforms only (rotation around z-axis)
>>> transforms = [np.eye(4), td.Transformed.rotation(np.pi/4, 2)]
>>> array = box.array(transforms=transforms)
>>> # Both None gives single instance of base geometry
>>> array = box.array()
"""
return GeometryArray(geometry=self, offsets=offsets, transforms=transforms)
""" Field and coordinate transformations """
[docs]
@staticmethod
def car_2_sph(x: float, y: float, z: float) -> tuple[float, float, float]:
"""Convert Cartesian to spherical coordinates.
Parameters
----------
x : float
x coordinate relative to ``local_origin``.
y : float
y coordinate relative to ``local_origin``.
z : float
z coordinate relative to ``local_origin``.
Returns
-------
tuple[float, float, float]
r, theta, and phi coordinates relative to ``local_origin``.
"""
r = np.sqrt(x**2 + y**2 + z**2)
theta = np.arccos(z / r)
phi = np.arctan2(y, x)
return r, theta, phi
[docs]
@staticmethod
def sph_2_car(r: float, theta: float, phi: float) -> tuple[float, float, float]:
"""Convert spherical to Cartesian coordinates.
Parameters
----------
r : float
radius.
theta : float
polar angle (rad) downward from x=y=0 line.
phi : float
azimuthal (rad) angle from y=z=0 line.
Returns
-------
tuple[float, float, float]
x, y, and z coordinates relative to ``local_origin``.
"""
r_sin_theta = r * np.sin(theta)
x = r_sin_theta * np.cos(phi)
y = r_sin_theta * np.sin(phi)
z = r * np.cos(theta)
return x, y, z
[docs]
@staticmethod
def sph_2_car_field(
f_r: float, f_theta: float, f_phi: float, theta: float, phi: float
) -> tuple[complex, complex, complex]:
"""Convert vector field components in spherical coordinates to cartesian.
Parameters
----------
f_r : float
radial component of the vector field.
f_theta : float
polar angle component of the vector fielf.
f_phi : float
azimuthal angle component of the vector field.
theta : float
polar angle (rad) of location of the vector field.
phi : float
azimuthal angle (rad) of location of the vector field.
Returns
-------
tuple[float, float, float]
x, y, and z components of the vector field in cartesian coordinates.
"""
sin_theta = np.sin(theta)
cos_theta = np.cos(theta)
sin_phi = np.sin(phi)
cos_phi = np.cos(phi)
f_x = f_r * sin_theta * cos_phi + f_theta * cos_theta * cos_phi - f_phi * sin_phi
f_y = f_r * sin_theta * sin_phi + f_theta * cos_theta * sin_phi + f_phi * cos_phi
f_z = f_r * cos_theta - f_theta * sin_theta
return f_x, f_y, f_z
[docs]
@staticmethod
def car_2_sph_field(
f_x: float, f_y: float, f_z: float, theta: float, phi: float
) -> tuple[complex, complex, complex]:
"""Convert vector field components in cartesian coordinates to spherical.
Parameters
----------
f_x : float
x component of the vector field.
f_y : float
y component of the vector fielf.
f_z : float
z component of the vector field.
theta : float
polar angle (rad) of location of the vector field.
phi : float
azimuthal angle (rad) of location of the vector field.
Returns
-------
tuple[float, float, float]
radial (s), elevation (theta), and azimuthal (phi) components
of the vector field in spherical coordinates.
"""
sin_theta = np.sin(theta)
cos_theta = np.cos(theta)
sin_phi = np.sin(phi)
cos_phi = np.cos(phi)
f_r = f_x * sin_theta * cos_phi + f_y * sin_theta * sin_phi + f_z * cos_theta
f_theta = f_x * cos_theta * cos_phi + f_y * cos_theta * sin_phi - f_z * sin_theta
f_phi = -f_x * sin_phi + f_y * cos_phi
return f_r, f_theta, f_phi
[docs]
@staticmethod
def kspace_2_sph(ux: float, uy: float, axis: Axis) -> tuple[float, float]:
"""Convert normalized k-space coordinates to angles.
Parameters
----------
ux : float
normalized kx coordinate.
uy : float
normalized ky coordinate.
axis : int
axis along which the observation plane is oriented.
Returns
-------
tuple[float, float]
theta and phi coordinates relative to ``local_origin``.
"""
phi_local = np.arctan2(uy, ux)
with np.errstate(invalid="ignore"):
theta_local = np.arcsin(np.sqrt(ux**2 + uy**2))
# Spherical coordinates rotation matrix reference:
# https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula#Matrix_notation
if axis == 2:
return theta_local, phi_local
x = np.cos(theta_local)
y = np.sin(theta_local) * np.cos(phi_local)
z = np.sin(theta_local) * np.sin(phi_local)
if axis == 1:
x, y, z = y, x, z
theta = np.arccos(z)
phi = np.arctan2(y, x)
return theta, phi
[docs]
@staticmethod
@verify_packages_import(["gdstk"])
def load_gds_vertices_gdstk(
gds_cell: Cell,
gds_layer: int,
gds_dtype: int | None = None,
gds_scale: PositiveFloat = 1.0,
) -> list[ArrayFloat2D]:
"""Load polygon vertices 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 polygon vertices
"""
# apply desired scaling and load the polygon vertices
if gds_dtype is not None:
# if both layer and datatype are specified, let gdstk do the filtering for better
# performance on large layouts
all_vertices = [
polygon.scale(gds_scale).points
for polygon in gds_cell.get_polygons(layer=gds_layer, datatype=gds_dtype)
]
else:
all_vertices = [
polygon.scale(gds_scale).points
for polygon in gds_cell.get_polygons()
if polygon.layer == gds_layer
]
# make sure something got loaded, otherwise error
if not all_vertices:
raise Tidy3dKeyError(
f"Couldn't load gds_cell, no vertices found at gds_layer={gds_layer} "
f"with specified gds_dtype={gds_dtype}."
)
return all_vertices
[docs]
@staticmethod
@verify_packages_import(["gdstk"])
def from_gds(
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",
merge_adjacent: bool = False,
) -> Geometry:
"""Import a ``gdstk.Cell`` and extrude it into a GeometryGroup.
Parameters
----------
gds_cell : gdstk.Cell
``gdstk.Cell`` containing 2D geometric data.
axis : int
Integer index defining the extrusion axis: 0 (x), 1 (y), or 2 (z).
slab_bounds: tuple[float, float]
Minimal and maximal positions of the extruded 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 (positive) or erosion (negative) amount to be applied to the original polygons.
sidewall_angle : float = 0
Angle of the extrusion sidewalls, away from the vertical direction, in radians. Positive
(negative) values result in slabs larger (smaller) at the base than at the top.
reference_plane : PlanePosition = "middle"
Reference position of the (dilated/eroded) polygons along the slab axis. One of
``"middle"`` (polygons correspond to the center of the slab bounds), ``"bottom"``
(minimal slab bound position), or ``"top"`` (maximal slab bound position). This value
has no effect if ``sidewall_angle == 0``.
merge_adjacent : bool = False
Merge polygons that become fractured into multiple adjacent GDS polygons, for example
due to the GDS vertex limit. Enable to import those fragments as a single merged shape.
Returns
-------
:class:`~tidy3d.Geometry`
Geometries created from the 2D data.
"""
import gdstk
if not isinstance(gds_cell, gdstk.Cell):
# Check if it might be a gdstk cell but gdstk is not found (should be caught by decorator)
# or if it's an entirely different type.
if "gdstk" in gds_cell.__class__.__name__.lower():
raise Tidy3dImportError(
"Module 'gdstk' not found. It is required to import gdstk cells."
)
raise Tidy3dImportError("Argument 'gds_cell' must be an instance of 'gdstk.Cell'.")
def iter_import_shapes(shape: Shapely) -> Iterable[Shapely]:
if shape.is_empty:
return
if shape.geom_type in {"MultiPolygon", "GeometryCollection"}:
for subshape in shape.geoms:
yield from iter_import_shapes(subshape)
else:
yield shape
def cleaned_shape(vertices: NDArray, consolidated_logger: Any) -> Shapely | None:
shape = shapely.set_precision(shapely.Polygon(vertices).buffer(0), POLY_GRID_SIZE)
if shape.is_empty:
consolidated_logger.warning(
"A GDS polygon collapsed during topology cleanup in "
"'Geometry.from_gds()' and will be skipped."
)
return None
return shape
geometries = []
with log as consolidated_logger:
gds_loader_fn = Geometry.load_gds_vertices_gdstk
all_vertices = gds_loader_fn(gds_cell, gds_layer, gds_dtype, gds_scale)
if merge_adjacent:
shapes = []
for vertices in all_vertices:
shape = cleaned_shape(vertices, consolidated_logger)
if shape is not None:
shapes.append(shape)
if len(shapes) > 1:
shapes = [shapely.set_precision(shapely.union_all(shapes), POLY_GRID_SIZE)]
import_shapes = (
import_shape for shape in shapes for import_shape in iter_import_shapes(shape)
)
else:
import_shapes = (
import_shape
for vertices in all_vertices
for shape in [cleaned_shape(vertices, consolidated_logger)]
if shape is not None
for import_shape in iter_import_shapes(shape)
)
for import_shape in import_shapes:
try:
geometries.append(
from_shapely(
import_shape,
axis,
slab_bounds,
dilation,
sidewall_angle,
reference_plane,
)
)
except ValidationError as error:
consolidated_logger.warning(str(error))
except Tidy3dError as error:
consolidated_logger.warning(str(error))
if not geometries:
raise SetupError(
"Couldn't import any valid geometries from 'gds_cell' at "
f"gds_layer={gds_layer} with specified gds_dtype={gds_dtype}. "
"All polygons were skipped during cleanup or failed conversion."
)
return geometries[0] if len(geometries) == 1 else GeometryGroup(geometries=geometries)
[docs]
@staticmethod
def from_shapely(
shape: Shapely,
axis: Axis,
slab_bounds: tuple[float, float],
dilation: float = 0.0,
sidewall_angle: float = 0,
reference_plane: PlanePosition = "middle",
) -> Geometry:
"""Convert a shapely primitive into a geometry instance by extrusion.
Parameters
----------
shape : shapely.geometry.base.BaseGeometry
Shapely primitive to be converted. It must be a linear ring, a polygon or a collection
of any of those.
axis : int
Integer index defining the extrusion axis: 0 (x), 1 (y), or 2 (z).
slab_bounds: tuple[float, float]
Minimal and maximal positions of the extruded slab along ``axis``.
dilation : float
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 extrusion sidewalls, away from the vertical direction, in radians. Positive
(negative) values result in slabs larger (smaller) at the base than at the top.
reference_plane : PlanePosition = "middle"
Reference position of the (dilated/eroded) polygons along the slab axis. One of
``"middle"`` (polygons correspond to the center of the slab bounds), ``"bottom"``
(minimal slab bound position), or ``"top"`` (maximal slab bound position). This value
has no effect if ``sidewall_angle == 0``.
Returns
-------
:class:`~tidy3d.Geometry`
Geometry extruded from the 2D data.
"""
return from_shapely(shape, axis, slab_bounds, dilation, sidewall_angle, reference_plane)
[docs]
@verify_packages_import(["gdstk"])
def to_gdstk(
self,
x: float | None = None,
y: float | None = None,
z: float | None = None,
gds_layer: NonNegativeInt = 0,
gds_dtype: NonNegativeInt = 0,
) -> list:
"""Convert a Geometry object's planar slice to a .gds type polygon.
Parameters
----------
x : float = None
Position of plane in x direction, only one of x,y,z can be specified to define plane.
y : float = None
Position of plane in y direction, only one of x,y,z can be specified to define plane.
z : float = None
Position of plane in z direction, only one of x,y,z can be specified to define plane.
gds_layer : int = 0
Layer index to use for the shapes stored in the .gds file.
gds_dtype : int = 0
Data-type index to use for the shapes stored in the .gds file.
Return
------
List
List of `gdstk.Polygon`.
"""
import gdstk
shapes = self.intersections_plane(x=x, y=y, z=z)
polygons = []
for shape in shapes:
for vertices in vertices_from_shapely(shape):
if len(vertices) == 1:
polygons.append(gdstk.Polygon(vertices[0], gds_layer, gds_dtype))
else:
polygons.extend(
gdstk.boolean(
vertices[:1],
vertices[1:],
"not",
layer=gds_layer,
datatype=gds_dtype,
)
)
return polygons
[docs]
@verify_packages_import(["gdstk"])
def to_gds(
self,
cell: Cell,
x: float | None = None,
y: float | None = None,
z: float | None = None,
gds_layer: NonNegativeInt = 0,
gds_dtype: NonNegativeInt = 0,
) -> None:
"""Append a Geometry object's planar slice to a .gds cell.
Parameters
----------
cell : ``gdstk.Cell``
Cell object to which the generated polygons are added.
x : float = None
Position of plane in x direction, only one of x,y,z can be specified to define plane.
y : float = None
Position of plane in y direction, only one of x,y,z can be specified to define plane.
z : float = None
Position of plane in z direction, only one of x,y,z can be specified to define plane.
gds_layer : int = 0
Layer index to use for the shapes stored in the .gds file.
gds_dtype : int = 0
Data-type index to use for the shapes stored in the .gds file.
"""
import gdstk
if not isinstance(cell, gdstk.Cell):
if "gdstk" in cell.__class__.__name__.lower():
raise Tidy3dImportError(
"Module 'gdstk' not found. It is required to export shapes to gdstk cells."
)
raise Tidy3dImportError("Argument 'cell' must be an instance of 'gdstk.Cell'.")
polygons = self.to_gdstk(x=x, y=y, z=z, gds_layer=gds_layer, gds_dtype=gds_dtype)
if polygons:
cell.add(*polygons)
[docs]
@verify_packages_import(["gdstk"])
def to_gds_file(
self,
fname: PathLike,
x: float | None = None,
y: float | None = None,
z: float | None = None,
gds_layer: NonNegativeInt = 0,
gds_dtype: NonNegativeInt = 0,
gds_cell_name: str = "MAIN",
gds_precision: PositiveFloat = 1e-3,
) -> None:
"""Export a Geometry object's planar slice to a .gds file.
Parameters
----------
fname : PathLike
Full path to the .gds file to save the :class:`~tidy3d.Geometry` slice to.
x : float = None
Position of plane in x direction, only one of x,y,z can be specified to define plane.
y : float = None
Position of plane in y direction, only one of x,y,z can be specified to define plane.
z : float = None
Position of plane in z direction, only one of x,y,z can be specified to define plane.
gds_layer : int = 0
Layer index to use for the shapes stored in the .gds file.
gds_dtype : int = 0
Data-type index to use for the shapes stored in the .gds file.
gds_cell_name : str = 'MAIN'
Name of the cell created in the .gds file to store the geometry.
gds_precision : float = 1e-3
Coordinate precision for the written GDS file in micrometers. The default matches
the gdstk default of ``1e-9`` meters. If the requested precision is too fine for the
written slice coordinates, export raises :class:`.SetupError`. The minimum safe value
scales with the maximum absolute written planar coordinate as
``max_abs_coord / (2**31 - 1)``.
"""
try:
import gdstk
except ImportError as e:
raise Tidy3dImportError(
format_chained_exception_message(
"Python module 'gdstk' not found. To export geometries to .gds files, "
"please install it",
e,
)
) from e
polygons = self.to_gdstk(
x=x,
y=y,
z=z,
gds_layer=gds_layer,
gds_dtype=gds_dtype,
)
gds_precision = self._validate_gds_precision(
polygons=polygons,
gds_precision=gds_precision,
context="Geometry.to_gds_file()",
)
library = gdstk.Library(unit=1e-6, precision=gds_precision * 1e-6)
cell = library.new_cell(gds_cell_name)
if polygons:
cell.add(*polygons)
fname = pathlib.Path(fname)
fname.parent.mkdir(parents=True, exist_ok=True)
library.write_gds(fname)
def _compute_derivatives(self, derivative_info: DerivativeInfo) -> AutogradFieldMap:
"""Compute the adjoint derivatives for this object."""
raise NotImplementedError(f"Can't compute derivative for 'Geometry': '{type(self)}'.")
def _resolve_autograd_route(self, field_path: tuple[Any, ...]) -> AutogradRoute:
"""Resolve and validate one traced geometry path for adjoint routing."""
return validate_traced_path(
parameter_kind="geometry",
owner_kind="geometry type",
owner_name=type(self).__name__,
field_path=field_path,
supported_paths=self._traced_supported_paths,
supported_parameters=type(self)._traced_autograd_supported_parameters(),
)
def _as_union(self) -> list[Geometry]:
"""Return a list of geometries that, united, make up the given geometry."""
if isinstance(self, GeometryGroup):
return self.geometries
if isinstance(self, ClipOperation) and self.operation == "union":
return (self.geometry_a, self.geometry_b)
return (self,)
def __add__(self, other: int | Geometry) -> Self | GeometryGroup:
"""Union of geometries"""
# This allows the user to write sum(geometries...) with the default start=0
if isinstance(other, int):
return self
if not isinstance(other, Geometry):
return NotImplemented
return GeometryGroup(geometries=self._as_union() + other._as_union())
def __radd__(self, other: int | Geometry) -> Self | GeometryGroup:
"""Union of geometries"""
# This allows the user to write sum(geometries...) with the default start=0
if isinstance(other, int):
return self
if not isinstance(other, Geometry):
return NotImplemented
return GeometryGroup(geometries=other._as_union() + self._as_union())
def __or__(self, other: Geometry) -> GeometryGroup:
"""Union of geometries"""
if not isinstance(other, Geometry):
return NotImplemented
return GeometryGroup(geometries=self._as_union() + other._as_union())
def __mul__(self, other: Geometry) -> ClipOperation:
"""Intersection of geometries"""
if not isinstance(other, Geometry):
return NotImplemented
return ClipOperation(operation="intersection", geometry_a=self, geometry_b=other)
def __and__(self, other: Geometry) -> ClipOperation:
"""Intersection of geometries"""
if not isinstance(other, Geometry):
return NotImplemented
return ClipOperation(operation="intersection", geometry_a=self, geometry_b=other)
def __sub__(self, other: Geometry) -> ClipOperation:
"""Difference of geometries"""
if not isinstance(other, Geometry):
return NotImplemented
return ClipOperation(operation="difference", geometry_a=self, geometry_b=other)
def __xor__(self, other: Geometry) -> ClipOperation:
"""Symmetric difference of geometries"""
if not isinstance(other, Geometry):
return NotImplemented
return ClipOperation(operation="symmetric_difference", geometry_a=self, geometry_b=other)
def __pos__(self) -> Self:
"""No op"""
return self
def __neg__(self) -> ClipOperation:
"""Opposite of a geometry"""
return ClipOperation(
operation="difference", geometry_a=Box(size=(inf, inf, inf)), geometry_b=self
)
def __invert__(self) -> ClipOperation:
"""Opposite of a geometry"""
return ClipOperation(
operation="difference", geometry_a=Box(size=(inf, inf, inf)), geometry_b=self
)
""" Abstract subclasses """
[docs]
class Centered(Geometry, ABC):
"""Geometry with a well defined center."""
center: TracedCoordinate = Field(
(0.0, 0.0, 0.0),
title="Center",
description="Center of object in x, y, and z.",
json_schema_extra={"units": MICROMETER},
)
@field_validator("center")
@classmethod
def _center_not_inf(cls, val: tuple[float, float, float]) -> tuple[float, float, float]:
"""Make sure center is not infinitiy."""
if any(np.isinf(v) for v in val):
raise ValidationError("center can not contain td.inf terms.")
return val
[docs]
class SimplePlaneIntersection(Geometry, ABC):
"""A geometry where intersections with an axis aligned plane may be computed efficiently."""
[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.
Checks special cases before relying on the complete computation.
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. If ``None``, uses
high-quality visualization settings.
section_tolerance_2d : bool = False
See :meth:`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>`_.
"""
# Check if normal is a special case, where the normal is aligned with an axis.
if np.sum(np.isclose(normal, 0.0)) == 2:
axis = np.argmax(np.abs(normal)).item()
coord = "xyz"[axis]
kwargs = {coord: origin[axis]}
section = self.intersections_plane(
cleanup=cleanup,
quad_segs=quad_segs,
section_tolerance_2d=section_tolerance_2d,
**kwargs,
)
# Apply transformation in the plane by removing row and column
to_2D_in_plane = np.delete(np.delete(to_2D, 2, 0), axis, 1)
def transform(p_array: NDArray) -> NDArray:
return np.dot(
np.hstack((p_array, np.ones((p_array.shape[0], 1)))), to_2D_in_plane.T
)[:, :2]
transformed_section = shapely.transform(section, transformation=transform)
return transformed_section
# Otherwise compute the arbitrary intersection
return self._do_intersections_tilted_plane(
normal=normal, origin=origin, to_2D=to_2D, quad_segs=quad_segs
)
@abstractmethod
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.
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>`_.
"""
[docs]
class Planar(SimplePlaneIntersection, Geometry, ABC):
"""Geometry with one ``axis`` that is slab-like with thickness ``height``."""
axis: Axis = Field(
2,
title="Axis",
description="Specifies dimension of the planar axis (0,1,2) -> (x,y,z).",
)
sidewall_angle: TracedFloat = Field(
0.0,
title="Sidewall angle",
description="Angle of the sidewall. "
"``sidewall_angle=0`` (default) specifies a vertical wall; "
"``0<sidewall_angle<np.pi/2`` specifies a shrinking cross section "
"along the ``axis`` direction; "
"and ``-np.pi/2<sidewall_angle<0`` specifies an expanding cross section "
"along the ``axis`` direction.",
json_schema_extra={"units": RADIAN},
)
reference_plane: PlanePosition = Field(
"middle",
title="Reference plane for cross section",
description="The position of the plane where the supplied cross section are "
"defined. The plane is perpendicular to the ``axis``. "
"The plane is located at the ``bottom``, ``middle``, or ``top`` of the "
"geometry with respect to the axis. "
"E.g. if ``axis=1``, ``bottom`` refers to the negative side of the y-axis, and "
"``top`` refers to the positive side of the y-axis.",
)
[docs]
@field_validator("sidewall_angle")
@classmethod
def validate_angle(cls, val: float) -> float:
lower_bound = -np.pi / 2
upper_bound = np.pi / 2
if (val <= lower_bound) or (val >= upper_bound):
# u03C0 is unicode for pi
raise ValidationError(f"Sidewall angle ({val}) must be between -π/2 and π/2 rad.")
return val
@property
@abstractmethod
def center_axis(self) -> float:
"""Gets the position of the center of the geometry in the out of plane dimension."""
@property
@abstractmethod
def length_axis(self) -> float:
"""Gets the length of the geometry along the out of plane dimension."""
@property
def finite_length_axis(self) -> float:
"""Gets the length of the geometry along the out of plane dimension.
If the length is td.inf, return ``LARGE_NUMBER``
"""
return min(self.length_axis, LARGE_NUMBER)
@property
def reference_axis_pos(self) -> float:
"""Coordinate along the slab axis at the reference plane.
Returns the axis coordinate corresponding to the selected
reference_plane:
- "bottom": lower bound of slab_bounds
- "middle": center_axis
- "top": upper bound of slab_bounds
"""
if self.reference_plane == "bottom":
return self.slab_bounds[0]
if self.reference_plane == "top":
return self.slab_bounds[1]
# default to middle
return self.center_axis
[docs]
def intersections_plane(
self,
x: float | None = None,
y: float | None = None,
z: float | None = None,
cleanup: bool = True,
quad_segs: int | None = None,
section_tolerance_2d: bool = False,
) -> list[Shapely]:
"""Returns shapely geometry at plane specified by one non None value of x,y,z.
Parameters
----------
x : float
Position of plane in x direction, only one of x,y,z can be specified to define plane.
y : float
Position of plane in y direction, only one of x,y,z can be specified to define plane.
z : float
Position of plane in z direction, only one of x,y,z can be specified to define plane.
cleanup : bool = True
If True, removes extremely small features from each polygon's boundary.
quad_segs : Optional[int] = None
Number of segments used to discretize circular shapes. If ``None``, uses
high-quality visualization settings.
section_tolerance_2d : bool = False
See :meth:`Geometry.intersections_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>``.
"""
axis, position = self.parse_xyz_kwargs(x=x, y=y, z=z)
use_2d_tolerance = (
section_tolerance_2d
and axis == self.axis
and np.isclose(self.length_axis, 0.0, rtol=fp_eps, atol=fp_eps)
)
if not self.intersects_axis_position(axis, position, section_tolerance_2d=use_2d_tolerance):
return []
if use_2d_tolerance and np.isclose(
position, self.reference_axis_pos, rtol=fp_eps, atol=fp_eps
):
position = self.reference_axis_pos
if axis == self.axis:
return self._intersections_normal(position, quad_segs=quad_segs)
return self._intersections_side(position, axis)
@abstractmethod
def _intersections_normal(self, z: float, quad_segs: int | None = None) -> list:
"""Find shapely geometries intersecting planar geometry with axis normal to slab.
Parameters
----------
z : float
Position along the axis normal to slab
quad_segs : Optional[int] = None
Number of segments used to discretize circular shapes. If ``None``, uses
high-quality visualization settings.
Returns
-------
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>`_.
"""
@abstractmethod
def _intersections_side(self, position: float, axis: Axis) -> list[Shapely]:
"""Find shapely geometries intersecting planar geometry with axis orthogonal to plane.
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>`_.
"""
def _order_axis(self, axis: int) -> int:
"""Order the axis as if self.axis is along z-direction.
Parameters
----------
axis : int
Integer index into the structure's planar axis.
Returns
-------
int
New index of axis.
"""
axis_index = [0, 1]
axis_index.insert(self.axis, 2)
return axis_index[axis]
def _order_by_axis(self, plane_val: Any, axis_val: Any, axis: int) -> tuple[Any, Any]:
"""Orders a value in the plane and value along axis in correct (x,y) order for plotting.
Note: sometimes if axis=1 and we compute cross section values orthogonal to axis,
they can either be x or y in the plots.
This function allows one to figure out the ordering.
Parameters
----------
plane_val : Any
The value in the planar coordinate.
axis_val : Any
The value in the ``axis`` coordinate.
axis : int
Integer index into the structure's planar axis.
Returns
-------
``(Any, Any)``
The two planar coordinates in this new coordinate system.
"""
vals = 3 * [plane_val]
vals[self.axis] = axis_val
_, (val_x, val_y) = self.pop_axis(vals, axis=axis)
return val_x, val_y
@cached_property
def _tanq(self) -> float:
"""Value of ``tan(sidewall_angle)``.
The (possibliy infinite) geometry offset is given by ``_tanq * length_axis``.
"""
return np.tan(self.sidewall_angle)
[docs]
class Circular(Geometry):
"""Geometry with circular characteristics (specified by a radius)."""
radius: NonNegativeFloat = Field(
title="Radius",
description="Radius of geometry.",
json_schema_extra={"units": MICROMETER},
)
@field_validator("radius")
@classmethod
def _radius_not_inf(cls, val: float) -> float:
"""Make sure center is not infinitiy."""
if np.isinf(val):
raise ValidationError("radius can not be 'td.inf'.")
return val
def _intersect_dist(self, position: float, z0: float) -> float:
"""Distance between points on circle at z=position where center of circle at z=z0.
Parameters
----------
position : float
position along z.
z0 : float
center of circle in z.
Returns
-------
float
Distance between points on the circle intersecting z=z, if no points, ``None``.
"""
dz = np.abs(z0 - position)
if dz > self.radius:
return None
return 2 * np.sqrt(self.radius**2 - dz**2)
"""Primitive classes"""
[docs]
class Box(SimplePlaneIntersection, Centered):
"""Rectangular prism.
Also base class for :class:`.Simulation`, :class:`Monitor`, and :class:`Source`.
Example
-------
>>> b = Box(center=(1,2,3), size=(2,2,2))
"""
_traced_supported_paths: ClassVar[tuple[PathType, ...]] = traced_paths(
"center",
"size",
*indexed_traced_paths("center", 3),
*indexed_traced_paths("size", 3),
)
size: TracedSize = Field(
title="Size",
description="Size in x, y, and z directions.",
json_schema_extra={"units": MICROMETER},
)
[docs]
@classmethod
def from_bounds(cls, rmin: Coordinate, rmax: Coordinate, **kwargs: Any) -> Self:
"""Constructs a :class:`~tidy3d.Box` from minimum and maximum coordinate bounds
Parameters
----------
rmin : tuple[float, float, float]
(x, y, z) coordinate of the minimum values.
rmax : tuple[float, float, float]
(x, y, z) coordinate of the maximum values.
Example
-------
>>> b = Box.from_bounds(rmin=(-1, -2, -3), rmax=(3, 2, 1))
"""
center = tuple(cls._get_center(pt_min, pt_max) for pt_min, pt_max in zip(rmin, rmax))
size = tuple((pt_max - pt_min) for pt_min, pt_max in zip(rmin, rmax))
return cls(center=center, size=size, **kwargs)
@cached_property
def _normal_axis(self) -> Axis:
"""Axis normal to the Box. Errors if box is not planar."""
if self.size.count(0.0) != 1:
raise ValidationError(
f"Tried to get 'normal_axis' of 'Box' that is not planar. Given 'size={self.size}.'"
)
return self.size.index(0.0)
@staticmethod
def _surface_keys(size: Size) -> tuple[list[str], set[int]]:
"""Return the canonical surface keys and indices dropped for infinite dimensions."""
surface_keys = [coord + direction for coord in "xyz" for direction in "-+"]
del_idx = {
2 * idx + offset for idx, _size in enumerate(size) if _size == inf for offset in (0, 1)
}
surface_keys = [key for idx, key in enumerate(surface_keys) if idx not in del_idx]
return surface_keys, del_idx
[docs]
@classmethod
def surfaces(cls, size: Size, center: Coordinate, **kwargs: Any) -> list[Self]:
"""Returns a list of 6 :class:`~tidy3d.Box` instances corresponding to each surface of a 3D volume.
The output surfaces are stored in the order [x-, x+, y-, y+, z-, z+], where x, y, and z
denote which axis is perpendicular to that surface, while "-" and "+" denote the direction
of the normal vector of that surface. If a name is provided, each output surface's name
will be that of the provided name appended with the above symbols. E.g., if the provided
name is "box", the x+ surfaces's name will be "box_x+".
Parameters
----------
size : tuple[float, float, float]
Size of object in x, y, and z directions.
center : tuple[float, float, float]
Center of object in x, y, and z.
Example
-------
>>> b = Box.surfaces(size=(1, 2, 3), center=(3, 2, 1))
"""
if any(s == 0.0 for s in size):
raise SetupError(
"Can't generate surfaces for the given object because it has zero volume."
)
bounds = Box(center=center, size=size).bounds
# Set up geometry data and names for each surface:
centers = [list(center) for _ in range(6)]
sizes = [list(size) for _ in range(6)]
surface_index = 0
for dim_index in range(3):
for min_max_index in range(2):
new_center = centers[surface_index]
new_size = sizes[surface_index]
new_center[dim_index] = bounds[min_max_index][dim_index]
new_size[dim_index] = 0.0
centers[surface_index] = new_center
sizes[surface_index] = new_size
surface_index += 1
surface_keys, del_idx = cls._surface_keys(size)
name_base = kwargs.pop("name", "")
kwargs.pop("normal_dir", None)
def del_items(items: Iterable, indices: set[int]) -> list:
"""Delete list items at indices."""
return [i for j, i in enumerate(items) if j not in indices]
centers = del_items(centers, del_idx)
sizes = del_items(sizes, del_idx)
names = [name_base + "_" + surface_key for surface_key in surface_keys]
normal_dirs = [surface_key[-1] for surface_key in surface_keys]
surfaces = []
for _cent, _size, _name, _normal_dir in zip(centers, sizes, names, normal_dirs):
if "normal_dir" in cls.model_fields:
kwargs["normal_dir"] = _normal_dir
if "name" in cls.model_fields:
kwargs["name"] = _name
surface = cls(center=_cent, size=_size, **kwargs)
surfaces.append(surface)
return surfaces
[docs]
@classmethod
def surfaces_with_exclusion(cls, size: Size, center: Coordinate, **kwargs: Any) -> list[Self]:
"""Returns a list of 6 :class:`~tidy3d.Box` instances corresponding to each surface of a 3D volume.
The output surfaces are stored in the order [x-, x+, y-, y+, z-, z+], where x, y, and z
denote which axis is perpendicular to that surface, while "-" and "+" denote the direction
of the normal vector of that surface. If a name is provided, each output surface's name
will be that of the provided name appended with the above symbols. E.g., if the provided
name is "box", the x+ surfaces's name will be "box_x+". If ``kwargs`` contains an
``exclude_surfaces`` parameter, the returned list of surfaces will not include the excluded
surfaces. Otherwise, the behavior is identical to that of ``surfaces()``.
Parameters
----------
size : tuple[float, float, float]
Size of object in x, y, and z directions.
center : tuple[float, float, float]
Center of object in x, y, and z.
Example
-------
>>> b = Box.surfaces_with_exclusion(
... size=(1, 2, 3), center=(3, 2, 1), exclude_surfaces=["x-"]
... )
"""
exclude_surfaces = kwargs.pop("exclude_surfaces", None)
surfaces = cls.surfaces(size=size, center=center, **kwargs)
if exclude_surfaces:
surface_keys, _ = cls._surface_keys(size)
exclude_surfaces = set(exclude_surfaces)
surfaces = [
surf
for surf, surface_key in zip(surfaces, surface_keys)
if surface_key not in exclude_surfaces
]
return surfaces
@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 Box 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
(x0, y0, z0), (x1, y1, z1) = self.bounds
vertices = [
(x0, y0, z0), # 0
(x0, y0, z1), # 1
(x0, y1, z0), # 2
(x0, y1, z1), # 3
(x1, y0, z0), # 4
(x1, y0, z1), # 5
(x1, y1, z0), # 6
(x1, y1, z1), # 7
]
faces = [
(0, 1, 3, 2), # -x
(4, 6, 7, 5), # +x
(0, 4, 5, 1), # -y
(2, 3, 7, 6), # +y
(0, 2, 6, 4), # -z
(1, 5, 7, 3), # +z
]
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
[docs]
def intersections_plane(
self,
x: float | None = None,
y: float | None = None,
z: float | None = None,
cleanup: bool = True,
quad_segs: int | None = None,
section_tolerance_2d: bool = False,
) -> list[Shapely]:
"""Returns shapely geometry at plane specified by one non None value of x,y,z.
Parameters
----------
x : float = None
Position of plane in x direction, only one of x,y,z can be specified to define plane.
y : float = None
Position of plane in y direction, only one of x,y,z can be specified to define plane.
z : float = None
Position of plane in z direction, only one of x,y,z can be specified to define plane.
cleanup : bool = True
If True, removes extremely small features from each polygon's boundary.
quad_segs : Optional[int] = None
Number of segments used to discretize circular shapes. Not used for Box geometry.
section_tolerance_2d : bool = False
See :meth:`Geometry.intersections_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>`_.
"""
axis, position = self.parse_xyz_kwargs(x=x, y=y, z=z)
use_2d_tolerance = section_tolerance_2d and np.isclose(
self.size[axis], 0.0, rtol=fp_eps, atol=fp_eps
)
if not self.intersects_axis_position(axis, position, section_tolerance_2d=use_2d_tolerance):
return []
z0, (x0, y0) = self.pop_axis(self.center, axis=axis)
Lz, (Lx, Ly) = self.pop_axis(self.size, axis=axis)
if use_2d_tolerance and np.isclose(position, z0, rtol=fp_eps, atol=fp_eps):
position = z0
dz = np.abs(z0 - position)
if dz > Lz / 2 + fp_eps:
return []
minx = x0 - Lx / 2
miny = y0 - Ly / 2
maxx = x0 + Lx / 2
maxy = y0 + Ly / 2
# handle case where the box vertices are identical
if np.isclose(minx, maxx) and np.isclose(miny, maxy):
return [self.make_shapely_point(minx, miny)]
return [self.make_shapely_box(minx, miny, maxx, maxy)]
[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.
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)
x0, y0, z0 = self.center
Lx, Ly, Lz = self.size
dist_x = np.abs(x - x0)
dist_y = np.abs(y - y0)
dist_z = np.abs(z - z0)
return (dist_x <= Lx / 2) * (dist_y <= Ly / 2) * (dist_z <= Lz / 2)
[docs]
def intersections_with(
self,
other: Geometry,
cleanup: bool = True,
quad_segs: int | None = None,
section_tolerance_2d: bool = False,
) -> list[Shapely]:
"""Returns list of shapely geometries representing the intersections of the geometry with
this 2D box.
Parameters
----------
other : :class:`~tidy3d.Geometry`
Geometry to intersect with.
cleanup : bool = True
If True, removes extremely small features from each polygon's boundary.
quad_segs : Optional[int] = None
Number of segments used to discretize circular shapes. If ``None``, uses
high-quality visualization settings.
section_tolerance_2d : bool = False
See :meth:`Geometry.intersections_plane`.
Returns
-------
list[shapely.geometry.base.BaseGeometry]
List of 2D shapes that intersect this 2D box.
For more details refer to
`Shapely's Documentation <https://shapely.readthedocs.io/en/stable/project.html>`_.
"""
# Verify 2D
if self.size.count(0.0) != 1:
raise ValidationError(
"Intersections with other geometry are only calculated from a 2D box."
)
# Don't bother if the geometry doesn't intersect the self at all.
# Plotting opts into the section-aware precheck so zero-thickness 2D shapes
# that are only off by ``fp_eps`` still reach ``intersections_plane()`` below.
normal_ind = self.size.index(0.0)
if section_tolerance_2d:
if not other.intersects_axis_position(
normal_ind,
self.center[normal_ind],
section_tolerance_2d=True,
):
return []
elif not other.intersects(self):
return []
# get list of Shapely shapes that intersect at the self
dim = "xyz"[normal_ind]
pos = self.center[normal_ind]
xyz_kwargs = {dim: pos}
shapes_plane = other.intersections_plane(
cleanup=cleanup,
quad_segs=quad_segs,
section_tolerance_2d=section_tolerance_2d,
**xyz_kwargs,
)
# intersect all shapes with the input self
bs_min, bs_max = (self.pop_axis(bounds, axis=normal_ind)[1] for bounds in self.bounds)
shapely_box = self.make_shapely_box(bs_min[0], bs_min[1], bs_max[0], bs_max[1])
shapely_box = Geometry.evaluate_inf_shape(shapely_box)
return [Geometry.evaluate_inf_shape(shape) & shapely_box for shape in shapes_plane]
[docs]
def slightly_enlarged_copy(self) -> Box:
"""Box size slightly enlarged around machine precision."""
size = [increment_float(orig_length, 1) for orig_length in self.size]
return self.updated_copy(size=size)
[docs]
def padded_copy(
self,
x: tuple[pydantic.NonNegativeFloat, pydantic.NonNegativeFloat] | None = None,
y: tuple[pydantic.NonNegativeFloat, pydantic.NonNegativeFloat] | None = None,
z: tuple[pydantic.NonNegativeFloat, pydantic.NonNegativeFloat] | None = None,
) -> Box:
"""Created a padded copy of a :class:`~tidy3d.Box` instance.
Parameters
----------
x : Optional[tuple[pydantic.NonNegativeFloat, pydantic.NonNegativeFloat]] = None
Padding sizes at the left and right boundaries of the box along x-axis.
y : Optional[tuple[pydantic.NonNegativeFloat, pydantic.NonNegativeFloat]] = None
Padding sizes at the left and right boundaries of the box along y-axis.
z : Optional[tuple[pydantic.NonNegativeFloat, pydantic.NonNegativeFloat]] = None
Padding sizes at the left and right boundaries of the box along z-axis.
Returns
-------
Box
Padded instance of :class:`~tidy3d.Box`.
"""
# Validate that padding values are non-negative
for axis_name, axis_padding in zip(("x", "y", "z"), (x, y, z)):
if axis_padding is not None:
if not isinstance(axis_padding, (tuple, list)) or len(axis_padding) != 2:
raise ValueError(f"Padding for {axis_name}-axis must be a tuple of two values.")
if any(p < 0 for p in axis_padding):
raise ValueError(
f"Padding values for {axis_name}-axis must be non-negative. Got {axis_padding}."
)
rmin, rmax = self.bounds
def bound_array(arrs: ArrayLike, idx: int) -> NDArray:
return np.array([(a[idx] if a is not None else 0) for a in arrs])
# parse padding sizes for simulation
drmin = bound_array((x, y, z), 0)
drmax = bound_array((x, y, z), 1)
rmin = np.array(rmin) - drmin
rmax = np.array(rmax) + drmax
return Box.from_bounds(rmin=rmin, rmax=rmax)
@cached_property
def bounds(self) -> Bound:
"""Returns bounding box min and max coordinates.
Returns
-------
tuple[float, float, float], tuple[float, float float]
Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
"""
size = self.size
center = self.center
coord_min = tuple(c - s / 2 for (s, c) in zip(size, center))
coord_max = tuple(c + s / 2 for (s, c) in zip(size, center))
return (coord_min, coord_max)
@cached_property
def geometry(self) -> Box:
""":class:`~tidy3d.Box` representation of self (used for subclasses of Box).
Returns
-------
:class:`~tidy3d.Box`
Instance of :class:`~tidy3d.Box` representing self's geometry.
"""
return Box(center=self.center, size=self.size)
@cached_property
def zero_dims(self) -> list[Axis]:
"""A list of axes along which the :class:`~tidy3d.Box` is zero-sized."""
return [dim for dim, size in enumerate(self.size) if size == 0]
@cached_property
def _normal_2dmaterial(self) -> Axis:
"""Get the normal to the given geometry, checking that it is a 2D geometry."""
if np.count_nonzero(self.size) != 2:
raise ValidationError(
"'Medium2D' requires exactly one of the 'Box' dimensions to have size zero."
)
return self.size.index(0)
def _update_from_bounds(self, bounds: tuple[float, float], axis: Axis) -> Box:
"""Returns an updated geometry which has been transformed to fit within ``bounds``
along the ``axis`` direction."""
new_center = list(self.center)
new_center[axis] = (bounds[0] + bounds[1]) / 2
new_size = list(self.size)
new_size[axis] = bounds[1] - bounds[0]
return self.updated_copy(center=tuple(new_center), size=tuple(new_size))
def _plot_arrow(
self,
direction: tuple[float, float, float],
x: float | None = None,
y: float | None = None,
z: float | None = None,
color: str | None = None,
alpha: float | None = None,
bend_radius: float | None = None,
bend_axis: Axis = None,
both_dirs: bool = False,
ax: Ax = None,
arrow_base: Coordinate = None,
) -> Ax:
"""Adds an arrow to the axis if with options if certain conditions met.
Parameters
----------
direction: tuple[float, float, float]
Normalized vector describing the arrow direction.
x : float = None
Position of plotting plane in x direction.
y : float = None
Position of plotting plane in y direction.
z : float = None
Position of plotting plane in z direction.
color : str = None
Color of the arrow.
alpha : float = None
Opacity of the arrow (0, 1)
bend_radius : float = None
Radius of curvature for this arrow.
bend_axis : Axis = None
Axis of curvature of ``bend_radius``.
both_dirs : bool = False
If True, plots an arrow pointing in direction and one in -direction.
arrow_base : :class:`.Coordinate` = None
Custom base of the arrow. Uses the geometry's center if not provided.
Returns
-------
matplotlib.axes._subplots.Axes
The matplotlib axes with the arrow added.
"""
from matplotlib import patches
from tidy3d.components.viz.styles import arrow_style
plot_axis, _ = self.parse_xyz_kwargs(x=x, y=y, z=z)
_, (dx, dy) = self.pop_axis(direction, axis=plot_axis)
# conditions to check to determine whether to plot arrow, taking into account the
# possibility of a custom arrow base
arrow_intersecting_plane = (
len(self.intersections_plane(x=x, y=y, z=z, section_tolerance_2d=True)) > 0
)
center = self.center
if arrow_base:
arrow_intersecting_plane = arrow_intersecting_plane and any(
a == b for a, b in zip(arrow_base, [x, y, z])
)
center = arrow_base
_, (dx, dy) = self.pop_axis(direction, axis=plot_axis)
components_in_plane = any(not np.isclose(component, 0) for component in (dx, dy))
# plot if arrow in plotting plane and some non-zero component can be displayed.
if arrow_intersecting_plane and components_in_plane:
_, (x0, y0) = self.pop_axis(center, axis=plot_axis)
# Reasonable value for temporary arrow size. The correct size and direction
# have to be calculated after all transforms have been set. That is why we
# use a callback to do these calculations only at the drawing phase.
xmin, xmax = ax.get_xlim()
ymin, ymax = ax.get_ylim()
v_x = (xmax - xmin) / 10
v_y = (ymax - ymin) / 10
directions = (1.0, -1.0) if both_dirs else (1.0,)
for sign in directions:
arrow = patches.FancyArrowPatch(
(x0, y0),
(x0 + v_x, y0 + v_y),
arrowstyle=arrow_style(),
color=color,
alpha=alpha,
zorder=np.inf,
)
# Don't draw this arrow until it's been reshaped
arrow.set_visible(False)
callback = self._arrow_shape_cb(
arrow, (x0, y0), (dx, dy), sign, bend_radius if bend_axis == plot_axis else None
)
callback_id = ax.figure.canvas.mpl_connect("draw_event", callback)
# Store a reference to the callback because mpl_connect does not.
arrow.set_shape_cb = (callback_id, callback)
ax.add_patch(arrow)
return ax
@staticmethod
def _arrow_shape_cb(
arrow: FancyArrowPatch,
pos: tuple[float, float],
direction: ArrayLike,
sign: float,
bend_radius: float | None,
) -> Callable[[Event], None]:
from matplotlib import patches
def _cb(event: Event) -> None:
# We only want to set the shape once, so we disconnect ourselves
event.canvas.mpl_disconnect(arrow.set_shape_cb[0])
transform = arrow.axes.transData.transform
scale_x = transform((1, 0))[0] - transform((0, 0))[0]
scale_y = transform((0, 1))[1] - transform((0, 0))[1]
scale = max(scale_x, scale_y) # <-- Hack: This is a somewhat arbitrary choice.
arrow_length = ARROW_LENGTH * event.canvas.figure.get_dpi() / scale
if bend_radius:
v_norm = (direction[0] ** 2 + direction[1] ** 2) ** 0.5
vx_norm = direction[0] / v_norm
vy_norm = direction[1] / v_norm
bend_angle = -sign * arrow_length / bend_radius
t_x = 1 - np.cos(bend_angle)
t_y = np.sin(bend_angle)
v_x = -bend_radius * (vx_norm * t_y - vy_norm * t_x)
v_y = -bend_radius * (vx_norm * t_x + vy_norm * t_y)
tangent_angle = np.arctan2(direction[1], direction[0])
arrow.set_connectionstyle(
patches.ConnectionStyle.Angle3(
angleA=180 / np.pi * tangent_angle,
angleB=180 / np.pi * (tangent_angle + bend_angle),
)
)
else:
v_x = sign * arrow_length * direction[0]
v_y = sign * arrow_length * direction[1]
arrow.set_positions(pos, (pos[0] + v_x, pos[1] + v_y))
arrow.set_visible(True)
arrow.draw(event.renderer)
return _cb
def _volume(self, bounds: Bound) -> float:
"""Returns object's volume within given bounds."""
volume = 1
for axis in range(3):
min_bound = max(self.bounds[0][axis], bounds[0][axis])
max_bound = min(self.bounds[1][axis], bounds[1][axis])
volume *= max_bound - min_bound
return volume
def _surface_area(self, bounds: Bound) -> float:
"""Returns object's surface area within given bounds."""
min_bounds = list(self.bounds[0])
max_bounds = list(self.bounds[1])
in_bounds_factor = [2, 2, 2]
length = [0, 0, 0]
for axis in (0, 1, 2):
if min_bounds[axis] < bounds[0][axis]:
min_bounds[axis] = bounds[0][axis]
in_bounds_factor[axis] -= 1
if max_bounds[axis] > bounds[1][axis]:
max_bounds[axis] = bounds[1][axis]
in_bounds_factor[axis] -= 1
length[axis] = max_bounds[axis] - min_bounds[axis]
return (
length[0] * length[1] * in_bounds_factor[2]
+ length[1] * length[2] * in_bounds_factor[0]
+ length[2] * length[0] * in_bounds_factor[1]
)
""" Autograd code """
def _compute_derivatives(self, derivative_info: DerivativeInfo) -> AutogradFieldMap:
"""Compute the adjoint derivatives for this object."""
# get gradients w.r.t. each of the 6 faces (in normal direction)
vjps_faces = self._derivative_faces(derivative_info=derivative_info)
# post-process these values to give the gradients w.r.t. center and size
vjps_center_size = self._derivatives_center_size(vjps_faces=vjps_faces)
# store only the gradients asked for in 'field_paths'
derivative_map = {}
for field_path in derivative_info.paths:
field_name, *index = field_path
if field_name in vjps_center_size:
# if the vjp calls for a specific index into the tuple
if index and len(index) == 1:
index = int(index[0])
if field_path not in derivative_map:
derivative_map[field_path] = vjps_center_size[field_name][index]
# otherwise, just grab the whole array
else:
derivative_map[field_path] = vjps_center_size[field_name]
return derivative_map
@staticmethod
def _derivatives_center_size(vjps_faces: Bound) -> dict[str, Coordinate]:
"""Derivatives with respect to the ``center`` and ``size`` fields in the ``Box``."""
vjps_faces_min, vjps_faces_max = np.array(vjps_faces)
# post-process min and max face gradients into center and size
vjp_center = vjps_faces_max - vjps_faces_min
vjp_size = (vjps_faces_min + vjps_faces_max) / 2.0
return {
"center": tuple(vjp_center.tolist()),
"size": tuple(vjp_size.tolist()),
}
def _derivative_faces(self, derivative_info: DerivativeInfo) -> Bound:
"""Derivative with respect to normal position of 6 faces of ``Box``."""
axes_to_compute = (0, 1, 2)
if len(derivative_info.paths[0]) > 1:
axes_to_compute = tuple(info[1] for info in derivative_info.paths)
# change in permittivity between inside and outside
vjp_faces = np.zeros((2, 3))
for min_max_index, _ in enumerate((0, -1)):
for axis in axes_to_compute:
vjp_face = self._derivative_face(
min_max_index=min_max_index,
axis_normal=axis,
derivative_info=derivative_info,
)
# record vjp for this face
vjp_faces[min_max_index, axis] = vjp_face
return vjp_faces
def _derivative_face(
self,
min_max_index: int,
axis_normal: Axis,
derivative_info: DerivativeInfo,
) -> float:
"""Compute the derivative w.r.t. shifting a face in the normal direction."""
interpolators = derivative_info.interpolators or derivative_info.create_interpolators()
_, axis_perp = self.pop_axis((0, 1, 2), axis=axis_normal)
# First, check if the face is outside the simulation domain in which case set the
# face gradient to 0.
bounds_normal, _ = self.pop_axis(np.array(derivative_info.bounds).T, axis=axis_normal)
coord_normal_face = bounds_normal[min_max_index]
if min_max_index == 0:
if coord_normal_face < derivative_info.simulation_bounds[0][axis_normal]:
return 0.0
else:
if coord_normal_face > derivative_info.simulation_bounds[1][axis_normal]:
return 0.0
intersect_min, intersect_max = map(np.asarray, derivative_info.bounds_intersect)
extents = intersect_max - intersect_min
_, intersect_min_perp = self.pop_axis(np.array(intersect_min), axis=axis_normal)
_, intersect_max_perp = self.pop_axis(np.array(intersect_max), axis=axis_normal)
is_2d_map = []
for axis_idx in range(3):
if axis_idx == axis_normal:
continue
is_2d_map.append(np.isclose(extents[axis_idx], 0.0))
if np.all(is_2d_map):
return 0.0
is_2d = np.any(is_2d_map)
# Build point grid
adaptive_spacing = derivative_info.adaptive_vjp_spacing()
def spacing_to_grid_points(
spacing: float, min_coord: float, max_coord: float
) -> NDArray[float]:
N = np.maximum(3, 1 + int((max_coord - min_coord) / spacing))
points = np.linspace(min_coord, max_coord, N)
centers = 0.5 * (points[0:-1] + points[1:])
return centers
def verify_integration_interval(bound: tuple[float, float]) -> bool:
# assume the bounds should not be equal or else this integration interval
# would be the flat dimension of a 2D geometry.
return bound[1] > bound[0]
def compute_integration_weight(grid_points: NDArray[float]) -> float:
grid_spacing = grid_points[1] - grid_points[0]
if grid_spacing == 0.0:
integration_weight = 1.0 / len(grid_points)
else:
integration_weight = grid_points[1] - grid_points[0]
return integration_weight
if is_2d:
# build 1D grid for sampling points along the face, which is an edge in the 2D case
zero_dim = np.where(is_2d_map)[0][0]
# zero dim is one of the perpendicular directions, so the other perpendicular direction
# is the nonzero dimension
nonzero_dim = 1 - zero_dim
# clip at simulation bounds for integration dimension
integration_bounds_perp = (
intersect_min_perp[nonzero_dim],
intersect_max_perp[nonzero_dim],
)
if not verify_integration_interval(integration_bounds_perp):
return 0.0
grid_points_linear = spacing_to_grid_points(
adaptive_spacing, integration_bounds_perp[0], integration_bounds_perp[1]
)
integration_weight = compute_integration_weight(grid_points_linear)
grid_points = np.repeat(np.expand_dims(grid_points_linear.copy(), 1), 3, axis=1)
# set up grid points to pass into evaluate_gradient_at_points
grid_points[:, axis_perp[nonzero_dim]] = grid_points_linear
grid_points[:, axis_perp[zero_dim]] = intersect_min_perp[zero_dim]
grid_points[:, axis_normal] = coord_normal_face
else:
# build 3D grid for sampling points along the face
# clip at simulation bounds for each integration dimension
integration_bounds_perp = (
(intersect_min_perp[0], intersect_max_perp[0]),
(intersect_min_perp[1], intersect_max_perp[1]),
)
if not np.all([verify_integration_interval(b) for b in integration_bounds_perp]):
return 0.0
grid_points_perp_1 = spacing_to_grid_points(
adaptive_spacing, integration_bounds_perp[0][0], integration_bounds_perp[0][1]
)
grid_points_perp_2 = spacing_to_grid_points(
adaptive_spacing, integration_bounds_perp[1][0], integration_bounds_perp[1][1]
)
integration_weight = compute_integration_weight(
grid_points_perp_1
) * compute_integration_weight(grid_points_perp_2)
mesh_perp1, mesh_perp2 = np.meshgrid(grid_points_perp_1, grid_points_perp_2)
zip_perp_coords = np.array(list(zip(mesh_perp1.flatten(), mesh_perp2.flatten())))
grid_points = np.pad(zip_perp_coords.copy(), ((0, 0), (1, 0)), mode="constant")
# set up grid points to pass into evaluate_gradient_at_points
grid_points[:, axis_perp[0]] = zip_perp_coords[:, 0]
grid_points[:, axis_perp[1]] = zip_perp_coords[:, 1]
grid_points[:, axis_normal] = coord_normal_face
normals = np.zeros_like(grid_points)
perps1 = np.zeros_like(grid_points)
perps2 = np.zeros_like(grid_points)
normals[:, axis_normal] = -1 if (min_max_index == 0) else 1
perps1[:, axis_perp[0]] = 1
perps2[:, axis_perp[1]] = 1
gradient_at_points = derivative_info.evaluate_gradient_at_points(
spatial_coords=grid_points,
normals=normals,
perps1=perps1,
perps2=perps2,
interpolators=interpolators,
)
vjp_value = np.sum(integration_weight * np.real(gradient_at_points))
return vjp_value
"""Compound subclasses"""
[docs]
class ClipOperation(Geometry):
"""Class representing the result of a set operation between geometries."""
operation: ClipOperationType = Field(
title="Operation Type",
description="Operation to be performed between geometries.",
)
geometry_a: discriminated_union(GeometryType) = Field(
title="Geometry A",
description="First operand for the set operation. It can be any geometry type, including "
":class:`GeometryGroup`.",
)
geometry_b: discriminated_union(GeometryType) = Field(
title="Geometry B",
description="Second operand for the set operation. It can also be any geometry type.",
)
[docs]
@staticmethod
def to_polygon_list(base_geometry: Shapely, cleanup: bool = False) -> list[Shapely]:
"""Return a list of valid polygons from a shapely geometry, discarding points, lines, and
empty polygons, and empty triangles within polygons.
Parameters
----------
base_geometry : shapely.geometry.base.BaseGeometry
Base geometry for inspection.
cleanup: bool = False
If True, removes extremely small features from each polygon's boundary.
This is useful for removing artifacts from 2D plots displayed to the user.
Returns
-------
list[shapely.geometry.base.BaseGeometry]
Valid polygons retrieved from ``base geometry``.
"""
unfiltered_geoms = []
if base_geometry.geom_type == "GeometryCollection":
unfiltered_geoms = [
p
for geom in base_geometry.geoms
for p in ClipOperation.to_polygon_list(geom, cleanup)
]
if base_geometry.geom_type == "MultiPolygon":
unfiltered_geoms = [p for p in base_geometry.geoms if not p.is_empty]
if base_geometry.geom_type == "Polygon" and not base_geometry.is_empty:
unfiltered_geoms = [base_geometry]
geoms = []
if cleanup:
# Optional: "clean" each of the polygons (by removing extremely small or thin features).
for geom in unfiltered_geoms:
geom_clean = cleanup_shapely_object(geom)
if geom_clean.geom_type == "Polygon":
geoms.append(geom_clean)
if geom_clean.geom_type == "MultiPolygon":
geoms += [p for p in geom_clean.geoms if not p.is_empty]
# Ignore other types of shapely objects (points and lines)
else:
geoms = unfiltered_geoms
return geoms
@property
def _shapely_operation(self) -> Callable[[Shapely, Shapely], Shapely]:
"""Return a Shapely function equivalent to this operation."""
result = _shapely_operations.get(self.operation, None)
if not result:
raise ValueError(
"'operation' must be one of 'union', 'intersection', 'difference', or "
"'symmetric_difference'."
)
return result
@property
def _bit_operation(self) -> Callable[[Any, Any], Any]:
"""Return a function equivalent to this operation using bit operators."""
result = _bit_operations.get(self.operation, None)
if not result:
raise ValueError(
"'operation' must be one of 'union', 'intersection', 'difference', or "
"'symmetric_difference'."
)
return result
[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. If ``None``, uses
high-quality visualization settings.
Returns
-------
list[shapely.geometry.base.BaseGeometry]
List of 2D shapes that intersect plane.
For more details refer to
`Shapely's Documentation <https://shapely.readthedocs.io/en/stable/project.html>`_.
"""
a = self.geometry_a.intersections_tilted_plane(
normal,
origin,
to_2D,
cleanup=cleanup,
quad_segs=quad_segs,
section_tolerance_2d=section_tolerance_2d,
)
b = self.geometry_b.intersections_tilted_plane(
normal,
origin,
to_2D,
cleanup=cleanup,
quad_segs=quad_segs,
section_tolerance_2d=section_tolerance_2d,
)
geom_a = shapely.unary_union([Geometry.evaluate_inf_shape(g) for g in a])
geom_b = shapely.unary_union([Geometry.evaluate_inf_shape(g) for g in b])
return ClipOperation.to_polygon_list(
self._shapely_operation(geom_a, geom_b),
cleanup=cleanup,
)
[docs]
def intersections_plane(
self,
x: float | None = None,
y: float | None = None,
z: float | None = None,
cleanup: bool = True,
quad_segs: int | None = None,
section_tolerance_2d: bool = False,
) -> list[Shapely]:
"""Returns list of shapely geometries at plane specified by one non-None value of x,y,z.
Parameters
----------
x : float = None
Position of plane in x direction, only one of x,y,z can be specified to define plane.
y : float = None
Position of plane in y direction, only one of x,y,z can be specified to define plane.
z : float = None
Position of plane in z direction, only one of x,y,z can be specified to define plane.
cleanup : bool = True
If True, removes extremely small features from each polygon's boundary.
quad_segs : Optional[int] = None
Number of segments used to discretize circular shapes. If ``None``, uses
high-quality visualization settings.
section_tolerance_2d : bool = False
See :meth:`Geometry.intersections_plane`.
Returns
-------
list[shapely.geometry.base.BaseGeometry]
List of 2D shapes that intersect plane.
For more details refer to
`Shapely's Documentaton <https://shapely.readthedocs.io/en/stable/project.html>`_.
"""
a = self.geometry_a.intersections_plane(
x,
y,
z,
cleanup=cleanup,
quad_segs=quad_segs,
section_tolerance_2d=section_tolerance_2d,
)
b = self.geometry_b.intersections_plane(
x,
y,
z,
cleanup=cleanup,
quad_segs=quad_segs,
section_tolerance_2d=section_tolerance_2d,
)
geom_a = shapely.unary_union([Geometry.evaluate_inf_shape(g) for g in a])
geom_b = shapely.unary_union([Geometry.evaluate_inf_shape(g) for g in b])
return ClipOperation.to_polygon_list(
self._shapely_operation(geom_a, geom_b),
cleanup=cleanup,
)
@cached_property
def bounds(self) -> Bound:
"""Returns bounding box min and max coordinates.
Returns
-------
tuple[float, float, float], tuple[float, float float]
Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
"""
# Overestimates
if self.operation == "difference":
result = self.geometry_a.bounds
elif self.operation == "intersection":
bounds = (self.geometry_a.bounds, self.geometry_b.bounds)
result = (
tuple(max(b[i] for b, _ in bounds) for i in range(3)),
tuple(min(b[i] for _, b in bounds) for i in range(3)),
)
if any(result[0][i] > result[1][i] for i in range(3)):
result = ((0, 0, 0), (0, 0, 0))
else:
bounds = (self.geometry_a.bounds, self.geometry_b.bounds)
result = (
tuple(min(b[i] for b, _ in bounds) for i in range(3)),
tuple(max(b[i] for _, b in bounds) for i in range(3)),
)
return result
[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.
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.
"""
inside_a = self.geometry_a.inside(x, y, z)
inside_b = self.geometry_b.inside(x, y, z)
return self._bit_operation(inside_a, inside_b)
[docs]
def inside_meshgrid(
self, x: NDArray[float], y: NDArray[float], z: NDArray[float]
) -> NDArray[bool]:
"""Faster way to check ``self.inside`` on a meshgrid. The input arrays are assumed sorted.
Parameters
----------
x : np.ndarray[float]
1D array of point positions in x direction.
y : np.ndarray[float]
1D array of point positions in y direction.
z : np.ndarray[float]
1D array of point positions in z direction.
Returns
-------
np.ndarray[bool]
Array with shape ``(x.size, y.size, z.size)``, which is ``True`` for every
point that is inside the geometry.
"""
inside_a = self.geometry_a.inside_meshgrid(x, y, z)
inside_b = self.geometry_b.inside_meshgrid(x, y, z)
return self._bit_operation(inside_a, inside_b)
def _volume(self, bounds: Bound) -> float:
"""Returns object's volume within given bounds."""
# Overestimates
if self.operation == "intersection":
return min(self.geometry_a.volume(bounds), self.geometry_b.volume(bounds))
if self.operation == "difference":
return self.geometry_a.volume(bounds)
return self.geometry_a.volume(bounds) + self.geometry_b.volume(bounds)
def _surface_area(self, bounds: Bound) -> float:
"""Returns object's surface area within given bounds."""
# Overestimates
return self.geometry_a.surface_area(bounds) + self.geometry_b.surface_area(bounds)
@cached_property
def _normal_2dmaterial(self) -> Axis:
"""Get the normal to the given geometry, checking that it is a 2D geometry."""
normal_a = self.geometry_a._normal_2dmaterial
normal_b = self.geometry_b._normal_2dmaterial
if normal_a != normal_b:
raise ValidationError(
"'Medium2D' requires both geometries in the 'ClipOperation' to "
"have exactly one dimension with zero size in common."
)
plane_position_a = self.geometry_a.bounds[0][normal_a]
plane_position_b = self.geometry_b.bounds[0][normal_b]
if plane_position_a != plane_position_b:
raise ValidationError(
"'Medium2D' requires both geometries in the 'ClipOperation' to be co-planar."
)
return normal_a
def _update_from_bounds(self, bounds: tuple[float, float], axis: Axis) -> ClipOperation:
"""Returns an updated geometry which has been transformed to fit within ``bounds``
along the ``axis`` direction."""
new_geom_a = self.geometry_a._update_from_bounds(bounds=bounds, axis=axis)
new_geom_b = self.geometry_b._update_from_bounds(bounds=bounds, axis=axis)
return self.updated_copy(geometry_a=new_geom_a, geometry_b=new_geom_b)
def _resolve_autograd_route(self, field_path: tuple[Any, ...]) -> AutogradRoute:
"""Resolve and validate one traced ClipOperation path for adjoint routing."""
return resolve_delegated_autograd_route(
parameter_kind="geometry",
owner_kind="geometry type",
owner_name=type(self).__name__,
field_path=field_path,
delegates={"geometry_a": self.geometry_a, "geometry_b": self.geometry_b},
supported_parameters=(
"geometry_a.<parameter>",
"geometry_b.<parameter>",
),
)
def _compute_derivatives(self, derivative_info: DerivativeInfo) -> AutogradFieldMap:
"""Compute adjoint derivatives by accumulating contributions from both operands."""
geometry_paths = {"geometry_a": [], "geometry_b": []}
for path in derivative_info.paths:
geometry_key, *sub_path = path
geometry_paths[geometry_key].append(tuple(sub_path))
if derivative_info.clipped_geometry is None:
raise ValidationError(
"ClipOperation derivative evaluation requires `clipped_geometry`."
)
geometry_map = {
"geometry_a": self.geometry_a,
"geometry_b": self.geometry_b,
}
grad_vjps = {}
# Reuse interpolation data for both operands to avoid duplicate setup.
interpolators = derivative_info.interpolators or derivative_info.create_interpolators()
with derivative_info.cache_min_spacing_from_permittivity():
for geometry_key, geometry in geometry_map.items():
paths = geometry_paths[geometry_key]
if not paths:
continue
geometry_info = derivative_info.updated_copy(
paths=paths,
bounds=geometry.bounds,
bounds_intersect=self.bounds_intersection(
geometry.bounds, derivative_info.simulation_bounds
),
deep=False,
interpolators=interpolators,
)
vjp_dict_geometry = geometry._compute_derivatives(geometry_info)
for geo_path, geo_vjp in vjp_dict_geometry.items():
full_path = (geometry_key, *geo_path)
if full_path in grad_vjps:
existing = grad_vjps[full_path]
if isinstance(existing, (list, tuple)) and isinstance(
geo_vjp, (list, tuple)
):
grad_vjps[full_path] = type(existing)(
x + y for x, y in zip(existing, geo_vjp)
)
else:
grad_vjps[full_path] = existing + geo_vjp
else:
grad_vjps[full_path] = geo_vjp
return grad_vjps
[docs]
class GeometryGroup(Geometry):
"""A collection of Geometry objects that can be called as a single geometry object."""
geometries: tuple[discriminated_union(GeometryType), ...] = Field(
title="Geometries",
description="Tuple of geometries in a single grouping. "
"Can provide significant performance enhancement in ``Structure`` when all geometries are "
"assigned the same medium.",
)
@field_validator("geometries")
@classmethod
def _geometries_not_empty(cls, val: tuple[GeometryType, ...]) -> tuple[GeometryType, ...]:
"""make sure geometries are not empty."""
if not len(val) > 0:
raise ValidationError("GeometryGroup.geometries must not be empty.")
return val
@cached_property
def bounds(self) -> Bound:
"""Returns bounding box min and max coordinates.
Returns
-------
tuple[float, float, float], tuple[float, float, float]
Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
"""
bounds = tuple(geometry.bounds for geometry in self.geometries)
return (
tuple(min(b[i] for b, _ in bounds) for i in range(3)),
tuple(max(b[i] for _, b in bounds) for i in range(3)),
)
[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. If ``None``, uses
high-quality visualization settings.
Returns
-------
list[shapely.geometry.base.BaseGeometry]
List of 2D shapes that intersect plane.
For more details refer to
`Shapely's Documentation <https://shapely.readthedocs.io/en/stable/project.html>`_.
"""
return [
intersection
for geometry in self.geometries
for intersection in geometry.intersections_tilted_plane(
normal,
origin,
to_2D,
cleanup=cleanup,
quad_segs=quad_segs,
section_tolerance_2d=section_tolerance_2d,
)
]
[docs]
def intersections_plane(
self,
x: float | None = None,
y: float | None = None,
z: float | None = None,
cleanup: bool = True,
quad_segs: int | None = None,
section_tolerance_2d: bool = False,
) -> list[Shapely]:
"""Returns list of shapely geometries at plane specified by one non-None value of x,y,z.
Parameters
----------
x : float = None
Position of plane in x direction, only one of x,y,z can be specified to define plane.
y : float = None
Position of plane in y direction, only one of x,y,z can be specified to define plane.
z : float = None
Position of plane in z direction, only one of x,y,z can be specified to define plane.
cleanup : bool = True
If True, removes extremely small features from each polygon's boundary.
quad_segs : Optional[int] = None
Number of segments used to discretize circular shapes. If ``None``, uses
high-quality visualization settings.
Returns
-------
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>`_.
"""
axis, position = self.parse_xyz_kwargs(x=x, y=y, z=z)
return [
intersection
for geometry in self.geometries
if geometry.intersects_axis_position(
axis, position, section_tolerance_2d=section_tolerance_2d
)
for intersection in geometry.intersections_plane(
x=x,
y=y,
z=z,
cleanup=cleanup,
quad_segs=quad_segs,
section_tolerance_2d=section_tolerance_2d,
)
]
[docs]
def intersects_axis_position(
self, axis: float, position: float, section_tolerance_2d: bool = False
) -> bool:
"""Whether self intersects plane specified by a given position along a normal axis.
Parameters
----------
axis : int = None
Axis normal to the plane.
position : float = None
Position of plane along the normal axis.
section_tolerance_2d : bool = False
See :meth:`Geometry.intersects_axis_position`.
Returns
-------
bool
Whether this geometry intersects the plane.
"""
return any(
geom.intersects_axis_position(axis, position, section_tolerance_2d=section_tolerance_2d)
for geom in self.geometries
)
[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.
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.
"""
individual_insides = (geometry.inside(x, y, z) for geometry in self.geometries)
return functools.reduce(lambda a, b: a | b, individual_insides)
[docs]
def inside_meshgrid(
self, x: NDArray[float], y: NDArray[float], z: NDArray[float]
) -> NDArray[bool]:
"""Faster way to check ``self.inside`` on a meshgrid. The input arrays are assumed sorted.
Parameters
----------
x : np.ndarray[float]
1D array of point positions in x direction.
y : np.ndarray[float]
1D array of point positions in y direction.
z : np.ndarray[float]
1D array of point positions in z direction.
Returns
-------
np.ndarray[bool]
Array with shape ``(x.size, y.size, z.size)``, which is ``True`` for every
point that is inside the geometry.
"""
individual_insides = (geom.inside_meshgrid(x, y, z) for geom in self.geometries)
return functools.reduce(lambda a, b: a | b, individual_insides)
def _volume(self, bounds: Bound) -> float:
"""Returns object's volume within given bounds."""
return sum(geometry.volume(bounds) for geometry in self.geometries)
def _surface_area(self, bounds: Bound) -> float:
"""Returns object's surface area within given bounds."""
return sum(geometry.surface_area(bounds) for geometry in self.geometries)
@cached_property
def _normal_2dmaterial(self) -> Axis:
"""Get the normal to the given geometry, checking that it is a 2D geometry."""
normals = {geom._normal_2dmaterial for geom in self.geometries}
if len(normals) != 1:
raise ValidationError(
"'Medium2D' requires all geometries in the 'GeometryGroup' to "
"share exactly one dimension with zero size."
)
normal = list(normals)[0]
positions = {geom.bounds[0][normal] for geom in self.geometries}
if len(positions) != 1:
raise ValidationError(
"'Medium2D' requires all geometries in the 'GeometryGroup' to be co-planar."
)
return normal
def _update_from_bounds(self, bounds: tuple[float, float], axis: Axis) -> GeometryGroup:
"""Returns an updated geometry which has been transformed to fit within ``bounds``
along the ``axis`` direction."""
new_geometries = tuple(
geometry._update_from_bounds(bounds=bounds, axis=axis) for geometry in self.geometries
)
return self.updated_copy(geometries=new_geometries)
def _resolve_autograd_route(self, field_path: tuple[Any, ...]) -> AutogradRoute:
"""Resolve and validate one traced GeometryGroup path for adjoint routing."""
if len(field_path) < 2 or field_path[0] != "geometries":
_raise_unsupported_traced_geometry_path(
type(self).__name__,
field_path,
supported_parameters=("geometries[index].<parameter>",),
)
index = field_path[1]
sub_path = field_path[2:]
try:
self.geometries[index]._resolve_autograd_route(sub_path)
except AdjointError as err:
raise_with_traced_path_context(
err,
parameter_kind="geometry",
local_path=sub_path,
full_path=field_path,
)
return AutogradRoute(local_path=field_path)
def _compute_derivatives(self, derivative_info: DerivativeInfo) -> AutogradFieldMap:
"""Compute the adjoint derivatives for this object."""
grad_vjps = {}
# create interpolators once for all geometries to avoid redundant field data conversions
interpolators = derivative_info.interpolators or derivative_info.create_interpolators()
with derivative_info.cache_min_spacing_from_permittivity():
for field_path in derivative_info.paths:
_, index, *geo_path = field_path
geo = self.geometries[index]
# pass pre-computed interpolators if available
geo_info = derivative_info.updated_copy(
paths=[tuple(geo_path)],
bounds=geo.bounds,
bounds_intersect=self.bounds_intersection(
geo.bounds, derivative_info.simulation_bounds
),
deep=False,
interpolators=interpolators,
)
vjp_dict_geo = geo._compute_derivatives(geo_info)
if len(vjp_dict_geo) != 1:
raise AssertionError("Got multiple gradients for single geometry field.")
grad_vjps[field_path] = vjp_dict_geo.popitem()[1]
return grad_vjps
[docs]
class GeometryArray(Geometry):
"""A geometry representing an array of copies of a base geometry, with optional offsets
and/or linear transformations applied to each copy.
Notes
-----
This class provides an efficient way to represent arrays of repeated geometries,
avoiding the need to create many individual geometry objects.
The instance pose for each copy is defined as: ``T(offsets[i]) @ L(transforms[i])``,
where ``T`` is a translation matrix and ``L`` is the linear transform. In other words,
the transform is applied first, then the translation.
- ``offsets`` represent all per-instance translation.
- ``transforms`` represent linear transforms only (rotation/reflection/scale/shear)
and must not contain translation. Use ``offsets`` for translations.
- If both ``offsets`` and ``transforms`` are ``None``, the array contains a single
instance of the base geometry at the origin.
- If both are provided, they must have the same length.
- Adjoint/autodiff is not currently supported for ``GeometryArray``.
Example
-------
>>> import tidy3d as td
>>> import numpy as np
>>> box = td.Box(size=(1, 1, 1))
>>> # Using offsets only:
>>> offsets = [[0, 0, 0], [2, 0, 0], [0, 2, 0], [2, 2, 0]]
>>> array = td.GeometryArray(geometry=box, offsets=offsets)
>>> # Or use the convenience method:
>>> array = box.array(offsets=offsets)
>>> # Using linear transforms only (rotation around z-axis):
>>> rot_0 = td.Transformed.rotation(0, 2) # no rotation
>>> rot_90 = td.Transformed.rotation(np.pi/2, 2) # 90 degree rotation
>>> array = td.GeometryArray(geometry=box, transforms=[rot_0, rot_90])
>>> # Both None gives single instance of base geometry:
>>> array = td.GeometryArray(geometry=box)
"""
geometry: discriminated_union(GeometryType) = Field(
...,
title="Geometry",
description="Base geometry to be repeated in the array.",
)
offsets: tuple[Coordinate, ...] | None = Field(
None,
title="Offsets",
description="A tuple of 3D coordinate offsets. Each offset translates the base "
"geometry (after any transform is applied) to create a copy. If not provided, no "
"additional translation is applied beyond any transforms.",
)
transforms: tuple[MatrixReal4x4, ...] | None = Field(
None,
title="Transforms",
description="A tuple of 4x4 linear-only transformation matrices "
"(rotation/reflection/scale/shear, no translation). Typical transforms can be "
"created using ``Transformed.rotation``, ``Transformed.reflection``, or ``Transformed.scaling``. "
"Each transform is applied to the base geometry before the corresponding offset translation. "
"If not provided, only translations from offsets are applied.",
)
_geometry_is_finite = assert_geometry_finite("geometry")
@field_validator("transforms")
@classmethod
def _validate_transforms(
cls, val: tuple[MatrixReal4x4, ...] | None
) -> tuple[MatrixReal4x4, ...] | None:
"""Validate that transforms are invertible, linear-only, and non-empty if provided."""
if val is None:
return val
# Must not be empty if provided
if len(val) < 1:
raise ValidationError("'transforms' must have at least one transform when provided.")
# Check each transform
for i, transform in enumerate(val):
# Check invertibility
check_transform_invertible(transform, index=i)
# Check linear-only (no translation)
transform_array = np.asarray(transform)
# Check translation column: transform[:3, 3] should be zero
translation = transform_array[:3, 3]
if not np.allclose(translation, 0, atol=LINEAR_TRANSFORM_TOL):
idx_msg = f"at index {i}"
raise ValidationError(
f"Transform {idx_msg} contains translation in [:3, 3] = {translation.tolist()}. "
"GeometryArray transforms must be linear-only (rotation/reflection/scale/shear). "
"Use the 'offsets' parameter for translations."
)
# Check bottom row: transform[3, :] should be [0, 0, 0, 1]
bottom_row = transform_array[3, :]
expected_bottom = np.array([0, 0, 0, 1])
if not np.allclose(bottom_row, expected_bottom, atol=LINEAR_TRANSFORM_TOL):
idx_msg = f"at index {i}"
raise ValidationError(
f"Transform {idx_msg} has invalid homogeneous form: [3, :] = {bottom_row.tolist()}. "
"Expected [0, 0, 0, 1]."
)
return val
@model_validator(mode="after")
def _validate_offsets_and_transforms(self) -> Self:
"""Validate offsets and transforms are consistent."""
offsets = self.offsets
transforms = self.transforms
# If offsets provided, must not be empty
if offsets is not None and len(offsets) < 1:
self._raise_validation_error_at_loc(
ValidationError("'offsets' must have at least one offset when provided."), "offsets"
)
# If both provided, lengths must match
if offsets is not None and transforms is not None:
if len(offsets) != len(transforms):
self._raise_validation_error_at_loc(
ValidationError(
f"Number of transforms ({len(transforms)}) must match "
f"number of offsets ({len(offsets)}) when both are provided."
),
"transforms",
)
return self
@cached_property
def num_geometries(self) -> int:
"""Number of geometries in the array."""
if self.offsets is not None:
return len(self.offsets)
if self.transforms is not None:
return len(self.transforms)
# Both None means single geometry (base geometry at origin)
return 1
@cached_property
def _all_transforms(self) -> np.ndarray:
"""Compute all 4x4 transforms for all geometries in a vectorized way.
Returns
-------
numpy.ndarray
Array of shape (num_geometries, 4, 4) containing the full transform
(rotation/scale + translation) for each geometry in the array.
"""
n = self.num_geometries
shape = (n, 4, 4)
# Get all transforms, defaulting to identity if not provided
if self.transforms is not None:
transforms = np.array(self.transforms)
else:
transforms = np.broadcast_to(Transformed.identity(), shape)
# Build translation matrices for all offsets
# translation matrix: [[1,0,0,x], [0,1,0,y], [0,0,1,z], [0,0,0,1]]
translations = np.broadcast_to(Transformed.identity(), shape).copy()
if self.offsets is not None:
translations[:, :3, 3] = self.offsets
# Apply transform, then translation: result = translation @ transform
return np.matmul(translations, transforms)
def _get_full_transform(self, index: int) -> MatrixReal4x4:
"""Get the full 4x4 transform for a geometry at given index (transform + translation)."""
return self._all_transforms[index]
@cached_property
def _transformed_geometries(self) -> list[Transformed]:
"""List of transformed geometries in the array."""
return [
Transformed(geometry=self.geometry, transform=transform)
for transform in self._all_transforms
]
@cached_property
def _geometry_group(self) -> GeometryGroup:
"""Return a GeometryGroup containing all transformed geometries in the array."""
return GeometryGroup(geometries=tuple(self._transformed_geometries))
@cached_property
def bounds(self) -> Bound:
"""Returns bounding box min and max coordinates.
Returns
-------
Tuple[float, float, float], Tuple[float, float, float]
Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
"""
return self._geometry_group.bounds
[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. If ``None``, uses
high-quality visualization settings.
section_tolerance_2d : bool = False
See :meth:`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 self._geometry_group.intersections_tilted_plane(
normal,
origin,
to_2D,
cleanup=cleanup,
quad_segs=quad_segs,
section_tolerance_2d=section_tolerance_2d,
)
[docs]
def intersections_plane(
self,
x: float | None = None,
y: float | None = None,
z: float | None = None,
cleanup: bool = True,
quad_segs: int | None = None,
section_tolerance_2d: bool = False,
) -> list[Shapely]:
"""Returns list of shapely geometries at plane specified by one non-None value of x,y,z.
Parameters
----------
x : float = None
Position of plane in x direction, only one of x,y,z can be specified to define plane.
y : float = None
Position of plane in y direction, only one of x,y,z can be specified to define plane.
z : float = None
Position of plane in z direction, only one of x,y,z can be specified to define plane.
cleanup : bool = True
If True, removes extremely small features from each polygon's boundary.
quad_segs : Optional[int] = None
Number of segments used to discretize circular shapes. If ``None``, uses
high-quality visualization settings.
section_tolerance_2d : bool = False
See :meth:`Geometry.intersections_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 self._geometry_group.intersections_plane(
x=x,
y=y,
z=z,
cleanup=cleanup,
quad_segs=quad_segs,
section_tolerance_2d=section_tolerance_2d,
)
[docs]
def intersects_axis_position(
self, axis: int, position: float, section_tolerance_2d: bool = False
) -> bool:
"""Whether self intersects plane specified by a given position along a normal axis.
Parameters
----------
axis : int = None
Axis normal to the plane.
position : float = None
Position of plane along the normal axis.
section_tolerance_2d : bool = False
See :meth:`Geometry.intersects_axis_position`.
Returns
-------
bool
Whether this geometry intersects the plane.
"""
return self._geometry_group.intersects_axis_position(
axis, position, section_tolerance_2d=section_tolerance_2d
)
[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:`Geometry`, and ``False`` otherwise.
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.
"""
return self._geometry_group.inside(x, y, z)
def _volume(self, bounds: Bound) -> float:
"""Returns object's volume within given bounds."""
return self._geometry_group._volume(bounds)
def _surface_area(self, bounds: Bound) -> float:
"""Returns object's surface area within given bounds."""
# Surface area cannot be reliably computed when non-trivial transforms are present
if self.transforms is not None:
log.warning("Surface area of transformed elements cannot be calculated.")
return None
# For pure translations, sum surface areas using local bounds for base geometry
total_area = 0.0
for geom in self._transformed_geometries:
# Transform bounds to local coordinate system
vertices = np.dot(geom.inverse, Transformed._vertices_from_bounds(bounds))[:3]
local_bounds = (tuple(vertices.min(axis=1)), tuple(vertices.max(axis=1)))
instance_area = self.geometry.surface_area(local_bounds)
if instance_area is None:
return None
total_area += instance_area
return total_area
@cached_property
def _normal_2dmaterial(self) -> Axis:
"""Get the normal to the given geometry, checking that it is a 2D geometry."""
return self._geometry_group._normal_2dmaterial
def _update_from_bounds(self, bounds: tuple[float, float], axis: Axis) -> GeometryGroup:
"""Returns an updated geometry which has been transformed to fit within ``bounds``
along the ``axis`` direction."""
return self._geometry_group._update_from_bounds(bounds=bounds, axis=axis)
def _compute_derivatives(self, derivative_info: DerivativeInfo) -> AutogradFieldMap:
"""Compute the adjoint derivatives for this object.
Raises
------
NotImplementedError
Adjoint/autodiff is not currently supported for GeometryArray.
"""
raise NotImplementedError(
"Adjoint is not currently supported for 'GeometryArray'.",
)
def cleanup_shapely_object(obj: Shapely, tolerance_ratio: float = POLY_TOLERANCE_RATIO) -> Shapely:
"""Remove small geometric features from the boundaries of a shapely object including
inward and outward spikes, thin holes, and thin connections between larger regions.
Parameters
----------
obj : shapely
a shapely object (typically a ``Polygon`` or a ``MultiPolygon``)
tolerance_ratio : float = ``POLY_TOLERANCE_RATIO``
Features on the boundaries of polygons will be discarded if they are smaller
or narrower than ``tolerance_ratio`` multiplied by the size of the object.
Returns
-------
Shapely
A new shapely object whose small features (eg. thin spikes or holes) are removed.
Notes
-----
This function does not attempt to delete overlapping, nearby, or collinear vertices.
To solve that problem, use ``shapely.simplify()`` afterwards.
"""
if _package_is_older_than("shapely", "2.1"):
log.warning("Versions of shapely prior to v2.1 may cause plot errors.", log_once=True)
return obj
if obj.is_empty:
return obj
centroid = obj.centroid
object_size = min(obj.bounds[2] - obj.bounds[0], obj.bounds[3] - obj.bounds[1])
if object_size == 0.0:
return shapely.Polygon([])
# To prevent numerical over- or underflow errors, subtract the centroid and rescale
normalized_obj = shapely.affinity.affine_transform(
obj,
matrix=[
1 / object_size,
0.0,
0.0,
1 / object_size,
-centroid.x / object_size,
-centroid.y / object_size,
],
)
# Important: Remove any self intersections beforehand using `shapely.make_valid()`.
valid_obj = shapely.make_valid(normalized_obj, method="structure", keep_collapsed=False)
# To get rid of small thin features, erode(shrink), dilate(expand), and erode again.
eroded_obj = shapely.buffer(
valid_obj,
distance=-tolerance_ratio,
cap_style="square",
quad_segs=3,
)
dilated_obj = shapely.buffer(
eroded_obj,
distance=2 * tolerance_ratio,
cap_style="square",
quad_segs=3,
)
cleaned_obj = dilated_obj
# Optional: Now shrink the polygon back to the original size.
cleaned_obj = shapely.buffer(
cleaned_obj,
distance=-tolerance_ratio,
cap_style="square",
quad_segs=3,
)
# Clean vertices of very close distances created during the erosion/dilation process.
# The distance value is heuristic.
cleaned_obj = cleaned_obj.simplify(POLY_DISTANCE_TOLERANCE, preserve_topology=True)
# Revert to the original scale and position.
rescaled_clean_obj = shapely.affinity.affine_transform(
cleaned_obj,
matrix=[
object_size,
0.0,
0.0,
object_size,
centroid.x,
centroid.y,
],
)
return rescaled_clean_obj
from .utils import GeometryType, from_shapely, vertices_from_shapely # noqa: E402
