From 0afb55216b1bbee4635010429e074560416e2920 Mon Sep 17 00:00:00 2001 From: grabowski Date: Wed, 15 Jul 2026 19:43:21 +0700 Subject: [PATCH] Fix swarm-review findings: empty-layer crashes, teardrop fills, Jmag - 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) --- fill_resistance/adaptive.py | 16 +++++++++---- fill_resistance/board_io.py | 8 ++++--- fill_resistance/geometry.py | 11 ++++++--- fill_resistance/raster.py | 11 +++++---- fill_resistance/solver.py | 46 ++++++++++++++++++++++++++++++------- tests/test_quadtree.py | 8 +++---- 6 files changed, 73 insertions(+), 27 deletions(-) diff --git a/fill_resistance/adaptive.py b/fill_resistance/adaptive.py index 3eef0a8..799b407 100644 --- a/fill_resistance/adaptive.py +++ b/fill_resistance/adaptive.py @@ -130,14 +130,13 @@ def run_solve_adaptive(problem: Problem, stack: RasterStack, cxg[offs[li]:offs[li + 1]] = cxl cyg[offs[li]:offs[li + 1]] = cyl sig_leaf = np.full(g_.n, sigmas[li]) - t_m = problem.layers[li].thickness_nm * 1e-9 s2d = sv._sigma_2d(stack, li, sigmas[li], sigma_buildup) fine = g_.size == 1 if s2d is not None and fine.any(): sig_leaf[fine] = s2d[g_.y0[fine], g_.x0[fine]] - # J reference thickness: same convention as the uniform grid - teq_leaves.append(sig_leaf * problem.rho_ohm_m if s2d is not None - else np.full(g_.n, t_m)) + # J reference thickness: conduction-equivalent copper (= geometric + # t at DC, skin-reduced at AC), same convention as the uniform grid + teq_leaves.append(sig_leaf * problem.rho_ohm_m) sig_leaves.append(sig_leaf) chainleaf = np.zeros(g_.n, dtype=bool) if stack.chain is not None and fine.any(): @@ -156,7 +155,7 @@ def run_solve_adaptive(problem: Problem, stack: RasterStack, ee.append(np.full(len(ia), li, dtype=np.int16)) if stack.chain_edges is not None and len(stack.chain_edges[0]): - ca, cb, cg, cl = stack.chain_edges + ca, cb, cg, cl, _ = stack.chain_edges alive = stack.masks.ravel()[ca] & stack.masks.ravel()[cb] if alive.any(): na = np.empty(len(ca), dtype=np.int64) @@ -241,6 +240,8 @@ def run_solve_adaptive(problem: Problem, stack: RasterStack, dead_barrels=dead_barrels) e_delta, e_axis, e_layer = e_delta[sel], e_axis[sel], e_layer[sel] for li in range(L): + if grids[li].n == 0: + continue ids = grids[li].id_grid kept_cells = (ids >= 0) & keepn[offs[li] + np.maximum(ids, 0)] stack.masks[li] &= kept_cells @@ -389,6 +390,8 @@ def run_solve_adaptive(problem: Problem, stack: RasterStack, # (fine regions stay plain copper = fully resolved) stack.mesh = np.zeros_like(stack.masks) for li in range(L): + if grids[li].n == 0: + continue ids = grids[li].id_grid b = np.zeros_like(stack.masks[li]) b[:, 1:] |= ids[:, 1:] != ids[:, :-1] @@ -437,6 +440,9 @@ def run_solve_adaptive(problem: Problem, stack: RasterStack, cellP = Pnode[offs[li]:offs[li + 1]] \ / (g_.size.astype(float) ** 2 * h_m * h_m) Parea[li][m] = np.maximum(cellP, 0.0)[ids[m]] + # chain cells accumulate no leaf-face currents (their links carry + # axis -1): overlay the true 1D link density + sv.overlay_chain_density(stack, problem.rho_ohm_m, V3, J3) timings["postprocess_s"] = time.perf_counter() - t0 return sv.Result( diff --git a/fill_resistance/board_io.py b/fill_resistance/board_io.py index dee3bef..7fb85a7 100644 --- a/fill_resistance/board_io.py +++ b/fill_resistance/board_io.py @@ -292,7 +292,9 @@ def gather_net_fills(board: Board) -> dict[str, dict[str, list[Polygon]]]: """net -> layer_name -> merged fill polygons (non-empty only).""" fills: dict[str, dict[str, list[Polygon]]] = {} for zone in board.get_zones(): - if zone.type != ZoneType.ZT_COPPER: + # teardrop fills are conducting copper too, but KiCad types them + # ZT_TEARDROP instead of ZT_COPPER + if zone.type not in (ZoneType.ZT_COPPER, ZoneType.ZT_TEARDROP): continue net = zone.net.name if zone.net is not None else "" for layer, polys in zone.filled_polygons.items(): @@ -392,8 +394,8 @@ def gather_mask_buildups(board: Board) -> dict[str, list[Polygon]]: def any_zone_unfilled(board: Board) -> bool: - return any(z.type == ZoneType.ZT_COPPER and not z.filled - for z in board.get_zones()) + return any(z.type in (ZoneType.ZT_COPPER, ZoneType.ZT_TEARDROP) + and not z.filled for z in board.get_zones()) def refill(board: Board) -> None: diff --git a/fill_resistance/geometry.py b/fill_resistance/geometry.py index f679813..f757035 100644 --- a/fill_resistance/geometry.py +++ b/fill_resistance/geometry.py @@ -159,10 +159,15 @@ class Problem: 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: - ring = seg.outline(100_000.0) # coarse tol: bbox only - xs.append(ring[:, 0]) - ys.append(ring[:, 1]) + # 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()) diff --git a/fill_resistance/raster.py b/fill_resistance/raster.py index c308908..270d930 100644 --- a/fill_resistance/raster.py +++ b/fill_resistance/raster.py @@ -40,8 +40,9 @@ class RasterStack: # (mask opening ∩ copper) chain: np.ndarray | None = None # bool (L, ny, nx): cells that are # copper only through a 1D trace chain - chain_edges: tuple | None = None # (a, b, g_dc, layer) arrays: explicit - # DC conductances of the chain links + chain_edges: tuple | None = None # (a, b, g_dc, layer, dl_m) arrays: + # explicit DC conductances and link + # lengths of the chain links thick_scale: np.ndarray | None = None # float (L, ny, nx): per-cell # copper-thickness factor (via # mouths: cap-thin or partially @@ -320,7 +321,7 @@ def _build_chains(stack: RasterStack, problem: Problem, h = stack.h_nm regular = stack.masks.copy() chain = np.zeros_like(stack.masks) - aa, bb, gg, ll = [], [], [], [] + aa, bb, gg, ll, dd = [], [], [], [], [] for li, seg in narrow: pts = seg.centerline(0.2 * h) d = np.hypot(*np.diff(pts, axis=0).T) @@ -352,12 +353,14 @@ def _build_chains(stack: RasterStack, problem: Problem, bb.append(li * plane + int(ci[k + 1]) * nx + int(cj[k + 1])) gg.append(g0 / dl) ll.append(li) + dd.append(dl) stack.chain = chain & ~regular stack.masks |= stack.chain stack.chain_edges = (np.asarray(aa, dtype=np.int64), np.asarray(bb, dtype=np.int64), np.asarray(gg, dtype=float), - np.asarray(ll, dtype=np.int64)) + np.asarray(ll, dtype=np.int64), + np.asarray(dd, dtype=float)) return len(aa) diff --git a/fill_resistance/solver.py b/fill_resistance/solver.py index c019a5a..47f3fea 100644 --- a/fill_resistance/solver.py +++ b/fill_resistance/solver.py @@ -233,7 +233,7 @@ def build_edges(stack: RasterStack, problem: Problem, sigmas: list[float], vv.append(np.array([vi], dtype=np.int32)) if stack.chain_edges is not None and len(stack.chain_edges[0]): - ca, cb, cg, cl = stack.chain_edges + ca, cb, cg, cl, _ = stack.chain_edges alive = stack.masks.ravel()[ca] & stack.masks.ravel()[cb] if alive.any(): # skin correction: scale like the layer's sheet conductance @@ -440,16 +440,16 @@ def _face_current_density(V2: np.ndarray, mask2: np.ndarray, sigma: float, sig2d: np.ndarray | None = None, rho: float | None = None) -> np.ndarray: """|J| (A/m^2) for one layer from face currents; V2 in volts. - With a per-cell conductance map (buildup), face currents use the - harmonic mean and J is referenced to the conductance-equivalent - copper thickness t_eq = sigma_cell * rho (equals the geometric t for - plain DC copper).""" + J is referenced to the conductance-equivalent copper thickness + t_eq = sigma_cell * rho: the geometric t for plain DC copper, the + skin-reduced conducting cross-section at AC. With a per-cell + conductance map (buildup), face currents use the harmonic mean.""" ny, nx = mask2.shape face_x = mask2[:, :-1] & mask2[:, 1:] face_y = mask2[:-1, :] & mask2[1:, :] if sig2d is None: wx = wy = sigma - teq = np.full((ny, nx), t_m) + teq = np.full((ny, nx), sigma * rho if rho is not None else t_m) else: wx = 2.0 * sig2d[:, :-1] * sig2d[:, 1:] / (sig2d[:, :-1] + sig2d[:, 1:]) wy = 2.0 * sig2d[:-1, :] * sig2d[1:, :] / (sig2d[:-1, :] + sig2d[1:, :]) @@ -468,6 +468,31 @@ def _face_current_density(V2: np.ndarray, mask2: np.ndarray, sigma: float, return Jmag +def overlay_chain_density(stack: RasterStack, rho: float, V3: np.ndarray, + J3: np.ndarray) -> None: + """Fill chain (sub-resolution trace) cells of J3 with the true 1D + link current density |dV| / (rho * dl), referenced to the + conduction-equivalent trace cross-section: the AC scaling of the + link conductance and of the cross-section cancel, so the expression + holds at any frequency. V3/J3 are the display-scaled (L, ny, nx) + maps; chain cells carry the max density of their attached links.""" + if stack.chain is None or stack.chain_edges is None \ + or not len(stack.chain_edges[0]): + return + ca, cb, _, _, cdl = stack.chain_edges + mflat = stack.masks.reshape(-1) + alive = mflat[ca] & mflat[cb] + if not alive.any(): + return + V3f = np.nan_to_num(V3.reshape(-1)) + Jl = np.abs(V3f[ca] - V3f[cb]) / (rho * cdl) + Jc = np.zeros(mflat.size) + np.maximum.at(Jc, ca[alive], Jl[alive]) + np.maximum.at(Jc, cb[alive], Jl[alive]) + fill = stack.chain & stack.masks + J3[fill] = Jc.reshape(J3.shape)[fill] + + def _equipotential_core(state: np.ndarray, edges: Edges): """Dirichlet solve on any node space (fine cells or leaves): state codes 0 off / 1 free / 2 V+ / 3 V-. Returns (Vflat_unit, R, I1, I2, @@ -719,17 +744,22 @@ def run_solve(problem: Problem, stack: RasterStack, e1: np.ndarray, parts2 or [], Ie, edges, e2.ravel(), s, i_test, contact_model, int(e2.sum())) - # embedded potential + per-layer current density @ I_test + # embedded potential + per-layer current density @ I_test; chain + # cells have no sheet faces in the model, so keep them out of the + # face computation and overlay their true 1D link density instead V3 = np.full((L, ny, nx), np.nan) V3[stack.masks] = Vflat.reshape(L, ny, nx)[stack.masks] * s + sheet = stack.masks if stack.chain is None \ + else stack.masks & ~stack.chain J3 = np.stack([ _face_current_density( - np.nan_to_num(V3[li]), stack.masks[li], sigmas[li], + np.nan_to_num(V3[li]), sheet[li], sigmas[li], h_m, problem.layers[li].thickness_nm * 1e-9, sig2d=_sigma_2d(stack, li, sigmas[li], sigma_buildup), rho=problem.rho_ohm_m) for li in range(L) ]) + overlay_chain_density(stack, problem.rho_ohm_m, V3, J3) timings["postprocess_s"] = time.perf_counter() - t0 return Result( diff --git a/tests/test_quadtree.py b/tests/test_quadtree.py index ba4408b..c751db0 100644 --- a/tests/test_quadtree.py +++ b/tests/test_quadtree.py @@ -56,13 +56,13 @@ def test_uniform_limit_graph_identical(): ny, nx = stack.shape2d flat_of_leaf = grid.y0.astype(np.int64) * nx + grid.x0 ia, ib, g = quadtree.leaf_edges(grid, p.sigma_s(0)) - ours = np.sort(np.stack([ + ours = np.stack([ np.minimum(flat_of_leaf[ia], flat_of_leaf[ib]), - np.maximum(flat_of_leaf[ia], flat_of_leaf[ib])], axis=1), axis=0) + np.maximum(flat_of_leaf[ia], flat_of_leaf[ib])], axis=1) edges = solver.build_edges(stack, p, [p.sigma_s(0)]) - ref = np.sort(np.stack([np.minimum(edges.a, edges.b), - np.maximum(edges.a, edges.b)], axis=1), axis=0) + ref = np.stack([np.minimum(edges.a, edges.b), + np.maximum(edges.a, edges.b)], axis=1) assert ours.shape == ref.shape assert np.array_equal(np.sort(ours.view("i8,i8"), order=["f0", "f1"], axis=0),