"""Geometry extruded from polygonal shapes."""
from __future__ import annotations
import math
from copy import copy
from functools import lru_cache
from typing import Optional, Union
import autograd.numpy as np
import pydantic.v1 as pydantic
import shapely
from autograd.tracer import getval, isbox
from numpy.polynomial.legendre import leggauss as _leggauss
from tidy3d.components.autograd import AutogradFieldMap, TracedVertices, get_static
from tidy3d.components.autograd.constants import (
EDGE_CLIP_TOLERANCE,
GAUSS_QUADRATURE_ORDER,
GRADIENT_DTYPE_FLOAT,
QUAD_SAMPLE_FRACTION,
)
from tidy3d.components.autograd.derivative_utils import DerivativeInfo
from tidy3d.components.autograd.types import TracedFloat
from tidy3d.components.base import cached_property, skip_if_fields_missing
from tidy3d.components.transformation import ReflectionFromPlane, RotationAroundAxis
from tidy3d.components.types import (
ArrayFloat1D,
ArrayFloat2D,
ArrayLike,
Axis,
Bound,
Coordinate,
MatrixReal4x4,
PlanePosition,
Shapely,
)
from tidy3d.constants import LARGE_NUMBER, MICROMETER, fp_eps
from tidy3d.exceptions import SetupError, Tidy3dImportError, ValidationError
from tidy3d.log import log
from tidy3d.packaging import verify_packages_import
from . import base, triangulation
# sampling polygon along dilation for validating polygon to be
# non self-intersecting during the entire dilation process
_N_SAMPLE_POLYGON_INTERSECT = 5
_IS_CLOSE_RTOL = np.finfo(float).eps
# Warn for too many divided polyslabs
_COMPLEX_POLYSLAB_DIVISIONS_WARN = 100
# Warn before triangulating large polyslabs due to inefficiency
_MAX_POLYSLAB_VERTICES_FOR_TRIANGULATION = 500
_MIN_POLYGON_AREA = fp_eps
@lru_cache(maxsize=128)
def leggauss(n):
"""Cached version of leggauss with dtype conversions."""
g, w = _leggauss(n)
return g.astype(GRADIENT_DTYPE_FLOAT, copy=False), w.astype(GRADIENT_DTYPE_FLOAT, copy=False)
[docs]
class PolySlab(base.Planar):
"""Polygon extruded with optional sidewall angle along axis direction.
Example
-------
>>> vertices = np.array([(0,0), (1,0), (1,1)])
>>> p = PolySlab(vertices=vertices, axis=2, slab_bounds=(-1, 1))
"""
slab_bounds: tuple[TracedFloat, TracedFloat] = pydantic.Field(
...,
title="Slab Bounds",
description="Minimum and maximum positions of the slab along axis dimension.",
units=MICROMETER,
)
dilation: float = pydantic.Field(
0.0,
title="Dilation",
description="Dilation of the supplied polygon by shifting each edge along its "
"normal outwards direction by a distance; a negative value corresponds to erosion.",
units=MICROMETER,
)
vertices: TracedVertices = pydantic.Field(
...,
title="Vertices",
description="List of (d1, d2) defining the 2 dimensional positions of the polygon "
"face vertices at the ``reference_plane``. "
"The index of dimension should be in the ascending order: e.g. if "
"the slab normal axis is ``axis=y``, the coordinate of the vertices will be in (x, z)",
units=MICROMETER,
)
[docs]
@staticmethod
def make_shapely_polygon(vertices: ArrayLike) -> shapely.Polygon:
"""Make a shapely polygon from some vertices, first ensures they are untraced."""
vertices = get_static(vertices)
return shapely.Polygon(vertices)
[docs]
@pydantic.validator("slab_bounds", always=True)
def slab_bounds_order(cls, val):
"""Maximum position of the slab should be no smaller than its minimal position."""
if val[1] < val[0]:
raise SetupError(
"Polyslab.slab_bounds must be specified in the order of "
"minimum and maximum positions of the slab along the axis. "
f"But now the maximum {val[1]} is smaller than the minimum {val[0]}."
)
return val
[docs]
@pydantic.validator("vertices", always=True)
def correct_shape(cls, val):
"""Makes sure vertices size is correct.
Make sure no intersecting edges.
"""
# overall shape of vertices
if val.shape[1] != 2:
raise SetupError(
"PolySlab.vertices must be a 2 dimensional array shaped (N, 2). "
f"Given array with shape of {val.shape}."
)
# make sure no polygon splitting, islands, 0 area
poly_heal = shapely.make_valid(cls.make_shapely_polygon(val))
if poly_heal.area < _MIN_POLYGON_AREA:
raise SetupError("The polygon almost collapses to a 1D curve.")
if not poly_heal.geom_type == "Polygon" or len(poly_heal.interiors) > 0:
raise SetupError(
"Polygon is self-intersecting, resulting in "
"polygon splitting or generation of holes/islands. "
"A general treatment to self-intersecting polygon will be available "
"in future releases."
)
return val
[docs]
@pydantic.validator("vertices", always=True)
@skip_if_fields_missing(["dilation"])
def no_complex_self_intersecting_polygon_at_reference_plane(cls, val, values):
"""At the reference plane, check if the polygon is self-intersecting.
There are two types of self-intersection that can occur during dilation:
1) the one that creates holes/islands, or splits polygons, or removes everything;
2) the one that does not.
For 1), we issue an error since it is yet to be supported;
For 2), we heal the polygon, and warn that the polygon has been cleaned up.
"""
# no need to validate anything here
if math.isclose(values["dilation"], 0):
return val
val_np = PolySlab._proper_vertices(val)
dist = values["dilation"]
# 0) fully eroded
if dist < 0 and dist < -PolySlab._maximal_erosion(val_np):
raise SetupError("Erosion value is too large. The polygon is fully eroded.")
# no edge events
if not PolySlab._edge_events_detection(val_np, dist, ignore_at_dist=False):
return val
poly_offset = PolySlab._shift_vertices(val_np, dist)[0]
if PolySlab._area(poly_offset) < fp_eps**2:
raise SetupError("Erosion value is too large. The polygon is fully eroded.")
# edge events
poly_offset = shapely.make_valid(cls.make_shapely_polygon(poly_offset))
# 1) polygon split or create holes/islands
if not poly_offset.geom_type == "Polygon" or len(poly_offset.interiors) > 0:
raise SetupError(
"Dilation/Erosion value is too large, resulting in "
"polygon splitting or generation of holes/islands. "
"A general treatment to self-intersecting polygon will be available "
"in future releases."
)
# case 2
log.warning(
"The dilation/erosion value is too large. resulting in a "
"self-intersecting polygon. "
"The vertices have been modified to make a valid polygon."
)
return val
[docs]
@pydantic.validator("vertices", always=True)
@skip_if_fields_missing(["sidewall_angle", "dilation", "slab_bounds", "reference_plane"])
def no_self_intersecting_polygon_during_extrusion(cls, val, values):
"""In this simple polyslab, we don't support self-intersecting polygons yet, meaning that
any normal cross section of the PolySlab cannot be self-intersecting. This part checks
if any self-interction will occur during extrusion with non-zero sidewall angle.
There are two types of self-intersection, known as edge events,
that can occur during dilation:
1) neighboring vertex-vertex crossing. This type of edge event can be treated with
``ComplexPolySlab`` which divides the polyslab into a list of simple polyslabs.
2) other types of edge events that can create holes/islands or split polygons.
To detect this, we sample _N_SAMPLE_POLYGON_INTERSECT cross sections to see if any creation
of polygons/holes, and changes in vertices number.
"""
# no need to validate anything here
if math.isclose(values["sidewall_angle"], 0):
return val
# apply dilation
poly_ref = PolySlab._proper_vertices(val)
if not math.isclose(values["dilation"], 0):
poly_ref = PolySlab._shift_vertices(poly_ref, values["dilation"])[0]
poly_ref = PolySlab._heal_polygon(poly_ref)
slab_min, slab_max = values["slab_bounds"]
slab_bounds = [getval(slab_min), getval(slab_max)]
# first, check vertex-vertex crossing at any point during extrusion
length = slab_bounds[1] - slab_bounds[0]
dist = [-length * np.tan(values["sidewall_angle"])]
# reverse the dilation value if it's defined on the top
if values["reference_plane"] == "top":
dist = [-dist[0]]
# for middle, both direction needs to be examined
elif values["reference_plane"] == "middle":
dist = [dist[0] / 2, -dist[0] / 2]
# capture vertex crossing events
max_thick = []
for dist_val in dist:
max_dist = PolySlab._neighbor_vertices_crossing_detection(poly_ref, dist_val)
if max_dist is not None:
max_thick.append(max_dist / abs(dist_val) * length)
if len(max_thick) > 0:
max_thick = min(max_thick)
raise SetupError(
"Sidewall angle or structure thickness is so large that the polygon "
"is self-intersecting during extrusion. "
f"Please either reduce structure thickness to be < {max_thick:.3e}, "
"or use our plugin 'ComplexPolySlab' to divide the complex polyslab "
"into a list of simple polyslabs."
)
# vertex-edge crossing event.
for dist_val in dist:
if PolySlab._edge_events_detection(poly_ref, dist_val):
raise SetupError(
"Sidewall angle or structure thickness is too large, "
"resulting in polygon splitting or generation of holes/islands. "
"A general treatment to self-intersecting polygon will be available "
"in future releases."
)
return val
[docs]
@classmethod
def from_gds(
cls,
gds_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",
) -> list[PolySlab]:
"""Import :class:`PolySlab` from a ``gdstk.Cell``.
Parameters
----------
gds_cell : gdstk.Cell
``gdstk.Cell`` containing 2D geometric data.
axis : int
Integer index into the polygon's slab axis. (0,1,2) -> (x,y,z).
slab_bounds: Tuple[float, float]
Minimum and maximum positions of the slab along ``axis``.
gds_layer : int
Layer index in the ``gds_cell``.
gds_dtype : int = None
Data-type index in the ``gds_cell``.
If ``None``, imports all data for this layer into the returned list.
gds_scale : float = 1.0
Length scale used in GDS file in units of MICROMETER.
For example, if gds file uses nanometers, set ``gds_scale=1e-3``.
Must be positive.
dilation : float = 0.0
Dilation of the polygon in the base by shifting each edge along its
normal outwards direction by a distance;
a negative value corresponds to erosion.
sidewall_angle : float = 0
Angle of the sidewall.
``sidewall_angle=0`` (default) specifies vertical wall,
while ``0<sidewall_angle<np.pi/2`` for the base to be larger than the top.
reference_plane : PlanePosition = "middle"
The position of the GDS layer. It can be at the ``bottom``, ``middle``,
or ``top`` of the PolySlab. E.g. if ``axis=1``, ``bottom`` refers to the
negative side of y-axis, and ``top`` refers to the positive side of y-axis.
Returns
-------
List[:class:`PolySlab`]
List of :class:`PolySlab` objects sharing ``axis`` and slab bound properties.
"""
all_vertices = PolySlab._load_gds_vertices(gds_cell, gds_layer, gds_dtype, gds_scale)
return [
cls(
vertices=verts,
axis=axis,
slab_bounds=slab_bounds,
dilation=dilation,
sidewall_angle=sidewall_angle,
reference_plane=reference_plane,
)
for verts in all_vertices
]
@staticmethod
def _load_gds_vertices(
gds_cell,
gds_layer: int,
gds_dtype: Optional[int] = None,
gds_scale: pydantic.PositiveFloat = 1.0,
) -> list[ArrayFloat2D]:
"""Import :class:`PolySlab` from a ``gdstk.Cell``.
Parameters
----------
gds_cell : gdstk.Cell
``gdstk.Cell`` containing 2D geometric data.
gds_layer : int
Layer index in the ``gds_cell``.
gds_dtype : int = None
Data-type index in the ``gds_cell``.
If ``None``, imports all data for this layer into the returned list.
gds_scale : float = 1.0
Length scale used in GDS file in units of MICROMETER.
For example, if gds file uses nanometers, set ``gds_scale=1e-3``.
Must be positive.
Returns
-------
List[ArrayFloat2D]
List of :class:`.ArrayFloat2D`
"""
import gdstk
gds_cell_class_name = str(gds_cell.__class__)
if not isinstance(gds_cell, gdstk.Cell):
if (
"gdstk" in gds_cell_class_name
): # Check if it might be a gdstk cell but gdstk is not found
raise Tidy3dImportError(
"Module 'gdstk' not found. It is required to import gdstk cells."
)
raise ValueError(
f"validate 'gds_cell' of type '{gds_cell_class_name}' "
"does not seem to be associated with 'gdstk' package "
"and therefore can't be loaded by Tidy3D."
)
all_vertices = base.Geometry.load_gds_vertices_gdstk(
gds_cell=gds_cell,
gds_layer=gds_layer,
gds_dtype=gds_dtype,
gds_scale=gds_scale,
)
# convert vertices into polyslabs
polygons = [PolySlab.make_shapely_polygon(vertices).buffer(0) for vertices in all_vertices]
polys_union = shapely.unary_union(polygons, grid_size=base.POLY_GRID_SIZE)
if polys_union.geom_type == "Polygon":
all_vertices = [np.array(polys_union.exterior.coords)]
elif polys_union.geom_type == "MultiPolygon":
all_vertices = [np.array(polygon.exterior.coords) for polygon in polys_union.geoms]
return all_vertices
@property
def center_axis(self) -> float:
"""Gets the position of the center of the geometry in the out of plane dimension."""
zmin, zmax = self.slab_bounds
if np.isneginf(zmin) and np.isposinf(zmax):
return 0.0
zmin = max(zmin, -LARGE_NUMBER)
zmax = min(zmax, LARGE_NUMBER)
return (zmax + zmin) / 2.0
@property
def length_axis(self) -> float:
"""Gets the length of the geometry along the out of plane dimension."""
zmin, zmax = self.slab_bounds
return zmax - zmin
@property
def finite_length_axis(self) -> float:
"""Gets the length of the PolySlab along the out of plane dimension.
First clips the slab bounds to LARGE_NUMBER and then returns difference.
"""
zmin, zmax = self.slab_bounds
zmin = max(zmin, -LARGE_NUMBER)
zmax = min(zmax, LARGE_NUMBER)
return zmax - zmin
@cached_property
def reference_polygon(self) -> np.ndarray:
"""The polygon at the reference plane.
Returns
-------
ArrayLike[float, float]
The vertices of the polygon at the reference plane.
"""
vertices = self._proper_vertices(self.vertices)
if math.isclose(self.dilation, 0):
return vertices
offset_vertices = self._shift_vertices(vertices, self.dilation)[0]
return self._heal_polygon(offset_vertices)
@cached_property
def middle_polygon(self) -> np.ndarray:
"""The polygon at the middle.
Returns
-------
ArrayLike[float, float]
The vertices of the polygon at the middle.
"""
dist = self._extrusion_length_to_offset_distance(self.finite_length_axis / 2)
if self.reference_plane == "bottom":
return self._shift_vertices(self.reference_polygon, dist)[0]
if self.reference_plane == "top":
return self._shift_vertices(self.reference_polygon, -dist)[0]
# middle case
return self.reference_polygon
@cached_property
def base_polygon(self) -> np.ndarray:
"""The polygon at the base, derived from the ``middle_polygon``.
Returns
-------
ArrayLike[float, float]
The vertices of the polygon at the base.
"""
if self.reference_plane == "bottom":
return self.reference_polygon
dist = self._extrusion_length_to_offset_distance(-self.finite_length_axis / 2)
return self._shift_vertices(self.middle_polygon, dist)[0]
@cached_property
def top_polygon(self) -> np.ndarray:
"""The polygon at the top, derived from the ``middle_polygon``.
Returns
-------
ArrayLike[float, float]
The vertices of the polygon at the top.
"""
if self.reference_plane == "top":
return self.reference_polygon
dist = self._extrusion_length_to_offset_distance(self.finite_length_axis / 2)
return self._shift_vertices(self.middle_polygon, dist)[0]
@cached_property
def _normal_2dmaterial(self) -> Axis:
"""Get the normal to the given geometry, checking that it is a 2D geometry."""
if self.slab_bounds[0] != self.slab_bounds[1]:
raise ValidationError("'Medium2D' requires the 'PolySlab' bounds to be equal.")
return self.axis
def _update_from_bounds(self, bounds: tuple[float, float], axis: Axis) -> PolySlab:
"""Returns an updated geometry which has been transformed to fit within ``bounds``
along the ``axis`` direction."""
if axis != self.axis:
raise ValueError(
f"'_update_from_bounds' may only be applied along axis '{self.axis}', "
f"but was given axis '{axis}'."
)
return self.updated_copy(slab_bounds=bounds)
@cached_property
def is_ccw(self) -> bool:
"""Is this ``PolySlab`` CCW-oriented?"""
return PolySlab._area(self.vertices) > 0
[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:`Geometry`, and ``False`` otherwise.
Note
----
For slanted sidewalls, this function only works if x, y, and z are arrays produced by a
``meshgrid call``, i.e. 3D arrays and each is constant along one axis.
Parameters
----------
x : np.ndarray[float]
Array of point positions in x direction.
y : np.ndarray[float]
Array of point positions in y direction.
z : np.ndarray[float]
Array of point positions in z direction.
Returns
-------
np.ndarray[bool]
``True`` for every point that is inside the geometry.
"""
self._ensure_equal_shape(x, y, z)
z, (x, y) = self.pop_axis((x, y, z), axis=self.axis)
z0 = self.center_axis
dist_z = np.abs(z - z0)
inside_height = dist_z <= (self.finite_length_axis / 2)
# avoid going into face checking if no points are inside slab bounds
if not np.any(inside_height):
return inside_height
# check what points are inside polygon cross section (face)
z_local = z - z0 # distance to the middle
dist = -z_local * self._tanq
if isinstance(x, np.ndarray):
inside_polygon = np.zeros_like(inside_height)
xs_slab = x[inside_height]
ys_slab = y[inside_height]
# vertical sidewall
if math.isclose(self.sidewall_angle, 0):
face_polygon = shapely.Polygon(self.reference_polygon).buffer(fp_eps)
shapely.prepare(face_polygon)
inside_polygon_slab = shapely.contains_xy(face_polygon, x=xs_slab, y=ys_slab)
inside_polygon[inside_height] = inside_polygon_slab
# slanted sidewall, offsetting vertices at each z
else:
# a helper function for moving axis
def _move_axis(arr):
return np.moveaxis(arr, source=self.axis, destination=-1)
def _move_axis_reverse(arr):
return np.moveaxis(arr, source=-1, destination=self.axis)
inside_polygon_axis = _move_axis(inside_polygon)
x_axis = _move_axis(x)
y_axis = _move_axis(y)
for z_i in range(z.shape[self.axis]):
if not _move_axis(inside_height)[0, 0, z_i]:
continue
vertices_z = self._shift_vertices(
self.middle_polygon, _move_axis(dist)[0, 0, z_i]
)[0]
face_polygon = shapely.Polygon(vertices_z).buffer(fp_eps)
shapely.prepare(face_polygon)
xs = x_axis[:, :, 0].flatten()
ys = y_axis[:, :, 0].flatten()
inside_polygon_slab = shapely.contains_xy(face_polygon, x=xs, y=ys)
inside_polygon_axis[:, :, z_i] = inside_polygon_slab.reshape(x_axis.shape[:2])
inside_polygon = _move_axis_reverse(inside_polygon_axis)
else:
vertices_z = self._shift_vertices(self.middle_polygon, dist)[0]
face_polygon = self.make_shapely_polygon(vertices_z).buffer(fp_eps)
point = shapely.Point(x, y)
inside_polygon = face_polygon.covers(point)
return inside_height * inside_polygon
@verify_packages_import(["trimesh"])
def _do_intersections_tilted_plane(
self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4
) -> 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.
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 len(self.base_polygon) > _MAX_POLYSLAB_VERTICES_FOR_TRIANGULATION:
log.warning(
f"Processing PolySlabs with over {_MAX_POLYSLAB_VERTICES_FOR_TRIANGULATION} vertices can be slow.",
log_once=True,
)
base_triangles = triangulation.triangulate(self.base_polygon)
top_triangles = (
base_triangles
if math.isclose(self.sidewall_angle, 0)
else triangulation.triangulate(self.top_polygon)
)
n = len(self.base_polygon)
faces = (
[[a, b, c] for c, b, a in base_triangles]
+ [[n + a, n + b, n + c] for a, b, c in top_triangles]
+ [(i, (i + 1) % n, n + i) for i in range(n)]
+ [((i + 1) % n, n + ((i + 1) % n), n + i) for i in range(n)]
)
x = np.hstack((self.base_polygon[:, 0], self.top_polygon[:, 0]))
y = np.hstack((self.base_polygon[:, 1], self.top_polygon[:, 1]))
z = np.hstack((np.full(n, self.slab_bounds[0]), np.full(n, self.slab_bounds[1])))
vertices = np.vstack(self.unpop_axis(z, (x, y), self.axis)).T
mesh = trimesh.Trimesh(vertices, faces)
section = mesh.section(plane_origin=origin, plane_normal=normal)
if section is None:
return []
path, _ = section.to_2D(to_2D=to_2D)
return path.polygons_full
def _intersections_normal(self, z: float):
"""Find shapely geometries intersecting planar geometry with axis normal to slab.
Parameters
----------
z : float
Position along the axis normal to slab.
Returns
-------
List[shapely.geometry.base.BaseGeometry]
List of 2D shapes that intersect plane.
For more details refer to
`Shapely's Documentation <https://shapely.readthedocs.io/en/stable/project.html>`_.
"""
if math.isclose(self.sidewall_angle, 0):
return [self.make_shapely_polygon(self.reference_polygon)]
z0 = self.center_axis
z_local = z - z0 # distance to the middle
dist = -z_local * self._tanq
vertices_z = self._shift_vertices(self.middle_polygon, dist)[0]
return [self.make_shapely_polygon(vertices_z)]
def _intersections_side(self, position, axis) -> list:
"""Find shapely geometries intersecting planar geometry with axis orthogonal to slab.
For slanted polyslab, the procedure is as follows,
1) Find out all z-coordinates where the plane will intersect directly with a vertex.
Denote the coordinates as (z_0, z_1, z_2, ... )
2) Find out all polygons that can be formed between z_i and z_{i+1}. There are two
types of polygons:
a) formed by the plane intersecting the edges
b) formed by the plane intersecting the vertices.
For either type, one needs to compute:
i) intersecting position
ii) angle between the plane and the intersecting edge
For a), both are straightforward to compute; while for b), one needs to compute
which edge the plane will slide into.
3) Looping through z_i, and merge all polygons. The partition by z_i is because once
the plane intersects the vertex, it can intersect with other edges during
the extrusion.
Parameters
----------
position : float
Position along ``axis``.
axis : int
Integer index into 'xyz' (0,1,2).
Returns
-------
List[shapely.geometry.base.BaseGeometry]
List of 2D shapes that intersect plane.
For more details refer to
`Shapely's Documentation <https://shapely.readthedocs.io/en/stable/project.html>`_.
"""
# find out all z_i where the plane will intersect the vertex
z0 = self.center_axis
z_base = z0 - self.finite_length_axis / 2
axis_ordered = self._order_axis(axis)
height_list = self._find_intersecting_height(position, axis_ordered)
polys = []
# looping through z_i to assemble the polygons
height_list = np.append(height_list, self.finite_length_axis)
h_base = 0.0
for h_top in height_list:
# length within between top and bottom
h_length = h_top - h_base
# coordinate of each subsection
z_min = z_base + h_base
z_max = np.inf if np.isposinf(h_top) else z_base + h_top
# for vertical sidewall, no need for complications
if math.isclose(self.sidewall_angle, 0):
ints_y, ints_angle = self._find_intersecting_ys_angle_vertical(
self.reference_polygon, position, axis_ordered
)
else:
# for slanted sidewall, move up by `fp_eps` in case vertices are degenerate at the base.
dist = -(h_base - self.finite_length_axis / 2 + fp_eps) * self._tanq
vertices = self._shift_vertices(self.middle_polygon, dist)[0]
ints_y, ints_angle = self._find_intersecting_ys_angle_slant(
vertices, position, axis_ordered
)
# make polygon with intersections and z axis information
for y_index in range(len(ints_y) // 2):
y_min = ints_y[2 * y_index]
y_max = ints_y[2 * y_index + 1]
minx, miny = self._order_by_axis(plane_val=y_min, axis_val=z_min, axis=axis)
maxx, maxy = self._order_by_axis(plane_val=y_max, axis_val=z_max, axis=axis)
if math.isclose(self.sidewall_angle, 0):
polys.append(self.make_shapely_box(minx, miny, maxx, maxy))
else:
angle_min = ints_angle[2 * y_index]
angle_max = ints_angle[2 * y_index + 1]
angle_min = np.arctan(np.tan(self.sidewall_angle) / np.sin(angle_min))
angle_max = np.arctan(np.tan(self.sidewall_angle) / np.sin(angle_max))
dy_min = h_length * np.tan(angle_min)
dy_max = h_length * np.tan(angle_max)
x1, y1 = self._order_by_axis(plane_val=y_min, axis_val=z_min, axis=axis)
x2, y2 = self._order_by_axis(plane_val=y_max, axis_val=z_min, axis=axis)
x3, y3 = self._order_by_axis(
plane_val=y_max - dy_max, axis_val=z_max, axis=axis
)
x4, y4 = self._order_by_axis(
plane_val=y_min + dy_min, axis_val=z_max, axis=axis
)
vertices = ((x1, y1), (x2, y2), (x3, y3), (x4, y4))
polys.append(self.make_shapely_polygon(vertices).buffer(0))
# update the base coordinate for the next subsection
h_base = h_top
# merge touching polygons
polys_union = shapely.unary_union(polys, grid_size=base.POLY_GRID_SIZE)
if polys_union.geom_type == "Polygon":
return [polys_union]
if polys_union.geom_type == "MultiPolygon":
return polys_union.geoms
# in other cases, just return the original unmerged polygons
return polys
def _find_intersecting_height(self, position: float, axis: int) -> np.ndarray:
"""Found a list of height where the plane will intersect with the vertices;
For vertical sidewall, just return np.array([]).
Assumes axis is handles so this function works on xy plane.
Parameters
----------
position : float
position along axis.
axis : int
Integer index into 'xyz' (0,1,2).
Returns
np.ndarray
Height (relative to the base) where the plane will intersect with vertices.
"""
if math.isclose(self.sidewall_angle, 0):
return np.array([])
# shift rate
dist = 1.0
shift_x, shift_y = PolySlab._shift_vertices(self.middle_polygon, dist)[2]
shift_val = shift_x if axis == 0 else shift_y
shift_val[np.isclose(shift_val, 0, rtol=_IS_CLOSE_RTOL)] = np.inf # for static vertices
# distance to the plane in the direction of vertex shifting
distance = self.middle_polygon[:, axis] - position
height = distance / self._tanq / shift_val + self.finite_length_axis / 2
height = np.unique(height)
# further filter very close ones
is_not_too_close = np.insert((np.diff(height) > fp_eps), 0, True)
height = height[is_not_too_close]
height = height[height > fp_eps]
height = height[height < self.finite_length_axis - fp_eps]
return height
def _find_intersecting_ys_angle_vertical(
self,
vertices: np.ndarray,
position: float,
axis: int,
exclude_on_vertices: bool = False,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Finds pairs of forward and backwards vertices where polygon intersects position at axis,
Find intersection point (in y) assuming straight line,and intersecting angle between plane
and edges. (For unslanted polyslab).
Assumes axis is handles so this function works on xy plane.
Parameters
----------
vertices : np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
position : float
position along axis.
axis : int
Integer index into 'xyz' (0,1,2).
exclude_on_vertices : bool = False
Whether to exclude those intersecting directly with the vertices.
Returns
-------
Union[np.ndarray, np.ndarray]
List of intersection points along y direction.
List of angles between plane and edges.
"""
vertices_axis = vertices
# flip vertices x,y for axis = y
if axis == 1:
vertices_axis = np.roll(vertices_axis, shift=1, axis=1)
# get the forward vertices
vertices_f = np.roll(vertices_axis, shift=-1, axis=0)
# x coordinate of the two sets of vertices
x_vertices_f, _ = vertices_f.T
x_vertices_axis, _ = vertices_axis.T
# find which segments intersect
f_left_to_intersect = x_vertices_f <= position
orig_right_to_intersect = x_vertices_axis > position
intersects_b = np.logical_and(f_left_to_intersect, orig_right_to_intersect)
f_right_to_intersect = x_vertices_f > position
orig_left_to_intersect = x_vertices_axis <= position
intersects_f = np.logical_and(f_right_to_intersect, orig_left_to_intersect)
# exclude vertices at the position if exclude_on_vertices is True
if exclude_on_vertices:
intersects_on = np.isclose(x_vertices_axis, position, rtol=_IS_CLOSE_RTOL)
intersects_f_on = np.isclose(x_vertices_f, position, rtol=_IS_CLOSE_RTOL)
intersects_both_off = np.logical_not(np.logical_or(intersects_on, intersects_f_on))
intersects_f &= intersects_both_off
intersects_b &= intersects_both_off
intersects_segment = np.logical_or(intersects_b, intersects_f)
iverts_b = vertices_axis[intersects_segment]
iverts_f = vertices_f[intersects_segment]
# intersecting positions and angles
ints_y = []
ints_angle = []
for vertices_f_local, vertices_b_local in zip(iverts_b, iverts_f):
x1, y1 = vertices_f_local
x2, y2 = vertices_b_local
slope = (y2 - y1) / (x2 - x1)
y = y1 + slope * (position - x1)
ints_y.append(y)
ints_angle.append(np.pi / 2 - np.arctan(np.abs(slope)))
ints_y = np.array(ints_y)
ints_angle = np.array(ints_angle)
sort_index = np.argsort(ints_y)
ints_y_sort = ints_y[sort_index]
ints_angle_sort = ints_angle[sort_index]
return ints_y_sort, ints_angle_sort
def _find_intersecting_ys_angle_slant(
self, vertices: np.ndarray, position: float, axis: int
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Finds pairs of forward and backwards vertices where polygon intersects position at axis,
Find intersection point (in y) assuming straight line,and intersecting angle between plane
and edges. (For slanted polyslab)
Assumes axis is handles so this function works on xy plane.
Parameters
----------
vertices : np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
position : float
position along axis.
axis : int
Integer index into 'xyz' (0,1,2).
Returns
-------
Union[np.ndarray, np.ndarray]
List of intersection points along y direction.
List of angles between plane and edges.
"""
vertices_axis = vertices.copy()
# flip vertices x,y for axis = y
if axis == 1:
vertices_axis = np.roll(vertices_axis, shift=1, axis=1)
# get the forward vertices
vertices_f = np.roll(vertices_axis, shift=-1, axis=0)
# get the backward vertices
vertices_b = np.roll(vertices_axis, shift=1, axis=0)
## First part, plane intersects with edges, same as vertical
ints_y, ints_angle = self._find_intersecting_ys_angle_vertical(
vertices, position, axis, exclude_on_vertices=True
)
ints_y = ints_y.tolist()
ints_angle = ints_angle.tolist()
## Second part, plane intersects directly with vertices
# vertices on the intersection
intersects_on = np.isclose(vertices_axis[:, 0], position, rtol=_IS_CLOSE_RTOL)
iverts_on = vertices_axis[intersects_on]
# position of the neighbouring vertices
iverts_b = vertices_b[intersects_on]
iverts_f = vertices_f[intersects_on]
# shift rate
dist = -np.sign(self.sidewall_angle)
shift_x, shift_y = self._shift_vertices(self.middle_polygon, dist)[2]
shift_val = shift_x if axis == 0 else shift_y
shift_val = shift_val[intersects_on]
for vertices_f_local, vertices_b_local, vertices_on_local, shift_local in zip(
iverts_f, iverts_b, iverts_on, shift_val
):
x_on, y_on = vertices_on_local
x_f, y_f = vertices_f_local
x_b, y_b = vertices_b_local
num_added = 0 # keep track the number of added vertices
slope = [] # list of slopes for added vertices
# case 1, shifting velocity is 0
if np.isclose(shift_local, 0, rtol=_IS_CLOSE_RTOL):
ints_y.append(y_on)
# Slope w.r.t. forward and backward should equal,
# just pick one of them.
slope.append((y_on - y_b) / (x_on - x_b))
ints_angle.append(np.pi / 2 - np.arctan(np.abs(slope[0])))
continue
# case 2, shifting towards backward direction
if (x_b - position) * shift_local < 0:
ints_y.append(y_on)
slope.append((y_on - y_b) / (x_on - x_b))
num_added += 1
# case 3, shifting towards forward direction
if (x_f - position) * shift_local < 0:
ints_y.append(y_on)
slope.append((y_on - y_f) / (x_on - x_f))
num_added += 1
# in case 2, and case 3, if just num_added = 1
if num_added == 1:
ints_angle.append(np.pi / 2 - np.arctan(np.abs(slope[0])))
# if num_added = 2, the order of the two new vertices needs to handled correctly;
# it should be sorted according to the -slope * moving direction
elif num_added == 2:
dressed_slope = [-s_i * shift_local for s_i in slope]
sort_index = np.argsort(np.array(dressed_slope))
sorted_slope = np.array(slope)[sort_index]
ints_angle.append(np.pi / 2 - np.arctan(np.abs(sorted_slope[0])))
ints_angle.append(np.pi / 2 - np.arctan(np.abs(sorted_slope[1])))
ints_y = np.array(ints_y)
ints_angle = np.array(ints_angle)
sort_index = np.argsort(ints_y)
ints_y_sort = ints_y[sort_index]
ints_angle_sort = ints_angle[sort_index]
return ints_y_sort, ints_angle_sort
@cached_property
def bounds(self) -> Bound:
"""Returns bounding box min and max coordinates. The dilation and slant angle are not
taken into account exactly for speed. Instead, the polygon may be slightly smaller than
the returned bounds, but it should always be fully contained.
Returns
-------
Tuple[float, float, float], Tuple[float, float float]
Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
"""
# check for the maximum possible contribution from dilation/slant on each side
max_offset = self.dilation
if not math.isclose(self.sidewall_angle, 0):
if self.reference_plane == "bottom":
max_offset += max(0, -self._tanq * self.finite_length_axis)
elif self.reference_plane == "top":
max_offset += max(0, self._tanq * self.finite_length_axis)
elif self.reference_plane == "middle":
max_offset += max(0, abs(self._tanq) * self.finite_length_axis / 2)
# special care when dilated
if max_offset > 0:
dilated_vertices = self._shift_vertices(
self._proper_vertices(self.vertices), max_offset
)[0]
xmin, ymin = np.amin(dilated_vertices, axis=0)
xmax, ymax = np.amax(dilated_vertices, axis=0)
else:
# otherwise, bounds are directly based on the supplied vertices
xmin, ymin = np.amin(self.vertices, axis=0)
xmax, ymax = np.amax(self.vertices, axis=0)
# get bounds in (local) z
zmin, zmax = self.slab_bounds
# rearrange axes
coords_min = self.unpop_axis(zmin, (xmin, ymin), axis=self.axis)
coords_max = self.unpop_axis(zmax, (xmax, ymax), axis=self.axis)
return (tuple(coords_min), tuple(coords_max))
def _extrusion_length_to_offset_distance(self, extrusion: float) -> float:
"""Convert extrusion length to offset distance."""
if math.isclose(self.sidewall_angle, 0):
return 0
return -extrusion * self._tanq
@staticmethod
def _area(vertices: np.ndarray) -> float:
"""Compute the signed polygon area (positive for CCW orientation).
Parameters
----------
vertices : np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
Returns
-------
float
Signed polygon area (positive for CCW orientation).
"""
vert_shift = np.roll(vertices, axis=0, shift=-1)
xs, ys = vertices.T
xs_shift, ys_shift = vert_shift.T
term1 = xs * ys_shift
term2 = ys * xs_shift
return np.sum(term1 - term2) * 0.5
@staticmethod
def _perimeter(vertices: np.ndarray) -> float:
"""Compute the polygon perimeter.
Parameters
----------
vertices : np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
Returns
-------
float
Polygon perimeter.
"""
vert_shift = np.roll(vertices, axis=0, shift=-1)
squared_diffs = (vertices - vert_shift) ** 2
# distance along each edge
dists = np.sqrt(squared_diffs.sum(axis=-1))
# total distance along all edges
return np.sum(dists)
@staticmethod
def _orient(vertices: np.ndarray) -> np.ndarray:
"""Return a CCW-oriented polygon.
Parameters
----------
vertices : np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
Returns
-------
np.ndarray
Vertices of a CCW-oriented polygon.
"""
return vertices if PolySlab._area(vertices) > 0 else vertices[::-1, :]
@staticmethod
def _remove_duplicate_vertices(vertices: np.ndarray) -> np.ndarray:
"""Remove redundant/identical nearest neighbour vertices.
Parameters
----------
vertices : np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
Returns
-------
np.ndarray
Vertices of polygon.
"""
vertices_f = np.roll(vertices, shift=-1, axis=0)
vertices_diff = np.linalg.norm(vertices - vertices_f, axis=1)
return vertices[~np.isclose(vertices_diff, 0, rtol=_IS_CLOSE_RTOL)]
@staticmethod
def _proper_vertices(vertices: ArrayFloat2D) -> np.ndarray:
"""convert vertices to np.array format,
removing duplicate neighbouring vertices,
and oriented in CCW direction.
Returns
-------
ArrayLike[float, float]
The vertices of the polygon for internal use.
"""
vertices_np = PolySlab.vertices_to_array(vertices)
return PolySlab._orient(PolySlab._remove_duplicate_vertices(vertices_np))
@staticmethod
def _edge_events_detection(
proper_vertices: np.ndarray, dilation: float, ignore_at_dist: bool = True
) -> bool:
"""Detect any edge events within the offset distance ``dilation``.
If ``ignore_at_dist=True``, the edge event at ``dist`` is ignored.
"""
# ignore the event that occurs right at the offset distance
if ignore_at_dist:
dilation -= fp_eps * dilation / abs(dilation)
# number of vertices before offsetting
num_vertices = proper_vertices.shape[0]
# 0) fully eroded?
if dilation < 0 and dilation < -PolySlab._maximal_erosion(proper_vertices):
return True
# sample at a few dilation values
dist_list = (
dilation
* np.linspace(0, 1, 1 + _N_SAMPLE_POLYGON_INTERSECT, dtype=GRADIENT_DTYPE_FLOAT)[1:]
)
for dist in dist_list:
# offset: we offset the vertices first, and then use shapely to make it proper
# in principle, one can offset with shapely.buffer directly, but shapely somehow
# automatically removes some vertices even though no change of topology.
poly_offset = PolySlab._shift_vertices(proper_vertices, dist)[0]
# flipped winding number
if PolySlab._area(poly_offset) < fp_eps**2:
return True
poly_offset = shapely.make_valid(PolySlab.make_shapely_polygon(poly_offset))
# 1) polygon split or create holes/islands
if not poly_offset.geom_type == "Polygon" or len(poly_offset.interiors) > 0:
return True
# 2) reduction in vertex number
offset_vertices = PolySlab._proper_vertices(list(poly_offset.exterior.coords))
if offset_vertices.shape[0] != num_vertices:
return True
# 3) some split polygon might fully disappear after the offset, but they
# can be detected if we offset back.
poly_offset_back = shapely.make_valid(
PolySlab.make_shapely_polygon(PolySlab._shift_vertices(offset_vertices, -dist)[0])
)
if poly_offset_back.geom_type == "MultiPolygon" or len(poly_offset_back.interiors) > 0:
return True
offset_back_vertices = list(poly_offset_back.exterior.coords)
if PolySlab._proper_vertices(offset_back_vertices).shape[0] != num_vertices:
return True
return False
@staticmethod
def _neighbor_vertices_crossing_detection(
vertices: np.ndarray, dist: float, ignore_at_dist: bool = True
) -> float:
"""Detect if neighboring vertices will cross after a dilation distance dist.
Parameters
----------
vertices : np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
dist : float
Distance to offset.
ignore_at_dist : bool, optional
whether to ignore the event right at ``dist`.
Returns
-------
float
the absolute value of the maximal allowed dilation
if there are any crossing, otherwise return ``None``.
"""
# ignore the event that occurs right at the offset distance
if ignore_at_dist:
dist -= fp_eps * dist / abs(dist)
edge_length, edge_reduction = PolySlab._edge_length_and_reduction_rate(vertices)
length_remaining = edge_length - edge_reduction * dist
if np.any(length_remaining < 0):
index_oversized = length_remaining < 0
max_dist = np.min(
np.abs(edge_length[index_oversized] / edge_reduction[index_oversized])
)
return max_dist
return None
[docs]
@staticmethod
def array_to_vertices(arr_vertices: np.ndarray) -> ArrayFloat2D:
"""Converts a numpy array of vertices to a list of tuples."""
return list(arr_vertices)
[docs]
@staticmethod
def vertices_to_array(vertices_tuple: ArrayFloat2D) -> np.ndarray:
"""Converts a list of tuples (vertices) to a numpy array."""
return np.array(vertices_tuple)
@cached_property
def interior_angle(self) -> ArrayFloat1D:
"""Angle formed inside polygon by two adjacent edges."""
def normalize(v):
return v / np.linalg.norm(v, axis=0)
vs_orig = self.reference_polygon.T
vs_next = np.roll(vs_orig, axis=-1, shift=-1)
vs_previous = np.roll(vs_orig, axis=-1, shift=+1)
asp = normalize(vs_next - vs_orig)
asm = normalize(vs_previous - vs_orig)
cos_angle = asp[0] * asm[0] + asp[1] * asm[1]
sin_angle = asp[0] * asm[1] - asp[1] * asm[0]
angle = np.arccos(cos_angle)
# concave angles
angle[sin_angle < 0] = 2 * np.pi - angle[sin_angle < 0]
return angle
@staticmethod
def _shift_vertices(
vertices: np.ndarray, dist
) -> tuple[np.ndarray, np.ndarray, tuple[np.ndarray, np.ndarray]]:
"""Shifts the vertices of a polygon outward uniformly by distances
`dists`.
Parameters
----------
np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
dist : float
Distance to offset.
Returns
-------
Tuple[np.ndarray, np.narray,Tuple[np.ndarray,np.ndarray]]
New polygon vertices;
and the shift of vertices in direction parallel to the edges.
Shift along x and y direction.
"""
if math.isclose(dist, 0):
return vertices, np.zeros(vertices.shape[0], dtype=float), None
def rot90(v):
"""90 degree rotation of 2d vector
vx -> vy
vy -> -vx
"""
vxs, vys = v
return np.stack((-vys, vxs), axis=0)
def cross(u, v):
return u[0] * v[1] - u[1] * v[0]
def normalize(v):
return v / np.linalg.norm(v, axis=0)
vs_orig = copy(vertices.T)
vs_next = np.roll(copy(vs_orig), axis=-1, shift=-1)
vs_previous = np.roll(copy(vs_orig), axis=-1, shift=+1)
asp = normalize(vs_next - vs_orig)
asm = normalize(vs_orig - vs_previous)
# the vertex shift is decomposed into parallel and perpendicular directions
perpendicular_shift = -dist
det = cross(asm, asp)
tan_half_angle = np.where(
np.isclose(det, 0, rtol=_IS_CLOSE_RTOL),
0.0,
cross(asm, rot90(asm - asp)) / (det + np.isclose(det, 0, rtol=_IS_CLOSE_RTOL)),
)
parallel_shift = dist * tan_half_angle
shift_total = perpendicular_shift * rot90(asm) + parallel_shift * asm
shift_x = shift_total[0, :]
shift_y = shift_total[1, :]
return (
np.swapaxes(vs_orig + shift_total, -2, -1),
parallel_shift,
(shift_x, shift_y),
)
@staticmethod
def _edge_length_and_reduction_rate(
vertices: np.ndarray,
) -> tuple[np.ndarray, np.ndarray]:
"""Edge length of reduction rate of each edge with unit offset length.
Parameters
----------
vertices : np.ndarray
Shape (N, 2) defining the polygon vertices in the xy-plane.
Returns
-------
Tuple[np.ndarray, np.narray]
edge length, and reduction rate
"""
# edge length
vs_orig = copy(vertices.T)
vs_next = np.roll(copy(vs_orig), axis=-1, shift=-1)
edge_length = np.linalg.norm(vs_next - vs_orig, axis=0)
# edge length remaining
dist = 1
parallel_shift = PolySlab._shift_vertices(vertices, dist)[1]
parallel_shift_p = np.roll(copy(parallel_shift), shift=-1)
edge_reduction = -(parallel_shift + parallel_shift_p)
return edge_length, edge_reduction
@staticmethod
def _maximal_erosion(vertices: np.ndarray) -> float:
"""The erosion value that reduces the length of
all edges to be non-positive.
"""
edge_length, edge_reduction = PolySlab._edge_length_and_reduction_rate(vertices)
ind_nonzero = abs(edge_reduction) > fp_eps
return -np.min(edge_length[ind_nonzero] / edge_reduction[ind_nonzero])
@staticmethod
def _heal_polygon(vertices: np.ndarray) -> np.ndarray:
"""heal a self-intersecting polygon."""
shapely_poly = PolySlab.make_shapely_polygon(vertices)
if shapely_poly.is_valid:
return vertices
if isbox(vertices):
raise NotImplementedError(
"The dilation caused damage to the polygon. "
"Automatically healing this is currently not supported when "
"differentiating w.r.t. the vertices. Try increasing the spacing "
"between vertices or reduce the amount of dilation."
)
# perform healing
poly_heal = shapely.make_valid(shapely_poly)
return PolySlab._proper_vertices(list(poly_heal.exterior.coords))
def _volume(self, bounds: Bound) -> float:
"""Returns object's volume within given bounds."""
z_min, z_max = self.slab_bounds
z_min = max(z_min, bounds[0][self.axis])
z_max = min(z_max, bounds[1][self.axis])
length = z_max - z_min
top_area = abs(self._area(self.top_polygon))
base_area = abs(self._area(self.base_polygon))
# https://mathworld.wolfram.com/PyramidalFrustum.html
return 1.0 / 3.0 * length * (top_area + base_area + np.sqrt(top_area * base_area))
def _surface_area(self, bounds: Bound) -> float:
"""Returns object's surface area within given bounds."""
area = 0
top_polygon = self.top_polygon
base_polygon = self.base_polygon
top_area = abs(self._area(top_polygon))
base_area = abs(self._area(base_polygon))
top_perim = self._perimeter(top_polygon)
base_perim = self._perimeter(base_polygon)
z_min, z_max = self.slab_bounds
if z_min < bounds[0][self.axis]:
z_min = bounds[0][self.axis]
else:
area += base_area
if z_max > bounds[1][self.axis]:
z_max = bounds[1][self.axis]
else:
area += top_area
length = z_max - z_min
area += 0.5 * (top_perim + base_perim) * length
return area
""" Autograd code """
def _compute_derivatives(self, derivative_info: DerivativeInfo) -> AutogradFieldMap:
"""
Return VJPs while handling several edge-cases:
- If the slab volume does not overlap the simulation, all grads are zero
(one warning is issued).
- Faces that lie completely outside the simulation give zero ``slab_bounds``
gradients; this includes the +/- inf cases.
- A 2d simulation collapses the surface integral to a line integral
"""
vjps: AutogradFieldMap = {}
sim_min, sim_max = map(np.asarray, derivative_info.bounds_intersect)
extents = sim_max - sim_min
is_2d = np.isclose(extents[self.axis], 0.0)
# early return if polyslab is not in simulation domain
slab_min, slab_max = self.slab_bounds
if (slab_max <= sim_min[self.axis]) or (slab_min >= sim_max[self.axis]):
log.warning(
"'PolySlab' lies completely outside the simulation domain.",
log_once=True,
)
for p in derivative_info.paths:
vjps[p] = np.zeros_like(self.vertices) if p == ("vertices",) else 0.0
return vjps
# create interpolators once for ALL derivative computations
# use provided interpolators if available to avoid redundant field data conversions
interpolators = derivative_info.interpolators or derivative_info.create_interpolators(
dtype=GRADIENT_DTYPE_FLOAT
)
for path in derivative_info.paths:
if path == ("vertices",):
vjps[path] = self._compute_derivative_vertices(
derivative_info, sim_min, sim_max, is_2d, interpolators
)
elif path[0] == "slab_bounds":
idx = path[1]
face_coord = self.slab_bounds[idx]
# face entirely outside -> gradient 0
if (
np.isinf(face_coord)
or face_coord < sim_min[self.axis]
or face_coord > sim_max[self.axis]
or is_2d
):
vjps[path] = 0.0
continue
v = self._compute_derivative_slab_bounds(derivative_info, idx, interpolators)
# outward-normal convention
if idx == 0:
v *= -1
vjps[path] = v
else:
raise ValueError(f"No derivative defined w.r.t. 'PolySlab' field '{path}'.")
return vjps
def _compute_derivative_slab_bounds(
self, derivative_info: DerivativeInfo, min_max_index: int, interpolators: dict
) -> float:
"""VJP for one of the two horizontal faces of a ``PolySlab``.
The face is discretized into a Cartesian grid of small planar patches
whose linear size does not exceed ``_VJP_SAMPLE_SPACING``. The adjoint surface
integral is evaluated on every retained patch; the resulting derivative
is split equally between the two vertices that bound the edge segment.
"""
# rmin/rmax over the geometry and simulation box
rmin, rmax = derivative_info.bounds_intersect
_, (r1_min, r2_min) = self.pop_axis(rmin, axis=self.axis)
_, (r1_max, r2_max) = self.pop_axis(rmax, axis=self.axis)
ax_val = self.slab_bounds[min_max_index]
# planar grid resolution, clipped to polygon bounding box
face_verts = self.base_polygon if min_max_index == 0 else self.top_polygon
face_poly = shapely.Polygon(face_verts).buffer(fp_eps)
# limit the patch grid to the face that lives inside the simulation box
poly_min_r1, poly_min_r2, poly_max_r1, poly_max_r2 = face_poly.bounds
r1_min = max(r1_min, poly_min_r1)
r1_max = min(r1_max, poly_max_r1)
r2_min = max(r2_min, poly_min_r2)
r2_max = min(r2_max, poly_max_r2)
if (r1_max <= r1_min) and (r2_max <= r2_min):
# the polygon does not intersect the current simulation slice
return 0.0
# re-compute the extents after clipping to the polygon bounds
extents = np.array([r1_max - r1_min, r2_max - r2_min])
# choose surface or line integral
integral_fun = (
self.compute_derivative_slab_bounds_line
if np.isclose(extents, 0).any()
else self.compute_derivative_slab_bounds_surface
)
return integral_fun(
derivative_info,
extents,
r1_min,
r1_max,
r2_min,
r2_max,
ax_val,
face_poly,
min_max_index,
interpolators,
)
[docs]
def compute_derivative_slab_bounds_line(
self,
derivative_info: DerivativeInfo,
extents: np.ndarray,
r1_min: float,
r1_max: float,
r2_min: float,
r2_max: float,
ax_val: float,
face_poly: shapely.Polygon,
min_max_index: int,
interpolators: dict,
) -> float:
"""Handle degenerate line cross-section case"""
line_dim = 1 if np.isclose(extents[0], 0) else 0
poly_min_r1, poly_min_r2, poly_max_r1, poly_max_r2 = face_poly.bounds
if line_dim == 0: # x varies, y is fixed
l_min = max(r1_min, poly_min_r1)
l_max = min(r1_max, poly_max_r1)
else: # y varies, x is fixed
l_min = max(r2_min, poly_min_r2)
l_max = min(r2_max, poly_max_r2)
length = l_max - l_min
if np.isclose(length, 0):
return 0.0
dx = derivative_info.adaptive_vjp_spacing()
n_seg = max(1, int(np.ceil(length / dx)))
coords = np.linspace(l_min, l_max, 2 * n_seg + 1, dtype=GRADIENT_DTYPE_FLOAT)[1::2]
# build XY coordinates and in-plane direction vectors
if line_dim == 0:
xy = np.column_stack((coords, np.full_like(coords, r2_min)))
dir_vec_plane = np.column_stack((np.ones_like(coords), np.zeros_like(coords)))
else:
xy = np.column_stack((np.full_like(coords, r1_min), coords))
dir_vec_plane = np.column_stack((np.zeros_like(coords), np.ones_like(coords)))
inside = shapely.contains_xy(face_poly, xy[:, 0], xy[:, 1])
if not inside.any():
return 0.0
xy = xy[inside]
dir_vec_plane = dir_vec_plane[inside]
n_pts = len(xy)
centers_xyz = self.unpop_axis_vect(np.full(n_pts, ax_val), xy)
areas = np.full(n_pts, length / n_seg) # patch length
normals_xyz = self.unpop_axis_vect(
np.full(n_pts, -1 if min_max_index == 0 else 1),
np.zeros_like(xy),
)
perps1_xyz = self.unpop_axis_vect(np.zeros(n_pts), dir_vec_plane)
perps2_xyz = self.unpop_axis_vect(np.zeros(n_pts), np.zeros_like(dir_vec_plane))
vjps = derivative_info.evaluate_gradient_at_points(
centers_xyz, normals_xyz, perps1_xyz, perps2_xyz, interpolators
)
return np.real(np.sum(vjps * areas)).item()
[docs]
def compute_derivative_slab_bounds_surface(
self,
derivative_info: DerivativeInfo,
extents: np.ndarray,
r1_min: float,
r1_max: float,
r2_min: float,
r2_max: float,
ax_val: float,
face_poly: shapely.Polygon,
min_max_index: int,
interpolators: dict,
) -> float:
"""2d surface integral on a Gauss quadrature grid"""
dx = derivative_info.adaptive_vjp_spacing()
# uniform grid would use n1 x n2 points
n1_uniform, n2_uniform = np.maximum(1, np.ceil(extents / dx).astype(int))
# use ~1/2 Gauss points in each direction for similar accuracy
n1 = max(2, n1_uniform // 2)
n2 = max(2, n2_uniform // 2)
g1, w1 = leggauss(n1)
g2, w2 = leggauss(n2)
coords1 = (0.5 * (r1_max - r1_min) * g1 + 0.5 * (r1_max + r1_min)).astype(
GRADIENT_DTYPE_FLOAT, copy=False
)
coords2 = (0.5 * (r2_max - r2_min) * g2 + 0.5 * (r2_max + r2_min)).astype(
GRADIENT_DTYPE_FLOAT, copy=False
)
r1_grid, r2_grid = np.meshgrid(coords1, coords2, indexing="ij")
r1_flat = r1_grid.flatten()
r2_flat = r2_grid.flatten()
pts = np.column_stack((r1_flat, r2_flat))
in_face = shapely.contains_xy(face_poly, pts[:, 0], pts[:, 1])
if not in_face.any():
return 0.0
xyz = self.unpop_axis_vect(
np.full(in_face.sum(), ax_val, dtype=GRADIENT_DTYPE_FLOAT), pts[in_face]
)
n_patches = xyz.shape[0]
normals_xyz = self.unpop_axis_vect(
np.full(n_patches, -1 if min_max_index == 0 else 1, dtype=GRADIENT_DTYPE_FLOAT),
np.zeros((n_patches, 2), dtype=GRADIENT_DTYPE_FLOAT),
)
perps1_xyz = self.unpop_axis_vect(
np.zeros(n_patches, dtype=GRADIENT_DTYPE_FLOAT),
np.column_stack(
(
np.ones(n_patches, dtype=GRADIENT_DTYPE_FLOAT),
np.zeros(n_patches, dtype=GRADIENT_DTYPE_FLOAT),
)
),
)
perps2_xyz = self.unpop_axis_vect(
np.zeros(n_patches, dtype=GRADIENT_DTYPE_FLOAT),
np.column_stack(
(
np.zeros(n_patches, dtype=GRADIENT_DTYPE_FLOAT),
np.ones(n_patches, dtype=GRADIENT_DTYPE_FLOAT),
)
),
)
w1_grid, w2_grid = np.meshgrid(w1, w2, indexing="ij")
weights_flat = (w1_grid * w2_grid).flatten()[in_face]
jacobian = 0.25 * (r1_max - r1_min) * (r2_max - r2_min)
# area-based correction for non-rectangular domains (e.g. concave polygon)
# for constant integrand, integral should equal polygon area
sum_weights = np.sum(weights_flat)
if sum_weights > 0:
area_correction = face_poly.area / (sum_weights * jacobian)
weights_flat = weights_flat * area_correction
vjps = derivative_info.evaluate_gradient_at_points(
xyz, normals_xyz, perps1_xyz, perps2_xyz, interpolators
)
return np.real(np.sum(vjps * weights_flat * jacobian)).item()
def _compute_derivative_vertices(
self,
derivative_info: DerivativeInfo,
sim_min: np.ndarray,
sim_max: np.ndarray,
is_2d: bool = False,
interpolators: Optional[dict] = None,
) -> np.ndarray:
"""VJP for the vertices of a ``PolySlab``.
Each side-wall is sliced in the **in-plane** direction (along the edge)
and, when applicable, in the **out-of-plane** direction (along the
extrusion axis) so that every rectangular surface patch is no larger
than ``_VJP_SAMPLE_SPACING`` in any dimension. The adjoint surface
integral is evaluated on every retained patch; the resulting force is
split (weighted) between the two vertices that bound the edge segment.
Special cases:
- Pure 2d simulations - integral collapses to a line; only one z-slice is taken.
- Partial overlap with the simulation box - patches falling
outside ``[sim_min, sim_max]`` are skipped.
- Degenerate edges (zero length) - ignored.
"""
vertices = np.asarray(self.vertices, dtype=GRADIENT_DTYPE_FLOAT)
next_v = np.roll(vertices, -1, axis=0)
edges = next_v - vertices
basis = self.edge_basis_vectors(edges)
dx = derivative_info.adaptive_vjp_spacing()
# compute z‐slices
z0 = max(self.slab_bounds[0], sim_min[self.axis])
z1 = min(self.slab_bounds[1], sim_max[self.axis])
# early return if no z-slices
if (not is_2d) and z1 <= z0:
return np.zeros_like(vertices)
if is_2d:
z_centers = np.array([self.center_axis], dtype=GRADIENT_DTYPE_FLOAT)
dz_surf = 1.0
else:
n_z = max(1, int(np.ceil((z1 - z0) / dx)))
dz = (z1 - z0) / n_z
dz_surf = dz / np.cos(self.sidewall_angle)
z_centers = np.linspace(z0 + dz / 2, z1 - dz / 2, n_z, dtype=GRADIENT_DTYPE_FLOAT)
vjp_per_vertex = np.zeros_like(vertices, dtype=float)
normals_2d = np.delete(basis["norm"], self.axis, axis=1)
# use provided interpolators or create them if not provided
if interpolators is None:
interpolators = derivative_info.create_interpolators(dtype=GRADIENT_DTYPE_FLOAT)
def compute_edge_clip_bounds(v0_3d, v1_3d, sim_min, sim_max):
"""
Compute parametric bounds [t_start, t_end] for the portion of edge
(v0_3d -> v1_3d) that lies within [sim_min, sim_max].
Returns:
- (t_start, t_end): parametric bounds, or None if edge is fully outside
"""
t_start, t_end = 0.0, 1.0
for dim in range(3):
v0_d, v1_d = v0_3d[dim], v1_3d[dim]
min_d, max_d = sim_min[dim], sim_max[dim]
if np.isclose(v0_d, v1_d):
if v0_d < min_d or v0_d > max_d:
return None
continue
t_min = (min_d - v0_d) / (v1_d - v0_d)
t_max = (max_d - v0_d) / (v1_d - v0_d)
if t_min > t_max:
t_min, t_max = t_max, t_min
t_start = max(t_start, t_min)
t_end = min(t_end, t_max)
if t_start >= t_end:
return None
# avoid numerical issues at boundaries
if t_end <= t_start + EDGE_CLIP_TOLERANCE:
return None
return (t_start, t_end)
# estimate total number of patches for pre-allocation
# this avoids repeated list.extend calls
estimated_patches = 0
n_edges = len(vertices)
for ei in range(n_edges):
v0, v1 = vertices[ei], next_v[ei]
L = np.linalg.norm(v1 - v0)
if not np.isclose(L, 0.0):
# estimate samples per edge (conservative estimate)
n_samples = max(1, int(np.ceil(L / dx) * 0.4)) # 40% of uniform for Gauss
estimated_patches += n_samples * len(z_centers)
# pre-allocate arrays with estimated size
# over-allocate by 20% to handle edge cases, will trim later
estimated_patches = int(estimated_patches * 1.2)
all_centers = np.empty((estimated_patches, 3), dtype=GRADIENT_DTYPE_FLOAT)
all_areas = np.empty(estimated_patches, dtype=GRADIENT_DTYPE_FLOAT)
all_normals = np.empty((estimated_patches, 3), dtype=GRADIENT_DTYPE_FLOAT)
all_perps1 = np.empty((estimated_patches, 3), dtype=GRADIENT_DTYPE_FLOAT)
all_perps2 = np.empty((estimated_patches, 3), dtype=GRADIENT_DTYPE_FLOAT)
all_edge_info = [] # keep as list since it stores tuples
# track actual number of patches
patch_idx = 0
def get_adaptive_samples(L, dx, t_start=0.0, t_end=1.0):
"""Get sampling points and weights along edge using optimal quadrature.
Parameters:
- L: Full edge length
- dx: Target spacing
- t_start, t_end: Parametric bounds for clipped edges (0 <= t_start < t_end <= 1)
Returns:
- s: Parametric coordinates in [t_start, t_end]
- weights: Quadrature weights (sum to t_end - t_start)
"""
# compute effective length for the clipped portion
L_effective = L * (t_end - t_start)
n_uniform = max(1, int(np.ceil(L_effective / dx)))
# Gauss quadrature is more accurate than uniform sampling
if n_uniform <= 3:
n_gauss = n_uniform # for very short edges, keep same number
else:
n_gauss = max(
2, int(n_uniform * QUAD_SAMPLE_FRACTION)
) # use 40% of points for longer edges
if n_gauss <= GAUSS_QUADRATURE_ORDER:
g, w = leggauss(n_gauss)
# map from [-1, 1] to [t_start, t_end]
s = (0.5 * (t_end - t_start) * g + 0.5 * (t_end + t_start)).astype(
GRADIENT_DTYPE_FLOAT, copy=False
)
weights = (w * 0.5 * (t_end - t_start)).astype(GRADIENT_DTYPE_FLOAT, copy=False)
else:
# for longer edges, use composite Gauss quadrature
pts_per_interval = GAUSS_QUADRATURE_ORDER
n_intervals = max(1, (n_gauss + pts_per_interval - 1) // pts_per_interval)
g_ref, w_ref = leggauss(pts_per_interval)
s_list = []
w_list = []
for i in range(n_intervals):
# sub-interval bounds in [0, 1]
a = i / n_intervals
b = (i + 1) / n_intervals
# map to [t_start, t_end]
a_mapped = t_start + a * (t_end - t_start)
b_mapped = t_start + b * (t_end - t_start)
# transform Gauss points to sub-interval [a_mapped, b_mapped]
s_interval = 0.5 * (b_mapped - a_mapped) * g_ref + 0.5 * (b_mapped + a_mapped)
w_interval = 0.5 * (b_mapped - a_mapped) * w_ref
s_list.extend(s_interval)
w_list.extend(w_interval)
s = np.array(s_list, dtype=GRADIENT_DTYPE_FLOAT)
weights = np.array(w_list, dtype=GRADIENT_DTYPE_FLOAT)
return s, weights
# pre-compute 2D bounding box for edge clipping optimization
sim_min_2d = np.delete(sim_min, self.axis)
sim_max_2d = np.delete(sim_max, self.axis)
for ei, (v0, v1) in enumerate(zip(vertices, next_v)):
edge_vec = v1 - v0
L = np.linalg.norm(edge_vec)
if np.isclose(L, 0.0):
continue
# check if edge is definitely inside 2D bounds
edge_min_2d = np.minimum(v0, v1)
edge_max_2d = np.maximum(v0, v1)
definitely_inside_2d = np.all(edge_min_2d >= sim_min_2d) and np.all(
edge_max_2d <= sim_max_2d
)
# process each z-slice with proper clipping
for zc in z_centers:
# fast path: if edge is definitely inside 2D bounds AND z is in bounds, skip clipping
if definitely_inside_2d and z0 <= zc <= z1:
t_start, t_end = 0.0, 1.0
else:
v0_3d = self.unpop_axis_vect(np.array([zc]), v0[None])[0]
v1_3d = self.unpop_axis_vect(np.array([zc]), v1[None])[0]
clip_bounds = compute_edge_clip_bounds(v0_3d, v1_3d, sim_min, sim_max)
if clip_bounds is None:
continue # edge is entirely outside bounds
t_start, t_end = clip_bounds
s_list, s_weights = get_adaptive_samples(L, dx, t_start, t_end)
pts2d = v0 + np.outer(s_list, edge_vec)
# patch areas include quadrature weights for accurate integration
patch_areas = L * s_weights * (dz_surf if not is_2d else 1.0)
# all points are within bounds by construction due to clipping
xyz = self.unpop_axis_vect(np.full(len(s_list), zc), pts2d)
n_patches = len(s_list)
# ensure we don't exceed pre-allocated size
if patch_idx + n_patches > estimated_patches:
# resize arrays if needed (rare case)
new_size = int((patch_idx + n_patches) * 1.5)
# at this stage, no other views into these arrays exist _yet_
# so `refcheck=False` should be safe
all_centers.resize((new_size, 3), refcheck=False)
all_areas.resize(new_size, refcheck=False)
all_normals.resize((new_size, 3), refcheck=False)
all_perps1.resize((new_size, 3), refcheck=False)
all_perps2.resize((new_size, 3), refcheck=False)
all_centers[patch_idx : patch_idx + n_patches] = xyz
all_areas[patch_idx : patch_idx + n_patches] = patch_areas
all_normals[patch_idx : patch_idx + n_patches] = basis["norm"][ei]
all_perps1[patch_idx : patch_idx + n_patches] = basis["perp1"][ei]
all_perps2[patch_idx : patch_idx + n_patches] = basis["perp2"][ei]
edge_norm_2d = normals_2d[ei]
all_edge_info.extend([(ei, s_val, edge_norm_2d) for s_val in s_list])
patch_idx += n_patches
if patch_idx > 0:
# trim arrays to actual size used
centers = all_centers[:patch_idx]
areas = all_areas[:patch_idx]
normals = all_normals[:patch_idx]
perps1 = all_perps1[:patch_idx]
perps2 = all_perps2[:patch_idx]
patch_vjps = derivative_info.evaluate_gradient_at_points(
centers, normals, perps1, perps2, interpolators
)
patch_vjps = (patch_vjps * areas).real
for vjp, (ei, s_val, edge_norm_2d) in zip(patch_vjps, all_edge_info):
# vertex weights: (1-s) for vertex i, s for vertex i+1
# quadrature weights are already included in areas via patch_areas
w0 = 1.0 - s_val
w1 = s_val
vjp_per_vertex[ei] += (w0 * vjp) * edge_norm_2d
vjp_per_vertex[(ei + 1) % len(vertices)] += (w1 * vjp) * edge_norm_2d
return vjp_per_vertex
[docs]
def edge_basis_vectors(
self,
edges: np.ndarray, # (N, 2)
) -> dict[str, np.ndarray]: # (N, 3)
"""Normalized basis vectors for ``normal`` direction, ``slab`` tangent direction and ``edge``."""
# ensure edges have consistent dtype
edges = edges.astype(GRADIENT_DTYPE_FLOAT, copy=False)
num_vertices, _ = edges.shape
zeros = np.zeros(num_vertices, dtype=GRADIENT_DTYPE_FLOAT)
ones = np.ones(num_vertices, dtype=GRADIENT_DTYPE_FLOAT)
# normalized vectors along edges
edges_norm_in_plane = self.normalize_vect(edges)
edges_norm_xyz = self.unpop_axis_vect(zeros, edges_norm_in_plane)
# normalized vectors from base of edges to tops of edges
cos_angle = np.cos(self.sidewall_angle, dtype=GRADIENT_DTYPE_FLOAT)
sin_angle = np.sin(self.sidewall_angle, dtype=GRADIENT_DTYPE_FLOAT)
slabs_axis_components = cos_angle * ones
# create axis_norm as array directly to avoid tuple->array conversion in np.cross
axis_norm = np.zeros(3, dtype=GRADIENT_DTYPE_FLOAT)
axis_norm[self.axis] = 1.0
slab_normal_xyz = -sin_angle * np.cross(edges_norm_xyz, axis_norm)
_, slab_normal_in_plane = self.pop_axis_vect(slab_normal_xyz)
slabs_norm_xyz = self.unpop_axis_vect(slabs_axis_components, slab_normal_in_plane)
# normalized vectors pointing in normal direction of edge
# cross yields inward normal when the extrusion axis is y, so negate once for axis==1
sign = (-1 if self.axis == 1 else 1) * (-1 if not self.is_ccw else 1)
normals_norm_xyz = sign * np.cross(edges_norm_xyz, slabs_norm_xyz)
return {
"norm": normals_norm_xyz,
"perp1": edges_norm_xyz,
"perp2": slabs_norm_xyz,
}
[docs]
def unpop_axis_vect(self, ax_coords: np.ndarray, plane_coords: np.ndarray) -> np.ndarray:
"""Combine coordinate along axis with coordinates on the plane tangent to the axis.
ax_coords.shape == [N]
plane_coords.shape == [N, 2]
return shape == [N, 3]
"""
n_pts = ax_coords.shape[0]
arr_xyz = np.zeros((n_pts, 3), dtype=ax_coords.dtype)
plane_axes = [i for i in range(3) if i != self.axis]
arr_xyz[:, self.axis] = ax_coords
arr_xyz[:, plane_axes] = plane_coords
return arr_xyz
[docs]
def pop_axis_vect(self, coord: np.ndarray) -> tuple[np.ndarray, tuple[np.ndarray, np.ndarray]]:
"""Combine coordinate along axis with coordinates on the plane tangent to the axis.
coord.shape == [N, 3]
return shape == ([N], [N, 2]
"""
arr_axis, arrs_plane = self.pop_axis(coord.T, axis=self.axis)
arrs_plane = np.array(arrs_plane).T
return arr_axis, arrs_plane
[docs]
@staticmethod
def normalize_vect(arr: np.ndarray) -> np.ndarray:
"""normalize an array shaped (N, d) along the `d` axis and return (N, 1)."""
norm = np.linalg.norm(arr, axis=-1, keepdims=True)
norm = np.where(norm == 0, 1, norm)
return arr / norm
[docs]
def translated(self, x: float, y: float, z: float) -> PolySlab:
"""Return a translated copy of this geometry.
Parameters
----------
x : float
Translation along x.
y : float
Translation along y.
z : float
Translation along z.
Returns
-------
:class:`PolySlab`
Translated copy of this ``PolySlab``.
"""
t_normal, t_plane = self.pop_axis((x, y, z), axis=self.axis)
translated_vertices = np.array(self.vertices) + np.array(t_plane)[None, :]
translated_slab_bounds = (self.slab_bounds[0] + t_normal, self.slab_bounds[1] + t_normal)
return self.updated_copy(vertices=translated_vertices, slab_bounds=translated_slab_bounds)
[docs]
def scaled(self, x: float = 1.0, y: float = 1.0, z: float = 1.0) -> PolySlab:
"""Return a scaled copy of this geometry.
Parameters
----------
x : float = 1.0
Scaling factor along x.
y : float = 1.0
Scaling factor along y.
z : float = 1.0
Scaling factor along z.
Returns
-------
:class:`Geometry`
Scaled copy of this geometry.
"""
scale_normal, scale_in_plane = self.pop_axis((x, y, z), axis=self.axis)
scaled_vertices = self.vertices * np.array(scale_in_plane)
scaled_slab_bounds = tuple(scale_normal * bound for bound in self.slab_bounds)
return self.updated_copy(vertices=scaled_vertices, slab_bounds=scaled_slab_bounds)
[docs]
def rotated(self, angle: float, axis: Union[Axis, Coordinate]) -> PolySlab:
"""Return a rotated copy of this geometry.
Parameters
----------
angle : float
Rotation angle (in radians).
axis : Union[int, Tuple[float, float, float]]
Axis of rotation: 0, 1, or 2 for x, y, and z, respectively, or a 3D vector.
Returns
-------
:class:`PolySlab`
Rotated copy of this ``PolySlab``.
"""
_, plane_axs = self.pop_axis([0, 1, 2], self.axis)
if (isinstance(axis, int) and axis == self.axis) or (
isinstance(axis, tuple) and all(axis[ax] == 0 for ax in plane_axs)
):
verts_3d = np.zeros((3, self.vertices.shape[0]))
verts_3d[plane_axs[0], :] = self.vertices[:, 0]
verts_3d[plane_axs[1], :] = self.vertices[:, 1]
rotation = RotationAroundAxis(angle=angle, axis=axis)
rotated_vertices = rotation.rotate_vector(verts_3d)
rotated_vertices = rotated_vertices[plane_axs, :].T
return self.updated_copy(vertices=rotated_vertices)
return super().rotated(angle=angle, axis=axis)
[docs]
def reflected(self, normal: Coordinate) -> PolySlab:
"""Return a reflected copy of this geometry.
Parameters
----------
normal : Tuple[float, float, float]
The 3D normal vector of the plane of reflection. The plane is assumed
to pass through the origin (0,0,0).
Returns
-------
-------
:class:`PolySlab`
Reflected copy of this ``PolySlab``.
"""
if math.isclose(normal[self.axis], 0):
_, plane_axs = self.pop_axis((0, 1, 2), self.axis)
verts_3d = np.zeros((3, self.vertices.shape[0]))
verts_3d[plane_axs[0], :] = self.vertices[:, 0]
verts_3d[plane_axs[1], :] = self.vertices[:, 1]
reflection = ReflectionFromPlane(normal=normal)
reflected_vertices = reflection.reflect_vector(verts_3d)
reflected_vertices = reflected_vertices[plane_axs, :].T
return self.updated_copy(vertices=reflected_vertices)
return super().reflected(normal=normal)
class ComplexPolySlabBase(PolySlab):
"""Interface for dividing a complex polyslab where self-intersecting polygon can
occur during extrusion. This class should not be used directly. Use instead
:class:`plugins.polyslab.ComplexPolySlab`."""
@pydantic.validator("vertices", always=True)
def no_self_intersecting_polygon_during_extrusion(cls, val, values):
"""Turn off the validation for this class."""
return val
@classmethod
def from_gds(
cls,
gds_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",
) -> list[PolySlab]:
"""Import :class:`.PolySlab` from a ``gdstk.Cell``.
Parameters
----------
gds_cell : gdstk.Cell
``gdstk.Cell`` containing 2D geometric data.
axis : int
Integer index into the polygon's slab axis. (0,1,2) -> (x,y,z).
slab_bounds: Tuple[float, float]
Minimum and maximum positions of the slab along ``axis``.
gds_layer : int
Layer index in the ``gds_cell``.
gds_dtype : int = None
Data-type index in the ``gds_cell``.
If ``None``, imports all data for this layer into the returned list.
gds_scale : float = 1.0
Length scale used in GDS file in units of MICROMETER.
For example, if gds file uses nanometers, set ``gds_scale=1e-3``.
Must be positive.
dilation : float = 0.0
Dilation of the polygon in the base by shifting each edge along its
normal outwards direction by a distance;
a negative value corresponds to erosion.
sidewall_angle : float = 0
Angle of the sidewall.
``sidewall_angle=0`` (default) specifies vertical wall,
while ``0<sidewall_angle<np.pi/2`` for the base to be larger than the top.
reference_plane : PlanePosition = "middle"
The position of the GDS layer. It can be at the ``bottom``, ``middle``,
or ``top`` of the PolySlab. E.g. if ``axis=1``, ``bottom`` refers to the
negative side of y-axis, and ``top`` refers to the positive side of y-axis.
Returns
-------
List[:class:`.PolySlab`]
List of :class:`.PolySlab` objects sharing ``axis`` and slab bound properties.
"""
# TODO: change for 2.0
# handle reference plane kwarg
all_vertices = PolySlab._load_gds_vertices(gds_cell, gds_layer, gds_dtype, gds_scale)
polyslabs = [
cls(
vertices=verts,
axis=axis,
slab_bounds=slab_bounds,
dilation=dilation,
sidewall_angle=sidewall_angle,
reference_plane=reference_plane,
)
for verts in all_vertices
]
return [sub_poly for sub_polys in polyslabs for sub_poly in sub_polys.sub_polyslabs]
@property
def geometry_group(self) -> base.GeometryGroup:
"""Divide a complex polyslab into a list of simple polyslabs, which
are assembled into a :class:`.GeometryGroup`.
Returns
-------
:class:`.GeometryGroup`
GeometryGroup for a list of simple polyslabs divided from the complex
polyslab.
"""
return base.GeometryGroup(geometries=self.sub_polyslabs)
@property
def sub_polyslabs(self) -> list[PolySlab]:
"""Divide a complex polyslab into a list of simple polyslabs.
Only neighboring vertex-vertex crossing events are treated in this
version.
Returns
-------
List[PolySlab]
A list of simple polyslabs.
"""
sub_polyslab_list = []
num_division_count = 0
# initialize sub-polyslab parameters
sub_polyslab_dict = self.dict(exclude={"type"}).copy()
if math.isclose(self.sidewall_angle, 0):
return [PolySlab.parse_obj(sub_polyslab_dict)]
sub_polyslab_dict.update({"dilation": 0}) # dilation accounted in setup
# initialize offset distance
offset_distance = 0
for dist_val in self._dilation_length:
dist_now = 0.0
vertices_now = self.reference_polygon
# constructing sub-polyslabs until reaching the base/top
while not math.isclose(dist_now, dist_val):
# bounds for sub-polyslabs assuming no self-intersection
slab_bounds = [
self._dilation_value_at_reference_to_coord(dist_now),
self._dilation_value_at_reference_to_coord(dist_val),
]
# 1) find out any vertices touching events between the current
# position to the base/top
max_dist = PolySlab._neighbor_vertices_crossing_detection(
vertices_now, dist_val - dist_now
)
# vertices touching events captured, update bounds for sub-polyslab
if max_dist is not None:
# max_dist doesn't have sign, so construct signed offset distance
offset_distance = max_dist * dist_val / abs(dist_val)
slab_bounds[1] = self._dilation_value_at_reference_to_coord(
dist_now + offset_distance
)
# 2) construct sub-polyslab
slab_bounds.sort() # for reference_plane=top/bottom, bounds need to be ordered
# direction of marching
reference_plane = "bottom" if dist_val / self._tanq < 0 else "top"
sub_polyslab_dict.update(
{
"slab_bounds": tuple(slab_bounds),
"vertices": vertices_now,
"reference_plane": reference_plane,
}
)
sub_polyslab_list.append(PolySlab.parse_obj(sub_polyslab_dict))
# Now Step 3
if max_dist is None:
break
dist_now += offset_distance
# new polygon vertices where collapsing vertices are removed but keep one
vertices_now = PolySlab._shift_vertices(vertices_now, offset_distance)[0]
vertices_now = PolySlab._remove_duplicate_vertices(vertices_now)
# all vertices collapse
if len(vertices_now) < 3:
break
# polygon collapse into 1D
if self.make_shapely_polygon(vertices_now).buffer(0).area < fp_eps:
break
vertices_now = PolySlab._orient(vertices_now)
num_division_count += 1
if num_division_count > _COMPLEX_POLYSLAB_DIVISIONS_WARN:
log.warning(
f"Too many self-intersecting events: the polyslab has been divided into "
f"{num_division_count} polyslabs; more than {_COMPLEX_POLYSLAB_DIVISIONS_WARN} may "
f"slow down the simulation."
)
return sub_polyslab_list
@property
def _dilation_length(self) -> list[float]:
"""dilation length from reference plane to the top/bottom of the polyslab."""
# for "bottom", only needs to compute the offset length to the top
dist = [self._extrusion_length_to_offset_distance(self.finite_length_axis)]
# reverse the dilation value if the reference plane is on the top
if self.reference_plane == "top":
dist = [-dist[0]]
# for middle, both directions
elif self.reference_plane == "middle":
dist = [dist[0] / 2, -dist[0] / 2]
return dist
def _dilation_value_at_reference_to_coord(self, dilation: float) -> float:
"""Compute the coordinate based on the dilation value to the reference plane."""
z_coord = -dilation / self._tanq + self.slab_bounds[0]
if self.reference_plane == "middle":
return z_coord + self.finite_length_axis / 2
if self.reference_plane == "top":
return z_coord + self.finite_length_axis
# bottom case
return z_coord
def intersections_tilted_plane(
self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4
) -> 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.
Returns
-------
List[shapely.geometry.base.BaseGeometry]
List of 2D shapes that intersect plane.
For more details refer to
`Shapely's Documentation <https://shapely.readthedocs.io/en/stable/project.html>`_.
"""
return [
shapely.unary_union(
[
base.Geometry.evaluate_inf_shape(shape)
for polyslab in self.sub_polyslabs
for shape in polyslab.intersections_tilted_plane(normal, origin, to_2D)
]
)
]
