min size feature autodetect

Author: dbochkov-flexcompute

CHANGELOG.md

@@ -13,6 +13,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - The method `Geometry.reflected` can be used to create a reflected copy of any geometry off a plane. As for other transformations, for efficiency, `reflected` `PolySlab` directly returns an updated `PolySlab` object rather than a `Transformed` object, except when the normal of the plane of reflection has a non-zero component along the slab axis, in which case `Transformed` is still returned.
 - Validation check for unit error in grid spacing.
 - Validation that when symmetry is imposed along a given axis, the boundary conditions on each side of the axis are identical.
+- Fields `convex_resolution`, `concave_resolution`, and `mixed_resolution` in `CornerFinderSpec` can be used to take into account the dimensions of autodetected convex, concave, or mixed geometric features when `dl_min` is automatically inferred during automatic grid generation.
 
 ### Changed
 - Supplying autograd-traced values to geometric fields (`center`, `size`) of simulations, monitors, and sources now logs a warning and falls back to the static value instead of erroring.

tests/test_components/test_layerrefinement.py

@@ -339,3 +339,66 @@ def count_grids_within_gap(sim_t):
     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
+
+
+def test_dl_min_from_smallest_feature():
+    structure = td.Structure(
+        geometry=td.PolySlab(
+            vertices=[
+                [0, 0],
+                [2, 0],
+                [2, 1],
+                [1, 1],
+                [1, 1.1],
+                [2, 1.1],
+                [2, 2],
+                [1, 2],
+                [1, 2.2],
+                [0.7, 2.2],
+                [0.7, 2],
+                [0, 2],
+            ],
+            slab_bounds=[-1, 1],
+            axis=2,
+        ),
+        medium=td.PECMedium(),
+    )
+
+    # check expected dl_min
+    layer_spec = td.LayerRefinementSpec(
+        axis=2,
+        size=(td.inf, td.inf, 2),
+        corner_finder=td.CornerFinderSpec(
+            convex_resolution=10,
+        ),
+    )
+    dl_min = layer_spec._dl_min_from_smallest_feature([structure])
+    assert np.allclose(0.3 / 10, dl_min)
+
+    layer_spec = td.LayerRefinementSpec(
+        axis=2,
+        size=(td.inf, td.inf, 2),
+        corner_finder=td.CornerFinderSpec(mixed_resolution=10),
+    )
+    dl_min = layer_spec._dl_min_from_smallest_feature([structure])
+    assert np.allclose(0.2 / 10, dl_min)
+
+    layer_spec = td.LayerRefinementSpec(
+        axis=2,
+        size=(td.inf, td.inf, 2),
+        corner_finder=td.CornerFinderSpec(
+            concave_resolution=10,
+        ),
+    )
+    dl_min = layer_spec._dl_min_from_smallest_feature([structure])
+    assert np.allclose(0.1 / 10, dl_min)
+
+    # check grid is generated succesfully
+    sim = td.Simulation(
+        size=(5, 5, 5),
+        structures=[structure],
+        grid_spec=td.GridSpec.auto(layer_refinement_specs=[layer_spec], wavelength=100 * td.C_0),
+        run_time=1e-20,
+    )
+
+    _ = sim.grid

tidy3d/components/grid/corner_finder.py

@@ -1,17 +1,17 @@
 """Find corners of structures on a 2D plane."""
 
-from typing import List, Literal, Optional
+from typing import List, Literal, Optional, Tuple
 
 import numpy as np
 import pydantic.v1 as pd
 
 from ...constants import inf
-from ..base import Tidy3dBaseModel
+from ..base import Tidy3dBaseModel, cached_property
 from ..geometry.base import Box, ClipOperation
 from ..geometry.utils import merging_geometries_on_plane
 from ..medium import PEC, LossyMetalMedium
 from ..structure import Structure
-from ..types import ArrayFloat2D, Axis
+from ..types import ArrayFloat1D, ArrayFloat2D, Axis
 
 CORNER_ANGLE_THRESOLD = 0.1 * np.pi
 
@@ -44,12 +44,44 @@ class CornerFinderSpec(Tidy3dBaseModel):
         "is below the threshold value based on Douglas-Peucker algorithm, the vertex is disqualified as a corner.",
     )
 
