More efficient mesh refinement in LayerRefinementSpec for interior-disjoint structures

Author: weiliangjin2021

CHANGELOG.md

@@ -14,6 +14,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added support for `tidy3d-extras`, an optional plugin that enables more accurate local mode solving via subpixel averaging.
 
 ### Changed
+- `LayerRefinementSpec` defaults to assuming structures made of different materials are interior-disjoint for more efficient mesh generation.
 
 ### Fixed
 - Stricter validation for `bend_radius` in mode simulations, preventing the bend center from coinciding with the simulation boundary.

tests/test_components/test_layerrefinement.py

@@ -288,6 +288,62 @@ def count_grids_within_layer(sim_t):
     )
 
 
+def test_grid_spec_with_layers_interior_disjoint():
+    """Test the application of layer_specs with interior_disjoint assumption to GridSpec."""
+
+    thickness = 1e-3
+    lumped_elements = []
+    # a PEC thin layer structure
+    box1 = td.Box(size=(thickness, 2, 2))
+    box2 = td.Box(center=(0, -1, 0), size=(thickness, 1, 1))
+    pec_str = td.Structure(geometry=box1 - box2, medium=td.PEC)
+    # a dielectric that overlaps with PEC structure, hence breaking interior-disjoint assumption
+    die_str = td.Structure(
+        geometry=td.Box(size=(thickness, 0.5, 0.5)),
+        medium=td.Medium(),
+    )
+    # layer spec assuming interior disjoint geometries
+    layer_disjoint = LayerRefinementSpec.from_structures(
+        [
+            pec_str,
+        ],
+        interior_disjoint_geometries=True,
+    )
+    # layer spec for general geometries
+    layer_general = layer_disjoint.updated_copy(interior_disjoint_geometries=False)
+
+    sim = td.Simulation(
+        size=(4, 4, 4),
+        grid_spec=td.GridSpec.auto(
+            min_steps_per_wvl=11, wavelength=1, layer_refinement_specs=[layer_disjoint]
+        ),
+        boundary_spec=td.BoundarySpec.pml(),
+        structures=[pec_str],
+        run_time=1e-12,
+    )
+    sim_general = sim.updated_copy(
+        grid_spec=sim.grid_spec.updated_copy(layer_refinement_specs=[layer_general])
+    )
+
+    # for simulations with interior-disjoint geometries, same corner finder results for both layer settings
+    def is_equal_snapping_points(plist1, plist2):
+        if len(plist1) != len(plist2):
+            return False
+        for p1, p2 in zip(plist1, plist2):
+            if p1 != p2:
+                return False
+        return True
+
+    assert is_equal_snapping_points(
+        sim.internal_snapping_points, sim_general.internal_snapping_points
+    )
+
+    # for simulations with overlapping geometries, assuming interior-disjoint misses some corners in this simulation
+    sim = sim.updated_copy(structures=[pec_str, die_str])
+    sim_general = sim_general.updated_copy(structures=[pec_str, die_str])
+    assert len(sim_general.internal_snapping_points) > len(sim.internal_snapping_points)
+
+
 @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."""
@@ -573,6 +629,7 @@ def test_gap_meshing():
                     corner_finder=None,
                     gap_meshing_iters=2,
                     dl_min_from_gap_width=True,
+                    interior_disjoint_geometries=False,
                 )
             ],
             min_steps_per_wvl=10,

tidy3d/components/geometry/utils.py

@@ -2,12 +2,14 @@
 
 from __future__ import annotations
 
+from collections import defaultdict
 from enum import Enum
 from math import isclose
 from typing import Any, Optional, Union
 
 import numpy as np
 import pydantic
+import shapely
 
 from tidy3d.components.autograd.utils import get_static
 from tidy3d.components.base import Tidy3dBaseModel
@@ -45,6 +47,7 @@ def merging_geometries_on_plane(
     geometries: list[GeometryType],
     plane: Box,
     property_list: list[Any],
+    interior_disjoint_geometries: bool = False,
 ) -> list[tuple[Any, Shapely]]:
     """Compute list of shapes on plane. Overlaps are removed or merged depending on
     provided property_list.
