more comments

Author: dbochkov-flexcompute

tidy3d/components/grid/grid_spec.py

@@ -1539,36 +1539,66 @@ def _find_vertical_intersections(self, x, y, vertices):
 
         cells_ij = []
         cells_dy = []
-        h_inds = np.argmax(x[:, None] > vertices[None, :, 0], axis=0)
-        h_inds[vertices[:, 0] > x[-1]] = len(x)
+
+        # for each polygon vertex find the index of the first grid line on the right
+        grid_lines_on_right = np.argmax(x[:, None] > vertices[None, :, 0], axis=0)
+        grid_lines_on_right[vertices[:, 0] > x[-1]] = len(x)
+        # once we know these indices then we can find grid lines intersected by the i-th
+        # segment of the polygon as
+        # [grid_lines_on_right[i], grid_lines_on_right[i+1]) for grid_lines_on_right[i] > grid_lines_on_right[i+1]
+        # or
+        # [grid_lines_on_right[i+1], grid_lines_on_right[i]) for grid_lines_on_right[i] < grid_lines_on_right[i+1]
+
+        # loop over segments of the polygon and determine in which cells and where exactly they cross grid lines
+        # v_beg and v_end are the starting and ending points of the segment
+        # ind_beg and ind_end are starting and ending indices of vertical grid lines that the segment intersects
+        # as described above
         for ind_beg, ind_end, v_beg, v_end in zip(
-            h_inds, np.roll(h_inds, -1), vertices, np.roll(vertices, axis=0, shift=-1)
+            grid_lines_on_right,
+            np.roll(grid_lines_on_right, -1),
+            vertices,
+            np.roll(vertices, axis=0, shift=-1),
         ):
-            # sort vertices in ascending order
+            # no intersections
+            if ind_end == ind_beg:
+                continue
+
+            # sort vertices in ascending order to make treatmeant unifrom
             if ind_beg > ind_end:
                 ind_beg, ind_end, v_beg, v_end = ind_end, ind_beg, v_end, v_beg
 
-            if ind_end > ind_beg:
-                # x coordinates are simply grid x points that are covered by each segment
-                h_x = x[ind_beg:ind_end]
+            # x coordinates are simply x coordinates of intersected vertical grid lines
+            intersections_x = x[ind_beg:ind_end]
+
+            # y coordinates can be found from line equation
+            intersections_y = v_beg[1] + (v_end[1] - v_beg[1]) / (v_end[0] - v_beg[0]) * (
+                intersections_x - v_beg[0]
+            )
 
-                # y coordinates can be found from line equation
-                h_y = v_beg[1] + (v_end[1] - v_beg[1]) / (v_end[0] - v_beg[0]) * (h_x - v_beg[0])
+            # however, some of the vertical lines might be crossed
+            # outside of computational domain
+            # so we need to see which ones are actually inside along y axis
+            inds_inside_grid = np.logical_and(intersections_y >= y[0], intersections_y <= y[-1])
 
-                # however, we need to see which ones are actually inside along y axis
-                inds_inside_grid = np.logical_and(h_y >= y[0], h_y <= y[-1])
+            intersections_y = intersections_y[inds_inside_grid]
 
-                h_y = h_y[inds_inside_grid]
+            # find i and j indices of cells which contain these intersections
 
-                # find cell j coordinates which contain these intersections
-                cell_is = np.arange(ind_beg, ind_end)[inds_inside_grid]
-                cell_js = np.searchsorted(y, h_y) - 1
+            # i indices are simply indices of crossed vertical grid lines
+            cell_is = np.arange(ind_beg, ind_end)[inds_inside_grid]
 
-                # find local dy
-                dy = h_y - y[cell_js]
+            # j indices can be computed by finding insertion indices
+            # of y coordinates of intersection points into array of y coordinates
+            # of the grid lines that preserve sorting
+            cell_js = np.searchsorted(y, intersections_y) - 1
 
-                cells_ij.append(np.transpose([cell_is, cell_js]))
-                cells_dy.append(dy)
+            # find local dy, that is, the distance between the intersection point
+            # and the bottom edge of the cell
+            dy = intersections_y - y[cell_js]
+
+            # record info
+            cells_ij.append(np.transpose([cell_is, cell_js]))
+            cells_dy.append(dy)
 
         if len(cells_ij) > 0:
             cells_ij = np.concatenate(cells_ij)
@@ -1580,26 +1610,47 @@ def _find_vertical_intersections(self, x, y, vertices):
 
     def _process_poly(self, x, y, vertices):
         """Detect intersection points of single polygon and grid lines."""
+
+        # find cells that contain intersections of vertical grid lines
+        # and locations of those intersections (along y axis)
         v_cells_ij, v_cells_dy = self._find_vertical_intersections(x, y, vertices)
+
+        # find cells that contain intersections of horizontal grid lines
+        # and locations of those intersections (along x axis)
         # reuse the same command but flip dimensions
-        h_cells_ij, h_cells_dx = self._find_vertical_intersections(y, x, vertices[:, [1, 0]])
+        h_cells_ij, h_cells_dx = self._find_vertical_intersections(y, x, np.flip(vertices, axis=1))
         if len(h_cells_ij) > 0:
             h_cells_ij = np.roll(h_cells_ij, axis=1, shift=1)
+
         return v_cells_ij, v_cells_dy, h_cells_ij, h_cells_dx
 
     def _process_slice(self, x, y, merged_geos):
         """Detect intersection points of geometries boundaries and grid lines."""
