From 3a9d515c5e9a38fa68cea77d37c3b930cddfb74b Mon Sep 17 00:00:00 2001 From: dbochkov-flexcompute Date: Tue, 13 May 2025 16:47:29 -0700 Subject: [PATCH] automatic iterative gap meshing --- CHANGELOG.md | 1 + tests/test_components/test_layerrefinement.py | 400 ++++++++++- tidy3d/components/grid/corner_finder.py | 73 +- tidy3d/components/grid/grid_spec.py | 680 +++++++++++++++++- tidy3d/components/simulation.py | 59 +- 5 files changed, 1185 insertions(+), 28 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a5ebc5a0e5..61652088ca 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - 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. +- `LayerRefinementSpec` now supports automatic thin gap meshing through fields `gap_meshing_iters` and `dl_min_from_gap_width`. ### 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. diff --git a/tests/test_components/test_layerrefinement.py b/tests/test_components/test_layerrefinement.py index dadc92602c..b0c6487e12 100644 --- a/tests/test_components/test_layerrefinement.py +++ b/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 @@ -402,3 +405,396 @@ def test_dl_min_from_smallest_feature(): ) _ = sim.grid + + +def test_gap_meshing(): + w = 1 + length = 10 + + l_shape_1 = td.Structure( + medium=td.PECMedium(), + geometry=td.PolySlab( + axis=2, + slab_bounds=[0, 2], + vertices=[ + [0, 0], + [length, 0], + [length, w], + [w, w], + [w, length], + [0, length], + ], + ), + ) + + 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 + length, w + gap], + [w + gap + length, w + gap + w], + [w + gap + w, w + gap + w], + [w + gap + w, w + gap + length], + [w + gap, w + gap + length], + [0, w + gap + length], + [0, gap + length], + [w + gap, gap + length], + ], + ), + ) + + # 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 * length, 1.2 * length, 2), + center=(0.5 * length, 0.5 * length, 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 + + # test ingored small feature + sim_size = (1, 1, 1) + box = td.Structure( + geometry=td.Box(size=(0.95, 0.2, 0.95)), + medium=td.PECMedium(), + ) + + reentry_gap = td.Structure( + geometry=td.Box(size=(0.3, 0.4, 0.3)) + .rotated(axis=1, angle=np.pi / 4) + .translated(x=0, y=0, z=0.5), + medium=td.Medium(), + ) + + aux = td.Structure( + geometry=td.Box(center=(0.0, 0, 0.5), size=(2, 0.4, 0.23)), + medium=td.Medium(), + ) + + sim = td.Simulation( + size=sim_size, + boundary_spec=td.BoundarySpec( + x=td.Boundary.periodic(), + y=td.Boundary.periodic(), + z=td.Boundary.periodic(), + ), + structures=[box, reentry_gap, aux], + grid_spec=td.GridSpec.auto( + layer_refinement_specs=[ + td.LayerRefinementSpec( + axis=1, + size=(td.inf, 0.2, td.inf), + corner_finder=None, + gap_meshing_iters=1, + dl_min_from_gap_width=True, + ) + ], + min_steps_per_wvl=10, + wavelength=1, + ), + run_time=1e-15, + ) + assert not any( + np.isclose(sim.grid.boundaries.x, reentry_gap.geometry.bounding_box.center[0], rtol=1e-2) + ) + # _, ax = plt.subplots(1, 1, figsize=(10, 10)) + # sim.plot(y=0, ax=ax) + # sim.plot_grid(y=0, ax=ax) + # plt.show() + + # test internal polygon + + gap_z = td.Structure( + geometry=td.Box(center=(0.3, 0, 0.1), size=(0.05, 0.4, 0.5)), + medium=td.Medium(), + ) + + gap_x = td.Structure( + geometry=td.Box(center=(0.03, 0, 0.07), size=(0.3, 0.4, 0.05)), + medium=td.Medium(), + ) + + sim = td.Simulation( + size=sim_size, + boundary_spec=td.BoundarySpec( + x=td.Boundary.periodic(), + y=td.Boundary.periodic(), + z=td.Boundary.periodic(), + ), + structures=[box, gap_x, gap_z], + grid_spec=td.GridSpec.auto( + layer_refinement_specs=[ + td.LayerRefinementSpec( + axis=1, + size=(td.inf, 0.2, td.inf), + corner_finder=None, + gap_meshing_iters=2, + dl_min_from_gap_width=True, + ) + ], + min_steps_per_wvl=10, + wavelength=1, + ), + run_time=1e-15, + ) + assert any(np.isclose(sim.grid.boundaries.x, gap_z.geometry.center[0], atol=1e-4)) + assert any(np.isclose(sim.grid.boundaries.z, gap_x.geometry.center[2], atol=1e-4)) + # _, ax = plt.subplots(1, 1, figsize=(10, 10)) + # sim.plot(y=0, ax=ax) + # sim.plot_grid(y=0, ax=ax) + # plt.show() + + # test gaps near pec/pmc + sim = td.Simulation( + size=sim_size, + boundary_spec=td.BoundarySpec( + x=td.Boundary.pec(), + y=td.Boundary.pml(), + z=td.Boundary.pmc(), + ), + structures=[box], + grid_spec=td.GridSpec.auto( + layer_refinement_specs=[ + td.LayerRefinementSpec( + axis=1, + size=(td.inf, 0.2, td.inf), + corner_finder=None, + gap_meshing_iters=1, + dl_min_from_gap_width=True, + ) + ], + min_steps_per_wvl=10, + wavelength=1, + ), + run_time=1e-15, + ) + + expected_grid_line = box.geometry.size[0] / 2 + (sim_size[0] - box.geometry.size[0]) / 4 + assert any(np.isclose(sim.grid.boundaries.x, expected_grid_line)) + assert any(np.isclose(sim.grid.boundaries.x, -expected_grid_line)) + + expected_grid_line = box.geometry.size[2] / 2 + (sim_size[2] - box.geometry.size[2]) / 4 + assert any(np.isclose(sim.grid.boundaries.z, expected_grid_line)) + assert any(np.isclose(sim.grid.boundaries.z, -expected_grid_line)) + + # _, ax = plt.subplots(1, 1, figsize=(10, 10)) + # sim.plot(y=0, ax=ax) + # sim.plot_grid(y=0, ax=ax) + # plt.show() + + # test limited size of layer spec + sim = td.Simulation( + size=sim_size, + boundary_spec=td.BoundarySpec( + x=td.Boundary.pec(), + y=td.Boundary.pml(), + z=td.Boundary.pmc(), + ), + structures=[box], + grid_spec=td.GridSpec.auto( + layer_refinement_specs=[ + td.LayerRefinementSpec( + axis=1, + size=(0.5, 0.2, 0.5), + center=(0.5, 0, -0.5), + corner_finder=None, + gap_meshing_iters=1, + dl_min_from_gap_width=True, + ) + ], + min_steps_per_wvl=10, + wavelength=1, + ), + run_time=1e-15, + ) + + expected_grid_line = box.geometry.size[0] / 2 + (sim_size[0] - box.geometry.size[0]) / 4 + assert any(np.isclose(sim.grid.boundaries.x, expected_grid_line)) + assert not any(np.isclose(sim.grid.boundaries.x, -expected_grid_line, atol=1e-2)) + + expected_grid_line = box.geometry.size[2] / 2 + (sim_size[2] - box.geometry.size[2]) / 4 + assert not any(np.isclose(sim.grid.boundaries.z, expected_grid_line, atol=1e-2)) + assert any(np.isclose(sim.grid.boundaries.z, -expected_grid_line)) + + # pretty much zero size layer spec + sim = td.Simulation( + size=sim_size, + boundary_spec=td.BoundarySpec( + x=td.Boundary.pec(), + y=td.Boundary.pml(), + z=td.Boundary.pmc(), + ), + structures=[box], + grid_spec=td.GridSpec.auto( + layer_refinement_specs=[ + td.LayerRefinementSpec( + axis=1, + size=(0.05, 0.2, 0.05), + center=(0.05, 0, 0.05), + corner_finder=None, + gap_meshing_iters=2, + dl_min_from_gap_width=True, + ) + ], + min_steps_per_wvl=10, + wavelength=1, + ), + run_time=1e-15, + ) + + expected_grid_line = box.geometry.size[0] / 2 + (sim_size[0] - box.geometry.size[0]) / 4 + assert not any(np.isclose(sim.grid.boundaries.x, expected_grid_line)) + assert not any(np.isclose(sim.grid.boundaries.x, -expected_grid_line, atol=1e-2)) + + expected_grid_line = box.geometry.size[2] / 2 + (sim_size[2] - box.geometry.size[2]) / 4 + assert not any(np.isclose(sim.grid.boundaries.z, expected_grid_line, atol=1e-2)) + assert not any(np.isclose(sim.grid.boundaries.z, -expected_grid_line)) + + # layer spec is outside + sim = td.Simulation( + size=sim_size, + boundary_spec=td.BoundarySpec( + x=td.Boundary.pec(), + y=td.Boundary.pml(), + z=td.Boundary.pmc(), + ), + structures=[box], + grid_spec=td.GridSpec.auto( + layer_refinement_specs=[ + td.LayerRefinementSpec( + axis=1, + size=(0.01, 0.2, 0.01), + center=(10.0, 0, 10.0), + corner_finder=None, + gap_meshing_iters=1, + dl_min_from_gap_width=True, + ) + ], + min_steps_per_wvl=10, + wavelength=1, + ), + run_time=1e-15, + ) + + expected_grid_line = box.geometry.size[0] / 2 + (sim_size[0] - box.geometry.size[0]) / 4 + assert not any(np.isclose(sim.grid.boundaries.x, expected_grid_line)) + assert not any(np.isclose(sim.grid.boundaries.x, -expected_grid_line, atol=1e-2)) + + expected_grid_line = box.geometry.size[2] / 2 + (sim_size[2] - box.geometry.size[2]) / 4 + assert not any(np.isclose(sim.grid.boundaries.z, expected_grid_line, atol=1e-2)) + assert not any(np.isclose(sim.grid.boundaries.z, -expected_grid_line)) + + # _, ax = plt.subplots(1, 1, figsize=(10, 10)) + # sim.plot(y=0, ax=ax) + # sim.plot_grid(y=0, ax=ax) + # plt.show() + + # test gap near periodic + sim = td.Simulation( + size=sim_size, + center=(0.025, -0.25, 0), + boundary_spec=td.BoundarySpec( + x=td.Boundary.periodic(), + y=td.Boundary.pec(), + z=td.Boundary.periodic(), + ), + structures=[box], + grid_spec=td.GridSpec.auto( + layer_refinement_specs=[ + td.LayerRefinementSpec( + axis=1, + size=(td.inf, 0.2, td.inf), + corner_finder=None, + gap_meshing_iters=1, + dl_min_from_gap_width=True, + ) + ], + min_steps_per_wvl=10, + wavelength=1, + ), + run_time=1e-15, + ) + + # _, ax = plt.subplots(1, 1, figsize=(10, 10)) + # sim.plot(y=0, ax=ax) + # sim.plot_grid(y=0, ax=ax) + # plt.show() + + assert any(np.isclose(sim.grid.boundaries.x, 0.5, atol=1e-4)) + assert any(np.isclose(sim.grid.boundaries.z, -0.5, atol=1e-4)) + + # test a thin strip near periodic + strip = td.Structure( + geometry=td.Box(center=(0, 0.5, 0), size=(0.2, 0.05, 0.6)), + medium=td.PECMedium(), + ) + + sim = td.Simulation( + size=sim_size, + boundary_spec=td.BoundarySpec( + x=td.Boundary.pec(), + y=td.Boundary.periodic(), + z=td.Boundary.pmc(), + ), + structures=[strip], + grid_spec=td.GridSpec.auto( + layer_refinement_specs=[ + td.LayerRefinementSpec( + axis=0, + size=(0.2, td.inf, td.inf), + corner_finder=None, + gap_meshing_iters=1, + dl_min_from_gap_width=True, + ) + ], + min_steps_per_wvl=10, + wavelength=1, + ), + run_time=1e-15, + ) + + assert any(np.isclose(sim.grid.boundaries.y, strip.geometry.center[1], atol=1e-4)) + # _, ax = plt.subplots(1, 1, figsize=(10, 10)) + # sim.plot(x=0, ax=ax) + # sim.plot_grid(x=0, ax=ax) + # plt.show() diff --git a/tidy3d/components/grid/corner_finder.py b/tidy3d/components/grid/corner_finder.py index 19aa0369cd..aa3c5347bf 100644 --- a/tidy3d/components/grid/corner_finder.py +++ b/tidy3d/components/grid/corner_finder.py @@ -1,6 +1,6 @@ """Find corners of structures on a 2D plane.""" -from typing import List, Literal, Optional, Tuple +from typing import Any, List, Literal, Optional, Tuple import numpy as np import pydantic.v1 as pd @@ -11,7 +11,7 @@ from ..geometry.utils import merging_geometries_on_plane from ..medium import PEC, LossyMetalMedium from ..structure import Structure -from ..types import ArrayFloat1D, ArrayFloat2D, Axis +from ..types import ArrayFloat1D, ArrayFloat2D, Axis, Shapely CORNER_ANGLE_THRESOLD = 0.1 * np.pi @@ -75,16 +75,16 @@ def _no_min_dl_override(self): ) ) - def _corners_and_convexity( - self, + @classmethod + def _merged_pec_on_plane( + cls, normal_axis: Axis, coord: float, structure_list: List[Structure], - 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`. - + center: Tuple[float, float] = [0, 0, 0], + size: Tuple[float, float, float] = [inf, inf, inf], + ) -> List[Tuple[Any, Shapely]]: + """On a 2D plane specified by axis = `normal_axis` and coordinate `coord`, merge geometries made of PEC. Parameters ---------- @@ -94,21 +94,23 @@ def _corners_and_convexity( 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. + center : Tuple[float, float] = [0, 0, 0] + 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) Returns ------- - ArrayFloat2D - Corner coordinates. + List[Tuple[Any, Shapely]] + List of shapes and their property value on the plane after merging. """ # Construct plane - center = [0, 0, 0] - size = [inf, inf, inf] - center[normal_axis] = coord - size[normal_axis] = 0 - plane = Box(center=center, size=size) + slice_center = list(center) + slice_size = list(size) + slice_center[normal_axis] = coord + slice_size[normal_axis] = 0 + plane = Box(center=slice_center, size=slice_size) # prepare geometry and medium list geometry_list = [structure.geometry for structure in structure_list] @@ -121,6 +123,41 @@ def _corners_and_convexity( # merge geometries merged_geos = merging_geometries_on_plane(geometry_list, plane, medium_list) + return merged_geos + + def _corners_and_convexity( + self, + normal_axis: Axis, + coord: float, + structure_list: List[Structure], + 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 PEC. + + + 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. + ravel : bool + Whether to put the resulting corners in a single list or per polygon. + + Returns + ------- + Tuple[ArrayFloat2D, ArrayFloat1D] + Corner coordinates and their convexity. + """ + + # merge geometries + merged_geos = self._merged_pec_on_plane( + normal_axis=normal_axis, coord=coord, structure_list=structure_list + ) + # corner finder corner_list = [] convexity_list = [] diff --git a/tidy3d/components/grid/grid_spec.py b/tidy3d/components/grid/grid_spec.py index 702e98af9d..17f41e088c 100644 --- a/tidy3d/components/grid/grid_spec.py +++ b/tidy3d/components/grid/grid_spec.py @@ -8,11 +8,11 @@ import numpy as np import pydantic.v1 as pd -from ...constants import C_0, MICROMETER, inf +from ...constants import C_0, MICROMETER, dp_eps, inf from ...exceptions import SetupError from ...log import log from ..base import Tidy3dBaseModel, cached_property, skip_if_fields_missing -from ..geometry.base import Box +from ..geometry.base import Box, ClipOperation from ..lumped_element import LumpedElementType from ..source.utils import SourceType from ..structure import MeshOverrideStructure, Structure, StructureType @@ -36,6 +36,9 @@ # Default refinement factor in GridRefinement when both dl and refinement_factor are not defined DEFAULT_REFINEMENT_FACTOR = 2 +# Tolerance for distinguishing pec/grid intersections +GAP_MESHING_TOL = 1e-3 + class GridSpec1d(Tidy3dBaseModel, ABC): """Abstract base class, defines 1D grid generation specifications.""" @@ -62,6 +65,9 @@ def make_coords( symmetry : Tuple[Symmetry, Symmetry, Symmetry] Reflection symmetry across a plane bisecting the simulation domain normal to each of the three axes. + periodic : bool + Apply periodic boundary condition or not. + Only relevant for autogrids. wavelength : float Free-space wavelength. num_pml_layers : Tuple[int, int] @@ -1037,6 +1043,20 @@ class LayerRefinementSpec(Box): "and the projection of the simulation domain overlaps.", ) + gap_meshing_iters: pd.NonNegativeInt = pd.Field( + 1, + title="Gap Meshing Iterations", + description="Number of recursive iterations for resolving thin gaps. " + "The underlying algorithm detects gaps contained in a single cell and places a snapping plane at the gaps's centers.", + ) + + dl_min_from_gap_width: bool = pd.Field( + True, + title="Set ``dl_min`` from Estimated Gap Width", + description="Take into account autodetected minimal PEC gap width when determining ``dl_min``. " + "This only applies if ``dl_min`` in ``AutoGrid`` specification is not set.", + ) + @pd.validator("axis", always=True) @skip_if_fields_missing(["size"]) def _finite_size_along_axis(cls, val, values): @@ -1057,6 +1077,8 @@ def from_layer_bounds( corner_snapping: bool = True, corner_refinement: GridRefinement = GridRefinement(), refinement_inside_sim_only: bool = True, + gap_meshing_iters: pd.NonNegativeInt = 1, + dl_min_from_gap_width: bool = True, ): """Constructs a :class:`LayerRefiementSpec` that is unbounded in inplane dimensions from bounds along layer thickness dimension. @@ -1082,6 +1104,10 @@ def from_layer_bounds( Inplane mesh refinement factor around corners. refinement_inside_sim_only : bool = True Apply refinement only to features inside simulation domain. + gap_meshing_iters : bool = True + Number of recursive iterations for resolving thin gaps. + dl_min_from_gap_width : bool = True + Take into account autodetected minimal PEC gap width when determining ``dl_min``. Example @@ -1103,6 +1129,8 @@ def from_layer_bounds( corner_snapping=corner_snapping, corner_refinement=corner_refinement, refinement_inside_sim_only=refinement_inside_sim_only, + gap_meshing_iters=gap_meshing_iters, + dl_min_from_gap_width=dl_min_from_gap_width, ) @classmethod @@ -1118,6 +1146,8 @@ def from_bounds( corner_snapping: bool = True, corner_refinement: GridRefinement = GridRefinement(), refinement_inside_sim_only: bool = True, + gap_meshing_iters: pd.NonNegativeInt = 1, + dl_min_from_gap_width: bool = True, ): """Constructs a :class:`LayerRefiementSpec` from minimum and maximum coordinate bounds. @@ -1145,6 +1175,10 @@ def from_bounds( Inplane mesh refinement factor around corners. refinement_inside_sim_only : bool = True Apply refinement only to features inside simulation domain. + gap_meshing_iters : bool = True + Number of recursive iterations for resolving thin gaps. + dl_min_from_gap_width : bool = True + Take into account autodetected minimal PEC gap width when determining ``dl_min``. Example @@ -1166,6 +1200,8 @@ def from_bounds( corner_snapping=corner_snapping, corner_refinement=corner_refinement, refinement_inside_sim_only=refinement_inside_sim_only, + gap_meshing_iters=gap_meshing_iters, + dl_min_from_gap_width=dl_min_from_gap_width, ) @classmethod @@ -1180,6 +1216,8 @@ def from_structures( corner_snapping: bool = True, corner_refinement: GridRefinement = GridRefinement(), refinement_inside_sim_only: bool = True, + gap_meshing_iters: pd.NonNegativeInt = 1, + dl_min_from_gap_width: bool = True, ): """Constructs a :class:`LayerRefiementSpec` from the bounding box of a list of structures. @@ -1205,6 +1243,10 @@ def from_structures( Inplane mesh refinement factor around corners. refinement_inside_sim_only : bool = True Apply refinement only to features inside simulation domain. + gap_meshing_iters : bool = True + Number of recursive iterations for resolving thin gaps. + dl_min_from_gap_width : bool = True + Take into account autodetected minimal PEC gap width when determining ``dl_min``. """ @@ -1226,6 +1268,8 @@ def from_structures( corner_snapping=corner_snapping, corner_refinement=corner_refinement, refinement_inside_sim_only=refinement_inside_sim_only, + gap_meshing_iters=gap_meshing_iters, + dl_min_from_gap_width=dl_min_from_gap_width, ) @cached_property @@ -1510,6 +1554,449 @@ 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, boundary + ) -> 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) + cells_ij = [] + # relative displacements of intersection from the bottom of the cell along y axis + cells_dy = [] + + # for each polygon vertex find the index of the first grid line on the right + grid_lines_on_right = np.argmax(grid_x_coords[:, None] >= poly_vertices[None, :, 0], axis=0) + grid_lines_on_right[poly_vertices[:, 0] >= grid_x_coords[-1]] = len(grid_x_coords) + # once we know these indices then we can find grid lines intersected by the i-th + # segment of the polygon as + # [grid_lines_on_right[i], grid_lines_on_right[i+1]) for grid_lines_on_right[i] > grid_lines_on_right[i+1] + # or + # [grid_lines_on_right[i+1], grid_lines_on_right[i]) for grid_lines_on_right[i] < grid_lines_on_right[i+1] + + # loop over segments of the polygon and determine in which cells and where exactly they cross grid lines + # v_beg and v_end are the starting and ending points of the segment + # ind_beg and ind_end are starting and ending indices of vertical grid lines that the segment intersects + # as described above + for ind_beg, ind_end, v_beg, v_end in zip( + grid_lines_on_right, + np.roll(grid_lines_on_right, -1), + poly_vertices, + np.roll(poly_vertices, axis=0, shift=-1), + ): + # no intersections + if ind_end == ind_beg: + continue + + # sort vertices in ascending order to make treatmeant unifrom + reverse = False + if ind_beg > ind_end: + reverse = True + ind_beg, ind_end, v_beg, v_end = ind_end, ind_beg, v_end, v_beg + + # x coordinates are simply x coordinates of intersected vertical grid lines + intersections_x = grid_x_coords[ind_beg:ind_end] + + # y coordinates can be found from line equation + intersections_y = v_beg[1] + (v_end[1] - v_beg[1]) / (v_end[0] - v_beg[0]) * ( + intersections_x - v_beg[0] + ) + + # 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 >= grid_y_coords[0], intersections_y <= grid_y_coords[-1] + ) + + intersections_y = intersections_y[inds_inside_grid] + + # find i and j indices of cells which contain these intersections + + # i indices are simply indices of crossed vertical grid lines + cell_is = np.arange(ind_beg, ind_end)[inds_inside_grid] + + # j indices can be computed by finding insertion indices + # of y coordinates of intersection points into array of y coordinates + # of the grid lines that preserve sorting + cell_js = np.searchsorted(grid_y_coords, intersections_y) - 1 + + # 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] + ) + + # preserve uniform ordering along perimeter of the polygon + if reverse: + cell_is = cell_is[::-1] + cell_js = cell_js[::-1] + dy = dy[::-1] + + # 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) + 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 = 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] + + # 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] + + # if we go beyond simulation domain boundary, either ignore + # or wrap periodically depending on boundary conditions + cells_ij_zero_side = points_to_duplicate_near_zero - np.array([0, 1]) + cells_zero_side_out = cells_ij_zero_side[:, 1] == -1 + if boundary[0] == "periodic": + cells_ij_zero_side[cells_zero_side_out, 1] = len(grid_y_coords) - 2 + else: + cells_ij_zero_side = cells_ij_zero_side[cells_zero_side_out == 0] + + cells_ij_one_side = points_to_duplicate_near_one + np.array([0, 1]) + cells_one_side_out = cells_ij_one_side[:, 1] == len(grid_y_coords) - 1 + if boundary[1] == "periodic": + cells_ij_one_side[cells_one_side_out, 1] = 0 + else: + cells_ij_one_side = cells_ij_one_side[cells_one_side_out == 0] + + cells_ij = np.concatenate( + [ + cells_ij, + cells_ij_zero_side, + cells_ij_one_side, + ] + ) + cells_dy = np.concatenate( + [ + cells_dy, + np.ones(len(cells_ij_zero_side)), + np.zeros(len(cells_ij_one_side)), + ] + ) + + return cells_ij, cells_dy + + def _process_poly( + self, grid_x_coords, grid_y_coords, poly_vertices, boundaries + ) -> 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, boundaries[1] + ) + + # 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), boundaries[0] + ) + 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, boundaries + ) -> Tuple[List[Tuple[int, int]], List[float], List[Tuple[int, int]], List[float]]: + """Detect intersection points of geometries boundaries and grid lines.""" + + # cells that contain intersections of vertical grid lines + v_cells_ij = [] + # relative locations of those intersections (along y axis) + v_cells_dy = [] + + # cells that contain intersections of horizontal grid lines + h_cells_ij = [] + # relative locations of those intersections (along x axis) + h_cells_dx = [] + + # for PEC and PMC boundary - treat them as PEC structure + # so that gaps are resolved near boundaries if any + nx = len(x) + ny = len(y) + + if boundaries[0][0] == "pec/pmc": + h_cells_ij.append(np.transpose([np.zeros(ny), np.arange(ny)]).astype(int)) + h_cells_dx.append(np.zeros(ny)) + + if boundaries[0][1] == "pec/pmc": + h_cells_ij.append(np.transpose([(nx - 2) * np.ones(ny), np.arange(ny)]).astype(int)) + h_cells_dx.append(np.ones(ny)) + + if boundaries[1][0] == "pec/pmc": + v_cells_ij.append(np.transpose([np.arange(nx), np.zeros(nx)]).astype(int)) + v_cells_dy.append(np.zeros(nx, dtype=int)) + + if boundaries[1][1] == "pec/pmc": + v_cells_ij.append(np.transpose([np.arange(nx), (ny - 2) * np.ones(nx)]).astype(int)) + v_cells_dy.append(np.ones(nx)) + + # 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: + poly = poly.normalize().buffer(0) + + # find intersections of a polygon with grid lines + # specifically: + # 0. cells that contain intersections of vertical grid lines + # 1. relative locations of those intersections along y axis + # 2. cells that contain intersections of horizontal grid lines + # 3. relative locations of those intersections along x axis + data = self._process_poly(x, y, np.array(poly.exterior.coords)[:-1], boundaries) + + 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], boundaries) + 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]) + + if len(v_cells_ij) > 0: + v_cells_ij = np.concatenate(v_cells_ij) + v_cells_dy = np.concatenate(v_cells_dy) + + if len(h_cells_ij) > 0: + h_cells_ij = np.concatenate(h_cells_ij) + 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 + ) -> 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 + resolve thin gaps and strips. + """ + min_gap_width = inf + + snapping_lines_y = [] + 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_and_disps = np.transpose([cell_linear_inds, relative_vert_disp]) + # 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 + cell_linear_inds_unique, counts = np.unique( + cell_linear_inds_and_disps_unique[:, 0], return_counts=True + ) + # when we count intersections we use linearized 2d index because we really + # need to count intersections in each cell separately + + # but when we need to decide about refinement, due to cartesian nature of grid + # we will need to consider all cells with a given j index at a time + + # so, let's compute j index for each cell in the unique list + cell_linear_inds_unique_j = cell_linear_inds_unique % len(grid_y_coords) + + # loop through all j rows that contain intersections + for ind_j in np.unique(cell_linear_inds_unique_j): + # we need to refine between two grid lines corresponding to index j + # if at least one cell with given j contains > 1 intersections + + # get all intersected cells with given j index + 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 + ] + ) + + # 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 + ) + # 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 + ] + + # 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]), + ) + + return snapping_lines_y, min_gap_width + + def _resolve_gaps( + self, structures: List[Structure], grid: Grid, boundaries: Tuple, center, size + ) -> Tuple[List[CoordinateOptional], float]: + """Detect underresolved gaps and place snapping lines in them. Also return the detected minimal gap width.""" + + # get x and y coordinates of grid lines + _, tan_dims = Box.pop_axis([0, 1, 2], self.axis) + x = grid.boundaries.to_list[tan_dims[0]] + y = grid.boundaries.to_list[tan_dims[1]] + + _, boundaries_tan = Box.pop_axis(boundaries, self.axis) + + # restrict to the size of layer spec + rmin, rmax = self.bounds + _, rmin = Box.pop_axis(rmin, self.axis) + _, rmax = Box.pop_axis(rmax, self.axis) + + new_coords = [] + new_boundaries = [] + for coord, cmin, cmax, bdry in zip([x, y], rmin, rmax, boundaries_tan): + if cmax <= coord[0] or cmin >= coord[-1]: + return [], inf + else: + if cmin < coord[0]: + ind_min = 0 + else: + ind_min = max(0, np.argmax(coord >= cmin) - 1) + + if cmax > coord[-1]: + ind_max = len(coord) - 1 + else: + ind_max = np.argmax(coord >= cmax) + + if ind_min >= ind_max - 1: + return [], inf + + new_coords.append(coord[ind_min : (ind_max + 1)]) + # ignore boundary conditions if we are not touching them + new_boundaries.append( + [ + None if ind_min > 0 else bdry[0], + None if ind_max < len(coord) - 1 else bdry[1], + ] + ) + + x, y = new_coords + + # restrict size of the plane where pec polygons are found in case of periodic boundary conditions + # this is to make sure gaps across periodic boundary conditions are resolved + # (if there is a PEC structure going into periodic boundary, now it will generate a grid line + # intersection next to that boundary and it will be propagated to the other side) + restricted_size_tan = [ + s * (1.0 - dp_eps) if b[0] == "periodic" else inf + for b, s in zip( + new_boundaries, + size, + ) + ] + restricted_size = Box.unpop_axis(size[self.axis], restricted_size_tan, self.axis) + + # 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, + center=center, + size=restricted_size, + ) + + # find intersections of pec polygons with grid lines + # specifically: + # 0. cells that contain intersections of vertical grid lines + # 1. relative locations of those intersections along y axis + # 2. cells that contain intersections of horizontal grid lines + # 3. relative locations of those intersections along x axis + v_cells_ij, v_cells_dy, h_cells_ij, h_cells_dx = self._process_slice( + x, y, plane_slice, new_boundaries + ) + + # generate horizontal snapping lines + 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 + 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(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 + ] + + return snapping_lines_x_3d + snapping_lines_y_3d, detected_gap_width + class GridSpec(Tidy3dBaseModel): """Collective grid specification for all three dimensions. @@ -1906,6 +2393,11 @@ def make_grid( lumped_elements: List[LumpedElementType] = (), internal_override_structures: List[MeshOverrideStructure] = None, internal_snapping_points: List[CoordinateOptional] = None, + boundary_types: Tuple[Tuple[str, str], Tuple[str, str], Tuple[str, str]] = [ + [None, None], + [None, None], + [None, None], + ], ) -> Grid: """Make the entire simulation grid based on some simulation parameters. @@ -1917,6 +2409,9 @@ def make_grid( symmetry : Tuple[Symmetry, Symmetry, Symmetry] Reflection symmetry across a plane bisecting the simulation domain normal to each of the three axes. + periodic: Tuple[bool, bool, bool] + Apply periodic boundary condition or not along each of the dimensions. + Only relevant for autogrids. sources : List[SourceType] List of sources. num_pml_layers : List[Tuple[float, float]] @@ -1927,6 +2422,186 @@ def make_grid( If `None`, recomputes internal override structures. internal_snapping_points : List[CoordinateOptional] If `None`, recomputes internal snapping points. + boundary_types : Tuple[Tuple[str, str], Tuple[str, str], Tuple[str, str]] = [[None, None], [None, None], [None, None]] + Type of boundary conditions along each dimension: "pec/pmc", "periodic", or + None for any other. This is relevant only for gap meshing. + + Returns + ------- + Grid: + Entire simulation grid. + """ + + grid, _ = self._make_grid_and_snapping_lines( + structures=structures, + symmetry=symmetry, + periodic=periodic, + sources=sources, + num_pml_layers=num_pml_layers, + lumped_elements=lumped_elements, + internal_override_structures=internal_override_structures, + internal_snapping_points=internal_snapping_points, + ) + + return grid + + def _make_grid_and_snapping_lines( + self, + structures: List[Structure], + symmetry: Tuple[Symmetry, Symmetry, Symmetry], + periodic: Tuple[bool, bool, bool], + sources: List[SourceType], + num_pml_layers: List[Tuple[pd.NonNegativeInt, pd.NonNegativeInt]], + lumped_elements: List[LumpedElementType] = (), + internal_override_structures: List[MeshOverrideStructure] = None, + internal_snapping_points: List[CoordinateOptional] = None, + boundary_types: Tuple[Tuple[str, str], Tuple[str, str], Tuple[str, str]] = [ + [None, None], + [None, None], + [None, None], + ], + ) -> Tuple[Grid, List[CoordinateOptional]]: + """Make the entire simulation grid based on some simulation parameters. + Also return snappiung point resulted from iterative gap meshing. + + Parameters + ---------- + structures : List[Structure] + List of structures present in the simulation. The first structure must be the + simulation geometry with the simulation background medium. + symmetry : Tuple[Symmetry, Symmetry, Symmetry] + Reflection symmetry across a plane bisecting the simulation domain + normal to each of the three axes. + periodic: Tuple[bool, bool, bool] + Apply periodic boundary condition or not along each of the dimensions. + Only relevant for autogrids. + sources : List[SourceType] + List of sources. + num_pml_layers : List[Tuple[float, float]] + List containing the number of absorber layers in - and + boundaries. + lumped_elements : List[LumpedElementType] + List of lumped elements. + internal_override_structures : List[MeshOverrideStructure] + If `None`, recomputes internal override structures. + internal_snapping_points : List[CoordinateOptional] + If `None`, recomputes internal snapping points. + boundary_types : Tuple[Tuple[str, str], Tuple[str, str], Tuple[str, str]] = [[None, None], [None, None], [None, None]] + Type of boundary conditions along each dimension: "pec/pmc", "periodic", or + None for any other. This is relevant only for gap meshing. + + Returns + ------- + Tuple[Grid, List[CoordinateOptional]]: + Entire simulation grid and snapping points generated during iterative gap meshing. + """ + + old_grid = self._make_grid_one_iteration( + structures=structures, + symmetry=symmetry, + periodic=periodic, + sources=sources, + num_pml_layers=num_pml_layers, + lumped_elements=lumped_elements, + internal_override_structures=internal_override_structures, + internal_snapping_points=internal_snapping_points, + ) + + sim_geometry = structures[0].geometry + + snapping_lines = [] + if len(self.layer_refinement_specs) > 0: + num_iters = max( + layer_spec.gap_meshing_iters for layer_spec in self.layer_refinement_specs + ) + + 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, + boundary_types, + center=sim_geometry.center, + size=sim_geometry.size, + ) + 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) + + if len(new_snapping_lines) == 0: + log.info( + "Grid is no longer changing. " + f"Stopping iterative gap meshing after {ind+1}/{num_iters} iterations." + ) + break + + snapping_lines = snapping_lines + new_snapping_lines + + new_grid = self._make_grid_one_iteration( + structures=structures, + symmetry=symmetry, + periodic=periodic, + sources=sources, + num_pml_layers=num_pml_layers, + lumped_elements=lumped_elements, + internal_override_structures=internal_override_structures, + internal_snapping_points=snapping_lines + internal_snapping_points, + dl_min_from_gaps=0.45 * min_gap_width, + ) + + same = old_grid == new_grid + + if same: + log.info( + "Grid is no longer changing. " + f"Stopping iterative gap meshing after {ind+1}/{num_iters} iterations." + ) + break + + old_grid = new_grid + + return old_grid, snapping_lines + + def _make_grid_one_iteration( + self, + structures: List[Structure], + symmetry: Tuple[Symmetry, Symmetry, Symmetry], + periodic: Tuple[bool, bool, bool], + sources: List[SourceType], + num_pml_layers: List[Tuple[pd.NonNegativeInt, pd.NonNegativeInt]], + lumped_elements: List[LumpedElementType] = (), + internal_override_structures: List[MeshOverrideStructure] = None, + internal_snapping_points: List[CoordinateOptional] = None, + dl_min_from_gaps: pd.PositiveFloat = inf, + ) -> Grid: + """Make the entire simulation grid based on some simulation parameters. + + Parameters + ---------- + structures : List[Structure] + List of structures present in the simulation. The first structure must be the + simulation geometry with the simulation background medium. + symmetry : Tuple[Symmetry, Symmetry, Symmetry] + Reflection symmetry across a plane bisecting the simulation domain + normal to each of the three axes. + periodic: Tuple[bool, bool, bool] + Apply periodic boundary condition or not along each of the dimensions. + Only relevant for autogrids. + sources : List[SourceType] + List of sources. + num_pml_layers : List[Tuple[float, float]] + List containing the number of absorber layers in - and + boundaries. + lumped_elements : List[LumpedElementType] + List of lumped elements. + internal_override_structures : List[MeshOverrideStructure] + If `None`, recomputes internal override structures. + internal_snapping_points : List[CoordinateOptional] + If `None`, recomputes internal snapping points. + dl_min_from_gaps : pd.PositiveFloat + Minimal grid size computed based on autodetected gaps. + Returns ------- @@ -2005,6 +2680,7 @@ def make_grid( sim_size, lumped_elements, ) + new_dl_min = min(new_dl_min, dl_min_from_gaps) for ind, grid in enumerate(grids_1d): if isinstance(grid, AutoGrid) and grid._undefined_dl_min: grids_1d[ind] = grid.updated_copy(dl_min=new_dl_min) diff --git a/tidy3d/components/simulation.py b/tidy3d/components/simulation.py index d89960434a..882aeeef4f 100644 --- a/tidy3d/components/simulation.py +++ b/tidy3d/components/simulation.py @@ -988,7 +988,7 @@ def plot_grid( edgecolor=kwargs["colors"], alpha=override_structures_alpha, ), - ] * 2 + ] * 3 plot_params[0] = plot_params[0].include_kwargs(edgecolor=kwargs["colors_internal"]) if self.grid_spec.auto_grid_used: @@ -1014,7 +1014,12 @@ def plot_grid( # Plot snapping points for points, plot_param in zip( - [self.internal_snapping_points, self.grid_spec.snapping_points], 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) @@ -1188,20 +1193,31 @@ def set_plot_params(boundary_edge, lim, side, thickness): # return plot_sim_3d(self, width=width, height=height) @cached_property - def grid(self) -> Grid: + def _grid_and_snapping_lines(self) -> Tuple[Grid, List[CoordinateOptional]]: """FDTD grid spatial locations and information. Returns ------- - :class:`.Grid` - :class:`.Grid` storing the spatial locations relevant to the simulation. + Tuple[:class:`.Grid`, List[CoordinateOptional]] + :class:`.Grid` storing the spatial locations relevant to the simulation + the list of snapping points generated during iterative gap meshing. """ # Add a simulation Box as the first structure structures = [Structure(geometry=self.geometry, medium=self.medium)] structures += self.static_structures - grid = self.grid_spec.make_grid( + # Get boundary types (relevant for gap meshing) + boundary_types = [[None] * 2] * 3 + + for dim, boundary in enumerate(self.boundary_spec.to_list): + for side, edge in enumerate(boundary): + if isinstance(edge, (PECBoundary, PMCBoundary)): + boundary_types[dim][side] = "pec/pmc" + elif isinstance(edge, (Periodic, BlochBoundary)): + boundary_types[dim][side] = "periodic" + + grid, lines = self.grid_spec._make_grid_and_snapping_lines( structures=structures, symmetry=self.symmetry, periodic=self._periodic, @@ -1210,12 +1226,43 @@ def grid(self) -> Grid: lumped_elements=self.lumped_elements, internal_snapping_points=self.internal_snapping_points, internal_override_structures=self.internal_override_structures, + boundary_types=boundary_types, ) + # This would AutoGrid the in-plane directions of the 2D materials + # return self._grid_corrections_2dmaterials(grid) + return grid, lines + + @cached_property + def grid(self) -> Grid: + """FDTD grid spatial locations and information. + + Returns + ------- + :class:`.Grid` + :class:`.Grid` storing the spatial locations relevant to the simulation. + """ + + grid, _ = self._grid_and_snapping_lines + # 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]: + """Snapping points resulted from iterative gap meshing. + + Returns + ------- + List[CoordinateOptional] + List of snapping lines resolving thin gaps and strips. + """ + + _, lines = self._grid_and_snapping_lines + + return lines + @cached_property def static_structures(self) -> list[Structure]: """Structures in simulation with all autograd tracers removed."""