"""Abstract base classes for geometry."""
from __future__ import annotations
import functools
import pathlib
from abc import ABC, abstractmethod
from collections.abc import Iterable, Sequence
from os import PathLike
from typing import TYPE_CHECKING, Any, Callable, Optional, Union
import autograd.numpy as np
import pydantic.v1 as pydantic
import shapely
from numpy._typing import ArrayLike, NDArray
from typing_extensions import Self
try:
from matplotlib import patches
except ImportError:
pass
from tidy3d.compat import _shapely_is_older_than
from tidy3d.components.autograd import (
AutogradFieldMap,
TracedCoordinate,
TracedFloat,
TracedSize,
get_static,
)
from tidy3d.components.autograd.derivative_utils import DerivativeInfo
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 (
ArrayFloat2D,
ArrayFloat3D,
Ax,
Axis,
Bound,
ClipOperationType,
Coordinate,
Coordinate2D,
LengthUnit,
MatrixReal4x4,
PlanePosition,
Shapely,
Size,
annotate_type,
)
from tidy3d.components.viz import (
ARROW_LENGTH,
PLOT_BUFFER,
PlotParams,
VisualizationSpec,
add_ax_if_none,
arrow_style,
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 (
SetupError,
Tidy3dError,
Tidy3dImportError,
Tidy3dKeyError,
ValidationError,
)
from tidy3d.log import log
from tidy3d.packaging import verify_packages_import
if TYPE_CHECKING:
from gdstk import Cell
from matplotlib.backend_bases import Event
from matplotlib.patches import FancyArrowPatch
POLY_GRID_SIZE = 1e-12
POLY_TOLERANCE_RATIO = 1e-12
POLY_DISTANCE_TOLERANCE = 8e-12
_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,
}
[docs]
class Geometry(Tidy3dBaseModel, ABC):
"""Abstract base class, defines where something exists in space."""
@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:`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: Optional[int] = 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.
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>`_.
"""
[docs]
def intersections_plane(
self,
x: Optional[float] = None,
y: Optional[float] = None,
z: Optional[float] = None,
cleanup: bool = True,
quad_segs: Optional[int] = None,
) -> list[Shapely]:
"""Returns list of shapely geometries at plane specified by one non-None value of x,y,z.
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)
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
)
[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:`Geometry` have intersecting `.bounds`.
Parameters
----------
other : :class:`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:`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: Optional[float] = None, y: Optional[float] = None, z: Optional[float] = 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) -> 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.
Returns
-------
bool
Whether this geometry intersects the plane.
"""
return self.bounds[0][axis] <= position <= self.bounds[1][axis]
@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:`Box` representation of the bounding box of a :class:`Geometry`.
Returns
-------
:class:`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:`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: Optional[float] = None,
y: Optional[float] = None,
z: Optional[float] = 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)
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: Optional[float] = None,
y: Optional[float] = None,
z: Optional[float] = 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
def _processed_coords(coords: Sequence[tuple[Any, ...]]) -> list[tuple[float, ...]]:
evaluated = Geometry._evaluate_inf(np.array(coords))
return [tuple(point) for point in evaluated.tolist()]
if shape.geom_type == "Polygon":
shell = _processed_coords(shape.exterior.coords)
holes = [_processed_coords(g.coords) for g in shape.interiors]
return shapely.Polygon(shell, holes)
if shape.geom_type in {"Point", "LineString", "LinearRing"}:
return shape.__class__(Geometry._evaluate_inf(np.array(shape.coords)))
if shape.geom_type in {
"MultiPoint",
"MultiLineString",
"MultiPolygon",
"GeometryCollection",
}:
return shape.__class__([Geometry.evaluate_inf_shape(g) for g in shape.geoms])
return shape
[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
[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:`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:`Geometry`
Scaled copy of this geometry.
"""
return Transformed(geometry=self, transform=Transformed.scaling(x, y, z))
[docs]
def rotated(self, angle: float, axis: Union[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:`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:`Geometry`
Reflected copy of this geometry.
"""
return Transformed(geometry=self, transform=Transformed.reflection(normal))
""" 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: Optional[int] = None,
gds_scale: pydantic.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: Optional[int] = None,
gds_scale: pydantic.PositiveFloat = 1.0,
dilation: float = 0.0,
sidewall_angle: float = 0,
reference_plane: PlanePosition = "middle",
) -> 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``.
Returns
-------
:class:`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'.")
gds_loader_fn = Geometry.load_gds_vertices_gdstk
geometries = []
with log as consolidated_logger:
for vertices in gds_loader_fn(gds_cell, gds_layer, gds_dtype, gds_scale):
# buffer(0) is necessary to merge self-intersections
shape = shapely.set_precision(shapely.Polygon(vertices).buffer(0), POLY_GRID_SIZE)
try:
geometries.append(
from_shapely(
shape, axis, slab_bounds, dilation, sidewall_angle, reference_plane
)
)
except pydantic.ValidationError as error:
consolidated_logger.warning(str(error))
except Tidy3dError as error:
consolidated_logger.warning(str(error))
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:`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: Optional[float] = None,
y: Optional[float] = None,
z: Optional[float] = None,
gds_layer: pydantic.NonNegativeInt = 0,
gds_dtype: pydantic.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: Optional[float] = None,
y: Optional[float] = None,
z: Optional[float] = None,
gds_layer: pydantic.NonNegativeInt = 0,
gds_dtype: pydantic.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: Optional[float] = None,
y: Optional[float] = None,
z: Optional[float] = None,
gds_layer: pydantic.NonNegativeInt = 0,
gds_dtype: pydantic.NonNegativeInt = 0,
gds_cell_name: str = "MAIN",
) -> 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:`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.
"""
try:
import gdstk
except ImportError as e:
raise Tidy3dImportError(
"Python module 'gdstk' not found. To export geometries to .gds "
"files, please install it."
) from e
library = gdstk.Library()
cell = library.new_cell(gds_cell_name)
self.to_gds(cell, x=x, y=y, z=z, gds_layer=gds_layer, gds_dtype=gds_dtype)
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 _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,)
[docs]
def __add__(self, other: Union[int, Geometry]) -> Union[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 # type: ignore[return-value]
return GeometryGroup(geometries=self._as_union() + other._as_union())
[docs]
def __radd__(self, other: Union[int, Geometry]) -> Union[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 # type: ignore[return-value]
return GeometryGroup(geometries=other._as_union() + self._as_union())
[docs]
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())
[docs]
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)
[docs]
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)
[docs]
def __sub__(self, other: Geometry) -> ClipOperation:
"""Difference of geometries"""
if not isinstance(other, Geometry):
return NotImplemented # type: ignore[return-value]
return ClipOperation(operation="difference", geometry_a=self, geometry_b=other)
[docs]
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)
[docs]
def __pos__(self) -> Self:
"""No op"""
return self
[docs]
def __neg__(self) -> ClipOperation:
"""Opposite of a geometry"""
return ClipOperation(
operation="difference", geometry_a=Box(size=(inf, inf, inf)), geometry_b=self
)
[docs]
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 = pydantic.Field(
(0.0, 0.0, 0.0),
title="Center",
description="Center of object in x, y, and z.",
units=MICROMETER,
)
@pydantic.validator("center", always=True)
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: Optional[int] = None,
) -> list[Shapely]:
"""Return a list of shapely geometries at the plane specified by normal and origin.
Checks special cases before relying on the complete computation.
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>`_.
"""
# 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, **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: Optional[int] = 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 = pydantic.Field(
2, title="Axis", description="Specifies dimension of the planar axis (0,1,2) -> (x,y,z)."
)
sidewall_angle: TracedFloat = pydantic.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.",
units=RADIAN,
)
reference_plane: PlanePosition = pydantic.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]
@pydantic.validator("sidewall_angle", always=True)
def validate_angle(cls, value: float) -> float:
lower_bound = -np.pi / 2
upper_bound = np.pi / 2
if (value <= lower_bound) or (value >= upper_bound):
# u03C0 is unicode for pi
raise ValidationError(f"Sidewall angle ({value}) must be between -Ο/2 and Ο/2 rad.")
return value
@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: Optional[float] = None,
y: Optional[float] = None,
z: Optional[float] = None,
cleanup: bool = True,
quad_segs: Optional[int] = None,
) -> list[Shapely]:
"""Returns shapely geometry at plane specified by one non None value of x,y,z.
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.
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)
if not self.intersects_axis_position(axis, position):
return []
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: Optional[int] = None) -> list:
"""Find shapely geometries intersecting planar geometry with axis normal to slab.
Parameters
----------
z : float
Position along the axis normal to slab
quad_segs : Optional[int] = None
Number of segments used to discretize circular shapes. If ``None``, uses
high-quality visualization settings.
Returns
-------
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: pydantic.NonNegativeFloat = pydantic.Field(
..., title="Radius", description="Radius of geometry.", units=MICROMETER
)
@pydantic.validator("radius", always=True)
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))
"""
size: TracedSize = pydantic.Field(
...,
title="Size",
description="Size in x, y, and z directions.",
units=MICROMETER,
)
[docs]
@classmethod
def from_bounds(cls, rmin: Coordinate, rmax: Coordinate, **kwargs: Any) -> Self:
"""Constructs a :class:`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)
[docs]
@classmethod
def surfaces(cls, size: Size, center: Coordinate, **kwargs: Any) -> list[Self]:
"""Returns a list of 6 :class:`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
name_base = kwargs.pop("name", "")
kwargs.pop("normal_dir", None)
names = []
normal_dirs = []
for coord in "xyz":
for direction in "-+":
surface_name = name_base + "_" + coord + direction
names.append(surface_name)
normal_dirs.append(direction)
# ignore surfaces that are infinitely far away
del_idx = []
for idx, _size in enumerate(size):
if _size == inf:
del_idx.append(idx)
del_idx = [[2 * i, 2 * i + 1] for i in del_idx]
del_idx = [item for sublist in del_idx for item in sublist]
def del_items(items: Iterable, indices: 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 = del_items(names, del_idx)
normal_dirs = del_items(normal_dirs, del_idx)
surfaces = []
for _cent, _size, _name, _normal_dir in zip(centers, sizes, names, normal_dirs):
if "normal_dir" in cls.__dict__["__fields__"]:
kwargs["normal_dir"] = _normal_dir
if "name" in cls.__dict__["__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:`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 "name" in cls.__dict__["__fields__"] and exclude_surfaces:
surfaces = [surf for surf in surfaces if surf.name[-2:] 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: Optional[int] = 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: Optional[float] = None,
y: Optional[float] = None,
z: Optional[float] = None,
cleanup: bool = True,
quad_segs: Optional[int] = None,
) -> list[Shapely]:
"""Returns shapely geometry at plane specified by one non None value of x,y,z.
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.
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)
if not self.intersects_axis_position(axis, position):
return []
z0, (x0, y0) = self.pop_axis(self.center, axis=axis)
Lz, (Lx, Ly) = self.pop_axis(self.size, axis=axis)
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:`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: Shapely, cleanup: bool = True, quad_segs: Optional[int] = None
) -> list[Shapely]:
"""Returns list of shapely geometries representing the intersections of the geometry with
this 2D box.
Parameters
----------
other : Shapely
Geometry to intersect with.
cleanup : bool = True
If True, removes extremely small features from each polygon's boundary.
quad_segs : Optional[int] = None
Number of segments used to discretize circular shapes. If ``None``, uses
high-quality visualization settings.
Returns
-------
List[shapely.geometry.base.BaseGeometry]
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."
)
# dont bother if the geometry doesn't intersect the self at all
if not other.intersects(self):
return []
# get list of Shapely shapes that intersect at the self
normal_ind = self.size.index(0.0)
dim = "xyz"[normal_ind]
pos = self.center[normal_ind]
xyz_kwargs = {dim: pos}
shapes_plane = other.intersections_plane(cleanup=cleanup, quad_segs=quad_segs, **xyz_kwargs)
# intersect all shapes with the input self
bs_min, bs_max = (self.pop_axis(bounds, axis=normal_ind)[1] for bounds in self.bounds)
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: Optional[tuple[pydantic.NonNegativeFloat, pydantic.NonNegativeFloat]] = None,
y: Optional[tuple[pydantic.NonNegativeFloat, pydantic.NonNegativeFloat]] = None,
z: Optional[tuple[pydantic.NonNegativeFloat, pydantic.NonNegativeFloat]] = None,
) -> Box:
"""Created a padded copy of a :class:`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:`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:`Box` representation of self (used for subclasses of Box).
Returns
-------
:class:`Box`
Instance of :class:`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:`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=new_center, size=new_size)
def _plot_arrow(
self,
direction: tuple[float, float, float],
x: Optional[float] = None,
y: Optional[float] = None,
z: Optional[float] = None,
color: Optional[str] = None,
alpha: Optional[float] = None,
bend_radius: Optional[float] = 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.
"""
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)) > 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]:
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, bounds_perp = 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)
sim_bounds_normal, sim_bounds_perp = self.pop_axis(
np.array(derivative_info.simulation_bounds).T, axis=axis_normal
)
# 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 = pydantic.Field(
...,
title="Operation Type",
description="Operation to be performed between geometries.",
)
geometry_a: annotate_type(GeometryType) = pydantic.Field(
...,
title="Geometry A",
description="First operand for the set operation. It can be any geometry type, including "
":class:`GeometryGroup`.",
)
geometry_b: annotate_type(GeometryType) = pydantic.Field(
...,
title="Geometry B",
description="Second operand for the set operation. It can also be any geometry type.",
)
@pydantic.validator("geometry_a", "geometry_b", always=True)
def _geometries_untraced(cls, val: GeometryType) -> GeometryType:
"""Make sure that ``ClipOperation`` geometries do not contain tracers."""
traced = val._strip_traced_fields()
if traced:
raise ValidationError(
f"{val.type} contains traced fields {list(traced.keys())}. Note that "
"'ClipOperation' does not currently support automatic differentiation."
)
return val
[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: Optional[int] = 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.
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
)
b = self.geometry_b.intersections_tilted_plane(
normal, origin, to_2D, cleanup=cleanup, quad_segs=quad_segs
)
geom_a = shapely.unary_union([Geometry.evaluate_inf_shape(g) for g in a])
geom_b = shapely.unary_union([Geometry.evaluate_inf_shape(g) for g in b])
return ClipOperation.to_polygon_list(
self._shapely_operation(geom_a, geom_b),
cleanup=cleanup,
)
[docs]
def intersections_plane(
self,
x: Optional[float] = None,
y: Optional[float] = None,
z: Optional[float] = None,
cleanup: bool = True,
quad_segs: Optional[int] = None,
) -> list[Shapely]:
"""Returns list of shapely geometries at plane specified by one non-None value of x,y,z.
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 Documentaton <https://shapely.readthedocs.io/en/stable/project.html>`_.
"""
a = self.geometry_a.intersections_plane(x, y, z, cleanup=cleanup, quad_segs=quad_segs)
b = self.geometry_b.intersections_plane(x, y, z, cleanup=cleanup, quad_segs=quad_segs)
geom_a = shapely.unary_union([Geometry.evaluate_inf_shape(g) for g in a])
geom_b = shapely.unary_union([Geometry.evaluate_inf_shape(g) for g in b])
return ClipOperation.to_polygon_list(
self._shapely_operation(geom_a, geom_b),
cleanup=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:`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)
[docs]
class GeometryGroup(Geometry):
"""A collection of Geometry objects that can be called as a single geometry object."""
geometries: tuple[annotate_type(GeometryType), ...] = pydantic.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.",
)
@pydantic.validator("geometries", always=True)
def _geometries_not_empty(
cls, val: tuple[annotate_type(GeometryType), ...]
) -> tuple[annotate_type(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: Optional[int] = 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.
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
)
]
[docs]
def intersections_plane(
self,
x: Optional[float] = None,
y: Optional[float] = None,
z: Optional[float] = None,
cleanup: bool = True,
quad_segs: Optional[int] = None,
) -> list[Shapely]:
"""Returns list of shapely geometries at plane specified by one non-None value of x,y,z.
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>`_.
"""
if not self.intersects_plane(x, y, z):
return []
return [
intersection
for geometry in self.geometries
for intersection in geometry.intersections_plane(
x=x, y=y, z=z, cleanup=cleanup, quad_segs=quad_segs
)
]
[docs]
def intersects_axis_position(self, axis: float, position: float) -> 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.
Returns
-------
bool
Whether this geometry intersects the plane.
"""
return any(geom.intersects_axis_position(axis, position) 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:`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."""
individual_volumes = (geometry.volume(bounds) for geometry in self.geometries)
return np.sum(individual_volumes)
def _surface_area(self, bounds: Bound) -> float:
"""Returns object's surface area within given bounds."""
individual_areas = (geometry.surface_area(bounds) for geometry in self.geometries)
return np.sum(individual_areas)
@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 = [
geometry._update_from_bounds(bounds=bounds, axis=axis) for geometry in self.geometries
]
return self.updated_copy(geometries=new_geometries)
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()
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
),
eps_approx=True,
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
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 _shapely_is_older_than("2.1"):
log.warning(
"Using old versions of the shapely library (prior to v2.1) may cause "
"plot errors. This can be solved by upgrading to Python 3.10 "
"(or later) and reinstalling Tidy3d.",
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([])
# In order to prevent numerical overflow or underflow errors, we first subtract
# the centroid and divide by (rescale) the size of the object so it is not too big.
normalized_obj = shapely.affinity.affine_transform(
# https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations
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( # This removes outward spikes
valid_obj,
distance=-tolerance_ratio,
cap_style="square", # (optional parameter to reduce computation time)
quad_segs=3, # (optional parameter to reduce computation time)
)
dilated_obj = shapely.buffer( # This removes inward spikes and tiny holes
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
