Coverage for klayout_pex/fastercap/fastercap_model_generator.py: 96%
662 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-17 17:24 +0000
1#
2# --------------------------------------------------------------------------------
3# SPDX-FileCopyrightText: 2024 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/martinjankoehler/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] = []
116 # layer, zstart, zstop
117 self.clayers: Dict[str, List[Tuple[kdb.Region, float, float]]] = {}
118 self.dlayers: Dict[str, List[Tuple[kdb.Region, float, float]]] = {}
120 info(f"DBU: {'%.12g' % self.dbu}")
121 info(f"Delaunay b: {'%.12g' % self.delaunay_b}")
122 info(f"Delaunay area_max: {'%.12g' % self.delaunay_amax}")
124 def add_material(self, name: str, k: float):
125 self.materials[name] = k
127 def add_dielectric(self,
128 material_name: str,
129 layer: kdb.Region,
130 z: float,
131 height: float):
132 if hasattr(layer, 'data'):
133 layer = layer.data
134 self._add_layer(name=material_name, layer=layer, is_dielectric=True, z=z, height=height)
136 def add_conductor(self,
137 net_name: str,
138 layer: kdb.Region,
139 z: float,
140 height: float):
141 if hasattr(layer, 'data'):
142 layer = layer.data
143 self._add_layer(name=net_name, layer=layer, is_dielectric=False, z=z, height=height)
145 def _norm2z(self, z: float) -> float:
146 return z * self.dbu
148 def _z2norm(self, z: float) -> float:
149 return math.floor(z / self.dbu + 1e-6)
151 def _add_layer(self,
152 name: str,
153 layer: kdb.Region,
154 z: float,
155 height: float,
156 is_dielectric: bool):
157 if is_dielectric and name not in self.materials:
158 raise ValueError(f"Unknown material {name} - did you use 'add_material'?")
160 zstart: float = z
161 zstop: float = zstart + height
163 if is_dielectric:
164 if name not in self.dlayers:
165 self.dlayers[name] = []
166 self.dlayers[name].append((layer, self._z2norm(zstart), self._z2norm(zstop)))
167 else:
168 if name not in self.clayers:
169 self.clayers[name] = []
170 self.clayers[name].append((layer, self._z2norm(zstart), self._z2norm(zstop)))
172 def generate(self) -> Optional[FasterCapModelGenerator]:
173 z: List[float] = []
174 for ll in (self.dlayers, self.clayers):
175 for k, v in ll.items():
176 for l in v:
177 z.extend((l[1], l[2]))
178 z = sorted([*{*z}]) # sort & uniq
179 if len(z) == 0:
180 return None
182 gen = FasterCapModelGenerator(dbu=self.dbu,
183 k_void=self.k_void,
184 delaunay_amax=self.delaunay_amax,
185 delaunay_b=self.delaunay_b,
186 materials=self.materials,
187 net_names=list(self.clayers.keys()))
188 for zcurr in z:
189 gen.next_z(self._norm2z(zcurr))
191 for nn, v in self.clayers.items():
192 for l in v:
193 if l[1] <= zcurr < l[2]:
194 gen.add_in(name=f"+{nn}", layer=l[0])
195 if l[1] < zcurr <= l[2]:
196 gen.add_out(name=f"+{nn}", layer=l[0])
197 for mn, v in self.dlayers.items():
198 for l in v:
199 if l[1] <= zcurr < l[2]:
200 gen.add_in(name=f"-{mn}", layer=l[0])
201 if l[1] < zcurr <= l[2]:
202 gen.add_out(name=f"-{mn}", layer=l[0])
204 gen.finish_z()
206 gen.finalize()
207 return gen
210@dataclass(frozen=True)
211class HDielKey:
212 outside: Optional[str]
213 inside: Optional[str]
215 def __str__(self) -> str:
216 return f"{self.outside or 'void'} <-> {self.inside or 'void'}"
218 @property
219 def topic(self) -> str:
220 return 'dielectric'
222 def reversed(self) -> HDielKey:
223 return HDielKey(self.inside, self.outside)
226@dataclass(frozen=True)
227class HCondKey:
228 net_name: str
229 outside: Optional[str]
231 def __str__(self) -> str:
232 return f"{self.outside or 'void'} <-> {self.net_name}"
234 @property
235 def topic(self) -> str:
236 return 'conductor'
239@dataclass(frozen=True)
240class VKey:
241 kk: HDielKey | HCondKey
242 p0: kdb.DPoint
243 de: kdb.DVector
246@dataclass(frozen=True)
247class Point:
248 x: float
249 y: float
250 z: float
252 def __sub__(self, other: Point) -> Point:
253 return Point(self.x - other.x, self.y - other.y, self.z - other.z)
255 def sq_length(self) -> float:
256 return self.x**2 + self.y**2 + self.z**2
258 def to_fastcap(self) -> str:
259 return '%.12g %.12g %.12g' % (self.x, self.y, self.z)
262def vector_product(a: Point, b: Point) -> Point:
263 vp = Point(
264 a.y * b.z - a.z * b.y,
265 a.z * b.x - a.x * b.z,
266 a.x * b.y - a.y * b.x
267 )
268 return vp
271def dot_product(a: Point, b: Point) -> float:
272 dp = a.x * b.x + a.y * b.y + a.z * b.z
273 return dp
276@dataclass(frozen=True)
277class Triangle:
278 p0: Point
279 p1: Point
280 p2: Point
282 def reversed(self) -> Triangle:
283 return Triangle(self.p2, self.p1, self.p0)
285 def outside_reference_point(self) -> Point:
286 v1 = self.p1 - self.p0
287 v2 = self.p2 - self.p0
288 vp = Point(v1.y * v2.z - v1.z * v2.y,
289 -v1.x * v2.z + v1.z * v2.x,
290 v1.x * v2.y - v1.y * v2.x)
291 vp_abs = math.sqrt(vp.x ** 2 + vp.y ** 2 + vp.z ** 2)
292 rp = Point(self.p0.x + vp.x / vp_abs,
293 self.p0.y + vp.y / vp_abs,
294 self.p0.z + vp.z / vp_abs)
295 return rp
297 def to_fastcap(self) -> str:
298 return ' '.join([p.to_fastcap() for p in (self.p0, self.p1, self.p2)])
300 def __len__(self):
301 return 3
303 def __getitem__(self, i) -> Point:
304 match i:
305 case 0: return self.p0
306 case 1: return self.p1
307 case 2: return self.p2
308 case _: raise IndexError("list index out of range")
311@dataclass(frozen=True)
312class Edge:
313 p0: Point
314 p1: Point
316 def vector_of_edge(self) -> Point:
317 return Point(
318 self.p1.x - self.p0.x,
319 self.p1.y - self.p0.y,
320 self.p1.z - self.p0.z
321 )
323 def reversed(self) -> Edge:
324 return Edge(self.p1, self.p0)
327@dataclass
328class FasterCapModelGenerator:
329 dbu: float
330 """Database unit"""
332 k_void: float
333 """Default dielectric of 'empty space'"""
335 delaunay_amax: float
336 """Maximum area parameter for the Delaunay triangulation"""
338 delaunay_b: float
339 """
340 The delaunay_b parameter for the Delaunay triangulation
341 corresponds to the minimum angle
342 and should be <=1 (b=sin(min_angle)*2).
343 I.e. b=1 -> min_angle=30 deg, b=0.5 -> min_angle~14.5 deg.
344 """
346 materials: Dict[str, float]
347 """Maps material name to dielectric k"""
349 net_names: List[str]
351 def __init__(self,
352 dbu: float,
353 k_void: float,
354 delaunay_amax: float,
355 delaunay_b: float,
356 materials: Dict[str, float],
357 net_names: List[str]):
358 self.k_void = k_void
359 self.delaunay_amax = delaunay_amax
360 self.delaunay_b = delaunay_b
361 self.dbu = dbu
362 self.materials = materials
363 self.net_names = net_names
365 self.z: Optional[float] = None
366 self.zz: Optional[float] = None
367 self.layers_in: Dict[str, kdb.Region] = {}
368 self.layers_out: Dict[str, kdb.Region] = {}
369 self.state: Dict[str, kdb.Region] = {}
370 self.current: Dict[str, List[kdb.Region]] = {}
371 self.diel_data: Dict[HDielKey, List[Triangle]] = {}
372 self.diel_vdata: Dict[VKey, kdb.Region] = {}
373 self.cond_data: Dict[HCondKey, List[Triangle]] = {}
374 self.cond_vdata: Dict[VKey, kdb.Region] = {}
376 def reset(self):
377 self.layers_in = {}
378 self.layers_out = {}
380 def add_in(self, name: str, layer: kdb.Region):
381 debug(f"add_in: {name} -> {layer}")
382 if name not in self.layers_in:
383 # NOTE: Calling kdb.Region([]) with the empty list enforces, that a flat layer is created.
384 # The following "+=" packs the polynomials into this flat layer and
385 # does not copy the hierarchical layer from the LayoutToNetlist database.
386 self.layers_in[name] = kdb.Region([])
387 self.layers_in[name] += layer
389 def add_out(self, name: str, layer: kdb.Region):
390 debug(f"add_out: {name} -> {layer}")
391 if name not in self.layers_out:
392 # NOTE: Calling kdb.Region([]) with the empty list enforces, that a flat layer is created.
393 # The following "+=" packs the polynomials into this flat layer and
394 # does not copy the hierarchical layer from the LayoutToNetlist database.
395 self.layers_out[name] = kdb.Region([])
396 self.layers_out[name] += layer
398 def finish_z(self):
399 debug(f"Finishing layer z={self.z}")
401 din: Dict[str, kdb.Region] = {}
402 dout: Dict[str, kdb.Region] = {}
403 all_in = kdb.Region()
404 all_out = kdb.Region()
405 all = kdb.Region()
406 all_cin: Optional[kdb.Region] = None
407 all_cout: Optional[kdb.Region] = None
409 for names, prefix in ((self.net_names, '+'), (self.materials.keys(), '-')):
410 for nn in names:
411 mk = prefix + nn
413 # compute merged events
414 if mk not in self.current:
415 self.current[mk] = []
416 current_before = self.current[mk][0].dup() if len(self.current[mk]) >= 1 else kdb.Region()
417 lin, lout, current = self._merge_events(pyra=self.current[mk],
418 lin=self.layers_in.get(mk, None),
419 lout=self.layers_out.get(mk, None))
420 debug(f"Merged events & status for {mk}:")
421 debug(f" in = {lin}")
422 debug(f" out = {lout}")
423 debug(f" state = {current}")
425 if mk not in self.state:
426 self.state[mk] = kdb.Region()
428 # legalize in and out events
429 lin_org = lin.dup()
430 lout_org = lout.dup()
431 lout &= self.state[mk]
432 lin -= all
433 lout += current & all_in
434 lin += current_before & all_out
435 lin -= lout_org
436 lout -= lin_org
438 # tracks the legalized horizontal cuts
439 self.state[mk] += lin
440 self.state[mk] -= lout
442 din[mk] = lin
443 dout[mk] = lout
445 debug(f"Legalized events & status for '{mk}':")
446 debug(f" in = {din[mk]}")
447 debug(f" out = {dout[mk]}")
448 debug(f" state = {self.state[mk]}")
450 all_in += lin
451 all_out += lout
452 all += self.state[mk]
454 if prefix == '+':
455 all_cin = all_in.dup()
456 all_cout = all_out.dup()
458 debug(f"All conductor region in: {all_cin}")
459 debug(f"All conductor region out: {all_cout}")
461 # check whether states are separated
462 a = reduce(lambda x, y: x+y, self.state.values())
463 for k, s in self.state.items():
464 r: kdb.Region = s - a
465 if not r.is_empty():
466 error(f"State region of {k} ({s}) is not contained entirely "
467 f"in remaining all state region ({a}) - this means there is an overlap")
468 a -= s
470 # Now we have legalized the in and out events
471 for mni in self.materials.keys():
472 lin = din.get(f"-{mni}", None)
473 if lin:
474 lin = lin.dup()
475 lin -= all_cout # handled with the conductor
476 for mno in self.materials.keys():
477 lout = dout.get(f"-{mno}", None)
478 if lout:
479 d: kdb.Region = lout & lin
480 if not d.is_empty():
481 self.generate_hdiel(below=mno, above=mni, layer=d)
482 lin -= lout
483 if not lin.is_empty():
484 self.generate_hdiel(below=None, above=mni, layer=lin)
486 for mno in self.materials.keys():
487 lout = dout.get(f"-{mno}", None)
488 if lout:
489 lout = lout.dup()
490 lout -= all_cin # handled with the conductor
491 for mni in self.materials.keys():
492 lin = din.get(f"-{mni}", None)
493 if lin:
494 lout -= lin
495 if not lout.is_empty():
496 self.generate_hdiel(below=mno, above=None, layer=lout)
498 for nn in self.net_names:
499 lin = din.get(f"+{nn}", None)
500 if lin:
501 lin = lin.dup()
502 for mno in self.materials.keys():
503 lout = dout.get(f"-{mno}", None)
504 if lout:
505 d = lout & lin
506 if not d.is_empty():
507 self.generate_hcond_in(net_name=nn, below=mno, layer=d)
508 lin -= lout
509 if not lin.is_empty():
510 self.generate_hcond_in(net_name=nn, below=None, layer=lin)
512 for nn in self.net_names:
513 lout = dout.get(f"+{nn}", None)
514 if lout:
515 lout = lout.dup()
516 lout -= all_cin # handled with the conductor
517 for mni in self.materials.keys():
518 lin = din.get(f"-{mni}", None)
519 if lin:
520 d = lout & lin
521 if not d.is_empty():
522 self.generate_hcond_out(net_name=nn, above=mni, layer=d)
523 lout -= lin
524 if not lout.is_empty():
525 self.generate_hcond_out(net_name=nn, above=None, layer=lout)
527 def next_z(self, z: float):
528 debug(f"Next layer {z}")
530 self.reset()
532 if self.z is None:
533 self.z = z
534 return
536 self.zz = z
538 all_cond = kdb.Region()
539 for nn in self.net_names:
540 mk = f"+{nn}"
541 if mk in self.state:
542 all_cond += self.state[mk]
543 all_cond = all_cond.edges()
545 for i, mni in enumerate(self.materials):
546 linside = self.state.get(f"-{mni}", None)
547 if linside:
548 linside = linside.edges()
549 linside -= all_cond # handled with the conductor
550 for o, mno in enumerate(self.materials):
551 if i != o:
552 loutside = self.state.get(f"-{mno}", None)
553 if loutside:
554 loutside = loutside.edges()
555 if o > i:
556 d = loutside & linside
557 for e in d:
558 # NOTE: we need to swap points as we started from "outside"
559 self.generate_vdiel(outside=mno, inside=mni, edge=e.swapped_points())
560 linside -= loutside
562 for e in linside:
563 self.generate_vdiel(outside=None, inside=mni, edge=e)
565 for nn in self.net_names:
566 mk = f"+{nn}"
567 linside = self.state.get(mk, None)
568 if linside:
569 linside = linside.edges()
570 for mno in self.materials:
571 loutside = self.state.get(f"-{mno}", None)
572 if loutside:
573 loutside = loutside.edges()
574 d = loutside & linside
575 for e in d:
576 # NOTE: we need to swap points as we started from "outside"
577 self.generate_vcond(net_name=nn, outside=mno, edge=e.swapped_points())
578 linside -= loutside
579 for e in linside:
580 self.generate_vcond(net_name=nn, outside=None, edge=e)
582 self.z = z
584 def generate_hdiel(self,
585 below: Optional[str],
586 above: Optional[str],
587 layer: kdb.Region):
588 k = HDielKey(below, above)
589 debug(f"Generating horizontal dielectric surface {k} as {layer}")
590 if k not in self.diel_data:
591 self.diel_data[k] = []
592 data = self.diel_data[k]
594 for t in layer.delaunay(self.delaunay_amax / self.dbu ** 2, self.delaunay_b):
595 # NOTE: normal is facing downwards (to "below")
596 pl = list(map(lambda pt: Point(pt.x * self.dbu, pt.y * self.dbu, self.z),
597 t.each_point_hull()))
598 tri = Triangle(*pl)
599 data.append(tri)
600 debug(f" {tri}")
602 def generate_v_surface(self,
603 kk: HDielKey | HCondKey,
604 edge: kdb.Edge) -> Tuple[VKey, kdb.Box]:
605 debug(f"Generating vertical {kk.topic} surface {kk} with edge {edge}")
607 el = math.sqrt(edge.sq_length())
608 de = kdb.DVector(edge.d().x / el, edge.d().y / el)
609 ne = kdb.DVector(edge.d().y / el, -edge.d().x / el)
610 p0 = ne * ne.sprod(kdb.DPoint(edge.p1) - kdb.DPoint()) + kdb.DPoint()
611 x1 = (edge.p1 - p0).sprod(de)
612 x2 = (edge.p2 - p0).sprod(de)
614 key = VKey(kk, p0, de)
615 surface = kdb.Box(x1,
616 math.floor(self.z / self.dbu + 0.5),
617 x2,
618 math.floor(self.zz / self.dbu + 0.5))
619 return key, surface
621 def generate_vdiel(self,
622 outside: Optional[str],
623 inside: Optional[str],
624 edge: kdb.Edge):
625 if edge.is_degenerate():
626 return
628 key, surface = self.generate_v_surface(HDielKey(outside, inside), edge)
629 if key not in self.diel_vdata:
630 self.diel_vdata[key] = kdb.Region()
632 self.diel_vdata[key].insert(surface)
634 def generate_hcond_in(self,
635 net_name: str,
636 below: Optional[str],
637 layer: kdb.Region):
638 k = HCondKey(net_name, below)
639 debug(f"Generating horizontal bottom conductor surface {k} as {layer}")
641 if k not in self.cond_data:
642 self.cond_data[k] = []
643 data = self.cond_data[k]
645 for t in layer.delaunay(self.delaunay_amax / self.dbu ** 2, self.delaunay_b):
646 # NOTE: normal is facing downwards (to "below")
647 pl = list(map(lambda pt: Point(pt.x * self.dbu, pt.y * self.dbu, self.z),
648 t.each_point_hull()))
649 tri = Triangle(*pl)
650 data.append(tri)
651 debug(f" {tri}")
653 def generate_hcond_out(self,
654 net_name: str,
655 above: Optional[str],
656 layer: kdb.Region):
657 k = HCondKey(net_name, above)
658 debug(f"Generating horizontal top conductor surface {k} as {layer}")
660 if k not in self.cond_data:
661 self.cond_data[k] = []
662 data = self.cond_data[k]
664 for t in layer.delaunay(self.delaunay_amax / self.dbu ** 2, self.delaunay_b):
665 # NOTE: normal is facing downwards (into conductor)
666 pl = list(map(lambda pt: Point(pt.x * self.dbu, pt.y * self.dbu, self.z),
667 t.each_point_hull()))
668 tri = Triangle(*pl)
669 # now it is facing outside (to "above")
670 tri = tri.reversed()
671 data.append(tri)
672 debug(f" {tri}")
674 def generate_vcond(self,
675 net_name: str,
676 outside: Optional[str],
677 edge: kdb.Edge):
678 if edge.is_degenerate():
679 return
681 key, surface = self.generate_v_surface(HCondKey(net_name, outside), edge)
682 if key not in self.cond_vdata:
683 self.cond_vdata[key] = kdb.Region()
685 self.cond_vdata[key].insert(surface)
687 def triangulate(self, p0: kdb.DPoint, de: kdb.DVector, region: kdb.Region, data: List[Triangle]):
688 def convert_point(pt: kdb.Point) -> Point:
689 pxy = (p0 + de * pt.x) * self.dbu
690 pz = pt.y * self.dbu
691 return Point(pxy.x, pxy.y, pz)
693 for t in region.delaunay(self.delaunay_amax / self.dbu ** 2, self.delaunay_b):
694 # NOTE: normal is facing outwards (to "left")
695 pl = list(map(convert_point, t.each_point_hull()))
696 tri = Triangle(*pl)
697 # now it is facing outside (to "above")
698 data.append(tri)
699 debug(f" {tri}")
701 def finalize(self):
702 for k, r in self.diel_vdata.items():
703 debug(f"Finishing vertical dielectric plane {k.kk} at {k.p0}/{k.de}")
705 if k.kk not in self.diel_data:
706 self.diel_data[k.kk] = []
707 data = self.diel_data[k.kk]
709 self.triangulate(p0=k.p0, de=k.de, region=r, data=data)
711 for k, r in self.cond_vdata.items():
712 debug(f"Finishing vertical conductor plane {k.kk} at {k.p0} / {k.de}")
714 if k.kk not in self.cond_data:
715 self.cond_data[k.kk] = []
716 data = self.cond_data[k.kk]
718 self.triangulate(p0=k.p0, de=k.de, region=r, data=data)
720 dk: Dict[HDielKey, List[Triangle]] = {}
722 for k in self.diel_data.keys():
723 kk = k.reversed()
724 if kk not in dk:
725 dk[k] = []
726 else:
727 debug(f"Combining dielectric surfaces {kk} with reverse")
729 for k, v in self.diel_data.items():
730 kk = k.reversed()
731 if kk in dk:
732 dk[kk] += list(map(lambda t: t.reversed(), v))
733 else:
734 dk[k] += v
736 self.diel_data = dk
738 def write_fastcap(self, output_dir_path: str, prefix: str) -> str:
739 max_filename_length: Optional[int] = None
740 try:
741 max_filename_length = os.pathconf(output_dir_path, 'PC_NAME_MAX')
742 except AttributeError:
743 pass # NOTE: windows does not support the os.pathconf attribute
745 lst_fn = os.path.join(output_dir_path, f"{prefix}.lst")
746 file_num = 0
747 lst_file: List[str] = [f"* k_void={'%.12g' % self.k_void}"]
749 for k, data in self.diel_data.items():
750 if len(data) == 0:
751 continue
753 file_num += 1
755 k_outside = self.materials[k.outside] if k.outside else self.k_void
756 k_inside = self.materials[k.inside] if k.inside else self.k_void
758 # lst_file.append(f"* Dielectric interface: outside={outside}, inside={inside}")
760 fn = f"{prefix}{file_num}_outside={k.outside or '(void)'}_inside={k.inside or '(void)'}.geo"
761 output_path = os.path.join(output_dir_path, fn)
762 self._write_fastercap_geo(output_path=output_path,
763 data=data,
764 cond_name=None,
765 cond_number=file_num,
766 rename_conductor=False)
768 # NOTE: for now, we compute the reference points for each triangle
769 # This is a FasterCap feature, reference point in the *.geo file (end of each T line)
770 rp_s = "0 0 0"
771 lst_file.append(f"D {fn} {'%.12g' % k_outside} {'%.12g' % k_inside} 0 0 0 {rp_s}")
773 #
774 # Feedback from FastFieldSolvers:
775 #
776 # - using the '+' trailing statements (conductor collation),
777 # only the same conductor should be collated
778 #
779 # - renaming different conductor numbers ('N' rule line) is not allowed (currently a bug)
780 # - Example: 1->VDD (1.geo) and 2->VDD (2.geo) is not possible
781 # - Both conductor *.geo files should have the same number
782 # - only the last conductor *.geo file should contain the 'N' rule
783 #
784 # - reference points
785 #
786 cond_data_grouped_by_net = defaultdict(list)
787 for k, data in self.cond_data.items():
788 if len(data) == 0:
789 continue
790 cond_data_grouped_by_net[k.net_name].append((k.outside, data))
792 cond_num = file_num
794 for nn, cond_list in cond_data_grouped_by_net.items():
795 cond_num += 1
796 last_cond_index = len(cond_list) - 1
797 for idx, (outside, data) in enumerate(cond_list):
798 file_num += 1
799 k_outside = self.materials[outside] if outside else self.k_void
801 outside = outside or '(void)'
802 # lst_file.append(f"* Conductor interface: outside={outside}, net={nn}")
803 fn = f"{prefix}{file_num}_outside={outside}_net={nn}.geo"
804 if max_filename_length is not None and len(fn) > max_filename_length:
805 warning(f"Unusual long net name detected: {nn}")
806 d = hashlib.md5(nn.encode('utf-8')).digest()
807 h = base64.urlsafe_b64encode(d).decode('utf-8').rstrip('=')
808 remaining_len = len(f"{prefix}_{file_num}_outside={outside}_net=.geo")
809 short_nn = nn[0: (max_filename_length - remaining_len - len(h) - 1)] + f"_{h}"
810 fn = f"{prefix}{file_num}_outside={outside}_net={short_nn}.geo"
811 output_path = os.path.join(output_dir_path, fn)
812 self._write_fastercap_geo(output_path=output_path,
813 data=data,
814 cond_number=cond_num,
815 cond_name=nn,
816 rename_conductor=(idx == last_cond_index))
817 collation_operator = '' if idx == last_cond_index else ' +'
818 lst_file.append(f"C {fn} {'%.12g' % k_outside} 0 0 0{collation_operator}")
820 subproc(lst_fn)
821 with open(lst_fn, "w") as f:
822 f.write('\n'.join(lst_file))
823 f.write('\n')
825 return lst_fn
827 @staticmethod
828 def _write_fastercap_geo(output_path: str,
829 data: List[Triangle],
830 cond_number: int,
831 cond_name: Optional[str],
832 rename_conductor: bool):
833 subproc(output_path)
834 with open(output_path, "w") as f:
835 f.write(f"0 GEO File\n")
836 for t in data:
837 f.write(f"T {cond_number}")
838 f.write(' ' + t.to_fastcap())
840 # compute a reference point "outside"
841 rp = t.outside_reference_point()
842 rp_s = rp.to_fastcap()
844 f.write(f" {rp_s}\n")
845 if cond_name and rename_conductor:
846 f.write(f"N {cond_number} {cond_name}\n")
848 def check(self):
849 info("Checking …")
850 errors = 0
852 for mn in self.materials.keys():
853 tris = self._collect_diel_tris(mn)
854 info(f"Material {mn} -> {len(tris)} triangles")
855 errors += self._check_tris(f"Material '{mn}'", tris)
857 for nn in self.net_names:
858 tris = self._collect_cond_tris(nn)
859 info(f"Net '{nn}' -> {len(tris)} triangles")
860 errors += self._check_tris(f"Net '{nn}'", tris)
862 if errors == 0:
863 info(" No errors found")
864 else:
865 info(f" {errors} error{'s' if errors >= 2 else ''} found")
867 def _check_tris(self, msg: str, triangles: List[Triangle]) -> int:
868 errors = 0
870 edge_set: Set[Edge] = set()
871 edges = self._normed_edges(triangles)
873 for e in edges:
874 if e in edge_set:
875 error(f"{msg}: duplicate edge {self._edge2s(e)}")
876 errors += 1
877 else:
878 edge_set.add(e)
880 self._split_edges(edge_set)
882 for e in edge_set:
883 if e.reversed() not in edge_set:
884 error(f"{msg}: edge {self._edge2s(e)} not connected with reverse edge (open surface)")
885 errors += 1
887 return errors
889 def _normed_edges(self, triangles: List[Triangle]) -> List[Edge]:
890 edges = []
892 def normed_dbu(p: Point):
893 return Point(*tuple(map(lambda c: math.floor(c / self.dbu + 0.5),
894 (p.x, p.y, p.z))))
896 for t in triangles:
897 for i in range(0, 3):
898 p1 = normed_dbu(t[i])
899 p2 = normed_dbu(t[(i + 1) % 3])
900 edges.append(Edge(p1, p2))
902 return edges
904 def _point2s(self, p: Point) -> str:
905 return f"(%.12g, %.12g, %.12g)" % (p.x * self.dbu, p.y * self.dbu, p.z * self.dbu)
907 def _edge2s(self, e: Edge) -> str:
908 return f"{self._point2s(e.p0)}-{self._point2s(e.p1)}"
910 @staticmethod
911 def _is_antiparallel(a: Point,
912 b: Point) -> bool:
913 vp = vector_product(a, b)
914 if abs(vp.sq_length()) > 0.5: # we got normalized!
915 return False
917 sp = dot_product(a, b)
918 return sp < 0
920 def _split_edges(self, edges: set[Edge]):
921 edges_by_p2: DefaultDict[Point, List[Edge]] = defaultdict(list)
922 edges_by_p1: DefaultDict[Point, List[Edge]] = defaultdict(list)
923 for e in edges:
924 edges_by_p2[e.p1].append(e)
925 edges_by_p1[e.p0].append(e)
927 i = 0
928 while True:
929 i += 1
930 subst: DefaultDict[Edge, List[Edge]] = defaultdict(list)
932 for e in edges:
933 ee = edges_by_p2.get(e.p0, [])
934 for eee in ee:
935 ve = e.vector_of_edge()
936 veee = eee.vector_of_edge()
937 if self._is_antiparallel(ve, veee) and \
938 (veee.sq_length() < ve.sq_length() - 0.5):
939 # There is a shorter edge antiparallel ->
940 # this means we need to insert a split point into e
941 subst[e] += [Edge(e.p0, eee.p0), Edge(eee.p0, e.p1)]
943 for e in edges:
944 ee = edges_by_p1.get(e.p1, [])
945 for eee in ee:
946 ve = e.vector_of_edge()
947 veee = eee.vector_of_edge()
948 if self._is_antiparallel(ve, veee) and \
949 (veee.sq_length() < ve.sq_length() - 0.5):
950 # There is a shorter edge antiparallel ->
951 # this means we need to insert a split point into e
952 subst[e] += [Edge(e.p0, eee.p1), Edge(eee.p1, e.p1)]
954 if len(subst) == 0:
955 break
957 for e, replacement in subst.items():
958 edges_by_p1[e.p0].remove(e)
959 edges_by_p2[e.p1].remove(e)
960 edges.remove(e)
961 for r in replacement:
962 edges.add(r)
963 edges_by_p1[r.p0].append(r)
964 edges_by_p2[r.p1].append(r)
966 def dump_stl(self, output_dir_path: str, prefix: str):
967 for mn in self.materials.keys():
968 tris = self._collect_diel_tris(mn)
969 output_path = os.path.join(output_dir_path, f"{prefix}diel_{mn}.stl")
970 self._write_as_stl(output_path, tris)
972 for nn in self.net_names:
973 tris = self._collect_cond_tris(nn)
974 output_path = os.path.join(output_dir_path, f"{prefix}cond_{nn}.stl")
975 self._write_as_stl(output_path, tris)
977 @staticmethod
978 def _write_as_stl(file_name: str,
979 tris: List[Triangle]):
980 if len(tris) == 0:
981 return
983 subproc(file_name)
984 with open(file_name, "w") as f:
985 f.write("solid stl\n")
986 for t in tris:
987 f.write(" facet normal 0 0 0\n")
988 f.write(" outer loop\n")
989 t = t.reversed()
990 for p in (t.p0, t.p1, t.p2):
991 f.write(f" vertex {p.to_fastcap()}\n")
992 f.write(" endloop\n")
993 f.write(" endfacet\n")
994 f.write("endsolid stl\n")
996 @staticmethod
997 def _merge_events(pyra: List[Optional[kdb.Region]],
998 lin: Optional[kdb.Region],
999 lout: Optional[kdb.Region]) -> Tuple[kdb.Region, kdb.Region, kdb.Region]:
1000 lin = lin.dup() if lin else kdb.Region()
1001 lout = lout.dup() if lout else kdb.Region()
1002 past = pyra[0].dup() if len(pyra) >= 1 else kdb.Region()
1004 for i in range(0, len(pyra)):
1005 ii = len(pyra) - i
1006 added: kdb.Region = lin & pyra[ii - 1]
1007 if not added.is_empty():
1008 if ii >= len(pyra):
1009 pyra.append(kdb.Region())
1010 assert len(pyra) == ii + 1
1011 pyra[ii] += added
1012 lin -= added
1014 if len(pyra) == 0:
1015 pyra.append(kdb.Region())
1016 pyra[0] += lin
1018 for i in range(0, len(pyra)):
1019 ii = len(pyra) - i
1020 removed: kdb.Region = lout & pyra[ii - 1]
1021 if not removed.is_empty():
1022 pyra[ii - 1] -= removed
1023 lout -= removed
1025 # compute merged events
1026 lin = pyra[0] - past
1027 lout = past - pyra[0]
1028 return lin, lout, pyra[0]
1030 def _collect_diel_tris(self, material_name: str) -> List[Triangle]:
1031 tris = []
1033 for k, v in self.diel_data.items():
1034 if material_name == k.outside:
1035 tris += v
1036 elif material_name == k.inside:
1037 tris += [t.reversed() for t in v]
1039 for k, v in self.cond_data.items():
1040 if material_name == k.outside:
1041 tris += v
1043 return tris
1045 def _collect_cond_tris(self, net_name: str) -> List[Triangle]:
1046 tris = []
1047 for k, v in self.cond_data.items():
1048 if k.net_name == net_name:
1049 tris += [t.reversed() for t in v]
1050 return tris