Coverage for klayout_pex/fastercap/fastercap_model_generator.py: 96%
690 statements
« prev ^ index » next coverage.py v7.17.2, created at 2026-11-06 13:14 +0000
1#
2# --------------------------------------------------------------------------------
3# SPDX-FileCopyrightText: 2024-2025 Martin Jan Köhler and Harald Pretl
4# Johannes Kepler University, Institute for Integrated Circuits.
5#
6# This file is part of KPEX
7# (see https://github.com/iic-jku/klayout-pex).
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#
25# A class providing a service for building FastCap2 or FasterCap models
26#
27# This class is used the following way:
28#
29# 1) Create a FasterCapModelBuilder object
30# Specify the default k value which is the k assumed for "empty space".
31# You can also specify a maximum area and the "b" parameter for the
32# triangulation. The b parameter corresponds to the minimum angle
33# and should be <=1 (b=sin(min_angle)*2).
34# I.e. b=1 -> min_angle=30 deg, b=0.5 -> min_angle~14.5 deg.
35#
36# 2) Add material definitions for the dielectrics
37# Each material definition consists of a k value and
38# a material name.
39#
40# 3) Add layers in the 2.5d view fashion
41# Each layer is a sheet in 3d space that is extruded in vertical
42# direction with the given start and stop z (or height)
43# The layer must be a DRC::Layer or RBA::Region object.
44#
45# Layers can be added in two ways:
46#
47# * As conductors: specify the net name
48#
49# * As dielectric layer: specify the material name
50#
51# The layers can intersect. The package resolves intersections
52# based on priority: conductors first, dielectrics according to
53# their position in the "materials" definition (first entries have
54# higher prio)
55#
56# 4) Generate a 3d model using "generate"
57# This method returns an object you can use to generate STL files
58# or FastCap files.
61from __future__ import annotations
63import base64
64from collections import defaultdict
65import hashlib
66import os
67from typing import *
68from dataclasses import dataclass
69from functools import reduce
70import math
72import klayout.db as kdb
74from ..log import (
75 debug,
76 info,
77 warning,
78 error,
79 subproc
80)
83@dataclass
84class FasterCapModelBuilder:
85 dbu: float
86 """Database unit"""
88 k_void: float
89 """Default dielectric of 'empty space'"""
91 delaunay_amax: float
92 """Maximum area parameter for the Delaunay triangulation"""
94 delaunay_b: float
95 """
96 The delaunay_b parameter for the Delaunay triangulation
97 corresponds to the minimum angle
98 and should be <=1 (b=sin(min_angle)*2).
99 I.e. b=1 -> min_angle=30 deg, b=0.5 -> min_angle~14.5 deg.
100 """
102 def __init__(self,
103 dbu: float,
104 k_void: float,
105 delaunay_amax: float = 0.0,
106 delaunay_b: float = 1.0,
107 ):
108 self.dbu = dbu
109 self.k_void = k_void
110 self.delaunay_amax = delaunay_amax
111 self.delaunay_b = delaunay_b
113 self.materials: Dict[str, float] = {}
114 self.net_names: List[str] = []
115 self.group_by_conductor: Dict[str, str] = {}
117 # layer, zstart, zstop
118 self.clayers: Dict[str, List[Tuple[kdb.Region, float, float]]] = {}
119 self.dlayers: Dict[str, List[Tuple[kdb.Region, float, float]]] = {}
121 info(f"DBU: {'%.12g' % self.dbu}")
122 info(f"Delaunay b: {'%.12g' % self.delaunay_b}")
123 info(f"Delaunay area_max: {'%.12g' % self.delaunay_amax}")
125 def add_material(self, name: str, k: float):
126 self.materials[name] = k
128 def add_dielectric(self,
129 material_name: str,
130 layer: kdb.Region,
131 z: float,
132 height: float):
133 if hasattr(layer, 'data'):
134 layer = layer.data
135 self._add_layer(name=material_name, layer=layer, is_dielectric=True, z=z, height=height)
137 def add_conductor(self,
138 net_name: str,
139 layer: kdb.Region,
140 z: float,
141 height: float,
142 group_name: Optional[str] = None):
143 """
144 :param net_name: the name of the conductor
145 :param group_name: of the conductors that touch, as they are parts of one net, which have no surfaces
146 between them (see FasterCapModelGenerator), by default the conductor is on its own
147 """
148 if hasattr(layer, 'data'):
149 layer = layer.data
150 self.group_by_conductor[net_name] = group_name or net_name
151 self._add_layer(name=net_name, layer=layer, is_dielectric=False, z=z, height=height)
153 def _norm2z(self, z: float) -> float:
154 return z * self.dbu
156 def _z2norm(self, z: float) -> float:
157 return math.floor(z / self.dbu + 1e-6)
159 def _add_layer(self,
160 name: str,
161 layer: kdb.Region,
162 z: float,
163 height: float,
164 is_dielectric: bool):
165 if is_dielectric and name not in self.materials:
166 raise ValueError(f"Unknown material {name} - did you use 'add_material'?")
168 zstart: float = z
169 zstop: float = zstart + height
171 if is_dielectric:
172 if name not in self.dlayers:
173 self.dlayers[name] = []
174 self.dlayers[name].append((layer, self._z2norm(zstart), self._z2norm(zstop)))
175 else:
176 if name not in self.clayers:
177 self.clayers[name] = []
178 self.clayers[name].append((layer, self._z2norm(zstart), self._z2norm(zstop)))
180 def generate(self) -> Optional[FasterCapModelGenerator]:
181 z: List[float] = []
182 for ll in (self.dlayers, self.clayers):
183 for k, v in ll.items():
184 for l in v:
185 z.extend((l[1], l[2]))
186 z = sorted([*{*z}]) # sort & uniq
187 if len(z) == 0:
188 return None
190 gen = FasterCapModelGenerator(dbu=self.dbu,
191 k_void=self.k_void,
192 delaunay_amax=self.delaunay_amax,
193 delaunay_b=self.delaunay_b,
194 materials=self.materials,
195 net_names=list(self.clayers.keys()),
196 group_by_conductor=self.group_by_conductor)
197 for zcurr in z:
198 gen.next_z(self._norm2z(zcurr))
200 for nn, v in self.clayers.items():
201 for l in v:
202 if l[1] <= zcurr < l[2]:
203 gen.add_in(name=f"+{nn}", layer=l[0])
204 if l[1] < zcurr <= l[2]:
205 gen.add_out(name=f"+{nn}", layer=l[0])
206 for mn, v in self.dlayers.items():
207 for l in v:
208 if l[1] <= zcurr < l[2]:
209 gen.add_in(name=f"-{mn}", layer=l[0])
210 if l[1] < zcurr <= l[2]:
211 gen.add_out(name=f"-{mn}", layer=l[0])
213 gen.finish_z()
215 gen.finalize()
216 return gen
219@dataclass(frozen=True)
220class HDielKey:
221 outside: Optional[str]
222 inside: Optional[str]
224 def __str__(self) -> str:
225 return f"{self.outside or 'void'} <-> {self.inside or 'void'}"
227 @property
228 def topic(self) -> str:
229 return 'dielectric'
231 def reversed(self) -> HDielKey:
232 return HDielKey(self.inside, self.outside)
235@dataclass(frozen=True)
236class HCondKey:
237 net_name: str
238 outside: Optional[str]
240 def __str__(self) -> str:
241 return f"{self.outside or 'void'} <-> {self.net_name}"
243 @property
244 def topic(self) -> str:
245 return 'conductor'
248@dataclass(frozen=True)
249class VKey:
250 kk: HDielKey | HCondKey
251 p0: kdb.DPoint
252 de: kdb.DVector
255@dataclass(frozen=True)
256class Point:
257 x: float
258 y: float
259 z: float
261 def __sub__(self, other: Point) -> Point:
262 return Point(self.x - other.x, self.y - other.y, self.z - other.z)
264 def sq_length(self) -> float:
265 return self.x**2 + self.y**2 + self.z**2
267 def to_fastcap(self) -> str:
268 return '%.12g %.12g %.12g' % (self.x, self.y, self.z)
271def vector_product(a: Point, b: Point) -> Point:
272 vp = Point(
273 a.y * b.z - a.z * b.y,
274 a.z * b.x - a.x * b.z,
275 a.x * b.y - a.y * b.x
276 )
277 return vp
280def dot_product(a: Point, b: Point) -> float:
281 dp = a.x * b.x + a.y * b.y + a.z * b.z
282 return dp
285@dataclass(frozen=True)
286class Triangle:
287 p0: Point
288 p1: Point
289 p2: Point
291 def reversed(self) -> Triangle:
292 return Triangle(self.p2, self.p1, self.p0)
294 def outside_reference_point(self) -> Point:
295 v1 = self.p1 - self.p0
296 v2 = self.p2 - self.p0
297 vp = Point(v1.y * v2.z - v1.z * v2.y,
298 -v1.x * v2.z + v1.z * v2.x,
299 v1.x * v2.y - v1.y * v2.x)
300 vp_abs = math.sqrt(vp.x ** 2 + vp.y ** 2 + vp.z ** 2)
301 rp = Point(self.p0.x + vp.x / vp_abs,
302 self.p0.y + vp.y / vp_abs,
303 self.p0.z + vp.z / vp_abs)
304 return rp
306 def to_fastcap(self) -> str:
307 return ' '.join([p.to_fastcap() for p in (self.p0, self.p1, self.p2)])
309 def __len__(self):
310 return 3
312 def __getitem__(self, i) -> Point:
313 match i:
314 case 0: return self.p0
315 case 1: return self.p1
316 case 2: return self.p2
317 case _: raise IndexError("list index out of range")
320@dataclass(frozen=True)
321class Edge:
322 p0: Point
323 p1: Point
325 def vector_of_edge(self) -> Point:
326 return Point(
327 self.p1.x - self.p0.x,
328 self.p1.y - self.p0.y,
329 self.p1.z - self.p0.z
330 )
332 def reversed(self) -> Edge:
333 return Edge(self.p1, self.p0)
336@dataclass
337class FasterCapModelGenerator:
338 dbu: float
339 """Database unit"""
341 k_void: float
342 """Default dielectric of 'empty space'"""
344 delaunay_amax: float
345 """Maximum area parameter for the Delaunay triangulation"""
347 delaunay_b: float
348 """
349 The delaunay_b parameter for the Delaunay triangulation
350 corresponds to the minimum angle
351 and should be <=1 (b=sin(min_angle)*2).
352 I.e. b=1 -> min_angle=30 deg, b=0.5 -> min_angle~14.5 deg.
353 """
355 materials: Dict[str, float]
356 """Maps material name to dielectric k"""
358 net_names: List[str]
360 group_by_conductor: Dict[str, str]
361 """The conductors that are parts of a net by the name of the net, the others are on their own"""
363 def __init__(self,
364 dbu: float,
365 k_void: float,
366 delaunay_amax: float,
367 delaunay_b: float,
368 materials: Dict[str, float],
369 net_names: List[str],
370 group_by_conductor: Optional[Dict[str, str]] = None):
371 self.k_void = k_void
372 self.delaunay_amax = delaunay_amax
373 self.delaunay_b = delaunay_b
374 self.dbu = dbu
375 self.materials = materials
376 self.net_names = net_names
377 self.group_by_conductor = group_by_conductor or {}
379 # NOTE: without surfaces between the conductors of a group, whose coincident panels would be singular.
380 # The charge they'd have is the one between the parts of a net, which no netlist has
381 conductors_by_group: Dict[str, List[str]] = defaultdict(list)
382 for conductor_name, group_name in self.group_by_conductor.items():
383 conductors_by_group[group_name].append(conductor_name)
384 self.conductors_by_group = conductors_by_group
386 self.z: Optional[float] = None
387 self.zz: Optional[float] = None
388 self.layers_in: Dict[str, kdb.Region] = {}
389 self.layers_out: Dict[str, kdb.Region] = {}
390 self.state: Dict[str, kdb.Region] = {}
391 self.current: Dict[str, List[kdb.Region]] = {}
392 self.diel_data: Dict[HDielKey, List[Triangle]] = {}
393 self.diel_vdata: Dict[VKey, kdb.Region] = {}
394 self.cond_data: Dict[HCondKey, List[Triangle]] = {}
395 self.cond_vdata: Dict[VKey, kdb.Region] = {}
397 def reset(self):
398 self.layers_in = {}
399 self.layers_out = {}
401 def other_conductors_of_group(self, conductor_name: str) -> List[str]:
402 group_name = self.group_by_conductor.get(conductor_name, conductor_name)
403 return [name for name in self.conductors_by_group.get(group_name, ()) if name != conductor_name]
405 def add_in(self, name: str, layer: kdb.Region):
406 debug(f"add_in: {name} -> {layer}")
407 if name not in self.layers_in:
408 # NOTE: Calling kdb.Region([]) with the empty list enforces, that a flat layer is created.
409 # The following "+=" packs the polynomials into this flat layer and
410 # does not copy the hierarchical layer from the LayoutToNetlist database.
411 self.layers_in[name] = kdb.Region([])
412 self.layers_in[name] += layer
414 def add_out(self, name: str, layer: kdb.Region):
415 debug(f"add_out: {name} -> {layer}")
416 if name not in self.layers_out:
417 # NOTE: Calling kdb.Region([]) with the empty list enforces, that a flat layer is created.
418 # The following "+=" packs the polynomials into this flat layer and
419 # does not copy the hierarchical layer from the LayoutToNetlist database.
420 self.layers_out[name] = kdb.Region([])
421 self.layers_out[name] += layer
423 def finish_z(self):
424 debug(f"Finishing layer z={self.z}")
426 din: Dict[str, kdb.Region] = {}
427 dout: Dict[str, kdb.Region] = {}
428 all_in = kdb.Region()
429 all_out = kdb.Region()
430 all = kdb.Region()
431 all_cin: Optional[kdb.Region] = None
432 all_cout: Optional[kdb.Region] = None
434 for names, prefix in ((self.net_names, '+'), (self.materials.keys(), '-')):
435 for nn in names:
436 mk = prefix + nn
438 # compute merged events
439 if mk not in self.current:
440 self.current[mk] = []
441 current_before = self.current[mk][0].dup() if len(self.current[mk]) >= 1 else kdb.Region()
442 lin, lout, current = self._merge_events(pyra=self.current[mk],
443 lin=self.layers_in.get(mk, None),
444 lout=self.layers_out.get(mk, None))
445 debug(f"Merged events & status for {mk}:")
446 debug(f" in = {lin}")
447 debug(f" out = {lout}")
448 debug(f" state = {current}")
450 if mk not in self.state:
451 self.state[mk] = kdb.Region()
453 # legalize in and out events
454 lin_org = lin.dup()
455 lout_org = lout.dup()
456 lout &= self.state[mk]
457 lin -= all
458 lout += current & all_in
459 lin += current_before & all_out
460 lin -= lout_org
461 lout -= lin_org
463 # tracks the legalized horizontal cuts
464 self.state[mk] += lin
465 self.state[mk] -= lout
467 din[mk] = lin
468 dout[mk] = lout
470 debug(f"Legalized events & status for '{mk}':")
471 debug(f" in = {din[mk]}")
472 debug(f" out = {dout[mk]}")
473 debug(f" state = {self.state[mk]}")
475 all_in += lin
476 all_out += lout
477 all += self.state[mk]
479 if prefix == '+':
480 all_cin = all_in.dup()
481 all_cout = all_out.dup()
483 debug(f"All conductor region in: {all_cin}")
484 debug(f"All conductor region out: {all_cout}")
486 # check whether states are separated
487 a = reduce(lambda x, y: x+y, self.state.values())
488 for k, s in self.state.items():
489 r: kdb.Region = s - a
490 if not r.is_empty():
491 error(f"State region of {k} ({s}) is not contained entirely "
492 f"in remaining all state region ({a}) - this means there is an overlap")
493 a -= s
495 # Now we have legalized the in and out events
496 for mni in self.materials.keys():
497 lin = din.get(f"-{mni}", None)
498 if lin:
499 lin = lin.dup()
500 lin -= all_cout # handled with the conductor
501 for mno in self.materials.keys():
502 lout = dout.get(f"-{mno}", None)
503 if lout:
504 d: kdb.Region = lout & lin
505 if not d.is_empty():
506 self.generate_hdiel(below=mno, above=mni, layer=d)
507 lin -= lout
508 if not lin.is_empty():
509 self.generate_hdiel(below=None, above=mni, layer=lin)
511 for mno in self.materials.keys():
512 lout = dout.get(f"-{mno}", None)
513 if lout:
514 lout = lout.dup()
515 lout -= all_cin # handled with the conductor
516 for mni in self.materials.keys():
517 lin = din.get(f"-{mni}", None)
518 if lin:
519 lout -= lin
520 if not lout.is_empty():
521 self.generate_hdiel(below=mno, above=None, layer=lout)
523 for nn in self.net_names:
524 lin = din.get(f"+{nn}", None)
525 if lin:
526 lin = lin.dup()
527 for other in self.other_conductors_of_group(nn):
528 lout = dout.get(f"+{other}", None)
529 if lout:
530 lin -= lout
531 for mno in self.materials.keys():
532 lout = dout.get(f"-{mno}", None)
533 if lout:
534 d = lout & lin
535 if not d.is_empty():
536 self.generate_hcond_in(net_name=nn, below=mno, layer=d)
537 lin -= lout
538 if not lin.is_empty():
539 self.generate_hcond_in(net_name=nn, below=None, layer=lin)
541 for nn in self.net_names:
542 lout = dout.get(f"+{nn}", None)
543 if lout:
544 lout = lout.dup()
545 lout -= all_cin # handled with the conductor
546 for mni in self.materials.keys():
547 lin = din.get(f"-{mni}", None)
548 if lin:
549 d = lout & lin
550 if not d.is_empty():
551 self.generate_hcond_out(net_name=nn, above=mni, layer=d)
552 lout -= lin
553 if not lout.is_empty():
554 self.generate_hcond_out(net_name=nn, above=None, layer=lout)
556 def next_z(self, z: float):
557 debug(f"Next layer {z}")
559 self.reset()
561 if self.z is None:
562 self.z = z
563 return
565 self.zz = z
567 all_cond = kdb.Region()
568 for nn in self.net_names:
569 mk = f"+{nn}"
570 if mk in self.state:
571 all_cond += self.state[mk]
572 all_cond = all_cond.edges()
574 for i, mni in enumerate(self.materials):
575 linside = self.state.get(f"-{mni}", None)
576 if linside:
577 linside = linside.edges()
578 linside -= all_cond # handled with the conductor
579 for o, mno in enumerate(self.materials):
580 if i != o:
581 loutside = self.state.get(f"-{mno}", None)
582 if loutside:
583 loutside = loutside.edges()
584 if o > i:
585 d = loutside & linside
586 for e in d:
587 # NOTE: we need to swap points as we started from "outside"
588 self.generate_vdiel(outside=mno, inside=mni, edge=e.swapped_points())
589 linside -= loutside
591 for e in linside:
592 self.generate_vdiel(outside=None, inside=mni, edge=e)
594 for nn in self.net_names:
595 mk = f"+{nn}"
596 linside = self.state.get(mk, None)
597 if linside:
598 linside = linside.edges()
599 for other in self.other_conductors_of_group(nn):
600 other_state = self.state.get(f"+{other}", None)
601 if other_state:
602 linside -= other_state.edges()
603 for mno in self.materials:
604 loutside = self.state.get(f"-{mno}", None)
605 if loutside:
606 loutside = loutside.edges()
607 d = loutside & linside
608 for e in d:
609 # NOTE: we need to swap points as we started from "outside"
610 self.generate_vcond(net_name=nn, outside=mno, edge=e.swapped_points())
611 linside -= loutside
612 for e in linside:
613 self.generate_vcond(net_name=nn, outside=None, edge=e)
615 self.z = z
617 def generate_hdiel(self,
618 below: Optional[str],
619 above: Optional[str],
620 layer: kdb.Region):
621 k = HDielKey(below, above)
622 debug(f"Generating horizontal dielectric surface {k} as {layer}")
623 if k not in self.diel_data:
624 self.diel_data[k] = []
625 data = self.diel_data[k]
627 for t in layer.delaunay(self.delaunay_amax / self.dbu ** 2, self.delaunay_b):
628 # NOTE: normal is facing downwards (to "below")
629 pl = list(map(lambda pt: Point(pt.x * self.dbu, pt.y * self.dbu, self.z),
630 t.each_point_hull()))
631 tri = Triangle(*pl)
632 data.append(tri)
633 debug(f" {tri}")
635 def generate_v_surface(self,
636 kk: HDielKey | HCondKey,
637 edge: kdb.Edge) -> Tuple[VKey, kdb.Box]:
638 debug(f"Generating vertical {kk.topic} surface {kk} with edge {edge}")
640 el = math.sqrt(edge.sq_length())
641 de = kdb.DVector(edge.d().x / el, edge.d().y / el)
642 ne = kdb.DVector(edge.d().y / el, -edge.d().x / el)
643 p0 = ne * ne.sprod(kdb.DPoint(edge.p1) - kdb.DPoint()) + kdb.DPoint()
644 x1 = (edge.p1 - p0).sprod(de)
645 x2 = (edge.p2 - p0).sprod(de)
647 key = VKey(kk, p0, de)
648 surface = kdb.Box(x1,
649 math.floor(self.z / self.dbu + 0.5),
650 x2,
651 math.floor(self.zz / self.dbu + 0.5))
652 return key, surface
654 def generate_vdiel(self,
655 outside: Optional[str],
656 inside: Optional[str],
657 edge: kdb.Edge):
658 if edge.is_degenerate():
659 return
661 key, surface = self.generate_v_surface(HDielKey(outside, inside), edge)
662 if key not in self.diel_vdata:
663 self.diel_vdata[key] = kdb.Region()
665 self.diel_vdata[key].insert(surface)
667 def generate_hcond_in(self,
668 net_name: str,
669 below: Optional[str],
670 layer: kdb.Region):
671 k = HCondKey(net_name, below)
672 debug(f"Generating horizontal bottom conductor surface {k} as {layer}")
674 if k not in self.cond_data:
675 self.cond_data[k] = []
676 data = self.cond_data[k]
678 for t in layer.delaunay(self.delaunay_amax / self.dbu ** 2, self.delaunay_b):
679 # NOTE: normal is facing downwards (to "below")
680 pl = list(map(lambda pt: Point(pt.x * self.dbu, pt.y * self.dbu, self.z),
681 t.each_point_hull()))
682 tri = Triangle(*pl)
683 data.append(tri)
684 debug(f" {tri}")
686 def generate_hcond_out(self,
687 net_name: str,
688 above: Optional[str],
689 layer: kdb.Region):
690 k = HCondKey(net_name, above)
691 debug(f"Generating horizontal top conductor surface {k} as {layer}")
693 if k not in self.cond_data:
694 self.cond_data[k] = []
695 data = self.cond_data[k]
697 for t in layer.delaunay(self.delaunay_amax / self.dbu ** 2, self.delaunay_b):
698 # NOTE: normal is facing downwards (into conductor)
699 pl = list(map(lambda pt: Point(pt.x * self.dbu, pt.y * self.dbu, self.z),
700 t.each_point_hull()))
701 tri = Triangle(*pl)
702 # now it is facing outside (to "above")
703 tri = tri.reversed()
704 data.append(tri)
705 debug(f" {tri}")
707 def generate_vcond(self,
708 net_name: str,
709 outside: Optional[str],
710 edge: kdb.Edge):
711 if edge.is_degenerate():
712 return
714 key, surface = self.generate_v_surface(HCondKey(net_name, outside), edge)
715 if key not in self.cond_vdata:
716 self.cond_vdata[key] = kdb.Region()
718 self.cond_vdata[key].insert(surface)
720 def triangulate(self, p0: kdb.DPoint, de: kdb.DVector, region: kdb.Region, data: List[Triangle]):
721 def convert_point(pt: kdb.Point) -> Point:
722 pxy = (p0 + de * pt.x) * self.dbu
723 pz = pt.y * self.dbu
724 return Point(pxy.x, pxy.y, pz)
726 for t in region.delaunay(self.delaunay_amax / self.dbu ** 2, self.delaunay_b):
727 # NOTE: normal is facing outwards (to "left")
728 pl = list(map(convert_point, t.each_point_hull()))
729 tri = Triangle(*pl)
730 # now it is facing outside (to "above")
731 data.append(tri)
732 debug(f" {tri}")
734 def finalize(self):
735 for k, r in self.diel_vdata.items():
736 debug(f"Finishing vertical dielectric plane {k.kk} at {k.p0}/{k.de}")
738 if k.kk not in self.diel_data:
739 self.diel_data[k.kk] = []
740 data = self.diel_data[k.kk]
742 self.triangulate(p0=k.p0, de=k.de, region=r, data=data)
744 for k, r in self.cond_vdata.items():
745 debug(f"Finishing vertical conductor plane {k.kk} at {k.p0} / {k.de}")
747 if k.kk not in self.cond_data:
748 self.cond_data[k.kk] = []
749 data = self.cond_data[k.kk]
751 self.triangulate(p0=k.p0, de=k.de, region=r, data=data)
753 dk: Dict[HDielKey, List[Triangle]] = {}
755 for k in self.diel_data.keys():
756 kk = k.reversed()
757 if kk not in dk:
758 dk[k] = []
759 else:
760 debug(f"Combining dielectric surfaces {kk} with reverse")
762 for k, v in self.diel_data.items():
763 kk = k.reversed()
764 if kk in dk:
765 dk[kk] += list(map(lambda t: t.reversed(), v))
766 else:
767 dk[k] += v
769 self.diel_data = dk
771 def write_fastcap(self, output_dir_path: str, prefix: str) -> List[str]:
772 """Return the paths written, with the primary list file first."""
773 max_filename_length: Optional[int] = None
774 try:
775 max_filename_length = os.pathconf(output_dir_path, 'PC_NAME_MAX')
776 except AttributeError:
777 pass # NOTE: windows does not support the os.pathconf attribute
779 lst_fn = os.path.join(output_dir_path, f"{prefix}.lst")
780 written: List[str] = []
781 file_num = 0
782 lst_file: List[str] = [f"* k_void={'%.12g' % self.k_void}"]
784 for k, data in self.diel_data.items():
785 if len(data) == 0:
786 continue
788 file_num += 1
790 k_outside = self.materials[k.outside] if k.outside else self.k_void
791 k_inside = self.materials[k.inside] if k.inside else self.k_void
793 # lst_file.append(f"* Dielectric interface: outside={outside}, inside={inside}")
795 fn = f"{prefix}{file_num}_outside={k.outside or '(void)'}_inside={k.inside or '(void)'}.geo"
796 output_path = os.path.join(output_dir_path, fn)
797 self._write_fastercap_geo(output_path=output_path,
798 data=data,
799 cond_name=None,
800 cond_number=file_num,
801 rename_conductor=False)
802 written.append(output_path)
804 # NOTE: for now, we compute the reference points for each triangle
805 # This is a FasterCap feature, reference point in the *.geo file (end of each T line)
806 rp_s = "0 0 0"
807 lst_file.append(f"D {fn} {'%.12g' % k_outside} {'%.12g' % k_inside} 0 0 0 {rp_s}")
809 #
810 # Feedback from FastFieldSolvers:
811 #
812 # - using the '+' trailing statements (conductor collation),
813 # only the same conductor should be collated
814 #
815 # - renaming different conductor numbers ('N' rule line) is not allowed (currently a bug)
816 # - Example: 1->VDD (1.geo) and 2->VDD (2.geo) is not possible
817 # - Both conductor *.geo files should have the same number
818 # - only the last conductor *.geo file should contain the 'N' rule
819 #
820 # - reference points
821 #
822 cond_data_grouped_by_net = defaultdict(list)
823 for k, data in self.cond_data.items():
824 if len(data) == 0:
825 continue
826 cond_data_grouped_by_net[k.net_name].append((k.outside, data))
828 cond_num = file_num
830 for nn, cond_list in cond_data_grouped_by_net.items():
831 cond_num += 1
832 last_cond_index = len(cond_list) - 1
833 for idx, (outside, data) in enumerate(cond_list):
834 file_num += 1
835 k_outside = self.materials[outside] if outside else self.k_void
837 outside = outside or '(void)'
838 # lst_file.append(f"* Conductor interface: outside={outside}, net={nn}")
839 fn = f"{prefix}{file_num}_outside={outside}_net={nn}.geo"
840 if max_filename_length is not None and len(fn) > max_filename_length:
841 warning(f"Unusual long net name detected: {nn}")
842 d = hashlib.md5(nn.encode('utf-8')).digest()
843 h = base64.urlsafe_b64encode(d).decode('utf-8').rstrip('=')
844 remaining_len = len(f"{prefix}_{file_num}_outside={outside}_net=.geo")
845 short_nn = nn[0: (max_filename_length - remaining_len - len(h) - 1)] + f"_{h}"
846 fn = f"{prefix}{file_num}_outside={outside}_net={short_nn}.geo"
847 output_path = os.path.join(output_dir_path, fn)
848 self._write_fastercap_geo(output_path=output_path,
849 data=data,
850 cond_number=cond_num,
851 cond_name=nn,
852 rename_conductor=(idx == last_cond_index))
853 written.append(output_path)
854 collation_operator = '' if idx == last_cond_index else ' +'
855 lst_file.append(f"C {fn} {'%.12g' % k_outside} 0 0 0{collation_operator}")
857 subproc(lst_fn)
858 with open(lst_fn, 'w', encoding='utf-8') as f:
859 f.write('\n'.join(lst_file))
860 f.write('\n')
862 return [lst_fn, *written]
864 @staticmethod
865 def _write_fastercap_geo(output_path: str,
866 data: List[Triangle],
867 cond_number: int,
868 cond_name: Optional[str],
869 rename_conductor: bool):
870 subproc(output_path)
871 with open(output_path, 'w', encoding='utf-8') as f:
872 f.write(f"0 GEO File\n")
873 for t in data:
874 f.write(f"T {cond_number}")
875 f.write(' ' + t.to_fastcap())
877 # compute a reference point "outside"
878 rp = t.outside_reference_point()
879 rp_s = rp.to_fastcap()
881 f.write(f" {rp_s}\n")
882 if cond_name and rename_conductor:
883 f.write(f"N {cond_number} {cond_name}\n")
885 def check(self):
886 info("Checking …")
887 errors = 0
889 for mn in self.materials.keys():
890 tris = self._collect_diel_tris(mn)
891 info(f"Material {mn} -> {len(tris)} triangles")
892 errors += self._check_tris(f"Material '{mn}'", tris)
894 # NOTE: the surface of a net is closed, not the ones of its parts
895 for group_name in dict.fromkeys(self.group_by_conductor.get(nn, nn) for nn in self.net_names):
896 tris = [t for nn in self.net_names if self.group_by_conductor.get(nn, nn) == group_name
897 for t in self._collect_cond_tris(nn)]
898 info(f"Net '{group_name}' -> {len(tris)} triangles")
899 errors += self._check_tris(f"Net '{group_name}'", tris)
901 if errors == 0:
902 info(" No errors found")
903 else:
904 info(f" {errors} error{'s' if errors >= 2 else ''} found")
906 def _check_tris(self, msg: str, triangles: List[Triangle]) -> int:
907 errors = 0
909 edge_set: Set[Edge] = set()
910 edges = self._normed_edges(triangles)
912 for e in edges:
913 if e in edge_set:
914 error(f"{msg}: duplicate edge {self._edge2s(e)}")
915 errors += 1
916 else:
917 edge_set.add(e)
919 self._split_edges(edge_set)
921 for e in edge_set:
922 if e.reversed() not in edge_set:
923 error(f"{msg}: edge {self._edge2s(e)} not connected with reverse edge (open surface)")
924 errors += 1
926 return errors
928 def _normed_edges(self, triangles: List[Triangle]) -> List[Edge]:
929 edges = []
931 def normed_dbu(p: Point):
932 return Point(*tuple(map(lambda c: math.floor(c / self.dbu + 0.5),
933 (p.x, p.y, p.z))))
935 for t in triangles:
936 for i in range(0, 3):
937 p1 = normed_dbu(t[i])
938 p2 = normed_dbu(t[(i + 1) % 3])
939 edges.append(Edge(p1, p2))
941 return edges
943 def _point2s(self, p: Point) -> str:
944 return f"(%.12g, %.12g, %.12g)" % (p.x * self.dbu, p.y * self.dbu, p.z * self.dbu)
946 def _edge2s(self, e: Edge) -> str:
947 return f"{self._point2s(e.p0)}-{self._point2s(e.p1)}"
949 @staticmethod
950 def _is_antiparallel(a: Point,
951 b: Point) -> bool:
952 vp = vector_product(a, b)
953 if abs(vp.sq_length()) > 0.5: # we got normalized!
954 return False
956 sp = dot_product(a, b)
957 return sp < 0
959 def _split_edges(self, edges: set[Edge]):
960 edges_by_p2: DefaultDict[Point, List[Edge]] = defaultdict(list)
961 edges_by_p1: DefaultDict[Point, List[Edge]] = defaultdict(list)
962 for e in edges:
963 edges_by_p2[e.p1].append(e)
964 edges_by_p1[e.p0].append(e)
966 i = 0
967 while True:
968 i += 1
969 subst: DefaultDict[Edge, List[Edge]] = defaultdict(list)
971 for e in edges:
972 ee = edges_by_p2.get(e.p0, [])
973 for eee in ee:
974 ve = e.vector_of_edge()
975 veee = eee.vector_of_edge()
976 if self._is_antiparallel(ve, veee) and \
977 (veee.sq_length() < ve.sq_length() - 0.5):
978 # There is a shorter edge antiparallel ->
979 # this means we need to insert a split point into e
980 subst[e] += [Edge(e.p0, eee.p0), Edge(eee.p0, e.p1)]
982 for e in edges:
983 ee = edges_by_p1.get(e.p1, [])
984 for eee in ee:
985 ve = e.vector_of_edge()
986 veee = eee.vector_of_edge()
987 if self._is_antiparallel(ve, veee) and \
988 (veee.sq_length() < ve.sq_length() - 0.5):
989 # There is a shorter edge antiparallel ->
990 # this means we need to insert a split point into e
991 subst[e] += [Edge(e.p0, eee.p1), Edge(eee.p1, e.p1)]
993 if len(subst) == 0:
994 break
996 for e, replacement in subst.items():
997 edges_by_p1[e.p0].remove(e)
998 edges_by_p2[e.p1].remove(e)
999 edges.remove(e)
1000 for r in replacement:
1001 edges.add(r)
1002 edges_by_p1[r.p0].append(r)
1003 edges_by_p2[r.p1].append(r)
1005 def dump_stl(self, output_dir_path: str, prefix: str) -> List[str]:
1006 """Return STL paths written; empty solids produce no file."""
1007 written: List[str] = []
1008 for mn in self.materials.keys():
1009 tris = self._collect_diel_tris(mn)
1010 output_path = os.path.join(output_dir_path, f"{prefix}diel_{mn}.stl")
1011 if self._write_as_stl(output_path, tris):
1012 written.append(output_path)
1014 for nn in self.net_names:
1015 tris = self._collect_cond_tris(nn)
1016 output_path = os.path.join(output_dir_path, f"{prefix}cond_{nn}.stl")
1017 if self._write_as_stl(output_path, tris):
1018 written.append(output_path)
1020 return written
1022 @staticmethod
1023 def _write_as_stl(file_name: str,
1024 tris: List[Triangle]) -> bool:
1025 if len(tris) == 0:
1026 return False
1028 subproc(file_name)
1029 with open(file_name, 'w', encoding='utf-8') as f:
1030 f.write("solid stl\n")
1031 for t in tris:
1032 f.write(" facet normal 0 0 0\n")
1033 f.write(" outer loop\n")
1034 t = t.reversed()
1035 for p in (t.p0, t.p1, t.p2):
1036 f.write(f" vertex {p.to_fastcap()}\n")
1037 f.write(" endloop\n")
1038 f.write(" endfacet\n")
1039 f.write("endsolid stl\n")
1040 return True
1042 @staticmethod
1043 def _merge_events(pyra: List[Optional[kdb.Region]],
1044 lin: Optional[kdb.Region],
1045 lout: Optional[kdb.Region]) -> Tuple[kdb.Region, kdb.Region, kdb.Region]:
1046 lin = lin.dup() if lin else kdb.Region()
1047 lout = lout.dup() if lout else kdb.Region()
1048 past = pyra[0].dup() if len(pyra) >= 1 else kdb.Region()
1050 for i in range(0, len(pyra)):
1051 ii = len(pyra) - i
1052 added: kdb.Region = lin & pyra[ii - 1]
1053 if not added.is_empty():
1054 if ii >= len(pyra):
1055 pyra.append(kdb.Region())
1056 assert len(pyra) == ii + 1
1057 pyra[ii] += added
1058 lin -= added
1060 if len(pyra) == 0:
1061 pyra.append(kdb.Region())
1062 pyra[0] += lin
1064 for i in range(0, len(pyra)):
1065 ii = len(pyra) - i
1066 removed: kdb.Region = lout & pyra[ii - 1]
1067 if not removed.is_empty():
1068 pyra[ii - 1] -= removed
1069 lout -= removed
1071 # compute merged events
1072 lin = pyra[0] - past
1073 lout = past - pyra[0]
1074 return lin, lout, pyra[0]
1076 def _collect_diel_tris(self, material_name: str) -> List[Triangle]:
1077 tris = []
1079 for k, v in self.diel_data.items():
1080 if material_name == k.outside:
1081 tris += v
1082 elif material_name == k.inside:
1083 tris += [t.reversed() for t in v]
1085 for k, v in self.cond_data.items():
1086 if material_name == k.outside:
1087 tris += v
1089 return tris
1091 def _collect_cond_tris(self, net_name: str) -> List[Triangle]:
1092 tris = []
1093 for k, v in self.cond_data.items():
1094 if k.net_name == net_name:
1095 tris += [t.reversed() for t in v]
1096 return tris