Coverage for plugins/klayout-pex-plugin-gmsh/src/klayout_pex_plugins/gmsh/mesh.py: 97%

1# 

2# -------------------------------------------------------------------------------- 

3# SPDX-FileCopyrightText: 2026 Martin Jan Köhler and Harald Pretl 

4# Johannes Kepler University, Institute for Integrated Circuits. 

5# 

6# This file is part of klayout-pex-plugin-gmsh 

7# (see https://github.com/iic-jku/klayout-pex-plugin-gmsh). 

8# 

9# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <http://www.gnu.org/licenses/>. 

21# SPDX-License-Identifier: GPL-3.0-or-later 

22# -------------------------------------------------------------------------------- 

23# 

24 

25""" 

26Disjoint prisms → gmsh volume mesh for electrostatics. 

27 

28The prisms are extruded with OpenCASCADE and fragmented, so that touching solids 

29share their faces and the mesh is conformal. Conductors are then removed: only 

30the dielectrics are meshed, and each net is the set of faces where it touches a 

31dielectric. Physical groups carry the names: 

32 

33* one volume group per dielectric, the background included 

34* one surface group per terminal: each net, each floating conductor, the ground plane 

35* one surface group ``outer_boundary`` for the faces of the domain box 

36 

37Tags count from 1 within each kind, so they can serve directly as Palace 

38attributes and Elmer body / boundary indices. Coordinates are in µm. 

39""" 

40 

41from __future__ import annotations 

42 

43import contextlib 

44from dataclasses import asdict, dataclass, field 

45from importlib import metadata 

46import json 

47import math 

48import os 

49from typing import TYPE_CHECKING, Any, Dict, Iterator, List, Optional 

50 

51from klayout_pex.log import debug, info 

52from klayout_pex.plugin_api.v1 import ExportError, ExporterUnavailable 

53 

54if TYPE_CHECKING: 

55 from klayout_pex_protobuf.kpex.pex25d.pex25d_scene_pb2 import PEX25DScene 

56 

57 from .options import GmshMeshOptions 

58 from .solids import Body, SceneSolids 

59 

60OUTER_BOUNDARY = 'outer_boundary' 

61 

62LENGTH_UNIT_M = 1e-6 

63"""Mesh coordinates are in µm.""" 

64 

65MAX_SIZE_PER_EDGE_SIZE = 20.0 

66"""Default largest element, relative to the size at conductor edges.""" 

67 

68SIZE_GROWTH = 0.5 

69"""How fast elements grow away from conductor edges: size increase per distance.""" 

70 

71MAX_EDGE_SAMPLES = 2000 

72 

73 

74@dataclass(frozen=True) 

75class PhysicalGroup: 

76 dim: int 

77 tag: int 

78 name: str 

79 permittivity: Optional[float] = None 

80 """Relative permittivity of a dielectric volume.""" 

81 

82 

83@dataclass 

84class MeshResult: 

85 """What was written, and which physical group is what.""" 

86 

87 msh_path: str 

88 groups_path: str 

89 dielectrics: List[PhysicalGroup] = field(default_factory=list) 

90 terminals: List[PhysicalGroup] = field(default_factory=list) 

91 outer_boundary: Optional[PhysicalGroup] = None 

92 mesh_size_edge_um: Optional[float] = None 

93 """Element size used at conductor edges.""" 

94 

95 mesh_size_max_um: Optional[float] = None 

96 """Largest element size used.""" 

97 

98 node_count: int = 0 

99 tetrahedron_count: int = 0 

100 terminal_faces: List[int] = field(default_factory=list, repr=False) 

101 """The model's faces of all terminals, for refining the mesh near conductors.""" 

102 

103 @property 

104 def paths(self) -> List[str]: 

105 return [self.msh_path, self.groups_path] 

106 

107 

108def write_mesh(scene: PEX25DScene, 

109 settings: GmshMeshOptions, 

110 output_dir_path: str, 

111 prefix: str = '') -> MeshResult: 

