Add LayerRefinementSpec for automatic mesh refinement in layered stru…

…cture

Snapping points with optional coordinate

Shift _filter_structures_plane to geometry/utils.py

Automatically sets dl_min if not set yet

Revise and improve plot_grid visualization

Change linestyle for override structures in plot_grid

CHANGELOG.md

@@ -9,6 +9,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 - Advanced option `dist_type` that allows `LinearLumpedElement` to be distributed across grid cells in different ways. For example, restricting the network portion to a single cell and using PEC wires to connect to the desired location of the terminals. 
+- New `layer_refinement_specs` field in `GridSpec` that takes a list of `LayerRefinementSpec` for automatic mesh refinement and snapping in layered structures. Structure corners on the cross section perpendicular to layer thickness direction can be automatically identified. Mesh is automatically snapped and refined around those corners.
+
+### Changed
+- The coordinate of snapping points in `GridSpec` can take value `None`, so that mesh can be selectively snapped only along certain dimensions.
 
 ## [2.8.0rc2] - 2025-01-28
 
@@ -35,6 +39,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 -`HeatChargeSimulation` supersedes `HeatSimulation`. Though both of the can be used for Heat simulation interchangeably, the latter has been deprecated and will disappear in the future. 
 - `FluidSpec` and `SolidSpec` are now deprecated in favor of `FluidMedium` and `SolidMedium`. Both can still be used interchangeably.
 - Exclude `FluxMonitor` frequencies from adjoint design region gradient monitors since we do not differentiate through `FluxData`
+- The coordinate of snapping points in `GridSpec` can take value `None`, so that mesh can be selectively snapped only along certain dimensions.
 
 ### Fixed
 - NumPy 2.1 compatibility issue where `numpy.float64` values passed to xarray interpolation would raise TypeError.

docs/api/discretization.rst

@@ -17,3 +17,6 @@ Discretization
    tidy3d.FieldGrid
    tidy3d.YeeGrid
    tidy3d.Grid
+   tidy3d.LayerRefinementSpec
+   tidy3d.GridRefinement
+   tidy3d.CornerFinderSpec

tests/test_components/test_grid_spec.py

@@ -349,6 +349,11 @@ def test_zerosize_dimensions():
 def test_custom_grid_boundaries():
     custom = td.CustomGridBoundaries(coords=np.linspace(-1, 1, 11))
     grid_spec = td.GridSpec(grid_x=custom, grid_y=custom, grid_z=custom)
+
+    # estimated minimal step size
+    estimated_dl = grid_spec.grid_x.estimated_min_dl(wavelength=1, structure_list=[], sim_size=[])
+    assert np.isclose(estimated_dl, 0.2)
+
     source = td.PointDipole(
         source_time=td.GaussianPulse(freq0=3e14, fwidth=1e14), polarization="Ex"
     )
@@ -400,6 +405,30 @@ def test_small_sim_with_min_steps_per_sim_size():
     assert sim.num_cells > 500
 
 
