PR comments, small feature filter bug, fix for boundary conditions

Author: dbochkov-flexcompute

tidy3d/components/grid/corner_finder.py

@@ -1,6 +1,6 @@
 """Find corners of structures on a 2D plane."""
 
-from typing import List, Literal, Optional, Tuple
+from typing import Any, List, Literal, Optional, Tuple
 
 import numpy as np
 import pydantic.v1 as pd
@@ -11,7 +11,7 @@
 from ..geometry.utils import merging_geometries_on_plane
 from ..medium import PEC, LossyMetalMedium
 from ..structure import Structure
-from ..types import ArrayFloat1D, ArrayFloat2D, Axis
+from ..types import ArrayFloat1D, ArrayFloat2D, Axis, Shapely
 
 CORNER_ANGLE_THRESOLD = 0.1 * np.pi
 
@@ -81,7 +81,9 @@ def _merged_pec_on_plane(
         normal_axis: Axis,
         coord: float,
         structure_list: List[Structure],
-    ) -> Tuple[ArrayFloat2D, ArrayFloat1D]:
+        center: Tuple[float, float] = [0, 0, 0],
+        size: Tuple[float, float, float] = [inf, inf, inf],
+    ) -> List[Tuple[Any, Shapely]]:
         """On a 2D plane specified by axis = `normal_axis` and coordinate `coord`, merge geometries made of PEC.
 
         Parameters
@@ -92,19 +94,23 @@ def _merged_pec_on_plane(
             Position of plane along the normal axis.
         structure_list : List[Structure]
             List of structures present in simulation.
+        center : Tuple[float, float] = [0, 0, 0]
+            Center of the 2D plane (coordinate along ``axis`` is ignored)
+        size : Tuple[float, float, float] = [inf, inf, inf]
+            Size of the 2D plane (size along ``axis`` is ignored)
 
         Returns
         -------
-        ArrayFloat2D
-            Corner coordinates.
+        List[Tuple[Any, Shapely]]
+            List of shapes and their property value on the plane after merging.
         """
 
         # Construct plane
-        center = [0, 0, 0]
-        size = [inf, inf, inf]
-        center[normal_axis] = coord
-        size[normal_axis] = 0
-        plane = Box(center=center, size=size)
+        slice_center = list(center)
+        slice_size = list(size)
+        slice_center[normal_axis] = coord
+        slice_size[normal_axis] = 0
+        plane = Box(center=slice_center, size=slice_size)
 
         # prepare geometry and medium list
         geometry_list = [structure.geometry for structure in structure_list]
