"""Concrete primitive geometrical objects."""
from __future__ import annotations
from math import isclose
from typing import TYPE_CHECKING, Any
import autograd.numpy as anp
import numpy as np
import shapely
from pydantic import Field, PrivateAttr, model_validator
from tidy3d.components.autograd import TracedSize1D, get_static
from tidy3d.components.base import cached_property
from tidy3d.components.geometry import base
from tidy3d.components.geometry.mesh import TriangleMesh
from tidy3d.components.geometry.polyslab import PolySlab
from tidy3d.config import config
from tidy3d.constants import LARGE_NUMBER, MICROMETER
from tidy3d.exceptions import SetupError, ValidationError
from tidy3d.log import log
from tidy3d.packaging import verify_packages_import
if TYPE_CHECKING:
from typing import Optional
from shapely.geometry.base import BaseGeometry
from tidy3d.compat import Self
from tidy3d.components.autograd import AutogradFieldMap
from tidy3d.components.autograd.derivative_utils import DerivativeInfo
from tidy3d.components.types import Axis, Bound, Coordinate, MatrixReal4x4, Shapely
# for sampling conical frustum in visualization
_N_SAMPLE_CURVE_SHAPELY = 40
# for shapely circular shapes discretization in visualization
_N_SHAPELY_QUAD_SEGS_VISUALIZATION = 200
# Default number of points to discretize polyslab in `Cylinder.to_polyslab()`
_N_PTS_CYLINDER_POLYSLAB = 51
_MAX_ICOSPHERE_SUBDIVISIONS = 7 # this would have 164K vertices and 328K faces
_DEFAULT_EDGE_FRACTION = 0.25
def _base_icosahedron() -> tuple[np.ndarray, np.ndarray]:
"""Return vertices and faces of a unit icosahedron."""
phi = (1.0 + np.sqrt(5.0)) / 2.0
vertices = np.array(
[
(-1, phi, 0),
(1, phi, 0),
(-1, -phi, 0),
(1, -phi, 0),
(0, -1, phi),
(0, 1, phi),
(0, -1, -phi),
(0, 1, -phi),
(phi, 0, -1),
(phi, 0, 1),
(-phi, 0, -1),
(-phi, 0, 1),
],
dtype=float,
)
vertices /= np.linalg.norm(vertices, axis=1)[:, None]
faces = np.array(
[
(0, 11, 5),
(0, 5, 1),
(0, 1, 7),
(0, 7, 10),
(0, 10, 11),
(1, 5, 9),
(5, 11, 4),
(11, 10, 2),
(10, 7, 6),
(7, 1, 8),
(3, 9, 4),
(3, 4, 2),
(3, 2, 6),
(3, 6, 8),
(3, 8, 9),
(4, 9, 5),
(2, 4, 11),
(6, 2, 10),
(8, 6, 7),
(9, 8, 1),
],
dtype=int,
)
return vertices, faces
_ICOSAHEDRON_VERTS, _ICOSAHEDRON_FACES = _base_icosahedron()
def discretization_wavelength(derivative_info: DerivativeInfo, geometry_label: str) -> float:
"""Choose reference wavelength for surface discretization."""
wvl0_min = derivative_info.wavelength_min
wvl_mat = wvl0_min / np.max([1.0, np.max(np.sqrt(abs(derivative_info.eps_in)))])
grid_cfg = config.adjoint
min_wvl_mat = grid_cfg.min_wvl_fraction * wvl0_min
if wvl_mat < min_wvl_mat:
log.warning(
f"The minimum wavelength inside the {geometry_label} material is {wvl_mat:.3e} μm, which would "
f"create a large number of discretization points for computing the gradient. "
f"To prevent performance degradation, the discretization wavelength has "
f"been clipped to {min_wvl_mat:.3e} μm.",
log_once=True,
)
return max(wvl_mat, min_wvl_mat)
[docs]
class Sphere(base.Centered, base.Circular):
"""Spherical geometry.
Example
-------
>>> b = Sphere(center=(1,2,3), radius=2)
"""
radius: TracedSize1D = Field(
title="Radius",
description="Radius of geometry.",
json_schema_extra={"units": MICROMETER},
)
_icosphere_cache: dict[int, tuple[np.ndarray, float]] = PrivateAttr(default_factory=dict)
[docs]
@verify_packages_import(["trimesh"])
def to_triangle_mesh(
self,
*,
max_edge_length: Optional[float] = None,
subdivisions: Optional[int] = None,
) -> TriangleMesh:
"""Approximate the sphere surface with a ``TriangleMesh``.
Parameters
----------
max_edge_length : float = None
Maximum edge length for triangulation in micrometers.
subdivisions : int = None
Number of subdivisions for icosphere generation.
Returns
-------
TriangleMesh
Triangle mesh approximation of the sphere surface.
"""
triangles, _ = self._triangulated_surface(
max_edge_length=max_edge_length, subdivisions=subdivisions
)
return TriangleMesh.from_triangles(triangles)
[docs]
def inside(
self, x: np.ndarray[float], y: np.ndarray[float], z: np.ndarray[float]
) -> np.ndarray[bool]:
"""For input arrays ``x``, ``y``, ``z`` of arbitrary but identical shape, return an array
with the same shape which is ``True`` for every point in zip(x, y, z) that is inside the
volume of the :class:`~tidy3d.Geometry`, and ``False`` otherwise.
Parameters
----------
x : np.ndarray[float]
Array of point positions in x direction.
y : np.ndarray[float]
Array of point positions in y direction.
z : np.ndarray[float]
Array of point positions in z direction.
Returns
-------
np.ndarray[bool]
``True`` for every point that is inside the geometry.
"""
self._ensure_equal_shape(x, y, z)
x0, y0, z0 = self.center
dist_x = np.abs(x - x0)
dist_y = np.abs(y - y0)
dist_z = np.abs(z - z0)
return (dist_x**2 + dist_y**2 + dist_z**2) <= (self.radius**2)
[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
``_N_SHAPELY_QUAD_SEGS_VISUALIZATION`` for high-quality visualization.
Returns
-------
list[shapely.geometry.base.BaseGeometry]
List of 2D shapes that intersect plane.
For more details refer to
`Shapely's Documentation <https://shapely.readthedocs.io/en/stable/project.html>`_.
"""
if quad_segs is None:
quad_segs = _N_SHAPELY_QUAD_SEGS_VISUALIZATION
normal = np.array(normal)
unit_normal = normal / (np.sum(normal**2) ** 0.5)
projection = np.dot(np.array(origin) - np.array(self.center), unit_normal)
if abs(projection) >= self.radius:
return []
radius = (self.radius**2 - projection**2) ** 0.5
center = np.array(self.center) + projection * unit_normal
v = np.zeros(3)
v[np.argmin(np.abs(unit_normal))] = 1
u = np.cross(unit_normal, v)
u /= np.sum(u**2) ** 0.5
v = np.cross(unit_normal, u)
angles = np.linspace(0, 2 * np.pi, quad_segs * 4 + 1)[:-1]
circ = center + np.outer(np.cos(angles), radius * u) + np.outer(np.sin(angles), radius * v)
vertices = np.dot(np.hstack((circ, np.ones((angles.size, 1)))), to_2D.T)
return [shapely.Polygon(vertices[:, :2])]
[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[BaseGeometry]:
"""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 x direction, only one of x,y,z can be specified to define plane.
z : float = None
Position of plane in x direction, only one of x,y,z can be specified to define plane.
cleanup : bool = True
If True, removes extremely small features from each polygon's boundary.
quad_segs : Optional[int] = None
Number of segments used to discretize circular shapes. If ``None``, uses
``_N_SHAPELY_QUAD_SEGS_VISUALIZATION`` for high-quality visualization.
Returns
-------
list[shapely.geometry.base.BaseGeometry]
List of 2D shapes that intersect plane.
For more details refer to
`Shapely's Documentation <https://shapely.readthedocs.io/en/stable/project.html>``.
"""
if quad_segs is None:
quad_segs = _N_SHAPELY_QUAD_SEGS_VISUALIZATION
axis, position = self.parse_xyz_kwargs(x=x, y=y, z=z)
if not self.intersects_axis_position(axis, position):
return []
z0, (x0, y0) = self.pop_axis(self.center, axis=axis)
intersect_dist = self._intersect_dist(position, z0)
if not intersect_dist:
return []
return [shapely.Point(x0, y0).buffer(0.5 * intersect_dist, quad_segs=quad_segs)]
@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)``.
"""
coord_min = tuple(c - self.radius for c in self.center)
coord_max = tuple(c + self.radius for c in self.center)
return (coord_min, coord_max)
def _volume(self, bounds: Bound) -> float:
"""Returns object's volume within given bounds."""
volume = 4.0 / 3.0 * np.pi * self.radius**3
# a very loose upper bound on how much of sphere is in bounds
for axis in range(3):
if self.center[axis] <= bounds[0][axis] or self.center[axis] >= bounds[1][axis]:
volume *= 0.5
return volume
def _surface_area(self, bounds: Bound) -> float:
"""Returns object's surface area within given bounds."""
area = 4.0 * np.pi * self.radius**2
# a very loose upper bound on how much of sphere is in bounds
for axis in range(3):
if self.center[axis] <= bounds[0][axis] or self.center[axis] >= bounds[1][axis]:
area *= 0.5
return area
[docs]
@classmethod
def unit_sphere_triangles(
cls,
*,
target_edge_length: Optional[float] = None,
subdivisions: Optional[int] = None,
) -> np.ndarray:
"""Return unit sphere triangles discretized via an icosphere."""
unit_tris = UNIT_SPHERE._unit_sphere_triangles(
target_edge_length=target_edge_length,
subdivisions=subdivisions,
copy_result=True,
)
return unit_tris
def _compute_derivatives(self, derivative_info: DerivativeInfo) -> AutogradFieldMap:
"""Compute adjoint derivatives using smooth sphere surface samples."""
valid_paths = {("radius",), *{("center", i) for i in range(3)}}
for path in derivative_info.paths:
if path not in valid_paths:
raise ValueError(
f"No derivative defined w.r.t. 'Sphere' field '{path}'. "
"Supported fields are 'radius' and 'center'."
)
if not derivative_info.paths:
return {}
grid_cfg = config.adjoint
radius = float(get_static(self.radius))
if radius == 0.0:
log.warning(
"Sphere gradients cannot be computed for zero radius; gradients are zero.",
log_once=True,
)
return dict.fromkeys(derivative_info.paths, 0.0)
wvl_mat = discretization_wavelength(derivative_info, "sphere")
target_edge = max(wvl_mat / grid_cfg.points_per_wavelength, np.finfo(float).eps)
triangles, _ = self._triangulated_surface(max_edge_length=target_edge)
triangles = triangles.astype(grid_cfg.gradient_dtype_float, copy=False)
sim_min, sim_max = (
np.asarray(arr, dtype=grid_cfg.gradient_dtype_float)
for arr in derivative_info.simulation_bounds
)
tol = config.adjoint.edge_clip_tolerance
sim_extents = sim_max - sim_min
collapsed_indices = np.flatnonzero(np.isclose(sim_extents, 0.0, atol=tol))
if collapsed_indices.size:
if collapsed_indices.size > 1:
return dict.fromkeys(derivative_info.paths, 0.0)
axis_idx = int(collapsed_indices[0])
plane_value = float(sim_min[axis_idx])
return self._compute_derivatives_collapsed_axis(
derivative_info=derivative_info,
axis_idx=axis_idx,
plane_value=plane_value,
)
trimesh_obj = TriangleMesh._triangles_to_trimesh(triangles)
vertices = np.asarray(trimesh_obj.vertices, dtype=grid_cfg.gradient_dtype_float)
center = np.asarray(self.center, dtype=grid_cfg.gradient_dtype_float)
verts_centered = vertices - center
norms = np.linalg.norm(verts_centered, axis=1, keepdims=True)
norms = np.where(norms == 0, 1, norms)
normals = verts_centered / norms
if vertices.size == 0:
return dict.fromkeys(derivative_info.paths, 0.0)
# get vertex weights
faces = np.asarray(trimesh_obj.faces, dtype=int)
face_areas = np.asarray(trimesh_obj.area_faces, dtype=grid_cfg.gradient_dtype_float)
weights = np.zeros(len(vertices), dtype=grid_cfg.gradient_dtype_float)
np.add.at(weights, faces[:, 0], face_areas / 3.0)
np.add.at(weights, faces[:, 1], face_areas / 3.0)
np.add.at(weights, faces[:, 2], face_areas / 3.0)
perp1, perp2 = self._tangent_basis_from_normals(normals)
valid_axes = np.abs(sim_max - sim_min) > tol
inside_mask = np.all(
vertices[:, valid_axes] >= (sim_min - tol)[valid_axes], axis=1
) & np.all(vertices[:, valid_axes] <= (sim_max + tol)[valid_axes], axis=1)
if not np.any(inside_mask):
return dict.fromkeys(derivative_info.paths, 0.0)
points = vertices[inside_mask]
normals_sel = normals[inside_mask]
perp1_sel = perp1[inside_mask]
perp2_sel = perp2[inside_mask]
weights_sel = weights[inside_mask]
interpolators = derivative_info.interpolators
if interpolators is None:
interpolators = derivative_info.create_interpolators(
dtype=grid_cfg.gradient_dtype_float
)
g = derivative_info.evaluate_gradient_at_points(
points,
normals_sel,
perp1_sel,
perp2_sel,
interpolators,
)
weighted = (weights_sel * g).real
grad_center = np.sum(weighted[:, None] * normals_sel, axis=0)
grad_radius = np.sum(weighted)
vjps: AutogradFieldMap = {}
for path in derivative_info.paths:
if path == ("radius",):
vjps[path] = float(grad_radius)
else:
_, idx = path
vjps[path] = float(grad_center[idx])
return vjps
def _compute_derivatives_collapsed_axis(
self,
derivative_info: DerivativeInfo,
axis_idx: int,
plane_value: float,
) -> AutogradFieldMap:
"""Delegate collapsed-axis gradients to a Cylinder cross section."""
tol = config.adjoint.edge_clip_tolerance
radius = float(self.radius)
center = np.asarray(self.center, dtype=float)
delta = plane_value - center[axis_idx]
radius_sq = radius**2 - delta**2
if radius_sq <= tol**2:
return dict.fromkeys(derivative_info.paths, 0.0)
radius_plane = float(np.sqrt(max(radius_sq, 0.0)))
if radius_plane <= tol:
return dict.fromkeys(derivative_info.paths, 0.0)
cyl_paths: set[tuple[str, int | None]] = set()
need_radius = False
for path in derivative_info.paths:
if path == ("radius",) or path == ("center", axis_idx):
cyl_paths.add(("radius",))
need_radius = True
elif path[0] == "center" and path[1] != axis_idx:
cyl_paths.add(("center", path[1]))
if not cyl_paths:
return dict.fromkeys(derivative_info.paths, 0.0)
cyl_center = center.copy()
cyl_center[axis_idx] = plane_value
cylinder = Cylinder(
center=tuple(cyl_center),
radius=radius_plane,
length=discretization_wavelength(derivative_info, "sphere") * 2.0,
axis=axis_idx,
)
bounds_min = list(cyl_center)
bounds_max = list(cyl_center)
for dim in range(3):
if dim == axis_idx:
continue
bounds_min[dim] = center[dim] - radius_plane
bounds_max[dim] = center[dim] + radius_plane
bounds = (tuple(bounds_min), tuple(bounds_max))
sim_min_arr, sim_max_arr = (
np.asarray(arr, dtype=float) for arr in derivative_info.simulation_bounds
)
intersect_min = tuple(max(bounds[0][i], sim_min_arr[i]) for i in range(3))
intersect_max = tuple(min(bounds[1][i], sim_max_arr[i]) for i in range(3))
if any(lo > hi for lo, hi in zip(intersect_min, intersect_max)):
return dict.fromkeys(derivative_info.paths, 0.0)
derivative_info_cyl = derivative_info.updated_copy(
paths=list(cyl_paths),
bounds=bounds,
bounds_intersect=(intersect_min, intersect_max),
)
vjps_cyl = cylinder._compute_derivatives(derivative_info_cyl)
result = dict.fromkeys(derivative_info.paths, 0.0)
vjp_radius = float(vjps_cyl.get(("radius",), 0.0)) if need_radius else 0.0
for path in derivative_info.paths:
if path == ("radius",):
result[path] = vjp_radius * (radius / radius_plane)
elif path == ("center", axis_idx):
result[path] = vjp_radius * (delta / radius_plane)
elif path[0] == "center" and path[1] != axis_idx:
result[path] = float(vjps_cyl.get(("center", path[1]), 0.0))
return result
def _edge_length_on_unit_sphere(
self, max_edge_length: Optional[float] = _DEFAULT_EDGE_FRACTION
) -> Optional[float]:
"""Convert ``max_edge_length`` in μm to unit-sphere coordinates."""
max_edge_length = _DEFAULT_EDGE_FRACTION if max_edge_length is None else max_edge_length
radius = float(self.radius)
if radius <= 0.0:
return None
return max_edge_length / radius
def _triangulated_surface(
self,
*,
max_edge_length: Optional[float] = None,
subdivisions: Optional[int] = None,
) -> tuple[np.ndarray, np.ndarray]:
"""Return physical and unit triangles for the surface discretization. Pass either max_edge_length or subdivisions."""
max_edge_length_unit = None
if subdivisions is None:
max_edge_length_unit = self._edge_length_on_unit_sphere(max_edge_length)
unit_tris = self._unit_sphere_triangles(
target_edge_length=max_edge_length_unit,
subdivisions=subdivisions,
copy_result=False,
)
radius = float(get_static(self.radius))
center = np.asarray(self.center, dtype=float)
dtype = config.adjoint.gradient_dtype_float
physical = radius * unit_tris + center
return physical.astype(dtype, copy=False), unit_tris.astype(dtype, copy=False)
def _unit_sphere_triangles(
self,
*,
target_edge_length: Optional[float] = None,
subdivisions: Optional[int] = None,
copy_result: bool = True,
) -> np.ndarray:
"""Return cached unit-sphere triangles with optional copying. Pass either target_edge_length or subdivisions."""
if target_edge_length is not None and subdivisions is not None:
raise ValueError("Specify either target_edge_length OR subdivisions, not both.")
if subdivisions is None:
subdivisions = self._subdivisions_for_edge(target_edge_length)
triangles, _ = self._icosphere_data(subdivisions)
return np.array(triangles, copy=copy_result)
def _subdivisions_for_edge(self, target_edge_length: Optional[float]) -> int:
if target_edge_length is None or target_edge_length <= 0.0:
return 0
for subdiv in range(_MAX_ICOSPHERE_SUBDIVISIONS + 1):
_, max_edge = self._icosphere_data(subdiv)
if max_edge <= target_edge_length:
return subdiv
log.warning(
f"Requested sphere mesh edge length {target_edge_length:.3e} μm requires more than "
f"{_MAX_ICOSPHERE_SUBDIVISIONS} subdivisions. "
"Clipping to the finest available mesh.",
log_once=True,
)
return _MAX_ICOSPHERE_SUBDIVISIONS
@staticmethod
def _tangent_basis_from_normals(normals: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
"""Construct orthonormal tangential bases for each normal vector (vectorized)."""
dtype = normals.dtype
tol = np.finfo(dtype).eps
# Normalize normals (in case they are not perfectly unit length).
n_norm = np.linalg.norm(normals, axis=1)
n = normals / np.maximum(n_norm, tol)[:, None]
# Pick a reference axis least aligned with each normal: argmin(|nx|,|ny|,|nz|).
ref_idx = np.argmin(np.abs(n), axis=1)
ref = np.zeros_like(n)
ref[np.arange(n.shape[0]), ref_idx] = 1.0
basis1 = np.cross(n, ref)
b1_norm = np.linalg.norm(basis1, axis=1)
basis1 = basis1 / np.maximum(b1_norm, tol)[:, None]
basis2 = np.cross(n, basis1)
b2_norm = np.linalg.norm(basis2, axis=1)
basis2 = basis2 / np.maximum(b2_norm, tol)[:, None]
return basis1, basis2
def _icosphere_data(self, subdivisions: int) -> tuple[np.ndarray, float]:
cache = self._icosphere_cache
if subdivisions in cache:
return cache[subdivisions]
vertices = np.asarray(_ICOSAHEDRON_VERTS, dtype=float)
faces = np.asarray(_ICOSAHEDRON_FACES, dtype=int)
if subdivisions > 0:
vertices = vertices.copy()
faces = faces.copy()
for _ in range(subdivisions):
vertices, faces = TriangleMesh.subdivide_faces(vertices, faces)
norms = np.linalg.norm(vertices, axis=1, keepdims=True)
norms = np.where(norms == 0.0, 1.0, norms)
vertices = vertices / norms
triangles = vertices[faces]
max_edge = self._max_edge_length(triangles)
cache[subdivisions] = (triangles, max_edge)
return triangles, max_edge
@staticmethod
def _max_edge_length(triangles: np.ndarray) -> float:
v = triangles
edges = np.stack(
[
v[:, 1] - v[:, 0],
v[:, 2] - v[:, 1],
v[:, 0] - v[:, 2],
],
axis=1,
)
return float(np.linalg.norm(edges, axis=2).max())
UNIT_SPHERE = Sphere(center=(0.0, 0.0, 0.0), radius=1.0)
[docs]
class Cylinder(base.Centered, base.Circular, base.Planar):
"""Cylindrical geometry with optional sidewall angle along axis
direction. When ``sidewall_angle`` is nonzero, the shape is a
conical frustum or a cone.
Example
-------
>>> c = Cylinder(center=(1,2,3), radius=2, length=5, axis=2)
See Also
--------
**Notebooks**
* `THz integrated demultiplexer/filter based on a ring resonator <../../../notebooks/THzDemultiplexerFilter.html>`_
* `Photonic crystal waveguide polarization filter <../../../notebooks/PhotonicCrystalWaveguidePolarizationFilter.html>`_
"""
# Provide more explanations on where radius is defined
radius: TracedSize1D = Field(
title="Radius",
description="Radius of geometry at the ``reference_plane``.",
json_schema_extra={"units": MICROMETER},
)
length: TracedSize1D = Field(
title="Length",
description="Defines thickness of cylinder along axis dimension.",
json_schema_extra={"units": MICROMETER},
)
@model_validator(mode="after")
def _only_middle_for_infinite_length_slanted_cylinder(self: Self) -> Self:
"""For a slanted cylinder of infinite length, ``reference_plane`` can only
be ``middle``; otherwise, the radius at ``center`` is either td.inf or 0.
"""
if isclose(self.sidewall_angle, 0) or not np.isinf(self.length):
return self
if self.reference_plane != "middle":
raise SetupError(
"For a slanted cylinder here is of infinite length, "
"defining the reference_plane other than 'middle' "
"leads to undefined cylinder behaviors near 'center'."
)
return self
[docs]
def to_polyslab(
self, num_pts_circumference: int = _N_PTS_CYLINDER_POLYSLAB, **kwargs: Any
) -> PolySlab:
"""Convert instance of ``Cylinder`` into a discretized version using ``PolySlab``.
Parameters
----------
num_pts_circumference : int = 51
Number of points in the circumference of the discretized polyslab.
**kwargs:
Extra keyword arguments passed to ``PolySlab()``, such as ``dilation``.
Returns
-------
PolySlab
Extruded polygon representing a discretized version of the cylinder.
"""
center_axis = self.center_axis
length_axis = self.length_axis
slab_bounds = (center_axis - length_axis / 2.0, center_axis + length_axis / 2.0)
if num_pts_circumference < 3:
raise ValueError("'PolySlab' from 'Cylinder' must have 3 or more radius points.")
_, (x0, y0) = self.pop_axis(self.center, axis=self.axis)
xs_, ys_ = self._points_unit_circle(num_pts_circumference=num_pts_circumference)
xs = x0 + self.radius * xs_
ys = y0 + self.radius * ys_
vertices = anp.stack((xs, ys), axis=-1)
return PolySlab(
vertices=vertices,
axis=self.axis,
slab_bounds=slab_bounds,
sidewall_angle=self.sidewall_angle,
reference_plane=self.reference_plane,
**kwargs,
)
def _points_unit_circle(
self, num_pts_circumference: int = _N_PTS_CYLINDER_POLYSLAB
) -> np.ndarray:
"""Set of x and y points for the unit circle when discretizing cylinder as a polyslab."""
angles = np.linspace(0, 2 * np.pi, num_pts_circumference, endpoint=False)
xs = np.cos(angles)
ys = np.sin(angles)
return np.stack((xs, ys), axis=0)
def _compute_derivatives(self, derivative_info: DerivativeInfo) -> AutogradFieldMap:
"""Compute the adjoint derivatives for this object."""
# compute circumference discretization
wvl_mat = discretization_wavelength(derivative_info, "cylinder")
circumference = 2 * np.pi * self.radius
wvls_in_circumference = circumference / wvl_mat
grid_cfg = config.adjoint
num_pts_circumference = int(np.ceil(grid_cfg.points_per_wavelength * wvls_in_circumference))
num_pts_circumference = max(3, num_pts_circumference)
# construct equivalent polyslab and compute the derivatives
polyslab = self.to_polyslab(num_pts_circumference=num_pts_circumference)
# build PolySlab derivative paths based on requested Cylinder paths
ps_paths = set()
for path in derivative_info.paths:
if path == ("length",):
ps_paths.update({("slab_bounds", 0), ("slab_bounds", 1)})
elif path == ("radius",):
ps_paths.add(("vertices",))
elif "center" in path:
_, center_index = path
_, (index_x, index_y) = self.pop_axis((0, 1, 2), axis=self.axis)
if center_index in (index_x, index_y):
ps_paths.add(("vertices",))
else:
ps_paths.update({("slab_bounds", 0), ("slab_bounds", 1)})
elif path == ("sidewall_angle",):
ps_paths.add(("sidewall_angle",))
# pass interpolators to PolySlab if available to avoid redundant conversions
update_kwargs = {
"paths": list(ps_paths),
"deep": False,
}
if derivative_info.interpolators is not None:
update_kwargs["interpolators"] = derivative_info.interpolators
derivative_info_polyslab = derivative_info.updated_copy(**update_kwargs)
vjps_polyslab = polyslab._compute_derivatives(derivative_info_polyslab)
vjps = {}
for path in derivative_info.paths:
if path == ("length",):
vjp_top = vjps_polyslab.get(("slab_bounds", 0), 0.0)
vjp_bot = vjps_polyslab.get(("slab_bounds", 1), 0.0)
vjps[path] = vjp_top - vjp_bot
elif path == ("radius",):
# transform polyslab vertices derivatives into radius derivative
xs_, ys_ = self._points_unit_circle(num_pts_circumference=num_pts_circumference)
if ("vertices",) not in vjps_polyslab:
vjps[path] = 0.0
else:
vjps_vertices_xs, vjps_vertices_ys = vjps_polyslab[("vertices",)].T
vjp_xs = np.sum(xs_ * vjps_vertices_xs)
vjp_ys = np.sum(ys_ * vjps_vertices_ys)
vjps[path] = vjp_xs + vjp_ys
elif "center" in path:
_, center_index = path
_, (index_x, index_y) = self.pop_axis((0, 1, 2), axis=self.axis)
if center_index == index_x:
if ("vertices",) not in vjps_polyslab:
vjps[path] = 0.0
else:
vjps_vertices_xs = vjps_polyslab[("vertices",)][:, 0]
vjps[path] = np.sum(vjps_vertices_xs)
elif center_index == index_y:
if ("vertices",) not in vjps_polyslab:
vjps[path] = 0.0
else:
vjps_vertices_ys = vjps_polyslab[("vertices",)][:, 1]
vjps[path] = np.sum(vjps_vertices_ys)
else:
vjp_top = vjps_polyslab.get(("slab_bounds", 0), 0.0)
vjp_bot = vjps_polyslab.get(("slab_bounds", 1), 0.0)
vjps[path] = vjp_top + vjp_bot
elif path == ("sidewall_angle",):
# direct mapping: cylinder angle equals polyslab angle
vjps[path] = vjps_polyslab.get(("sidewall_angle",), 0.0)
else:
raise NotImplementedError(
f"Differentiation with respect to 'Cylinder' '{path}' field not supported. "
"If you would like this feature added, please feel free to raise "
"an issue on the tidy3d front end repository."
)
return vjps
@property
def center_axis(self) -> Any:
"""Gets the position of the center of the geometry in the out of plane dimension."""
z0, _ = self.pop_axis(self.center, axis=self.axis)
return z0
@property
def length_axis(self) -> float:
"""Gets the length of the geometry along the out of plane dimension."""
return self.length
@cached_property
def _normal_2dmaterial(self) -> Axis:
"""Get the normal to the given geometry, checking that it is a 2D geometry."""
if self.length != 0:
raise ValidationError("'Medium2D' requires the 'Cylinder' length to be zero.")
return self.axis
def _update_from_bounds(self, bounds: tuple[float, float], axis: Axis) -> Cylinder:
"""Returns an updated geometry which has been transformed to fit within ``bounds``
along the ``axis`` direction."""
if axis != self.axis:
raise ValueError(
f"'_update_from_bounds' may only be applied along axis '{self.axis}', "
f"but was given axis '{axis}'."
)
new_center = list(self.center)
new_center[axis] = (bounds[0] + bounds[1]) / 2
new_length = bounds[1] - bounds[0]
return self.updated_copy(center=tuple(new_center), length=new_length)
@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. If ``None``, uses
``_N_SHAPELY_QUAD_SEGS_VISUALIZATION`` for high-quality visualization.
Returns
-------
list[shapely.geometry.base.BaseGeometry]
List of 2D shapes that intersect plane.
For more details refer to
`Shapely's Documentation <https://shapely.readthedocs.io/en/stable/project.html>`_.
"""
import trimesh
if quad_segs is None:
quad_segs = _N_SHAPELY_QUAD_SEGS_VISUALIZATION
z0, (x0, y0) = self.pop_axis(self.center, self.axis)
half_length = self.finite_length_axis / 2
z_top = z0 + half_length
z_bot = z0 - half_length
if np.isclose(self.sidewall_angle, 0):
r_top = self.radius
r_bot = self.radius
else:
r_top = self.radius_top
r_bot = self.radius_bottom
if r_top < 0 or np.isclose(r_top, 0):
r_top = 0
z_top = z0 + self._radius_z(z0) / self._tanq
elif r_bot < 0 or np.isclose(r_bot, 0):
r_bot = 0
z_bot = z0 + self._radius_z(z0) / self._tanq
angles = np.linspace(0, 2 * np.pi, quad_segs * 4 + 1)
if r_bot > 0:
x_bot = x0 + r_bot * np.cos(angles)
y_bot = y0 + r_bot * np.sin(angles)
x_bot[-1] = x0
y_bot[-1] = y0
else:
x_bot = np.array([x0])
y_bot = np.array([y0])
if r_top > 0:
x_top = x0 + r_top * np.cos(angles)
y_top = y0 + r_top * np.sin(angles)
x_top[-1] = x0
y_top[-1] = y0
else:
x_top = np.array([x0])
y_top = np.array([y0])
x = np.hstack((x_bot, x_top))
y = np.hstack((y_bot, y_top))
z = np.hstack((np.full_like(x_bot, z_bot), np.full_like(x_top, z_top)))
vertices = np.vstack(self.unpop_axis(z, (x, y), self.axis)).T
if x_bot.shape[0] == 1:
m = 1
n = x_top.shape[0] - 1
faces_top = [(m + n, m + i, m + (i + 1) % n) for i in range(n)]
faces_side = [(m + (i + 1) % n, m + i, 0) for i in range(n)]
faces = faces_top + faces_side
elif x_top.shape[0] == 1:
m = x_bot.shape[0]
n = m - 1
faces_bot = [(n, (i + 1) % n, i) for i in range(n)]
faces_side = [(i, (i + 1) % n, m) for i in range(n)]
faces = faces_bot + faces_side
else:
m = x_bot.shape[0]
n = m - 1
faces_bot = [(n, (i + 1) % n, i) for i in range(n)]
faces_top = [(m + n, m + i, m + (i + 1) % n) for i in range(n)]
faces_side_bot = [(i, (i + 1) % n, m + (i + 1) % n) for i in range(n)]
faces_side_top = [(m + (i + 1) % n, m + i, i) for i in range(n)]
faces = faces_bot + faces_top + faces_side_bot + faces_side_top
mesh = trimesh.Trimesh(vertices, faces)
section = mesh.section(plane_origin=origin, plane_normal=normal)
if section is None:
return []
path, _ = section.to_2D(to_2D=to_2D)
return path.polygons_full
def _intersections_normal(
self, z: float, quad_segs: Optional[int] = None
) -> list[BaseGeometry]:
"""Find shapely geometries intersecting cylindrical geometry with axis normal to slab.
Parameters
----------
z : float
Position along the axis normal to slab
quad_segs : Optional[int] = None
Number of segments used to discretize circular shapes. If ``None``, uses
``_N_SHAPELY_QUAD_SEGS_VISUALIZATION`` for high-quality visualization.
Returns
-------
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 quad_segs is None:
quad_segs = _N_SHAPELY_QUAD_SEGS_VISUALIZATION
static_self = self.to_static()
# radius at z
radius_offset = static_self._radius_z(z)
if radius_offset <= 0:
return []
_, (x0, y0) = self.pop_axis(static_self.center, axis=self.axis)
return [shapely.Point(x0, y0).buffer(radius_offset, quad_segs=quad_segs)]
def _intersections_side(self, position: float, axis: int) -> list[BaseGeometry]:
"""Find shapely geometries intersecting cylindrical geometry with axis orthogonal to length.
When ``sidewall_angle`` is nonzero, so that it's in fact a conical frustum or cone, the
cross section can contain hyperbolic curves. This is currently approximated by a polygon
of many vertices.
Parameters
----------
position : float
Position along axis direction.
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>`_.
"""
# position in the local coordinate of the cylinder
position_local = position - self.center[axis]
# no intersection
if abs(position_local) >= self.radius_max:
return []
# half of intersection length at the top and bottom
intersect_half_length_max = np.sqrt(self.radius_max**2 - position_local**2)
intersect_half_length_min = -LARGE_NUMBER
if abs(position_local) < self.radius_min:
intersect_half_length_min = np.sqrt(self.radius_min**2 - position_local**2)
# the vertices on the max side of top/bottom
# The two vertices are present in all scenarios.
vertices_max = [
self._local_to_global_side_cross_section([-intersect_half_length_max, 0], axis),
self._local_to_global_side_cross_section([intersect_half_length_max, 0], axis),
]
# Extending to a cone, the maximal height of the cone
h_cone = (
LARGE_NUMBER if isclose(self.sidewall_angle, 0) else self.radius_max / abs(self._tanq)
)
# The maximal height of the cross section
height_max = min(
(1 - abs(position_local) / self.radius_max) * h_cone, self.finite_length_axis
)
# more vertices to add for conical frustum shape
vertices_frustum_right = []
vertices_frustum_left = []
if not (isclose(position, self.center[axis]) or isclose(self.sidewall_angle, 0)):
# The y-coordinate for the additional vertices
y_list = height_max * np.linspace(0, 1, _N_SAMPLE_CURVE_SHAPELY)
# `abs()` to make sure np.sqrt(0-fp_eps) goes through
x_list = np.sqrt(
np.abs(self.radius_max**2 * (1 - y_list / h_cone) ** 2 - position_local**2)
)
for i in range(_N_SAMPLE_CURVE_SHAPELY):
vertices_frustum_right.append(
self._local_to_global_side_cross_section([x_list[i], y_list[i]], axis)
)
vertices_frustum_left.append(
self._local_to_global_side_cross_section(
[
-x_list[_N_SAMPLE_CURVE_SHAPELY - i - 1],
y_list[_N_SAMPLE_CURVE_SHAPELY - i - 1],
],
axis,
)
)
# the vertices on the min side of top/bottom
vertices_min = []
## termination at the top/bottom
if intersect_half_length_min > 0:
vertices_min.append(
self._local_to_global_side_cross_section(
[intersect_half_length_min, self.finite_length_axis], axis
)
)
vertices_min.append(
self._local_to_global_side_cross_section(
[-intersect_half_length_min, self.finite_length_axis], axis
)
)
## early termination
else:
vertices_min.append(self._local_to_global_side_cross_section([0, height_max], axis))
return [
shapely.Polygon(
vertices_max + vertices_frustum_right + vertices_min + vertices_frustum_left
)
]
[docs]
def inside(
self, x: np.ndarray[float], y: np.ndarray[float], z: np.ndarray[float]
) -> np.ndarray[bool]:
"""For input arrays ``x``, ``y``, ``z`` of arbitrary but identical shape, return an array
with the same shape which is ``True`` for every point in zip(x, y, z) that is inside the
volume of the :class:`~tidy3d.Geometry`, and ``False`` otherwise.
Parameters
----------
x : np.ndarray[float]
Array of point positions in x direction.
y : np.ndarray[float]
Array of point positions in y direction.
z : np.ndarray[float]
Array of point positions in z direction.
Returns
-------
np.ndarray[bool]
``True`` for every point that is inside the geometry.
"""
# radius at z
self._ensure_equal_shape(x, y, z)
z0, (x0, y0) = self.pop_axis(self.center, axis=self.axis)
z, (x, y) = self.pop_axis((x, y, z), axis=self.axis)
radius_offset = self._radius_z(z)
positive_radius = radius_offset > 0
dist_x = np.abs(x - x0)
dist_y = np.abs(y - y0)
dist_z = np.abs(z - z0)
inside_radius = (dist_x**2 + dist_y**2) <= (radius_offset**2)
inside_height = dist_z <= (self.finite_length_axis / 2)
return positive_radius * inside_radius * inside_height
@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)``.
"""
coord_min = [c - self.radius_max for c in self.center]
coord_max = [c + self.radius_max for c in self.center]
coord_min[self.axis] = self.center[self.axis] - self.length_axis / 2.0
coord_max[self.axis] = self.center[self.axis] + self.length_axis / 2.0
return (tuple(coord_min), tuple(coord_max))
def _volume(self, bounds: Bound) -> float:
"""Returns object's volume within given bounds."""
coord_min = max(self.bounds[0][self.axis], bounds[0][self.axis])
coord_max = min(self.bounds[1][self.axis], bounds[1][self.axis])
length = coord_max - coord_min
volume = np.pi * self.radius_max**2 * length
# a very loose upper bound on how much of the cylinder is in bounds
for axis in range(3):
if axis != self.axis:
if self.center[axis] <= bounds[0][axis] or self.center[axis] >= bounds[1][axis]:
volume *= 0.5
return volume
def _surface_area(self, bounds: Bound) -> float:
"""Returns object's surface area within given bounds."""
area = 0
coord_min = self.bounds[0][self.axis]
coord_max = self.bounds[1][self.axis]
if coord_min < bounds[0][self.axis]:
coord_min = bounds[0][self.axis]
else:
area += np.pi * self.radius_max**2
if coord_max > bounds[1][self.axis]:
coord_max = bounds[1][self.axis]
else:
area += np.pi * self.radius_max**2
length = coord_max - coord_min
area += 2.0 * np.pi * self.radius_max * length
# a very loose upper bound on how much of the cylinder is in bounds
for axis in range(3):
if axis != self.axis:
if self.center[axis] <= bounds[0][axis] or self.center[axis] >= bounds[1][axis]:
area *= 0.5
return area
@cached_property
def radius_bottom(self) -> float:
"""radius of bottom"""
return self._radius_z(self.center_axis - self.finite_length_axis / 2)
@cached_property
def radius_top(self) -> float:
"""radius of bottom"""
return self._radius_z(self.center_axis + self.finite_length_axis / 2)
@cached_property
def radius_max(self) -> float:
"""max(radius of top, radius of bottom)"""
return max(self.radius_bottom, self.radius_top)
@cached_property
def radius_min(self) -> float:
"""min(radius of top, radius of bottom). It can be negative for a large
sidewall angle.
"""
return min(self.radius_bottom, self.radius_top)
def _radius_z(self, z: float) -> float:
"""Compute the radius of the cross section at the position z.
Parameters
----------
z : float
Position along the axis normal to slab
"""
if isclose(self.sidewall_angle, 0):
return self.radius
radius_middle = self.radius
if self.reference_plane == "top":
radius_middle += self.finite_length_axis / 2 * self._tanq
elif self.reference_plane == "bottom":
radius_middle -= self.finite_length_axis / 2 * self._tanq
return radius_middle - (z - self.center_axis) * self._tanq
def _local_to_global_side_cross_section(self, coords: list[float], axis: int) -> list[float]:
"""Map a point (x,y) from local to global coordinate system in the
side cross section.
The definition of the local: y=0 lies at the base if ``sidewall_angle>=0``,
and at the top if ``sidewall_angle<0``; x=0 aligns with the corresponding
``self.center``. In both cases, y-axis is pointing towards the narrowing
direction of cylinder.
Parameters
----------
axis : int
Integer index into 'xyz' (0, 1, 2).
coords : list[float, float]
The value in the planar coordinate.
Returns
-------
Tuple[float, float]
The point in the global coordinate for plotting `_intersection_side`.
"""
# For negative sidewall angle, quantities along axis direction usually needs a flipped sign
axis_sign = 1
if self.sidewall_angle < 0:
axis_sign = -1
lx_offset, ly_offset = self._order_by_axis(
plane_val=coords[0],
axis_val=axis_sign * (-self.finite_length_axis / 2 + coords[1]),
axis=axis,
)
_, (x_center, y_center) = self.pop_axis(self.center, axis=axis)
return [x_center + lx_offset, y_center + ly_offset]