+def test_autogrid_estimated_dl():
+    """Test that estimiated minimal step size in AutoGrid works as expected."""
+    box = td.Structure(
+        geometry=td.Box(size=(1, 1, 1), center=(0, 0, 1)), medium=td.Medium(permittivity=4)
+    )
+    sim_size = [10, 10, 10]
+    wavelength = 1.0
+    grid_spec = td.AutoGrid(min_steps_per_wvl=10)
+
+    # decided by medium
+    estimated = grid_spec.estimated_min_dl(wavelength, [box], sim_size)
+    assert np.isclose(estimated, wavelength / 10 / 2)
+
+    # overridden by dl_min
+    grid_spec_min = td.AutoGrid(min_steps_per_wvl=10, dl_min=0.1)
+    estimated = grid_spec_min.estimated_min_dl(wavelength, [box], sim_size)
+    assert np.isclose(estimated, 0.1)
+
+    # decided by sim_size
+    sim_size = [0.1, 0.1, 0.1]
+    estimated = grid_spec.estimated_min_dl(wavelength, [box], sim_size)
+    assert np.isclose(estimated, 0.1 / 10)
+
+
 def test_quasiuniform_grid():
     """Test grid is quasi-uniform that adjusts to structure boundaries."""
     box = td.Structure(
@@ -420,7 +449,29 @@ def test_quasiuniform_grid():
     # add snapping points
     pos = 0.281
     snapping_points = [(pos, pos, pos)]
-    sim = sim.updated_copy(
+    sim2 = sim.updated_copy(
         grid_spec=td.GridSpec.quasiuniform(dl=0.1, snapping_points=snapping_points)
     )
-    assert any(np.isclose(sim.grid.boundaries.x, pos))
+    assert any(np.isclose(sim2.grid.boundaries.x, pos))
+
+    # override structures of dl=None takes no effect
+    sim3 = sim.updated_copy(
+        grid_spec=td.GridSpec.quasiuniform(
+            dl=0.1,
+            override_structures=[
+                td.MeshOverrideStructure(
+                    geometry=td.Box(size=(0.2, 0.2, 0.2)), dl=[None, None, None]
+                ),
+            ],
+        )
+    )
+    np.allclose(sim.grid.boundaries.x, sim3.grid.boundaries.x)
+
+    # a larger dl_min can override step size
+    grid_1d = td.QuasiUniformGrid(dl=0.1, dl_min=0.2)
+    sim4 = sim.updated_copy(grid_spec=td.GridSpec(grid_x=grid_1d, grid_y=grid_1d, grid_z=grid_1d))
+    assert np.all(sim4.grid.sizes.x > 0.11)
+
+    # 2d simulation
+    sim5 = sim.updated_copy(size=(0, 1, 1))
+    assert np.isclose(sim5.grid.sizes.x, 0.1)

tests/test_components/test_layerrefinement.py

@@ -0,0 +1,268 @@
+"""Tests 2d corner finder."""
+
+import numpy as np
+import pydantic.v1 as pydantic
+import pytest
+import tidy3d as td
+from tidy3d.components.grid.corner_finder import CornerFinderSpec
+from tidy3d.components.grid.grid_spec import GridRefinement, LayerRefinementSpec
+
+CORNER_FINDER = CornerFinderSpec()
+GRID_REFINEMENT = GridRefinement()
+LAYER_REFINEMENT = LayerRefinementSpec(axis=2, size=(td.inf, td.inf, 2))
+LAYER2D_REFINEMENT = LayerRefinementSpec(axis=2, size=(td.inf, td.inf, 0))
+
+
+def test_2dcorner_finder_filter_collinear_vertex():
+    """In corner finder, test that collinear vertices are filtered"""
+    # 2nd and 3rd vertices are on a collinear line
+    vertices = ((0, 0), (0.1, 0), (0.5, 0), (1, 0), (1, 1))
+    polyslab = td.PolySlab(vertices=vertices, axis=2, slab_bounds=[-1, 1])
+    structures = [td.Structure(geometry=polyslab, medium=td.PEC)]
+    corners = CORNER_FINDER.corners(normal_axis=2, coord=0, structure_list=structures)
+    assert len(corners) == 3
+
+    # if angle threshold is 0, collinear vertex will not be filtered
+    corner_finder = CORNER_FINDER.updated_copy(angle_threshold=0)
+    corners = corner_finder.corners(normal_axis=2, coord=0, structure_list=structures)
+    assert len(corners) == 5
+
+
+def test_2dcorner_finder_filter_nearby_vertex():
+    """In corner finder, test that vertices that are very close are filtered"""
+    # filter duplicate vertices
+    vertices = ((0, 0), (0, 0), (1e-4, -1e-4), (1, 0), (1, 1))
+    polyslab = td.PolySlab(vertices=vertices, axis=2, slab_bounds=[-1, 1])
+    structures = [td.Structure(geometry=polyslab, medium=td.PEC)]
+    corners = CORNER_FINDER.corners(normal_axis=2, coord=0, structure_list=structures)
+    assert len(corners) == 4
+
+    # filter very close vertices
+    corner_finder = CORNER_FINDER.updated_copy(distance_threshold=2e-4)
+    corners = corner_finder.corners(normal_axis=2, coord=0, structure_list=structures)
+    assert len(corners) == 3
+
+
+def test_2dcorner_finder_medium():
+    """No corner found if the medium is dielectric while asking to search for corner of metal."""
+    structures = [td.Structure(geometry=td.Box(size=(1, 1, 1)), medium=td.Medium())]
+    corners = CORNER_FINDER.corners(normal_axis=2, coord=0, structure_list=structures)
+    assert len(corners) == 0
+
+
+def test_2dcorner_finder_polygon_with_hole():
+    """Find corners related to interior holes of a polygon."""
+    structures = [
+        td.Structure(geometry=td.Box(size=(2, 2, 2)), medium=td.PEC),
+        td.Structure(geometry=td.Box(size=(1, 1, 1)), medium=td.Medium()),
+    ]
+    corners = CORNER_FINDER.corners(normal_axis=2, coord=0, structure_list=structures)
+    # 4 interior, 4 exterior
+    assert len(corners) == 8
+
+
+def test_gridrefinement():
+    """Test GradRefinement is working as expected."""
+
+    # generate override structures for z-axis
+    center = [None, None, 0]
+    grid_size_in_vaccum = 1
+    structure = GRID_REFINEMENT.override_structure(center, grid_size_in_vaccum)
+    assert not structure.shadow
+    for axis in range(2):
+        assert structure.dl[axis] is None
+        assert structure.geometry.size[axis] == td.inf
+    dl = grid_size_in_vaccum / GRID_REFINEMENT._refinement_factor
+    assert np.isclose(structure.dl[2], dl)
+    assert np.isclose(structure.geometry.size[2], dl * GRID_REFINEMENT.num_cells)
+
+    # explicitly define step size in refinement region that is smaller than that of refinement_factor
+    dl = 1
+    grid_refinement = GRID_REFINEMENT.updated_copy(dl=dl)
+    structure = grid_refinement.override_structure(center, grid_size_in_vaccum)
+    for axis in range(2):
+        assert structure.dl[axis] is None
+        assert structure.geometry.size[axis] == td.inf
+    assert np.isclose(structure.dl[2], dl)
+    assert np.isclose(structure.geometry.size[2], dl * GRID_REFINEMENT.num_cells)
+
+
+def test_layerrefinement():
+    """Test LayerRefinementSpec is working as expected."""
+
+    # size along axis must be inf
+    with pytest.raises(pydantic.ValidationError):
+        _ = LayerRefinementSpec(axis=0, size=(td.inf, 0, 0))
+
+    # classmethod
+    for axis in range(3):
+        layer = LayerRefinementSpec.from_layer_bounds(axis=axis, bounds=(0, 1))
+        assert layer.center[axis] == 0.5
+        assert layer.size[axis] == 1
+        assert layer.size[(axis + 1) % 3] == td.inf
+        assert layer.size[(axis + 2) % 3] == td.inf
+        assert layer._is_inplane_unbounded
+
+    layer = LayerRefinementSpec.from_bounds(axis=axis, rmin=(0, 0, 0), rmax=(1, 2, 3))
+    layer = LayerRefinementSpec.from_bounds(rmin=(0, 0, 0), rmax=(1, 2, 3))
+    assert layer.axis == 0
+    assert np.isclose(layer.length_axis, 1)
+    assert np.isclose(layer.center_axis, 0.5)
+    assert not layer._is_inplane_unbounded
+
+    # from structures
+    structures = [td.Structure(geometry=td.Box(size=(td.inf, 2, 3)), medium=td.Medium())]
+    layer = LayerRefinementSpec.from_structures(structures)
+    assert layer.axis == 1
+
+    with pytest.raises(pydantic.ValidationError):
+        structures = [
+            td.Structure(geometry=td.Box(size=(td.inf, td.inf, td.inf)), medium=td.Medium())
+        ]
+        layer = LayerRefinementSpec.from_structures(structures)
+
+    with pytest.raises(pydantic.ValidationError):
+        _ = LayerRefinementSpec.from_layer_bounds(axis=axis, bounds=(0, td.inf))
+    with pytest.raises(pydantic.ValidationError):
+        _ = LayerRefinementSpec.from_layer_bounds(axis=axis, bounds=(td.inf, 0))
+    with pytest.raises(pydantic.ValidationError):
+        _ = LayerRefinementSpec.from_layer_bounds(axis=axis, bounds=(-td.inf, 0))
+    with pytest.raises(pydantic.ValidationError):
+        _ = LayerRefinementSpec.from_layer_bounds(axis=axis, bounds=(1, -1))
+
+
+def test_layerrefinement_inplane_inside():
+    # inplane inside
+    layer = LayerRefinementSpec.from_layer_bounds(axis=2, bounds=(0, 1))
+    assert layer._inplane_inside([3e3, 4e4])
+    layer = LayerRefinementSpec(axis=1, size=(1, 0, 1))
+    assert layer._inplane_inside([0, 0])
+    assert not layer._inplane_inside([2, 0])
+
+
+def test_layerrefinement_snapping_points():
+    """Test snapping points for LayerRefinementSpec is working as expected."""
+
+    # snapping points for layer bounds
+    points = LAYER2D_REFINEMENT._snapping_points_along_axis
+    assert len(points) == 1
+    assert points[0] == (None, None, 0)
+
+    points = LAYER_REFINEMENT._snapping_points_along_axis
+    assert len(points) == 1
+    assert points[0] == (None, None, -1)
+
+
+def test_grid_spec_with_layers():
+    """Test the application of layer_specs to GridSpec."""
+
+    thickness = 1e-3
+    # 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 pin
+    pin_str = td.Structure(
+        geometry=td.Cylinder(axis=0, radius=0.1, length=1.1, center=(-0.5, 0, 0)), medium=td.PEC
+    )
+    layer = LayerRefinementSpec.from_structures(
+        [
+            pec_str,
+        ]
+    )
+    assert layer.axis == 0
+
+    sim = td.Simulation(
+        size=(4, 4, 4),
+        grid_spec=td.GridSpec.auto(
+            min_steps_per_wvl=11, wavelength=1, layer_refinement_specs=[layer]
+        ),
+        boundary_spec=td.BoundarySpec.pml(),
+        structures=[pec_str, pin_str],
+        run_time=1e-12,
+    )
+    # lower bound is snapped
+    assert any(np.isclose(sim.grid.boundaries.x, -thickness / 2))
+    # corner snapped
+    assert any(np.isclose(sim.grid.boundaries.y, -0.5))
+    assert any(np.isclose(sim.grid.boundaries.z, -0.5))
+    assert any(np.isclose(sim.grid.boundaries.z, 0.5))
+
+    # differnt laye parameters
+    def update_sim_with_newlayer(layer):
+        return sim.updated_copy(
+            grid_spec=td.GridSpec.auto(
+                min_steps_per_wvl=11, wavelength=1, layer_refinement_specs=[layer]
+            )
+        )
+
+    # bounds snapping
+    layer = LayerRefinementSpec.from_structures(
+        [
+            pec_str,
+        ],
+        bounds_snapping="bounds",
+    )
+    sim2 = update_sim_with_newlayer(layer)
+    assert any(np.isclose(sim2.grid.boundaries.x, -thickness / 2))
+    assert any(np.isclose(sim2.grid.boundaries.x, thickness / 2))
+
+    # layer thickness refinement
+    def count_grids_within_layer(sim_t):
+        float_relax = 1.001
+        x = sim_t.grid.boundaries.x
+        x = x[x >= -thickness / 2 * float_relax]
+        x = x[x <= thickness / 2 * float_relax]
+        return len(x)
+
+    layer = LayerRefinementSpec.from_structures(
+        [
+            pec_str,
+        ],
+        min_steps_along_axis=3,
+    )
+    sim2 = update_sim_with_newlayer(layer)
+    assert count_grids_within_layer(sim2) == 4
+
+    # layer thickness refinement + bounds refinement, but the latter is too coarse so that it's abandoned
+    layer = LayerRefinementSpec.from_structures(
+        [
+            pec_str,
+        ],
+        min_steps_along_axis=3,
+        bounds_refinement=td.GridRefinement(),
+    )
+    sim2 = update_sim_with_newlayer(layer)
+    assert count_grids_within_layer(sim2) == 4
+    # much finer refinement to be included
+    layer = LayerRefinementSpec.from_structures(
+        [
+            pec_str,
+        ],
+        min_steps_along_axis=3,
+        bounds_refinement=td.GridRefinement(dl=thickness / 10),
+    )
+    sim2 = update_sim_with_newlayer(layer)
+    assert count_grids_within_layer(sim2) > 10
+
+    # layer bounds refinement: combined into one structure when they overlap
+    layer = LayerRefinementSpec.from_structures(
+        [
+            pec_str,
+        ],
+        bounds_refinement=td.GridRefinement(dl=thickness * 1.1),
+        corner_finder=None,
+    )
+    sim2 = update_sim_with_newlayer(layer)
+    assert len(sim2.grid_spec.all_override_structures(list(sim2.structures), 1.0, sim2.size)) == 1
+
+    # separate when they don't overlap
+    layer = LayerRefinementSpec.from_structures(
+        [
+            pec_str,
+        ],
+        bounds_refinement=td.GridRefinement(dl=thickness * 0.9),
+        corner_finder=None,
+    )
+    sim2 = update_sim_with_newlayer(layer)
+    assert len(sim2.grid_spec.all_override_structures(list(sim2.structures), 1.0, sim2.size)) == 2

tidy3d/__init__.py

@@ -188,12 +188,15 @@
 from .components.geometry.mesh import TriangleMesh
 from .components.geometry.polyslab import PolySlab
 from .components.geometry.primitives import Cylinder, Sphere
+from .components.grid.corner_finder import CornerFinderSpec
 from .components.grid.grid import Coords, Coords1D, FieldGrid, Grid, YeeGrid
 from .components.grid.grid_spec import (
     AutoGrid,
     CustomGrid,
     CustomGridBoundaries,
+    GridRefinement,
     GridSpec,
+    LayerRefinementSpec,
     QuasiUniformGrid,
     UniformGrid,
 )
@@ -374,6 +377,9 @@ def set_logging_level(level: str) -> None:
     "CustomGrid",
     "AutoGrid",
     "CustomGridBoundaries",
+    "LayerRefinementSpec",
+    "GridRefinement",
+    "CornerFinderSpec",
     "Box",
     "Sphere",
     "Cylinder",