@@ -143,8 +149,8 @@ def _corners_and_convexity(
 
         Returns
         -------
-        ArrayFloat2D
-            Corner coordinates.
+        Tuple[ArrayFloat2D, ArrayFloat1D]
+            Corner coordinates and their convexity.
         """
 
         # merge geometries

tidy3d/components/grid/grid_spec.py

@@ -1052,8 +1052,9 @@ class LayerRefinementSpec(Box):
 
     dl_min_from_gap_width: bool = pd.Field(
         True,
-        title="Set ``dl_min`` from Gap Width",
-        description="Take into account autodetected minimal PEC gap width when determining ``dl_min``.",
+        title="Set ``dl_min`` from Estimated Gap Width",
+        description="Take into account autodetected minimal PEC gap width when determining ``dl_min``. "
+        "This only applies if ``dl_min`` in ``AutoGrid`` specification is not set.",
     )
 
     @pd.validator("axis", always=True)
@@ -1554,7 +1555,7 @@ def _override_structures_along_axis(
         return override_structures
 
     def _find_vertical_intersections(
-        self, grid_x_coords, grid_y_coords, poly_vertices
+        self, grid_x_coords, grid_y_coords, poly_vertices, boundary
     ) -> Tuple[List[Tuple[int, int]], List[float]]:
         """Detect intersection points of single polygon and vertical grid lines."""
 
@@ -1587,7 +1588,9 @@ def _find_vertical_intersections(
                 continue
 
             # sort vertices in ascending order to make treatmeant unifrom
+            reverse = False
             if ind_beg > ind_end:
+                reverse = True
                 ind_beg, ind_end, v_beg, v_end = ind_end, ind_beg, v_end, v_beg
 
             # x coordinates are simply x coordinates of intersected vertical grid lines
@@ -1623,6 +1626,12 @@ def _find_vertical_intersections(
                 grid_y_coords[cell_js + 1] - grid_y_coords[cell_js]
             )
 
+            # preserve uniform ordering along perimeter of the polygon
+            if reverse:
+                cell_is = cell_is[::-1]
+                cell_js = cell_js[::-1]
+                dy = dy[::-1]
+
             # record info
             cells_ij.append(np.transpose([cell_is, cell_js]))
             cells_dy.append(dy)
@@ -1665,39 +1674,55 @@ def _find_vertical_intersections(
             points_to_duplicate_near_zero = cells_ij[close_to_zero]
             points_to_duplicate_near_one = cells_ij[close_to_one]
 
+            # if we go beyond simulation domain boundary, either ignore
+            # or wrap periodically depending on boundary conditions
+            cells_ij_zero_side = points_to_duplicate_near_zero - np.array([0, 1])
+            cells_zero_side_out = cells_ij_zero_side[:, 1] == -1
+            if boundary[0] == "periodic":
+                cells_ij_zero_side[cells_zero_side_out, 1] = len(grid_y_coords) - 2
+            else:
+                cells_ij_zero_side = cells_ij_zero_side[cells_zero_side_out == 0]
+
+            cells_ij_one_side = points_to_duplicate_near_one + np.array([0, 1])
+            cells_one_side_out = cells_ij_one_side[:, 1] == len(grid_y_coords) - 1
+            if boundary[1] == "periodic":
+                cells_ij_one_side[cells_one_side_out, 1] = 0
+            else:
+                cells_ij_one_side = cells_ij_one_side[cells_one_side_out == 0]
+
             cells_ij = np.concatenate(
                 [
                     cells_ij,
-                    points_to_duplicate_near_zero - np.array([0, 1]),
-                    points_to_duplicate_near_one + np.array([0, 1]),
+                    cells_ij_zero_side,
+                    cells_ij_one_side,
                 ]
             )
             cells_dy = np.concatenate(
                 [
                     cells_dy,
-                    np.ones(len(points_to_duplicate_near_zero)),
-                    np.zeros(len(points_to_duplicate_near_one)),
+                    np.ones(len(cells_ij_zero_side)),
+                    np.zeros(len(cells_ij_one_side)),
                 ]
             )
 
         return cells_ij, cells_dy
 
     def _process_poly(
-        self, grid_x_coords, grid_y_coords, poly_vertices
+        self, grid_x_coords, grid_y_coords, poly_vertices, boundaries
     ) -> Tuple[List[Tuple[int, int]], List[float], List[Tuple[int, int]], List[float]]:
         """Detect intersection points of single polygon and grid lines."""
 
         # find cells that contain intersections of vertical grid lines
         # and relative locations of those intersections (along y axis)
         v_cells_ij, v_cells_dy = self._find_vertical_intersections(
-            grid_x_coords, grid_y_coords, poly_vertices
+            grid_x_coords, grid_y_coords, poly_vertices, boundaries[1]
         )
 
         # find cells that contain intersections of horizontal grid lines
         # and relative locations of those intersections (along x axis)
         # reuse the same command but flip dimensions
         h_cells_ij, h_cells_dx = self._find_vertical_intersections(
-            grid_y_coords, grid_x_coords, np.flip(poly_vertices, axis=1)
+            grid_y_coords, grid_x_coords, np.flip(poly_vertices, axis=1), boundaries[0]
         )
         if len(h_cells_ij) > 0:
             # flip dimensions back
@@ -1706,7 +1731,7 @@ def _process_poly(
         return v_cells_ij, v_cells_dy, h_cells_ij, h_cells_dx
 
     def _process_slice(
-        self, x, y, merged_geos
+        self, x, y, merged_geos, boundaries
     ) -> Tuple[List[Tuple[int, int]], List[float], List[Tuple[int, int]], List[float]]:
         """Detect intersection points of geometries boundaries and grid lines."""
 
@@ -1720,6 +1745,27 @@ def _process_slice(
         # relative locations of those intersections (along x axis)
         h_cells_dx = []
 
+        # for PEC and PMC boundary - treat them as PEC structure
+        # so that gaps are resolved near boundaries if any
+        nx = len(x)
+        ny = len(y)
+
+        if boundaries[0][0] == "pec/pmc":
+            h_cells_ij.append(np.transpose([np.zeros(ny), np.arange(ny)]).astype(int))
+            h_cells_dx.append(np.zeros(ny))
+
+        if boundaries[0][1] == "pec/pmc":
+            h_cells_ij.append(np.transpose([(nx - 2) * np.ones(ny), np.arange(ny)]).astype(int))
+            h_cells_dx.append(np.ones(ny))
+
+        if boundaries[1][0] == "pec/pmc":
+            v_cells_ij.append(np.transpose([np.arange(nx), np.zeros(nx)]).astype(int))
+            v_cells_dy.append(np.zeros(nx, dtype=int))
+
+        if boundaries[1][1] == "pec/pmc":
+            v_cells_ij.append(np.transpose([np.arange(nx), (ny - 2) * np.ones(nx)]).astype(int))
+            v_cells_dy.append(np.ones(nx))
+
         # loop over all shapes
         for mat, shapes in merged_geos:
             if not mat.is_pec:
@@ -1736,7 +1782,7 @@ def _process_slice(
                 # 1. relative locations of those intersections along y axis
                 # 2. cells that contain intersections of horizontal grid lines
                 # 3. relative locations of those intersections along x axis
-                data = self._process_poly(x, y, np.array(poly.exterior.coords)[:-1])
+                data = self._process_poly(x, y, np.array(poly.exterior.coords)[:-1], boundaries)
 
                 if len(data[0]) > 0:
                     v_cells_ij.append(data[0])
@@ -1748,7 +1794,7 @@ def _process_slice(
 
                 # in case the polygon has holes
                 for poly_inner in poly.interiors:
-                    data = self._process_poly(x, y, np.array(poly_inner.coords)[:-1])
+                    data = self._process_poly(x, y, np.array(poly_inner.coords)[:-1], boundaries)
                     if len(data[0]) > 0:
                         v_cells_ij.append(data[0])
                         v_cells_dy.append(data[1])
@@ -1845,27 +1891,82 @@ def _generate_horizontal_snapping_lines(
 
         return snapping_lines_y, min_gap_width
 
-    def _resolve_gaps(self, structures: List, grid: Grid) -> Tuple[List[CoordinateOptional], float]:
+    def _resolve_gaps(
+        self, structures: List[Structure], grid: Grid, boundaries: Tuple, center, size
+    ) -> Tuple[List[CoordinateOptional], float]:
         """Detect underresolved gaps and place snapping lines in them. Also return the detected minimal gap width."""
 
-        # get merged pec structures on plane
-        # note that we expect this function to also convert all LossyMetal's into PEC
-        plane_slice = CornerFinderSpec._merged_pec_on_plane(
-            coord=self.center_axis, normal_axis=self.axis, structure_list=structures
-        )
-
         # get x and y coordinates of grid lines
         _, tan_dims = Box.pop_axis([0, 1, 2], self.axis)
         x = grid.boundaries.to_list[tan_dims[0]]
         y = grid.boundaries.to_list[tan_dims[1]]
 
+        _, boundaries_tan = Box.pop_axis(boundaries, self.axis)
+
+        # restrict to the size of layer spec
+        rmin, rmax = self.bounds
+        _, rmin = Box.pop_axis(rmin, self.axis)
+        _, rmax = Box.pop_axis(rmax, self.axis)
+
+        new_coords = []
+        new_boundaries = []
+        for coord, cmin, cmax, bdry in zip([x, y], rmin, rmax, boundaries_tan):
+            if cmax <= coord[0] or cmin >= coord[-1]:
+                return [], inf
+            else:
+                if cmin < coord[0]:
+                    ind_min = 0
+                else:
+                    ind_min = max(0, np.argmax(coord >= cmin) - 1)
+
+                if cmax > coord[-1]:
+                    ind_max = len(coord) - 1
+                else:
+                    ind_max = np.argmax(coord >= cmax)
+
+                if ind_min >= ind_max:
+                    return [], inf
+
+                new_coords.append(coord[ind_min : (ind_max + 1)])
+                # ignore boundary conditions if we are not touching them
+                new_boundaries.append(
+                    [
+                        None if ind_min > 0 else bdry[0],
+                        None if ind_max < len(coord) - 1 else bdry[1],
+                    ]
+                )
+
+        x, y = new_coords
+
+        # restrict size of the plane where pec polygons are found in case of periodic boundary conditions
+        # this is to make sure gaps across periodic boundary conditions are resolved
+        # (if there is a PEC structure going into periodic boundary, now it will generate a grid line
+        # intersection next to that boundary and it will be propagated to the other side)
+        restricted_size_tan = [
+            s - GAP_MESHING_TOL / 2 if b[0] == "periodic" else inf
+            for b, s in zip(new_boundaries, size)
+        ]
+        restricted_size = Box.unpop_axis(size[self.axis], restricted_size_tan, self.axis)
+
+        # get merged pec structures on plane
+        # note that we expect this function to also convert all LossyMetal's into PEC
+        plane_slice = CornerFinderSpec._merged_pec_on_plane(
+            coord=self.center_axis,
+            normal_axis=self.axis,
+            structure_list=structures,
+            center=center,
+            size=restricted_size,
+        )
+
         # find intersections of pec polygons with grid lines
         # specifically:
         # 0. cells that contain intersections of vertical grid lines
         # 1. relative locations of those intersections along y axis
         # 2. cells that contain intersections of horizontal grid lines
         # 3. relative locations of those intersections along x axis
-        v_cells_ij, v_cells_dy, h_cells_ij, h_cells_dx = self._process_slice(x, y, plane_slice)
+        v_cells_ij, v_cells_dy, h_cells_ij, h_cells_dx = self._process_slice(
+            x, y, plane_slice, new_boundaries
+        )
 
         # generate horizontal snapping lines
         snapping_lines_y, min_gap_width_along_y = self._generate_horizontal_snapping_lines(
@@ -2289,6 +2390,11 @@ def make_grid(
         lumped_elements: List[LumpedElementType] = (),
         internal_override_structures: List[MeshOverrideStructure] = None,
         internal_snapping_points: List[CoordinateOptional] = None,
+        boundary_types: Tuple[Tuple[str, str], Tuple[str, str], Tuple[str, str]] = [
+            [None, None],
+            [None, None],
+            [None, None],
+        ],
     ) -> Grid:
         """Make the entire simulation grid based on some simulation parameters.
 
@@ -2313,6 +2419,9 @@ def make_grid(
             If `None`, recomputes internal override structures.
         internal_snapping_points : List[CoordinateOptional]
             If `None`, recomputes internal snapping points.
+        boundary_types : Tuple[Tuple[str, str], Tuple[str, str], Tuple[str, str]] = [[None, None], [None, None], [None, None]]
+            Type of boundary conditions along each dimension: "pec/pmc", "periodic", or
+            None for any other. This is relevant only for gap meshing.
 
         Returns
         -------
@@ -2343,6 +2452,11 @@ def _make_grid_and_snapping_lines(
         lumped_elements: List[LumpedElementType] = (),
         internal_override_structures: List[MeshOverrideStructure] = None,
         internal_snapping_points: List[CoordinateOptional] = None,
+        boundary_types: Tuple[Tuple[str, str], Tuple[str, str], Tuple[str, str]] = [
+            [None, None],
+            [None, None],
+            [None, None],
+        ],
     ) -> Tuple[Grid, List[CoordinateOptional]]:
         """Make the entire simulation grid based on some simulation parameters.
         Also return snappiung point resulted from iterative gap meshing.
@@ -2368,11 +2482,14 @@ def _make_grid_and_snapping_lines(
             If `None`, recomputes internal override structures.
         internal_snapping_points : List[CoordinateOptional]
             If `None`, recomputes internal snapping points.
+        boundary_types : Tuple[Tuple[str, str], Tuple[str, str], Tuple[str, str]] = [[None, None], [None, None], [None, None]]
+            Type of boundary conditions along each dimension: "pec/pmc", "periodic", or
+            None for any other. This is relevant only for gap meshing.
 
         Returns
         -------
-        Grid:
-            Entire simulation grid.
+        Tuple[Grid, List[CoordinateOptional]]:
+            Entire simulation grid and snapping points generated during iterative gap meshing.
         """
 
         old_grid = self._make_grid_one_iteration(
@@ -2386,6 +2503,8 @@ def _make_grid_and_snapping_lines(
             internal_snapping_points=internal_snapping_points,
         )
 
+        sim_geometry = structures[0].geometry
+
         snapping_lines = []
         if len(self.layer_refinement_specs) > 0:
             num_iters = max(
@@ -2398,7 +2517,11 @@ def _make_grid_and_snapping_lines(
                 for layer_spec in self.layer_refinement_specs:
                     if layer_spec.gap_meshing_iters > ind:
                         one_layer_snapping_lines, gap_width = layer_spec._resolve_gaps(
-                            structures, old_grid
+                            structures,
+                            old_grid,
+                            boundary_types,
+                            center=sim_geometry.center,
+                            size=sim_geometry.size,
                         )
                         new_snapping_lines = new_snapping_lines + one_layer_snapping_lines
                         if layer_spec.dl_min_from_gap_width: