formatting

Author: dbochkov-flexcompute

tests/test_components/test_layerrefinement.py

@@ -285,7 +285,8 @@ def count_grids_within_layer(sim_t):
     )
 
 
-def test_corner_refinement_outside_domain():
+@pytest.mark.parametrize("gap_meshing_iters", [0, 1])
+def test_corner_refinement_outside_domain(gap_meshing_iters):
     """Test the behavior of corner refinement if corners are outside the simulation domain."""
 
     # CPW waveguides that goes through the simulation domain along x-axis, so that the corners
@@ -313,6 +314,7 @@ def test_corner_refinement_outside_domain():
         axis=2,
         corner_refinement=td.GridRefinement(refinement_factor=5),
         refinement_inside_sim_only=True,
+        gap_meshing_iters=gap_meshing_iters,
     )
 
     sim = td.Simulation(
@@ -333,11 +335,12 @@ def count_grids_within_gap(sim_t):
         return len(y)
 
     # just 2 grids sampling the gap
-    assert count_grids_within_gap(sim) == 2
+    assert count_grids_within_gap(sim) == gap_meshing_iters + 2
 
     # 2) refined if corners outside simulation domain is accounted for.
     layer = layer.updated_copy(refinement_inside_sim_only=False)
     sim = sim.updated_copy(grid_spec=td.GridSpec.auto(wavelength=1, layer_refinement_specs=[layer]))
+
     assert count_grids_within_gap(sim) > 2
 
 

tidy3d/components/grid/grid_spec.py

@@ -1553,7 +1553,9 @@ def _override_structures_along_axis(
                 override_structures += refinement_structures
         return override_structures
 
-    def _find_vertical_intersections(self, grid_x_coords, grid_y_coords, poly_vertices) -> Tuple[List[Tuple[int, int]], List[float]]:
+    def _find_vertical_intersections(
+        self, grid_x_coords, grid_y_coords, poly_vertices
+    ) -> Tuple[List[Tuple[int, int]], List[float]]:
         """Detect intersection points of single polygon and vertical grid lines."""
 
         # indices of cells that contain intersection with grid lines (left edge of a cell)
@@ -1599,7 +1601,9 @@ def _find_vertical_intersections(self, grid_x_coords, grid_y_coords, poly_vertic
             # 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])
+            inds_inside_grid = np.logical_and(
+                intersections_y >= grid_y_coords[0], intersections_y <= grid_y_coords[-1]
+            )
 
             intersections_y = intersections_y[inds_inside_grid]
 
@@ -1615,60 +1619,94 @@ def _find_vertical_intersections(self, grid_x_coords, grid_y_coords, poly_vertic
 
             # find local dy, that is, the distance between the intersection point
             # and the bottom edge of the cell
-            dy = (intersections_y - grid_y_coords[cell_js]) / (grid_y_coords[cell_js + 1] - grid_y_coords[cell_js])
+            dy = (intersections_y - grid_y_coords[cell_js]) / (
+                grid_y_coords[cell_js + 1] - grid_y_coords[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)
-
-        if len(cells_dy) > 0:
             cells_dy = np.concatenate(cells_dy)
 
+            # Filter from re-entering subcell features. That is, we discard any consecutive
+            # intersections if they are crossing the same edge. This happens, for example,
+            # when a tiny feature pokes through an edge. This helps not to set dl_min
+            # to a very low value, and take into account only actual gaps and strips.
+
+            # To do that we use the fact that intersection points are recorded and stored
+            # in the order as they appear along the border of the polygon.
+
+            # first we calculate linearized indices of edges (cells) they cross
             linear_index = cells_ij[:, 0] * len(grid_y_coords) + cells_ij[:, 1]
 
+            # then look at the differences with next and previous neighbors
             fwd_diff = linear_index - np.roll(linear_index, -1)
-            bwd_diff = linear_index - np.roll(linear_index, 1)
+            bwd_diff = np.roll(fwd_diff, 1)
 
+            # an intersection point is not a part of a "re-entering subcell feature"
+            # if it doesn't cross the same edges as its neighbors
             valid = np.logical_and(fwd_diff != 0, bwd_diff != 0)
 
             cells_dy = cells_dy[valid]
             cells_ij = cells_ij[valid]
 
-            # duplicate crossings that are too close to cell boundaries
-            close_zero = cells_dy < GAP_MESHING_TOL
-            close_one = (1.0 - cells_dy) < GAP_MESHING_TOL
-
-            duplicate_ij_zero = cells_ij[close_zero]
-            duplicate_ij_one = cells_ij[close_one]
-
-            cells_ij = np.concatenate([cells_ij, duplicate_ij_zero - np.array([0, 1]), duplicate_ij_one + np.array([0, 1])])
-            cells_dy = np.concatenate([cells_dy, (1.0 - GAP_MESHING_TOL / 10) * np.ones(len(duplicate_ij_zero)), GAP_MESHING_TOL / 10 * np.ones(len(duplicate_ij_one))])
+            # Now we are duplicating intersection points very close to cell boundaries
+            # to corresponding adjacent cells. Basically, if we have a line crossing
+            # very close to a grid node, we consider that it crosses edges on both sides
+            # from that node. That is, this serves as a tolerance allowance.
+            # Note that duplicated intersections and their originals will be snapped to
+            # cell boundaries during quantization later.
+            close_to_zero = cells_dy < GAP_MESHING_TOL
+            close_to_one = (1.0 - cells_dy) < GAP_MESHING_TOL
+
+            points_to_duplicate_near_zero = cells_ij[close_to_zero]
+            points_to_duplicate_near_one = cells_ij[close_to_one]
+
+            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_dy = np.concatenate(
+                [
+                    cells_dy,
+                    np.ones(len(points_to_duplicate_near_zero)),
+                    np.zeros(len(points_to_duplicate_near_one)),
+                ]
+            )
 
         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
     ) -> 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)
+        v_cells_ij, v_cells_dy = self._find_vertical_intersections(
+            grid_x_coords, grid_y_coords, poly_vertices
+        )
 
         # 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))
+        h_cells_ij, h_cells_dx = self._find_vertical_intersections(
+            grid_y_coords, grid_x_coords, np.flip(poly_vertices, axis=1)
+        )
         if len(h_cells_ij) > 0:
+            # flip dimensions back
             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
+        self, x, y, merged_geos
     ) -> Tuple[List[Tuple[int, int]], List[float], List[Tuple[int, int]], List[float]]:
         """Detect intersection points of geometries boundaries and grid lines."""
 
@@ -1685,6 +1723,8 @@ def _process_slice(
         # loop over all shapes
         for mat, shapes in merged_geos:
             if not mat.is_pec:
+                # note that we expect LossyMetal's converted into PEC in merged_geos
+                # that is why we are not checking for that separately
                 continue
             polygon_list = ClipOperation.to_polygon_list(shapes)
             for poly in polygon_list:
@@ -1726,12 +1766,12 @@ def _process_slice(
             h_cells_dx = np.concatenate(h_cells_dx)
 
         return v_cells_ij, v_cells_dy, h_cells_ij, h_cells_dx
-    
+
     def _generate_horizontal_snapping_lines(
-            self, grid_y_coords, intersected_cells_ij, relative_vert_disp
+        self, grid_y_coords, intersected_cells_ij, relative_vert_disp
     ) -> Tuple[List[CoordinateOptional], float]:
         """Convert a list of intersections of vertical grid lines, given as coordinates of cells
-        and relative vertical displacement inside each cell, into locations of snapping lines that 
+        and relative vertical displacement inside each cell, into locations of snapping lines that
         resolve thin gaps and strips.
         """
         min_gap_width = inf
@@ -1740,9 +1780,11 @@ def _generate_horizontal_snapping_lines(
         if len(intersected_cells_ij) > 0:
             # quantize intersection locations
             relative_vert_disp = np.round(relative_vert_disp / GAP_MESHING_TOL).astype(int)
-            cell_linear_inds = intersected_cells_ij[:, 0] * len(grid_y_coords) + intersected_cells_ij[:, 1]
+            cell_linear_inds = (
+                intersected_cells_ij[:, 0] * len(grid_y_coords) + intersected_cells_ij[:, 1]
+            )
             cell_linear_inds_and_disps = np.transpose([cell_linear_inds, relative_vert_disp])
-            # remove duplicates 
+            # remove duplicates
             cell_linear_inds_and_disps_unique = np.unique(cell_linear_inds_and_disps, axis=0)
 
             # count intersections of vertical grid lines in each cell
@@ -1767,31 +1809,47 @@ def _generate_horizontal_snapping_lines(
                 j_selection = cell_linear_inds_unique_j == ind_j
                 # and number intersections in each of them
                 counts_j = counts[j_selection]
-                
+
                 # find cell with max intersections
                 max_count_el = np.argmax(counts_j)
                 max_count = counts_j[max_count_el]
                 if max_count > 1:
                     # get its linear index
                     target_cell_linear_ind = cell_linear_inds_unique[j_selection][max_count_el]
                     # look up relative positions of intersections in that cells
-                    target_disps = np.sort(cell_linear_inds_and_disps_unique[cell_linear_inds_and_disps_unique[:, 0] == target_cell_linear_ind, 1])
+                    target_disps = np.sort(
+                        cell_linear_inds_and_disps_unique[
+                            cell_linear_inds_and_disps_unique[:, 0] == target_cell_linear_ind, 1
+                        ]
+                    )
 
                     # place a snapping line between any two neighboring intersections (in relative units)
-                    relative_snap_lines_pos = 0.5 * (target_disps[1:] + target_disps[:-1]) * GAP_MESHING_TOL
+                    relative_snap_lines_pos = (
+                        0.5 * (target_disps[1:] + target_disps[:-1]) * GAP_MESHING_TOL
+                    )
                     # convert relative positions to absolute ones
-                    snapping_lines_y += [grid_y_coords[ind_j] + rel_pos * (grid_y_coords[ind_j + 1] - grid_y_coords[ind_j]) for rel_pos in relative_snap_lines_pos]
+                    snapping_lines_y += [
+                        grid_y_coords[ind_j]
+                        + rel_pos * (grid_y_coords[ind_j + 1] - grid_y_coords[ind_j])
+                        for rel_pos in relative_snap_lines_pos
+                    ]
 
                     # compute minimal gap/strip width
-                    min_gap_width_current = np.min(target_disps[1:] - target_disps[:-1]) * GAP_MESHING_TOL
-                    min_gap_width = min(min_gap_width, min_gap_width_current * (grid_y_coords[ind_j + 1] - grid_y_coords[ind_j]))
+                    min_gap_width_current = (
+                        np.min(target_disps[1:] - target_disps[:-1]) * GAP_MESHING_TOL
+                    )
+                    min_gap_width = min(
+                        min_gap_width,
+                        min_gap_width_current * (grid_y_coords[ind_j + 1] - grid_y_coords[ind_j]),
+                    )
 
         return snapping_lines_y, min_gap_width
 
     def _resolve_gaps(self, structures: List, grid: Grid) -> 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
         )
@@ -1810,16 +1868,28 @@ def _resolve_gaps(self, structures: List, grid: Grid) -> Tuple[List[CoordinateOp
         v_cells_ij, v_cells_dy, h_cells_ij, h_cells_dx = self._process_slice(x, y, plane_slice)
 
         # generate horizontal snapping lines
-        snapping_lines_y, min_gap_width_along_y = self._generate_horizontal_snapping_lines(y, v_cells_ij, v_cells_dy)
+        snapping_lines_y, min_gap_width_along_y = self._generate_horizontal_snapping_lines(
+            y, v_cells_ij, v_cells_dy
+        )
+        detected_gap_width = min_gap_width_along_y
 
         # generate vertical snapping lines
-        snapping_lines_x, min_gap_width_along_x = self._generate_horizontal_snapping_lines(x, np.roll(h_cells_ij, shift=1, axis=1), h_cells_dx)
+        if len(h_cells_ij) > 0:  # check, otherwise np.roll fails
+            snapping_lines_x, min_gap_width_along_x = self._generate_horizontal_snapping_lines(
+                x, np.roll(h_cells_ij, shift=1, axis=1), h_cells_dx
+            )
 
-        detected_gap_width = min(min_gap_width_along_y, min_gap_width_along_x)
+            detected_gap_width = min(detected_gap_width, min_gap_width_along_x)
+        else:
+            snapping_lines_x = []
 
         # convert snapping lines' coordinates into 3d coordinates
-        snapping_lines_y_3d = [Box.unpop_axis(Y, (None, None), axis=tan_dims[1]) for Y in snapping_lines_y]
-        snapping_lines_x_3d = [Box.unpop_axis(X, (None, None), axis=tan_dims[0]) for X in snapping_lines_x]
+        snapping_lines_y_3d = [
+            Box.unpop_axis(Y, (None, None), axis=tan_dims[1]) for Y in snapping_lines_y
+        ]
+        snapping_lines_x_3d = [
+            Box.unpop_axis(X, (None, None), axis=tan_dims[0]) for X in snapping_lines_x
+        ]
 
         return snapping_lines_x_3d + snapping_lines_y_3d, detected_gap_width
 
@@ -2405,7 +2475,7 @@ def _make_grid_one_iteration(
             If `None`, recomputes internal snapping points.
         dl_min_from_gaps : pd.PositiveFloat
             Minimal grid size computed based on autodetected gaps.
-            
+
 
         Returns
         -------

tidy3d/components/simulation.py

@@ -944,7 +944,12 @@ def plot_grid(
 
         # Plot snapping points
         for points, plot_param in zip(
-            [self.internal_snapping_points, self.grid_spec.snapping_points, self._gap_meshing_snapping_lines], plot_params
+            [
+                self.internal_snapping_points,
+                self.grid_spec.snapping_points,
+                self._gap_meshing_snapping_lines,
+            ],
+            plot_params,
         ):
             for point in points:
                 _, (x_point, y_point) = Geometry.pop_axis(point, axis=axis)
@@ -1161,7 +1166,6 @@ def grid(self) -> Grid:
         # This would AutoGrid the in-plane directions of the 2D materials
         # return self._grid_corrections_2dmaterials(grid)
         return grid
-    
 
     @cached_property
     def _gap_meshing_snapping_lines(self) -> List[CoordinateOptional]: