diff --git a/docs/api.md b/docs/api.md index e5fec09..d729090 100644 --- a/docs/api.md +++ b/docs/api.md @@ -178,6 +178,8 @@ ::: gsim.meep.Material +::: gsim.meep.MeepMeshSim + ::: gsim.meep.ModeSource ::: gsim.meep.ResolutionConfig diff --git a/nbs/mesh_ybranch.ipynb b/nbs/mesh_ybranch.ipynb new file mode 100644 index 0000000..7848041 --- /dev/null +++ b/nbs/mesh_ybranch.ipynb @@ -0,0 +1,268 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# GMSH Mesh Generation via MeepMeshSim\n", + "\n", + "This notebook demonstrates using `MeepMeshSim` to generate a GMSH mesh\n", + "from any gdsfactory component + LayerStack — mirroring Palace's\n", + "`mesh()` → `write_config()` workflow for photonic simulations.\n", + "\n", + "We use a UBC PDK Y-branch (photonic SOI) as an example." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "### Load component from UBC PDK" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "from ubcpdk import PDK, cells\n", + "\n", + "PDK.activate()\n", + "\n", + "# c = cells.ebeam_y_1550()\n", + "# c = cells.ring_single()\n", + "c = cells.coupler()\n", + "\n", + "dirname = \"./sim-data-mesh-----coupler\"\n", + "\n", + "c" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "### Configure MeepMeshSim" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "from gsim.common.stack.extractor import Layer, LayerStack\n", + "from gsim.meep import MeepMeshSim\n", + "\n", + "# Build a LayerStack for the SOI process\n", + "stack = LayerStack(pdk_name=\"ubcpdk\")\n", + "\n", + "stack.layers[\"core\"] = Layer(\n", + " name=\"core\",\n", + " gds_layer=(1, 0),\n", + " zmin=0.0,\n", + " zmax=0.22,\n", + " thickness=0.22,\n", + " material=\"si\",\n", + " layer_type=\"dielectric\",\n", + ")\n", + "\n", + "stack.dielectrics = [\n", + " {\"name\": \"box\", \"material\": \"SiO2\", \"zmin\": -3.0, \"zmax\": 0.0},\n", + " {\"name\": \"clad\", \"material\": \"SiO2\", \"zmin\": 0.0, \"zmax\": 1.8},\n", + "]\n", + "\n", + "stack.materials = {\n", + " \"si\": {\"type\": \"dielectric\", \"permittivity\": 11.7},\n", + " \"SiO2\": {\"type\": \"dielectric\", \"permittivity\": 2.1},\n", + "}\n", + "\n", + "# Configure simulation\n", + "sim = MeepMeshSim()\n", + "sim.geometry(component=c, stack=stack)\n", + "sim.materials = {\"si\": 3.47, \"SiO2\": 1.44}\n", + "sim.source(port=\"o1\", wavelength=1.55, wavelength_span=0.1, num_freqs=11)\n", + "sim.monitors = [\"o1\", \"o2\", \"o3\"]\n", + "sim.domain(pml=1.0, margin=0.5)\n", + "sim.solver(resolution=32)" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "### Generate mesh and write config" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "result = sim.mesh(dirname, refined_mesh_size=0.2, max_mesh_size=1.0)\n", + "config_path = sim.write_config()\n", + "\n", + "print(f\"Mesh file: {result.mesh_path}\")\n", + "print(f\"File size: {result.mesh_path.stat().st_size / 1024:.0f} KB\")\n", + "print(f\"Config: {config_path}\")" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "### Mesh summary" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "import meshio\n", + "import numpy as np\n", + "\n", + "stats = result.mesh_stats\n", + "groups = result.groups\n", + "m = meshio.read(result.mesh_path)\n", + "tag_to_name = {tag: name for name, (tag, _) in m.field_data.items()}\n", + "\n", + "# --- Header ---\n", + "size_kb = result.mesh_path.stat().st_size / 1024\n", + "print(f\"Mesh: {result.mesh_path.name} ({size_kb:.0f} KB)\")\n", + "print(f\"Nodes: {stats.get('nodes', '?'):,} Tets: {stats.get('tetrahedra', '?'):,}\")\n", + "print()\n", + "\n", + "# --- Quality ---\n", + "q = stats.get(\"quality\", {})\n", + "sicn = stats.get(\"sicn\", {})\n", + "edge = stats.get(\"edge_length\", {})\n", + "print(\"Quality\")\n", + "print(\n", + " f\" Shape (gamma): min={q.get('min', '?')} mean={q.get('mean', '?')} max={q.get('max', '?')}\"\n", + ")\n", + "print(\n", + " f\" SICN: min={sicn.get('min', '?')} mean={sicn.get('mean', '?')} invalid={sicn.get('invalid', '?')}\"\n", + ")\n", + "print(f\" Edge length: min={edge.get('min', '?')} max={edge.get('max', '?')}\")\n", + "\n", + "# Edge ratio from actual mesh\n", + "for cells in m.cells:\n", + " if cells.type == \"tetra\":\n", + " pts = m.points[cells.data]\n", + " pairs = [(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)]\n", + " all_edges = np.stack(\n", + " [np.linalg.norm(pts[:, i] - pts[:, j], axis=1) for i, j in pairs]\n", + " )\n", + " ratios = all_edges.max(axis=0) / np.maximum(all_edges.min(axis=0), 1e-15)\n", + " bad = int(np.sum(ratios > 20))\n", + " print(\n", + " f\" Edge ratio: mean={ratios.mean():.1f} max={ratios.max():.1f} bad(>20)={bad}\"\n", + " )\n", + "print()\n", + "\n", + "# --- Bounding box ---\n", + "bb = stats.get(\"bbox\", {})\n", + "print(f\"Bounding box\")\n", + "print(\n", + " f\" x: [{bb.get('xmin', 0):.2f}, {bb.get('xmax', 0):.2f}] ({bb.get('xmax', 0) - bb.get('xmin', 0):.2f} um)\"\n", + ")\n", + "print(\n", + " f\" y: [{bb.get('ymin', 0):.2f}, {bb.get('ymax', 0):.2f}] ({bb.get('ymax', 0) - bb.get('ymin', 0):.2f} um)\"\n", + ")\n", + "print(\n", + " f\" z: [{bb.get('zmin', 0):.2f}, {bb.get('zmax', 0):.2f}] ({bb.get('zmax', 0) - bb.get('zmin', 0):.2f} um)\"\n", + ")\n", + "print()\n", + "\n", + "# --- Physical groups ---\n", + "print(\"Physical groups\")\n", + "for cells, phys in zip(m.cells, m.cell_data[\"gmsh:physical\"]):\n", + " if cells.type != \"tetra\":\n", + " continue\n", + " for tag, name in sorted(tag_to_name.items(), key=lambda x: x[1]):\n", + " mask = phys == tag\n", + " n = int(np.sum(mask))\n", + " if n == 0:\n", + " continue\n", + " pts = m.points[cells.data[mask].ravel()]\n", + " zlo, zhi = pts[:, 2].min(), pts[:, 2].max()\n", + " mat = \"\"\n", + " if name in groups.get(\"layer_volumes\", {}):\n", + " mat = f\" (material: {groups['layer_volumes'][name].get('material', '?')})\"\n", + " print(f\" {name:16s} {n:>7,} tets z=[{zlo:.3f}, {zhi:.3f}]{mat}\")\n", + "\n", + "if groups.get(\"port_surfaces\"):\n", + " print()\n", + " print(\"Port surfaces\")\n", + " for name, info in groups[\"port_surfaces\"].items():\n", + " c = info[\"center\"]\n", + " print(\n", + " f\" {name:16s} center=({c[0]:.2f}, {c[1]:.2f}) width={info['width']:.2f} z=[{info['z_range'][0]:.3f}, {info['z_range'][1]:.3f}]\"\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "### Visualize the mesh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "# Full mesh wireframe\n", + "sim.plot_mesh(interactive=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "# Show only the waveguide core\n", + "sim.plot_mesh(show_groups=[\"core\", \"port_\"], interactive=False)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "gsim", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/nbs/mesh_ybranch_split_clad.ipynb b/nbs/mesh_ybranch_split_clad.ipynb new file mode 100644 index 0000000..4bba18a --- /dev/null +++ b/nbs/mesh_ybranch_split_clad.ipynb @@ -0,0 +1,268 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# GMSH Mesh — Separate Box and Clad Groups\n", + "\n", + "Same as `mesh_ybranch.ipynb` but uses different material names for\n", + "box and clad so they appear as separate physical groups in the mesh." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "### Load component from UBC PDK" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "from ubcpdk import PDK, cells\n", + "\n", + "PDK.activate()\n", + "\n", + "# c = cells.ebeam_y_1550()\n", + "# c = cells.ring_single()\n", + "c = cells.coupler()\n", + "\n", + "dirname = \"./sim-data-mesh-----coupler\"\n", + "\n", + "c" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "### Configure MeepMeshSim" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "from gsim.common.stack.extractor import Layer, LayerStack\n", + "from gsim.meep import MeepMeshSim\n", + "\n", + "# Build a LayerStack for the SOI process\n", + "stack = LayerStack(pdk_name=\"ubcpdk\")\n", + "\n", + "stack.layers[\"core\"] = Layer(\n", + " name=\"core\",\n", + " gds_layer=(1, 0),\n", + " zmin=0.0,\n", + " zmax=0.22,\n", + " thickness=0.22,\n", + " material=\"si\",\n", + " layer_type=\"dielectric\",\n", + ")\n", + "\n", + "# Use different material names for box vs clad so they get separate\n", + "# physical groups in the mesh (both are SiO2 but named differently).\n", + "stack.dielectrics = [\n", + " {\"name\": \"box\", \"material\": \"SiO2_box\", \"zmin\": -3.0, \"zmax\": 0.0},\n", + " {\"name\": \"clad\", \"material\": \"SiO2_clad\", \"zmin\": 0.0, \"zmax\": 1.8},\n", + "]\n", + "\n", + "stack.materials = {\n", + " \"si\": {\"type\": \"dielectric\", \"permittivity\": 11.7},\n", + " \"SiO2_box\": {\"type\": \"dielectric\", \"permittivity\": 2.1},\n", + " \"SiO2_clad\": {\"type\": \"dielectric\", \"permittivity\": 2.1},\n", + "}\n", + "\n", + "# Configure simulation\n", + "sim = MeepMeshSim()\n", + "sim.geometry(component=c, stack=stack)\n", + "sim.materials = {\"si\": 3.47, \"SiO2_box\": 1.44, \"SiO2_clad\": 1.44}\n", + "sim.source(port=\"o1\", wavelength=1.55, wavelength_span=0.1, num_freqs=11)\n", + "sim.monitors = [\"o1\", \"o2\", \"o3\"]\n", + "sim.domain(pml=1.0, margin=0.5)\n", + "sim.solver(resolution=32)" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "### Generate mesh and write config" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "result = sim.mesh(dirname, refined_mesh_size=0.2, max_mesh_size=1.0)\n", + "config_path = sim.write_config()\n", + "\n", + "print(f\"Mesh file: {result.mesh_path}\")\n", + "print(f\"File size: {result.mesh_path.stat().st_size / 1024:.0f} KB\")\n", + "print(f\"Config: {config_path}\")" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "### Mesh summary" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "import meshio\n", + "import numpy as np\n", + "\n", + "stats = result.mesh_stats\n", + "groups = result.groups\n", + "m = meshio.read(result.mesh_path)\n", + "tag_to_name = {tag: name for name, (tag, _) in m.field_data.items()}\n", + "\n", + "# --- Header ---\n", + "size_kb = result.mesh_path.stat().st_size / 1024\n", + "print(f\"Mesh: {result.mesh_path.name} ({size_kb:.0f} KB)\")\n", + "print(f\"Nodes: {stats.get('nodes', '?'):,} Tets: {stats.get('tetrahedra', '?'):,}\")\n", + "print()\n", + "\n", + "# --- Quality ---\n", + "q = stats.get(\"quality\", {})\n", + "sicn = stats.get(\"sicn\", {})\n", + "edge = stats.get(\"edge_length\", {})\n", + "print(\"Quality\")\n", + "print(\n", + " f\" Shape (gamma): min={q.get('min', '?')} mean={q.get('mean', '?')} max={q.get('max', '?')}\"\n", + ")\n", + "print(\n", + " f\" SICN: min={sicn.get('min', '?')} mean={sicn.get('mean', '?')} invalid={sicn.get('invalid', '?')}\"\n", + ")\n", + "print(f\" Edge length: min={edge.get('min', '?')} max={edge.get('max', '?')}\")\n", + "\n", + "# Edge ratio from actual mesh\n", + "for cells in m.cells:\n", + " if cells.type == \"tetra\":\n", + " pts = m.points[cells.data]\n", + " pairs = [(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)]\n", + " all_edges = np.stack(\n", + " [np.linalg.norm(pts[:, i] - pts[:, j], axis=1) for i, j in pairs]\n", + " )\n", + " ratios = all_edges.max(axis=0) / np.maximum(all_edges.min(axis=0), 1e-15)\n", + " bad = int(np.sum(ratios > 20))\n", + " print(\n", + " f\" Edge ratio: mean={ratios.mean():.1f} max={ratios.max():.1f} bad(>20)={bad}\"\n", + " )\n", + "print()\n", + "\n", + "# --- Bounding box ---\n", + "bb = stats.get(\"bbox\", {})\n", + "print(f\"Bounding box\")\n", + "print(\n", + " f\" x: [{bb.get('xmin', 0):.2f}, {bb.get('xmax', 0):.2f}] ({bb.get('xmax', 0) - bb.get('xmin', 0):.2f} um)\"\n", + ")\n", + "print(\n", + " f\" y: [{bb.get('ymin', 0):.2f}, {bb.get('ymax', 0):.2f}] ({bb.get('ymax', 0) - bb.get('ymin', 0):.2f} um)\"\n", + ")\n", + "print(\n", + " f\" z: [{bb.get('zmin', 0):.2f}, {bb.get('zmax', 0):.2f}] ({bb.get('zmax', 0) - bb.get('zmin', 0):.2f} um)\"\n", + ")\n", + "print()\n", + "\n", + "# --- Physical groups ---\n", + "print(\"Physical groups\")\n", + "for cells, phys in zip(m.cells, m.cell_data[\"gmsh:physical\"]):\n", + " if cells.type != \"tetra\":\n", + " continue\n", + " for tag, name in sorted(tag_to_name.items(), key=lambda x: x[1]):\n", + " mask = phys == tag\n", + " n = int(np.sum(mask))\n", + " if n == 0:\n", + " continue\n", + " pts = m.points[cells.data[mask].ravel()]\n", + " zlo, zhi = pts[:, 2].min(), pts[:, 2].max()\n", + " mat = \"\"\n", + " if name in groups.get(\"layer_volumes\", {}):\n", + " mat = f\" (material: {groups['layer_volumes'][name].get('material', '?')})\"\n", + " print(f\" {name:16s} {n:>7,} tets z=[{zlo:.3f}, {zhi:.3f}]{mat}\")\n", + "\n", + "if groups.get(\"port_surfaces\"):\n", + " print()\n", + " print(\"Port surfaces\")\n", + " for name, info in groups[\"port_surfaces\"].items():\n", + " c = info[\"center\"]\n", + " print(\n", + " f\" {name:16s} center=({c[0]:.2f}, {c[1]:.2f}) width={info['width']:.2f} z=[{info['z_range'][0]:.3f}, {info['z_range'][1]:.3f}]\"\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "### Visualize the mesh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "# Full mesh wireframe\n", + "sim.plot_mesh(interactive=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "# Show only the waveguide core\n", + "sim.plot_mesh(show_groups=[\"core\", \"port_\"], interactive=False)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "gsim", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/gsim/common/mesh/__init__.py b/src/gsim/common/mesh/__init__.py new file mode 100644 index 0000000..024116f --- /dev/null +++ b/src/gsim/common/mesh/__init__.py @@ -0,0 +1,36 @@ +"""Generic GMSH mesh generation for any LayerStack. + +Provides solver-agnostic mesh generation from gdsfactory components +and LayerStack definitions. Used by Palace internally; also available +for standalone mesh inspection or other solvers. + +Usage: + from gsim.common.mesh import generate_mesh, GeometryData, MeshResult + + result = generate_mesh(component, stack, output_dir="./mesh_output") +""" + +from __future__ import annotations + +from gsim.common.mesh import gmsh_utils +from gsim.common.mesh.generator import collect_mesh_stats, generate_mesh +from gsim.common.mesh.geometry import ( + GeometryData, + add_dielectrics, + add_layer_volumes, + extract_geometry, + get_layer_info, +) +from gsim.common.mesh.types import MeshResult + +__all__ = [ + "GeometryData", + "MeshResult", + "add_dielectrics", + "add_layer_volumes", + "collect_mesh_stats", + "extract_geometry", + "generate_mesh", + "get_layer_info", + "gmsh_utils", +] diff --git a/src/gsim/common/mesh/generator.py b/src/gsim/common/mesh/generator.py new file mode 100644 index 0000000..4694dca --- /dev/null +++ b/src/gsim/common/mesh/generator.py @@ -0,0 +1,504 @@ +"""Generic GMSH mesh generator for any LayerStack. + +Orchestrates geometry extraction, volume creation, dielectric boxes, +fragmentation, physical group assignment, refinement, and meshing. +Solver-agnostic — produces a .msh file usable by Palace, other solvers, +or standalone mesh inspection. +""" + +from __future__ import annotations + +import logging +from pathlib import Path +from typing import TYPE_CHECKING + +import gmsh + +from gsim.common.mesh import gmsh_utils +from gsim.common.mesh.geometry import ( + add_dielectrics, + add_layer_volumes, + extract_geometry, +) +from gsim.common.mesh.types import MeshResult + +if TYPE_CHECKING: + from gsim.common.stack import LayerStack + +logger = logging.getLogger(__name__) + + +def generate_mesh( + component, + stack: LayerStack, + output_dir: str | Path, + *, + model_name: str = "mesh", + refined_mesh_size: float = 5.0, + max_mesh_size: float = 300.0, + margin: float = 50.0, + air_margin: float = 50.0, + include_airbox: bool = True, + include_ports: bool = True, + mesh_scale: float | None = None, + show_gui: bool = False, +) -> MeshResult: + """Generate a generic GMSH mesh from a gdsfactory component and LayerStack. + + Steps: + 1. Extract polygon geometry from the component + 2. Extrude all layer polygons into 3D volumes + 3. Add dielectric background boxes and optionally airbox + 3b. (Optional) Add port surfaces from component ports + 4. Fragment all geometry for conformal meshing + 5. Assign physical groups (one per material, plus outer boundary) + 5b. Assign port physical groups + 6. Set up mesh refinement near layer volume boundaries + 7. Generate 3D mesh and write .msh file + + Args: + component: gdsfactory Component + stack: LayerStack with layer and material definitions + output_dir: Directory for output files + model_name: Base name for output files (default: "mesh") + refined_mesh_size: Fine mesh size near layer boundaries (um) + max_mesh_size: Coarse mesh size in air/dielectric (um) + margin: XY margin around design (um) + air_margin: Air box margin around dielectric envelope (um) + include_airbox: Whether to add a surrounding airbox volume. + Needed for RF (Palace); not needed for photonics (Meep). + include_ports: Whether to create port surfaces from component + ports as physical groups. Default ``True``. + mesh_scale: Scale factor applied to mesh coordinates after generation. + E.g. ``1000.0`` converts um → nm. Default ``None`` (no scaling). + show_gui: Show gmsh GUI during meshing + + Returns: + MeshResult with mesh path, statistics, and physical group info + """ + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + msh_path = output_dir / f"{model_name}.msh" + + # 1. Extract geometry + logger.info("Extracting geometry...") + geometry = extract_geometry(component, stack) + logger.info(" Polygons: %s", len(geometry.polygons)) + logger.info(" Bbox: %s", geometry.bbox) + + # Initialize gmsh + gmsh.initialize() + gmsh.option.setNumber("General.Verbosity", 3) + + model_tag = f"{model_name}_mesh" + if model_tag in gmsh.model.list(): + gmsh.model.setCurrent(model_tag) + gmsh.model.remove() + gmsh.model.add(model_tag) + + kernel = gmsh.model.occ + + try: + # 2. Extrude all layer polygons into volumes + logger.info("Adding layer volumes...") + layer_volume_tags = add_layer_volumes(kernel, geometry, stack) + + # 3. Add dielectric boxes (and optionally airbox) + logger.info("Adding dielectrics...") + dielectric_tags = add_dielectrics( + kernel, + geometry, + stack, + margin, + air_margin, + include_airbox=include_airbox, + ) + + # 4. Fragment all geometry + logger.info("Fragmenting geometry...") + geom_dimtags, geom_map = gmsh_utils.fragment_all(kernel) + + # 5. Assign physical groups + logger.info("Assigning physical groups...") + groups = _assign_generic_groups( + kernel, + layer_volume_tags, + dielectric_tags, + geom_dimtags, + geom_map, + stack, + ) + + # 5b. Find and tag port surfaces on existing volume boundaries + if include_ports: + logger.info("Identifying port surfaces...") + _assign_port_groups(component, stack, groups) + + # 6. Mesh refinement + logger.info("Setting up mesh refinement...") + _setup_generic_refinement(kernel, groups, refined_mesh_size, max_mesh_size) + + # Show GUI if requested + if show_gui: + gmsh.fltk.run() + + # 7. Generate mesh and write + logger.info("Generating 3D mesh...") + gmsh.model.mesh.generate(3) + + if mesh_scale is not None: + s = mesh_scale + gmsh.model.mesh.affineTransform([s, 0, 0, 0, 0, s, 0, 0, 0, 0, s, 0]) + + mesh_stats = collect_mesh_stats() + + gmsh.option.setNumber("Mesh.Binary", 0) + gmsh.option.setNumber("Mesh.SaveAll", 0) + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + gmsh.write(str(msh_path)) + + logger.info("Mesh saved: %s", msh_path) + + finally: + gmsh.clear() + gmsh.finalize() + + return MeshResult( + mesh_path=msh_path, + mesh_stats=mesh_stats, + groups=groups, + ) + + +def _assign_generic_groups( + kernel, + layer_volume_tags: dict[str, list[int]], + dielectric_tags: dict[str, list[int]], + geom_dimtags: list, + geom_map: list, + stack: LayerStack, +) -> dict: + """Assign physical groups for the generic mesh. + + Creates one volume physical group per material name (merging layers + that share the same material) plus an outer boundary surface group + from the airbox. + + Returns: + {"volumes": {material: {"phys_group": int, "tags": [int]}}, + "layer_volumes": {layer_name: {"phys_group": int, "tags": [int]}}, + "outer_boundary": {"phys_group": int, "tags": [int]}} + """ + groups: dict = { + "volumes": {}, + "layer_volumes": {}, + "outer_boundary": {}, + } + + # --- Layer volume groups first (to know which tags to exclude from dielectrics) --- + layer_volume_tag_set: set[int] = set() + for layer_name, tags in layer_volume_tags.items(): + new_tags = gmsh_utils.get_tags_after_fragment( + tags, geom_dimtags, geom_map, dimension=3 + ) + if new_tags: + layer_volume_tag_set.update(new_tags) + # Look up material for this layer to also merge into volume groups + layer = stack.layers.get(layer_name) + material_name = layer.material if layer else layer_name + + phys_group = gmsh_utils.assign_physical_group(3, new_tags, layer_name) + groups["layer_volumes"][layer_name] = { + "phys_group": phys_group, + "tags": new_tags, + "material": material_name, + } + + # --- Dielectric / airbox volume groups (exclude layer volume fragments) --- + for material_name, tags in dielectric_tags.items(): + new_tags = gmsh_utils.get_tags_after_fragment( + tags, geom_dimtags, geom_map, dimension=3 + ) + if new_tags: + # Exclude fragments that belong to layer volumes to avoid + # double-assignment (e.g. waveguide core inside an oxide box) + new_tags = [t for t in new_tags if t not in layer_volume_tag_set] + if new_tags: + phys_group = gmsh_utils.assign_physical_group(3, new_tags, material_name) + groups["volumes"][material_name] = { + "phys_group": phys_group, + "tags": new_tags, + } + + # --- Outer boundary from airbox --- + if "airbox" in groups["volumes"]: + airbox_tags = groups["volumes"]["airbox"]["tags"] + if airbox_tags: + try: + _, simulation_boundary = kernel.getSurfaceLoops(airbox_tags[0]) + if simulation_boundary: + boundary_tags = list(next(iter(simulation_boundary))) + phys_group = gmsh_utils.assign_physical_group( + 2, boundary_tags, "outer_boundary" + ) + groups["outer_boundary"] = { + "phys_group": phys_group, + "tags": boundary_tags, + } + except Exception: + logger.warning("Could not extract outer boundary from airbox") + + kernel.synchronize() + return groups + + +def _assign_port_groups( + component, + stack: LayerStack, + groups: dict, +) -> None: + """Find existing volume boundary faces at port locations. + + Assign them as physical groups. + + After fragmentation, the layer volumes have boundary faces at each + component port. This function identifies those faces by matching + their bounding box to the expected port rectangle, then assigns + them to named surface physical groups (``port_``). + + Modifies *groups* in-place, adding a ``"port_surfaces"`` key. + """ + from gsim.common.mesh.geometry import _resolve_port_layer + + # Collect all boundary surfaces of post-fragment layer volumes + layer_boundary_surfs: set[int] = set() + for layer_info in groups["layer_volumes"].values(): + for vtag in layer_info["tags"]: + try: + boundary = gmsh.model.getBoundary([(3, vtag)], oriented=False) + for _, stag in boundary: + layer_boundary_surfs.add(stag) + except Exception: + pass + + if not layer_boundary_surfs: + return + + # Cache bounding boxes for all boundary surfaces + surf_bboxes: dict[int, tuple[float, ...]] = {} + import contextlib + + for stag in layer_boundary_surfs: + with contextlib.suppress(Exception): + surf_bboxes[stag] = gmsh.model.getBoundingBox(2, stag) + + tol = 0.02 # bounding-box match tolerance (um) + + groups["port_surfaces"] = {} + + for port in component.ports: + resolved = _resolve_port_layer(port, component, stack) + if resolved is None: + continue + layer_name, zmin, zmax = resolved + + cx, cy = float(port.center[0]), float(port.center[1]) + hw = float(port.width) / 2 + angle = float(port.orientation) % 360 + + # Expected port rectangle bounds + if angle < 45 or angle >= 315 or (135 <= angle < 225): + # East/West → YZ plane at x=cx + exp = (cx, cy - hw, zmin, cx, cy + hw, zmax) + else: + # North/South → XZ plane at y=cy + exp = (cx - hw, cy, zmin, cx + hw, cy, zmax) + + matching_tags: list[int] = [] + for stag, bbox in surf_bboxes.items(): + sxmin, symin, szmin, sxmax, symax, szmax = bbox + if ( + abs(sxmin - exp[0]) < tol + and abs(symin - exp[1]) < tol + and abs(szmin - exp[2]) < tol + and abs(sxmax - exp[3]) < tol + and abs(symax - exp[4]) < tol + and abs(szmax - exp[5]) < tol + ): + matching_tags.append(stag) + + if matching_tags: + phys_group = gmsh_utils.assign_physical_group( + 2, matching_tags, f"port_{port.name}" + ) + groups["port_surfaces"][port.name] = { + "phys_group": phys_group, + "tags": matching_tags, + "center": [cx, cy], + "width": float(port.width), + "orientation": angle, + "layer": layer_name, + "z_range": [zmin, zmax], + } + logger.info( + "Port '%s': phys_group=%d (%d face(s)) on layer '%s'", + port.name, + phys_group, + len(matching_tags), + layer_name, + ) + else: + logger.warning( + "Port '%s': no matching boundary face found at center=(%.3f, %.3f)", + port.name, + cx, + cy, + ) + + gmsh.model.occ.synchronize() + + +def _setup_generic_refinement( + kernel, + groups: dict, + refined_cellsize: float, + max_cellsize: float, +) -> None: + """Set up mesh refinement near layer volume boundaries.""" + boundary_lines: list[int] = [] + + # Collect boundary lines from all layer volumes + for layer_info in groups["layer_volumes"].values(): + for vol_tag in layer_info["tags"]: + # Get surfaces bounding this volume + surfs = gmsh.model.getBoundary([(3, vol_tag)], oriented=False) + for _dim, surf_tag in surfs: + try: + lines = gmsh_utils.get_boundary_lines(surf_tag, kernel) + boundary_lines.extend(lines) + except Exception: + pass + + field_ids: list[int] = [] + if boundary_lines: + field_id = gmsh_utils.setup_mesh_refinement( + boundary_lines, refined_cellsize, max_cellsize + ) + field_ids.append(field_id) + + if field_ids: + gmsh_utils.finalize_mesh_fields(field_ids) + + +def collect_mesh_stats() -> dict: + """Collect mesh statistics from gmsh after mesh generation. + + Must be called while gmsh is initialized and the mesh is generated. + + Returns: + Dict with mesh statistics including: + - bbox: Bounding box coordinates + - nodes: Number of nodes + - elements: Total element count + - tetrahedra: Tet count + - quality: Shape quality metrics (gamma) + - sicn: Signed Inverse Condition Number + - edge_length: Min/max edge lengths + - groups: Physical group info + """ + stats: dict = {} + + # Get bounding box + try: + xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.getBoundingBox(-1, -1) + stats["bbox"] = { + "xmin": xmin, + "ymin": ymin, + "zmin": zmin, + "xmax": xmax, + "ymax": ymax, + "zmax": zmax, + } + except Exception: + pass + + # Get node count + try: + node_tags, _, _ = gmsh.model.mesh.getNodes() + stats["nodes"] = len(node_tags) + except Exception: + pass + + # Get element counts and collect tet tags for quality + tet_tags = [] + try: + element_types, element_tags, _ = gmsh.model.mesh.getElements() + total_elements = sum(len(tags) for tags in element_tags) + stats["elements"] = total_elements + + # Count tetrahedra (type 4) and save tags + for etype, tags in zip(element_types, element_tags, strict=False): + if etype == 4: # 4-node tetrahedron + stats["tetrahedra"] = len(tags) + tet_tags = list(tags) + except Exception: + pass + + # Get mesh quality for tetrahedra + if tet_tags: + # Gamma: inscribed/circumscribed radius ratio (shape quality) + try: + qualities = gmsh.model.mesh.getElementQualities(tet_tags, "gamma") + if len(qualities) > 0: + stats["quality"] = { + "min": round(min(qualities), 3), + "max": round(max(qualities), 3), + "mean": round(sum(qualities) / len(qualities), 3), + } + except Exception: + pass + + # SICN: Signed Inverse Condition Number (negative = invalid element) + try: + sicn = gmsh.model.mesh.getElementQualities(tet_tags, "minSICN") + if len(sicn) > 0: + sicn_min = min(sicn) + invalid_count = sum(1 for s in sicn if s < 0) + stats["sicn"] = { + "min": round(sicn_min, 3), + "mean": round(sum(sicn) / len(sicn), 3), + "invalid": invalid_count, + } + except Exception: + pass + + # Edge lengths + try: + min_edges = gmsh.model.mesh.getElementQualities(tet_tags, "minEdge") + max_edges = gmsh.model.mesh.getElementQualities(tet_tags, "maxEdge") + if len(min_edges) > 0 and len(max_edges) > 0: + stats["edge_length"] = { + "min": round(min(min_edges), 3), + "max": round(max(max_edges), 3), + } + except Exception: + pass + + # Get physical groups with tags + try: + phys_groups: dict[str, list] = {"volumes": [], "surfaces": []} + for dim, tag in gmsh.model.getPhysicalGroups(): + name = gmsh.model.getPhysicalName(dim, tag) + entry = {"name": name, "tag": tag} + if dim == 3: + phys_groups["volumes"].append(entry) + elif dim == 2: + phys_groups["surfaces"].append(entry) + stats["groups"] = phys_groups + except Exception: + pass + + return stats + + +__all__ = ["collect_mesh_stats", "generate_mesh"] diff --git a/src/gsim/common/mesh/geometry.py b/src/gsim/common/mesh/geometry.py new file mode 100644 index 0000000..2e8d194 --- /dev/null +++ b/src/gsim/common/mesh/geometry.py @@ -0,0 +1,351 @@ +"""Geometry extraction and creation for generic GMSH mesh generation. + +This module handles extracting polygons from gdsfactory components +and creating 3D geometry in gmsh. Solver-agnostic — no Palace or Meep imports. +""" + +from __future__ import annotations + +import logging +import math +from dataclasses import dataclass +from typing import TYPE_CHECKING + +from gsim.common.mesh import gmsh_utils + +if TYPE_CHECKING: + from gsim.common.stack import LayerStack + +logger = logging.getLogger(__name__) + + +@dataclass +class GeometryData: + """Container for geometry data extracted from component.""" + + polygons: list # List of (layer_num, pts_x, pts_y, holes) tuples + bbox: tuple[float, float, float, float] # (xmin, ymin, xmax, ymax) + layer_bboxes: dict # layer_num -> (xmin, ymin, xmax, ymax) + + +def extract_geometry(component, stack: LayerStack) -> GeometryData: + """Extract polygon geometry from a gdsfactory component. + + Args: + component: gdsfactory Component + stack: LayerStack for layer mapping + + Returns: + GeometryData with polygons and bounding boxes + """ + polygons = [] + global_bbox = [math.inf, math.inf, -math.inf, -math.inf] + layer_bboxes = {} + + # Get polygons from component + polygons_by_index = component.get_polygons() + + # Build layer_index -> GDS tuple mapping + layout = component.kcl.layout + index_to_gds = {} + for layer_index in range(layout.layers()): + if layout.is_valid_layer(layer_index): + info = layout.get_info(layer_index) + index_to_gds[layer_index] = (info.layer, info.datatype) + + # Build GDS tuple -> layer number mapping + gds_to_layernum = {} + for layer_data in stack.layers.values(): + gds_tuple = tuple(layer_data.gds_layer) + gds_to_layernum[gds_tuple] = gds_tuple[0] + + # Convert polygons + for layer_index, polys in polygons_by_index.items(): + gds_tuple = index_to_gds.get(layer_index) + if gds_tuple is None: + continue + + layernum = gds_to_layernum.get(gds_tuple) + if layernum is None: + continue + + for poly in polys: + # Convert klayout polygon to lists (nm -> um) + points = list(poly.each_point_hull()) + if len(points) < 3: + continue + + pts_x = [pt.x / 1000.0 for pt in points] + pts_y = [pt.y / 1000.0 for pt in points] + + # Extract holes from polygon + holes = [] + for hole_idx in range(poly.holes()): + hole_pts = list(poly.each_point_hole(hole_idx)) + if len(hole_pts) >= 3: + hx = [pt.x / 1000.0 for pt in hole_pts] + hy = [pt.y / 1000.0 for pt in hole_pts] + holes.append((hx, hy)) + + polygons.append((layernum, pts_x, pts_y, holes)) + + # Update bounding boxes + xmin, xmax = min(pts_x), max(pts_x) + ymin, ymax = min(pts_y), max(pts_y) + + global_bbox[0] = min(global_bbox[0], xmin) + global_bbox[1] = min(global_bbox[1], ymin) + global_bbox[2] = max(global_bbox[2], xmax) + global_bbox[3] = max(global_bbox[3], ymax) + + if layernum not in layer_bboxes: + layer_bboxes[layernum] = [xmin, ymin, xmax, ymax] + else: + bbox = layer_bboxes[layernum] + bbox[0] = min(bbox[0], xmin) + bbox[1] = min(bbox[1], ymin) + bbox[2] = max(bbox[2], xmax) + bbox[3] = max(bbox[3], ymax) + + return GeometryData( + polygons=polygons, + bbox=(global_bbox[0], global_bbox[1], global_bbox[2], global_bbox[3]), + layer_bboxes=layer_bboxes, + ) + + +def get_layer_info(stack: LayerStack, gds_layer: int) -> dict | None: + """Get layer info from stack by GDS layer number. + + Args: + stack: LayerStack with layer definitions + gds_layer: GDS layer number + + Returns: + Dict with layer info or None if not found + """ + for name, layer in stack.layers.items(): + if layer.gds_layer[0] == gds_layer: + return { + "name": name, + "zmin": layer.zmin, + "zmax": layer.zmax, + "thickness": layer.zmax - layer.zmin, + "material": layer.material, + "type": layer.layer_type, + } + return None + + +def add_layer_volumes( + kernel, + geometry: GeometryData, + stack: LayerStack, +) -> dict[str, list[int]]: + """Extrude ALL layer polygons into 3D volumes. + + Unlike Palace's ``add_metals()`` which only handles conductor/via layers + and creates shell surfaces for conductors, this function processes every + layer type and always creates true volumes. + + Args: + kernel: gmsh OCC kernel + geometry: Extracted geometry data + stack: LayerStack with layer definitions + + Returns: + Dict mapping layer_name to list of volume tags. + """ + layer_volumes: dict[str, list[int]] = {} + + # Group polygons by GDS layer number + polygons_by_layer: dict[int, list[tuple]] = {} + for layernum, pts_x, pts_y, holes in geometry.polygons: + if layernum not in polygons_by_layer: + polygons_by_layer[layernum] = [] + polygons_by_layer[layernum].append((pts_x, pts_y, holes)) + + # Process each layer + for layernum, polys in polygons_by_layer.items(): + layer_info = get_layer_info(stack, layernum) + if layer_info is None: + continue + + layer_name = layer_info["name"] + zmin = layer_info["zmin"] + thickness = layer_info["thickness"] + + if thickness <= 0: + continue + + if layer_name not in layer_volumes: + layer_volumes[layer_name] = [] + + for pts_x, pts_y, holes in polys: + vol_tag = gmsh_utils.extrude_polygon( + kernel, pts_x, pts_y, zmin, thickness, holes=holes + ) + if vol_tag is not None: + layer_volumes[layer_name].append(vol_tag) + + kernel.removeAllDuplicates() + kernel.synchronize() + + return layer_volumes + + +def add_dielectrics( + kernel, + geometry: GeometryData, + stack: LayerStack, + margin: float, + air_margin: float, + *, + include_airbox: bool = True, +) -> dict: + """Add dielectric boxes and optionally an airbox to gmsh. + + Args: + kernel: gmsh OCC kernel + geometry: Extracted geometry data + stack: LayerStack with dielectric definitions + margin: XY margin around design (um) + air_margin: Air box margin (um) + include_airbox: Whether to add a surrounding airbox volume. + Required for RF (Palace) simulations; typically not needed + for photonic (Meep) simulations. + + Returns: + Dict with material_name -> list of volume_tags + """ + dielectric_tags: dict[str, list[int]] = {} + + # Get overall geometry bounds + xmin, ymin, xmax, ymax = geometry.bbox + xmin -= margin + ymin -= margin + xmax += margin + ymax += margin + + # Track overall z range + z_min_all = math.inf + z_max_all = -math.inf + + # Merge adjacent dielectrics of the same material into single boxes + # to avoid coincident faces and sliver tets at shared z-boundaries. + # E.g. box (SiO2, -3→0) + clad (SiO2, 0→1.8) → one box (SiO2, -3→1.8). + merged: dict[str, list[tuple[float, float]]] = {} + for dielectric in stack.dielectrics: + material = dielectric["material"] + d_zmin = dielectric["zmin"] + d_zmax = dielectric["zmax"] + if material not in merged: + merged[material] = [] + merged[material].append((d_zmin, d_zmax)) + + # For each material, sort ranges and merge overlapping/adjacent ones + for material, ranges in merged.items(): + ranges.sort() + combined: list[tuple[float, float]] = [ranges[0]] + for lo, hi in ranges[1:]: + prev_lo, prev_hi = combined[-1] + if lo <= prev_hi + 1e-6: # adjacent or overlapping + combined[-1] = (prev_lo, max(prev_hi, hi)) + else: + combined.append((lo, hi)) + merged[material] = combined + + # Create one GMSH box per merged z-range. + # Adjacent boxes of different materials share a face at the same z; + # removeAllDuplicates() after creation merges these shared faces so + # fragment() works cleanly without sliver tets. + for material, ranges in merged.items(): + dielectric_tags[material] = [] + for d_zmin, d_zmax in ranges: + z_min_all = min(z_min_all, d_zmin) + z_max_all = max(z_max_all, d_zmax) + + box_tag = gmsh_utils.create_box( + kernel, + xmin, + ymin, + d_zmin, + xmax, + ymax, + d_zmax, + ) + dielectric_tags[material].append(box_tag) + + kernel.removeAllDuplicates() + kernel.synchronize() + + # Add surrounding airbox (needed for RF, not for photonics) + if include_airbox: + air_xmin = xmin - air_margin + air_ymin = ymin - air_margin + air_xmax = xmax + air_margin + air_ymax = ymax + air_margin + air_zmin = z_min_all - air_margin + air_zmax = z_max_all + air_margin + + airbox_tag = kernel.addBox( + air_xmin, + air_ymin, + air_zmin, + air_xmax - air_xmin, + air_ymax - air_ymin, + air_zmax - air_zmin, + ) + dielectric_tags["airbox"] = [airbox_tag] + + kernel.synchronize() + + return dielectric_tags + + +def _resolve_port_layer( + port, + component, + stack: LayerStack, +) -> tuple[str, float, float] | None: + """Resolve a gdsfactory port's layer to a stack layer name and z-range. + + Args: + port: gdsfactory Port object + component: gdsfactory Component (needed for klayout layer lookup) + stack: LayerStack + + Returns: + ``(layer_name, zmin, zmax)`` or ``None`` if no matching stack layer. + """ + # Build GDS tuple → stack layer lookup + gds_to_layer = {} + for name, layer in stack.layers.items(): + gds_to_layer[tuple(layer.gds_layer)] = (name, layer) + + # port.layer is a klayout layer index (int); convert to GDS tuple + port_layer_raw = port.layer + if isinstance(port_layer_raw, int): + layout = component.kcl.layout + if not layout.is_valid_layer(port_layer_raw): + return None + info = layout.get_info(port_layer_raw) + port_gds = (info.layer, info.datatype) + else: + port_gds = tuple(port_layer_raw) + + match = gds_to_layer.get(port_gds) + if match is None: + return None + + layer_name, layer = match + return (layer_name, layer.zmin, layer.zmax) + + +__all__ = [ + "GeometryData", + "add_dielectrics", + "add_layer_volumes", + "extract_geometry", + "get_layer_info", +] diff --git a/src/gsim/common/mesh/gmsh_utils.py b/src/gsim/common/mesh/gmsh_utils.py new file mode 100644 index 0000000..7658eff --- /dev/null +++ b/src/gsim/common/mesh/gmsh_utils.py @@ -0,0 +1,510 @@ +"""Generic GMSH utility functions for mesh generation. + +Provides low-level helpers for creating geometry, fragmenting, +assigning physical groups, and setting up mesh refinement in GMSH. +These are solver-agnostic — no Palace or Meep imports. +""" + +from __future__ import annotations + +import math + +import gmsh +import numpy as np + + +def create_box( + kernel, + xmin: float, + ymin: float, + zmin: float, + xmax: float, + ymax: float, + zmax: float, + meshseed: float = 0, +) -> int: + """Create a 3D box volume in gmsh. + + Args: + kernel: gmsh.model.occ kernel + xmin, ymin, zmin: minimum coordinates + xmax, ymax, zmax: maximum coordinates + meshseed: mesh seed size at corners (0 = auto) + + Returns: + Volume tag of created box + """ + if meshseed == 0: + # Use simple addBox + return kernel.addBox(xmin, ymin, zmin, xmax - xmin, ymax - ymin, zmax - zmin) + + # Create box with explicit mesh seed at corners + pt1 = kernel.addPoint(xmin, ymin, zmin, meshseed, -1) + pt2 = kernel.addPoint(xmin, ymax, zmin, meshseed, -1) + pt3 = kernel.addPoint(xmax, ymax, zmin, meshseed, -1) + pt4 = kernel.addPoint(xmax, ymin, zmin, meshseed, -1) + + line1 = kernel.addLine(pt1, pt2, -1) + line2 = kernel.addLine(pt2, pt3, -1) + line3 = kernel.addLine(pt3, pt4, -1) + line4 = kernel.addLine(pt4, pt1, -1) + linetaglist = [line1, line2, line3, line4] + + curvetag = kernel.addCurveLoop(linetaglist, tag=-1) + surfacetag = kernel.addPlaneSurface([curvetag], tag=-1) + returnval = kernel.extrude([(2, surfacetag)], 0, 0, zmax - zmin) + volumetag = returnval[1][1] + + return volumetag + + +def create_polygon_surface( + kernel, + pts_x: list[float], + pts_y: list[float], + z: float, + meshseed: float = 0, + holes: list[tuple[list[float], list[float]]] | None = None, +) -> int | None: + """Create a planar surface from polygon vertices at z height. + + Args: + kernel: gmsh.model.occ kernel + pts_x: list of x coordinates + pts_y: list of y coordinates + z: z coordinate of the surface + meshseed: mesh seed size at vertices (0 = auto) + holes: list of (hole_pts_x, hole_pts_y) tuples for interior holes + + Returns: + Surface tag, or None if polygon is invalid + """ + numvertices = len(pts_x) + if numvertices < 3: + return None + + linetaglist = [] + vertextaglist = [] + + # Create vertices + for v in range(numvertices): + vertextag = kernel.addPoint(pts_x[v], pts_y[v], z, meshseed, -1) + vertextaglist.append(vertextag) + + # Create lines connecting vertices + for v in range(numvertices): + pt_start = vertextaglist[v] + pt_end = vertextaglist[(v + 1) % numvertices] + try: + linetag = kernel.addLine(pt_start, pt_end, -1) + linetaglist.append(linetag) + except Exception: + pass # Skip degenerate lines + + if len(linetaglist) < 3: + return None + + # Create outer curve loop + curvetag = kernel.addCurveLoop(linetaglist, tag=-1) + + # Create hole curve loops + all_loops = [curvetag] + for hx, hy in holes or []: + hole_lines = [] + hole_verts = [] + for v in range(len(hx)): + vtag = kernel.addPoint(hx[v], hy[v], z, meshseed, -1) + hole_verts.append(vtag) + for v in range(len(hx)): + try: + ltag = kernel.addLine(hole_verts[v], hole_verts[(v + 1) % len(hx)], -1) + hole_lines.append(ltag) + except Exception: + pass + if len(hole_lines) >= 3: + hole_loop = kernel.addCurveLoop(hole_lines, tag=-1) + all_loops.append(hole_loop) + + surfacetag = kernel.addPlaneSurface(all_loops, tag=-1) + + return surfacetag + + +def extrude_polygon( + kernel, + pts_x: list[float], + pts_y: list[float], + zmin: float, + thickness: float, + meshseed: float = 0, + holes: list[tuple[list[float], list[float]]] | None = None, +) -> int | None: + """Create an extruded polygon volume (for vias, metals). + + Args: + kernel: gmsh.model.occ kernel + pts_x: list of x coordinates + pts_y: list of y coordinates + zmin: base z coordinate + thickness: extrusion height + meshseed: mesh seed size at vertices + holes: list of (hole_pts_x, hole_pts_y) tuples for interior holes + + Returns: + Volume tag if thickness > 0, surface tag if thickness == 0, or None if invalid + """ + surfacetag = create_polygon_surface( + kernel, pts_x, pts_y, zmin, meshseed, holes=holes + ) + if surfacetag is None: + return None + + if thickness > 0: + result = kernel.extrude([(2, surfacetag)], 0, 0, thickness) + # result[1] contains the volume (dim=3, tag) + return result[1][1] + + return surfacetag + + +def create_port_rectangle( + kernel, + xmin: float, + ymin: float, + zmin: float, + xmax: float, + ymax: float, + zmax: float, + meshseed: float = 0, +) -> int: + """Create a rectangular surface for a port. + + Handles both horizontal (z-plane) and vertical port surfaces. + + Args: + kernel: gmsh.model.occ kernel + xmin, ymin, zmin: minimum coordinates + xmax, ymax, zmax: maximum coordinates + meshseed: mesh seed size at corners + + Returns: + Surface tag of created port rectangle + """ + # Determine port orientation + dx = xmax - xmin + dz = zmax - zmin + + if dz < 1e-6: + # Horizontal port (in xy plane) + pt1 = kernel.addPoint(xmin, ymin, zmin, meshseed, -1) + pt2 = kernel.addPoint(xmin, ymax, zmin, meshseed, -1) + pt3 = kernel.addPoint(xmax, ymax, zmin, meshseed, -1) + pt4 = kernel.addPoint(xmax, ymin, zmin, meshseed, -1) + elif dx < 1e-6: + # Vertical port in yz plane + pt1 = kernel.addPoint(xmin, ymin, zmin, meshseed, -1) + pt2 = kernel.addPoint(xmin, ymax, zmin, meshseed, -1) + pt3 = kernel.addPoint(xmin, ymax, zmax, meshseed, -1) + pt4 = kernel.addPoint(xmin, ymin, zmax, meshseed, -1) + else: + # Vertical port in xz plane + pt1 = kernel.addPoint(xmin, ymin, zmin, meshseed, -1) + pt2 = kernel.addPoint(xmin, ymin, zmax, meshseed, -1) + pt3 = kernel.addPoint(xmax, ymin, zmax, meshseed, -1) + pt4 = kernel.addPoint(xmax, ymin, zmin, meshseed, -1) + + line1 = kernel.addLine(pt1, pt2, -1) + line2 = kernel.addLine(pt2, pt3, -1) + line3 = kernel.addLine(pt3, pt4, -1) + line4 = kernel.addLine(pt4, pt1, -1) + linetaglist = [line1, line2, line3, line4] + + curvetag = kernel.addCurveLoop(linetaglist, tag=-1) + surfacetag = kernel.addPlaneSurface([curvetag], tag=-1) + + return surfacetag + + +def fragment_all(kernel) -> tuple[list, list]: + """Fragment all geometry to ensure conformal mesh at intersections. + + Args: + kernel: gmsh.model.occ kernel + + Returns: + (geom_dimtags, geom_map) - original dimtags and mapping to new tags + """ + geom_dimtags = [x for x in kernel.getEntities() if x[0] in (2, 3)] + _, geom_map = kernel.fragment(geom_dimtags, []) + kernel.synchronize() + return geom_dimtags, geom_map + + +def get_tags_after_fragment( + original_tags: list[int], + geom_dimtags: list, + geom_map: list, + dimension: int = 2, +) -> list[int]: + """Get new tags after fragmenting, given original tags. + + Tags change after gmsh fragment operation. This function maps + original tags to their new values using the fragment mapping. + + Args: + original_tags: list of tags before fragmenting + geom_dimtags: list of all original dimtags before fragmenting + geom_map: mapping from fragment() function + dimension: dimension for tags (2=surfaces, 3=volumes) + + Returns: + List of new tags after fragmenting + """ + if isinstance(original_tags, int): + original_tags = [original_tags] + + indices = [ + i + for i, x in enumerate(geom_dimtags) + if x[0] == dimension and (x[1] in original_tags) + ] + raw = [geom_map[i] for i in indices] + flat = [item for sublist in raw for item in sublist] + newtags = [s[-1] for s in flat] + + return newtags + + +def assign_physical_group( + dim: int, + tags: list[int], + name: str, +) -> int: + """Assign tags to a physical group with a name. + + Args: + dim: dimension (2=surfaces, 3=volumes) + tags: list of entity tags + name: physical group name + + Returns: + Physical group tag + """ + if not tags: + return -1 + phys_group = gmsh.model.addPhysicalGroup(dim, tags, tag=-1) + gmsh.model.setPhysicalName(dim, phys_group, name) + return phys_group + + +def get_surface_normal(surface_tag: int) -> np.ndarray: + """Get the normal vector of a surface. + + Args: + surface_tag: surface entity tag + + Returns: + Normal vector as numpy array [nx, ny, nz] + """ + # Get the boundary of the surface + boundary_lines = gmsh.model.getBoundary([(2, surface_tag)], oriented=True) + + # Get points from these lines + points = [] + seen_points = set() + + for _dim, line_tag in boundary_lines: + line_points = gmsh.model.getBoundary([(1, line_tag)], oriented=True) + for _pdim, ptag in line_points: + if ptag not in seen_points: + coord = gmsh.model.getValue(0, ptag, []) + points.append(np.array(coord)) + seen_points.add(ptag) + if len(points) == 3: + break + if len(points) == 3: + break + + if len(points) < 3: + return np.array([0, 0, 1]) # Default to z-normal + + # Compute surface normal using cross product + v1 = points[1] - points[0] + v2 = points[2] - points[0] + normal = np.cross(v1, v2) + norm = np.linalg.norm(normal) + if norm > 0: + normal = normal / norm + return normal + + +def is_vertical_surface(surface_tag: int) -> bool: + """Check if a surface is vertical (not in xy plane). + + Args: + surface_tag: surface entity tag + + Returns: + True if surface is vertical (z component of normal is ~0) + """ + normal = get_surface_normal(surface_tag) + n = normal[2] + if not np.isnan(n): + return int(abs(n)) == 0 + return False + + +def get_volumes_at_z_range( + zmin: float, + zmax: float, + delta: float = 0.001, +) -> list[tuple[int, int]]: + """Get all volumes within a z-coordinate range. + + Args: + zmin: minimum z coordinate + zmax: maximum z coordinate + delta: tolerance for z comparison + + Returns: + List of (dim, tag) tuples for volumes in the z range + """ + volumes_in_bbox = gmsh.model.getEntitiesInBoundingBox( + -math.inf, + -math.inf, + zmin - delta / 2, + math.inf, + math.inf, + zmax + delta / 2, + 3, + ) + + volume_list = [] + for volume in volumes_in_bbox: + volume_tag = volume[1] + _, _, vzmin, _, _, vzmax = gmsh.model.getBoundingBox(3, volume_tag) + if ( + abs(vzmin - (zmin - delta / 2)) < delta + and abs(vzmax - (zmax + delta / 2)) < delta + ): + volume_list.append(volume) + + return volume_list + + +def get_surfaces_at_z(z: float, delta: float = 0.001) -> list[tuple[int, int]]: + """Get all surfaces at a specific z coordinate. + + Args: + z: z coordinate + delta: tolerance for z comparison + + Returns: + List of (dim, tag) tuples for surfaces at z + """ + return gmsh.model.getEntitiesInBoundingBox( + -math.inf, + -math.inf, + z - delta / 2, + math.inf, + math.inf, + z + delta / 2, + 2, + ) + + +def get_boundary_lines(surface_tag: int, kernel) -> list[int]: + """Get all boundary line tags of a surface. + + Args: + surface_tag: surface entity tag + kernel: gmsh.model.occ kernel + + Returns: + List of curve/line tags forming the surface boundary + """ + _clt, ct = kernel.getCurveLoops(surface_tag) + lines = [] + for curvetag in ct: + lines.extend(curvetag) + return lines + + +def setup_mesh_refinement( + boundary_line_tags: list[int], + refined_cellsize: float, + max_cellsize: float, +) -> int: + """Set up mesh refinement near boundary lines. + + Args: + boundary_line_tags: list of curve tags for refinement + refined_cellsize: mesh size near boundaries + max_cellsize: mesh size far from boundaries + + Returns: + Field ID for the minimum field + """ + # Distance field from boundary curves + gmsh.model.mesh.field.add("Distance", 1) + gmsh.model.mesh.field.setNumbers(1, "CurvesList", boundary_line_tags) + gmsh.model.mesh.field.setNumber(1, "Sampling", 200) + + # Threshold field for gradual size transition + gmsh.model.mesh.field.add("Threshold", 2) + gmsh.model.mesh.field.setNumber(2, "InField", 1) + gmsh.model.mesh.field.setNumber(2, "SizeMin", refined_cellsize) + gmsh.model.mesh.field.setNumber(2, "SizeMax", max_cellsize) + gmsh.model.mesh.field.setNumber(2, "DistMin", 0) + gmsh.model.mesh.field.setNumber(2, "DistMax", max_cellsize) + + return 2 + + +def setup_box_refinement( + field_id: int, + xmin: float, + ymin: float, + zmin: float, + xmax: float, + ymax: float, + zmax: float, + size_in: float, + size_out: float, +) -> None: + """Set up box-based mesh refinement. + + Args: + field_id: field ID to use + xmin, ymin, zmin: box minimum coordinates + xmax, ymax, zmax: box maximum coordinates + size_in: mesh size inside box + size_out: mesh size outside box + """ + gmsh.model.mesh.field.add("Box", field_id) + gmsh.model.mesh.field.setNumber(field_id, "VIn", size_in) + gmsh.model.mesh.field.setNumber(field_id, "VOut", size_out) + gmsh.model.mesh.field.setNumber(field_id, "XMin", xmin) + gmsh.model.mesh.field.setNumber(field_id, "XMax", xmax) + gmsh.model.mesh.field.setNumber(field_id, "YMin", ymin) + gmsh.model.mesh.field.setNumber(field_id, "YMax", ymax) + gmsh.model.mesh.field.setNumber(field_id, "ZMin", zmin) + gmsh.model.mesh.field.setNumber(field_id, "ZMax", zmax) + + +def finalize_mesh_fields(field_ids: list[int]) -> None: + """Finalize mesh fields by setting up minimum field. + + Args: + field_ids: list of field IDs to combine + """ + min_field_id = max(field_ids) + 1 + gmsh.model.mesh.field.add("Min", min_field_id) + gmsh.model.mesh.field.setNumbers(min_field_id, "FieldsList", field_ids) + gmsh.model.mesh.field.setAsBackgroundMesh(min_field_id) + + # Disable other mesh size sources + gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) + gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) + gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) + gmsh.option.setNumber("Mesh.Algorithm", 5) # Delaunay algorithm diff --git a/src/gsim/common/mesh/types.py b/src/gsim/common/mesh/types.py new file mode 100644 index 0000000..657a78c --- /dev/null +++ b/src/gsim/common/mesh/types.py @@ -0,0 +1,26 @@ +"""Generic mesh data types for GMSH mesh generation.""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from pathlib import Path + + +@dataclass +class MeshResult: + """Result from generic mesh generation. + + Attributes: + mesh_path: Path to the generated .msh file. + mesh_stats: Mesh statistics (nodes, elements, quality, etc.). + groups: Physical group info mapping: + {"volumes": {material: {"phys_group": int, "tags": [int]}}, + "outer_boundary": {"phys_group": int, "tags": [int]}} + """ + + mesh_path: Path + mesh_stats: dict = field(default_factory=dict) + groups: dict = field(default_factory=dict) + + +__all__ = ["MeshResult"] diff --git a/src/gsim/meep/__init__.py b/src/gsim/meep/__init__.py index 4bcd605..8f490f5 100644 --- a/src/gsim/meep/__init__.py +++ b/src/gsim/meep/__init__.py @@ -19,6 +19,7 @@ """ from gsim.gcloud import RunResult, register_result_parser +from gsim.meep.mesh_sim import MeepMeshSim from gsim.meep.models import ( FDTD, Domain, @@ -54,6 +55,7 @@ def _parse_meep_result(run_result: RunResult) -> SParameterResult: "DomainConfig", "Geometry", "Material", + "MeepMeshSim", "ModeSource", "ResolutionConfig", "SParameterResult", diff --git a/src/gsim/meep/mesh_config.py b/src/gsim/meep/mesh_config.py new file mode 100644 index 0000000..5c3241d --- /dev/null +++ b/src/gsim/meep/mesh_config.py @@ -0,0 +1,162 @@ +"""Mesh config writer for GMSH-based Meep simulations. + +Serializes the full photonic problem definition (mesh reference, materials, +source, monitors, domain, solver settings) to a ``mesh_config.json`` file +alongside the ``.msh`` mesh file. + +All spatial values are written in nanometers (length_scale = 1e-9). +""" + +from __future__ import annotations + +import json +import logging +from pathlib import Path +from typing import TYPE_CHECKING, Any + +if TYPE_CHECKING: + from gsim.common.mesh.types import MeshResult + from gsim.meep.models.api import FDTD, Domain, Material, ModeSource + +logger = logging.getLogger(__name__) + +# um → nm conversion factor (hardcoded for now) +_UM_TO_NM = 1000.0 + + +def _serialize_materials(materials: dict[str, float | Material]) -> dict[str, Any]: + """Convert materials dict to JSON-serializable form.""" + from gsim.meep.models.api import Material + + out: dict[str, Any] = {} + for name, val in materials.items(): + if isinstance(val, (int, float)): + out[name] = {"refractive_index": float(val), "extinction_coeff": 0.0} + elif isinstance(val, Material): + out[name] = {"refractive_index": val.n, "extinction_coeff": val.k} + else: + out[name] = val + return out + + +def _serialize_mesh_groups(groups: dict) -> dict[str, Any]: + """Extract serializable mesh group info (drop internal gmsh tags).""" + out: dict[str, Any] = {} + + if "volumes" in groups: + out["volumes"] = { + name: { + "phys_group": info["phys_group"], + "material": name, + } + for name, info in groups["volumes"].items() + } + + if "layer_volumes" in groups: + out["layer_volumes"] = { + name: { + "phys_group": info["phys_group"], + "material": info.get("material", name), + } + for name, info in groups["layer_volumes"].items() + } + + if groups.get("outer_boundary"): + out["outer_boundary"] = {"phys_group": groups["outer_boundary"]["phys_group"]} + + if "port_surfaces" in groups: + out["port_surfaces"] = {} + for name, info in groups["port_surfaces"].items(): + out["port_surfaces"][name] = { + "phys_group": info["phys_group"], + "center": [c * _UM_TO_NM for c in info["center"]], + "width": info["width"] * _UM_TO_NM, + "orientation": info["orientation"], + "layer": info["layer"], + "z_range": [z * _UM_TO_NM for z in info["z_range"]], + } + + return out + + +def _serialize_source(source: ModeSource) -> dict[str, Any]: + """Serialize source settings.""" + return { + "port": source.port, + "wavelength": source.wavelength * _UM_TO_NM, + "wavelength_span": source.wavelength_span * _UM_TO_NM, + "num_freqs": source.num_freqs, + } + + +def _serialize_domain(domain: Domain) -> dict[str, Any]: + """Serialize domain settings.""" + return { + "pml": domain.pml * _UM_TO_NM, + "margin_xy": domain.margin * _UM_TO_NM, + "margin_z_above": domain.margin_z_above * _UM_TO_NM, + "margin_z_below": domain.margin_z_below * _UM_TO_NM, + } + + +def _serialize_solver(solver: FDTD) -> dict[str, Any]: + """Serialize solver settings.""" + return { + "nanometers_per_cell": 1000.0 / solver.resolution, + "stopping": { + "mode": solver.stopping, + "threshold": solver.stopping_threshold, + }, + "wall_time_max": solver.wall_time_max, + } + + +def _serialize_mesh_stats(mesh_stats: dict) -> dict[str, Any]: + """Extract key mesh statistics.""" + out: dict[str, Any] = {} + if "nodes" in mesh_stats: + out["nodes"] = mesh_stats["nodes"] + if "tetrahedra" in mesh_stats: + out["tetrahedra"] = mesh_stats["tetrahedra"] + return out + + +def write_mesh_config( + mesh_result: MeshResult, + materials: dict[str, float | Material], + source: ModeSource, + monitors: list[str], + domain: Domain, + solver: FDTD, + output_dir: Path, +) -> Path: + """Write mesh_config.json alongside the .msh file. + + Args: + mesh_result: Result from ``generate_mesh()``. + materials: Material name → refractive index or Material object. + source: Mode source configuration. + monitors: List of monitor port names. + domain: Computational domain settings. + solver: FDTD solver settings. + output_dir: Directory to write the config file. + + Returns: + Path to the written ``mesh_config.json``. + """ + config = { + "length_scale": 1e-9, # mesh coordinates in nm; hardcoded for now + "mesh_filename": mesh_result.mesh_path.name, + "materials": _serialize_materials(materials), + "mesh_groups": _serialize_mesh_groups(mesh_result.groups), + "source": _serialize_source(source), + "monitors": list(monitors), + "domain": _serialize_domain(domain), + "solver": _serialize_solver(solver), + "mesh_stats": _serialize_mesh_stats(mesh_result.mesh_stats), + } + + config_path = Path(output_dir) / "mesh_config.json" + config_path.write_text(json.dumps(config, indent=2)) + logger.info("Mesh config written: %s", config_path) + return config_path diff --git a/src/gsim/meep/mesh_sim.py b/src/gsim/meep/mesh_sim.py new file mode 100644 index 0000000..7dba72f --- /dev/null +++ b/src/gsim/meep/mesh_sim.py @@ -0,0 +1,311 @@ +"""GMSH-based mesh generation for photonic simulations. + +Mirrors Palace's DrivenSim pattern: configure → mesh() → write_config(). +Produces a GMSH ``.msh`` + ``mesh_config.json`` that fully describes a +photonic problem. This is an alternative to Meep's existing prism-based +workflow. + +Example:: + + from gsim.meep import MeepMeshSim + + sim = MeepMeshSim() + sim.geometry(component=c, z_crop="auto") + sim.materials = {"si": 3.47, "SiO2": 1.44} + sim.source(port="o1", wavelength=1.55, wavelength_span=0.1, num_freqs=11) + sim.monitors = ["o1", "o2", "o3"] + sim.domain(pml=1.0, margin=0.5) + sim.solver(resolution=32) + + result = sim.mesh("./mesh-ybranch") + sim.write_config() + sim.plot_mesh(interactive=False) +""" + +from __future__ import annotations + +import logging +from pathlib import Path +from typing import Any + +from pydantic import BaseModel, ConfigDict, Field, PrivateAttr, field_validator + +from gsim.meep.models.api import ( + FDTD, + Domain, + Geometry, + Material, + ModeSource, +) + +logger = logging.getLogger(__name__) + + +class MeepMeshSim(BaseModel): + """GMSH-based mesh generation for photonic problems. + + Collects photonic simulation parameters and generates a GMSH mesh + + config JSON that together fully describe the problem. + + Uses the same API models as :class:`gsim.meep.Simulation` (``Geometry``, + ``ModeSource``, ``Domain``, ``FDTD``, ``Material``). + """ + + model_config = ConfigDict( + validate_assignment=True, + arbitrary_types_allowed=True, + ) + + geometry: Geometry = Field(default_factory=Geometry) + materials: dict[str, float | Material] = Field(default_factory=dict) + source: ModeSource = Field(default_factory=ModeSource) + monitors: list[str] = Field(default_factory=list) + domain: Domain = Field(default_factory=Domain) + solver: FDTD = Field(default_factory=FDTD) + + _stack_kwargs: dict[str, Any] = PrivateAttr(default_factory=dict) + _output_dir: Path | None = PrivateAttr(default=None) + _last_mesh_result: Any = PrivateAttr(default=None) + + # ------------------------------------------------------------------------- + # Validators + # ------------------------------------------------------------------------- + + @field_validator("materials", mode="before") + @classmethod + def _normalize_materials( + cls, + v: dict[str, float | Material | dict], + ) -> dict[str, float | Material]: + """Accept float shorthand: ``{"si": 3.47}`` → ``Material(n=3.47)``.""" + out: dict[str, float | Material] = {} + for name, val in v.items(): + if isinstance(val, (int, float)): + out[name] = Material(n=float(val)) + elif isinstance(val, dict): + out[name] = Material(**val) + else: + out[name] = val + return out + + # ------------------------------------------------------------------------- + # Internal: stack resolution + # ------------------------------------------------------------------------- + + def _ensure_stack(self) -> None: + """Lazily resolve the layer stack if not yet built.""" + if self.geometry.stack is not None: + return + + from gsim.common.stack import get_stack + + if self._stack_kwargs: + yaml_path = self._stack_kwargs.pop("yaml_path", None) + self.geometry.stack = get_stack(yaml_path=yaml_path, **self._stack_kwargs) + self._stack_kwargs["yaml_path"] = yaml_path + else: + self.geometry.stack = get_stack() + + # ------------------------------------------------------------------------- + # Internal: z-crop (reuses same logic as Simulation) + # ------------------------------------------------------------------------- + + def _apply_z_crop(self) -> None: + """Apply z-crop to the stack if geometry.z_crop is set.""" + if self.geometry.z_crop is None: + return + + from gsim.common.stack.extractor import Layer, LayerStack + from gsim.meep.ports import _find_highest_n_layer + + stack = self.geometry.stack + if stack is None: + raise ValueError("No stack configured for z-crop.") + + ref: Layer | None = None + if self.geometry.z_crop == "auto": + ref, best_n = _find_highest_n_layer(stack) + if ref is None or best_n <= 1.5: + raise ValueError( + "Could not auto-detect core layer (no layer with n > 1.5). " + "Set geometry.z_crop to an explicit layer name." + ) + else: + layer_name = self.geometry.z_crop + if layer_name not in stack.layers: + raise ValueError( + f"Layer '{layer_name}' not found. " + f"Available: {list(stack.layers.keys())}" + ) + ref = stack.layers[layer_name] + + z_lo = ref.zmin - self.domain.margin_z_below + z_hi = ref.zmax + self.domain.margin_z_above + + cropped: dict[str, Layer] = {} + for name, layer in stack.layers.items(): + if layer.zmax <= z_lo or layer.zmin >= z_hi: + continue + new_zmin = max(layer.zmin, z_lo) + new_zmax = min(layer.zmax, z_hi) + cropped[name] = layer.model_copy( + update={ + "zmin": new_zmin, + "zmax": new_zmax, + "thickness": new_zmax - new_zmin, + } + ) + + cropped_dielectrics = [] + for diel in stack.dielectrics: + if diel["zmax"] <= z_lo or diel["zmin"] >= z_hi: + continue + cropped_dielectrics.append( + { + **diel, + "zmin": max(diel["zmin"], z_lo), + "zmax": min(diel["zmax"], z_hi), + } + ) + + self.geometry.stack = LayerStack( + pdk_name=stack.pdk_name, + units=stack.units, + layers=cropped, + materials=stack.materials, + dielectrics=cropped_dielectrics, + simulation=stack.simulation, + ) + self.geometry.z_crop = None + + # ------------------------------------------------------------------------- + # mesh() + # ------------------------------------------------------------------------- + + def mesh( + self, + output_dir: str | Path, + *, + refined_mesh_size: float = 0.05, + max_mesh_size: float = 1.0, + margin: float | None = None, + air_margin: float = 2.0, + include_airbox: bool = False, + ) -> Any: + """Generate a GMSH mesh for this simulation. + + Creates ``mesh.msh`` inside *output_dir*. + + Args: + output_dir: Directory for mesh output files. + refined_mesh_size: Fine mesh size near waveguide boundaries (um). + max_mesh_size: Coarse mesh size in cladding/substrate (um). + margin: XY margin around geometry (um). Defaults to + ``domain.margin + domain.pml``. + air_margin: Airbox margin around dielectric envelope (um). + include_airbox: Whether to add a surrounding airbox volume. + Not needed for photonic simulations (default ``False``). + + Returns: + :class:`~gsim.common.mesh.types.MeshResult` with mesh path, + statistics, and physical group info. + + Raises: + ValueError: If no component or stack is configured. + """ + from gsim.common.mesh import generate_mesh + + if self.geometry.component is None: + raise ValueError( + "No component set. Assign sim.geometry(component=...) first." + ) + + self._ensure_stack() + if self.geometry.stack is None: + raise ValueError("Stack resolution failed.") + + self._apply_z_crop() + + if margin is None: + margin = self.domain.margin + self.domain.pml + + output_dir = Path(output_dir) + self._output_dir = output_dir + + result = generate_mesh( + component=self.geometry.component, + stack=self.geometry.stack, + output_dir=output_dir, + model_name="mesh", + refined_mesh_size=refined_mesh_size, + max_mesh_size=max_mesh_size, + margin=margin, + air_margin=air_margin, + include_airbox=include_airbox, + mesh_scale=1000.0, # um → nm + ) + + self._last_mesh_result = result + logger.info("Mesh generated: %s", result.mesh_path) + return result + + # ------------------------------------------------------------------------- + # write_config() + # ------------------------------------------------------------------------- + + def write_config(self) -> Path: + """Write ``mesh_config.json`` alongside the mesh. + + Must be called after :meth:`mesh`. + + Returns: + Path to the written config file. + + Raises: + ValueError: If :meth:`mesh` has not been called. + """ + from gsim.meep.mesh_config import write_mesh_config + + if self._last_mesh_result is None or self._output_dir is None: + raise ValueError("Call mesh() first to generate a mesh.") + + return write_mesh_config( + mesh_result=self._last_mesh_result, + materials=self.materials, + source=self.source, + monitors=self.monitors, + domain=self.domain, + solver=self.solver, + output_dir=self._output_dir, + ) + + # ------------------------------------------------------------------------- + # plot_mesh() + # ------------------------------------------------------------------------- + + def plot_mesh( + self, + show_groups: list[str] | None = None, + interactive: bool = True, + ) -> None: + """Visualize the generated mesh. + + Delegates to :func:`gsim.viz.plot_mesh`. + + Args: + show_groups: Physical group names to display (None = all). + interactive: Show interactive 3D viewer. + + Raises: + ValueError: If :meth:`mesh` has not been called. + """ + from gsim.viz import plot_mesh + + if self._last_mesh_result is None: + raise ValueError("Call mesh() first to generate a mesh.") + + plot_mesh( + self._last_mesh_result.mesh_path, + show_groups=show_groups, + interactive=interactive, + ) diff --git a/src/gsim/palace/mesh/config_generator.py b/src/gsim/palace/mesh/config_generator.py index 581faa2..9e88ab1 100644 --- a/src/gsim/palace/mesh/config_generator.py +++ b/src/gsim/palace/mesh/config_generator.py @@ -9,7 +9,7 @@ from pathlib import Path from typing import TYPE_CHECKING -import gmsh +from gsim.common.mesh.generator import collect_mesh_stats if TYPE_CHECKING: from gsim.common.stack import LayerStack @@ -231,117 +231,6 @@ def generate_palace_config( return config_path -def collect_mesh_stats() -> dict: - """Collect mesh statistics from gmsh after mesh generation. - - Must be called while gmsh is initialized and the mesh is generated. - - Returns: - Dict with mesh statistics including: - - bbox: Bounding box coordinates - - nodes: Number of nodes - - elements: Total element count - - tetrahedra: Tet count - - quality: Shape quality metrics (gamma) - - sicn: Signed Inverse Condition Number - - edge_length: Min/max edge lengths - - groups: Physical group info - """ - stats = {} - - # Get bounding box - try: - xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.getBoundingBox(-1, -1) - stats["bbox"] = { - "xmin": xmin, - "ymin": ymin, - "zmin": zmin, - "xmax": xmax, - "ymax": ymax, - "zmax": zmax, - } - except Exception: - pass - - # Get node count - try: - node_tags, _, _ = gmsh.model.mesh.getNodes() - stats["nodes"] = len(node_tags) - except Exception: - pass - - # Get element counts and collect tet tags for quality - tet_tags = [] - try: - element_types, element_tags, _ = gmsh.model.mesh.getElements() - total_elements = sum(len(tags) for tags in element_tags) - stats["elements"] = total_elements - - # Count tetrahedra (type 4) and save tags - for etype, tags in zip(element_types, element_tags, strict=False): - if etype == 4: # 4-node tetrahedron - stats["tetrahedra"] = len(tags) - tet_tags = list(tags) - except Exception: - pass - - # Get mesh quality for tetrahedra - if tet_tags: - # Gamma: inscribed/circumscribed radius ratio (shape quality) - try: - qualities = gmsh.model.mesh.getElementQualities(tet_tags, "gamma") - if len(qualities) > 0: - stats["quality"] = { - "min": round(min(qualities), 3), - "max": round(max(qualities), 3), - "mean": round(sum(qualities) / len(qualities), 3), - } - except Exception: - pass - - # SICN: Signed Inverse Condition Number (negative = invalid element) - try: - sicn = gmsh.model.mesh.getElementQualities(tet_tags, "minSICN") - if len(sicn) > 0: - sicn_min = min(sicn) - invalid_count = sum(1 for s in sicn if s < 0) - stats["sicn"] = { - "min": round(sicn_min, 3), - "mean": round(sum(sicn) / len(sicn), 3), - "invalid": invalid_count, - } - except Exception: - pass - - # Edge lengths - try: - min_edges = gmsh.model.mesh.getElementQualities(tet_tags, "minEdge") - max_edges = gmsh.model.mesh.getElementQualities(tet_tags, "maxEdge") - if len(min_edges) > 0 and len(max_edges) > 0: - stats["edge_length"] = { - "min": round(min(min_edges), 3), - "max": round(max(max_edges), 3), - } - except Exception: - pass - - # Get physical groups with tags - try: - groups = {"volumes": [], "surfaces": []} - for dim, tag in gmsh.model.getPhysicalGroups(): - name = gmsh.model.getPhysicalName(dim, tag) - entry = {"name": name, "tag": tag} - if dim == 3: - groups["volumes"].append(entry) - elif dim == 2: - groups["surfaces"].append(entry) - stats["groups"] = groups - except Exception: - pass - - return stats - - def write_config( mesh_result, stack: LayerStack, diff --git a/src/gsim/palace/mesh/generator.py b/src/gsim/palace/mesh/generator.py index cf71299..d503e16 100644 --- a/src/gsim/palace/mesh/generator.py +++ b/src/gsim/palace/mesh/generator.py @@ -14,9 +14,10 @@ import gmsh +from gsim.common.mesh.generator import collect_mesh_stats + from . import gmsh_utils from .config_generator import ( - collect_mesh_stats, generate_palace_config, write_config, ) diff --git a/src/gsim/palace/mesh/geometry.py b/src/gsim/palace/mesh/geometry.py index 59201bd..424de0a 100644 --- a/src/gsim/palace/mesh/geometry.py +++ b/src/gsim/palace/mesh/geometry.py @@ -1,140 +1,32 @@ """Geometry extraction and creation for Palace mesh generation. -This module handles extracting polygons from gdsfactory components -and creating 3D geometry in gmsh. +Generic geometry helpers (GeometryData, extract_geometry, get_layer_info, +add_dielectrics) are provided by ``gsim.common.mesh.geometry`` and +re-exported here for backward compatibility. + +This module keeps Palace-specific functions: ``add_metals`` and ``add_ports``. """ from __future__ import annotations import math -from dataclasses import dataclass from typing import TYPE_CHECKING -from . import gmsh_utils +from gsim.common.mesh import gmsh_utils + +# Re-export generic helpers so existing Palace code keeps working +from gsim.common.mesh.geometry import ( + GeometryData, + add_dielectrics, + extract_geometry, + get_layer_info, +) if TYPE_CHECKING: from gsim.common.stack import LayerStack from gsim.palace.ports.config import PalacePort -@dataclass -class GeometryData: - """Container for geometry data extracted from component.""" - - polygons: list # List of (layer_num, pts_x, pts_y, holes) tuples - bbox: tuple[float, float, float, float] # (xmin, ymin, xmax, ymax) - layer_bboxes: dict # layer_num -> (xmin, ymin, xmax, ymax) - - -def extract_geometry(component, stack: LayerStack) -> GeometryData: - """Extract polygon geometry from a gdsfactory component. - - Args: - component: gdsfactory Component - stack: LayerStack for layer mapping - - Returns: - GeometryData with polygons and bounding boxes - """ - polygons = [] - global_bbox = [math.inf, math.inf, -math.inf, -math.inf] - layer_bboxes = {} - - # Get polygons from component - polygons_by_index = component.get_polygons() - - # Build layer_index -> GDS tuple mapping - layout = component.kcl.layout - index_to_gds = {} - for layer_index in range(layout.layers()): - if layout.is_valid_layer(layer_index): - info = layout.get_info(layer_index) - index_to_gds[layer_index] = (info.layer, info.datatype) - - # Build GDS tuple -> layer number mapping - gds_to_layernum = {} - for layer_data in stack.layers.values(): - gds_tuple = tuple(layer_data.gds_layer) - gds_to_layernum[gds_tuple] = gds_tuple[0] - - # Convert polygons - for layer_index, polys in polygons_by_index.items(): - gds_tuple = index_to_gds.get(layer_index) - if gds_tuple is None: - continue - - layernum = gds_to_layernum.get(gds_tuple) - if layernum is None: - continue - - for poly in polys: - # Convert klayout polygon to lists (nm -> um) - points = list(poly.each_point_hull()) - if len(points) < 3: - continue - - pts_x = [pt.x / 1000.0 for pt in points] - pts_y = [pt.y / 1000.0 for pt in points] - - # Extract holes from polygon - holes = [] - for hole_idx in range(poly.holes()): - hole_pts = list(poly.each_point_hole(hole_idx)) - if len(hole_pts) >= 3: - hx = [pt.x / 1000.0 for pt in hole_pts] - hy = [pt.y / 1000.0 for pt in hole_pts] - holes.append((hx, hy)) - - polygons.append((layernum, pts_x, pts_y, holes)) - - # Update bounding boxes - xmin, xmax = min(pts_x), max(pts_x) - ymin, ymax = min(pts_y), max(pts_y) - - global_bbox[0] = min(global_bbox[0], xmin) - global_bbox[1] = min(global_bbox[1], ymin) - global_bbox[2] = max(global_bbox[2], xmax) - global_bbox[3] = max(global_bbox[3], ymax) - - if layernum not in layer_bboxes: - layer_bboxes[layernum] = [xmin, ymin, xmax, ymax] - else: - bbox = layer_bboxes[layernum] - bbox[0] = min(bbox[0], xmin) - bbox[1] = min(bbox[1], ymin) - bbox[2] = max(bbox[2], xmax) - bbox[3] = max(bbox[3], ymax) - - return GeometryData( - polygons=polygons, - bbox=(global_bbox[0], global_bbox[1], global_bbox[2], global_bbox[3]), - layer_bboxes=layer_bboxes, - ) - - -def get_layer_info(stack: LayerStack, gds_layer: int) -> dict | None: - """Get layer info from stack by GDS layer number. - - Args: - stack: LayerStack with layer definitions - gds_layer: GDS layer number - - Returns: - Dict with layer info or None if not found - """ - for name, layer in stack.layers.items(): - if layer.gds_layer[0] == gds_layer: - return { - "name": name, - "zmin": layer.zmin, - "zmax": layer.zmax, - "thickness": layer.zmax - layer.zmin, - "material": layer.material, - "type": layer.layer_type, - } - return None - - def add_metals( kernel, geometry: GeometryData, @@ -445,6 +337,7 @@ def build_entities( return entities + def add_ports( kernel, ports: list[PalacePort], diff --git a/src/gsim/palace/mesh/gmsh_utils.py b/src/gsim/palace/mesh/gmsh_utils.py index f461368..84c2b32 100644 --- a/src/gsim/palace/mesh/gmsh_utils.py +++ b/src/gsim/palace/mesh/gmsh_utils.py @@ -1,18 +1,24 @@ -"""Gmsh utility functions for Palace mesh generation. +"""Gmsh utility functions — backward-compatible re-export from common. -Extracted from gds2palace/util_simulation_setup.py and adapted -to work directly with palace-api data structures. +All generic GMSH helpers now live in ``gsim.common.mesh.gmsh_utils``. +This shim keeps existing ``from gsim.palace.mesh import gmsh_utils`` working. + +Palace-specific helpers (``Entity``, ``run_boolean_pipeline``) are defined +here because they encode Palace mesh-order conventions and physical-group +labelling that don't belong in the solver-agnostic common module. """ from __future__ import annotations -import math - import gmsh -import numpy as np + +# Re-export all generic helpers from common +from gsim.common.mesh.gmsh_utils import * # noqa: F403 +from gsim.common.mesh.gmsh_utils import is_vertical_surface + # --------------------------------------------------------------------------- -# Minimalistic meshwell-style entity + boolean pipeline +# Palace-specific: Entity + boolean pipeline # --------------------------------------------------------------------------- @@ -172,500 +178,3 @@ def run_boolean_pipeline(entities: list[Entity]) -> dict[str, int]: pg_map[name] = pg_tag return pg_map - - -def create_box( - kernel, - xmin: float, - ymin: float, - zmin: float, - xmax: float, - ymax: float, - zmax: float, - meshseed: float = 0, -) -> int: - """Create a 3D box volume in gmsh. - - Args: - kernel: gmsh.model.occ kernel - xmin, ymin, zmin: minimum coordinates - xmax, ymax, zmax: maximum coordinates - meshseed: mesh seed size at corners (0 = auto) - - Returns: - Volume tag of created box - """ - if meshseed == 0: - # Use simple addBox - return kernel.addBox(xmin, ymin, zmin, xmax - xmin, ymax - ymin, zmax - zmin) - - # Create box with explicit mesh seed at corners - pt1 = kernel.addPoint(xmin, ymin, zmin, meshseed, -1) - pt2 = kernel.addPoint(xmin, ymax, zmin, meshseed, -1) - pt3 = kernel.addPoint(xmax, ymax, zmin, meshseed, -1) - pt4 = kernel.addPoint(xmax, ymin, zmin, meshseed, -1) - - line1 = kernel.addLine(pt1, pt2, -1) - line2 = kernel.addLine(pt2, pt3, -1) - line3 = kernel.addLine(pt3, pt4, -1) - line4 = kernel.addLine(pt4, pt1, -1) - linetaglist = [line1, line2, line3, line4] - - curvetag = kernel.addCurveLoop(linetaglist, tag=-1) - surfacetag = kernel.addPlaneSurface([curvetag], tag=-1) - returnval = kernel.extrude([(2, surfacetag)], 0, 0, zmax - zmin) - volumetag = returnval[1][1] - - return volumetag - - -def create_polygon_surface( - kernel, - pts_x: list[float], - pts_y: list[float], - z: float, - meshseed: float = 0, - holes: list[tuple[list[float], list[float]]] | None = None, -) -> int | None: - """Create a planar surface from polygon vertices at z height. - - Args: - kernel: gmsh.model.occ kernel - pts_x: list of x coordinates - pts_y: list of y coordinates - z: z coordinate of the surface - meshseed: mesh seed size at vertices (0 = auto) - holes: list of (hole_pts_x, hole_pts_y) tuples for interior holes - - Returns: - Surface tag, or None if polygon is invalid - """ - numvertices = len(pts_x) - if numvertices < 3: - return None - - linetaglist = [] - vertextaglist = [] - - # Create vertices - for v in range(numvertices): - vertextag = kernel.addPoint(pts_x[v], pts_y[v], z, meshseed, -1) - vertextaglist.append(vertextag) - - # Create lines connecting vertices - for v in range(numvertices): - pt_start = vertextaglist[v] - pt_end = vertextaglist[(v + 1) % numvertices] - try: - linetag = kernel.addLine(pt_start, pt_end, -1) - linetaglist.append(linetag) - except Exception: - pass # Skip degenerate lines - - if len(linetaglist) < 3: - return None - - # Create outer curve loop - curvetag = kernel.addCurveLoop(linetaglist, tag=-1) - - # Create hole curve loops - all_loops = [curvetag] - for hx, hy in holes or []: - hole_lines = [] - hole_verts = [] - for v in range(len(hx)): - vtag = kernel.addPoint(hx[v], hy[v], z, meshseed, -1) - hole_verts.append(vtag) - for v in range(len(hx)): - try: - ltag = kernel.addLine(hole_verts[v], hole_verts[(v + 1) % len(hx)], -1) - hole_lines.append(ltag) - except Exception: - pass - if len(hole_lines) >= 3: - hole_loop = kernel.addCurveLoop(hole_lines, tag=-1) - all_loops.append(hole_loop) - - surfacetag = kernel.addPlaneSurface(all_loops, tag=-1) - - return surfacetag - - -def extrude_polygon( - kernel, - pts_x: list[float], - pts_y: list[float], - zmin: float, - thickness: float, - meshseed: float = 0, - holes: list[tuple[list[float], list[float]]] | None = None, -) -> int | None: - """Create an extruded polygon volume (for vias, metals). - - Args: - kernel: gmsh.model.occ kernel - pts_x: list of x coordinates - pts_y: list of y coordinates - zmin: base z coordinate - thickness: extrusion height - meshseed: mesh seed size at vertices - holes: list of (hole_pts_x, hole_pts_y) tuples for interior holes - - Returns: - Volume tag if thickness > 0, surface tag if thickness == 0, or None if invalid - """ - surfacetag = create_polygon_surface( - kernel, pts_x, pts_y, zmin, meshseed, holes=holes - ) - if surfacetag is None: - return None - - if thickness > 0: - result = kernel.extrude([(2, surfacetag)], 0, 0, thickness) - # result[1] contains the volume (dim=3, tag) - return result[1][1] - - return surfacetag - - -def create_port_rectangle( - kernel, - xmin: float, - ymin: float, - zmin: float, - xmax: float, - ymax: float, - zmax: float, - meshseed: float = 0, -) -> int: - """Create a rectangular surface for a port. - - Handles both horizontal (z-plane) and vertical port surfaces. - - Args: - kernel: gmsh.model.occ kernel - xmin, ymin, zmin: minimum coordinates - xmax, ymax, zmax: maximum coordinates - meshseed: mesh seed size at corners - - Returns: - Surface tag of created port rectangle - """ - # Determine port orientation - dx = xmax - xmin - dz = zmax - zmin - - if dz < 1e-6: - # Horizontal port (in xy plane) - pt1 = kernel.addPoint(xmin, ymin, zmin, meshseed, -1) - pt2 = kernel.addPoint(xmin, ymax, zmin, meshseed, -1) - pt3 = kernel.addPoint(xmax, ymax, zmin, meshseed, -1) - pt4 = kernel.addPoint(xmax, ymin, zmin, meshseed, -1) - elif dx < 1e-6: - # Vertical port in yz plane - pt1 = kernel.addPoint(xmin, ymin, zmin, meshseed, -1) - pt2 = kernel.addPoint(xmin, ymax, zmin, meshseed, -1) - pt3 = kernel.addPoint(xmin, ymax, zmax, meshseed, -1) - pt4 = kernel.addPoint(xmin, ymin, zmax, meshseed, -1) - else: - # Vertical port in xz plane - pt1 = kernel.addPoint(xmin, ymin, zmin, meshseed, -1) - pt2 = kernel.addPoint(xmin, ymin, zmax, meshseed, -1) - pt3 = kernel.addPoint(xmax, ymin, zmax, meshseed, -1) - pt4 = kernel.addPoint(xmax, ymin, zmin, meshseed, -1) - - line1 = kernel.addLine(pt1, pt2, -1) - line2 = kernel.addLine(pt2, pt3, -1) - line3 = kernel.addLine(pt3, pt4, -1) - line4 = kernel.addLine(pt4, pt1, -1) - linetaglist = [line1, line2, line3, line4] - - curvetag = kernel.addCurveLoop(linetaglist, tag=-1) - surfacetag = kernel.addPlaneSurface([curvetag], tag=-1) - - return surfacetag - - -def fragment_all(kernel) -> tuple[list, list]: - """Fragment all geometry to ensure conformal mesh at intersections. - - Args: - kernel: gmsh.model.occ kernel - - Returns: - (geom_dimtags, geom_map) - original dimtags and mapping to new tags - """ - geom_dimtags = [x for x in kernel.getEntities() if x[0] in (2, 3)] - _, geom_map = kernel.fragment(geom_dimtags, []) - kernel.synchronize() - return geom_dimtags, geom_map - - -def get_tags_after_fragment( - original_tags: list[int], - geom_dimtags: list, - geom_map: list, - dimension: int = 2, -) -> list[int]: - """Get new tags after fragmenting, given original tags. - - Tags change after gmsh fragment operation. This function maps - original tags to their new values using the fragment mapping. - - Args: - original_tags: list of tags before fragmenting - geom_dimtags: list of all original dimtags before fragmenting - geom_map: mapping from fragment() function - dimension: dimension for tags (2=surfaces, 3=volumes) - - Returns: - List of new tags after fragmenting - """ - if isinstance(original_tags, int): - original_tags = [original_tags] - - indices = [ - i - for i, x in enumerate(geom_dimtags) - if x[0] == dimension and (x[1] in original_tags) - ] - raw = [geom_map[i] for i in indices] - flat = [item for sublist in raw for item in sublist] - newtags = [s[-1] for s in flat] - - return newtags - - -def assign_physical_group( - dim: int, - tags: list[int], - name: str, -) -> int: - """Assign tags to a physical group with a name. - - Args: - dim: dimension (2=surfaces, 3=volumes) - tags: list of entity tags - name: physical group name - - Returns: - Physical group tag - """ - if not tags: - return -1 - phys_group = gmsh.model.addPhysicalGroup(dim, tags, tag=-1) - gmsh.model.setPhysicalName(dim, phys_group, name) - return phys_group - - -def get_surface_normal(surface_tag: int) -> np.ndarray: - """Get the normal vector of a surface. - - Args: - surface_tag: surface entity tag - - Returns: - Normal vector as numpy array [nx, ny, nz] - """ - # Get the boundary of the surface - boundary_lines = gmsh.model.getBoundary([(2, surface_tag)], oriented=True) - - # Get points from these lines - points = [] - seen_points = set() - - for _dim, line_tag in boundary_lines: - line_points = gmsh.model.getBoundary([(1, line_tag)], oriented=True) - for _pdim, ptag in line_points: - if ptag not in seen_points: - coord = gmsh.model.getValue(0, ptag, []) - points.append(np.array(coord)) - seen_points.add(ptag) - if len(points) == 3: - break - if len(points) == 3: - break - - if len(points) < 3: - return np.array([0, 0, 1]) # Default to z-normal - - # Compute surface normal using cross product - v1 = points[1] - points[0] - v2 = points[2] - points[0] - normal = np.cross(v1, v2) - norm = np.linalg.norm(normal) - if norm > 0: - normal = normal / norm - return normal - - -def is_vertical_surface(surface_tag: int) -> bool: - """Check if a surface is vertical (not in xy plane). - - Args: - surface_tag: surface entity tag - - Returns: - True if surface is vertical (z component of normal is ~0) - """ - normal = get_surface_normal(surface_tag) - n = normal[2] - if not np.isnan(n): - return int(abs(n)) == 0 - return False - - -def get_volumes_at_z_range( - zmin: float, - zmax: float, - delta: float = 0.001, -) -> list[tuple[int, int]]: - """Get all volumes within a z-coordinate range. - - Args: - zmin: minimum z coordinate - zmax: maximum z coordinate - delta: tolerance for z comparison - - Returns: - List of (dim, tag) tuples for volumes in the z range - """ - volumes_in_bbox = gmsh.model.getEntitiesInBoundingBox( - -math.inf, - -math.inf, - zmin - delta / 2, - math.inf, - math.inf, - zmax + delta / 2, - 3, - ) - - volume_list = [] - for volume in volumes_in_bbox: - volume_tag = volume[1] - _, _, vzmin, _, _, vzmax = gmsh.model.getBoundingBox(3, volume_tag) - if ( - abs(vzmin - (zmin - delta / 2)) < delta - and abs(vzmax - (zmax + delta / 2)) < delta - ): - volume_list.append(volume) - - return volume_list - - -def get_surfaces_at_z(z: float, delta: float = 0.001) -> list[tuple[int, int]]: - """Get all surfaces at a specific z coordinate. - - Args: - z: z coordinate - delta: tolerance for z comparison - - Returns: - List of (dim, tag) tuples for surfaces at z - """ - return gmsh.model.getEntitiesInBoundingBox( - -math.inf, - -math.inf, - z - delta / 2, - math.inf, - math.inf, - z + delta / 2, - 2, - ) - - -def get_boundary_lines(surface_tag: int, kernel) -> list[int]: - """Get all boundary line tags of a surface. - - Args: - surface_tag: surface entity tag - kernel: gmsh.model.occ kernel - - Returns: - List of curve/line tags forming the surface boundary - """ - _clt, ct = kernel.getCurveLoops(surface_tag) - lines = [] - for curvetag in ct: - lines.extend(curvetag) - return lines - - -def setup_mesh_refinement( - boundary_line_tags: list[int], - refined_cellsize: float, - max_cellsize: float, -) -> int: - """Set up mesh refinement near boundary lines. - - Args: - boundary_line_tags: list of curve tags for refinement - refined_cellsize: mesh size near boundaries - max_cellsize: mesh size far from boundaries - - Returns: - Field ID for the minimum field - """ - # Distance field from boundary curves - gmsh.model.mesh.field.add("Distance", 1) - gmsh.model.mesh.field.setNumbers(1, "CurvesList", boundary_line_tags) - gmsh.model.mesh.field.setNumber(1, "Sampling", 200) - - # Threshold field for gradual size transition - gmsh.model.mesh.field.add("Threshold", 2) - gmsh.model.mesh.field.setNumber(2, "InField", 1) - gmsh.model.mesh.field.setNumber(2, "SizeMin", refined_cellsize) - gmsh.model.mesh.field.setNumber(2, "SizeMax", max_cellsize) - gmsh.model.mesh.field.setNumber(2, "DistMin", 0) - gmsh.model.mesh.field.setNumber(2, "DistMax", max_cellsize) - - return 2 - - -def setup_box_refinement( - field_id: int, - xmin: float, - ymin: float, - zmin: float, - xmax: float, - ymax: float, - zmax: float, - size_in: float, - size_out: float, -) -> None: - """Set up box-based mesh refinement. - - Args: - field_id: field ID to use - xmin, ymin, zmin: box minimum coordinates - xmax, ymax, zmax: box maximum coordinates - size_in: mesh size inside box - size_out: mesh size outside box - """ - gmsh.model.mesh.field.add("Box", field_id) - gmsh.model.mesh.field.setNumber(field_id, "VIn", size_in) - gmsh.model.mesh.field.setNumber(field_id, "VOut", size_out) - gmsh.model.mesh.field.setNumber(field_id, "XMin", xmin) - gmsh.model.mesh.field.setNumber(field_id, "XMax", xmax) - gmsh.model.mesh.field.setNumber(field_id, "YMin", ymin) - gmsh.model.mesh.field.setNumber(field_id, "YMax", ymax) - gmsh.model.mesh.field.setNumber(field_id, "ZMin", zmin) - gmsh.model.mesh.field.setNumber(field_id, "ZMax", zmax) - - -def finalize_mesh_fields(field_ids: list[int]) -> None: - """Finalize mesh fields by setting up minimum field. - - Args: - field_ids: list of field IDs to combine - """ - min_field_id = max(field_ids) + 1 - gmsh.model.mesh.field.add("Min", min_field_id) - gmsh.model.mesh.field.setNumbers(min_field_id, "FieldsList", field_ids) - gmsh.model.mesh.field.setAsBackgroundMesh(min_field_id) - - # Disable other mesh size sources - gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) - gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) - gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) - gmsh.option.setNumber("Mesh.Algorithm", 5) # Delaunay algorithm diff --git a/src/gsim/palace/mesh/groups.py b/src/gsim/palace/mesh/groups.py index 360254d..4e997dc 100644 --- a/src/gsim/palace/mesh/groups.py +++ b/src/gsim/palace/mesh/groups.py @@ -9,7 +9,7 @@ import logging from typing import TYPE_CHECKING -from . import gmsh_utils +from gsim.common.mesh import gmsh_utils if TYPE_CHECKING: from gsim.common.stack import LayerStack diff --git a/src/gsim/viz.py b/src/gsim/viz.py index cd2b716..16ce2a7 100644 --- a/src/gsim/viz.py +++ b/src/gsim/viz.py @@ -96,27 +96,57 @@ def _plot_wireframe( mesh = pv.read(msh_path) plotter = _make_plotter(interactive) + # Determine which groups to display if show_groups: ids = [ tag for tag, name in group_map.items() if any(p in name for p in show_groups) ] - colors = ["red", "blue", "green", "orange", "purple", "cyan"] - for i, gid in enumerate(ids): - subset = mesh.extract_cells(mesh.cell_data["gmsh:physical"] == gid) - if subset.n_cells > 0: - plotter.add_mesh( - subset, - style="wireframe", - color=colors[i % len(colors)], - line_width=1, - label=group_map.get(gid, str(gid)), - ) - if ids: - plotter.add_legend() # type: ignore[call-arg] else: - plotter.add_mesh(mesh, style="wireframe", color="black", line_width=1) # type: ignore[arg-type] + # Show all groups + ids = list(group_map.keys()) + + # Find the largest volume group to render transparently + largest_vol_id = None + largest_vol_count = 0 + if not show_groups: + for gid in ids: + subset = mesh.extract_cells(mesh.cell_data["gmsh:physical"] == gid) + if subset.n_cells > largest_vol_count: + largest_vol_count = subset.n_cells + largest_vol_id = gid + + colors = ["red", "blue", "green", "orange", "purple", "cyan"] + color_idx = 0 + for gid in ids: + subset = mesh.extract_cells(mesh.cell_data["gmsh:physical"] == gid) + if subset.n_cells == 0: + continue + name = group_map.get(gid, str(gid)) + color = colors[color_idx % len(colors)] + color_idx += 1 + + if not show_groups and gid == largest_vol_id: + # Render the largest volume (e.g. cladding) as transparent + # wireframe so inner structures are visible + plotter.add_mesh( + subset, + style="wireframe", + color=color, + opacity=0.15, + line_width=1, + label=name, + ) + else: + plotter.add_mesh( + subset, + style="wireframe", + color=color, + line_width=1, + label=name, + ) + plotter.add_legend() _finish(plotter, msh_path, output=output, interactive=interactive) diff --git a/tests/common/__init__.py b/tests/common/__init__.py new file mode 100644 index 0000000..65309b3 --- /dev/null +++ b/tests/common/__init__.py @@ -0,0 +1 @@ +"""Tests for gsim.common modules.""" diff --git a/tests/common/test_mesh.py b/tests/common/test_mesh.py new file mode 100644 index 0000000..6bc0c80 --- /dev/null +++ b/tests/common/test_mesh.py @@ -0,0 +1,221 @@ +"""Tests for gsim.common.mesh module.""" + +from __future__ import annotations + +import pytest + +from gsim.common.mesh.geometry import GeometryData, get_layer_info +from gsim.common.mesh.types import MeshResult + + +class TestGeometryData: + """Test GeometryData dataclass.""" + + def test_create_empty(self): + gd = GeometryData(polygons=[], bbox=(0, 0, 10, 10), layer_bboxes={}) + assert gd.polygons == [] + assert gd.bbox == (0, 0, 10, 10) + assert gd.layer_bboxes == {} + + def test_create_with_polygons(self): + polys = [(1, [0.0, 1.0, 1.0], [0.0, 0.0, 1.0], [])] + gd = GeometryData( + polygons=polys, + bbox=(0, 0, 1, 1), + layer_bboxes={1: [0, 0, 1, 1]}, + ) + assert len(gd.polygons) == 1 + assert gd.polygons[0][0] == 1 + assert 1 in gd.layer_bboxes + + +class TestMeshResult: + """Test MeshResult dataclass.""" + + def test_create_minimal(self, tmp_path): + msh = tmp_path / "test.msh" + result = MeshResult(mesh_path=msh) + assert result.mesh_path == msh + assert result.mesh_stats == {} + assert result.groups == {} + + def test_create_with_data(self, tmp_path): + result = MeshResult( + mesh_path=tmp_path / "test.msh", + mesh_stats={"nodes": 100}, + groups={"volumes": {"SiO2": {"phys_group": 1, "tags": [1, 2]}}}, + ) + assert result.mesh_stats["nodes"] == 100 + assert "SiO2" in result.groups["volumes"] + + +class TestGetLayerInfo: + """Test get_layer_info with a mock LayerStack.""" + + @pytest.fixture + def mock_stack(self): + """Create a minimal LayerStack-like object for testing.""" + from gsim.common.stack.extractor import Layer, LayerStack + + stack = LayerStack(pdk_name="test") + stack.layers["metal1"] = Layer( + name="metal1", + gds_layer=(10, 0), + zmin=0.0, + zmax=0.5, + thickness=0.5, + material="copper", + layer_type="conductor", + ) + stack.layers["oxide"] = Layer( + name="oxide", + gds_layer=(20, 0), + zmin=-1.0, + zmax=0.0, + thickness=1.0, + material="SiO2", + layer_type="dielectric", + ) + return stack + + def test_found(self, mock_stack): + info = get_layer_info(mock_stack, 10) + assert info is not None + assert info["name"] == "metal1" + assert info["material"] == "copper" + assert info["type"] == "conductor" + assert info["zmin"] == 0.0 + assert info["thickness"] == 0.5 + + def test_not_found(self, mock_stack): + info = get_layer_info(mock_stack, 999) + assert info is None + + def test_dielectric_layer(self, mock_stack): + info = get_layer_info(mock_stack, 20) + assert info is not None + assert info["name"] == "oxide" + assert info["type"] == "dielectric" + + +class TestExtractGeometry: + """Test extract_geometry with a real gdsfactory component.""" + + @pytest.fixture + def simple_component_and_stack(self): + """Create a simple gdsfactory component with matching stack.""" + import gdsfactory as gf + + from gsim.common.stack.extractor import Layer, LayerStack + + gf.gpdk.PDK.activate() + c = gf.Component() + c.add_polygon( + [(0, 0), (10, 0), (10, 5), (0, 5)], + layer=(1, 0), + ) + + stack = LayerStack(pdk_name="test") + stack.layers["metal1"] = Layer( + name="metal1", + gds_layer=(1, 0), + zmin=0.0, + zmax=0.5, + thickness=0.5, + material="copper", + layer_type="conductor", + ) + return c, stack + + def test_extract_geometry(self, simple_component_and_stack): + from gsim.common.mesh.geometry import extract_geometry + + c, stack = simple_component_and_stack + gd = extract_geometry(c, stack) + + assert len(gd.polygons) == 1 + layernum, pts_x, _pts_y, holes = gd.polygons[0] + assert layernum == 1 + assert len(pts_x) >= 3 + assert holes == [] + # Bbox should match the polygon (in um) + assert gd.bbox[0] == pytest.approx(0.0, abs=0.01) + assert gd.bbox[2] == pytest.approx(10.0, abs=0.01) # 10um + + +class TestAddLayerVolumes: + """Test add_layer_volumes with gmsh.""" + + @pytest.fixture + def geometry_and_stack(self): + """Create geometry data and stack for testing.""" + from gsim.common.stack.extractor import Layer, LayerStack + + gd = GeometryData( + polygons=[ + (1, [0.0, 10.0, 10.0, 0.0], [0.0, 0.0, 5.0, 5.0], []), + ], + bbox=(0, 0, 10, 5), + layer_bboxes={1: [0, 0, 10, 5]}, + ) + stack = LayerStack(pdk_name="test") + stack.layers["metal1"] = Layer( + name="metal1", + gds_layer=(1, 0), + zmin=0.0, + zmax=0.5, + thickness=0.5, + material="copper", + layer_type="conductor", + ) + return gd, stack + + def test_add_layer_volumes(self, geometry_and_stack): + import gmsh + + from gsim.common.mesh.geometry import add_layer_volumes + + gd, stack = geometry_and_stack + + gmsh.initialize() + gmsh.model.add("test_volumes") + kernel = gmsh.model.occ + + try: + result = add_layer_volumes(kernel, gd, stack) + assert "metal1" in result + assert len(result["metal1"]) == 1 + assert isinstance(result["metal1"][0], int) + finally: + gmsh.clear() + gmsh.finalize() + + +class TestCollectMeshStats: + """Test collect_mesh_stats.""" + + def test_collect_stats_keys(self): + import gmsh + + from gsim.common.mesh.generator import collect_mesh_stats + + gmsh.initialize() + gmsh.model.add("test_stats") + kernel = gmsh.model.occ + + try: + # Create a simple box and mesh it + kernel.addBox(0, 0, 0, 10, 10, 10) + kernel.synchronize() + gmsh.model.mesh.generate(3) + + stats = collect_mesh_stats() + assert "nodes" in stats + assert "elements" in stats + assert "tetrahedra" in stats + assert "bbox" in stats + assert stats["nodes"] > 0 + assert stats["tetrahedra"] > 0 + finally: + gmsh.clear() + gmsh.finalize() diff --git a/tests/meep/test_mesh_sim.py b/tests/meep/test_mesh_sim.py new file mode 100644 index 0000000..27ed858 --- /dev/null +++ b/tests/meep/test_mesh_sim.py @@ -0,0 +1,367 @@ +"""Tests for MeepMeshSim and mesh config writer.""" + +from __future__ import annotations + +import json + +import pytest + +from gsim.common.mesh.types import MeshResult +from gsim.meep import FDTD, Domain, Material, MeepMeshSim, ModeSource +from gsim.meep.mesh_config import ( + _serialize_domain, + _serialize_materials, + _serialize_mesh_groups, + _serialize_solver, + _serialize_source, + write_mesh_config, +) + +# --------------------------------------------------------------------------- +# MeepMeshSim construction and attribute setting +# --------------------------------------------------------------------------- + + +class TestMeepMeshSimConstruction: + """Test MeepMeshSim creation and attribute updates.""" + + def test_defaults(self): + sim = MeepMeshSim() + assert sim.geometry.component is None + assert sim.materials == {} + assert sim.source.port is None + assert sim.monitors == [] + assert sim.domain.pml == 1.0 + assert sim.solver.resolution == 32 + + def test_callable_setters(self): + sim = MeepMeshSim() + sim.source(port="o1", wavelength=1.31) + assert sim.source.port == "o1" + assert sim.source.wavelength == 1.31 + + def test_domain_setter(self): + sim = MeepMeshSim() + sim.domain(pml=0.5, margin=0.3) + assert sim.domain.pml == 0.5 + assert sim.domain.margin == 0.3 + + def test_solver_setter(self): + sim = MeepMeshSim() + sim.solver(resolution=64) + assert sim.solver.resolution == 64 + + def test_materials_float_shorthand(self): + sim = MeepMeshSim() + sim.materials = {"si": 3.47, "SiO2": 1.44} + assert isinstance(sim.materials["si"], Material) + assert sim.materials["si"].n == 3.47 + assert sim.materials["SiO2"].n == 1.44 + + def test_materials_dict_shorthand(self): + sim = MeepMeshSim() + sim.materials = {"si": {"n": 3.47, "k": 0.01}} + assert isinstance(sim.materials["si"], Material) + assert sim.materials["si"].k == 0.01 + + def test_monitors_assignment(self): + sim = MeepMeshSim() + sim.monitors = ["o1", "o2", "o3"] + assert sim.monitors == ["o1", "o2", "o3"] + + +# --------------------------------------------------------------------------- +# MeepMeshSim error handling +# --------------------------------------------------------------------------- + + +class TestMeepMeshSimErrors: + """Test error conditions.""" + + def test_mesh_without_component_raises(self, tmp_path): + sim = MeepMeshSim() + with pytest.raises(ValueError, match="No component set"): + sim.mesh(tmp_path) + + def test_write_config_without_mesh_raises(self): + sim = MeepMeshSim() + with pytest.raises(ValueError, match="Call mesh"): + sim.write_config() + + def test_plot_mesh_without_mesh_raises(self): + sim = MeepMeshSim() + with pytest.raises(ValueError, match="Call mesh"): + sim.plot_mesh() + + +# --------------------------------------------------------------------------- +# Config serialization helpers +# --------------------------------------------------------------------------- + + +class TestSerializeMaterials: + """Test material serialization.""" + + def test_float_values(self): + result = _serialize_materials({"si": 3.47}) + assert result["si"]["refractive_index"] == 3.47 + assert result["si"]["extinction_coeff"] == 0.0 + + def test_material_values(self): + result = _serialize_materials({"si": Material(n=3.47, k=0.01)}) + assert result["si"]["refractive_index"] == 3.47 + assert result["si"]["extinction_coeff"] == 0.01 + + +class TestSerializeMeshGroups: + """Test mesh group serialization.""" + + def test_full_groups(self): + groups = { + "volumes": {"SiO2": {"phys_group": 1, "tags": [1, 2]}}, + "layer_volumes": {"core": {"phys_group": 3, "tags": [5], "material": "si"}}, + "outer_boundary": {"phys_group": 4, "tags": [10, 11]}, + } + result = _serialize_mesh_groups(groups) + assert result["volumes"]["SiO2"]["phys_group"] == 1 + assert result["volumes"]["SiO2"]["material"] == "SiO2" + assert "tags" not in result["volumes"]["SiO2"] + assert result["layer_volumes"]["core"]["material"] == "si" + assert result["outer_boundary"]["phys_group"] == 4 + + def test_empty_outer_boundary(self): + groups = {"volumes": {}, "layer_volumes": {}, "outer_boundary": {}} + result = _serialize_mesh_groups(groups) + assert "outer_boundary" not in result + + +class TestSerializeSource: + """Test source serialization.""" + + def test_default_source(self): + source = ModeSource(port="o1") + result = _serialize_source(source) + assert result["port"] == "o1" + assert result["wavelength"] == 1.55 + assert result["num_freqs"] == 11 + + +class TestSerializeDomain: + """Test domain serialization.""" + + def test_default_domain(self): + domain = Domain() + result = _serialize_domain(domain) + assert result["pml"] == 1.0 + assert result["margin_xy"] == 0.5 + + +class TestSerializeSolver: + """Test solver serialization.""" + + def test_default_solver(self): + solver = FDTD() + result = _serialize_solver(solver) + assert result["resolution"] == 32 + assert result["stopping"]["mode"] == "field_decay" + + +# --------------------------------------------------------------------------- +# write_mesh_config +# --------------------------------------------------------------------------- + + +class TestWriteMeshConfig: + """Test writing mesh_config.json.""" + + def test_writes_valid_json(self, tmp_path): + mesh_result = MeshResult( + mesh_path=tmp_path / "test.msh", + mesh_stats={"nodes": 1000, "tetrahedra": 5000}, + groups={ + "volumes": {"SiO2": {"phys_group": 1, "tags": [1]}}, + "layer_volumes": { + "core": {"phys_group": 2, "tags": [3], "material": "si"} + }, + "outer_boundary": {"phys_group": 3, "tags": [4]}, + }, + ) + + config_path = write_mesh_config( + mesh_result=mesh_result, + materials={"si": Material(n=3.47), "SiO2": Material(n=1.44)}, + source=ModeSource( + port="o1", wavelength=1.55, wavelength_span=0.1, num_freqs=11 + ), + monitors=["o1", "o2"], + domain=Domain(pml=1.0, margin=0.5), + solver=FDTD(resolution=32), + output_dir=tmp_path, + ) + + assert config_path.exists() + data = json.loads(config_path.read_text()) + + assert data["mesh_filename"] == "test.msh" + assert data["materials"]["si"]["refractive_index"] == 3.47 + assert data["source"]["port"] == "o1" + assert data["monitors"] == ["o1", "o2"] + assert data["domain"]["pml"] == 1.0 + assert data["solver"]["resolution"] == 32 + assert data["mesh_stats"]["nodes"] == 1000 + assert data["mesh_stats"]["tetrahedra"] == 5000 + assert data["mesh_groups"]["volumes"]["SiO2"]["phys_group"] == 1 + assert data["mesh_groups"]["volumes"]["SiO2"]["material"] == "SiO2" + assert data["mesh_groups"]["layer_volumes"]["core"]["material"] == "si" + + +# --------------------------------------------------------------------------- +# End-to-end: component + stack → .msh + mesh_config.json +# --------------------------------------------------------------------------- + + +class TestMeepMeshSimEndToEnd: + """End-to-end test with a real gdsfactory component and stack.""" + + @pytest.fixture + def sim_with_component(self): + """Create a MeepMeshSim with a simple waveguide component and stack.""" + import gdsfactory as gf + + from gsim.common.stack.extractor import Layer, LayerStack + + gf.gpdk.PDK.activate() + + # Simple waveguide + c = gf.Component() + c.add_polygon( + [(0, 0), (10, 0), (10, 0.5), (0, 0.5)], + layer=(1, 0), + ) + + stack = LayerStack(pdk_name="test") + stack.layers["core"] = Layer( + name="core", + gds_layer=(1, 0), + zmin=0.0, + zmax=0.22, + thickness=0.22, + material="si", + layer_type="dielectric", + ) + stack.dielectrics = [ + {"name": "box", "material": "SiO2", "zmin": -1.0, "zmax": 0.0}, + {"name": "clad", "material": "SiO2", "zmin": 0.0, "zmax": 1.0}, + ] + stack.materials = { + "si": {"type": "dielectric", "permittivity": 11.7}, + "SiO2": {"type": "dielectric", "permittivity": 2.1}, + } + + sim = MeepMeshSim() + sim.geometry(component=c, stack=stack) + sim.materials = {"si": 3.47, "SiO2": 1.44} + sim.source(port="o1", wavelength=1.55) + sim.monitors = ["o1", "o2"] + sim.domain(pml=1.0, margin=0.5) + sim.solver(resolution=32) + + return sim + + def test_mesh_generates_msh(self, sim_with_component, tmp_path): + sim = sim_with_component + result = sim.mesh(tmp_path) + + assert result.mesh_path.exists() + assert result.mesh_path.name == "mesh.msh" + assert result.mesh_stats.get("nodes", 0) > 0 + assert result.mesh_stats.get("tetrahedra", 0) > 0 + assert "volumes" in result.groups + assert "layer_volumes" in result.groups + + def test_write_config_after_mesh(self, sim_with_component, tmp_path): + sim = sim_with_component + sim.mesh(tmp_path) + config_path = sim.write_config() + + assert config_path.exists() + assert config_path.name == "mesh_config.json" + + data = json.loads(config_path.read_text()) + assert data["mesh_filename"] == "mesh.msh" + assert "si" in data["materials"] + assert data["source"]["wavelength"] == 1.55 + assert data["monitors"] == ["o1", "o2"] + assert data["solver"]["resolution"] == 32 + + def test_mesh_margin_defaults_to_domain(self, sim_with_component, tmp_path): + """Verify mesh margin defaults to domain.margin + domain.pml.""" + sim = sim_with_component + sim.domain(pml=0.5, margin=0.3) + result = sim.mesh(tmp_path) + # Should succeed — margin = 0.3 + 0.5 = 0.8 + assert result.mesh_path.exists() + + def test_dielectric_slices_have_consistent_xy(self, sim_with_component, tmp_path): + """Box and clad z-regions of SiO2 should have the same XY extent. + + Verifies that dielectric layers are not offset from each other + in the XY plane — regression test for the alternating-offset bug. + """ + import meshio + import numpy as np + + sim = sim_with_component + result = sim.mesh(tmp_path) + m = meshio.read(result.mesh_path) + + # Collect all tet nodes per physical group + tag_to_name = {tag: name for name, (tag, _) in m.field_data.items()} + group_points: dict[str, np.ndarray] = {} + for cells, phys in zip(m.cells, m.cell_data["gmsh:physical"], strict=False): + if cells.type != "tetra": + continue + for tag, name in tag_to_name.items(): + mask = phys == tag + if not np.any(mask): + continue + pts = m.points[cells.data[mask].ravel()] + if name in group_points: + group_points[name] = np.vstack([group_points[name], pts]) + else: + group_points[name] = pts + + # There should be both SiO2 (dielectric) and core (layer) groups + assert "SiO2" in group_points + assert "core" in group_points + + # SiO2 spans both box (-1→0) and clad (0→1) — check z range + sio2_pts = group_points["SiO2"] + tol = 0.01 + assert sio2_pts[:, 2].min() < -1.0 + tol, "SiO2 should extend below z=0 (box)" + assert sio2_pts[:, 2].max() > 1.0 - tol, "SiO2 should extend above z=0 (clad)" + + # Split SiO2 into box region (z < -0.1) and clad region (z > 0.3) + # to avoid the core transition zone + box_mask = sio2_pts[:, 2] < -0.1 + clad_mask = sio2_pts[:, 2] > 0.3 + assert np.any(box_mask), "No SiO2 nodes in box region" + assert np.any(clad_mask), "No SiO2 nodes in clad region" + + box_pts = sio2_pts[box_mask] + clad_pts = sio2_pts[clad_mask] + + # Compare XY bounding boxes of box vs clad regions + xy_tol = 0.1 # allow for mesh discretization + assert abs(box_pts[:, 0].min() - clad_pts[:, 0].min()) < xy_tol, ( + f"X min: box={box_pts[:, 0].min():.3f} vs clad={clad_pts[:, 0].min():.3f}" + ) + assert abs(box_pts[:, 0].max() - clad_pts[:, 0].max()) < xy_tol, ( + f"X max: box={box_pts[:, 0].max():.3f} vs clad={clad_pts[:, 0].max():.3f}" + ) + assert abs(box_pts[:, 1].min() - clad_pts[:, 1].min()) < xy_tol, ( + f"Y min: box={box_pts[:, 1].min():.3f} vs clad={clad_pts[:, 1].min():.3f}" + ) + assert abs(box_pts[:, 1].max() - clad_pts[:, 1].max()) < xy_tol, ( + f"Y max: box={box_pts[:, 1].max():.3f} vs clad={clad_pts[:, 1].max():.3f}" + )