feat(meshing): cache merged_geos in gap refinement

feat(meshing): gap refinement only insert snapping points

Author: weiliangjin2021

tests/test_components/test_layerrefinement.py

@@ -104,19 +104,19 @@ def test_layerrefinement():
         assert layer.size[axis] == 1
         assert layer.size[(axis + 1) % 3] == td.inf
         assert layer.size[(axis + 2) % 3] == td.inf
-        assert not layer._is_inplane_bounded
+        assert not layer._is_inplane_bounded(layer)
 
     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 layer._is_inplane_bounded
+    assert layer._is_inplane_bounded(layer)
 
     # from structures
     structures = [td.Structure(geometry=td.Box(size=(td.inf, 2, 3)), medium=td.Medium())]
     layer = LayerRefinementSpec.from_structures(structures)
-    assert layer._is_inplane_bounded
+    assert layer._is_inplane_bounded(layer)
     assert layer.axis == 1
 
     with pytest.raises(pydantic.ValidationError):
@@ -138,11 +138,11 @@ def test_layerrefinement():
 def test_layerrefinement_inplane_inside():
     # inplane inside
     layer = LayerRefinementSpec.from_layer_bounds(axis=2, bounds=(0, 1))
-    assert not layer._is_inplane_bounded
-    assert layer._inplane_inside([3e3, 4e4])
+    assert not layer._is_inplane_bounded(layer)
+    assert layer._inplane_inside(layer, [3e3, 4e4])
     layer = LayerRefinementSpec(axis=1, size=(1, 0, 1))
-    assert layer._inplane_inside([0, 0])
-    assert not layer._inplane_inside([2, 0])
+    assert layer._inplane_inside(layer, [0, 0])
+    assert not layer._inplane_inside(layer, [2, 0])
 
 
 def test_layerrefinement_snapping_points():
@@ -263,7 +263,11 @@ def count_grids_within_layer(sim_t):
     assert (
         len(
             sim2.grid_spec.all_override_structures(
-                list(sim2.structures), 1.0, sim2.size, lumped_elements
+                list(sim2.structures),
+                1.0,
+                lumped_elements,
+                sim2._internal_layerrefinement_boundary_types,
+                sim2.bounds,
             )
         )
         == 1
@@ -281,7 +285,11 @@ def count_grids_within_layer(sim_t):
     assert (
         len(
             sim2.grid_spec.all_override_structures(
-                list(sim2.structures), 1.0, sim2.size, lumped_elements
+                list(sim2.structures),
+                1.0,
+                lumped_elements,
+                sim2._internal_layerrefinement_boundary_types,
+                sim2.bounds,
             )
         )
         == 2
@@ -425,7 +433,11 @@ def test_dl_min_from_smallest_feature():
         ),
         medium=td.PECMedium(),
     )
-
+    boundary_types = [[None] * 2] * 3
+    sim_bounds = [
+        [-td.inf] * 3,
+        [td.inf] * 3,
+    ]
     # check expected dl_min
     layer_spec = td.LayerRefinementSpec(
         axis=2,
@@ -434,15 +446,15 @@ def test_dl_min_from_smallest_feature():
             convex_resolution=10,
         ),
     )
-    dl_min = layer_spec._dl_min_from_smallest_feature([structure])
+    dl_min = layer_spec._dl_min_from_smallest_feature([structure], sim_bounds, boundary_types)
     assert np.allclose(0.3 / 10, dl_min)
 
     layer_spec = td.LayerRefinementSpec(
         axis=2,
         size=(td.inf, td.inf, 2),
         corner_finder=td.CornerFinderSpec(mixed_resolution=10),
     )
-    dl_min = layer_spec._dl_min_from_smallest_feature([structure])
+    dl_min = layer_spec._dl_min_from_smallest_feature([structure], sim_bounds, boundary_types)
     assert np.allclose(0.2 / 10, dl_min)
 
     layer_spec = td.LayerRefinementSpec(
@@ -452,7 +464,7 @@ def test_dl_min_from_smallest_feature():
             concave_resolution=10,
         ),
     )
-    dl_min = layer_spec._dl_min_from_smallest_feature([structure])
+    dl_min = layer_spec._dl_min_from_smallest_feature([structure], sim_bounds, boundary_types)
     assert np.allclose(0.1 / 10, dl_min)
 
     # check grid is generated succesfully