@@ -57,6 +60,8 @@ def merging_geometries_on_plane(
         Plane specification.
     property_list : List = None
         Property value for each structure.
+    interior_disjoint_geometries: bool = False
+        If ``True``, geometries of different properties on the plane must not be overlapping.
 
     Returns
     -------
@@ -78,6 +83,20 @@ def merging_geometries_on_plane(
         for shape in shapes_plane:
             shapes.append((prop, shape, shape.bounds))
 
+    if interior_disjoint_geometries:
+        # No need to consider overlapping. We simply group shapes by property, and union_all
+        # shapes of the same property.
+        shapes_by_prop = defaultdict(list)
+        for prop, shape, _ in shapes:
+            shapes_by_prop[prop].append(shape)
+        # union shapes of same property
+        results = []
+        for prop, shapes in shapes_by_prop.items():
+            unionized = shapely.union_all(shapes).buffer(0).normalize()
+            if not unionized.is_empty:
+                results.append((prop, unionized))
+        return results
+
     background_shapes = []
     for prop, shape, bounds in shapes:
         minx, miny, maxx, maxy = bounds

tidy3d/components/grid/corner_finder.py

@@ -85,6 +85,7 @@ def _merged_pec_on_plane(
         structure_list: list[Structure],
         center: tuple[float, float] = [0, 0, 0],
         size: tuple[float, float, float] = [inf, inf, inf],
+        interior_disjoint_geometries: bool = False,
     ) -> list[tuple[Any, Shapely]]:
         """On a 2D plane specified by axis = `normal_axis` and coordinate `coord`, merge geometries made of PEC.
 
@@ -100,6 +101,8 @@ def _merged_pec_on_plane(
             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)
+        interior_disjoint_geometries: bool = False
+            If ``True``, geometries on the plane must not be overlapping.
 
         Returns
         -------
@@ -123,7 +126,9 @@ def _merged_pec_on_plane(
             PEC if (mat.is_pec or isinstance(mat, LossyMetalMedium)) else mat for mat in medium_list
         ]
         # merge geometries
-        merged_geos = merging_geometries_on_plane(geometry_list, plane, medium_list)
+        merged_geos = merging_geometries_on_plane(
+            geometry_list, plane, medium_list, interior_disjoint_geometries
+        )
 
         return merged_geos
 
@@ -133,6 +138,7 @@ def _corners_and_convexity(
         coord: float,
         structure_list: list[Structure],
         ravel: bool,
+        interior_disjoint_geometries: bool = False,
     ) -> tuple[ArrayFloat2D, ArrayFloat1D]:
         """On a 2D plane specified by axis = `normal_axis` and coordinate `coord`, find out corners of merged
         geometries made of PEC.
@@ -148,6 +154,8 @@ def _corners_and_convexity(
             List of structures present in simulation.
         ravel : bool
             Whether to put the resulting corners in a single list or per polygon.
+        interior_disjoint_geometries: bool = False
+            If ``True``, geometries made of different materials on the plane must not be overlapping.
 
         Returns
         -------
@@ -157,7 +165,10 @@ def _corners_and_convexity(
 
         # merge geometries
         merged_geos = self._merged_pec_on_plane(
-            normal_axis=normal_axis, coord=coord, structure_list=structure_list
+            normal_axis=normal_axis,
+            coord=coord,
+            structure_list=structure_list,
+            interior_disjoint_geometries=interior_disjoint_geometries,
         )
 
         # corner finder
@@ -183,18 +194,22 @@ def _corners_and_convexity(
                     )
                     corner_list.append(corners_xy)
                     convexity_list.append(corners_convexity)
+        return self._ravel_corners_and_convexity(ravel, corner_list, convexity_list)
 
+    def _ravel_corners_and_convexity(
+        self, ravel: bool, corner_list, convexity_list
+    ) -> tuple[ArrayFloat2D, ArrayFloat1D]:
+        """Whether to put the resulting corners in a single list or per polygon."""
         if ravel and len(corner_list) > 0:
-            corner_list = np.concatenate(corner_list)
-            convexity_list = np.concatenate(convexity_list)
-
+            return np.concatenate(corner_list), np.concatenate(convexity_list)
         return corner_list, convexity_list
 
     def corners(
         self,
         normal_axis: Axis,
         coord: float,
         structure_list: list[Structure],
+        interior_disjoint_geometries: bool = False,
     ) -> ArrayFloat2D:
         """On a 2D plane specified by axis = `normal_axis` and coordinate `coord`, find out corners of merged
         geometries made of `medium`.
@@ -208,15 +223,20 @@ def corners(
             Position of plane along the normal axis.
         structure_list : List[Structure]
             List of structures present in simulation.
-
+        interior_disjoint_geometries: bool = False
+            If ``True``, geometries made of different materials on the plane must not be overlapping.
         Returns
         -------
         ArrayFloat2D
             Corner coordinates.
         """
 
         corner_list, _ = self._corners_and_convexity(
-            normal_axis=normal_axis, coord=coord, structure_list=structure_list, ravel=True
+            normal_axis=normal_axis,
+            coord=coord,
+            structure_list=structure_list,
+            ravel=True,
+            interior_disjoint_geometries=interior_disjoint_geometries,
         )
         return corner_list