112 """ 

113 Mesh ``scene`` and write ``<prefix>mesh.msh`` (gmsh 2.2, ASCII), plus 

114 ``<prefix>mesh.json`` naming each physical group. 

115 """ 

116 try: 

117 import gmsh 

118 except ImportError as exc: 

119 raise ExporterUnavailable(f"The gmsh Python module is not installed: {exc}") from exc 

120 from .solids import scene_solids 

121 

122 solids = scene_solids(scene, settings.field_margin_um) 

123 

124 os.makedirs(output_dir_path, exist_ok=True) 

125 result = MeshResult(msh_path=os.path.join(output_dir_path, f"{prefix}mesh.msh"), 

126 groups_path=os.path.join(output_dir_path, f"{prefix}mesh.json")) 

127 with _session(gmsh, settings.verbose): 

128 _build(gmsh, solids, result) 

129 _generate(gmsh, settings, solids, result) 

130 gmsh.write(result.msh_path) 

131 

132 with open(result.groups_path, 'w', encoding='utf-8') as file: 

133 json.dump(_groups_document(result), file, indent=2) 

134 file.write('\n') 

135 info(f"Meshed {len(result.dielectrics)} dielectric(s) and {len(result.terminals)} " 

136 f"terminal(s): {result.node_count} nodes, {result.tetrahedron_count} tetrahedra, " 

137 f"{result.mesh_size_edge_um} µm at conductor edges, up to {result.mesh_size_max_um} µm") 

138 return result 

139 

140 

141@contextlib.contextmanager 

142def _session(gmsh: Any, verbose: bool) -> Iterator[None]: 

143 """A model of our own; gmsh itself is finalized only if it was not running.""" 

144 owner = not gmsh.isInitialized() 

145 if owner: 

146 gmsh.initialize(readConfigFiles=False) 

147 gmsh.option.setNumber('General.Terminal', 1 if verbose else 0) 

148 gmsh.model.add('klayout-pex-plugin-gmsh') 

149 try: 

150 yield 

151 finally: 

152 gmsh.model.remove() 

153 if owner: 

154 gmsh.finalize() 

155 

156 

157def _build(gmsh: Any, solids: SceneSolids, result: MeshResult): 

158 from .solids import BodyKind 

159 

160 occ = gmsh.model.occ 

161 body_of_volume: Dict[int, Body] = {} 

162 for prism in solids.prisms: 

163 z = solids.um(prism.zlow) 

164 height = solids.um(prism.zhigh - prism.zlow) 

165 for polygon in prism.region.each(): 

166 for volume in _extrude(occ, polygon, z, height, solids.grid_um): 

167 body_of_volume[volume] = prism.body 

168 debug(f"{len(body_of_volume)} volumes before fragmenting") 

169 

170 inputs = [(3, volume) for volume in body_of_volume] 

171 if len(inputs) > 1: 

172 # The prisms are disjoint, so each fragment comes from exactly one of them. 

173 _, children = occ.fragment(inputs[:1], inputs[1:]) 

174 fragments: Dict[int, Body] = {} 

175 for (_, parent), parent_children in zip(inputs, children): 

176 for _, child in parent_children: 

177 if fragments.setdefault(child, body_of_volume[parent]) != body_of_volume[parent]: 

178 raise ExportError(f"Solids of '{fragments[child].name}' and " 

179 f"'{body_of_volume[parent].name}' overlap") 

180 body_of_volume = fragments 

181 occ.synchronize() 

182 

183 terminal_faces: Dict[Body, List[int]] = {} 

184 outer_faces: List[int] = [] 

185 for _, face in gmsh.model.getEntities(2): 

186 volumes, _ = gmsh.model.getAdjacencies(2, face) 

187 bodies = [body_of_volume[volume] for volume in volumes] 

188 conductors = [body for body in bodies if body.kind == BodyKind.CONDUCTOR] 

189 if len(bodies) == 1 and not conductors: 

190 outer_faces.append(face) 

191 elif len(bodies) == 2 and len(conductors) == 1: 

