Snapping points with optional coordinate

Add LayerRefinementSpec for automatic mesh refinement

Shift _filter_structures_plane to geometry/utils.py

Automatically sets dl_min if not set yet

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(

tests/test_components/test_layerrefinement.py

@@ -0,0 +1,123 @@
+"""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_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_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_gridrefinement():
+    """Test GradRefinement is working as expected."""
+
+    # no override grid step information
+    with pytest.raises(pydantic.ValidationError):
+        _ = GridRefinement(dl=None, refinement_factor=None)
+
+    # 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 = 0.01
+    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_unbounded_layer(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
+
+    layer = LayerRefinementSpec.from_bounds(axis=axis, rmin=(0, 0, 0), rmax=(1, 2, 3))
+
+    with pytest.raises(pydantic.ValidationError):
+        _ = LayerRefinementSpec.from_unbounded_layer(axis=axis, bounds=(0, td.inf))
+    with pytest.raises(pydantic.ValidationError):
+        _ = LayerRefinementSpec.from_unbounded_layer(axis=axis, bounds=(td.inf, 0))
+    with pytest.raises(pydantic.ValidationError):
+        _ = LayerRefinementSpec.from_unbounded_layer(axis=axis, bounds=(-td.inf, 0))
+    with pytest.raises(pydantic.ValidationError):
+        _ = LayerRefinementSpec.from_unbounded_layer(axis=axis, bounds=(1, -1))
+
+
+def test_layerrefinement_inplane_inside():
+    # inplane inside
+    layer = LayerRefinementSpec.from_unbounded_layer(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) == 2
+#     assert points[0] == (None, None, -1)
+#     assert points[1] == (None, None, 1)

tidy3d/__init__.py

@@ -122,12 +122,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,
 )
@@ -329,6 +332,9 @@ def set_logging_level(level: str) -> None:
     "CustomGrid",
     "AutoGrid",
     "CustomGridBoundaries",
+    "LayerRefinementSpec",
+    "GridRefinement",
+    "CornerFinderSpec",
     "Box",
     "Sphere",
     "Cylinder",

tidy3d/components/geometry/utils.py

@@ -4,13 +4,13 @@
 
 from enum import Enum
 from math import isclose
-from typing import Optional, Tuple, Union
+from typing import Any, List, Optional, Tuple, Union
 
 import numpy as np
 import pydantic as pydantic
 
 from ...constants import fp_eps
-from ...exceptions import Tidy3dError
+from ...exceptions import SetupError, Tidy3dError
 from ..base import Tidy3dBaseModel
 from ..geometry.base import Box
 from ..grid.grid import Grid
@@ -30,6 +30,81 @@
 ]
 
 
+def merging_geometries_on_plane(
+    geometries: List[GeometryType],
+    plane: Box,
+    property_list: List[Any],
+) -> List[Tuple[Any, Shapely]]:
+    """Compute list of shapes on plane. Overlaps are removed or merged depending on
+    provided property_list.
+
+    Parameters
+    ----------
+    geometries : List[GeometryType]
+        List of structures to filter on the plane.
+    plane : Box
+        Plane specification.
+    property_list : List = None
+        Property value for each structure.
+
+    Returns
+    -------
+    List[Tuple[Any, shapely]]
+        List of shapes and their property value on the plane after merging.
+    """
+
+    if len(geometries) != len(property_list):
+        raise SetupError(
+            "Number of provided property values is not equal to the number of geometries."
+        )
+
+    shapes = []
+    for geo, prop in zip(geometries, property_list):
+        # get list of Shapely shapes that intersect at the plane
+        shapes_plane = plane.intersections_with(geo)
+
+        # Append each of them and their property information to the list of shapes
+        for shape in shapes_plane:
+            shapes.append((prop, shape, shape.bounds))
+
+    background_shapes = []
+    for prop, shape, bounds in shapes:
+        minx, miny, maxx, maxy = bounds
+
+        # loop through background_shapes (note: all background are non-intersecting or merged)
+        for index, (_prop, _shape, _bounds) in enumerate(background_shapes):
+            _minx, _miny, _maxx, _maxy = _bounds
+
+            # do a bounding box check to see if any intersection to do anything about
+            if minx > _maxx or _minx > maxx or miny > _maxy or _miny > maxy:
+                continue
+
+            # look more closely to see if intersected.
+            if shape.disjoint(_shape):
+                continue
+
+            # different prop, remove intersection from background shape
+            if prop != _prop:
+                diff_shape = (_shape - shape).buffer(0).normalize()
+                # mark background shape for removal if nothing left
+                if diff_shape.is_empty or len(diff_shape.bounds) == 0:
+                    background_shapes[index] = None
+                background_shapes[index] = (_prop, diff_shape, diff_shape.bounds)
+            # same prop, unionize shapes and mark background shape for removal
+            else:
+                shape = (shape | _shape).buffer(0).normalize()
+                background_shapes[index] = None
+
+        # after doing this with all background shapes, add this shape to the background
+        background_shapes.append((prop, shape, shape.bounds))
+
+        # remove any existing background shapes that have been marked as 'None'
+        background_shapes = [b for b in background_shapes if b is not None]
+
+    # filter out any remaining None or empty shapes (shapes with area completely removed)
+    return [(prop, shape) for (prop, shape, _) in background_shapes if shape]
+
+
 def flatten_groups(
     *geometries: GeometryType,
     flatten_nonunion_type: bool = False,