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

@@ -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
+        # J reference thickness: conduction-equivalent copper (= geometric
-        teq_leaves.append(sig_leaf * problem.rho_ohm_m if s2d is not None
+        # t at DC, skin-reduced at AC), same convention as the uniform grid
-                          else np.full(g_.n, t_m))
+        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(

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 "<no net>"
         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
+    return any(z.type in (ZoneType.ZT_COPPER, ZoneType.ZT_TEARDROP)
-               for z in board.get_zones())
+               and not z.filled for z in board.get_zones())
 
 
 def refill(board: Board) -> None:

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
+            # exact stroke bbox: centerline extrema + half width (round
-            xs.append(ring[:, 0])
+            # caps); a chord-tessellated outline undershoots arc and cap
-            ys.append(ring[:, 1])
+            # 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())

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
+    chain_edges: tuple | None = None    # (a, b, g_dc, layer, dl_m) arrays:
-                                        # DC conductances of the chain links
+                                        # 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)
 
 

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
+    J is referenced to the conductance-equivalent copper thickness
-    harmonic mean and J is referenced to the conductance-equivalent
+    t_eq = sigma_cell * rho: the geometric t for plain DC copper, the
-    copper thickness t_eq = sigma_cell * rho (equals the geometric t for
+    skin-reduced conducting cross-section at AC. With a per-cell
-    plain DC copper)."""
+    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(

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),
+    ref = np.stack([np.minimum(edges.a, edges.b),
-                            np.maximum(edges.a, edges.b)], axis=1), axis=0)
+                    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),