192 terminal_faces.setdefault(conductors[0], []).append(face) 

193 elif len(conductors) == 2 and conductors[0] != conductors[1]: 

194 raise ExportError(f"Conductors '{conductors[0].name}' and " 

195 f"'{conductors[1].name}' touch, which shorts them") 

196 

197 # Only the dielectrics are meshed. Removing the conductors keeps their faces 

198 # as the dielectrics' boundaries; faces left without a volume are dropped. 

199 occ.remove([(3, volume) for volume, body in body_of_volume.items() 

200 if body.kind == BodyKind.CONDUCTOR]) 

201 occ.synchronize() 

202 dangling = [(2, face) for _, face in gmsh.model.getEntities(2) 

203 if len(gmsh.model.getAdjacencies(2, face)[0]) == 0] 

204 if dangling: 

205 occ.remove(dangling, recursive=True) 

206 occ.synchronize() 

207 

208 for body in solids.bodies: 

209 volumes = [volume for volume, owner in body_of_volume.items() 

210 if owner == body and body.kind == BodyKind.DIELECTRIC] 

211 if volumes: 

212 result.dielectrics.append(_physical_group( 

213 gmsh, 3, volumes, len(result.dielectrics) + 1, body.name, body.permittivity)) 

214 for body in solids.bodies: 

215 if terminal_faces.get(body): 

216 result.terminals.append(_physical_group( 

217 gmsh, 2, terminal_faces[body], len(result.terminals) + 1, body.name)) 

218 result.terminal_faces.extend(terminal_faces[body]) 

219 if any(group.name == OUTER_BOUNDARY for group in result.terminals): 

220 raise ExportError(f"A net is named '{OUTER_BOUNDARY}', the name of the " 

221 f"domain boundary group") 

222 if outer_faces: 

223 result.outer_boundary = _physical_group( 

224 gmsh, 2, outer_faces, len(result.terminals) + 1, OUTER_BOUNDARY) 

225 

226 

227def _extrude(occ: Any, polygon: Any, z: float, height: float, scale: float) -> List[int]: 

228 loops = [_curve_loop(occ, polygon.each_point_hull(), z, scale)] 

229 loops.extend(_curve_loop(occ, polygon.each_point_hole(hole), z, scale) 

230 for hole in range(polygon.holes())) 

231 surface = occ.addPlaneSurface(loops) 

232 return [tag for dim, tag in occ.extrude([(2, surface)], 0, 0, height) if dim == 3] 

233 

234 

235def _curve_loop(occ: Any, points: Any, z: float, scale: float) -> int: 

236 # gmsh wants the holes wound like the hull, KLayout winds them the other 

237 # way; opposite windings make OCC add a hole's area instead of cutting it. 

238 points = [(point.x, point.y) for point in points] 

239 twice_area = sum(x0 * y1 - x1 * y0 

240 for (x0, y0), (x1, y1) in zip(points, points[1:] + points[:1])) 

241 if twice_area < 0: 

242 points.reverse() 

243 tags = [occ.addPoint(x * scale, y * scale, z) for x, y in points] 

244 lines = [occ.addLine(start, end) for start, end in zip(tags, tags[1:] + tags[:1])] 

245 return occ.addCurveLoop(lines) 

246 

247 

248def _physical_group(gmsh: Any, dim: int, entities: List[int], tag: int, name: str, 

249 permittivity: Optional[float] = None) -> PhysicalGroup: 

250 gmsh.model.addPhysicalGroup(dim, entities, tag) 

251 gmsh.model.setPhysicalName(dim, tag, name) 

252 return PhysicalGroup(dim, tag, name, permittivity) 

253 

254 

255def _generate(gmsh: Any, settings: GmshMeshOptions, solids: SceneSolids, result: MeshResult): 

256 gmsh.option.setNumber('Mesh.MshFileVersion', 2.2) 

257 gmsh.option.setNumber('Mesh.Binary', 0) 

258 gmsh.option.setNumber('Mesh.SaveAll', 0) 

259 