-    def corners(
+    concave_resolution: Optional[pd.PositiveInt] = pd.Field(
+        None,
+        title="Concave Region Resolution.",
+        description="Specifies number of steps to use for determining `dl_min` based on concave featues."
+        "If set to ``None``, then the corresponding `dl_min` reduction is not applied.",
+    )
+
+    convex_resolution: Optional[pd.PositiveInt] = pd.Field(
+        None,
+        title="Convex Region Resolution.",
+        description="Specifies number of steps to use for determining `dl_min` based on convex featues."
+        "If set to ``None``, then the corresponding `dl_min` reduction is not applied.",
+    )
+
+    mixed_resolution: Optional[pd.PositiveInt] = pd.Field(
+        None,
+        title="Mixed Region Resolution.",
+        description="Specifies number of steps to use for determining `dl_min` based on mixed featues."
+        "If set to ``None``, then the corresponding `dl_min` reduction is not applied.",
+    )
+
+    @cached_property
+    def _no_min_dl_override(self):
+        return all(
+            (
+                self.concave_resolution is None,
+                self.convex_resolution is None,
+                self.mixed_resolution is None,
+            )
+        )
+
+    def _corners_and_convexity(
         self,
         normal_axis: Axis,
         coord: float,
         structure_list: List[Structure],
-    ) -> ArrayFloat2D:
+        ravel: bool,
+    ) -> Tuple[ArrayFloat2D, ArrayFloat1D]:
         """On a 2D plane specified by axis = `normal_axis` and coordinate `coord`, find out corners of merged
         geometries made of `medium`.
 
@@ -62,6 +94,8 @@ def corners(
             Position of plane along the normal axis.
         structure_list : List[Structure]
             List of structures present in simulation.
+        ravel : bool
+            Whether to put the resulting corners in a single list or per polygon.
 
         Returns
         -------
@@ -89,6 +123,7 @@ def corners(
 
         # corner finder
         corner_list = []
+        convexity_list = []
         for mat, shapes in merged_geos:
             if self.medium != "all" and mat.is_pec != (self.medium == "metal"):
                 continue
@@ -97,17 +132,59 @@ def corners(
                 poly = poly.normalize().buffer(0)
                 if self.distance_threshold is not None:
                     poly = poly.simplify(self.distance_threshold, preserve_topology=True)
-                corner_list.append(self._filter_collinear_vertices(list(poly.exterior.coords)))
+                corners_xy, corners_convexity = self._filter_collinear_vertices(
+                    list(poly.exterior.coords)
+                )
+                corner_list.append(corners_xy)
+                convexity_list.append(corners_convexity)
                 # in case the polygon has holes
                 for poly_inner in poly.interiors:
-                    corner_list.append(self._filter_collinear_vertices(list(poly_inner.coords)))
+                    corners_xy, corners_convexity = self._filter_collinear_vertices(
+                        list(poly_inner.coords)
+                    )
+                    corner_list.append(corners_xy)
+                    convexity_list.append(corners_convexity)
 
-        if len(corner_list) > 0:
+        if ravel and len(corner_list) > 0:
             corner_list = np.concatenate(corner_list)
+            convexity_list = np.concatenate(convexity_list)
+
+        return corner_list, convexity_list
+
+    def corners(
+        self,
+        normal_axis: Axis,
+        coord: float,
+        structure_list: List[Structure],
+    ) -> ArrayFloat2D:
+        """On a 2D plane specified by axis = `normal_axis` and coordinate `coord`, find out corners of merged
+        geometries made of `medium`.
+
+
+        Parameters
+        ----------
+        normal_axis : Axis
+            Axis normal to the 2D plane.
+        coord : float
+            Position of plane along the normal axis.
+        structure_list : List[Structure]
+            List of structures present in simulation.
+
+        Returns
+        -------
+        ArrayFloat2D
+            Corner coordinates.
+        """
+
+        corner_list, _ = self._corners_and_convexity(
+            normal_axis=normal_axis, coord=coord, structure_list=structure_list, ravel=True
+        )
         return corner_list
 
-    def _filter_collinear_vertices(self, vertices: ArrayFloat2D) -> ArrayFloat2D:
-        """Filter collinear vertices of a polygon, and return corners.
+    def _filter_collinear_vertices(
+        self, vertices: ArrayFloat2D
+    ) -> Tuple[ArrayFloat2D, ArrayFloat1D]:
+        """Filter collinear vertices of a polygon, and return corners locations and their convexity.
 
         Parameters
         ----------
@@ -119,6 +196,8 @@ def _filter_collinear_vertices(self, vertices: ArrayFloat2D) -> ArrayFloat2D:
         -------
         ArrayFloat2D
             Corner coordinates.
+        ArrayFloat1D
+            Convexity of corners: True for outer corners, False for inner corners.
         """
 
         def normalize(v):
@@ -136,5 +215,12 @@ def normalize(v):
         inner_product = np.where(inner_product > 1, 1, inner_product)
         inner_product = np.where(inner_product < -1, -1, inner_product)
         angle = np.arccos(inner_product)
+        num_vs = len(vs_orig)
+        cross_product = np.cross(
+            np.hstack([unit_next, np.zeros((num_vs, 1))]),
+            np.hstack([unit_previous, np.zeros((num_vs, 1))]),
+            axis=-1,
+        )
+        convexity = cross_product[:, 2] < 0
         ind_filter = angle <= np.pi - self.angle_threshold
-        return vs_orig[ind_filter]
+        return vs_orig[ind_filter], convexity[ind_filter]

tidy3d/components/grid/grid_spec.py

@@ -1262,7 +1262,7 @@ def _unpop_axis(self, ax_coord: float, plane_coord: Any) -> CoordinateOptional:
         """
         return self.unpop_axis(ax_coord, [plane_coord, plane_coord], self.axis)
 
-    def suggested_dl_min(self, grid_size_in_vacuum: float) -> float:
+    def suggested_dl_min(self, grid_size_in_vacuum: float, structures: List[Structure]) -> float:
         """Suggested lower bound of grid step size for this layer.
 
         Parameters