+
+        # cells that contain intersections of vertical grid lines
         v_cells_ij = []
+        # locations of those intersections (along y axis)
         v_cells_dy = []
+
+        # cells that contain intersections of horizontal grid lines
         h_cells_ij = []
+        # locations of those intersections (along x axis)
         h_cells_dx = []
 
+        # loop over all shapes
         for mat, shapes in merged_geos:
             if not mat.is_pec:
                 continue
             polygon_list = ClipOperation.to_polygon_list(shapes)
             for poly in polygon_list:
                 poly = poly.normalize().buffer(0)
+
+                # find intersections of a polygon with grid lines
+                # specifically:
+                # 0. cells that contain intersections of vertical grid lines
+                # 1. locations of those intersections along y axis
+                # 2. cells that contain intersections of horizontal grid lines
+                # 3. locations of those intersections along x axis
                 data = self._process_poly(x, y, np.array(poly.exterior.coords)[:-1])
 
                 if len(data[0]) > 0:
@@ -1643,48 +1694,88 @@ def _resolve_gaps(self, structures: List, grid: Grid):
             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]]
 
-        # find in which cells pec boundaries cross grid lines and locations of intersections
+        # find intersections of pec polygons with grid lines
+        # specifically:
+        # 0. cells that contain intersections of vertical grid lines
+        # 1. locations of those intersections along y axis
+        # 2. cells that contain intersections of horizontal grid lines
+        # 3. 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)
 
         detected_gap_width = inf
 
-        # count crossings of vertical grid lines
+        # generate horizontal snapping lines
         snapping_lines_y = []
         if len(v_cells_ij) > 0:
-            # add number of intersections
+            # count intersections of vertical grid lines in each cell
             ind_unique, counts = np.unique(
                 v_cells_ij[:, 0] * len(y) + v_cells_ij[:, 1], return_counts=True
             )
+            # when we count intersections we use linearized 2d index because we really
+            # need to count intersections in each cell separately
+
+            # but when we need to decide about refinement, due to cartesian nature of grid
+            # we will need to consider all cells with a given j index at a time
+
+            # so, let's compute j index for each cell in the inique list
             ind_unique_j = ind_unique % len(y)
+
+            # loop through all j rows that contain intersections
             for ind_j in np.unique(ind_unique_j):
+                # we need to refine between two grid lines corresponding to index j
+                # if at least one cell with given j contains > 1 intersections
                 max_count = np.max(counts[ind_unique_j == ind_j])
-                # for those where we have more than one crossing, place a snapping line in the middle
                 if max_count > 1:
+                    # This could use some improvement in future:
+                    # currently we just calculate the position of the lowest intersection
+                    # and the highest intersection in row j
+                    # and just place a snapping line in between.
                     min_dy = np.min(v_cells_dy[v_cells_ij[:, 1] == ind_j])
                     max_dy = np.max(v_cells_dy[v_cells_ij[:, 1] == ind_j])
                     snapping_lines_y.append(y[ind_j] + 0.5 * (min_dy + max_dy))
                     detected_gap_width = min(detected_gap_width, max_dy - min_dy)
+                    # This is expected to work fine when we have a simple gap/strip passing
+                    # through a number of cells. But not optimal if we have slanted/curved
+                    # boundaries and/or more than 2 intersections.
 
-        # count crossings of horizontal grid lines
+        # generate vertical snapping lines
         snapping_lines_x = []
         if len(h_cells_ij) > 0:
-            # add number of intersections
+            # count intersections of horizontal grid lines in each cell
             ind_unique, counts = np.unique(
                 h_cells_ij[:, 0] * len(y) + h_cells_ij[:, 1], return_counts=True
             )
+            # when we count intersections we use linearized 2d index because we really
+            # need to count intersections in each cell separately
+
+            # but when we need to decide about refinement, due to cartesian nature of grid
+            # we will need to consider all cells with a given i index at a time
+
+            # so, let's compute i index for each cell in the inique list
             ind_unique_i = ind_unique // len(y)
+
+            # loop through all i columns that contain intersections
             for ind_i in np.unique(ind_unique_i):
+                # we need to refine between two grid lines corresponding to index i
+                # if at least one cell with given i contains > 1 intersections
                 max_count = np.max(counts[ind_unique_i == ind_i])
-                # for those where we have more than one crossing, place a snapping line in the middle
                 if max_count > 1:
+                    # This could use some improvement in future:
+                    # currently we just calculate the position of the leftmost intersection
+                    # and the rightmost intersection in column i
+                    # and just place a snapping line in between.
                     min_dx = np.min(h_cells_dx[h_cells_ij[:, 0] == ind_i])
                     max_dx = np.max(h_cells_dx[h_cells_ij[:, 0] == ind_i])
                     snapping_lines_x.append(x[ind_i] + 0.5 * (min_dx + max_dx))
                     detected_gap_width = min(detected_gap_width, max_dx - min_dx)
+                    # This is expected to work fine when we have a simple gap/strip passing
+                    # through a number of cells. But not optimal if we have slanted/curved
+                    # boundaries and/or more than 2 intersections.
 
         return [Box.unpop_axis(X, (None, None), axis=tan_dims[0]) for X in snapping_lines_x] + [
             Box.unpop_axis(Y, (None, None), axis=tan_dims[1]) for Y in snapping_lines_y