@@ -517,7 +529,8 @@ def test_gap_meshing():
             layer_refinement_specs=[
                 td.LayerRefinementSpec(
                     axis=2,
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=num_iters,
                     size=[td.inf, td.inf, 2],
                     center=[0, 0, 1],
@@ -583,7 +596,8 @@ def test_gap_meshing():
                 td.LayerRefinementSpec(
                     axis=1,
                     size=(td.inf, 0.2, td.inf),
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=1,
                     dl_min_from_gap_width=True,
                 )
@@ -626,7 +640,8 @@ def test_gap_meshing():
                 td.LayerRefinementSpec(
                     axis=1,
                     size=(td.inf, 0.2, td.inf),
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=2,
                     dl_min_from_gap_width=True,
                     interior_disjoint_geometries=False,
@@ -658,7 +673,8 @@ def test_gap_meshing():
                 td.LayerRefinementSpec(
                     axis=1,
                     size=(td.inf, 0.2, td.inf),
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=1,
                     dl_min_from_gap_width=True,
                 )
@@ -697,7 +713,8 @@ def test_gap_meshing():
                     axis=1,
                     size=(0.5, 0.2, 0.5),
                     center=(0.5, 0, -0.5),
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=1,
                     dl_min_from_gap_width=True,
                 )
@@ -731,7 +748,8 @@ def test_gap_meshing():
                     axis=1,
                     size=(0.05, 0.2, 0.05),
                     center=(0.05, 0, 0.05),
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=2,
                     dl_min_from_gap_width=True,
                 )
@@ -765,7 +783,8 @@ def test_gap_meshing():
                     axis=1,
                     size=(0.01, 0.2, 0.01),
                     center=(10.0, 0, 10.0),
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=1,
                     dl_min_from_gap_width=True,
                 )
@@ -804,7 +823,8 @@ def test_gap_meshing():
                 td.LayerRefinementSpec(
                     axis=1,
                     size=(td.inf, 0.2, td.inf),
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=1,
                     dl_min_from_gap_width=True,
                 )
@@ -842,7 +862,8 @@ def test_gap_meshing():
                 td.LayerRefinementSpec(
                     axis=0,
                     size=(0.2, td.inf, td.inf),
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=1,
                     dl_min_from_gap_width=True,
                 )

tidy3d/components/geometry/base.py

@@ -31,6 +31,7 @@
 from tidy3d.components.autograd.derivative_utils import DerivativeInfo
 from tidy3d.components.base import Tidy3dBaseModel, cached_property
 from tidy3d.components.geometry.bound_ops import bounds_intersection, bounds_union
+from tidy3d.components.geometry.float_utils import increment_float
 from tidy3d.components.transformation import ReflectionFromPlane, RotationAroundAxis
 from tidy3d.components.types import (
     ArrayFloat2D,
@@ -2173,6 +2174,11 @@ def intersections_with(self, other: Shapely) -> list[Shapely]:
         shapely_box = Geometry.evaluate_inf_shape(shapely_box)
         return [Geometry.evaluate_inf_shape(shape) & shapely_box for shape in shapes_plane]
 
+    def slightly_enlarged_copy(self) -> Box:
+        """Box size slightly enlarged around machine precision."""
+        size = [increment_float(orig_length, 1) for orig_length in self.size]
+        return self.updated_copy(size=size)
+
     def padded_copy(
         self,
         x: Optional[tuple[pydantic.NonNegativeFloat, pydantic.NonNegativeFloat]] = None,

tidy3d/components/geometry/float_utils.py

@@ -0,0 +1,31 @@
+"""Utilities for float manipulation."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from tidy3d.constants import inf
+
+
+def increment_float(val: float, sign: int) -> float:
+    """Applies a small positive or negative shift as though `val` is a 32bit float
+    using numpy.nextafter, but additionally handles some corner cases.
+    """
+    # Infinity is left unchanged
+    if val == inf or val == -inf:
+        return val
+
+    if sign >= 0:
+        sign = 1
+    else:
+        sign = -1
+
+    # Avoid small increments within subnormal values
+    if np.abs(val) <= np.finfo(np.float32).tiny:
+        return val + sign * np.finfo(np.float32).tiny
+
+    # Numpy seems to skip over the increment from -0.0 and +0.0
+    # which is different from c++
+    val_inc = np.nextafter(val, sign * inf, dtype=np.float32)
+
+    return np.float32(val_inc)

tidy3d/components/geometry/utils_2d.py

@@ -8,36 +8,13 @@
 import shapely
 
 from tidy3d.components.geometry.base import Box, ClipOperation, Geometry, GeometryGroup
+from tidy3d.components.geometry.float_utils import increment_float
 from tidy3d.components.geometry.polyslab import _MIN_POLYGON_AREA, PolySlab
 from tidy3d.components.grid.grid import Grid
 from tidy3d.components.scene import Scene
 from tidy3d.components.structure import Structure
 from tidy3d.components.types import Axis, Shapely
-from tidy3d.constants import fp_eps, inf
-
-
-def increment_float(val: float, sign: int) -> float:
-    """Applies a small positive or negative shift as though `val` is a 32bit float
-    using numpy.nextafter, but additionally handles some corner cases.
-    """
-    # Infinity is left unchanged
-    if val == inf or val == -inf:
-        return val
-
-    if sign >= 0:
-        sign = 1
-    else:
-        sign = -1
-
-    # Avoid small increments within subnormal values
-    if np.abs(val) <= np.finfo(np.float32).tiny:
-        return val + sign * np.finfo(np.float32).tiny
-
-    # Numpy seems to skip over the increment from -0.0 and +0.0
-    # which is different from c++
-    val_inc = np.nextafter(val, sign * inf, dtype=np.float32)
-
-    return np.float32(val_inc)
+from tidy3d.constants import fp_eps
 
 
 def snap_coordinate_to_grid(grid: Grid, center: float, axis: Axis) -> float: