fix(Medium2D): fix subdivision when merged geometry is a polygon with holes

CHANGELOG.md

@@ -81,6 +81,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Fixed DRC parsing for quoted categories in klayout plugin.
 - Removed mode solver warnings about evaluating permittivity of a `Medium2D`.
 - Maximum number of grid points in an EME simulation is now based solely on transverse grid points. Maximum number of EME cells is unchanged.
+- Fixed handling of polygons with holes (interiors) in `subdivide()` function. The function now properly converts polygons with interiors into `PolySlab` geometries using subtraction operations with `GeometryGroup`.
 
 ### Removed
 - Removed deprecated `use_complex_fields` parameter from `TwoPhotonAbsorption` and `KerrNonlinearity`. Parameters `beta` and `n2` are now real-valued only, as is `n0` if specified.

tests/test_components/test_geometry.py

@@ -1158,6 +1158,71 @@ def test_subdivide():
     subdivisions = subdivide(geom=overlapping_boxes, structures=[background_structure, box_sliver])
 
 
+def test_subdivide_geometry_group_with_polygon_holes():
+    """Test that unionized geometry containing a hole works correctly."""
+    mm = 1000.0
+
+    # Create four boxes arranged to form a cross pattern with a square hole in the middle
+    box_a = td.Box.from_bounds((-10 * mm, -10 * mm, 0 * mm), (-5 * mm, 5 * mm, 0 * mm))
+    box_b = td.Box.from_bounds((-10 * mm, 5 * mm, 0 * mm), (5 * mm, 10 * mm, 0 * mm))
+    box_c = td.Box.from_bounds((5 * mm, -5 * mm, 0 * mm), (10 * mm, 10 * mm, 0 * mm))
+    box_d = td.Box.from_bounds((-5 * mm, -10 * mm, 0 * mm), (10 * mm, -5 * mm, 0 * mm))
+
+    geom_group = td.GeometryGroup(geometries=(box_a, box_b, box_c, box_d))
+    geom_structures_group = [td.Structure(geometry=geom_group, medium=td.PEC2D)]
+
+    feed_pin_bottom = -2 * mm
+    feed_pin_top = 0 * mm
+    feed_pin_center = 0.5 * (feed_pin_top + feed_pin_bottom)
+    feed_pin_length = feed_pin_top - feed_pin_bottom
+    feed_center_x = -7.5 * mm
+    feed_center_y = -7.5 * mm
+    rfeed = 1.0 * mm
+
+    feed_pin = td.Structure(
+        geometry=td.Cylinder(
+            center=(feed_center_x, feed_center_y, feed_pin_center),
+            radius=rfeed,
+            length=feed_pin_length,
+            axis=2,
+        ),
+        medium=td.PECMedium(),
+    )
+
+    structures_list = [feed_pin, *geom_structures_group]
+
+    freq = 1500 * 1e6
+    dl = (td.C_0 / freq) / 300.0
+
+    mesh_overrides = [
+        td.MeshOverrideStructure(
+            geometry=td.Box(
+                center=(0, 0, 0 * mm),
+                size=(20 * mm, 20 * mm, 6 * mm),
+            ),
+            dl=[dl, dl, dl],
+        )
+    ]
+
+    sim = td.Simulation(
+        size=[100 * mm, 100 * mm, 30 * mm],
+        grid_spec=td.GridSpec.auto(
+            min_steps_per_wvl=20,
+            wavelength=td.C_0 / freq,
+            override_structures=mesh_overrides,
+        ),
+        structures=structures_list,
+        run_time=1e-13,
+    )
+
+    contains_difference_operation = False
+    for structure in sim._finalized.structures:
+        geo = structure.geometry
+        if isinstance(geo, td.ClipOperation) and geo.operation == "difference":
+            contains_difference_operation = True
+    assert contains_difference_operation
+
+
 @pytest.mark.parametrize("snap_location", [SnapLocation.Boundary, SnapLocation.Center])
 @pytest.mark.parametrize(
     "snap_behavior",

tidy3d/components/geometry/base.py

@@ -5,7 +5,7 @@
 import functools
 import pathlib
 from abc import ABC, abstractmethod
-from collections.abc import Iterable
+from collections.abc import Iterable, Sequence
 from os import PathLike
 from typing import TYPE_CHECKING, Any, Callable, Optional, Union
 
@@ -721,11 +721,14 @@ def evaluate_inf_shape(shape: Shapely) -> Shapely:
         if not any(np.isinf(b) for b in shape.bounds):
             return shape
 
+        def _processed_coords(coords: Sequence[tuple[Any, ...]]) -> list[tuple[float, ...]]:
+            evaluated = Geometry._evaluate_inf(np.array(coords))
+            return [tuple(point) for point in evaluated.tolist()]
+
         if shape.geom_type == "Polygon":
-            return shapely.Polygon(
-                Geometry._evaluate_inf(np.array(shape.exterior.coords)),
-                [Geometry._evaluate_inf(np.array(g.coords)) for g in shape.interiors],
-            )
+            shell = _processed_coords(shape.exterior.coords)
+            holes = [_processed_coords(g.coords) for g in shape.interiors]
+            return shapely.Polygon(shell, holes)
         if shape.geom_type in {"Point", "LineString", "LinearRing"}:
             return shape.__class__(Geometry._evaluate_inf(np.array(shape.coords)))
         if shape.geom_type in {

tidy3d/components/geometry/utils_2d.py

@@ -7,7 +7,7 @@
 import numpy as np
 import shapely
 
-from tidy3d.components.geometry.base import Box, ClipOperation, Geometry
+from tidy3d.components.geometry.base import Box, ClipOperation, Geometry, GeometryGroup
 from tidy3d.components.geometry.polyslab import _MIN_POLYGON_AREA, PolySlab
 from tidy3d.components.grid.grid import Grid
 from tidy3d.components.scene import Scene
@@ -122,10 +122,28 @@ def subdivide(
 
     """
 
-    def shapely_to_polyslab(polygon: shapely.Polygon, axis: Axis, center: float) -> PolySlab:
-        xx, yy = polygon.exterior.coords.xy
-        vertices = list(zip(xx, yy))
-        return PolySlab(slab_bounds=(center, center), vertices=vertices, axis=axis)
+    def shapely_to_polyslab(polygon: shapely.Polygon, axis: Axis, center: float) -> Geometry:
+        def ring_vertices(ring: shapely.LinearRing) -> list[tuple[float, float]]:
+            xx, yy = ring.coords.xy
+            return list(zip(xx, yy))
+
+        polyslab = PolySlab(
+            slab_bounds=(center, center),
+            vertices=ring_vertices(polygon.exterior),
+            axis=axis,
+        )
+        if len(polygon.interiors) == 0:
+            return polyslab
+
+        interiors = [
+            PolySlab(
+                slab_bounds=(center, center),
+                vertices=ring_vertices(interior),
+                axis=axis,
+            )
+            for interior in polygon.interiors
+        ]
+        return polyslab - GeometryGroup(geometries=interiors)
 
     def to_multipolygon(shapely_geometry: Shapely) -> shapely.MultiPolygon:
         return shapely.MultiPolygon(ClipOperation.to_polygon_list(shapely_geometry))