Source code for gplugins.gmsh.uz_xsection_mesh

# type: ignore
from __future__ import annotations

from collections.abc import Sequence

import gdsfactory as gf
import numpy as np
from gdsfactory.config import get_number_of_cores
from gdsfactory.technology import LayerStack
from gdsfactory.typings import ComponentOrReference
from meshwell.gmsh_entity import GMSH_entity
from meshwell.model import Model
from shapely.geometry import LineString, MultiPolygon, Point, Polygon
from shapely.ops import unary_union

from gplugins.common.utils.parse_layer_stack import (
    list_unique_layer_stack_z,
    order_layer_stack,
)
from gplugins.gmsh.define_polysurfaces import define_polysurfaces
from gplugins.gmsh.parse_component import (
    process_buffers,
)
from gplugins.gmsh.parse_gds import cleanup_component, to_polygons


[docs] def get_u_bounds_polygons( polygons: MultiPolygon | list[Polygon], xsection_bounds: tuple[tuple[float, float], tuple[float, float]], u_offset: float = 0.0, ): """Performs the bound extraction given a (Multi)Polygon or [Polygon] and cross-sectional line coordinates. Args: polygons: shapely MultiPolygon or list of Polygons. xsection_bounds: ( (x1,y1), (x2,y2) ), with x1,y1 beginning point of cross-sectional line and x2,y2 the end. u_offset: amount to offset the returned polygons in the lateral dimension. Returns: list of bounding box coordinates (u1,u2)) in xsection line coordinates (distance from xsection_bounds[0]). """ line = LineString(xsection_bounds) linestart = Point(xsection_bounds[0]) return_list = [] for polygon in polygons.geoms if hasattr(polygons, "geoms") else [polygons]: initial_settings = np.seterr() np.seterr(invalid="ignore") intersection = polygon.intersection(line) np.seterr(**initial_settings) if intersection: for entry in ( intersection.geoms if hasattr(intersection, "geoms") else [intersection] ): bounds = entry.bounds p1 = Point([bounds[0], bounds[1]]) p2 = Point([bounds[2], bounds[3]]) return_list.append( [ linestart.distance(p1) + u_offset, linestart.distance(p2) + u_offset, ] ) return return_list
[docs] def get_u_bounds_layers( layer_polygons_dict: dict[tuple(str, str, str), MultiPolygon], xsection_bounds: tuple[tuple[float, float], tuple[float, float]], ): """Given a layer_polygons_dict and two coordinates (x1,y1), (x2,y2). Computes the bounding box(es) of each layer in the xsection coordinate system (u). Args: layer_polygons_dict: dict containing layernames: shapely polygons pairs xsection_bounds: ( (x1,y1), (x2,y2) ), with x1,y1 beginning point of cross-sectional line and x2,y2 the end. Returns: Dict containing layer(list pairs, with list a list of bounding box coordinates (u1,u2)) in xsection line coordinates. """ bounds_dict = {} for layername, ( gds_layername, next_layername, polygons, next_polygons, ) in layer_polygons_dict.items(): bounds_dict[layername] = [] bounds = get_u_bounds_polygons(polygons, xsection_bounds) next_bounds = get_u_bounds_polygons(next_polygons, xsection_bounds) if bounds: bounds_dict[layername] = ( gds_layername, next_layername, bounds, next_bounds, ) return bounds_dict
[docs] def get_uz_bounds_layers( layer_polygons_dict: dict[str, tuple[str, MultiPolygon, MultiPolygon]], xsection_bounds: tuple[tuple[float, float], tuple[float, float]], layer_stack: LayerStack, u_offset: float = 0.0, z_bounds: tuple[float, float] | None = None, ): """Given a component and layer stack, computes the bounding box(es). For each layer in the xsection coordinate system (u,z). Args: layer_polygons_dict: dict containing layernames: shapely polygons pairs. xsection_bounds: ( (x1,y1), (x2,y2) ), with x1,y1 beginning point of cross-sectional line and x2,y2 the end. layer_stack: LayerStack object containing layer information. u_offset: amount to offset the returned polygons in the lateral dimension. z_bounds: optional tuple containing the zmin and zmax of the simulation region. Returns: Dict containing layer: polygon pairs, with (u1,u2) in xsection line coordinates """ if z_bounds is not None: z_min_sim = z_bounds[0] z_max_sim = z_bounds[1] else: z_min_sim = -np.inf z_max_sim = np.inf # Get in-plane cross-sections inplane_bounds_dict = get_u_bounds_layers(layer_polygons_dict, xsection_bounds) outplane_bounds_dict = {} layer_dict = layer_stack.to_dict() # Remove empty entries inplane_bounds_dict = { key: value for (key, value) in inplane_bounds_dict.items() if value } for layername, ( gds_layername, next_layername, inplane_bounds_list, next_inplane_bounds_list, ) in inplane_bounds_dict.items(): outplane_polygons_list = [] for inplane_bounds, next_inplane_bounds in zip( inplane_bounds_list, next_inplane_bounds_list ): zmin = layer_dict[layername]["zmin"] zmax = layer_dict[next_layername]["zmin"] # Get bounding box umin_zmin = np.min(inplane_bounds) + u_offset umax_zmin = np.max(inplane_bounds) + u_offset umin_zmax = np.min(next_inplane_bounds) + u_offset umax_zmax = np.max(next_inplane_bounds) + u_offset # If layer thickness are negative, then zmin is actually zmax if zmin > zmax: foo = zmax zmax = zmin zmin = foo foo = umin_zmin umin_zmin = umin_zmax umin_zmax = foo foo = umax_zmax umax_zmax = umax_zmin umax_zmin = foo if zmax < z_min_sim or zmin > z_max_sim: # The whole polygon is outside of the simulation region continue # Check if we need to clip on top or bottom if zmin < z_min_sim: # Need to clip the bottom # Do a linear interpolation of the u coordinates umin_zmin = np.interp(z_min_sim, [zmin, zmax], [umin_zmin, umin_zmax]) umax_zmin = np.interp(z_min_sim, [zmin, zmax], [umax_zmin, umax_zmax]) zmin = z_min_sim if zmax > z_max_sim: # Need to clip the top umin_zmax = np.interp(z_max_sim, [zmin, zmax], [umin_zmin, umin_zmax]) umax_zmax = np.interp(z_max_sim, [zmin, zmax], [umax_zmin, umax_zmax]) zmax = z_max_sim points = [ [umin_zmin, zmin], [umin_zmax, zmax], [umax_zmax, zmax], [umax_zmin, zmin], ] outplane_polygons_list.append(Polygon(points)) outplane_bounds_dict[layername] = (gds_layername, outplane_polygons_list) return outplane_bounds_dict
[docs] def uz_xsection_mesh( component: ComponentOrReference, xsection_bounds: tuple[tuple[float, float], tuple[float, float]], layer_stack: LayerStack, layer_physical_map: dict, layer_meshbool_map: dict, resolutions: dict | None = None, default_characteristic_length: float = 0.5, background_tag: str | None = None, background_padding: Sequence[float, float, float, float, float, float] = (2.0,) * 6, background_mesh_order: int | float = 2**63 - 1, global_scaling: float = 1, global_scaling_premesh: float = 1, global_2D_algorithm: int = 6, filename: str | None = None, round_tol: int = 4, simplify_tol: float = 1e-4, u_offset: float = 0.0, left_right_periodic_bcs: bool = False, verbosity: int | None = 0, n_threads: int = get_number_of_cores(), gmsh_version: float | None = None, interface_delimiter: str = "___", background_remeshing_file=None, optimization_flags: tuple[tuple[str, int]] | None = None, **kwargs, ): """Mesh uz cross-section of component along line u = [[x1,y1] , [x2,y2]]. Args: component (Component): gdsfactory component to mesh xsection_bounds (Tuple): Tuple [[x1,y1] , [x2,y2]] parametrizing the line u layer_stack (LayerStack): gdsfactory LayerStack to parse layer_physical_map: map layer names to physical names layer_meshbool_map: map layer names to mesh_bool (True: mesh the prisms, False: don't mesh) resolutions (Dict): Pairs {"layername": {"resolution": float, "distance": "float}} to roughly control mesh refinement default_characteristic_length (float): default characteristic length to use in the mesh. background_tag (str): name of the background layer to add (default: no background added) background_padding (Tuple): [xleft, ydown, xright, yup] distances to add to the components and to fill with background_tag background_mesh_order (int, float): mesh order to assign to the background round_tol: during gds --> mesh conversion cleanup, number of decimal points at which to round the gdsfactory/shapely points before introducing to gmsh global_scaling: global scaling factor to apply to the mesh. global_scaling_premesh: global scaling factor to apply to the mesh before meshing. global_2D_algorithm: gmsh 2D meshing algorithm to use. filename: filename to save the mesh to. 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 left_right_periodic_bcs: if True, makes the left and right simulation edge meshes periodic verbosity: verbosity level of gmsh n_threads: number of threads to use in gmsh. gmsh_version: gmsh version to use. interface_delimiter: delimiter to use in the interface names background_remeshing_file: .pos file to use as a remeshing field. Overrides resolutions if not None. optimization_flags: list of tuples containing optimization flags to pass to gmsh. kwargs: additional arguments to pass to the meshing function. """ # Fuse and cleanup polygons of same layer in case user overlapped them layer_polygons_dict = cleanup_component( component, layer_stack, round_tol, simplify_tol ) # GDS polygons to simulation polygons buffered_layer_polygons_dict, buffered_layer_stack = process_buffers( layer_polygons_dict, layer_stack ) # simulation polygons to u-z coordinates along cross-sectional line z_bounds = kwargs.get("z_bounds", None) bounds_dict = get_uz_bounds_layers( buffered_layer_polygons_dict, xsection_bounds, buffered_layer_stack, u_offset, z_bounds=z_bounds, ) # u-z coordinates to gmsh-friendly polygons # Remove terminal layers and merge polygons layer_order = order_layer_stack(layer_stack) # gds layers shapes = {} for layername in layer_order: current_shapes = [] for gds_name, bounds in bounds_dict.values(): if gds_name == layername: layer_shapes = list(bounds) current_shapes.append(MultiPolygon(to_polygons(layer_shapes))) shapes[layername] = unary_union(MultiPolygon(to_polygons(current_shapes))) # Define polysurfaces model = Model(n_threads=n_threads) polysurfaces_list = define_polysurfaces( polygons_dict=shapes, layer_stack=layer_stack, model=model, scale_factor=global_scaling_premesh, resolutions=resolutions, layer_physical_map=layer_physical_map, layer_meshbool_map=layer_meshbool_map, ) # Add background polygon if background_tag is not None: # shapes[background_tag] = bounds.buffer(background_padding[0]) # bounds = unary_union(list(shapes.values())).bounds zs = list_unique_layer_stack_z(buffered_layer_stack) zmin = np.min(zs) - background_padding[1] zmax = np.max(zs) + background_padding[3] umin = -1 * background_padding[0] + u_offset umax = ( np.linalg.norm(np.array(xsection_bounds[1]) - np.array(xsection_bounds[0])) + background_padding[2] + u_offset ) background_box = GMSH_entity( gmsh_function=model.occ.add_rectangle, gmsh_function_kwargs={ "x": umin, "y": zmin, "z": 0, "dx": umax - umin, "dy": zmax - zmin, }, dimension=2, model=model, mesh_order=background_mesh_order, physical_name=background_tag, ) polysurfaces_list.append(background_box) # Create interface surfaces and boundaries minu = np.inf minz = np.inf maxu = -np.inf maxz = -np.inf if left_right_periodic_bcs: # Figure out bbox of simulation for polygon in shapes.values(): if not polygon.is_empty: minx, miny, maxx, maxy = polygon.bounds minu = min(minx, minu) minz = min(miny, minz) maxu = max(maxx, maxu) maxz = max(maxy, maxz) minu += u_offset maxu += u_offset if background_tag: minu -= background_padding[0] maxu += background_padding[2] minz -= background_padding[1] maxz += background_padding[3] # Create boundary rectangles left_line_rectangle = GMSH_entity( gmsh_function=model.occ.add_rectangle, gmsh_function_kwargs={ "x": minu - 1, "y": minz, "z": 0, "dx": 1, "dy": maxz - minz, }, dimension=2, model=model, mesh_order=0, physical_name="left_line", mesh_bool=False, ) right_line_rectangle = GMSH_entity( gmsh_function=model.occ.add_rectangle, gmsh_function_kwargs={ "x": maxu, "y": minz, "z": 0, "dx": 1, "dy": maxz - minz, }, dimension=2, model=model, mesh_order=0, physical_name="right_line", mesh_bool=False, ) # Give meshwell all possible boundary interfaces periodic_entities = [ ( f"left_line{interface_delimiter}{entity.physical_name}", f"right_line{interface_delimiter}{entity.physical_name}", ) for entity in polysurfaces_list ] polysurfaces_list.append(left_line_rectangle) polysurfaces_list.append(right_line_rectangle) # Mesh return model.mesh( entities_list=polysurfaces_list, default_characteristic_length=default_characteristic_length, periodic_entities=periodic_entities if left_right_periodic_bcs else None, global_scaling=global_scaling, global_2D_algorithm=global_2D_algorithm, gmsh_version=gmsh_version, filename=filename, verbosity=verbosity, interface_delimiter=interface_delimiter, background_remeshing_file=background_remeshing_file, optimization_flags=optimization_flags, )
if __name__ == "__main__": from gdsfactory.pdk import get_layer_stack c = gf.Component() waveguide = c << gf.get_component(gf.components.straight_pin(length=10, taper=None)) undercut = c << gf.get_component( gf.components.rectangle( size=(5.0, 5.0), layer="UNDERCUT", centered=True, ) ) undercut.dmove((4, 0)) c.show() filtered_layer_stack = LayerStack( layers={ k: get_layer_stack().layers[k] for k in ( "slab90", "core", "via_contact", "undercut", "box", # "substrate", "clad", "metal1", ) # "slab90", "via_contact")#"via_contact") # "slab90", "core" } ) resolutions = {} resolutions["si"] = {"resolution": 0.02, "distance": 0.1} resolutions["siInt"] = {"resolution": 0.0003, "distance": 0.1} # resolutions["sio2"] = {"resolution": 0.5, "distance": 1} resolutions["Aluminum"] = {"resolution": 0.1, "distance": 1} # resolutions = {} # resolutions["core"] = {"resolution": 0.02, "distance": 0.1} # resolutions["coreInt"] = {"resolution": 0.0001, "distance": 0.1} # resolutions["slab90"] = {"resolution": 0.05, "distance": 0.1} # # resolutions["sio2"] = {"resolution": 0.5, "distance": 1} # resolutions["via_contact"] = {"resolution": 0.1, "distance": 1} geometry = uz_xsection_mesh( c, [(4, -15), (4, 15)], filtered_layer_stack, resolutions=resolutions, layer_physical_map={}, layer_meshbool_map={}, background_tag="Oxide", background_padding=(0, 0, 0, 0), filename="mesh.msh", round_tol=3, simplify_tol=1e-3, u_offset=-15, left_right_periodic_bcs=True, ) import meshio mesh_from_file = meshio.read("mesh.msh") def create_mesh(mesh, cell_type, prune_z=False): cells = mesh.get_cells_type(cell_type) cell_data = mesh.get_cell_data("gmsh:physical", cell_type) points = mesh.points[:, :2] if prune_z else mesh.points return meshio.Mesh( points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data]}, ) line_mesh = create_mesh(mesh_from_file, "line", prune_z=True) meshio.write("facet_mesh.xdmf", line_mesh) triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True) meshio.write("mesh.xdmf", triangle_mesh)