260 edge_size = settings.mesh_size_edge_um 

261 if edge_size is None and solids.conductor_thickness: 

262 edge_size = solids.um(solids.conductor_thickness) 

263 max_size = settings.mesh_size_max_um 

264 if max_size is None and edge_size is not None: 

265 max_size = edge_size * MAX_SIZE_PER_EDGE_SIZE 

266 if edge_size is not None and result.terminal_faces: 

267 _refine_at_conductor_edges(gmsh, result.terminal_faces, edge_size, max_size) 

268 result.mesh_size_edge_um = edge_size 

269 result.mesh_size_max_um = max_size 

270 

271 if settings.mesh_size_min_um is not None: 

272 gmsh.option.setNumber('Mesh.MeshSizeMin', settings.mesh_size_min_um) 

273 if max_size is not None: 

274 gmsh.option.setNumber('Mesh.MeshSizeMax', max_size) 

275 gmsh.model.mesh.generate(3) 

276 

277 result.node_count = len(gmsh.model.mesh.getNodes()[0]) 

278 element_types, element_tags, _ = gmsh.model.mesh.getElements(3) 

279 result.tetrahedron_count = sum(len(tags) for element_type, tags 

280 in zip(element_types, element_tags) if element_type == 4) 

281 

282 

283def _refine_at_conductor_edges(gmsh: Any, faces: List[int], edge_size: float, 

284 max_size: float): 

285 """ 

286 The field concentrates at conductor edges: elements are ``edge_size`` there and 

287 grow with the distance from the nearest edge, up to ``max_size``. 

288 """ 

289 curves = sorted({curve for _, curve in gmsh.model.getBoundary( 

290 [(2, face) for face in faces], combined=False, oriented=False)}) 

291 longest = max(gmsh.model.occ.getMass(1, curve) for curve in curves) 

292 

293 fields = gmsh.model.mesh.field 

294 distance = fields.add('Distance') 

295 fields.setNumbers(distance, 'CurvesList', curves) 

296 # Sample the longest edge at the element size, so the distance is not bumpy. 

297 fields.setNumber(distance, 'Sampling', min(MAX_EDGE_SAMPLES, max( 

298 2, math.ceil(longest / edge_size) + 1))) 

299 threshold = fields.add('Threshold') 

300 fields.setNumber(threshold, 'InField', distance) 

301 fields.setNumber(threshold, 'SizeMin', edge_size) 

302 fields.setNumber(threshold, 'SizeMax', max_size) 

303 fields.setNumber(threshold, 'DistMin', 0.0) 

304 fields.setNumber(threshold, 'DistMax', (max_size - edge_size) / SIZE_GROWTH) 

305 fields.setAsBackgroundMesh(threshold) 

306 

307 # The field alone decides; sizes from points or boundaries would override it. 

308 gmsh.option.setNumber('Mesh.MeshSizeExtendFromBoundary', 0) 

309 gmsh.option.setNumber('Mesh.MeshSizeFromPoints', 0) 

310 gmsh.option.setNumber('Mesh.MeshSizeFromCurvature', 0) 

311 

312 

313def _groups_document(result: MeshResult) -> Dict[str, Any]: 

314 def entry(group: PhysicalGroup) -> Dict[str, Any]: 

315 document = asdict(group) 

316 document.pop('dim') 

317 if document['permittivity'] is None: 

318 document.pop('permittivity') 

319 return document 

320 

321 return { 

322 'generator': f"klayout-pex-plugin-gmsh {metadata.version('klayout-pex-plugin-gmsh')}", 

323 'mesh': os.path.basename(result.msh_path), 

324 'length_unit_m': LENGTH_UNIT_M, 

325 'mesh_size_edge_um': result.mesh_size_edge_um, 

326 'mesh_size_max_um': result.mesh_size_max_um, 

327 'dielectrics': [entry(group) for group in result.dielectrics], 

328 'terminals': [entry(group) for group in result.terminals], 

329 'outer_boundary': entry(result.outer_boundary) if result.outer_boundary else None, 

330 }