tests

Author: dbochkov-flexcompute

tests/test_components/test_layerrefinement.py

@@ -402,3 +402,89 @@ def test_dl_min_from_smallest_feature():
     )
 
     _ = sim.grid
+
+
+def test_gap_meshing():
+    w = 1
+    l = 10
+
+    l_shape_1 = td.Structure(
+        medium=td.PECMedium(),
+        geometry=td.PolySlab(
+            axis=2,
+            slab_bounds=[0, 2],
+            vertices=[
+                [0, 0],
+                [l, 0],
+                [l, w],
+                [w, w],
+                [w, l],
+                [0, l],
+            ],
+        ),
+    )
+
+    gap = 0.1
+    l_shape_2 = td.Structure(
+        medium=td.PECMedium(),
+        geometry=td.PolySlab(
+            axis=2,
+            slab_bounds=[0, 2],
+            vertices=[
+                [w + gap, w + gap],
+                [w + gap + l, w + gap],
+                [w + gap + l, w + gap + w],
+                [w + gap + w, w + gap + w],
+                [w + gap + w, w + gap + l],
+                [w + gap, w + gap + l],
+                [0, w + gap + l],
+                [0, gap + l],
+                [w + gap, gap + l],
+            ],
+        ),
+    )
+
+    # ax = l_shape_1.plot(z=1)
+    # l_shape_2.plot(z=1, ax=ax)
+
+    for num_iters in range(2):
+        grid_spec = td.GridSpec(
+            grid_x=td.AutoGrid(min_steps_per_wvl=10),
+            grid_y=td.AutoGrid(min_steps_per_wvl=10),
+            grid_z=td.AutoGrid(min_steps_per_wvl=10),
+            layer_refinement_specs=[
+                td.LayerRefinementSpec(
+                    axis=2,
+                    corner_finder=None,
+                    gap_meshing_iters=num_iters,
+                    size=[td.inf, td.inf, 2],
+                    center=[0, 0, 1],
+                )
+            ],
+            wavelength=7,
+        )
+
+        sim = td.Simulation(
+            structures=[l_shape_1, l_shape_2],
+            grid_spec=grid_spec,
+            size=(1.2 * l, 1.2 * l, 2),
+            center=(0.5 * l, 0.5 * l, 0),
+            run_time=1e-15,
+        )
+
+        # ax = sim.plot(z=1)
+        # sim.plot_grid(z=1, ax=ax)
+        # ax.set_xlim([0, 2])
+        # ax.set_ylim([0, 2])
+
+        resolved_x = np.any(
+            np.logical_and(sim.grid.boundaries.x > w, sim.grid.boundaries.x < w + gap)
+        )
+        resolved_y = np.any(
+            np.logical_and(sim.grid.boundaries.y > w, sim.grid.boundaries.y < w + gap)
+        )
+
+        if num_iters == 0:
+            assert (not resolved_x) and (not resolved_y)
+        else:
+            assert resolved_x and resolved_y

tidy3d/components/grid/corner_finder.py

