fix polygons at gdstk level before shapely meshing preprocessing

gdsfactory/simulation/gmsh/mesh.py

@@ -19,7 +19,7 @@
 )
 
 
-def define_entities(model, shapes_dict: OrderedDict):
+def define_entities(model, shapes_dict: OrderedDict, atol=1e-3):
     """Adds the polygons and lines from a "shapes_dict" as physical entities in the pygmsh model "model".
 
     Args:
@@ -33,7 +33,7 @@ def define_entities(model, shapes_dict: OrderedDict):
     # Break up lines and polygon edges so that plane is tiled with no partially overlapping line segments
     polygons_broken_dict, lines_broken_dict = break_geometry(shapes_dict)
     # Instantiate shapely to gmsh map
-    meshtracker = MeshTracker(model=model, atol=1e-3)
+    meshtracker = MeshTracker(model=model, atol=atol)
     # Add lines, reusing line segments
     model, meshtracker = add_lines(model, meshtracker, lines_broken_dict)
     # Add surfaces, reusing lines
@@ -126,6 +126,7 @@ def mesh_from_polygons(
     global_meshsize_array: Optional[np.array] = None,
     global_meshsize_interpolant_func: Optional[callable] = NearestNDInterpolator,
     verbosity: Optional[bool] = False,
+    atol: Optional[float] = 1e-4,
 ):
     """Return a 2D mesh from an ordered dict of shapely polygons.
 
@@ -140,6 +141,7 @@ def mesh_from_polygons(
         global_meshsize_array: array [x,y,z,lc] defining local mesh sizes. Not used if None
         global_meshsize_interpolant: interpolation function for array [x,y,z,lc]. Default scipy.interpolate.NearestNDInterpolator
         verbosity: boolean, gmsh stdout as it meshes
+        atol: tolerance used to establish equivalency between vertices
     """
     global_meshsize_callback_bool = global_meshsize_array is not None
 
@@ -153,7 +155,7 @@ def mesh_from_polygons(
     (
         model,
         meshtracker,
-    ) = define_entities(model, shapes_dict)
+    ) = define_entities(model, shapes_dict, atol)
 
     # Synchronize
     model.synchronize()

gdsfactory/simulation/gmsh/meshtracker.py

@@ -9,7 +9,7 @@
 
 
 class MeshTracker:
-    def __init__(self, model, atol=1e-3):
+    def __init__(self, model, atol=1e-4):
         """Map between shapely and gmsh.
 
         Shapely is useful for built-in geometry equivalencies and extracting orientation, instead of doing it manually

gdsfactory/simulation/gmsh/parse_gds.py

@@ -1,11 +1,12 @@
 """Preprocessing involving mostly the GDS polygons."""
 from __future__ import annotations
 
+import gdsfactory as gf
 import shapely
 from shapely.geometry import LineString, MultiLineString, MultiPolygon, Polygon
 
 
-def round_coordinates(geom, ndigits=5):
+def round_coordinates(geom, ndigits=4):
     """Round coordinates to n_digits to eliminate floating point errors."""
 
     def _round_coords(x, y, z=None):
@@ -20,9 +21,21 @@ def _round_coords(x, y, z=None):
     return shapely.ops.transform(_round_coords, geom)
 
 
-def fuse_polygons(component, layername, layer, round_tol=5, simplify_tol=1e-5):
+def fuse_polygons(
+    component, layername, layer, round_tol=4, simplify_tol=1e-4, offset_tol=None
+):
     """Take all polygons from a layer, and returns a single (Multi)Polygon shapely object."""
     layer_component = component.extract(layer)
+
+    # gdstk union before shapely conversion helps with ill-formed polygons
+    offset_tol = offset_tol or gf.get_active_pdk().grid_size
+    layer_component = gf.geometry.offset(
+        layer_component, distance=offset_tol, precision=1e-6, layer=layer
+    )
+    layer_component = gf.geometry.offset(
+        layer_component, distance=-offset_tol, precision=1e-6, layer=layer
+    )
+
     shapely_polygons = [
         round_coordinates(shapely.geometry.Polygon(polygon), round_tol)
         for polygon in layer_component.get_polygons()

gdsfactory/simulation/gmsh/uz_xsection_mesh.py

@@ -166,9 +166,10 @@ def uz_xsection_mesh(
     extra_shapes_dict: Optional[OrderedDict] = None,
     merge_by_material: Optional[bool] = False,
     interface_surfaces: Optional[Dict[str, Tuple(float, float)]] = None,
-    round_tol: int = 3,
-    simplify_tol: float = 1e-2,
+    round_tol: int = 4,
+    simplify_tol: float = 1e-4,
     u_offset: float = 0.0,
+    atol: Optional[float] = 1e-5,
     **kwargs,
 ):
     """Mesh uz cross-section of component along line u = [[x1,y1] , [x2,y2]].
@@ -191,6 +192,7 @@ def uz_xsection_mesh(
         round_tol: during gds --> mesh conversion cleanup, number of decimal points at which to round the gdsfactory/shapely points before introducing to gmsh
         simplify_tol: during gds --> mesh conversion cleanup, shapely "simplify" tolerance (make it so all points are at least separated by this amount)
         u_offset: quantity to add to the "u" coordinates, useful to convert back to x or y if parallel to those axes
+        atol: tolerance used to establish equivalency between vertices
     """
     interface_surfaces = interface_surfaces or {}
 
@@ -275,6 +277,7 @@ def uz_xsection_mesh(
         default_resolution_max=default_resolution_max,
         global_meshsize_array=global_meshsize_array,
         global_meshsize_interpolant_func=global_meshsize_interpolant_func,
+        atol=atol,
     )
 
 

gdsfactory/simulation/gmsh/xy_xsection_mesh.py

@@ -37,8 +37,9 @@ def xy_xsection_mesh(
     global_meshsize_interpolant_func: Optional[callable] = NearestNDInterpolator,
     extra_shapes_dict: Optional[OrderedDict] = None,
     merge_by_material: Optional[bool] = False,
-    round_tol: int = 3,
-    simplify_tol: float = 1e-2,
+    round_tol: int = 4,
+    simplify_tol: float = 1e-4,
+    atol: Optional[float] = 1e-5,
 ):
     """Mesh xy cross-section of component at height z.
 
@@ -59,6 +60,7 @@ def xy_xsection_mesh(
         merge_by_material: boolean, if True will merge polygons from layers with the same layer.material. Physical keys will be material in this case.
         round_tol: during gds --> mesh conversion cleanup, number of decimal points at which to round the gdsfactory/shapely points before introducing to gmsh
         simplify_tol: during gds --> mesh conversion cleanup, shapely "simplify" tolerance (make it so all points are at least separated by this amount)
+        atol: tolerance used to establish equivalency between vertices
     """
     # Fuse and cleanup polygons of same layer in case user overlapped them
     layer_polygons_dict = cleanup_component(
@@ -113,6 +115,7 @@ def xy_xsection_mesh(
         default_resolution_max=default_resolution_max,
         global_meshsize_array=global_meshsize_array,
         global_meshsize_interpolant_func=global_meshsize_interpolant_func,
+        atol=atol,
     )
 
 

gdsfactory/simulation/gmsh/xyz_mesh.py

@@ -209,7 +209,7 @@ def xyz_mesh(
     verbosity: Optional[bool] = False,
     override_volumes: Optional[Dict] = None,
     round_tol: int = 3,
-    simplify_tol: float = 1e-2,
+    simplify_tol: float = 1e-3,
 ):
     """Full 3D mesh of component.