Equalize vertex processing in Geometry and PolySlab

Use the same grid size for both when loading vertices from GDSII, and
apply dilations in the same manner to avoid unexpected differences due
to operation order.

Signed-off-by: Lucas Heitzmann Gabrielli <lucas@flexcompute.com>

tests/test_components/test_geometry.py

@@ -561,7 +561,7 @@ def test_from_gds():
         cell, 2, (0, 1), gds_layer=0, dilation=-0.5, sidewall_angle=0.5, reference_plane="bottom"
     )
     assert len(geo.intersections_plane(z=0)) == 2
-    assert len(geo.intersections_plane(z=1)) == 5
+    assert len(geo.intersections_plane(z=1)) == 1
 
 
 def test_custom_surface_geometry(tmp_path):

tidy3d/components/geometry/base.py

@@ -23,6 +23,7 @@
 from ...exceptions import SetupError, ValidationError, Tidy3dKeyError, Tidy3dError
 from ...constants import MICROMETER, LARGE_NUMBER, RADIAN, inf, fp_eps
 
+POLY_GRID_SIZE = 1e-12
 
 class Geometry(Tidy3dBaseModel, ABC):
     """Abstract base class, defines where something exists in space."""
@@ -921,12 +922,16 @@ def from_gds(
         geometries = []
         with log as consolidated_logger:
             for vertices in gds_loader_fn(gds_cell, gds_layer, gds_dtype, gds_scale):
-                # buffer(0) is necessary to merge self-intersections before dilation/erosion
-                shape = shapely.Polygon(vertices).buffer(0).buffer(dilation)
+                # buffer(0) is necessary to merge self-intersections
+                shape = shapely.set_precision(shapely.Polygon(vertices).buffer(0), POLY_GRID_SIZE)
                 try:
                     geometries.append(
-                        from_shapely(shape, axis, slab_bounds, sidewall_angle, reference_plane)
+                        from_shapely(
+                            shape, axis, slab_bounds, dilation, sidewall_angle, reference_plane
+                        )
                     )
+                except pydantic.ValidationError as error:
+                    consolidated_logger.warning(str(error))
                 except Tidy3dError as error:
                     consolidated_logger.warning(str(error))
         return GeometryGroup(geometries=geometries)
@@ -936,6 +941,7 @@ def from_shapely(
         shape: Shapely,
         axis: Axis,
         slab_bounds: Tuple[float, float],
+        dilation: float = 0.0,
         sidewall_angle: float = 0,
         reference_plane: PlanePosition = "middle",
     ) -> Geometry:
@@ -950,6 +956,9 @@ def from_shapely(
             Integer index defining the extrusion axis: 0 (x), 1 (y), or 2 (z).
         slab_bounds: Tuple[float, float]
             Minimal and maximal positions of the extruded slab along ``axis``.
+        dilation : float
+            Dilation of the polygon in the base by shifting each edge along its normal outwards
+            direction by a distance; a negative value corresponds to erosion.
         sidewall_angle : float = 0
             Angle of the extrusion sidewalls, away from the vertical direction, in radians. Positive
             (negative) values result in slabs larger (smaller) at the base than at the top.
@@ -964,7 +973,7 @@ def from_shapely(
         :class:`Geometry`
             Geometry extruded from the 2D data.
         """
-        return from_shapely(shape, axis, slab_bounds, sidewall_angle, reference_plane)
+        return from_shapely(shape, axis, slab_bounds, dilation, sidewall_angle, reference_plane)
 
     def _as_union(self) -> List[Geometry]:
         """Return a list of geometries that, united, make up the given geometry."""

tidy3d/components/geometry/polyslab.py

@@ -24,9 +24,6 @@
 
 _IS_CLOSE_RTOL = np.finfo(float).eps
 
-# polygon merge
-_POLY_GRID_SIZE = 1e-12
-
 # Warn for too many divided polyslabs
 _COMPLEX_POLYSLAB_DIVISIONS_WARN = 100
 
@@ -335,7 +332,7 @@ def _load_gds_vertices(
 
         # convert vertices into polyslabs
         polygons = [shapely.Polygon(vertices).buffer(0) for vertices in all_vertices]
-        polys_union = shapely.unary_union(polygons, grid_size=_POLY_GRID_SIZE)
+        polys_union = shapely.unary_union(polygons, grid_size=base.POLY_GRID_SIZE)
 
         if polys_union.geom_type == "Polygon":
             all_vertices = [np.array(polys_union.exterior.coords)]
@@ -643,7 +640,7 @@ def _intersections_side(self, position, axis) -> list:
             h_base = h_top
 
         # merge touching polygons
-        polys_union = shapely.unary_union(polys, grid_size=_POLY_GRID_SIZE)
+        polys_union = shapely.unary_union(polys, grid_size=base.POLY_GRID_SIZE)
         if polys_union.geom_type == "Polygon":
             return [polys_union]
         if polys_union.geom_type == "MultiPolygon":

tidy3d/components/geometry/utils.py

@@ -22,10 +22,12 @@
 ]
 
 
+# pylint:disable=too-many-arguments
 def from_shapely(
     shape: Shapely,
     axis: Axis,
     slab_bounds: Tuple[float, float],
+    dilation: float = 0.0,
     sidewall_angle: float = 0,
     reference_plane: PlanePosition = "middle",
 ) -> base.Geometry:
@@ -40,6 +42,9 @@ def from_shapely(
         Integer index defining the extrusion axis: 0 (x), 1 (y), or 2 (z).
     slab_bounds: Tuple[float, float]
         Minimal and maximal positions of the extruded slab along ``axis``.
+    dilation : float
+        Dilation of the polygon in the base by shifting each edge along its normal outwards
+        direction by a distance; a negative value corresponds to erosion.
     sidewall_angle : float = 0
         Angle of the extrusion sidewalls, away from the vertical direction, in radians. Positive
         (negative) values result in slabs larger (smaller) at the base than at the top.
@@ -58,23 +63,27 @@ def from_shapely(
         if sidewall_angle == 0:
             return polyslab.PolySlab(
                 vertices=shape.coords[:-1],
-                slab_bounds=slab_bounds,
                 axis=axis,
+                slab_bounds=slab_bounds,
+                dilation=dilation,
                 reference_plane=reference_plane,
             )
         group = polyslab.ComplexPolySlabBase(
             vertices=shape.coords[:-1],
-            slab_bounds=slab_bounds,
             axis=axis,
+            slab_bounds=slab_bounds,
+            dilation=dilation,
             sidewall_angle=sidewall_angle,
             reference_plane=reference_plane,
         ).geometry_group
         return group.geometries[0] if len(group.geometries) == 1 else group
 
     if shape.geom_type == "Polygon":
-        exterior = from_shapely(shape.exterior, axis, slab_bounds, sidewall_angle, reference_plane)
+        exterior = from_shapely(
+            shape.exterior, axis, slab_bounds, dilation, sidewall_angle, reference_plane
+        )
         interior = [
-            from_shapely(hole, axis, slab_bounds, -sidewall_angle, reference_plane)
+            from_shapely(hole, axis, slab_bounds, -dilation, -sidewall_angle, reference_plane)
             for hole in shape.interiors
         ]
         if len(interior) == 0:
@@ -85,7 +94,7 @@ def from_shapely(
     if shape.geom_type in {"MultiPolygon", "GeometryCollection"}:
         return base.GeometryGroup(
             geometries=[
-                from_shapely(geo, axis, slab_bounds, sidewall_angle, reference_plane)
+                from_shapely(geo, axis, slab_bounds, dilation, sidewall_angle, reference_plane)
                 for geo in shape.geoms
             ]
         )