@@ -1293,6 +1293,12 @@ def suggested_dl_min(self, grid_size_in_vacuum: float) -> float:
         # inplane dimension
         if self.corner_finder is not None and self.corner_refinement is not None:
             dl_min = min(dl_min, self.corner_refinement._grid_size(grid_size_in_vacuum))
+
+        # min feature size
+        if self.corner_finder is not None and not self.corner_finder._no_min_dl_override:
+            dl_suggested = self._dl_min_from_smallest_feature(structures)
+            dl_min = min(dl_min, dl_suggested)
+
         return dl_min
 
     def generate_snapping_points(self, structure_list: List[Structure]) -> List[CoordinateOptional]:
@@ -1329,22 +1335,80 @@ def _inplane_inside(self, point: ArrayFloat2D) -> bool:
         )
         return self.inside(point_3d[0], point_3d[1], point_3d[2])
 
-    def _corners(self, structure_list: List[Structure]) -> List[CoordinateOptional]:
-        """Inplane corners in 3D coordinate."""
+    def _corners_and_convexity_2d(
+        self, structure_list: List[Structure], ravel: bool
+    ) -> List[CoordinateOptional]:
+        """Raw inplane corners and their convexity."""
         if self.corner_finder is None:
-            return []
+            return [], []
 
         # filter structures outside the layer
         structures_intersect = structure_list
         if self._is_inplane_bounded:
             structures_intersect = [s for s in structure_list if self.intersects(s.geometry)]
-        inplane_points = self.corner_finder.corners(
-            self.axis, self.center_axis, structures_intersect
+        inplane_points, convexity = self.corner_finder._corners_and_convexity(
+            self.axis, self.center_axis, structures_intersect, ravel
         )
 
         # filter corners outside the inplane bounds
-        if self._is_inplane_bounded:
-            inplane_points = [point for point in inplane_points if self._inplane_inside(point)]
+        if self._is_inplane_bounded and len(inplane_points) > 0:
+            # flatten temporary list of arrays for faster processing
+            if not ravel:
+                split_inds = np.cumsum([len(pts) for pts in inplane_points])[:-1]
+                inplane_points = np.concatenate(inplane_points)
+                convexity = np.concatenate(convexity)
+            inds = [self._inplane_inside(point) for point in inplane_points]
+            inplane_points = inplane_points[inds]
+            convexity = convexity[inds]
+            if not ravel:
+                inplane_points = np.split(inplane_points, split_inds)
+                convexity = np.split(convexity, split_inds)
+
+        return inplane_points, convexity
+
+    def _dl_min_from_smallest_feature(self, structure_list: List[Structure]):
+        """Calculate `dl_min` suggestion based on smallest feature size."""
+
+        inplane_points, convexity = self._corners_and_convexity_2d(
+            structure_list=structure_list, ravel=False
+        )
+
+        dl_min = inf
+
+        if self.corner_finder is None or self.corner_finder._no_min_dl_override:
+            return dl_min
+
+        finder = self.corner_finder
+
+        for points, conv in zip(inplane_points, convexity):
+            conv_nei = np.roll(conv, -1)
+            lengths = np.linalg.norm(points - np.roll(points, axis=0, shift=-1), axis=-1)
+
+            if finder.convex_resolution is not None:
+                convex_features = np.logical_and(conv, conv_nei)
+                if np.any(convex_features):
+                    min_convex_size = np.min(lengths[convex_features])
+                    dl_min = min(dl_min, min_convex_size / finder.convex_resolution)
+
+            if finder.concave_resolution is not None:
+                concave_features = np.logical_not(np.logical_or(conv, conv_nei))
+                if np.any(concave_features):
+                    min_concave_size = np.min(lengths[concave_features])
+                    dl_min = min(dl_min, min_concave_size / finder.concave_resolution)
+
+            if finder.mixed_resolution is not None:
+                mixed_features = np.logical_xor(conv, conv_nei)
+                if np.any(mixed_features):
+                    min_mixed_size = np.min(lengths[mixed_features])
+                    dl_min = min(dl_min, min_mixed_size / finder.mixed_resolution)
+
+        return dl_min
+
+    def _corners(self, structure_list: List[Structure]) -> List[CoordinateOptional]:
+        """Inplane corners in 3D coordinate."""
+        inplane_points, _ = self._corners_and_convexity_2d(
+            structure_list=structure_list, ravel=True
+        )
 
         # convert 2d points to 3d
         return [
@@ -1817,7 +1881,7 @@ def _dl_min(
         if self.layer_refinement_used:
             min_vacuum_dl = self._min_vacuum_dl_in_autogrid(wavelength, sim_size)
             for layer in self.layer_refinement_specs:
-                min_dl = min(min_dl, layer.suggested_dl_min(min_vacuum_dl))
+                min_dl = min(min_dl, layer.suggested_dl_min(min_vacuum_dl, structures))
         # from lumped elements
         for lumped_element in lumped_elements:
             for override_structure in lumped_element.to_mesh_overrides():