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_layerrefinement.py

@@ -0,0 +1,102 @@
+"""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, bounds=(-1, 1))
+LAYER2D_REFINEMENT = LayerRefinementSpec(axis=2, bounds=(0, 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)
+    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."""
+
+    # wrong bounds order
+    with pytest.raises(pydantic.ValidationError):
+        _ = LayerRefinementSpec(axis=0, bounds=(1, 0))
+
+    # bounds must be finite
+    with pytest.raises(pydantic.ValidationError):
+        _ = LayerRefinementSpec(axis=0, bounds=(-np.inf, 0))
+    with pytest.raises(pydantic.ValidationError):
+        _ = LayerRefinementSpec(axis=0, bounds=(0, np.inf))
+    with pytest.raises(pydantic.ValidationError):
+        _ = LayerRefinementSpec(axis=0, bounds=(-np.inf, np.inf))
+
+
+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

@@ -124,12 +124,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,
     UniformGrid,
 )
 from .components.heat_charge.boundary import (
@@ -328,6 +331,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,

tidy3d/components/grid/corner_finder.py

@@ -0,0 +1,137 @@
+"""Find corners of structures on a 2D plane."""
+
+from typing import List, Literal, Optional
+
+import numpy as np
+import pydantic.v1 as pd
+
+from ...constants import inf
+from ..base import Tidy3dBaseModel
+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
+
+CORNER_ANGLE_THRESOLD = 0.1 * np.pi
+
+
+class CornerFinderSpec(Tidy3dBaseModel):
+    """Specification for corner detection on a 2D plane."""
+
+    medium: Literal["metal", "dielectric", "all"] = pd.Field(
+        "metal",
+        title="Material To Be Considered For Refinement",
+        description="Apply mesh refinement to structures made of ``medium``, "
+        "which can take value ``metal`` for PEC and lossy metal, ``dielectric`` "
+        "for non-metallic materials, and ``all`` for all materials.",
+    )
+
+    angle_threshold: float = pd.Field(
+        CORNER_ANGLE_THRESOLD,
+        title="Angle Threshold In Defining Corner",
+        description="Classify a vertex as a corner if the angle spanned by the two edges "
+        "connecting the two neighboring vertices is larger than the supplementary angle of "
+        "the threshold value.",
+        ge=0,
+        lt=np.pi,
+    )
+
+    distance_threshold: Optional[pd.PositiveFloat] = pd.Field(
+        None,
+        title="Distance Threshold In Defining Corner",
+        description="If not ``None``, discard the corner if its distance to neighbors "
+        "are below the threshold value, based on Douglas-Peucker algorithm.",
+    )
+
+    def corners(
+        self,
+        normal_axis: Axis,
+        coord: float,
+        structure_list: List[Structure],
+    ) -> ArrayFloat2D:
+        """On a 2D plane specified by axis = `normal_axis` and coordinate at `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.
+        """
+
+        # Construct plane
+        center = [0, 0, 0]
+        size = [inf, inf, inf]
+        center[normal_axis] = coord
+        size[normal_axis] = 0
+        plane = Box(center=center, size=size)
+
+        geometry_list = [
+            structure.geometry for structure in structure_list if isinstance(structure, Structure)
+        ]
+        # For metal, we don't distinguish between LossyMetal and PEC,
+        # so they'll be merged to PEC. Other materials are considered as dielectric.
+        medium_list = (
+            structure.medium for structure in structure_list if isinstance(structure, Structure)
+        )
+        medium_list = [
+            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)
+
+        # corner finder here
+        corner_list = []
+        for mat, shapes in merged_geos:
+            if self.medium != "all" and mat.is_pec != (self.medium == "metal"):
+                continue
+            polygon_list = ClipOperation.to_polygon_list(shapes)
+            for poly in polygon_list:
+                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)))
+                # in case the polygon has holes
+                for poly_inner in poly.interiors:
+                    corner_list.append(self._filter_collinear_vertices(list(poly_inner.coords)))
+        return np.concatenate(corner_list)
+
+    def _filter_collinear_vertices(self, vertices: ArrayFloat2D) -> ArrayFloat2D:
+        """Filter collinear vertices of a polygon, and return corners.
+
+        Parameters
+        ----------
+        vertices : ArrayFloat2D
+            Polygon vertices from shapely.Polygon. The last vertex is identical to the 1st
+            vertex to make a valid polygon.
+
+        Returns
+        -------
+        ArrayFloat2D
+            Corner coordinates.
+        """
+
+        def normalize(v):
+            return v / np.linalg.norm(v, axis=-1)[:, np.newaxis]
+
+        # drop the last vertex, which is identical to the 1st one.
+        vs_orig = np.array(vertices[:-1])
+        # compute unit vector to next and previous vertex
+        vs_next = np.roll(vs_orig, axis=0, shift=-1)
+        vs_previous = np.roll(vs_orig, axis=0, shift=+1)
+        unit_next = normalize(vs_next - vs_orig)
+        unit_previous = normalize(vs_previous - vs_orig)
+        # angle
+        angle = np.arccos(np.sum(unit_next * unit_previous, axis=-1))
+        ind_filter = angle <= np.pi - self.angle_threshold
+        return vs_orig[ind_filter]