From 3bcfa27ec00ef7452dce177f7916df619a80128e Mon Sep 17 00:00:00 2001 From: Weiliang Jin Date: Wed, 9 Mar 2022 16:53:21 -0800 Subject: [PATCH] Slanted PolySlab implementation More Validators for PolySlab to make sure the base polygon is valid Validator to make sure no vertices crossing during extrusion with sidewall Change "inside" function in geometry to include boundaries Implement itersection_side for slanted PolySlab --- .pylintrc | 2 +- test_local.sh | 1 + tests/test_sidewall.py | 264 ++++++++++++ tidy3d/components/geometry.py | 787 ++++++++++++++++++++++++++++++---- 4 files changed, 969 insertions(+), 85 deletions(-) create mode 100644 tests/test_sidewall.py diff --git a/.pylintrc b/.pylintrc index a98464ac1..04bb03c9e 100644 --- a/.pylintrc +++ b/.pylintrc @@ -1,7 +1,7 @@ [MASTER] extension-pkg-allow-list=pydantic ignore=material_library.py, plugins -good-names=ax, im, Lx, Ly, Lz, x0, y0, z0, x, y, z, u, v, w, f, t, y1, y2, x1, x2, xs, ys, zs, Ax, Nx, Ny, Nz, N, dl, rr, E, H, xx, yy, zz, dx, dy, Jx, Jy, Jz, Ex, Ey, Ez, Mx, My, Mz, Hx, Hy, Hz, dz, e, fp, dt, a, c, kx, ky, kz, k0, Jx_k, Jy_k, My_k, Mx_k, N_theta, N_phi, L_theta, L_phi, ux, uy, uz +good-names=ax, im, Lx, Ly, Lz, x0, y0, z0, x, y, z, u, v, w, f, t, y1, y2, x1, x2, x3, x4, y3, y4, xs, ys, zs, Ax, Nx, Ny, Nz, N, dl, rr, E, H, xx, yy, zz, dx, dy, Jx, Jy, Jz, Ex, Ey, Ez, Mx, My, Mz, Hx, Hy, Hz, dz, e, fp, dt, a, c, kx, ky, kz, k0, Jx_k, Jy_k, My_k, Mx_k, N_theta, N_phi, L_theta, L_phi, ux, uy, uz [BASIC] diff --git a/test_local.sh b/test_local.sh index e0749f0da..08435c29f 100755 --- a/test_local.sh +++ b/test_local.sh @@ -7,5 +7,6 @@ pytest -ra tests/test_grid.py pytest -ra tests/test_IO.py pytest -ra tests/test_material_library.py pytest -ra tests/test_plugins.py +pytest -ra tests/test_sidewall.py pytest --doctest-modules tidy3d/components --ignore=tidy3d/components/base.py \ No newline at end of file diff --git a/tests/test_sidewall.py b/tests/test_sidewall.py new file mode 100644 index 000000000..58712d066 --- /dev/null +++ b/tests/test_sidewall.py @@ -0,0 +1,264 @@ +from typing import Dict +import pytest +import numpy as np +import pydantic +from shapely.geometry import Polygon, Point + +import tidy3d as td +from tidy3d.log import ValidationError, SetupError + + +def setup_polyslab(vertices, dilation, angle, bounds): + s = td.PolySlab( + vertices=vertices, + slab_bounds=bounds, + axis=2, + dilation=dilation, + sidewall_angle=angle, + ) + return s + + +def minimal_edge_length(vertices): + """compute the minimal edge length in a polygon""" + vs = vertices.T.copy() + vsp = np.roll(vs.copy(), axis=-1, shift=-1) + edge_length = np.linalg.norm(vsp - vs, axis=0) + return np.min(edge_length) + + +def convert_valid_polygon(vertices): + """Given vertices that might have intersecting eduges, generate + vertices of a valid polygon + """ + poly = Polygon(vertices).buffer(1e-3) # make sure no intersecting edges + if type(poly) is not Polygon: + poly = poly[0] + + vertices_n = np.array(poly.exterior.coords[:]) + return vertices_n + + +# default values +bounds = (0, 0.5) +dilation = 0.0 +angle = 0 + + +def test_remove_duplicate(): + """ + Make sure redundant neighboring vertices are removed + """ + vertices = np.random.random((10, 2)) + vertices[0, :] = vertices[9, :] + vertices[1, :] = vertices[0, :] + vertices[5, :] = vertices[6, :] + + vertices = td.PolySlab._remove_duplicate_vertices(vertices) + assert vertices.shape[0] == 7 + + +def test_valid_polygon(): + """No intersecting edges""" + + # area = 0 + vertices = ((0, 0), (1, 0), (2, 0)) + with pytest.raises(SetupError) as e_info: + s = setup_polyslab(vertices, dilation, angle, bounds) + + # only two points + vertices = ((0, 0), (1, 0), (1, 0)) + with pytest.raises(SetupError) as e_info: + s = setup_polyslab(vertices, dilation, angle, bounds) + + # intersecting edges + vertices = ((0, 0), (1, 0), (1, 1), (0, 1), (0.5, -1)) + + with pytest.raises(SetupError) as e_info: + s = setup_polyslab(vertices, dilation, angle, bounds) + + +def test_crossing_square(): + """ + Vertices crossing detection for a simple square + """ + vertices = ((0, 0), (1, 0), (1, -1), (0, -1)) + dilation = 0.0 + angle = np.pi / 4 + s = setup_polyslab(vertices, dilation, angle, bounds) + + # dilation too significant + dilation = -1.1 + angle = 0 + with pytest.raises(SetupError) as e_info: + s = setup_polyslab(vertices, dilation, angle, bounds) + + # angle too large + dilation = 0 + angle = np.pi / 3 + with pytest.raises(SetupError) as e_info: + s = setup_polyslab(vertices, dilation, angle, bounds) + + # combines both + dilation = -0.1 + angle = np.pi / 4 + with pytest.raises(SetupError) as e_info: + s = setup_polyslab(vertices, dilation, angle, bounds) + + +def test_max_erosion_polygon(): + """ + Maximal erosion distance validation + """ + N = 10 # number of vertices + for i in range(50): + vertices = convert_valid_polygon(np.random.random((N, 2)) * 10) + + dilation = 0 + angle = 0 + bounds = (0, 0.5) + s = setup_polyslab(vertices, dilation, angle, bounds) + + # compute maximal allowed erosion distance + _, max_dist = s._crossing_detection(s.base_polygon, -100) + + # verify it is indeed maximal allowed + dilation = -max_dist + 1e-10 + # avoid vertex-edge crossing case + try: + s = setup_polyslab(vertices, dilation, angle, bounds) + except: + continue + assert np.isclose(minimal_edge_length(s.base_polygon), 0, atol=1e-5) + + # verify it is indeed maximal allowed + dilation = 0.0 + bounds = (0, max_dist - 1e-10) + angle = np.pi / 4 + + # avoid vertex-edge crossing case + try: + s = setup_polyslab(vertices, dilation, angle, bounds) + except: + continue + assert np.isclose(minimal_edge_length(s.top_polygon), 0, atol=1e-5) + + +def test_shift_height(): + """Make sure a list of height where the plane will intersect with the vertices + works properly + """ + N = 10 # number of vertices + Lx = 10.0 + for i in range(50): + vertices = convert_valid_polygon(np.random.random((N, 2)) * Lx) + dilation = 0 + angle = 0 + bounds = (0, 1) + s = setup_polyslab(vertices, dilation, angle, bounds) + # set up proper thickness + _, max_dist = s._crossing_detection(s.base_polygon, -100) + dilation = 0.0 + bounds = (0, max_dist * 0.99) + angle = np.pi / 4 + # avoid vertex-edge crossing case + try: + s = setup_polyslab(vertices, dilation, angle, bounds) + except: + continue + + for axis in (0, 1): + position = np.random.random(1)[0] * Lx - Lx / 2 + height = s._find_intersecting_height(position, axis) + for h in height: + bounds = (0, h) + s = setup_polyslab(vertices, dilation, angle, bounds) + diff = s.top_polygon[:, axis] - position + assert np.any(np.isclose(diff, 0)) == True + + +def test_intersection_with_inside(): + """Make sure intersection produces the same result as inside""" + + N = 10 # number of vertices + Lx = 10 # maximal length in x,y direction + for i in range(50): + vertices = convert_valid_polygon(np.random.random((N, 2)) * Lx) + vertices = np.array(vertices).astype("float32") + dilation = 0 + angle = 0 + bounds = (0, 1) + s = setup_polyslab(vertices, dilation, angle, bounds) + + # set up proper thickness + _, max_dist = s._crossing_detection(s.base_polygon, -100) + dilation = 0.0 + bounds = (0, np.float32(max_dist * 0.95)) + angle = np.pi / 4 + # avoid vertex-edge crossing case + try: + s = setup_polyslab(vertices, dilation, angle, bounds) + except: + continue + + ### side x + xp = np.random.random(1)[0] * 2 * Lx - Lx + yp = np.random.random(10) * 2 * Lx - Lx + zp = np.random.random(10) * (bounds[1] - bounds[0]) + bounds[0] + xp = np.float32(xp) + yp = np.float32(yp) + zp = np.float32(zp) + shape_intersect = s.intersections(x=xp) + + for i in range(len(yp)): + for j in range(len(zp)): + # inside + res_inside = s.inside(xp, yp[i], zp[j]) + # intersect + res_inter = False + for shape in shape_intersect: + if shape.covers(Point(yp[i], zp[j])): + res_inter = True + # if res_inter != res_inside: + # print(repr(vertices)) + # print(repr(s.base_polygon)) + # print(bounds) + # print(xp, yp[i], zp[j]) + assert res_inter == res_inside + + ### side y + xp = np.random.random(10) * 2 * Lx - Lx + yp = np.random.random(1)[0] * 2 * Lx - Lx + zp = np.random.random(10) * (bounds[1] - bounds[0]) + bounds[0] + xp = np.float32(xp) + yp = np.float32(yp) + zp = np.float32(zp) + shape_intersect = s.intersections(y=yp) + + for i in range(len(xp)): + for j in range(len(zp)): + # inside + res_inside = s.inside(xp[i], yp, zp[j]) + # intersect + res_inter = False + for shape in shape_intersect: + if shape.covers(Point(xp[i], zp[j])): + res_inter = True + assert res_inter == res_inside + + ### norm z + xp = np.random.random(10) * 2 * Lx - Lx + yp = np.random.random(10) * 2 * Lx - Lx + zp = np.random.random(1)[0] * (bounds[1] - bounds[0]) + bounds[0] + shape_intersect = s.intersections(z=zp) + + for i in range(len(xp)): + for j in range(len(yp)): + # inside + res_inside = s.inside(xp[i], yp[j], zp) + # intersect + res_inter = False + for shape in shape_intersect: + if shape.covers(Point(xp[i], yp[j])): + res_inter = True + assert res_inter == res_inside diff --git a/tidy3d/components/geometry.py b/tidy3d/components/geometry.py index 74daba4c5..96e0cd326 100644 --- a/tidy3d/components/geometry.py +++ b/tidy3d/components/geometry.py @@ -16,7 +16,11 @@ from .viz import add_ax_if_none, equal_aspect from .viz import PLOT_BUFFER, ARROW_LENGTH_FACTOR, ARROW_WIDTH_FACTOR, MAX_ARROW_WIDTH_FACTOR from ..log import Tidy3dKeyError, SetupError, ValidationError -from ..constants import MICROMETER, LARGE_NUMBER +from ..constants import MICROMETER, LARGE_NUMBER, RADIAN + +# for sampling polygon in slanted polyslab along z-direction for +# validating polygon to be non_intersecting. +_N_SAMPLE_POLYGON_INTERSECT = 100 class Geometry(Tidy3dBaseModel, ABC): @@ -502,13 +506,18 @@ def intersections(self, x: float = None, y: float = None, z: float = None): z0, _ = self.pop_axis(self.center, axis=self.axis) if (position < z0 - self.length / 2) or (position > z0 + self.length / 2): return [] - return self._intersections_normal() + return self._intersections_normal(position) return self._intersections_side(position, axis) @abstractmethod - def _intersections_normal(self) -> list: + def _intersections_normal(self, z: float) -> list: """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] @@ -709,7 +718,7 @@ def inside(self, x, y, z) -> bool: dist_x = np.abs(x - x0) dist_y = np.abs(y - y0) dist_z = np.abs(z - z0) - return (dist_x < Lx / 2) * (dist_y < Ly / 2) * (dist_z < Lz / 2) + return (dist_x <= Lx / 2) * (dist_y <= Ly / 2) * (dist_z <= Lz / 2) @property def bounds(self) -> Bound: @@ -917,9 +926,14 @@ class Cylinder(Circular, Planar): units=MICROMETER, ) - def _intersections_normal(self): + def _intersections_normal(self, z: float): """Find shapely geometries intersecting cylindrical geometry with axis normal to slab. + Parameters + ---------- + z : float + Position along the axis normal to slab + Returns ------- List[shapely.geometry.base.BaseGeometry] @@ -983,7 +997,7 @@ def inside(self, x, y, z) -> bool: dist_y = np.abs(y - y0) dist_z = np.abs(z - z0) inside_radius = (dist_x**2 + dist_y**2) <= (self.radius**2) - inside_height = dist_z < (self.length / 2) + inside_height = dist_z <= (self.length / 2) return inside_radius * inside_height @property @@ -1003,7 +1017,7 @@ def bounds(self): class PolySlab(Planar): - """Polygon with constant thickness (slab) along axis direction. + """Polygon extruded with optional sidewall angle along axis direction. Example ------- @@ -1018,10 +1032,30 @@ class PolySlab(Planar): units=MICROMETER, ) + dilation: float = pydantic.Field( + 0.0, + title="Dilation", + description="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.", + units=MICROMETER, + ) + + sidewall_angle: float = pydantic.Field( + 0.0, + title="Sidewall angle", + description="Angle of the sidewall. " + "``sidewall_angle=0`` (default) specifies vertical wall, " + "while ``0 List["PolySlab"]: """Import :class:`PolySlab` from a ``gdspy.Cell``. @@ -1097,6 +1244,14 @@ def from_gds( # pylint:disable=too-many-arguments 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 float: + """ + tan(sidewall_angle). _tanq*height gives the offset value + """ + return np.tan(self.sidewall_angle) + + @property + def base_polygon(self) -> tidynumpy: + """The polygon at the base after potential dilation operation. + + Returns + ------- + tidynumpy + The vertices of the polygon at the base. + """ + + return self._shift_vertices(self._proper_vertices(self.vertices), self.dilation)[0] + + @property + def top_polygon(self) -> tidynumpy: + """The polygon at the top after potential dilation and sidewall operation. + + Returns + ------- + tidynumpy + The vertices of the polygon at the top. + """ + + dist = -self.length * self._tanq + return self._shift_vertices(self.base_polygon, dist)[0] def inside(self, x, y, z) -> bool: # pylint:disable=too-many-locals """Returns True if point ``(x,y,z)`` inside volume of geometry. + For slanted polyslab and x/y/z to be np.ndarray, a loop over z-axis + is performed to find out the offsetted polygon at each z-coordinate. Parameters ---------- @@ -1142,36 +1340,65 @@ def inside(self, x, y, z) -> bool: # pylint:disable=too-many-locals bool Whether point ``(x,y,z)`` is inside geometry. """ + z0, _ = self.pop_axis(self.center, axis=self.axis) dist_z = np.abs(z - z0) - inside_height = dist_z < (self.length / 2) + inside_height = dist_z <= (self.length / 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) - face_polygon = Polygon(self.vertices) + z_local = z - z0 + self.length / 2 # distance to the base + dist = -z_local * self._tanq + + def contains_pointwise(face_polygon): + def fun_contain(xy_point): + point = Point(xy_point) + return face_polygon.covers(point) + + return fun_contain + if isinstance(x, np.ndarray): inside_polygon = np.zeros_like(inside_height) xs_slab = x[inside_height] ys_slab = y[inside_height] - def contains_pointwise(xy_point): - point = Point(xy_point) - return face_polygon.contains(point) - - contains_vectorized = np.vectorize(contains_pointwise, signature="(n)->()") - points_stacked = np.stack((xs_slab, ys_slab), axis=1) - inside_polygon_slab = contains_vectorized(points_stacked) - inside_polygon[inside_height] = inside_polygon_slab + # vertical sidewall + if np.isclose(self.sidewall_angle, 0): + face_polygon = Polygon(self.base_polygon) + fun_contain = contains_pointwise(face_polygon) + contains_vectorized = np.vectorize(fun_contain, signature="(n)->()") + points_stacked = np.stack((xs_slab, ys_slab), axis=1) + inside_polygon_slab = contains_vectorized(points_stacked) + inside_polygon[inside_height] = inside_polygon_slab + # slanted sidewall, offsetting vertices at each z + else: + for z_i in range(z.shape[2]): + if not inside_height[0, 0, z_i]: + continue + vertices_z = self._shift_vertices(self.base_polygon, dist[0, 0, z_i])[0] + face_polygon = Polygon(vertices_z) + fun_contain = contains_pointwise(face_polygon) + contains_vectorized = np.vectorize(fun_contain, signature="(n)->()") + points_stacked = np.stack((x[:, :, 0].flatten(), y[:, :, 0].flatten()), axis=1) + inside_polygon_slab = contains_vectorized(points_stacked) + inside_polygon[:, :, z_i] = inside_polygon_slab.reshape(x.shape[:2]) else: + vertices_z = self._shift_vertices(self.base_polygon, dist)[0] + face_polygon = Polygon(vertices_z) point = Point(x, y) - inside_polygon = face_polygon.contains(point) + inside_polygon = face_polygon.covers(point) return inside_height * inside_polygon - def _intersections_normal(self): - """Find shapely geometries intersecting planar geometry with axis normal to slab + 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 ------- @@ -1180,15 +1407,35 @@ def _intersections_normal(self): For more details refer to `Shapely's Documentaton `_. """ - return [Polygon(self.vertices)] + z0, _ = self.pop_axis(self.center, axis=self.axis) + z_local = z - z0 + self.length / 2 # distance to the base + dist = -z_local * self._tanq + vertices_z = self._shift_vertices(self.base_polygon, dist)[0] + return [Polygon(vertices_z)] def _intersections_side(self, position, axis) -> list: # pylint:disable=too-many-locals - """Find shapely geometries intersecting planar geometry with axis orthogonal to slab + """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`` + Position along ``axis``. axis : int Integer index into 'xyz' (0,1,2). @@ -1200,97 +1447,289 @@ def _intersections_side(self, position, axis) -> list: # pylint:disable=too-man `Shapely's Documentaton `_. """ + # find out all z_i where the plane will intersect the vertex z0, _ = self.pop_axis(self.center, axis=self.axis) - z_min, z_max = z0 - self.length / 2, z0 + self.length / 2 - iverts_b, iverts_f = self._find_intersecting_vertices(position, axis) - ints_y = self._find_intersecting_ys(iverts_b, iverts_f, position) - - # make polygon with intersections and z axis information + z_base = z0 - self.length / 2 + height_list = self._find_intersecting_height(position, axis) polys = [] - 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) - polys.append(box(minx=minx, miny=miny, maxx=maxx, maxy=maxy)) + + # looping through z_i to assemble the polygons + height = 0.0 + height_list = np.append(height_list, self.length) + for h_local in height_list: + dist = -height * self._tanq + vertices = self._shift_vertices(self.base_polygon, dist)[0] + z_min, z_max = z_base + height, z_base + height + h_local + # for vertical sidewall, no need for complications + if np.isclose(self.sidewall_angle, 0): + ints_y, ints_angle = self._find_intersecting_ys_angle_vertical( + vertices, position, axis + ) + else: + ints_y, ints_angle = self._find_intersecting_ys_angle_slant( + vertices, position, axis + ) + + # 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 np.isclose(self.sidewall_angle, 0): + polys.append(box(minx=minx, miny=miny, maxx=maxx, maxy=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_local * np.tan(angle_min) + dy_max = h_local * 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) + + if y_max - y_min <= dy_min + dy_max: + # intersect before reaching top of polygon + # make triangle + h_mid = (y_max - y_min) / (dy_min + dy_max) * h_local + z_mid = z_min + h_mid + y_mid = y_min + dy_min / h_local * h_mid + x3, y3 = self._order_by_axis(plane_val=y_mid, axis_val=z_mid, axis=axis) + vertices = ((x1, y1), (x2, y2), (x3, y3)) + polys.append(Polygon(vertices)) + else: + 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(Polygon(vertices)) + + height += h_local return polys - def _find_intersecting_vertices( - self, position: float, axis: int - ) -> Tuple[np.ndarray, np.ndarray]: - """Finds pairs of forward and backwards vertices where polygon intersects position at axis. - Assumes axis is handles so this function works on xy plane. + 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 + position along axis. axis : int Integer index into 'xyz' (0,1,2). Returns - np.ndarray, np.ndarray - Backward (xy) vertices and forward (xy) vertices. + np.ndarray + Height (relative to the base) where the plane will intersect with vertices. """ + if np.isclose(self.sidewall_angle, 0): + return np.array([]) + + vertices = self.base_polygon.copy() + + # shift rate + dist = 1.0 + shift_x, shift_y = PolySlab._shift_vertices(vertices, dist)[2] + shift_val = shift_x if axis == 0 else shift_y + shift_val[np.isclose(shift_val, 0)] = np.inf # for static vertices + + # distance to the plane in the direction of vertex shifting + distance = vertices[:, axis] - position + height = distance / self._tanq / shift_val + height = np.unique(height) + if self.sidewall_angle < 0: + height *= -1 + + height = height[height > 0] + height = height[height < self.length] + return height + + def _find_intersecting_ys_angle_vertical( # pylint:disable=too-many-locals + 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 unslanted polyslab). + Assumes axis is handles so this function works on xy plane. - vertices_b = np.array(self.vertices) + 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() # if the first coordinate refers to bounds, need to flip the vertices x,y if (axis == 2) or ((self.axis == 2) and (axis == 1)): - vertices_b = np.roll(vertices_b, shift=1, axis=1) + vertices_axis = np.roll(vertices_axis, shift=1, axis=1) # get the forward vertices - vertices_f = np.roll(vertices_b, shift=1, axis=0) + vertices_f = np.roll(vertices_axis, shift=-1, axis=0) # find which segments intersect - intersects_b = np.logical_and((vertices_f[:, 0] <= position), (vertices_b[:, 0] > position)) - intersects_f = np.logical_and((vertices_b[:, 0] <= position), (vertices_f[:, 0] > position)) + intersects_b = np.logical_and( + (vertices_f[:, 0] <= position), (vertices_axis[:, 0] > position) + ) + intersects_f = np.logical_and( + (vertices_axis[:, 0] <= position), (vertices_f[:, 0] > position) + ) intersects_segment = np.logical_or(intersects_b, intersects_f) - iverts_b = vertices_b[intersects_segment] + iverts_b = vertices_axis[intersects_segment] iverts_f = vertices_f[intersects_segment] - return iverts_b, iverts_f + # 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))) - @staticmethod - def _find_intersecting_ys( - iverts_b: np.ndarray, iverts_f: np.ndarray, position: float - ) -> List[float]: - """For each intersecting segment, find intersection point (in y) assuming straight line. + 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( # pylint:disable=too-many-locals, too-many-statements + 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 ---------- - iverts_b : np.ndarray - Backward (x,y) vertices. - iverts_f : np.ndarray - Forward (x,y) vertices. + vertices : np.ndarray + Shape (N, 2) defining the polygon vertices in the xy-plane. position : float - Position along coordinate x. + position along axis. + axis : int + Integer index into 'xyz' (0,1,2). Returns ------- - List[float] + Union[np.ndarray, np.ndarray] List of intersection points along y direction. + List of angles between plane and edges. """ - ints_y = [] - for (vertices_f, vertices_b) in zip(iverts_b, iverts_f): - x1, y1 = vertices_f - x2, y2 = vertices_b - slope = (y2 - y1) / (x2 - x1) - y = y1 + slope * (position - x1) - ints_y.append(y) - ints_y.sort() - return ints_y + vertices_axis = vertices.copy() + # if the first coordinate refers to bounds, need to flip the vertices x,y + if (axis == 2) or ((self.axis == 2) and (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) + 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) + 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(vertices, 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 + # case 1, neighbouring vertices on opposite side + if (x_b - position) * (x_f - position) < 0: + ints_y.append(y_on) + # vertices keep on the plane + if np.isclose(shift_local, 0): + ints_angle.append(np.pi / 2) + else: + x1, y1 = x_b, y_b + # decide which edge the plane will slide into + if shift_local * (x_f - position) < 0: + x1, y1, = ( + x_f, + y_f, + ) + slope = (y_on - y1) / (x_on - x1) + ints_angle.append(np.pi / 2 - np.arctan(np.abs(slope))) + # case 2, neightbouring vertices on the same side: + elif (x_b - position) * (x_f - position) > 0: + # this vertex is no longer relevant + if (x_f - position) * shift_local > 0: + continue + # the vertex will split into two + ints_y.append(y_on) + ints_y.append(y_on) + slope = (y_on - y_b) / (x_on - x_b) + ints_angle.append(np.pi / 2 - np.arctan(np.abs(slope))) + slope = (y_on - y_f) / (x_on - x_f) + ints_angle.append(np.pi / 2 - np.arctan(np.abs(slope))) + # case 3, one of neightbouring vertices on the plane too: + else: + # now only for anglepi/2 + continue + + 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 @property def _bounds(self): # get the min and max points in polygon plane - xpoints = tuple(c[0] for c in self.vertices) - ypoints = tuple(c[1] for c in self.vertices) - xmin, xmax = min(xpoints), max(xpoints) - ymin, ymax = min(ypoints), max(ypoints) + xpoints_base = tuple(c[0] for c in self.base_polygon) + ypoints_base = tuple(c[1] for c in self.base_polygon) + xpoints_top = tuple(c[0] for c in self.top_polygon) + ypoints_top = tuple(c[1] for c in self.top_polygon) + xmin = min(min(xpoints_base), min(xpoints_top)) + ymin = min(min(ypoints_base), min(ypoints_top)) + xmax = max(max(xpoints_base), max(xpoints_top)) + ymax = max(max(ypoints_base), max(ypoints_top)) z0, _ = self.pop_axis(self.center, axis=self.axis) zmin = z0 - self.length / 2 zmax = z0 + self.length / 2 @@ -1298,6 +1737,186 @@ def _bounds(self): coords_max = self.unpop_axis(zmax, (xmax, ymax), axis=self.axis) return (tuple(coords_min), tuple(coords_max)) + @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.copy(), axis=0, shift=-1) + term1 = vertices[:, 0] * vert_shift[:, 1] + term2 = vertices[:, 1] * vert_shift[:, 0] + + return np.sum(term1 - term2) * 0.5 + + @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. + """ + if PolySlab._area(vertices) > 0: + return vertices + + return 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.copy(), shift=-1, axis=0) + vertices_diff = np.linalg.norm(vertices - vertices_f, axis=1) + return vertices[~np.isclose(vertices_diff, 0)] + + @staticmethod + def _crossing_detection(vertices: np.ndarray, dist: float) -> Tuple[bool, float]: + """Detect if 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. + + Returns + ------- + Tuple[bool,float] + True if there are any crossings; + if True, return the maximal allowed dilation. + """ + + # edge length + vs_orig = vertices.T.copy() + vs_next = np.roll(vs_orig.copy(), axis=-1, shift=-1) + edge_length = np.linalg.norm(vs_next - vs_orig, axis=0) + + # edge length remaining + parallel_shift = PolySlab._shift_vertices(vertices, dist)[1] + parallel_shift_p = np.roll(parallel_shift.copy(), shift=-1) + edge_reduction = -(parallel_shift + parallel_shift_p) + length_remaining = edge_length - edge_reduction + + if np.any(length_remaining < 0): + index_oversized = length_remaining < 0 + max_dist = np.min(edge_length[index_oversized] / edge_reduction[index_oversized]) + max_dist *= np.abs(dist) + return True, max_dist + return False, None + + @staticmethod + def array_to_vertices(arr_vertices: np.ndarray) -> Vertices: + """Converts a numpy array of vertices to a list of tuples.""" + return [(x, y) for (x, y) in arr_vertices] # pylint:disable=unnecessary-comprehension + + @staticmethod + def _proper_vertices(vertices: Vertices) -> np.ndarray: + """convert vertices to np.array format, + removing duplicate neighbouring vertices, + and oriented in CCW direction. + + Returns + ------- + tidynumpy + 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 vertices_to_array(vertices_tuple: Vertices) -> np.ndarray: + """Converts a list of tuples (vertices) to a numpy array.""" + return np.array(vertices_tuple) + + @staticmethod + def _shift_vertices( # pylint:disable=too-many-locals + 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 + ---------- + vertices : 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 np.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 np.cross(u, v, axis=0) + + def normalize(v): + return v / np.linalg.norm(v, axis=0) + + vs_orig = vertices.T.copy() + vs_next = np.roll(vs_orig.copy(), axis=-1, shift=-1) + vs_previous = np.roll(vs_orig.copy(), 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), 0.0, cross(asm, rot90(asm - asp)) / (det + np.isclose(det, 0)) + ) + 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) + # geometries that can be used to define structures. GeometryFields = (Box, Sphere, Cylinder, PolySlab)