@@ -119,7 +119,6 @@ def _merged_pec_on_plane(
 
         return merged_geos
 
-
     def _corners_and_convexity(
         self,
         normal_axis: Axis,
@@ -149,7 +148,9 @@ def _corners_and_convexity(
         """
 
         # merge geometries
-        merged_geos = self._merged_pec_on_plane(normal_axis=normal_axis, coord=coord, structure_list=structure_list)
+        merged_geos = self._merged_pec_on_plane(
+            normal_axis=normal_axis, coord=coord, structure_list=structure_list
+        )
 
         # corner finder
         corner_list = []

tidy3d/components/grid/grid_spec.py

@@ -1533,36 +1533,37 @@ def _override_structures_along_axis(
             if refinement_structures[0].dl[self.axis] <= dl:
                 override_structures += refinement_structures
         return override_structures
-    
+
     def _find_vertical_intersections(self, x, y, vertices):
         """Detect intersection points of single polygon and vertical grid lines."""
 
         cells_ij = []
         cells_dy = []
         h_inds = np.argmax(x[:, None] > vertices[None, :, 0], axis=0)
         h_inds[vertices[:, 0] > x[-1]] = len(x)
-        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)):
+        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)
+        ):
             # sort vertices in ascending order
             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]
-                
+
                 # 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, 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])
-        
+
                 h_y = h_y[inds_inside_grid]
-        
+
                 # 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
-        
+
                 # find local dy
                 dy = h_y - y[cell_js]
 
@@ -1576,7 +1577,7 @@ def _find_vertical_intersections(self, x, y, vertices):
             cells_dy = np.concatenate(cells_dy)
 
         return cells_ij, cells_dy
-    
+
     def _process_poly(self, x, y, vertices):
         """Detect intersection points of single polygon and grid lines."""
         v_cells_ij, v_cells_dy = self._find_vertical_intersections(x, y, vertices)
@@ -1585,7 +1586,7 @@ def _process_poly(self, x, y, vertices):
         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."""
         v_cells_ij = []
@@ -1600,15 +1601,15 @@ def _process_slice(self, x, y, merged_geos):
             for poly in polygon_list:
                 poly = poly.normalize().buffer(0)
                 data = self._process_poly(x, y, np.array(poly.exterior.coords)[:-1])
-                
+
                 if len(data[0]) > 0:
                     v_cells_ij.append(data[0])
                     v_cells_dy.append(data[1])
 
                 if len(data[2]) > 0:
                     h_cells_ij.append(data[2])
                     h_cells_dx.append(data[3])
-                
+
                 # in case the polygon has holes
                 for poly_inner in poly.interiors:
                     data = self._process_poly(x, y, np.array(poly_inner.coords)[:-1])
@@ -1633,13 +1634,14 @@ def _process_slice(self, x, y, merged_geos):
             h_cells_dx = np.concatenate(h_cells_dx)
 
         return v_cells_ij, v_cells_dy, h_cells_ij, h_cells_dx
-    
 
     def _resolve_gaps(self, structures: List, grid: Grid):
         """Detect underresolved gaps and place snapping lines in them. Also return the detected minimal gap width."""
 
         # get merged pec structures on plane
-        plane_slice = CornerFinderSpec._merged_pec_on_plane(coord=self.center_axis, normal_axis=self.axis, structure_list=structures)
+        plane_slice = CornerFinderSpec._merged_pec_on_plane(
+            coord=self.center_axis, normal_axis=self.axis, structure_list=structures
+        )
 
         _, tan_dims = Box.pop_axis([0, 1, 2], self.axis)
         x = grid.boundaries.to_list[tan_dims[0]]
@@ -1654,7 +1656,9 @@ def _resolve_gaps(self, structures: List, grid: Grid):
         snapping_lines_y = []
         if len(v_cells_ij) > 0:
             # add number of intersections
-            ind_unique, counts = np.unique(v_cells_ij[:, 0] * len(y) + v_cells_ij[:, 1], return_counts=True)
+            ind_unique, counts = np.unique(
+                v_cells_ij[:, 0] * len(y) + v_cells_ij[:, 1], return_counts=True
+            )
             ind_unique_j = ind_unique % len(y)
             for ind_j in np.unique(ind_unique_j):
                 max_count = np.max(counts[ind_unique_j == ind_j])
@@ -1669,7 +1673,9 @@ def _resolve_gaps(self, structures: List, grid: Grid):
         snapping_lines_x = []
         if len(h_cells_ij) > 0:
             # add number of intersections
-            ind_unique, counts = np.unique(h_cells_ij[:, 0] * len(y) + h_cells_ij[:, 1], return_counts=True)
+            ind_unique, counts = np.unique(
+                h_cells_ij[:, 0] * len(y) + h_cells_ij[:, 1], return_counts=True
+            )
             ind_unique_i = ind_unique // len(y)
             for ind_i in np.unique(ind_unique_i):
                 max_count = np.max(counts[ind_unique_i == ind_i])
@@ -1680,7 +1686,9 @@ def _resolve_gaps(self, structures: List, grid: Grid):
                     snapping_lines_x.append(x[ind_i] + 0.5 * (min_dx + max_dx))
                     detected_gap_width = min(detected_gap_width, max_dy - min_dy)
 
-        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], detected_gap_width
+        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
+        ], detected_gap_width
 
 
 class GridSpec(Tidy3dBaseModel):
@@ -2067,7 +2075,7 @@ def get_wavelength(self, sources: List[SourceType]) -> float:
             wavelength = self.wavelength_from_sources(sources)
             log.info(f"Auto meshing using wavelength {wavelength:1.4f} defined from sources.")
         return wavelength
-    
+
     def make_grid(
         self,
         structures: List[Structure],
@@ -2118,15 +2126,19 @@ def make_grid(
         )
 
         if len(self.layer_refinement_specs) > 0:
-            num_iters = max(layer_spec.gap_meshing_iters for layer_spec in self.layer_refinement_specs)
-            
+            num_iters = max(
+                layer_spec.gap_meshing_iters for layer_spec in self.layer_refinement_specs
+            )
+
             snapping_lines = []
             min_gap_width = inf
             for ind in range(num_iters):
                 new_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)
+                        one_layer_snapping_lines, gap_width = layer_spec._resolve_gaps(
+                            structures, old_grid
+                        )
                         new_snapping_lines = new_snapping_lines + one_layer_snapping_lines
                         if layer_spec.dl_min_from_gap_width:
                             min_gap_width = min(min_gap_width, gap_width)
@@ -2158,7 +2170,6 @@ def make_grid(
                 old_grid = new_grid
 
         return old_grid
-    
 
     def _make_grid_one_iteration(
         self,