grabowski
0afb55216b
- adaptive: skip layers with zero quadtree leaves in the connectivity restriction and mesh-boundary loops (IndexError on boards where a selected layer has no copper) - board_io: accept ZT_TEARDROP zones as conducting copper (KiCad types teardrop fills ZT_TEARDROP, never ZT_COPPER, so they were dropped) - geometry: copper_bbox uses the exact stroke bbox (centerline extrema + half width) instead of a 100 um chord tessellation that could undershoot arc/cap extrema past the raster guard margin - solver/adaptive: reference |J| to the conduction-equivalent thickness sigma*rho in every branch (the uniform branch used geometric t, so AC plots changed scale ~rs_ratio depending on unrelated per-cell maps) - solver/adaptive/raster: chain cells no longer show phantom sheet-face currents; store dl per chain link and overlay the true 1D density |dV|/(rho*dl) (exact at any frequency: AC scaling of link conductance and cross-section cancels) - test_quadtree: compare edge lists pair-for-pair (the independent column sort destroyed endpoint association)
454 lines
17 KiB
Python
454 lines
17 KiB
Python
"""Plain geometry data model. No kipy imports here.
|
|
|
|
Everything is int64 nanometers in KiCad board coordinates (y grows down);
|
|
z grows from the board top surface downwards through the stackup.
|
|
Problem is the complete solver input and doubles as the JSON dump schema,
|
|
so the whole pipeline downstream of board_io runs without KiCad.
|
|
|
|
Schema v2 is multi-layer: per-layer fills at stackup depths, linked by
|
|
via/through-pad barrels. v1 dumps (single layer, no vias) still load.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import json
|
|
import math
|
|
from dataclasses import dataclass, field
|
|
from pathlib import Path
|
|
|
|
import numpy as np
|
|
|
|
JSON_SCHEMA_VERSION = 5
|
|
|
|
|
|
@dataclass(frozen=True)
|
|
class Rect:
|
|
x0: int
|
|
y0: int
|
|
x1: int
|
|
y1: int
|
|
layer_name: str
|
|
|
|
@classmethod
|
|
def normalized(cls, xa: int, ya: int, xb: int, yb: int, layer_name: str) -> "Rect":
|
|
return cls(min(xa, xb), min(ya, yb), max(xa, xb), max(ya, yb), layer_name)
|
|
|
|
@property
|
|
def width(self) -> int:
|
|
return self.x1 - self.x0
|
|
|
|
@property
|
|
def height(self) -> int:
|
|
return self.y1 - self.y0
|
|
|
|
|
|
@dataclass
|
|
class Polygon:
|
|
outline: np.ndarray # (N, 2) int64 nm, open ring
|
|
holes: list[np.ndarray] = field(default_factory=list)
|
|
|
|
|
|
@dataclass
|
|
class LayerFill:
|
|
layer_name: str
|
|
thickness_nm: int
|
|
z_nm: int # copper center depth from board top
|
|
polygons: list[Polygon]
|
|
|
|
|
|
@dataclass
|
|
class SurfaceBuildup:
|
|
"""Solder (plus optional added copper) sitting on an outer copper
|
|
layer inside solder-mask openings (zones on F.Mask/B.Mask)."""
|
|
layer_name: str # copper layer it sits on
|
|
polygons: list[Polygon]
|
|
|
|
|
|
@dataclass
|
|
class TrackSeg:
|
|
"""One trace segment: straight ((2, 2) points) or arc ((3, 2)
|
|
start/mid/end points). Kept as centerline + width so the raster can
|
|
decide per run: wide traces are rasterized from their outline,
|
|
traces narrower than TRACK_1D_FACTOR grid cells become exact 1D
|
|
resistor chains along the centerline."""
|
|
layer_name: str
|
|
points: np.ndarray # (2|3, 2) int64 nm
|
|
width_nm: int
|
|
|
|
def outline(self, tol_nm: float) -> np.ndarray:
|
|
if len(self.points) == 3:
|
|
return arc_band_ring(self.points[0], self.points[1],
|
|
self.points[2], self.width_nm, tol_nm)
|
|
return capsule_ring(int(self.points[0][0]), int(self.points[0][1]),
|
|
int(self.points[1][0]), int(self.points[1][1]),
|
|
self.width_nm, tol_nm)
|
|
|
|
def centerline(self, tol_nm: float) -> np.ndarray:
|
|
"""(N, 2) float polyline along the trace center, start to end."""
|
|
if len(self.points) == 3:
|
|
pts = arc_points(self.points[0], self.points[1], self.points[2],
|
|
tol_nm)
|
|
return np.vstack([pts, self.points[2][None, :]]).astype(float)
|
|
return self.points.astype(float)
|
|
|
|
|
|
@dataclass
|
|
class Electrode:
|
|
"""One PART of a current-injection terminal: a drawn rectangle or a
|
|
selected pad. A terminal (V+ or V-) is a LIST of parts, all merged
|
|
into one equipotential contact (externally bonded). `polygons`
|
|
(board nm) is the exact copper shape when known (pads); None means
|
|
the rectangle itself is the shape. `contact` = 'all' or a layer
|
|
name: which included layers this part touches."""
|
|
rect: Rect # bounding box (labels/summary)
|
|
contact: str = "all"
|
|
polygons: list[Polygon] | None = None
|
|
label: str = "rect"
|
|
|
|
|
|
@dataclass
|
|
class ViaLink:
|
|
"""A conductive barrel (via or plated through-hole pad) linking copper
|
|
layers whose z lies within [z_top_nm, z_bot_nm]."""
|
|
x: int
|
|
y: int
|
|
drill_nm: int
|
|
z_top_nm: int
|
|
z_bot_nm: int
|
|
kind: str = "via" # "via" | "pad"
|
|
pad_nm: int = 0 # pad/annular diameter; 0 = unknown
|
|
|
|
def spans(self, z_nm: int) -> bool:
|
|
return self.z_top_nm - 1 <= z_nm <= self.z_bot_nm + 1
|
|
|
|
def barrel_resistance(self, length_nm: int, rho_ohm_m: float,
|
|
plating_nm: int) -> float:
|
|
"""Barrel segment resistance over length_nm: thin-wall annulus of
|
|
plating around the drill."""
|
|
area_m2 = math.pi * (self.drill_nm * 1e-9) * (plating_nm * 1e-9)
|
|
return rho_ohm_m * (length_nm * 1e-9) / area_m2
|
|
|
|
|
|
@dataclass
|
|
class Problem:
|
|
board_path: str
|
|
net_name: str
|
|
rho_ohm_m: float
|
|
plating_nm: int
|
|
layers: list[LayerFill] # sorted by z_nm (top first)
|
|
vias: list[ViaLink]
|
|
electrodes1: list[Electrode] # V+ terminal parts (merged)
|
|
electrodes2: list[Electrode] # V- terminal parts (merged)
|
|
thickness_source: str = "stackup"
|
|
buildups: list[SurfaceBuildup] = field(default_factory=list)
|
|
solder_thickness_nm: int = 50_000
|
|
solder_rho_ohm_m: float = 1.32e-7
|
|
extra_cu_nm: int = 0
|
|
tracks: list[TrackSeg] = field(default_factory=list)
|
|
vias_capped: bool = True # filled+capped vias: thin cap
|
|
cap_plating_nm: int = 15_000 # over outer-layer mouths;
|
|
# False = open mouths
|
|
|
|
@property
|
|
def layer_names(self) -> list[str]:
|
|
return [l.layer_name for l in self.layers]
|
|
|
|
def sigma_s(self, layer_index: int) -> float:
|
|
"""Sheet conductance of one layer [S per square]."""
|
|
return (self.layers[layer_index].thickness_nm * 1e-9) / self.rho_ohm_m
|
|
|
|
def copper_bbox(self) -> tuple[int, int, int, int]:
|
|
xs = [p.outline[:, 0] for l in self.layers for p in l.polygons]
|
|
ys = [p.outline[:, 1] for l in self.layers for p in l.polygons]
|
|
tol = 1_000.0
|
|
for seg in self.tracks:
|
|
# exact stroke bbox: centerline extrema + half width (round
|
|
# caps); a chord-tessellated outline undershoots arc and cap
|
|
# extrema by up to its sagitta tolerance
|
|
pts = seg.centerline(tol)
|
|
r = seg.width_nm / 2.0 + tol
|
|
xs.append(np.array([pts[:, 0].min() - r, pts[:, 0].max() + r]))
|
|
ys.append(np.array([pts[:, 1].min() - r, pts[:, 1].max() + r]))
|
|
x = np.concatenate(xs)
|
|
y = np.concatenate(ys)
|
|
return int(x.min()), int(y.min()), int(x.max()), int(y.max())
|
|
|
|
|
|
def _arc_params(start, mid, end) -> tuple[float, float, float, float, float] | None:
|
|
"""Circle through three points: (cx, cy, r, a0, sweep) with a0 the
|
|
start angle and sweep signed; None if the points are collinear."""
|
|
sx, sy = float(start[0]), float(start[1])
|
|
mx, my = float(mid[0]), float(mid[1])
|
|
ex, ey = float(end[0]), float(end[1])
|
|
|
|
d = 2.0 * (sx * (my - ey) + mx * (ey - sy) + ex * (sy - my))
|
|
chord = math.hypot(ex - sx, ey - sy)
|
|
if abs(d) < 1e-9 * max(chord, 1.0):
|
|
return None
|
|
cx = ((sx**2 + sy**2) * (my - ey) + (mx**2 + my**2) * (ey - sy)
|
|
+ (ex**2 + ey**2) * (sy - my)) / d
|
|
cy = ((sx**2 + sy**2) * (ex - mx) + (mx**2 + my**2) * (sx - ex)
|
|
+ (ex**2 + ey**2) * (mx - sx)) / d
|
|
r = math.hypot(sx - cx, sy - cy)
|
|
|
|
a0 = math.atan2(sy - cy, sx - cx)
|
|
a1 = math.atan2(my - cy, mx - cx)
|
|
a2 = math.atan2(ey - cy, ex - cx)
|
|
two_pi = 2.0 * math.pi
|
|
d01 = (a1 - a0) % two_pi
|
|
d02 = (a2 - a0) % two_pi
|
|
sweep = d02 if d01 <= d02 else d02 - two_pi
|
|
return cx, cy, r, a0, sweep
|
|
|
|
|
|
def _n_arc_segments(sweep_abs: float, r: float, tol_nm: float) -> int:
|
|
"""Segments needed to keep the sagitta of each chord <= tol_nm."""
|
|
tol = min(tol_nm, 0.999 * r)
|
|
dtheta_max = 2.0 * math.acos(1.0 - tol / r)
|
|
return max(2, int(math.ceil(sweep_abs / dtheta_max)))
|
|
|
|
|
|
def arc_points(start, mid, end, tol_nm: float) -> np.ndarray:
|
|
"""Tessellate a start/mid/end arc into points from start (inclusive)
|
|
to end (exclusive), max sagitta <= tol_nm. Collinear input degrades
|
|
to just the start point (straight segment)."""
|
|
params = _arc_params(start, mid, end)
|
|
if params is None:
|
|
return np.array([[start[0], start[1]]], dtype=np.int64)
|
|
cx, cy, r, a0, sweep = params
|
|
n = _n_arc_segments(abs(sweep), r, tol_nm)
|
|
ks = np.arange(n)
|
|
angs = a0 + sweep * ks / n
|
|
pts = np.stack([cx + r * np.cos(angs), cy + r * np.sin(angs)], axis=1)
|
|
return np.round(pts).astype(np.int64)
|
|
|
|
|
|
def capsule_ring(x1: int, y1: int, x2: int, y2: int, width_nm: int,
|
|
tol_nm: float) -> np.ndarray:
|
|
"""Outline (open ring, int64 nm) of a straight track segment: a
|
|
rectangle with semicircular end caps; a circle for a zero-length
|
|
segment. Cap sagitta <= tol_nm."""
|
|
r = width_nm / 2.0
|
|
dx, dy = float(x2 - x1), float(y2 - y1)
|
|
length = math.hypot(dx, dy)
|
|
n = _n_arc_segments(math.pi, r, tol_nm)
|
|
if length < 1.0:
|
|
angs = np.linspace(0.0, 2.0 * math.pi, 2 * n, endpoint=False)
|
|
pts = np.stack([x1 + r * np.cos(angs), y1 + r * np.sin(angs)],
|
|
axis=1)
|
|
return np.round(pts).astype(np.int64)
|
|
ux, uy = dx / length, dy / length
|
|
a0 = math.atan2(ux, -uy) # angle of the left normal
|
|
ks = np.arange(n + 1)
|
|
cap2 = a0 - ks * math.pi / n # +normal -> -normal, around end
|
|
cap1 = a0 - (ks + n) * math.pi / n # -normal -> +normal, around start
|
|
pts = np.concatenate([
|
|
np.stack([x2 + r * np.cos(cap2), y2 + r * np.sin(cap2)], axis=1),
|
|
np.stack([x1 + r * np.cos(cap1), y1 + r * np.sin(cap1)], axis=1),
|
|
])
|
|
return np.round(pts).astype(np.int64)
|
|
|
|
|
|
def arc_band_ring(start, mid, end, width_nm: int, tol_nm: float) -> np.ndarray:
|
|
"""Outline of an arc track: the annular band of the given width
|
|
around the start/mid/end centerline, with semicircular end caps.
|
|
Collinear input degrades to the straight capsule."""
|
|
params = _arc_params(start, mid, end)
|
|
if params is None:
|
|
return capsule_ring(start[0], start[1], end[0], end[1], width_nm,
|
|
tol_nm)
|
|
cx, cy, r, a0, sweep = params
|
|
w2 = width_nm / 2.0
|
|
router = r + w2
|
|
rinner = max(r - w2, 0.0)
|
|
sgn = 1.0 if sweep >= 0 else -1.0
|
|
a1 = a0 + sweep
|
|
m = _n_arc_segments(abs(sweep), router, tol_nm)
|
|
ncap = _n_arc_segments(math.pi, w2, tol_nm)
|
|
ks = np.arange(m + 1)
|
|
|
|
th = a0 + sweep * ks / m # outer arc, start -> end
|
|
parts = [np.stack([cx + router * np.cos(th),
|
|
cy + router * np.sin(th)], axis=1)]
|
|
ex_, ey_ = cx + r * math.cos(a1), cy + r * math.sin(a1)
|
|
ca = a1 + sgn * math.pi * np.arange(1, ncap) / ncap # end cap, bulges
|
|
parts.append(np.stack([ex_ + w2 * np.cos(ca), # along exit tangent
|
|
ey_ + w2 * np.sin(ca)], axis=1))
|
|
if rinner > 0:
|
|
th = a1 - sweep * ks / m # inner arc, end -> start
|
|
parts.append(np.stack([cx + rinner * np.cos(th),
|
|
cy + rinner * np.sin(th)], axis=1))
|
|
else:
|
|
parts.append(np.array([[cx, cy]])) # band swallows the center
|
|
sx_, sy_ = cx + r * math.cos(a0), cy + r * math.sin(a0)
|
|
ca = a0 + math.pi + sgn * math.pi * np.arange(1, ncap) / ncap
|
|
parts.append(np.stack([sx_ + w2 * np.cos(ca), # start cap, bulges
|
|
sy_ + w2 * np.sin(ca)], axis=1)) # backwards
|
|
return np.round(np.concatenate(parts)).astype(np.int64)
|
|
|
|
|
|
def linearize_ring(nodes: list, tol_nm: float) -> np.ndarray:
|
|
"""nodes: list of ('pt', (x, y)) or ('arc', (start, mid, end)) tuples,
|
|
already in board nm. Returns an (N, 2) int64 open ring."""
|
|
parts = []
|
|
for kind, data in nodes:
|
|
if kind == "pt":
|
|
parts.append(np.array([[data[0], data[1]]], dtype=np.int64))
|
|
elif kind == "arc":
|
|
parts.append(arc_points(data[0], data[1], data[2], tol_nm))
|
|
else:
|
|
raise ValueError(f"unknown polyline node kind: {kind}")
|
|
ring = np.concatenate(parts, axis=0)
|
|
if len(ring) > 1 and (ring[0] == ring[-1]).all():
|
|
ring = ring[:-1]
|
|
return ring
|
|
|
|
|
|
# --- JSON dump / load -------------------------------------------------------
|
|
|
|
def _poly_to_json(p: Polygon) -> dict:
|
|
return {"outline": p.outline.tolist(), "holes": [h.tolist() for h in p.holes]}
|
|
|
|
|
|
def _poly_from_json(d: dict) -> Polygon:
|
|
return Polygon(outline=np.asarray(d["outline"], dtype=np.int64),
|
|
holes=[np.asarray(h, dtype=np.int64) for h in d["holes"]])
|
|
|
|
|
|
def _electrode_to_json(e: Electrode) -> dict:
|
|
return {
|
|
"rect": vars(e.rect) | {},
|
|
"contact": e.contact,
|
|
"label": e.label,
|
|
"polygons": (None if e.polygons is None
|
|
else [_poly_to_json(poly) for poly in e.polygons]),
|
|
}
|
|
|
|
|
|
def _electrode_from_json(d: dict) -> Electrode:
|
|
return Electrode(
|
|
rect=_rect_from_json(d["rect"]),
|
|
contact=d.get("contact", "all"),
|
|
label=d.get("label", "rect"),
|
|
polygons=(None if d.get("polygons") is None
|
|
else [_poly_from_json(pd) for pd in d["polygons"]]),
|
|
)
|
|
|
|
|
|
def problem_to_json(p: Problem) -> dict:
|
|
return {
|
|
"schema_version": JSON_SCHEMA_VERSION,
|
|
"board_path": p.board_path,
|
|
"net_name": p.net_name,
|
|
"rho_ohm_m": p.rho_ohm_m,
|
|
"plating_nm": p.plating_nm,
|
|
"thickness_source": p.thickness_source,
|
|
"electrodes1": [_electrode_to_json(e) for e in p.electrodes1],
|
|
"electrodes2": [_electrode_to_json(e) for e in p.electrodes2],
|
|
"layers": [
|
|
{
|
|
"layer_name": l.layer_name,
|
|
"thickness_nm": l.thickness_nm,
|
|
"z_nm": l.z_nm,
|
|
"polygons": [_poly_to_json(poly) for poly in l.polygons],
|
|
}
|
|
for l in p.layers
|
|
],
|
|
"vias": [vars(v) | {} for v in p.vias],
|
|
"tracks": [
|
|
{"layer_name": s.layer_name, "points": s.points.tolist(),
|
|
"width_nm": s.width_nm}
|
|
for s in p.tracks
|
|
],
|
|
"buildups": [
|
|
{"layer_name": b.layer_name,
|
|
"polygons": [_poly_to_json(poly) for poly in b.polygons]}
|
|
for b in p.buildups
|
|
],
|
|
"solder_thickness_nm": p.solder_thickness_nm,
|
|
"solder_rho_ohm_m": p.solder_rho_ohm_m,
|
|
"extra_cu_nm": p.extra_cu_nm,
|
|
"vias_capped": p.vias_capped,
|
|
"cap_plating_nm": p.cap_plating_nm,
|
|
}
|
|
|
|
|
|
def _rect_from_json(rd: dict) -> Rect:
|
|
return Rect(int(rd["x0"]), int(rd["y0"]), int(rd["x1"]), int(rd["y1"]),
|
|
rd["layer_name"])
|
|
|
|
|
|
def problem_from_json(d: dict) -> Problem:
|
|
version = d.get("schema_version", 1)
|
|
if version == 1:
|
|
# v1: single layer, no vias, rect electrodes
|
|
return Problem(
|
|
board_path=d["board_path"],
|
|
net_name=d["net_name"],
|
|
rho_ohm_m=float(d["rho_ohm_m"]),
|
|
plating_nm=18_000,
|
|
layers=[LayerFill(
|
|
layer_name=d["layer_name"],
|
|
thickness_nm=int(d["thickness_nm"]),
|
|
z_nm=0,
|
|
polygons=[_poly_from_json(pd) for pd in d["polygons"]],
|
|
)],
|
|
vias=[],
|
|
electrodes1=[Electrode(rect=_rect_from_json(d["rect1"]))],
|
|
electrodes2=[Electrode(rect=_rect_from_json(d["rect2"]))],
|
|
thickness_source=d.get("thickness_source", "unknown"),
|
|
)
|
|
return Problem(
|
|
board_path=d["board_path"],
|
|
net_name=d["net_name"],
|
|
rho_ohm_m=float(d["rho_ohm_m"]),
|
|
plating_nm=int(d["plating_nm"]),
|
|
layers=[
|
|
LayerFill(
|
|
layer_name=ld["layer_name"],
|
|
thickness_nm=int(ld["thickness_nm"]),
|
|
z_nm=int(ld["z_nm"]),
|
|
polygons=[_poly_from_json(pd) for pd in ld["polygons"]],
|
|
)
|
|
for ld in d["layers"]
|
|
],
|
|
vias=[
|
|
ViaLink(x=int(vd["x"]), y=int(vd["y"]), drill_nm=int(vd["drill_nm"]),
|
|
z_top_nm=int(vd["z_top_nm"]), z_bot_nm=int(vd["z_bot_nm"]),
|
|
kind=vd.get("kind", "via"),
|
|
pad_nm=int(vd.get("pad_nm", 0)))
|
|
for vd in d["vias"]
|
|
],
|
|
electrodes1=(
|
|
[_electrode_from_json(ed) for ed in d["electrodes1"]]
|
|
if version >= 3 else [_electrode_from_json(d["electrode1"])]),
|
|
electrodes2=(
|
|
[_electrode_from_json(ed) for ed in d["electrodes2"]]
|
|
if version >= 3 else [_electrode_from_json(d["electrode2"])]),
|
|
thickness_source=d.get("thickness_source", "unknown"),
|
|
buildups=[
|
|
SurfaceBuildup(
|
|
layer_name=bd["layer_name"],
|
|
polygons=[_poly_from_json(pd) for pd in bd["polygons"]])
|
|
for bd in d.get("buildups", [])
|
|
],
|
|
solder_thickness_nm=int(d.get("solder_thickness_nm", 50_000)),
|
|
solder_rho_ohm_m=float(d.get("solder_rho_ohm_m", 1.32e-7)),
|
|
extra_cu_nm=int(d.get("extra_cu_nm", 0)),
|
|
tracks=[
|
|
TrackSeg(layer_name=td["layer_name"],
|
|
points=np.asarray(td["points"], dtype=np.int64),
|
|
width_nm=int(td["width_nm"]))
|
|
for td in d.get("tracks", []) # <= v4: baked into polygons
|
|
],
|
|
vias_capped=bool(d.get("vias_capped", True)),
|
|
cap_plating_nm=int(d.get("cap_plating_nm", 15_000)),
|
|
)
|
|
|
|
|
|
def save_problem(p: Problem, path: Path) -> None:
|
|
path.write_text(json.dumps(problem_to_json(p)), encoding="utf-8")
|
|
|
|
|
|
def load_problem(path: Path) -> Problem:
|
|
return problem_from_json(json.loads(Path(path).read_text(encoding="utf-8")))
|