Release notes come from a file in the repo, not the tag

Deriving the body from the tag annotation does not survive CI: the
release action does not fall back to it (v1.2.0 published empty), and
checkout leaves the tag lightweight, so %(contents) yields the commit
message instead - which is what v1.2.1 first published. Keep the notes
at docs/release-notes/v<version>.md, pass that as body_path, and fail
the run when it is missing.

.gitea/workflows/build-pcm.yml

@@ -35,11 +35,28 @@ jobs:
             dist/*.zip
             dist/metadata-registry.json
 
-      - name: Create release with the zip attached
+      # The release body comes from a file in the repo: the action does
+      # not fall back to the tag annotation (v1.2.0 published empty), and
+      # reading the annotation here is unreliable - checkout leaves the
+      # tag lightweight, so %(contents) yields the commit message instead.
+      - name: Check the release notes exist
+        if: startsWith(github.ref, 'refs/tags/v')
+        run: |
+          notes="docs/release-notes/${GITHUB_REF_NAME}.md"
+          if [ ! -s "$notes" ]; then
+            echo "$notes is missing or empty - write the release notes" \
+                 "before tagging" >&2
+            exit 1
+          fi
+          cat "$notes"
+
+      - name: Create release with the zip, registry metadata and figures
         if: startsWith(github.ref, 'refs/tags/v')
         uses: akkuman/gitea-release-action@b8d9144f302c68610911db1aaf722708d5c02d94 # v1
         with:
+          body_path: docs/release-notes/${{ github.ref_name }}.md
           files: |
             dist/*.zip
             dist/metadata-registry.json
+            docs/img/*.png
           token: ${{ secrets.GITEA_TOKEN }}

.gitignore

@@ -3,3 +3,11 @@ __pycache__/
 *.pyc
 .pytest_cache/
 dist/
+
+# local AI-tooling artifacts, never publish
+.claude/
+.claude-flow/
+.swarm/
+.mcp.json
+CLAUDE.md
+ruvector.db

README.md

@@ -6,11 +6,23 @@ between two contacts, **single- or multi-layer**: the chosen net's fills
 solved as coupled finite-difference sheets linked by the net's **via
 and through-hole-pad barrels** (18 µm plating, configurable). At a user-set **frequency** the exact 1D foil/barrel
 skin-effect correction is applied (AC results are a rigorous lower
-bound — see *Model & limits*). Shows per-layer rasterized maps,
+bound; see *Model & limits*). Shows per-layer rasterized maps,
-potential, current density, and **power density**, reports **per-via
+potential, current density, and **power density**, and reports **per-via
 currents** (via ampacity!) and total dissipation at a **selectable test
 current**. PNGs + a text summary are saved per run.
 
+![Current density on a two-layer demo net](docs/img/demo-current.png)
+*Real output on a synthetic two-layer net: current from a soldered
+THT-pad contact (V+, injected at the drill-wall ring) squeezes past a
+notch in the F.Cu pour, transfers through the stitching-via field into
+the B.Cu pour and leaves at the V− lug. Per-via currents and the
+hottest via are reported.*
+
+![Potential on the two-layer demo net](docs/img/demo-potential.png)
+*The matching potential map with equipotential contour lines: they
+bunch where the field is strongest — nearly the whole 8.7 mV drop
+happens around the notch on F.Cu.*
+
 Uses the KiCad **IPC API** (`kicad-python` / `kipy`), not the deprecated
 SWIG API. Requires KiCad **10.0.1+**.
 
@@ -23,7 +35,7 @@ SWIG API. Requires KiCad **10.0.1+**.
    on Windows or `/usr/bin/python3` on Linux (after a 9→10 upgrade it
    can point at KiCad 9).
 3. **Deploy** (dev checkout; end users install the PCM zip instead, see
-   *Packaging*):
+   *Packaging / publishing*):
    ```powershell
    powershell -ExecutionPolicy Bypass -File deploy.ps1        # junction (dev)
    powershell -ExecutionPolicy Bypass -File deploy.ps1 -Mode Copy
@@ -35,18 +47,30 @@ SWIG API. Requires KiCad **10.0.1+**.
    ```
 4. **Restart KiCad**; first load builds the plugin venv (numpy, scipy,
    matplotlib, PySide6 — takes minutes; the Ω button appears when done).
-   If stuck: Preferences → Plugins → *Recreate Plugin Environment*.
+   If stuck: in the PCB editor, Preferences → *PCB Editor → Action
+   Plugins*, **right-click** the plugin's row → *Recreate Plugin
+   Environment* (context menu only — there is no button). Manual
+   equivalent: delete
+   `%LOCALAPPDATA%\kicad\10.0\python-environments\th.co.b4l.fill-resistance`
+   and restart KiCad.
 
 ## Usage
 
 1. Mark the current-injection terminals. Each terminal may have
-   **multiple parts** (all merged into one externally-bonded contact):
+   **multiple parts** (all merged into one externally bonded contact):
    - **V+ rectangles on `User.1`**, **V− rectangles on `User.2`**
      (marker layers, configurable via `ELECTRODE_POS_LAYER` /
      `ELECTRODE_NEG_LAYER`), any number per side, axis-aligned;
-   - **pads** (real copper shape; through-hole pad contacts all layers,
+   - **pads and vias** (SMD pad: real copper shape on its own layer;
-     SMD pad its own layer) — selected pads fill a side that has no
+     through-hole pads and vias become **barrel contacts**: the current
-     rectangles;
+     enters at the drill wall on every spanned layer, see below).
+     Selected pads/vias fill the side that has **no rectangles**, so
+     mixing both kinds is the everyday workflow: e.g. select **one
+     rectangle on `User.1`** (V+) **plus any number of pads / THT
+     holes** (Ctrl-click) — the pads together form the V− terminal
+     (a connector's pin group, a via cluster, …). All selected
+     pads/vias go to that one side; if both marker layers already
+     provide rectangles, selecting pads on top is an error;
    - legacy: exactly 2 selected contacts with no marker rectangles still
      works; empty selection scans the whole board's marker layers.
 2. **Select the contacts**, click the **Fill Resistance** Ω button.
@@ -55,13 +79,28 @@ SWIG API. Requires KiCad **10.0.1+**.
    ("All selected layers" = bolted-lug/through contact), the **test
    current**, and optionally a grid cell size. Multiple layers are coupled
    through the net's via/pad barrels automatically.
-4. Read R / voltage drop / total power in the figure titles and status
+4. Wait for the solve. Depending on board size, included layers, cell
+   size and your hardware it can take **considerable time** — large
+   multi-layer pours at fine cell sizes may run for minutes (on our
+   test setup a typical real-board run finishes in ≈ 8 s). Then read
+   R / voltage drop / total power in the figure titles and status
    bar. Outputs land in `<board dir>\fill_res_results\<timestamp>\`:
    per-layer `1_raster_map` / `2_potential` / `3_current_density` /
    `4_power_density` PNGs, `summary.txt` (incl. the busiest vias with
    per-via current and dissipation, and the **current through each
    injection area** — computed flux with the equipotential model,
    prescribed area share with the uniform model), `geometry_dump.json`.
+5. **Experimental — overlays inside KiCad** (dialog checkbox, default
+   off; KiCad ≥ 10.0.1): after the solve, the per-layer **|J| heatmaps
+   are pushed into the open board** as unlocked reference images on
+   `User.9`…`User.12` (`OVERLAY_LAYERS`; enable them in Board Setup),
+   copper layers mapped in stackup order, top first. Toggle them in the
+   Appearance panel like any layer; opaque over copper, transparent
+   elsewhere, cold end lifted so it stays visible on the dark canvas.
+   Reference images never plot to gerbers. Every push **replaces all
+   reference images on those layers**, so don't store unrelated images
+   there. Also available headless:
+   `python tools/kicad_heatmap_overlay.py --net X --amps 10`.
 
 ## Model & limits
 
@@ -76,20 +115,56 @@ SWIG API. Requires KiCad **10.0.1+**.
   capped" checkbox** (default on, `VIAS_CAPPED`) the mouth carries a
   thin copper cap (`CAP_PLATING_UM = 15`, fab spec) on the **outer**
   layers and is an open hole on inner layers; unchecked, mouths are open
-  holes everywhere. Layer-to-layer the cap never matters at DC (it is in
+  holes everywhere. The fab caps only small vias: drills above the
-  parallel with the annular-ring contact, not in series) — the checkbox
+  dialog's **"capped up to drill"** threshold (default
+  `CAP_MAX_DRILL_MM = 0.5`) keep open mouths even with capping
+  selected. Layer-to-layer the cap never matters at DC (it is in
+  parallel with the annular-ring contact, not in series), so the checkbox
   only affects in-plane conduction across outer-layer mouths. Sub-cell
   mouths scale their cells' sheet conductance by the true covered
   fraction (4×4 supersampling), so coarse grids see the correct small
-  perturbation instead of a whole-cell hole. THT-pad copper and drills
+  perturbation instead of a whole-cell hole. Barrels are gathered in
-  remain outside the model; at f > 0 the thickness scaling is applied
+  **single-layer runs too** (drill mouths perforate a lone plane).
-  multiplicatively to the skin-corrected sheet conductance
+  **THT pads are fully modeled**: their exact copper shapes (incl.
+  oblong pads, fetched from KiCad; the outer shape stands in for inner
+  rings) are stamped onto every included layer, and every **populated**
+  pad carries its full **soldered joint** on its SOLDER side (opposite
+  the component; the component-side pad face stays bare). The hole
+  holds the **component lead** (a cylinder of drill −
+  `THT_LEAD_CLEARANCE_MM`, resistivity `THT_LEAD_RHO_OHM_M`, copper by
+  default; raise it for brass/steel leads) **plus solder** in the
+  remaining annulus, both in parallel with the plating. The filled
+  hole also conducts **in-plane on every spanned layer** (component
+  side and inner layers included): the mouth keeps its copper and
+  additionally carries the plug — lead disc plus solder bore — as
+  conduction-equivalent copper of the **full hole depth** (the pin
+  continues beyond both mouths, so each layer sees the whole plug
+  cross-section). The joint is side-symmetric except for the solder:
+  coat and cone on the solder side only. On the raster map these
+  mouths render in a darker tin color. The pad face
+  gets the average-thickness solder coat (exact pad shape) and the
+  protruding-lead cone (see barrel contacts below; on oblong pads the
+  cone tapers to the pad's short dimension). **Slotted (oval) holes**
+  keep their true stadium shape: the barrel wall, drill mouth, contact
+  ring and lead cone all follow the slot (rotated with the pad), and
+  the barrel conducts over the slot's real perimeter/bore area — not a
+  circle of the slot's long dimension. Whether a hole is a via or
+  a THT pad, the owning footprint's side, and its **Do not populate**
+  flag are all read from KiCad. **DNP pads** get an **open hole** and
+  a plating-only barrel, no joint. At f > 0 the thickness scaling is
+  applied multiplicatively to the skin-corrected sheet conductance
   (approximation). Per layer a barrel attaches to
   the fill cell under it, or to the nearest copper cell within the pad
   footprint plus one grid cell — fills joined by **thermal-relief
   spokes** still connect; wider antipads do not, and the barrel bridges
   the layers above/below with the full barrel length. Barrels that reach
   fill on fewer than two layers carry no current and are reported.
+
+  ![How drilled holes are modeled — cross-section](docs/img/hole-model.png)
+  *The four hole types: capped small via, open large via, populated THT
+  pad with its full solder joint (lead ∥ solder ∥ plating in the hole,
+  one-sided pad coat, protruding-lead solder cone), and a DNP THT pad.*
+
 - The net's **traces** (straight and arc tracks, exact outline polygons
   incl. rounded ends) conduct together with the fills — dialog checkbox,
   on by default (`INCLUDE_TRACKS`). Traces narrower than
@@ -98,8 +173,13 @@ SWIG API. Requires KiCad **10.0.1+**.
   series resistance carries no discretization error and no cell-size
   tuning is needed for thin traces. 1D-modeled traces show potential,
   power density, and |J| (the true in-trace density from the link
-  currents, |ΔV|/(ρ·Δl)). Pad copper other than the selected
+  currents, |ΔV|/(ρ·Δl)). Pad copper is part of the conductor: THT pad
-  contacts is still **not** part of the conductor model.
+  shapes are stamped on every included layer (see above), **SMD** pad
+  shapes on their own layer (`INCLUDE_SMD_PADS`) — pads are the
+  junctions where traces and thermal-relief spokes actually meet, so
+  without them a multi-track junction necks down to the accidental
+  overlap of the track ends. Dead-end pads (component terminals) are
+  dropped with the other copper not connected to both contacts.
 - **Solder buildup on mask openings** (dialog checkbox, **off by
   default**; `INCLUDE_MASK_BUILDUP`): zones drawn on `F.Mask`/`B.Mask`
   are treated as mask openings that collect `SOLDER_THICKNESS_UM`
@@ -110,6 +190,29 @@ SWIG API. Requires KiCad **10.0.1+**.
   interface faces use harmonic-mean conductances. Buildup areas render
   tin-gray on the raster map; |J| in them is referenced to the
   conductance-equivalent copper thickness.
+- **Barrel contacts**: a selected **via or through-hole pad** injects at
+  the **drill-wall ring** on every layer the barrel spans — the current
+  physically enters through the lead/wire soldered into the hole, so
+  the spreading resistance across the pad and surrounding pour is part
+  of the result (both contact models; verified against
+  R = ρ/(π·t)·acosh(d/2a) for two circular contacts on a sheet).
+  Slotted holes inject along the stadium-shaped slot wall. A
+  soldered **THT joint** additionally assumes the **hole is filled with
+  solder** (core in parallel with the plating) and the **pad face on
+  the solder side carries an average-thickness solder coat**
+  (`SOLDER_THICKNESS_UM`, 50 µm) over the modeled copper under the pad
+  shape — the solder side is the side opposite the component (taken
+  from the owning footprint; assumed `B.Cu` if it cannot be found),
+  and the component-side pad face stays bare. There the **clipped
+  lead protrudes** `THT_LEAD_PROTRUSION_MM` (1.5 mm, 0 = off)
+  and a **solder cone** wraps
+  it: full protrusion height at the drill wall, tapering linearly to
+  zero at the pad edge, applied as extra conduction-equivalent copper
+  per cell. The tall solder column at the wall pulls the joint
+  vicinity to lead potential — equivalent to extending the barrel wall
+  vertically — while the taper carries the radial spreading. To model
+  a probe pressed onto the pad face instead, draw a marker rectangle
+  over the pad.
 - **Contact models** (dialog / `CONTACT_MODEL`): default **uniform
   injection** — a conductor pressed on top feeds the current orthogonally
   with uniform surface density, so |J| ramps across the contact area
@@ -120,13 +223,19 @@ SWIG API. Requires KiCad **10.0.1+**.
   terminals (e.g. planes joined only through the bolted lugs), only the
   equipotential model is well-defined; the uniform model stops with an
   error instead of prescribing an arbitrary split.
+
+  ![The two contact models bracket a real contact](docs/img/contact-models.png)
+  *|J| around the same 3×3 mm contact under both models: the ideal
+  bonded lug crowds the current at the contact edges (no in-sheet
+  current inside an equipotential region), the pressed conductor ramps
+  it across the contact area.*
 - Fields are reported at the dialog's test current; power scales with I².
 - **Skin effect (f > 0)**: per-layer effective sheet resistance from the
   exact 1D foil-diffusion solution `Zs = τρ·coth(τt)`, `τ = (1+j)/δ`
   (`SKIN_SIDES = 1` in config: plane facing a return plane; `2` =
   isolated foil), and the analogous correction for the 18 µm barrel wall.
   Enter one frequency per run (e.g. a switching harmonic, with its RMS
-  amplitude as the test current) — suffixes `k`/`M` accepted.
+  amplitude as the test current); suffixes `k`/`M` are accepted.
   **Caveat:** only through-thickness crowding is modeled. Lateral
   (proximity-effect) redistribution needs a magneto-quasistatic solver
   and is not captured — since the resistance-driven distribution is the
@@ -138,15 +247,15 @@ SWIG API. Requires KiCad **10.0.1+**.
   not the geometric foil thickness.
 - 5-point FDM per layer on an auto-sized shared grid (~2 M fine cells
   with the uniform grid; ~8 M with the adaptive grid, whose unknown
-  count no longer scales with them). Direct sparse solve up to 500 k
+  count no longer scales with the fine-cell count). Direct sparse solve
-  unknowns, AMG-preconditioned CG (pyamg) above — Jacobi-CG if pyamg is
+  up to 500 k unknowns, AMG-preconditioned CG (pyamg) above (Jacobi-CG
-  missing. Discretization error typically ≲ 2 % at defaults — halve the
+  if pyamg is missing). Discretization error typically ≲ 2 % at
-  cell size and compare to judge convergence.
+  defaults; halve the cell size and compare to judge convergence.
 - **Adaptive cells** (dialog checkbox, **on by default**;
   `ADAPTIVE_CELLS`):
   solves on a 2:1-balanced quadtree — fine cells at copper boundaries,
   electrodes, traces, via mouths and buildup, blocks up to
-  `ADAPTIVE_MAX_CELL_UM` (2 mm) in plane interiors (`ADAPTIVE_GUARD`
+  `ADAPTIVE_MAX_CELL_UM` (1 mm) in plane interiors (`ADAPTIVE_GUARD`
   sets the clearance a block needs to grow). The **minimum element size
   is the grid cell size itself** (auto / dialog / `CELL_UM_OVERRIDE`);
   the uniform limit reproduces the normal grid exactly. Large
@@ -158,25 +267,40 @@ SWIG API. Requires KiCad **10.0.1+**.
   uniform grid ≲ 0.03 %, with the power-balance identity intact. All
   fields are expanded back to the fine grid for the maps and reports.
 
+  ![Rasterized map with the adaptive mesh overlay](docs/img/demo-raster.png)
+  *The raster map of the demo net: quadtree leaves drawn on the copper
+  (fine at boundaries, electrodes, via mouths and pads; coarse blocks
+  in plane interiors), the tin-gray solder coat of the THT-pad contact
+  P1, and the via field with its pad copper.*
+
+**Measured vs. computed**: we tested the plugin on a few real boards
+against a UT3513+ micro-ohm meter; the measured resistances were within
+±20 % of the computed values. We attribute the deviation to
+imperfections of the testing setup (probe placement and probe contact
+resistance vs. the ideal modeled contacts) and to manufacturing
+inaccuracies — actual copper and plating thicknesses routinely deviate
+from nominal. Relative comparisons between layout variants are
+accordingly more trustworthy than absolute numbers.
+
 ## Offline / development
 
 Every run writes `geometry_dump.json`; re-solve without KiCad:
 
 ```powershell
-.venv\Scripts\python.exe -m fill_resistance.standalone dump.json `
+uv run python -m fill_resistance.standalone dump.json `
     [--current 40] [--cell-um 50] [--layers F.Cu,In1.Cu] [--no-show] `
     [--out DIR] [--force-iterative]
 ```
 
-Dev environment, tests, headless extraction (Windows shown; on
+Dev environment, tests, headless extraction — [uv](https://docs.astral.sh/uv/)
-Linux/macOS use `.venv/bin/python`):
+manages the venv from `pyproject.toml`/`uv.lock` (`requirements.txt`
+stays: KiCad builds the plugin's runtime venv from it):
 
 ```powershell
-uv venv --python 3.11 .venv
+uv sync                                              # one-time env setup
-uv pip install --python .venv\Scripts\python.exe kicad-python numpy scipy pyamg matplotlib pytest
+uv run pytest -q                                     # incl. exact analytic cases
-.venv\Scripts\python.exe -m pytest tests -q          # incl. exact analytic cases
+uv run python tools/api_probe.py                     # IPC API probe vs live KiCad
-.venv\Scripts\python.exe tools\api_probe.py          # IPC API probe vs live KiCad
+uv run python -m fill_resistance.board_io dump.json [NET]  # extract only
-.venv\Scripts\python.exe -m fill_resistance.board_io dump.json [NET]  # extract only
 ```
 
 ## Packaging / publishing
@@ -188,7 +312,10 @@ filled in. To publish: upload the zip to a release, set `download_url`
 (and the `homepage` resource in `metadata.json`), then submit the
 registry copy as `packages/th.co.b4l.fill-resistance/metadata.json` in a
 merge request to <https://gitlab.com/kicad/addons/metadata>. Icons are
-regenerated with `python tools/gen_icons.py`.
+regenerated with `python tools/gen_icons.py`; the README figures in
+`docs/img/` with `uv run python tools/gen_readme_figs.py`
+(real solver output on small synthetic boards, plus the hand-drawn
+hole cross-section).
 
 ## License
 
@@ -197,7 +324,9 @@ GPL-3.0-or-later — see [LICENSE](LICENSE).
 ## Troubleshooting
 
 - **No toolbar button**: venv still building (wait), or build failed →
-  *Recreate Plugin Environment*; check the interpreter path (setup 2).
+  *Recreate Plugin Environment* (right-click the plugin's row in
+  Preferences → *PCB Editor → Action Plugins*); check the interpreter
+  path (setup 2).
 - **"Could not connect to KiCad's IPC API"**: API server not enabled, or
   KiCad not running (no headless mode in KiCad 10).
 - **"KiCad is busy"**: a modal dialog is open in KiCad — close it, rerun.
@@ -205,3 +334,23 @@ GPL-3.0-or-later — see [LICENSE](LICENSE).
   best-effort); PNGs are always saved regardless.
 - **Result seems too low/high**: remember the model is fills + barrels
   only, with ideal contacts; measure electrode-to-electrode.
+
+## LLM disclaimer
+
+This plugin was developed with an LLM: Anthropic's **Claude** (Claude
+Code, model Claude Fable 5). The physics model, solver, tests, tooling
+and this documentation (including the figures; all but the hand-drawn
+hole cross-section are generated by the solver itself) were written by
+the model, feature by feature, under human direction and review
+(janik / B4L); most commits carry a `Co-Authored-By: Claude` trailer.
+
+What keeps this honest: the test suite pins the numerics to exact
+analytic references (strip and annulus resistances, the acosh spreading
+resistance of two circular contacts, skin-effect limits, power-balance
+identities) and to convergence/regression checks; run it with
+`uv run pytest`. Real boards were measured against a UT3513+ micro-ohm
+meter (see *Measured vs. computed* above). Nevertheless, an LLM wrote
+this: read *Model & limits*
+critically, treat surprising numbers with the usual engineering
+suspicion, and cross-check against a hand estimate before trusting a
+result with hardware. Bug reports are very welcome.

deploy.ps1

@@ -30,7 +30,8 @@ if ($Mode -eq 'Junction') {
     New-Item -ItemType Junction -Path $dst -Target $src | Out-Null
     Write-Host "junction created: $dst -> $src"
 } else {
-    $exclude = @('.venv', '.git', 'tests', 'tools', '__pycache__', '.pytest_cache')
+    $exclude = @('.venv', '.git', 'tests', 'tools', '__pycache__', '.pytest_cache',
+                 'pyproject.toml', 'uv.lock')
     New-Item -ItemType Directory -Force $dst | Out-Null
     Get-ChildItem $src -Force | Where-Object { $exclude -notcontains $_.Name } |
         ForEach-Object { Copy-Item $_.FullName -Destination $dst -Recurse -Force }

docs/RELEASING.md

@@ -0,0 +1,75 @@
+# Releasing a new version
+
+The PCM addon zip is built by CI (`.gitea/workflows/build-pcm.yml`).
+Every push to `main` builds it as a downloadable artifact; pushing a
+`v<version>` tag additionally creates a Gitea release with the zip
+attached. The release job checks that the tag matches `metadata.json`
+and that the release notes exist, and fails on either mismatch.
+
+## Steps
+
+1. **Bump the version** in `metadata.json` — the single entry in
+   `versions` (plain `MAJOR.MINOR.PATCH`, no `v` prefix; the PCM schema
+   rejects anything else):
+
+   ```json
+   "versions": [
+       {
+           "version": "1.0.2",
+           "status": "stable",
+           "kicad_version": "10.0",
+           "runtime": "ipc"
+       }
+   ]
+   ```
+
+2. **Write the release notes** at `docs/release-notes/v<version>.md`.
+   This file becomes the release description verbatim; the job fails if
+   it is missing or empty (the release action publishes empty notes
+   rather than falling back to the tag message, which is how v1.2.0
+   shipped with a blank description). Say what changed for a user of
+   the previous version — in particular, whether results move for an
+   unchanged board.
+
+3. **Commit, tag, push** (tag = `v` + the manifest version):
+
+   ```powershell
+   git add metadata.json docs/release-notes/v1.0.2.md
+   git commit -m "Release 1.0.2"
+   git tag v1.0.2
+   git push
+   git push origin v1.0.2
+   ```
+
+4. **Verify**: the Actions run for the tag builds
+   `th.co.b4l.fill-resistance_<version>.zip` and publishes it at
+   <https://git.b4l.co.th/B4L/kicad-zone-resistance/releases>, together
+   with `metadata-registry.json`. The zip installs directly via
+   Plugin and Content Manager → *Install from File*.
+
+## Publishing to the official KiCad registry (optional)
+
+The attached `metadata-registry.json` already carries the release
+`download_url`, `download_sha256` and sizes. Submit it as
+`packages/th.co.b4l.fill-resistance/metadata.json` in a merge request
+to <https://gitlab.com/kicad/addons/metadata>. The registry keeps every
+published version: append the new entry to the `versions` array of the
+registry copy instead of replacing the previous one (the repo's own
+`metadata.json` only ever holds the current version —
+`tools/build_package.py` reads `versions[0]`).
+
+## Local build (no CI)
+
+```powershell
+python tools/build_package.py    # writes dist/<identifier>_<version>.zip
+```
+
+Pure stdlib — no venv needed. `dist/` is gitignored.
+
+## CI prerequisites (one-time, server side)
+
+- Actions enabled for the repo (Settings → Actions unit).
+- A runner registered with the `ubuntu-latest` label; the default
+  act_runner image works — the build needs only Python 3.
+- Workflow actions are pinned to commit SHAs; when bumping them, update
+  the SHA and the trailing version comment together.

docs/release-notes/v1.2.1.md

@@ -0,0 +1,36 @@
+Bug-fix release. Results are unchanged from 1.2.0 for a board that
+solves cleanly; the fixes are in the in-KiCad overlay push, pad copper
+selection and error reporting.
+
+Note for anyone coming from 1.1.0 or earlier: 1.2.0 changed the physics
+model (exact SMD and THT pad copper, populated THT holes conducting as
+their solder plug and lead, slotted holes as true stadiums) and fixed an
+adaptive barrel-refinement bug that could make via-field results read up
+to ~13% low. Numbers for an unchanged board differ from 1.1.0 - re-run
+any board you track across versions.
+
+Fixed:
+- Overlay push: a locked reference image silently survived removal and a
+  new one was stacked on top of it. KiCad reports the failure per item
+  while the overall request still reads OK; it is now checked, and the
+  layer is reported and skipped instead.
+- Overlay push: a run covering fewer layers than the previous one left
+  the earlier solve's heatmap on the unused slots, where it read as
+  current. Those slots are now cleared.
+- Overlay push: the whole push is one commit, so a single undo reverts
+  it rather than just the last layer.
+- Through-hole pad copper was always read from F.Cu even when the joint
+  protrudes on B.Cu, mis-sizing the modelled solder coat for pads sized
+  differently per copper layer. The solder side is now probed first.
+- A failure before the output directory existed - a broken plugin
+  Python environment, typically - reported nothing at all on screen.
+  The error figure now falls back to the temp directory.
+- Pads sitting on no single copper layer are noted rather than silently
+  skipped, and the frequency field keeps its specific rejection reason
+  ("1,500" is a thousands separator, "-5" is negative) as the other
+  numeric fields already did.
+
+The KiCad overlay push remains experimental and opt-in (off by default).
+It writes reference images to User.9-User.12 and replaces what is on
+those layers.
+

fill_resistance/adaptive.py

@@ -94,6 +94,7 @@ def run_solve_adaptive(problem: Problem, stack: RasterStack,
 
     # --- leaves per layer -------------------------------------------------
     t0 = time.perf_counter()
+    links, dead_barrels = sv._barrel_links(stack, problem)
     keep = e1 | e2
     if stack.chain is not None:
         keep |= stack.chain
@@ -101,6 +102,13 @@ def run_solve_adaptive(problem: Problem, stack: RasterStack,
         keep |= stack.buildup
     if stack.thick_scale is not None:
         keep |= stack.thick_scale != 1.0
+    # pin every barrel attachment cell fine: a point-like barrel
+    # injection into a coarse leaf makes the whole leaf equipotential
+    # and deletes the local spreading resistance (via fields read up
+    # to ~13% low otherwise); the guard ring then grades around it
+    for _vi, la, ia_, ja_, lb, ib_, jb_, _r in links:
+        keep[la, ia_, ja_] = True
+        keep[lb, ib_, jb_] = True
     mb = _max_block(stack.h_nm)
     grids = [quadtree.build_leaves(stack.masks[li], keep_fine=keep[li],
                                    max_block=mb,
@@ -181,7 +189,6 @@ def run_solve_adaptive(problem: Problem, stack: RasterStack,
             xx.append(np.full(k, -1, dtype=np.int8))
             ee.append(np.full(k, -1, dtype=np.int16))
 
-    links, dead_barrels = sv._barrel_links(stack, problem)
     for vi, la, ia_, ja_, lb, ib_, jb_, r_dc in links:
         na = offs[la] + grids[la].id_grid[ia_, ja_]
         nb = offs[lb] + grids[lb].id_grid[ib_, jb_]

fill_resistance/board_io.py

@@ -6,12 +6,13 @@ KiCad to extract without the dialog (all layers of the net, defaults).
 """
 from __future__ import annotations
 
+import math
 from dataclasses import dataclass, field
 from pathlib import Path
 
 from kipy import KiCad
 from kipy.board import Board
-from kipy.board_types import ArcTrack, BoardRectangle, Pad
+from kipy.board_types import ArcTrack, BoardRectangle, Pad, Via
 from kipy.proto.board.board_pb2 import BoardStackupLayerType
 from kipy.proto.board.board_types_pb2 import ZoneType
 from kipy.util.board_layer import (canonical_name, is_copper_layer,
@@ -22,7 +23,9 @@ import numpy as np
 from . import config
 from .errors import ApiVersionError, CandidateError, SelectionError
 from .geometry import (Electrode, LayerFill, Polygon, Problem, Rect,
-                       SurfaceBuildup, TrackSeg, ViaLink, linearize_ring)
+                       SurfaceBuildup, TrackSeg, ViaLink,
+                       contact_solder_buildups, linearize_ring,
+                       tht_joint_buildups)
 
 MASK_TO_COPPER = {"F.Mask": "F.Cu", "B.Mask": "B.Cu"}
 
@@ -137,11 +140,33 @@ def _convert_poly(poly_with_holes) -> Polygon:
                    holes=[ring(h) for h in poly_with_holes.holes])
 
 
-def _pad_drill_nm(pad_or_via) -> int:
+def _drill_info(pad_or_via) -> tuple[int, int, int]:
+    """(width_nm, slot_dx_nm, slot_dy_nm) of a padstack drill. Round
+    holes: (diameter, 0, 0). Slotted (oblong) holes: width is the
+    NARROW dimension, (slot_dx, slot_dy) the board-frame offset from
+    the drill center to each end-cap center of the slot. The slot
+    follows the pad rotation (KiCad rotates CCW with y down:
+    x' = x cos + y sin, y' = y cos - x sin)."""
     try:
-        return int(pad_or_via.padstack.drill.diameter.x)
+        d = pad_or_via.padstack.drill.diameter
+        dx, dy = int(d.x), int(d.y)
     except Exception:
-        return 0
+        return 0, 0, 0
+    if dx <= 0 or dy <= 0 or dx == dy:
+        return max(dx, 0), 0, 0
+    half = (max(dx, dy) - min(dx, dy)) / 2.0
+    try:
+        th = math.radians(pad_or_via.padstack.angle.degrees)
+    except Exception:
+        th = 0.0
+    ux, uy = (1.0, 0.0) if dx > dy else (0.0, 1.0)
+    return (min(dx, dy),
+            int(round(half * (ux * math.cos(th) + uy * math.sin(th)))),
+            int(round(half * (uy * math.cos(th) - ux * math.sin(th)))))
+
+
+def _pad_drill_nm(pad_or_via) -> int:
+    return _drill_info(pad_or_via)[0]
 
 
 def _pad_default_contact(pad: Pad) -> str:
@@ -157,11 +182,18 @@ def _pad_default_contact(pad: Pad) -> str:
     return "all"
 
 
-def _pad_polygons(board: Board, pad: Pad, contact: str) -> list[Polygon] | None:
+def _pad_polygons(board: Board, pad: Pad, contact: str,
+                  prefer: str | None = None) -> list[Polygon] | None:
+    """Exact pad copper. The first probed layer that has a shape wins, so
+    `prefer` (the solder side of a THT joint) must be tried before the
+    F.Cu/B.Cu fallback: KiCad allows a different pad size per copper
+    layer, and the solder coat is sized from this shape."""
     layer_ids = []
-    if contact != "all":
+    for name in (contact if contact != "all" else None, prefer):
+        if not name:
+            continue
         try:
-            layer_ids.append(layer_from_canonical_name(contact))
+            layer_ids.append(layer_from_canonical_name(name))
         except Exception:
             pass
     for name in ("F.Cu", "B.Cu"):
@@ -179,7 +211,42 @@ def _pad_polygons(board: Board, pad: Pad, contact: str) -> list[Polygon] | None:
     return None
 
 
-def _to_electrode(board: Board, item) -> Electrode:
+def _footprint_pad_map(footprints) -> dict:
+    """(x, y, number) -> owning FootprintInstance. Footprint pads are
+    stored with absolute positions, so the lookup is exact."""
+    out = {}
+    for fp in footprints or []:
+        try:
+            for fpad in fp.definition.pads:
+                out[(fpad.position.x, fpad.position.y, fpad.number)] = fp
+        except Exception:
+            continue
+    return out
+
+
+def _pad_owner(pad: Pad, pad_map: dict):
+    return pad_map.get((pad.position.x, pad.position.y, pad.number))
+
+
+def _tht_protrusion_side(pad: Pad, pad_map: dict, quiet: bool = False) -> str:
+    """Outer layer where the clipped THT lead protrudes (tent + solder
+    cone): the side OPPOSITE the component. Unknown owner -> assume the
+    component sits on F.Cu (lead tents on B.Cu)."""
+    fp = _pad_owner(pad, pad_map)
+    if fp is not None:
+        try:
+            side = canonical_name(fp.layer)
+            return "F.Cu" if side == "B.Cu" else "B.Cu"
+        except Exception:
+            pass
+    if not quiet:
+        print(f"note: no footprint found for pad {pad.number} - assuming "
+              f"its lead protrudes on B.Cu")
+    return "B.Cu"
+
+
+def _to_electrode(board: Board, item, stackup: StackupInfo | None = None,
+                  pad_map: dict | None = None) -> Electrode:
     if isinstance(item, BoardRectangle):
         tl, br = item.top_left, item.bottom_right
         rect = Rect.normalized(tl.x, tl.y, br.x, br.y,
@@ -188,6 +255,22 @@ def _to_electrode(board: Board, item) -> Electrode:
         cy = (rect.y0 + rect.y1) / 2e6
         return Electrode(rect=rect, contact="all",
                          label=f"rect({cx:.1f},{cy:.1f})")
+    if isinstance(item, Via):
+        via: Via = item
+        x, y = via.position.x, via.position.y
+        drill = int(via.drill_diameter or 0) or _pad_drill_nm(via)
+        if drill <= 0:
+            raise SelectionError(
+                f"Selected via at ({x / 1e6:.2f}, {y / 1e6:.2f}) mm has no "
+                f"drill diameter - cannot use it as a contact.")
+        pad_nm = _padstack_pad_nm(via)
+        r = max(pad_nm, drill) // 2
+        rect = Rect.normalized(x - r, y - r, x + r, y + r, "via")
+        return Electrode(
+            rect=rect, contact="all", label=f"via({x / 1e6:.1f},{y / 1e6:.1f})",
+            drill_nm=drill, pad_nm=pad_nm, center=(x, y),
+            barrel_z=(_padstack_span(via.padstack, stackup)
+                      if stackup is not None else None))
     # Pad
     pad: Pad = item
     contact = _pad_default_contact(pad)
@@ -197,37 +280,57 @@ def _to_electrode(board: Board, item) -> Electrode:
     if box is None:
         raise SelectionError(f"Could not get the bounding box of {label}.")
     rect = _box2_to_rect(box, "pad")
+    drill, slot_dx, slot_dy = _drill_info(pad)
+    prot = _tht_protrusion_side(pad, pad_map or {}) if drill > 0 else None
     return Electrode(rect=rect, contact=contact,
-                     polygons=_pad_polygons(board, pad, contact), label=label)
+                     polygons=_pad_polygons(board, pad, contact, prefer=prot),
+                     label=label,
+                     # through-hole pad: current enters at the soldered
+                     # barrel; the joint is solder-filled + pad-coated,
+                     # with a solder cone around the protruding lead
+                     drill_nm=drill, pad_nm=_padstack_pad_nm(pad),
+                     pad_min_nm=_padstack_pad_min_nm(pad),
+                     slot_dx_nm=slot_dx, slot_dy_nm=slot_dy,
+                     center=(pad.position.x, pad.position.y),
+                     solder=drill > 0, protrusion_side=prot)
 
 
-def _net_hint_of(pads: list[Pad]) -> str | None:
+def _net_hint_of(items: list) -> str | None:
-    for pad in pads:
+    for item in items:
-        if pad.net is not None:
+        if item.net is not None:
-            return pad.net.name
+            return item.net.name
     return None
 
 
-def get_electrodes(board: Board
+def get_electrodes(board: Board, stackup: StackupInfo | None = None
                    ) -> tuple[list[Electrode], list[Electrode], str | None]:
     """Terminals from the selection. Each terminal may have MULTIPLE parts
     (all merged into one externally-bonded contact):
 
     - rectangles on ELECTRODE_POS_LAYER -> V+ parts, on ELECTRODE_NEG_LAYER
-      -> V- parts; selected pads fill a side that has no rectangles;
+      -> V- parts; selected pads/vias fill a side that has no rectangles;
     - no marker rectangles selected: legacy mode, exactly 2 items
-      (rects/pads, any layer) -> one part each;
+      (rects/pads/vias, any layer) -> one part each;
     - empty selection: board-wide scan of both marker layers.
+
+    Selected vias and through-hole pads become BARREL contacts: current
+    enters at the drill-wall ring (the soldered lead/wire), not the pad
+    face. Draw a marker rectangle over the pad instead to model a probe
+    pressed onto the pad face.
     """
     pos_l = config.ELECTRODE_POS_LAYER
     neg_l = config.ELECTRODE_NEG_LAYER
     scheme = (f"Draw V+ rectangle(s) on {pos_l} and V- rectangle(s) on "
-              f"{neg_l} (axis-aligned), and/or select pads for a side "
+              f"{neg_l} (axis-aligned), and/or select pads/vias for a side "
               f"without rectangles.")
 
     selection = list(board.get_selection())
     rects = [s for s in selection if isinstance(s, BoardRectangle)]
-    pads = [s for s in selection if isinstance(s, Pad)]
+    pads = [s for s in selection if isinstance(s, (Pad, Via))]
+    # protrusion-side lookup needs the owning footprints (THT pads only)
+    pad_map = (_footprint_pad_map(board.get_footprints())
+               if any(isinstance(s, Pad) and _pad_drill_nm(s) > 0
+                      for s in pads) else {})
 
     if not selection:
         allr = [s for s in board.get_shapes() if isinstance(s, BoardRectangle)]
@@ -258,12 +361,13 @@ def get_electrodes(board: Board
         es2 = [_to_electrode(board, r) for r in neg]
         if pads and es1 and es2:
             raise SelectionError(
-                f"Cannot assign the {len(pads)} selected pad(s): both marker "
+                f"Cannot assign the {len(pads)} selected pad(s)/via(s): both "
-                f"layers already provide rectangles. Use pads only for a "
+                f"marker layers already provide rectangles. Use pads/vias "
-                f"side that has none."
+                f"only for a side that has none."
             )
         if pads:
-            pad_parts = [_to_electrode(board, p) for p in pads]
+            pad_parts = [_to_electrode(board, p, stackup, pad_map)
+                         for p in pads]
             if not es1:
                 es1 = pad_parts
             else:
@@ -277,12 +381,13 @@ def get_electrodes(board: Board
 
     items = rects + pads
     if len(items) == 2:
-        return ([_to_electrode(board, items[0])],
+        return ([_to_electrode(board, items[0], stackup, pad_map)],
-                [_to_electrode(board, items[1])], _net_hint_of(pads))
+                [_to_electrode(board, items[1], stackup, pad_map)],
+                _net_hint_of(pads))
     raise SelectionError(
         f"The selection has {len(rects)} rectangle(s) (none on the marker "
-        f"layers) and {len(pads)} pad(s); without marker layers exactly 2 "
+        f"layers) and {len(pads)} pad(s)/via(s); without marker layers "
-        f"contacts are needed.\n{scheme}"
+        f"exactly 2 contacts are needed.\n{scheme}"
     )
 
 
@@ -416,6 +521,18 @@ def _padstack_pad_nm(item) -> int:
         return 0
 
 
+def _padstack_pad_min_nm(item) -> int:
+    """Smallest dimension of the (largest) copper pad of a padstack; 0
+    if unknown. Bounds the lead-cone taper on oblong pads: the cone
+    stays within the inscribed circle."""
+    try:
+        sizes = [min(int(l.size.x), int(l.size.y))
+                 for l in item.padstack.copper_layers]
+        return max(sizes) if sizes else 0
+    except Exception:
+        return 0
+
+
 def _padstack_span(padstack, stackup: StackupInfo) -> tuple[int, int]:
     """(z_top, z_bot) of the barrel; falls back to the full stack."""
     try:
@@ -444,19 +561,209 @@ def gather_barrels(board: Board, net_name: str,
                                z_bot_nm=z_bot, kind="via",
                                pad_nm=_padstack_pad_nm(via)))
     if config.INCLUDE_TH_PADS:
-        for pad in board.get_pads():
+        net_pads = [pad for pad in board.get_pads()
-            if pad.net is None or pad.net.name != net_name:
+                    if pad.net is not None and pad.net.name == net_name
-                continue
+                    and _pad_drill_nm(pad) > 0]
-            drill = _pad_drill_nm(pad)
+        # populated (non-DNP) THT pads carry a soldered joint: filled
-            if drill <= 0:
+        # hole + coat + lead cone on the side opposite the component
-                continue
+        pad_map = (_footprint_pad_map(board.get_footprints())
-            barrels.append(ViaLink(x=pad.position.x, y=pad.position.y,
+                   if net_pads else {})
-                                   drill_nm=drill, z_top_nm=-1,
+        unknown = 0
-                                   z_bot_nm=stackup.z_bot_nm + 1, kind="pad",
+        for pad in net_pads:
-                                   pad_nm=_padstack_pad_nm(pad)))
+            fp = _pad_owner(pad, pad_map)
+            unknown += fp is None
+            populated = True
+            if fp is not None:
+                try:
+                    populated = not fp.attributes.do_not_populate
+                except Exception:
+                    pass
+            drill, slot_dx, slot_dy = _drill_info(pad)
+            barrels.append(ViaLink(
+                x=pad.position.x, y=pad.position.y,
+                drill_nm=drill, z_top_nm=-1,
+                z_bot_nm=stackup.z_bot_nm + 1, kind="pad",
+                pad_nm=_padstack_pad_nm(pad),
+                pad_min_nm=_padstack_pad_min_nm(pad),
+                slot_dx_nm=slot_dx, slot_dy_nm=slot_dy,
+                solder_filled=populated,
+                protrusion_side=(_tht_protrusion_side(pad, pad_map,
+                                                      quiet=True)
+                                 if populated else None)))
+        if unknown:
+            print(f"note: {unknown} THT pad(s) without an identifiable "
+                  f"footprint - assumed populated, leads on B.Cu")
     return barrels
 
 
+def gather_smd_pad_copper(board: Board, net_name: str
+                          ) -> dict[str, list[Polygon]]:
+    """layer name -> exact copper shape(s) of every SMD (undrilled) pad
+    on the net. Pads are junctions: traces and thermal-relief spokes
+    meet ON the pad copper, and without it the junction necks down to
+    the accidental overlap of the track ends - or is severed outright.
+    Dead-end pads (component terminals) become floating islands that
+    the solver's connectivity restriction drops. One API call per pad;
+    pads whose copper layer cannot be determined are skipped."""
+    shapes: dict[str, list[Polygon]] = {}
+    for pad in board.get_pads():
+        if pad.net is None or pad.net.name != net_name \
+                or _pad_drill_nm(pad) > 0:
+            continue
+        layer = _pad_default_contact(pad)      # SMD: its own copper layer
+        if layer == "all":
+            # zero or >1 copper layers (custom padstack): no single layer
+            # to stamp it on. Say so - a silent skip loses a real junction
+            print(f"note: pad {pad.number}@{net_name} sits on no single "
+                  f"copper layer - its pad copper is not modelled")
+            continue
+        polys = _pad_polygons(board, pad, layer)
+        if polys:
+            shapes.setdefault(layer, []).extend(polys)
+    return shapes
+
+
+def gather_tht_pad_copper(board: Board, net_name: str
+                          ) -> dict[tuple[int, int], list[Polygon]]:
+    """(x, y) -> exact copper shape(s) of every drilled (THT) pad on the
+    net. The annular-ring copper conducts on every layer the barrel
+    spans, so build_problem stamps these onto each included layer. One
+    API call per pad; the outer-layer shape stands in for the inner
+    rings (approximation - inner rings are usually the same or
+    smaller)."""
+    shapes: dict[tuple[int, int], list[Polygon]] = {}
+    for pad in board.get_pads():
+        if pad.net is None or pad.net.name != net_name \
+                or _pad_drill_nm(pad) <= 0:
+            continue
+        polys = _pad_polygons(board, pad, "all")
+        if polys:
+            shapes[(pad.position.x, pad.position.y)] = polys
+    return shapes
+
+
+# --- in-KiCad result overlays (EXPERIMENTAL) ---------------------------------
+
+# KiCad sizes reference images as pixels * (1 inch / PPI) * image_scale
+# and assumes 300 PPI for PNGs without a density chunk (BITMAP_BASE)
+OVERLAY_PIX_NM = 25.4e6 / 300
+
+
+def _create_reference_image(board: Board, ref) -> None:
+    """create_items with the per-item status surfaced (kipy <= 0.7.1
+    swallows it and returns an empty wrapper on failure)."""
+    from kipy.proto.common.commands.editor_commands_pb2 import (
+        CreateItems, CreateItemsResponse)
+    from kipy.util import pack_any
+
+    cmd = CreateItems()
+    cmd.header.document.CopyFrom(board._doc)
+    cmd.items.append(pack_any(ref.proto))
+    result = board._kicad.send(cmd, CreateItemsResponse).created_items[0]
+    if result.status.code != 1:                 # 1 = ISC_OK
+        raise RuntimeError(
+            f"KiCad rejected the image (status {result.status.code}) "
+            f"{result.status.error_message or ''} - is the layer enabled "
+            f"in Board Setup? (KiCad >= 10.0.1 required)")
+
+
+def remove_overlays(board: Board, layer) -> int:
+    """Remove every reference image on the given layer; returns count.
+
+    remove_items with the per-item status surfaced: kipy discards the
+    DeleteItemsResponse, and its own proto warns the overall status "may
+    return IRS_OK even if no items were deleted" - a locked image comes
+    back IDS_IMMUTABLE. Unchecked, the stale image survives and the new
+    one is stacked on top of it instead of replacing it."""
+    from kipy.proto.common.commands.editor_commands_pb2 import (
+        DeleteItems, DeleteItemsResponse, ItemDeletionStatus)
+
+    ours = [r for r in board.get_reference_images() if r.layer == layer]
+    if not ours:
+        return 0
+
+    cmd = DeleteItems()
+    cmd.header.document.CopyFrom(board._doc)
+    cmd.item_ids.extend([r.id for r in ours])
+    results = board._kicad.send(cmd, DeleteItemsResponse).deleted_items
+
+    stuck = [r for r in results
+             if r.status not in (ItemDeletionStatus.IDS_OK,
+                                 ItemDeletionStatus.IDS_NONEXISTENT)]
+    if stuck:
+        locked = sum(1 for r in stuck
+                     if r.status == ItemDeletionStatus.IDS_IMMUTABLE)
+        raise RuntimeError(
+            f"{len(stuck)} existing overlay image(s) could not be removed"
+            + (f" ({locked} locked)" if locked else "")
+            + " - unlock them in KiCad, or delete them by hand, then run "
+              "again (a new image would otherwise stack on top).")
+    return len(results)
+
+
+def push_result_overlays(board: Board, stack, result,
+                         lock: bool = False) -> None:
+    """EXPERIMENTAL: the solved |J| of every included copper layer as an
+    unlocked reference image on config.OVERLAY_LAYERS (stackup order,
+    top first; existing images there are replaced, and slots this run
+    does not write are cleared so no stale heatmap is left behind).
+    The whole push is one commit, so a single undo reverts it. Editor-
+    only - reference images never plot. Per-layer failures are reported
+    and skipped, never fatal to the run."""
+    from kipy.board_types import ReferenceImage
+    from kipy.geometry import Vector2
+
+    from .overlay import heatmap_png
+
+    names = stack.layer_names
+    pairs = list(zip(names, config.OVERLAY_LAYERS))
+    if len(names) > len(config.OVERLAY_LAYERS):
+        print(f"overlays: more copper layers than slots - "
+              f"{', '.join(names[len(config.OVERLAY_LAYERS):])} skipped")
+    ny, nx = stack.shape2d
+    w_nm, h_nm = nx * stack.h_nm, ny * stack.h_nm
+
+    commit = board.begin_commit() if hasattr(board, "begin_commit") else None
+    done = False
+    try:
+        # a narrower run than last time writes fewer slots; whatever the
+        # zip above left out still holds the previous solve's heatmap and
+        # would read as current, so clear it
+        for dest_name in config.OVERLAY_LAYERS[len(pairs):]:
+            try:
+                if remove_overlays(board, layer_from_canonical_name(dest_name)):
+                    print(f"overlay: cleared stale {dest_name}")
+            except Exception as e:
+                print(f"overlay: clearing stale {dest_name} failed: {e}")
+
+        for src, dest_name in pairs:
+            try:
+                dest = layer_from_canonical_name(dest_name)
+                png = heatmap_png(result.Jmag * 1e-6, names.index(src))
+                remove_overlays(board, dest)
+                ref = ReferenceImage()
+                ref.layer = dest
+                ref.position = Vector2.from_xy(round(stack.x0_nm + w_nm / 2),
+                                               round(stack.y0_nm + h_nm / 2))
+                ref.image_scale = w_nm / (nx * OVERLAY_PIX_NM)
+                ref.image_data = png
+                ref.locked = lock
+                _create_reference_image(board, ref)
+                print(f"overlay: |J| of {src} -> {dest_name} "
+                      f"({len(png) / 1024:.0f} kB)")
+            except Exception as e:
+                print(f"overlay: {src} -> {dest_name} failed: {e}")
+        if commit is not None:
+            board.push_commit(commit, "Fill Resistance |J| overlays")
+        done = True
+    finally:
+        if commit is not None and not done:
+            try:
+                board.drop_commit(commit)
+            except Exception:
+                pass
+
+
 # --- top level ----------------------------------------------------------------
 
 def build_problem(board: Board, net: str, layer_names: list[str],
@@ -465,7 +772,8 @@ def build_problem(board: Board, net: str, layer_names: list[str],
                   buildups: dict[str, list[Polygon]] | None = None,
                   extra_cu_um: float | None = None,
                   tracks: dict | None = None,
-                  vias_capped: bool | None = None) -> Problem:
+                  vias_capped: bool | None = None,
+                  cap_max_drill_mm: float | None = None) -> Problem:
     per_layer = fills.get(net, {})
     per_layer_tracks = (tracks or {}).get(net, {})
     layers = []
@@ -490,7 +798,32 @@ def build_problem(board: Board, net: str, layer_names: list[str],
             f"Net {net} has no fill on any of the selected layers "
             f"({', '.join(layer_names)})."
         )
-    vias = gather_barrels(board, net, stackup) if len(layers) > 1 else []
+    # barrels matter on a single layer too: via rings + drill mouths
+    # perforate the plane, THT joints locally stiffen it
+    vias = gather_barrels(board, net, stackup)
+    # THT pad copper is part of the conductor: stamp the exact pad
+    # shapes onto every included layer (the barrel spans the stack)
+    pad_shapes = (gather_tht_pad_copper(board, net)
+                  if any(v.kind == "pad" for v in vias) else {})
+    if pad_shapes:
+        extra = [poly for polys in pad_shapes.values() for poly in polys]
+        for layer in layers:
+            layer.polygons = list(layer.polygons) + extra
+        print(f"{len(pad_shapes)} THT pad shape(s) stamped on every "
+              f"included layer")
+    # SMD pad copper too: pads are the junctions where traces/spokes
+    # meet (also gives selected SMD-pad contacts their real copper)
+    smd_shapes = (gather_smd_pad_copper(board, net)
+                  if config.INCLUDE_SMD_PADS else {})
+    if smd_shapes:
+        n = 0
+        for layer in layers:
+            polys = smd_shapes.get(layer.layer_name, [])
+            if polys:
+                layer.polygons = list(layer.polygons) + polys
+                n += len(polys)
+        if n:
+            print(f"{n} SMD pad shape(s) stamped on their layers")
     included = {l.layer_name for l in layers}
     buildup_list = [
         SurfaceBuildup(layer_name=name, polygons=polys)
@@ -502,7 +835,7 @@ def build_problem(board: Board, net: str, layer_names: list[str],
           + (f", solder buildup on "
              f"{', '.join(b.layer_name for b in buildup_list)}"
              if buildup_list else ""))
-    return Problem(
+    problem = Problem(
         board_path=board.name or "",
         net_name=net,
         rho_ohm_m=config.RHO_CU_OHM_M,
@@ -522,7 +855,34 @@ def build_problem(board: Board, net: str, layer_names: list[str],
         vias_capped=(vias_capped if vias_capped is not None
                      else config.VIAS_CAPPED),
         cap_plating_nm=int(config.CAP_PLATING_UM * 1000),
+        cap_max_drill_nm=int((cap_max_drill_mm if cap_max_drill_mm is not None
+                              else config.CAP_MAX_DRILL_MM) * 1e6),
+        tht_protrusion_nm=int(config.THT_LEAD_PROTRUSION_MM * 1e6),
+        tht_lead_clearance_nm=int(config.THT_LEAD_CLEARANCE_MM * 1e6),
+        tht_lead_rho_ohm_m=config.THT_LEAD_RHO_OHM_M,
     )
+    solder_layers = contact_solder_buildups(problem)
+    if solder_layers:
+        sides = sorted({e.protrusion_side
+                        for e in problem.electrodes1 + problem.electrodes2
+                        if e.solder and e.protrusion_side})
+        cone = (f", {config.THT_LEAD_PROTRUSION_MM:g} mm lead + solder cone "
+                f"on {', '.join(sides)}"
+                if sides and problem.tht_protrusion_nm > 0 else "")
+        print(f"THT contact(s): solder-filled hole + "
+              f"{config.SOLDER_THICKNESS_UM:g} um average solder coat on the "
+              f"pad face ({', '.join(solder_layers)}){cone}")
+    tht_joint_buildups(problem, pad_shapes)
+    n_joint = sum(1 for v in problem.vias
+                  if v.kind == "pad" and v.solder_filled)
+    n_dnp = sum(1 for v in problem.vias
+                if v.kind == "pad" and not v.solder_filled)
+    if n_joint or n_dnp:
+        print(f"{n_joint} populated THT pad joint(s): lead + solder in the "
+              f"hole, coat + cone on the solder side"
+              + (f"; {n_dnp} DNP pad(s): open hole, plating-only"
+                 if n_dnp else ""))
+    return problem
 
 
 if __name__ == "__main__":
@@ -533,7 +893,7 @@ if __name__ == "__main__":
     out = Path(sys.argv[1]) if len(sys.argv) > 1 else Path("geometry_dump.json")
     _, board = connect()
     stackup = get_stackup_info(board)
-    es1, es2, net_hint = get_electrodes(board)
+    es1, es2, net_hint = get_electrodes(board, stackup)
     if any_zone_unfilled(board):
         refill(board)
     fills = gather_net_fills(board)

fill_resistance/config.py

@@ -28,7 +28,30 @@ VIAS_CAPPED = True              # filled + capped vias (dialog checkbox):
                                 # Ring/pad copper of vias is modeled either
                                 # way; THT-pad copper/drills are not.
 CAP_PLATING_UM = 15.0           # cap plating thickness (fab spec)
-INCLUDE_TH_PADS = True          # plated through-hole pads stitch layers too
+CAP_MAX_DRILL_MM = 0.5          # fab caps only small vias: drills above this
+                                # stay open even with VIAS_CAPPED
+                                # (dialog-settable)
+INCLUDE_SMD_PADS = True         # the net's SMD pad copper conducts too (exact
+                                # shapes on the pad's layer): pads are the
+                                # junctions where traces/spokes meet, and
+                                # selected pad contacts get their real copper.
+                                # Dead-end pads are dropped as floating islands
+INCLUDE_TH_PADS = True          # plated through-hole pads stitch layers too;
+                                # their holes are modeled solder-filled (a
+                                # soldered component lead), so the solder core
+                                # conducts in parallel with the plating
+THT_LEAD_PROTRUSION_MM = 1.5    # clipped THT lead protrusion on the side
+                                # opposite the component: a solder cone of
+                                # this height at the drill wall (tapering to
+                                # zero at the pad edge) wraps the lead of
+                                # every populated THT pad. 0 = no cones
+THT_LEAD_CLEARANCE_MM = 0.25    # hole diameter minus lead diameter (fab
+                                # rule): a lead cylinder of drill - this
+                                # conducts inside every solder-filled hole
+THT_LEAD_RHO_OHM_M = 1.68e-8    # lead material resistivity: copper leads/
+                                # wires; brass ~6.4e-8, phosphor bronze
+                                # ~1.1e-7, copper-clad steel higher - raise
+                                # this if your components use such leads
 SKIN_SIDES = 1                  # skin-effect field config: 1 = plane facing a
                                 # return plane (conservative), 2 = isolated foil
 
@@ -55,6 +78,20 @@ ELECTRODE_POS_LAYER = "User.1"  # rectangles on this layer mark V+ contact parts
 ELECTRODE_NEG_LAYER = "User.2"  # rectangles on this layer mark V- contact parts
 ALWAYS_REFILL = False           # refill zones even if KiCad says they are filled
 
+# --- In-KiCad result overlays (EXPERIMENTAL) ---
+PUSH_OVERLAYS = False           # after solving, push the per-layer |J|
+                                # heatmaps into the open board as unlocked
+                                # reference images (editor-only, never
+                                # plotted); dialog-toggleable
+OVERLAY_LAYERS = ("User.9", "User.10", "User.11", "User.12")
+                                # copper layers map here in stackup order
+                                # (top first); existing reference images on
+                                # these layers are REPLACED on every push;
+                                # each must be enabled in Board Setup
+OVERLAY_ALPHA = 255             # overlay opacity over copper (0-255);
+                                # translucency washes out over bright
+                                # copper - toggle the User layer instead
+
 # --- Adaptive grid ---
 ADAPTIVE_CELLS = True           # solve on a 2:1-balanced quadtree: fine at
                                 # copper boundaries/electrodes/features,

fill_resistance/dialog.py

@@ -37,7 +37,9 @@ class Selection:
     extra_cu_um: float = 0.0
     include_tracks: bool = True
     vias_capped: bool = True
-    adaptive: bool = False
+    cap_max_drill_mm: float = 0.5
+    adaptive: bool = True
+    push_overlays: bool = False   # EXPERIMENTAL in-KiCad |J| overlays
 
 
 class _Dialog(QDialog):
@@ -73,9 +75,14 @@ class _Dialog(QDialog):
         self.capped_check.setChecked(config.VIAS_CAPPED)
         form.addRow("Vias:", self.capped_check)
 
+        self.cap_drill_edit = QLineEdit(f"{config.CAP_MAX_DRILL_MM:g}")
+        self.cap_drill_edit.setEnabled(config.VIAS_CAPPED)
+        self.capped_check.toggled.connect(self.cap_drill_edit.setEnabled)
+        form.addRow("Capped up to drill [mm]:", self.cap_drill_edit)
+
         self.adaptive_check = QCheckBox(
             "adaptive cells (coarsen plane interiors; faster on large "
-            "boards, corrected to ≲0.1 % of the uniform grid)")
+            "boards, corrected to ≲0.03 % of the uniform grid)")
         self.adaptive_check.setChecked(config.ADAPTIVE_CELLS)
         form.addRow("Grid:", self.adaptive_check)
 
@@ -115,6 +122,14 @@ class _Dialog(QDialog):
         self.extracu_edit.setEnabled(bool(buildup_layers))
         form.addRow("Extra Cu in openings [µm]:", self.extracu_edit)
 
+        first, last = config.OVERLAY_LAYERS[0], config.OVERLAY_LAYERS[-1]
+        self.overlay_check = QCheckBox(
+            f"experimental: push per-layer |J| heatmaps into the board as "
+            f"reference images on {first}..{last} (replaces images there; "
+            f"layers must be enabled in Board Setup)")
+        self.overlay_check.setChecked(config.PUSH_OVERLAYS)
+        form.addRow("Overlays:", self.overlay_check)
+
         buttons = QDialogButtonBox(QDialogButtonBox.Ok | QDialogButtonBox.Cancel)
         buttons.accepted.connect(self._try_accept)
         buttons.rejected.connect(self.reject)
@@ -181,10 +196,13 @@ class _Dialog(QDialog):
             raise ValueError("Check at least one layer.")
 
         def number(edit: QLineEdit, name: str) -> float:
+            text = edit.text().strip()
             try:
-                return float(edit.text().strip().replace(",", "."))
+                return float(skin.normalize_decimal(text))
-            except ValueError:
+            except ValueError as exc:
-                raise ValueError(f"{name}: '{edit.text()}' is not a number.")
+                if "separator" in str(exc):
+                    raise ValueError(f"{name}: {exc}")
+                raise ValueError(f"{name}: '{text}' is not a number.")
 
         current = number(self.current_edit, "Test current")
         if current <= 0:
@@ -196,7 +214,11 @@ class _Dialog(QDialog):
                 raise ValueError("Cell size must be > 0 µm.")
         try:
             freq = skin.parse_frequency(self.freq_edit.text())
-        except ValueError:
+        except ValueError as exc:
+            # as in number(): keep parse_frequency's own explanation for
+            # the inputs it rejects deliberately, not just "unparseable"
+            if any(k in str(exc) for k in ("separator", "negative")):
+                raise ValueError(f"Frequency: {exc}")
             raise ValueError(
                 f"Frequency: cannot parse '{self.freq_edit.text()}' "
                 f"(examples: 0, 142k, 1.5M).")
@@ -205,6 +227,11 @@ class _Dialog(QDialog):
             extra_cu = number(self.extracu_edit, "Extra Cu")
             if extra_cu < 0:
                 raise ValueError("Extra Cu must be ≥ 0 µm.")
+        cap_max_drill = config.CAP_MAX_DRILL_MM
+        if self.capped_check.isChecked():
+            cap_max_drill = number(self.cap_drill_edit, "Capped up to drill")
+            if cap_max_drill <= 0:
+                raise ValueError("Capped-up-to drill must be > 0 mm.")
 
         def contact(box: QComboBox) -> str:
             t = box.currentText()
@@ -222,7 +249,9 @@ class _Dialog(QDialog):
                          extra_cu_um=extra_cu,
                          include_tracks=self.tracks_check.isChecked(),
                          vias_capped=self.capped_check.isChecked(),
-                         adaptive=self.adaptive_check.isChecked())
+                         cap_max_drill_mm=cap_max_drill,
+                         adaptive=self.adaptive_check.isChecked(),
+                         push_overlays=self.overlay_check.isChecked())
 
     def _try_accept(self) -> None:
         try:

fill_resistance/geometry.py

@@ -17,7 +17,7 @@ from pathlib import Path
 
 import numpy as np
 
-JSON_SCHEMA_VERSION = 5
+JSON_SCHEMA_VERSION = 6
 
 
 @dataclass(frozen=True)
@@ -93,16 +93,44 @@ class TrackSeg:
 
 @dataclass
 class Electrode:
-    """One PART of a current-injection terminal: a drawn rectangle or a
+    """One PART of a current-injection terminal: a drawn rectangle, a
-    selected pad. A terminal (V+ or V-) is a LIST of parts, all merged
+    selected pad, or a selected via. A terminal (V+ or V-) is a LIST of
-    into one equipotential contact (externally bonded). `polygons`
+    parts, all merged into one equipotential contact (externally
-    (board nm) is the exact copper shape when known (pads); None means
+    bonded). `polygons` (board nm) is the exact copper shape when known
-    the rectangle itself is the shape. `contact` = 'all' or a layer
+    (pads); None means the rectangle itself is the shape. `contact` =
-    name: which included layers this part touches."""
+    'all' or a layer name: which included layers this part touches.
+
+    drill_nm > 0 marks a BARREL contact (selected via or through-hole
+    pad): the current physically enters through the plated barrel (the
+    lead/wire soldered into the hole), so the contact cells are the
+    copper ring at the drill wall, not the whole pad face. `solder`
+    additionally models a soldered THT joint: the hole is filled with
+    solder and the pad face on the SOLDER side (protrusion_side,
+    opposite the component) carries an average-thickness solder coat
+    (Problem.solder_thickness_nm over `polygons`) plus the
+    protruding-lead cone."""
     rect: Rect                                # bounding box (labels/summary)
     contact: str = "all"
     polygons: list[Polygon] | None = None
     label: str = "rect"
+    drill_nm: int = 0                         # >0: barrel contact (slotted
+                                              # holes: the slot WIDTH)
+    pad_nm: int = 0                           # pad diameter (search bound;
+                                              # largest dimension if oblong)
+    pad_min_nm: int = 0                       # smallest pad dimension (cone
+                                              # taper bound); 0 = pad_nm
+    slot_dx_nm: int = 0                       # slotted (oblong) hole: offset
+    slot_dy_nm: int = 0                       # from `center` to each end-cap
+                                              # center of the slot, board
+                                              # frame; (0, 0) = round drill
+    center: tuple[int, int] | None = None     # drill center; None = rect center
+    barrel_z: tuple[int, int] | None = None   # (z_top, z_bot); None = full stack
+    solder: bool = False                      # soldered THT joint (see above)
+    protrusion_side: str | None = None        # outer layer where the clipped
+                                              # lead protrudes (opposite the
+                                              # component): a solder cone
+                                              # wraps it there, see
+                                              # Problem.tht_protrusion_nm
 
 
 @dataclass
@@ -111,21 +139,55 @@ class ViaLink:
     layers whose z lies within [z_top_nm, z_bot_nm]."""
     x: int
     y: int
-    drill_nm: int
+    drill_nm: int                             # slotted holes: the slot WIDTH
     z_top_nm: int
     z_bot_nm: int
     kind: str = "via"                         # "via" | "pad"
     pad_nm: int = 0                           # pad/annular diameter; 0 = unknown
+                                              # (oblong pads: LARGEST dimension,
+                                              # used as a search bound)
+    pad_min_nm: int = 0                       # smallest pad dimension (bounds
+                                              # the lead-cone taper on oblong
+                                              # pads); 0 = same as pad_nm
+    slot_dx_nm: int = 0                       # slotted (oblong) hole: offset
+    slot_dy_nm: int = 0                       # from (x, y) to each end-cap
+                                              # center of the slot, board
+                                              # frame; (0, 0) = round drill
+    solder_filled: bool = False               # populated THT pad: the hole
+                                              # holds lead + solder (in parallel
+                                              # with the plating); False for
+                                              # vias and DNP footprints
+    protrusion_side: str | None = None        # populated THT pad: outer layer
+                                              # where the clipped lead tents
+                                              # (solder cone), opposite the
+                                              # component side
 
     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:
+                          plating_nm: int,
+                          solder_rho_ohm_m: float | None = None,
+                          lead_nm: float = 0,
+                          lead_rho_ohm_m: float | None = None) -> float:
         """Barrel segment resistance over length_nm: thin-wall annulus of
-        plating around the drill."""
+        plating around the drill (slotted holes: thin wall around the
-        area_m2 = math.pi * (self.drill_nm * 1e-9) * (plating_nm * 1e-9)
+        stadium-shaped slot). With solder_rho_ohm_m the hole holds a
-        return rho_ohm_m * (length_nm * 1e-9) / area_m2
+        soldered THT joint: the component lead (a cylinder of lead_nm
+        diameter, resistivity lead_rho_ohm_m) and the solder filling the
+        remaining bore conduct in parallel with the plating."""
+        ext = 2.0 * math.hypot(self.slot_dx_nm, self.slot_dy_nm) * 1e-9
+        wall = math.pi * (self.drill_nm * 1e-9) + 2.0 * ext
+        ga = wall * (plating_nm * 1e-9) / rho_ohm_m
+                                              # conductance-area [m^2/ohm-m]
+        if solder_rho_ohm_m is not None:
+            r_core = max(self.drill_nm / 2.0 - plating_nm, 0.0) * 1e-9
+            r_lead = min(lead_nm * 1e-9 / 2.0, r_core)
+            if lead_rho_ohm_m is not None and r_lead > 0:
+                ga += math.pi * r_lead * r_lead / lead_rho_ohm_m
+            ga += (math.pi * r_core * r_core + 2.0 * r_core * ext
+                   - math.pi * r_lead * r_lead) / solder_rho_ohm_m
+        return (length_nm * 1e-9) / ga
 
 
 @dataclass
@@ -147,6 +209,23 @@ class Problem:
     vias_capped: bool = True                  # filled+capped vias: thin cap
     cap_plating_nm: int = 15_000              # over outer-layer mouths;
                                               # False = open mouths
+    cap_max_drill_nm: int = 500_000           # fab caps only small vias:
+                                              # drills above this stay open
+                                              # even with vias_capped
+    tht_protrusion_nm: int = 1_500_000        # clipped THT lead protrusion:
+                                              # a solder cone of this height
+                                              # at the drill wall (tapering
+                                              # to zero at the pad edge)
+                                              # wraps the lead on each solder
+                                              # contact's protrusion_side;
+                                              # 0 disables the cones
+    tht_lead_clearance_nm: int = 250_000      # hole minus lead diameter (fab
+                                              # rule): the lead cylinder of
+                                              # drill - this conducts inside
+                                              # every solder-filled hole
+    tht_lead_rho_ohm_m: float = 1.68e-8       # lead material resistivity
+                                              # (copper; brass ~6.4e-8,
+                                              # copper-clad steel higher)
 
     @property
     def layer_names(self) -> list[str]:
@@ -173,6 +252,101 @@ class Problem:
         return int(x.min()), int(y.min()), int(x.max()), int(y.max())
 
 
+def contact_solder_buildups(problem: Problem) -> list[str]:
+    """Soldered THT-joint contacts: the pad face on the SOLDER side (the
+    protrusion side, opposite the component - the component-side face
+    stays bare) is covered in solder of average thickness
+    solder_thickness_nm. Adds one SurfaceBuildup there for every
+    `solder` electrode's pad shape (the buildup machinery intersects
+    with actual copper at raster time). Returns the affected layer
+    names. Called once when the problem is built."""
+    included = {l.layer_name for l in problem.layers}
+    touched = []
+    for e in problem.electrodes1 + problem.electrodes2:
+        if not e.solder or not e.polygons \
+                or e.protrusion_side not in included:
+            continue
+        problem.buildups.append(
+            SurfaceBuildup(layer_name=e.protrusion_side,
+                           polygons=list(e.polygons)))
+        touched.append(e.protrusion_side)
+    return sorted(set(touched))
+
+
+def slot_distance(xg, yg, dx_nm: int, dy_nm: int):
+    """Distance from points (xg, yg) (numpy-broadcastable, coordinates
+    RELATIVE to the hole center) to a slotted hole's axis - the segment
+    (-dx, -dy)..(+dx, +dy) between the end-cap centers. The slot wall
+    sits at distance width/2. Round drills (dx = dy = 0) reduce to the
+    plain radius, so callers need no special case."""
+    if dx_nm == 0 and dy_nm == 0:
+        return np.hypot(xg, yg)
+    l2 = float(dx_nm) * dx_nm + float(dy_nm) * dy_nm
+    t = np.clip((xg * dx_nm + yg * dy_nm) / l2, -1.0, 1.0)
+    return np.hypot(xg - t * dx_nm, yg - t * dy_nm)
+
+
+def _disc_polygon(x_nm: float, y_nm: float, r_nm: float,
+                  n: int = 32) -> Polygon:
+    th = np.linspace(0.0, 2.0 * math.pi, n, endpoint=False)
+    return Polygon(outline=np.round(np.stack(
+        [x_nm + r_nm * np.cos(th), y_nm + r_nm * np.sin(th)],
+        axis=1)).astype(np.int64))
+
+
+def _capsule_polygon(x_nm: float, y_nm: float, dx_nm: float, dy_nm: float,
+                     r_nm: float, n: int = 16) -> Polygon:
+    """Stadium: two half-circle caps of radius r_nm centered at
+    (x +- dx, y +- dy), joined by straight flanks."""
+    a0 = math.atan2(dy_nm, dx_nm)
+    th = np.linspace(-0.5 * math.pi, 0.5 * math.pi, n) + a0
+    cap1 = np.stack([x_nm + dx_nm + r_nm * np.cos(th),
+                     y_nm + dy_nm + r_nm * np.sin(th)], axis=1)
+    cap2 = np.stack([x_nm - dx_nm + r_nm * np.cos(th + math.pi),
+                     y_nm - dy_nm + r_nm * np.sin(th + math.pi)], axis=1)
+    return Polygon(outline=np.round(np.vstack([cap1, cap2])).astype(np.int64))
+
+
+def tht_joint_buildups(problem: Problem,
+                       shapes: dict | None = None) -> list[str]:
+    """Solder coat of the net's populated STITCHING through-hole pads
+    (ViaLink kind 'pad' with solder_filled), on the pad's SOLDER side
+    (the protrusion side, opposite the component; the component-side
+    face stays bare). `shapes` maps (x, y) to the exact pad polygons
+    (fetched from KiCad); pads without one fall back to a pad-diameter
+    disc. Contact pads are skipped: contact_solder_buildups already
+    coats them with the exact pad shape. Returns the affected layer
+    names."""
+    included = {l.layer_name for l in problem.layers}
+    contacts = {e.center for e in problem.electrodes1 + problem.electrodes2
+                if e.drill_nm > 0 and e.center is not None}
+    touched = []
+    for v in problem.vias:
+        if v.kind != "pad" or not v.solder_filled \
+                or (v.x, v.y) in contacts \
+                or v.protrusion_side not in included:
+            continue
+        polys = (shapes or {}).get((v.x, v.y))
+        if polys is None:
+            if v.pad_nm <= v.drill_nm:
+                continue
+            # oblong pads: never coat past the pad - a capsule along the
+            # slot axis, or the inscribed disc when the axis is unknown
+            w = v.pad_min_nm or v.pad_nm
+            hl = math.hypot(v.slot_dx_nm, v.slot_dy_nm)
+            if hl > 0.0 and v.pad_nm > w:
+                s = (v.pad_nm - w) / 2.0 / hl
+                polys = [_capsule_polygon(v.x, v.y, v.slot_dx_nm * s,
+                                          v.slot_dy_nm * s, w / 2.0)]
+            else:
+                polys = [_disc_polygon(v.x, v.y, w / 2.0)]
+        problem.buildups.append(
+            SurfaceBuildup(layer_name=v.protrusion_side,
+                           polygons=list(polys)))
+        touched.append(v.protrusion_side)
+    return sorted(set(touched))
+
+
 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."""
@@ -321,6 +495,15 @@ def _electrode_to_json(e: Electrode) -> dict:
         "label": e.label,
         "polygons": (None if e.polygons is None
                      else [_poly_to_json(poly) for poly in e.polygons]),
+        "drill_nm": e.drill_nm,
+        "pad_nm": e.pad_nm,
+        "pad_min_nm": e.pad_min_nm,
+        "slot_dx_nm": e.slot_dx_nm,
+        "slot_dy_nm": e.slot_dy_nm,
+        "center": (None if e.center is None else list(e.center)),
+        "barrel_z": (None if e.barrel_z is None else list(e.barrel_z)),
+        "solder": e.solder,
+        "protrusion_side": e.protrusion_side,
     }
 
 
@@ -331,6 +514,17 @@ def _electrode_from_json(d: dict) -> Electrode:
         label=d.get("label", "rect"),
         polygons=(None if d.get("polygons") is None
                   else [_poly_from_json(pd) for pd in d["polygons"]]),
+        drill_nm=int(d.get("drill_nm", 0)),
+        pad_nm=int(d.get("pad_nm", 0)),
+        pad_min_nm=int(d.get("pad_min_nm", 0)),
+        slot_dx_nm=int(d.get("slot_dx_nm", 0)),
+        slot_dy_nm=int(d.get("slot_dy_nm", 0)),
+        center=(None if d.get("center") is None
+                else (int(d["center"][0]), int(d["center"][1]))),
+        barrel_z=(None if d.get("barrel_z") is None
+                  else (int(d["barrel_z"][0]), int(d["barrel_z"][1]))),
+        solder=bool(d.get("solder", False)),
+        protrusion_side=d.get("protrusion_side"),
     )
 
 
@@ -369,6 +563,10 @@ def problem_to_json(p: Problem) -> dict:
         "extra_cu_nm": p.extra_cu_nm,
         "vias_capped": p.vias_capped,
         "cap_plating_nm": p.cap_plating_nm,
+        "cap_max_drill_nm": p.cap_max_drill_nm,
+        "tht_protrusion_nm": p.tht_protrusion_nm,
+        "tht_lead_clearance_nm": p.tht_lead_clearance_nm,
+        "tht_lead_rho_ohm_m": p.tht_lead_rho_ohm_m,
     }
 
 
@@ -415,7 +613,14 @@ def problem_from_json(d: dict) -> Problem:
             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)))
+                    pad_nm=int(vd.get("pad_nm", 0)),
+                    pad_min_nm=int(vd.get("pad_min_nm", 0)),
+                    slot_dx_nm=int(vd.get("slot_dx_nm", 0)),
+                    slot_dy_nm=int(vd.get("slot_dy_nm", 0)),
+                    # older dumps: every THT pad counted as solder-filled
+                    solder_filled=bool(vd.get(
+                        "solder_filled", vd.get("kind", "via") == "pad")),
+                    protrusion_side=vd.get("protrusion_side"))
             for vd in d["vias"]
         ],
         electrodes1=(
@@ -442,6 +647,10 @@ def problem_from_json(d: dict) -> Problem:
         ],
         vias_capped=bool(d.get("vias_capped", True)),
         cap_plating_nm=int(d.get("cap_plating_nm", 15_000)),
+        cap_max_drill_nm=int(d.get("cap_max_drill_nm", 500_000)),
+        tht_protrusion_nm=int(d.get("tht_protrusion_nm", 1_500_000)),
+        tht_lead_clearance_nm=int(d.get("tht_lead_clearance_nm", 250_000)),
+        tht_lead_rho_ohm_m=float(d.get("tht_lead_rho_ohm_m", 1.68e-8)),
     )
 
 

fill_resistance/main.py

@@ -18,9 +18,21 @@ from .errors import CandidateError, UserFacingError
 
 def _fail(message: str, outdir) -> None:
     print(f"ERROR: {message}")
-    from . import plots
+    try:
-    fig = plots.fig_error(message)
+        if outdir is None:
-    plots.save_and_show([(fig, "error")], outdir)
+            # A failure before the run has an output directory (a broken
+            # plugin environment throws on import) would otherwise save
+            # no PNG - and with no GUI toolkit, plots falls back to
+            # opening the saved PNGs, so the figure would never be shown
+            # either. Exactly the case the docstring promises to cover.
+            import tempfile
+            from pathlib import Path
+            outdir = Path(tempfile.gettempdir()) / "fill-resistance-error"
+        from . import plots
+        fig = plots.fig_error(message)
+        plots.save_and_show([(fig, "error")], outdir)
+    except Exception:                     # reporting must not mask the fault
+        traceback.print_exc()
     sys.exit(1)
 
 
@@ -33,7 +45,7 @@ def main() -> None:
         try:
             kicad, board = board_io.connect()
             stackup = board_io.get_stackup_info(board)
-            es1, es2, net_hint = board_io.get_electrodes(board)
+            es1, es2, net_hint = board_io.get_electrodes(board, stackup)
             if board_io.any_zone_unfilled(board) or config.ALWAYS_REFILL:
                 board_io.refill(board)
             fills = board_io.gather_net_fills(board)
@@ -95,15 +107,21 @@ def main() -> None:
                 buildups=(buildups if selection.include_buildup else None),
                 extra_cu_um=selection.extra_cu_um,
                 tracks=(tracks if selection.include_tracks else None),
-                vias_capped=selection.vias_capped)
+                vias_capped=selection.vias_capped,
+                cap_max_drill_mm=selection.cap_max_drill_mm)
             outdir = report.make_output_dir(board_io.board_dir(board))
         except ApiError as e:
             raise UserFacingError(f"KiCad API error: {e}")
 
         report.write_geometry_dump(outdir, problem)
+        overlay_cb = None
+        if selection.push_overlays:
+            def overlay_cb(stack, result):
+                board_io.push_result_overlays(board, stack, result)
         pipeline.run(problem, outdir, show=True, i_test=selection.current_a,
                      freq_hz=selection.freq_hz,
-                     contact_model=selection.contact_model)
+                     contact_model=selection.contact_model,
+                     overlay=overlay_cb)
     except UserFacingError as e:
         _fail(str(e), outdir)
     except Exception:

fill_resistance/overlay.py

@@ -0,0 +1,64 @@
+"""Rendering for the experimental in-KiCad result overlays: a solved
+field (|J|) as an RGBA PNG, one pixel per grid cell, transparent where
+there is no copper. The pushing side (ReferenceImages via the IPC API)
+lives in board_io; this module stays KiCad-free so it is testable
+headless.
+"""
+from __future__ import annotations
+
+import io
+
+import numpy as np
+
+from . import config
+
+# the colormap's near-black bottom must stay distinguishable from
+# KiCad's dark canvas (matplotlib figures sit on a light background
+# instead), so the log scale starts this far up the colormap
+FLOOR = 0.18
+
+
+def heatmap_png(data3: np.ndarray, li: int, alpha: int | None = None,
+                bleed: bool = True) -> bytes:
+    """One layer of a field (e.g. |J|, NaN = no copper) as opaque-over-
+    copper RGBA PNG bytes. Color scale matches the plugin's log figure
+    (global vmax across layers). `bleed` extends the edge color one
+    pixel outward at half opacity: the raster mask covers cells whose
+    CENTER is inside the copper, so without it the overlay stops half a
+    cell short of the outline KiCad draws."""
+    import matplotlib
+    from PIL import Image
+    from scipy import ndimage
+
+    if alpha is None:
+        alpha = config.OVERLAY_ALPHA
+    if not np.isfinite(data3).any():
+        raise ValueError("field is empty - nothing to overlay")
+    vmax = float(np.nanmax(data3))
+    if vmax <= 0:
+        raise ValueError("field is empty - nothing to overlay")
+    vmin = vmax / config.CURRENT_DYNAMIC_RANGE
+    d = np.clip(data3[li], vmin, vmax)
+    if config.LOG_CURRENT_SCALE:
+        u = (np.log(d) - np.log(vmin)) / (np.log(vmax) - np.log(vmin))
+    else:
+        u = d / vmax
+    u = FLOOR + (1.0 - FLOOR) * u
+    cmap = matplotlib.colormaps[config.CMAP_CURRENT]
+    rgba = (cmap(np.nan_to_num(u)) * 255).astype(np.uint8)
+    copper = ~np.isnan(data3[li])
+    rgba[..., 3] = np.where(copper, alpha, 0)
+
+    if bleed and copper.any() and not copper.all():
+        ring = ndimage.binary_dilation(
+            copper, structure=np.ones((3, 3), dtype=bool)) & ~copper
+        iy, ix = ndimage.distance_transform_edt(
+            ~copper, return_distances=False, return_indices=True)
+        rgba[ring, :3] = rgba[iy[ring], ix[ring], :3]
+        rgba[ring, 3] = alpha // 2
+
+    buf = io.BytesIO()
+    # no dpi metadata: KiCad assumes its 300 PPI default, which the
+    # pusher's scale computation relies on
+    Image.fromarray(rgba, "RGBA").save(buf, format="PNG")
+    return buf.getvalue()

fill_resistance/pipeline.py

@@ -12,7 +12,9 @@ from .solver import Result
 
 def run(problem: Problem, outdir: Path | None, show: bool = True,
         i_test: float | None = None, freq_hz: float = 0.0,
-        contact_model: str | None = None) -> Result:
+        contact_model: str | None = None, overlay=None) -> Result:
+    """overlay: optional callback(stack, result) run after the solve
+    (EXPERIMENTAL in-KiCad overlays); its failures are non-fatal."""
     if i_test is None:
         i_test = config.TEST_CURRENT_A
     if i_test <= 0:
@@ -43,6 +45,12 @@ def run(problem: Problem, outdir: Path | None, show: bool = True,
         report.write_summary(outdir, problem, stack, result)
     print(report.result_line(result, problem, stack))
 
+    if overlay is not None:
+        try:
+            overlay(stack, result)
+        except Exception as e:
+            print(f"overlay push failed: {e}")
+
     figs = [
         (plots.fig_raster(stack, e1, e2, problem, result), "1_raster_map"),
         (plots.fig_potential(result, stack, e1, e2, problem), "2_potential"),

fill_resistance/plots.py

@@ -50,7 +50,9 @@ _COPPER = "#c98b4e"
 _E1_COLOR = "#c8385a"
 _E2_COLOR = "#2f6fb0"
 _VIA_COLOR = "#2d6b45"
+_PAD_COLOR = "#5b4a8a"   # THT pad barrels (kind='pad'), violet-ink
 _SOLDER = "#9aa3ad"      # tin-gray: solder buildup areas
+_PLUG = "#6e7885"        # darker tin: solder-filled THT holes (lead + plug)
 _MESH = "#a56c33"        # darker copper: adaptive leaf boundaries
 _INK = "#3a3a3a"
 _GRID_INK = "#b8b4ae"
@@ -187,10 +189,14 @@ def _electrode_labels(ax, stack, e1_l, e2_l):
 
 
 def _via_markers(ax, problem, layer):
-    xs = [v.x * 1e-6 for v in problem.vias if v.spans(layer.z_nm)]
+    """One dot per barrel spanning the layer: vias green, THT pad
-    ys = [v.y * 1e-6 for v in problem.vias if v.spans(layer.z_nm)]
+    barrels violet (same joint markers, different physics)."""
-    if xs:
+    for kind, color in (("via", _VIA_COLOR), ("pad", _PAD_COLOR)):
-        ax.plot(xs, ys, ".", ms=2.5, color=_VIA_COLOR, alpha=0.7)
+        pts = [(v.x * 1e-6, v.y * 1e-6) for v in problem.vias
+               if v.kind == kind and v.spans(layer.z_nm)]
+        if pts:
+            xs, ys = zip(*pts)
+            ax.plot(xs, ys, ".", ms=2.5, color=color, alpha=0.7)
 
 
 def area_tag(sign: str, index: int) -> str:
@@ -219,8 +225,9 @@ def _injection_area_labels(ax, li, layer_name, problem, result):
 
 def fig_raster(stack, e1, e2, problem, result=None):
     cmap = ListedColormap([_BG, _COPPER, _E1_COLOR, _E2_COLOR, _SOLDER,
-                           _MESH])
+                           _MESH, _PLUG])
     has_buildup = stack.buildup is not None and stack.buildup.any()
+    has_plug = stack.plug is not None and stack.plug.any()
     has_mesh = stack.mesh is not None and stack.mesh.any()
 
     def paint(ax, li):
@@ -228,11 +235,13 @@ def fig_raster(stack, e1, e2, problem, result=None):
         codes[stack.masks[li]] = 1
         if has_buildup:
             codes[stack.buildup[li]] = 4
+        if has_plug:
+            codes[stack.plug[li]] = 6
         if has_mesh:
             codes[stack.mesh[li]] = 5
         codes[e1[li]] = 2
         codes[e2[li]] = 3
-        ax.imshow(codes, cmap=cmap, vmin=0, vmax=5, origin="upper",
+        ax.imshow(codes, cmap=cmap, vmin=0, vmax=6, origin="upper",
                   extent=stack.extent_mm(), interpolation="nearest")
         _via_markers(ax, problem, problem.layers[li])
         if result is not None and (result.part_currents1
@@ -243,8 +252,12 @@ def fig_raster(stack, e1, e2, problem, result=None):
             _electrode_labels(ax, stack, e1[li], e2[li])
 
     def finalize(fig, rows):
-        handles = [Patch(fc=_COPPER, label="copper"),
+        kinds = {v.kind for v in problem.vias}
-                   Patch(fc=_VIA_COLOR, label="vias")]
+        handles = [Patch(fc=_COPPER, label="copper")]
+        if "via" in kinds or not kinds:
+            handles.append(Patch(fc=_VIA_COLOR, label="vias"))
+        if "pad" in kinds:
+            handles.append(Patch(fc=_PAD_COLOR, label="THT pad barrels"))
         if has_mesh:
             handles.append(Patch(fc=_MESH,
                                  label="adaptive mesh (coarse leaves)"))
@@ -255,6 +268,9 @@ def fig_raster(stack, e1, e2, problem, result=None):
                       f"({problem.solder_thickness_nm / 1000:.0f} µm"
                       + (f" + {problem.extra_cu_nm / 1000:.0f} µm Cu"
                          if problem.extra_cu_nm else "") + ")"))
+        if has_plug:
+            handles.append(Patch(
+                fc=_PLUG, label="solder-filled THT hole (lead + solder)"))
         if result is not None and (result.part_currents1
                                    or result.part_currents2):
             entries = ([("+", _E1_COLOR, i, amps)
@@ -416,7 +432,7 @@ def fig_power(result, stack, e1, e2, problem):
 def fig_error(message: str):
     fig, ax = plt.subplots(figsize=(9, 4.5), layout="constrained")
     ax.axis("off")
-    ax.set_title("Fill Resistance — ERROR", color="#b02a2a",
+    ax.set_title("Fill Resistance - ERROR", color="#b02a2a",
                  fontsize=14, fontweight="bold", loc="left")
     wrapped = "\n".join(
         textwrap.fill(line, width=90) for line in message.splitlines()

fill_resistance/raster.py

@@ -23,7 +23,7 @@ from scipy import ndimage
 
 from . import config
 from .errors import ElectrodeError, GridSizeError
-from .geometry import Electrode, Problem, Rect
+from .geometry import Electrode, Problem, Rect, slot_distance
 
 # 4-connectivity: matches the in-plane 5-point stencil of the solver
 _STRUCT4 = ndimage.generate_binary_structure(2, 1)
@@ -46,7 +46,16 @@ class RasterStack:
     thick_scale: np.ndarray | None = None   # float (L, ny, nx): per-cell
                                             # copper-thickness factor (via
                                             # mouths: cap-thin or partially
-                                            # drilled cells); None = all 1
+                                            # drilled cells; folded-in cone
+                                            # and plug extras); None = all 1
+    t_extra_nm: np.ndarray | None = None    # float (L, ny, nx): additive
+                                            # conduction-equivalent copper
+                                            # (lead cones + hole plugs),
+                                            # folded into thick_scale at the
+                                            # end of rasterize_stack
+    plug: np.ndarray | None = None      # bool (L, ny, nx): solder-filled THT
+                                        # hole mouths (lead + solder plug,
+                                        # drawn on the raster map)
     mesh: np.ndarray | None = None      # bool (L, ny, nx): adaptive leaf
                                         # boundaries (drawn on the raster map)
 
@@ -232,9 +241,89 @@ def rasterize_stack(problem: Problem, h_nm: float) -> RasterStack:
                     _paint_ring(stack, hole, False, pmask)
                 stack.buildup[li] |= pmask
         stack.buildup &= stack.masks        # solder wets exposed copper only
+
+    _paint_lead_fillets(stack, problem)
+
+    if stack.t_extra_nm is not None:
+        # cones + plugs are ADDITIVE conduction-equivalent copper; fold
+        # them into the multiplicative per-cell scale once (multiplying
+        # per contribution would overstate cells carrying both)
+        if stack.thick_scale is None:
+            stack.thick_scale = np.ones(stack.masks.shape)
+        for li, layer in enumerate(problem.layers):
+            stack.thick_scale[li] *= np.where(
+                stack.masks[li],
+                1.0 + stack.t_extra_nm[li] / layer.thickness_nm, 1.0)
     return stack
 
 
+def _paint_lead_fillets(stack: RasterStack, problem: Problem) -> None:
+    """Protruding THT leads (barrel contacts AND the net's populated
+    stitching through-hole pads): the clipped lead sticks
+    tht_protrusion_nm out of the hole on the side opposite the
+    component, wrapped by a solder cone - full protrusion height at
+    the drill wall, tapering linearly to zero at the pad edge. Modeled
+    as extra conduction-equivalent copper (stack.t_extra_nm, ADDITIVE
+    with the hole plug, folded into thick_scale by rasterize_stack): the
+    tall solder column next to the wall pulls those cells to lead
+    potential (equivalent to extending the barrel wall vertically), the
+    taper carries the radial spreading. At f > 0 the factor multiplies
+    the skin-corrected sheet conductance, like the via mouths
+    (approximation)."""
+    H = problem.tht_protrusion_nm
+    if H <= 0:
+        return
+    ny, nx = stack.shape2d
+    h = stack.h_nm
+    index = {name: li for li, name in enumerate(stack.layer_names)}
+
+    # one cone per joint: contact electrodes first (exact data), then the
+    # net's populated stitching THT pads, skipping the contacts' barrels
+    jobs = []
+    seen = set()
+    for e in problem.electrodes1 + problem.electrodes2:
+        if e.drill_nm <= 0:
+            continue
+        if e.center is not None:
+            x, y = e.center
+        else:
+            x = (e.rect.x0 + e.rect.x1) / 2.0
+            y = (e.rect.y0 + e.rect.y1) / 2.0
+        seen.add((int(x), int(y)))
+        if e.solder and e.protrusion_side:
+            # oblong pads: taper from the (slot) wall to the inscribed
+            # dimension (conservative)
+            jobs.append((x, y, e.drill_nm, e.pad_min_nm or e.pad_nm,
+                         e.protrusion_side, e.slot_dx_nm, e.slot_dy_nm))
+    for v in problem.vias:
+        if v.kind == "pad" and v.solder_filled and v.protrusion_side \
+                and (v.x, v.y) not in seen:
+            jobs.append((v.x, v.y, v.drill_nm, v.pad_min_nm or v.pad_nm,
+                         v.protrusion_side, v.slot_dx_nm, v.slot_dy_nm))
+
+    for x, y, drill_nm, pad_nm, side, sdx, sdy in jobs:
+        li = index.get(side)
+        if li is None or pad_nm <= drill_nm:
+            continue
+        ra, rb = drill_nm / 2.0, pad_nm / 2.0
+        ex, ey = rb + abs(sdx), rb + abs(sdy)
+        j0 = max(0, math.floor((x - ex - stack.x0_nm) / h))
+        j1 = min(nx, math.floor((x + ex - stack.x0_nm) / h) + 1)
+        i0 = max(0, math.floor((y - ey - stack.y0_nm) / h))
+        i1 = min(ny, math.floor((y + ey - stack.y0_nm) / h) + 1)
+        if i0 >= i1 or j0 >= j1:
+            continue
+        xs = stack.x0_nm + (np.arange(j0, j1) + 0.5) * h - x
+        ys = stack.y0_nm + (np.arange(i0, i1) + 0.5) * h - y
+        r = slot_distance(xs[None, :], ys[:, None], sdx, sdy)
+        t_sn = H * np.clip((rb - r) / (rb - ra), 0.0, 1.0)
+        t_eq = t_sn * (problem.rho_ohm_m / problem.solder_rho_ohm_m)
+        if stack.t_extra_nm is None:
+            stack.t_extra_nm = np.zeros(stack.masks.shape)
+        m = stack.masks[li, i0:i1, j0:j1]
+        stack.t_extra_nm[li, i0:i1, j0:j1] += np.where(m, t_eq, 0.0)
+
+
 def _via_span(problem: Problem, via) -> list[int]:
     return [li for li, layer in enumerate(problem.layers)
             if via.spans(layer.z_nm)]
@@ -268,39 +357,79 @@ def _apply_via_mouths(stack: RasterStack, problem: Problem) -> None:
     """Drill-mouth treatment, area-weighted per cell (4x4 supersampling):
     capped vias carry a cap_plating-thin copper cap over the mouth on the
     OUTER layers, uncapped vias (and inner layers either way) get an open
-    hole. Fully swallowed cells leave the mask; partially covered cells
+    hole. The fab caps only small vias: drills above cap_max_drill_nm
-    keep a thickness-scaled sheet conductance via stack.thick_scale."""
+    stay open even with vias_capped. THT pad mouths: populated pads are
+    solder-filled - the mouth keeps its copper and additionally carries
+    the PLUG (the component lead plus the solder filling the bore) as
+    in-plane conduction-equivalent copper of the FULL hole depth on
+    EVERY spanned layer (the pin continues beyond both mouths, so each
+    layer sees the whole plug cross-section); the joint is then
+    side-symmetric except for the solder: the solder-side coat and cone
+    come on top, additively (see _paint_lead_fillets).
+    DNP pad holes are cut open on every layer. Fully swallowed cells
+    leave the mask; partially covered cells keep a thickness-scaled
+    sheet conductance via stack.thick_scale."""
     ny, nx = stack.shape2d
     h = stack.h_nm
     outer = {li for li, n in enumerate(stack.layer_names)
              if n in ("F.Cu", "B.Cu")}
     sub = (np.arange(4) + 0.5) / 4.0
     for via in problem.vias:
-        if via.kind != "via" or via.drill_nm <= 0:
+        if via.drill_nm <= 0:
             continue
+        plugged = via.kind == "pad" and via.solder_filled
         r = via.drill_nm / 2.0
-        j0 = max(0, math.floor((via.x - r - stack.x0_nm) / h))
+        ex, ey = r + abs(via.slot_dx_nm), r + abs(via.slot_dy_nm)
-        j1 = min(nx, math.floor((via.x + r - stack.x0_nm) / h) + 1)
+        j0 = max(0, math.floor((via.x - ex - stack.x0_nm) / h))
-        i0 = max(0, math.floor((via.y - r - stack.y0_nm) / h))
+        j1 = min(nx, math.floor((via.x + ex - stack.x0_nm) / h) + 1)
-        i1 = min(ny, math.floor((via.y + r - stack.y0_nm) / h) + 1)
+        i0 = max(0, math.floor((via.y - ey - stack.y0_nm) / h))
+        i1 = min(ny, math.floor((via.y + ey - stack.y0_nm) / h) + 1)
         if i0 >= i1 or j0 >= j1:
             continue
         xs = stack.x0_nm + (np.arange(j0, j1)[:, None] + sub[None, :]) * h \
             - via.x
         ys = stack.y0_nm + (np.arange(i0, i1)[:, None] + sub[None, :]) * h \
             - via.y
-        cov = ((ys[:, None, :, None] ** 2 + xs[None, :, None, :] ** 2)
+        cov = (slot_distance(xs[None, :, None, :], ys[:, None, :, None],
-               <= r * r).mean(axis=(2, 3))
+                             via.slot_dx_nm, via.slot_dy_nm)
+               <= r).mean(axis=(2, 3))
         if not (cov > 0).any():
             continue                            # mouth far smaller than h
+        span = _via_span(problem, via)
+
+        if plugged:
+            # lead cylinder + solder bore: the pin continues beyond BOTH
+            # mouths (component body / clipped stickout), so every
+            # spanned layer sees the FULL plug depth for lateral
+            # spreading - no per-layer split
+            r_lead = max(via.drill_nm - problem.tht_lead_clearance_nm,
+                         0) / 2.0
+            cov_lead = (np.hypot(xs[None, :, None, :],
+                                 ys[:, None, :, None])
+                        <= r_lead).mean(axis=(2, 3))
+            t_sn = problem.rho_ohm_m / problem.solder_rho_ohm_m
+            t_pb = problem.rho_ohm_m / problem.tht_lead_rho_ohm_m
+            depth = max(float(via.z_bot_nm - via.z_top_nm), 0.0)
+            t_eq = depth * (cov_lead * t_pb + (cov - cov_lead) * t_sn)
+            if stack.t_extra_nm is None:
+                stack.t_extra_nm = np.zeros(stack.masks.shape)
+            if stack.plug is None:
+                stack.plug = np.zeros_like(stack.masks)
+            for li in span:
+                m = stack.masks[li, i0:i1, j0:j1]
+                stack.t_extra_nm[li, i0:i1, j0:j1] += np.where(m, t_eq, 0.0)
+                stack.plug[li, i0:i1, j0:j1] |= m & (cov > 0.5)
+            continue
+
         if stack.thick_scale is None:
             stack.thick_scale = np.ones(stack.masks.shape)
-        for li in _via_span(problem, via):
+        for li in span:
-            if problem.vias_capped and li in outer:
+            if via.kind == "via" and problem.vias_capped and li in outer \
+                    and via.drill_nm <= problem.cap_max_drill_nm:
                 ratio = min(problem.cap_plating_nm
                             / problem.layers[li].thickness_nm, 1.0)
             else:
-                ratio = 0.0
+                ratio = 0.0                 # open hole (also DNP THT holes)
             s = 1.0 - cov * (1.0 - ratio)
             gone = s <= 1e-9
             stack.masks[li, i0:i1, j0:j1] &= ~gone
@@ -399,26 +528,90 @@ def _electrode_cells2d(stack: RasterStack, e: Electrode) -> np.ndarray:
     return _rect_cells(stack, e.rect)
 
 
+def _barrel_ring2d(stack: RasterStack, e: Electrode,
+                   mask2d: np.ndarray) -> np.ndarray:
+    """Contact cells of a barrel electrode on one layer: the copper ring
+    at the drill wall (cell centers within one cell of radius drill/2;
+    slotted holes: within one cell of the stadium-shaped slot wall),
+    where the lead/wire soldered into the hole actually meets the layer.
+    If rasterization or an antipad leaves no copper there, fall back to
+    the nearest copper ring within the pad footprint (+1 cell of slop) -
+    the same search bound as the solver's barrel attachment."""
+    ny, nx = stack.shape2d
+    h = stack.h_nm
+    if e.center is not None:
+        x, y = e.center
+    else:
+        x = (e.rect.x0 + e.rect.x1) / 2.0
+        y = (e.rect.y0 + e.rect.y1) / 2.0
+    r = e.drill_nm / 2.0
+    rw = max(e.pad_nm, e.drill_nm + 300_000) / 2.0 + h
+    ex, ey = rw + abs(e.slot_dx_nm), rw + abs(e.slot_dy_nm)
+    out = np.zeros((ny, nx), dtype=bool)
+    j0 = max(0, math.floor((x - ex - stack.x0_nm) / h))
+    j1 = min(nx, math.floor((x + ex - stack.x0_nm) / h) + 1)
+    i0 = max(0, math.floor((y - ey - stack.y0_nm) / h))
+    i1 = min(ny, math.floor((y + ey - stack.y0_nm) / h) + 1)
+    if i0 >= i1 or j0 >= j1:
+        return out
+    xs = stack.x0_nm + (np.arange(j0, j1) + 0.5) * h - x
+    ys = stack.y0_nm + (np.arange(i0, i1) + 0.5) * h - y
+    d = slot_distance(xs[None, :], ys[:, None], e.slot_dx_nm, e.slot_dy_nm)
+    m = mask2d[i0:i1, j0:j1]
+    ring = m & (np.abs(d - r) <= h)
+    if not ring.any():
+        dc = np.where(m & (d <= rw), d, np.inf)
+        dmin = dc.min()
+        if np.isfinite(dmin):
+            ring = dc <= dmin + h              # e.g. thermal-spoke tips
+    out[i0:i1, j0:j1] = ring
+    return out
+
+
+def _part_mask3d(stack: RasterStack, problem: Problem,
+                 el: Electrode) -> np.ndarray:
+    """(L, ny, nx) contact cells of one electrode part: the barrel-wall
+    ring on every spanned layer for via/THT-pad contacts, else the
+    part's shape ∩ copper on its contact layer(s)."""
+    part = np.zeros_like(stack.masks)
+    if el.drill_nm > 0:
+        for li, name in enumerate(stack.layer_names):
+            if el.contact not in ("all", name):
+                continue
+            if el.barrel_z is not None:
+                z = problem.layers[li].z_nm
+                if not (el.barrel_z[0] - 1 <= z <= el.barrel_z[1] + 1):
+                    continue
+            part[li] = _barrel_ring2d(stack, el, stack.masks[li])
+        return part
+    cells2d = _electrode_cells2d(stack, el)
+    for li, name in enumerate(stack.layer_names):
+        if el.contact in ("all", name):
+            part[li] = cells2d & stack.masks[li]
+    return part
+
+
 def electrode_masks(stack: RasterStack, problem: Problem
                     ) -> tuple[np.ndarray, np.ndarray]:
     """Terminal mask = OR over its parts; part = shape ∩ copper on the
-    part's contact layer(s). contact 'all' = every included layer (bolted
+    part's contact layer(s), or the barrel-wall ring for via/THT-pad
-    lug / through pad); a layer name = that layer only. Every part must
+    contacts (current enters through the soldered barrel, not the pad
-    individually land on copper (clear feedback). V+/V- must not overlap;
+    face). contact 'all' = every included layer (bolted lug / through
-    touching is checked later, only for the equipotential contact model."""
+    pad); a layer name = that layer only. Every part must individually
+    land on copper (clear feedback). V+/V- must not overlap; touching is
+    checked later, only for the equipotential contact model."""
     def build(parts: list[Electrode], which: str) -> np.ndarray:
         e = np.zeros_like(stack.masks)
         for el in parts:
-            cells2d = _electrode_cells2d(stack, el)
+            part = _part_mask3d(stack, problem, el)
-            part = np.zeros_like(stack.masks)
-            for li, name in enumerate(stack.layer_names):
-                if el.contact == "all" or el.contact == name:
-                    part[li] = cells2d & stack.masks[li]
             if not part.any():
+                where = ("near its barrel (drill-wall ring / pad footprint)"
+                         if el.drill_nm > 0 else
+                         "(or is smaller than one grid cell)")
                 raise ElectrodeError(
                     f"A {which} contact part ({el.label}) does not overlap "
                     f"any copper of the selected fill on contact layer(s) "
-                    f"'{el.contact}' (or is smaller than one grid cell)."
+                    f"'{el.contact}' {where}."
                 )
             e |= part
         if not e.any():
@@ -446,11 +639,7 @@ def electrode_partition(stack: RasterStack, problem: Problem
         out = []
         claimed = np.zeros_like(stack.masks)
         for el in parts:
-            cells2d = _electrode_cells2d(stack, el)
+            m = _part_mask3d(stack, problem, el)
-            m = np.zeros_like(stack.masks)
-            for li, name in enumerate(stack.layer_names):
-                if el.contact == "all" or el.contact == name:
-                    m[li] = cells2d & stack.masks[li]
             m &= ~claimed
             claimed |= m
             out.append((el.label, m))

fill_resistance/skin.py

@@ -21,6 +21,7 @@ from __future__ import annotations
 
 import cmath
 import math
+import re
 
 MU0 = 4e-7 * math.pi
 
@@ -58,11 +59,25 @@ def resistance_factor(thickness_m: float, freq_hz: float,
             / (rho_ohm_m / thickness_m))
 
 
+def normalize_decimal(text: str) -> str:
+    """Accept a European decimal comma ('1,5' -> '1.5'); reject
+    thousands-separator commas ('1,500' would silently become 1.5,
+    a 1000x error that propagates unnoticed into the result)."""
+    if "," in text:
+        if "." in text or text.count(",") > 1 \
+                or re.search(r",\d{3}(?=\D|$)", text):
+            raise ValueError(
+                f"ambiguous comma in '{text}': use '.' as the decimal "
+                "separator and no thousands separators")
+        text = text.replace(",", ".")
+    return text
+
+
 def parse_frequency(text: str) -> float:
     """'0', '100k', '1.5M', '142500' -> Hz; empty -> 0 (DC).
-    Raises ValueError on unparseable or negative input (a typo silently
+    Raises ValueError on unparseable, ambiguous or negative input (a
-    becoming DC would mislabel the result)."""
+    typo silently becoming DC would mislabel the result)."""
-    t = text.strip().lower().replace(",", ".").removesuffix("hz").strip()
+    t = normalize_decimal(text.strip().lower()).removesuffix("hz").strip()
     if not t:
         return 0.0
     mult = 1.0

fill_resistance/solver.py

@@ -47,7 +47,7 @@ from scipy.sparse import linalg as sla
 
 from . import config, skin
 from .errors import ConnectivityError, ElectrodeError, SolverError
-from .geometry import Problem
+from .geometry import Problem, slot_distance
 from .raster import RasterStack, electrodes_touch
 
 
@@ -156,12 +156,14 @@ def _barrel_links(stack: RasterStack, problem: Problem
         span = [li for li, layer in enumerate(problem.layers)
                 if via.spans(layer.z_nm)]
         r_nm = max(via.pad_nm, via.drill_nm + 300_000) / 2.0 + h
-        win = int(r_nm // h) + 1
+        win_j = int((r_nm + abs(via.slot_dx_nm)) // h) + 1
-        i0, i1 = max(0, i - win), min(ny, i + win + 1)
+        win_i = int((r_nm + abs(via.slot_dy_nm)) // h) + 1
-        j0, j1 = max(0, j - win), min(nx, j + win + 1)
+        i0, i1 = max(0, i - win_i), min(ny, i + win_i + 1)
+        j0, j1 = max(0, j - win_j), min(nx, j + win_j + 1)
         xs = stack.x0_nm + (np.arange(j0, j1) + 0.5) * h - via.x
         ys = stack.y0_nm + (np.arange(i0, i1) + 0.5) * h - via.y
-        d2 = ys[:, None] ** 2 + xs[None, :] ** 2
+        d2 = slot_distance(xs[None, :], ys[:, None],
+                           via.slot_dx_nm, via.slot_dy_nm) ** 2
         d2 = np.where(d2 <= r_nm * r_nm, d2, np.inf)
         present = []                            # (layer, i, j) per layer
         for li in span:
@@ -178,8 +180,16 @@ def _barrel_links(stack: RasterStack, problem: Problem
             length = problem.layers[lb].z_nm - problem.layers[la].z_nm
             if length <= 0:
                 continue
-            r_dc = via.barrel_resistance(length, problem.rho_ohm_m,
+            # populated THT pads carry a soldered component lead: lead
-                                         problem.plating_nm)
+            # cylinder (drill minus the fab clearance) + solder annulus
+            # in parallel with the plating (DNP pads and vias stay
+            # plating-only)
+            r_dc = via.barrel_resistance(
+                length, problem.rho_ohm_m, problem.plating_nm,
+                solder_rho_ohm_m=(problem.solder_rho_ohm_m
+                                  if via.solder_filled else None),
+                lead_nm=max(via.drill_nm - problem.tht_lead_clearance_nm, 0),
+                lead_rho_ohm_m=problem.tht_lead_rho_ohm_m)
             links.append((vi, la, ia, ja, lb, ib, jb, r_dc))
     return links, dead
 

fill_resistance/standalone.py

@@ -42,6 +42,10 @@ def main(argv=None) -> int:
     ap.add_argument("--uncapped", action="store_true",
                     help="treat vias as uncapped (open drill mouths on "
                          "all layers)")
+    ap.add_argument("--cap-max-drill", type=float, default=None,
+                    metavar="MM",
+                    help="cap only vias with drill <= this [mm]; larger "
+                         "drills stay open (default: from the dump)")
     ap.add_argument("--extra-cu-um", type=float, default=None,
                     help="override the added copper in mask openings [um]")
     ap.add_argument("--force-iterative", action="store_true",
@@ -68,6 +72,8 @@ def main(argv=None) -> int:
         problem.buildups = []
     if args.uncapped:
         problem.vias_capped = False
+    if args.cap_max_drill is not None:
+        problem.cap_max_drill_nm = int(args.cap_max_drill * 1e6)
     if args.extra_cu_um is not None:
         problem.extra_cu_nm = int(args.extra_cu_um * 1000)
     if args.layers:

metadata.json

@@ -2,7 +2,7 @@
     "$schema": "https://go.kicad.org/pcm/schemas/v2",
     "name": "Fill Resistance",
     "description": "DC/AC resistance of copper zone fills and traces between two contacts, single- or multi-layer with via coupling; current and power density maps.",
-    "description_full": "Computes the DC or AC resistance of copper zone fills and traces between two contacts (marker rectangles on User.1/User.2 and/or selected pads), single- or multi-layer: the chosen net's fills and tracks are solved as coupled finite-difference sheets linked by the net's via and through-hole-pad barrels; traces narrower than the grid become exact 1D resistor chains, and an adaptive multi-resolution grid (fine at features, coarse plane interiors, deferred-corrected) keeps large boards fast.\n\nShows per-layer rasterized maps, potential, current density and power density, reports per-via currents (via ampacity) and total dissipation at a selectable test current. At a user-set frequency the exact 1D foil/barrel skin-effect correction is applied (AC results are a rigorous lower bound). PNGs, a text summary and a re-solvable geometry dump are saved per run.\n\nNote: the first load builds the plugin's Python environment (numpy, scipy, pyamg, matplotlib, PySide6) and can take several minutes.",
+    "description_full": "Computes the DC or AC resistance of copper zone fills and traces between two contacts (marker rectangles on User.1/User.2 and/or selected pads/vias), single- or multi-layer: the chosen net's fills and tracks are solved as coupled finite-difference sheets linked by the net's via and through-hole-pad barrels. Selected vias/THT pads inject at the drill-wall barrel, and every populated THT hole carries its full solder joint (component lead, solder fill, one-sided pad coat and protruding-lead cone) with exact pad shapes and do-not-populate flags read from KiCad, conducting in-plane as its solder plug and lead on every layer it spans. Every net pad's exact copper shape is stamped on the layers it sits on, SMD as well as through-hole, and oblong (slotted) holes are modelled as their true stadium shape rather than an approximating circle. Traces narrower than the grid become exact 1D resistor chains, and an adaptive multi-resolution grid (fine at features, coarse plane interiors, deferred-corrected) keeps large boards fast.\n\nShows per-layer rasterized maps, potential, current density and power density, reports per-via currents (via ampacity) and total dissipation at a selectable test current. At a user-set frequency the exact 1D foil/barrel skin-effect correction is applied (AC results are a rigorous lower bound). PNGs, a text summary and a re-solvable geometry dump are saved per run.\n\nNote: the first load builds the plugin's Python environment (numpy, scipy, pyamg, matplotlib, PySide6) and can take several minutes.",
     "identifier": "th.co.b4l.fill-resistance",
     "type": "plugin",
     "author": {
@@ -17,7 +17,7 @@
     },
     "versions": [
         {
-            "version": "1.0.1",
+            "version": "1.2.1",
             "status": "stable",
             "kicad_version": "10.0",
             "runtime": "ipc"

pyproject.toml

@@ -0,0 +1,28 @@
+# Development environment only (uv sync / uv run). The KiCad plugin
+# manager builds the runtime venv itself from requirements.txt — keep
+# the dependency list there in sync with [project.dependencies].
+[project]
+name = "fill-resistance"
+version = "1.2.1"
+description = "DC/AC resistance of copper zone fills and traces between two contacts (KiCad 10 plugin)"
+license = "GPL-3.0-or-later"
+requires-python = ">=3.11"
+dependencies = [
+    "kicad-python>=0.7.0",
+    "numpy",
+    "scipy",
+    "pyamg",
+    "matplotlib",
+    "PySide6",
+]
+
+[dependency-groups]
+dev = [
+    "pytest",
+]
+
+[tool.uv]
+package = false
+
+[tool.pytest.ini_options]
+testpaths = ["tests"]

tests/test_adaptive.py

@@ -166,6 +166,32 @@ def test_part_currents_and_ac(monkeypatch):
     assert ada.rs_ratios == ref.rs_ratios
 
 
+def test_stitching_pad_mid_plane_close(monkeypatch):
+    """A solder-filled THT stitching pad mid-pour leaves no keep-fine
+    marker of its own (mouth not cut, thick_scale untouched, no copper
+    boundary nearby): without the barrel-attachment pinning its links
+    landed in a coarse equipotential leaf and the local spreading
+    resistance vanished - R read ~20% low on this exact case."""
+    sq = [(0, 0), (40, 0), (40, 40), (0, 40)]
+
+    def prob():
+        p = make_multilayer(
+            [[(sq, [])], [(sq, [])]],
+            rect1_mm=(0, 15, 2, 25), rect2_mm=(38, 15, 40, 25),
+            contact1="L0", contact2="L1",
+            vias_mm=[(20, 20)], gap_mm=1.6, drill_mm=1.0)
+        v = p.vias[0]
+        v.kind = "pad"
+        v.pad_nm = int(1.8 * NM)
+        v.solder_filled = True
+        return p
+
+    ref = _run(prob(), 0.15, adaptive=False, monkeypatch=monkeypatch)
+    ada = _run(prob(), 0.15, adaptive=True, monkeypatch=monkeypatch)
+    assert ada.n_free < 0.4 * ref.n_free       # pour still coarsens
+    assert ada.R_ohm == pytest.approx(ref.R_ohm, rel=2e-3)
+
+
 def test_auto_cell_size_finer_with_adaptive(monkeypatch):
     """The auto sizer affords a larger fine-cell budget (finer h) when
     the adaptive grid is on."""

tests/test_barrel_contacts.py

@@ -0,0 +1,542 @@
+"""Barrel (via / through-hole pad) contact tests: current enters at the
+drill-wall ring, not the pad face, and soldered THT joints carry a
+solder-filled hole plus an average-thickness solder coat on the pad."""
+import math
+
+import numpy as np
+import pytest
+
+from fill_resistance import raster, solver
+from fill_resistance.geometry import (Electrode, Polygon, ViaLink,
+                                      contact_solder_buildups, load_problem,
+                                      problem_from_json, problem_to_json,
+                                      save_problem, tht_joint_buildups)
+from tests.util import NM, make_multilayer, make_problem, rect_mm, ring_mm
+
+PLATE20 = [(0, 0), (20, 0), (20, 20), (0, 20)]
+
+
+def _barrel(x_mm, y_mm, drill_mm, pad_mm=0.0, solder=False, polygons=None):
+    r = max(pad_mm, drill_mm) / 2
+    return Electrode(
+        rect=rect_mm((x_mm - r, y_mm - r, x_mm + r, y_mm + r)),
+        contact="all", label=f"via({x_mm},{y_mm})",
+        drill_nm=int(drill_mm * NM), pad_nm=int(pad_mm * NM),
+        center=(int(x_mm * NM), int(y_mm * NM)), solder=solder,
+        polygons=polygons)
+
+
+def _disc(x_mm, y_mm, r_mm, n=64) -> Polygon:
+    ang = np.linspace(0, 2 * np.pi, n, endpoint=False)
+    return Polygon(outline=ring_mm(
+        [(x_mm + r_mm * np.cos(a), y_mm + r_mm * np.sin(a)) for a in ang]))
+
+
+def _solve(p, h_mm, model="equipotential"):
+    stack = raster.rasterize_stack(p, h_mm * NM)
+    e1, e2 = raster.electrode_masks(stack, p)
+    return solver.run_solve(p, stack, e1, e2, 1.0, contact_model=model), stack
+
+
+def test_ring_cells_at_drill_wall():
+    """The contact cells of a barrel electrode form a ring at the drill
+    wall (one-cell tolerance), not the pad face."""
+    p = make_problem([(PLATE20, [])],
+                     rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+    p.electrodes1 = [_barrel(10, 10, drill_mm=1.0, pad_mm=1.6)]
+    stack = raster.rasterize_stack(p, 0.1 * NM)
+    e1, _ = raster.electrode_masks(stack, p)
+    ii, jj = np.nonzero(e1[0])
+    xs = stack.x0_nm + (jj + 0.5) * stack.h_nm - 10 * NM
+    ys = stack.y0_nm + (ii + 0.5) * stack.h_nm - 10 * NM
+    d = np.hypot(xs, ys)
+    assert len(ii) >= 8
+    assert (np.abs(d - 0.5 * NM) <= stack.h_nm + 1).all()
+    # far fewer cells than the full 1.6 mm pad disc
+    assert len(ii) < 0.5 * math.pi * (0.8 * NM / stack.h_nm) ** 2
+
+
+def test_two_barrel_contacts_match_acosh():
+    """Two equipotential circular contacts of radius a, centers d apart,
+    on a large sheet: R = rho/(pi t) * acosh(d / 2a). The barrel-ring
+    contact must reproduce the analytic spreading resistance."""
+    t_um, rho = 70.0, 1.68e-8
+    plate = [(0, 0), (80, 0), (80, 60), (0, 60)]
+    p = make_problem([(plate, [])], rect1_mm=(0, 0, 1, 1),
+                     rect2_mm=(79, 59, 80, 60), t_um=t_um, rho=rho)
+    p.electrodes1 = [_barrel(30, 30, drill_mm=2.0)]
+    p.electrodes2 = [_barrel(50, 30, drill_mm=2.0)]
+    res, _ = _solve(p, 0.15)
+    r_ref = rho / (math.pi * t_um * 1e-6) * math.acosh(20e-3 / (2 * 1e-3))
+    assert res.R_ohm == pytest.approx(r_ref, rel=0.08)
+
+
+def test_barrel_includes_pad_spreading_resistance():
+    """Injecting at the barrel wall (0.5 mm ring) sees the spreading
+    resistance the whole-pad-face contact (2.4 mm equipotential disc)
+    short-circuits: R_barrel > R_pad_face."""
+    p1 = make_problem([(PLATE20, [])],
+                      rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+    p1.electrodes1 = [_barrel(10, 10, drill_mm=1.0, pad_mm=2.4)]
+    r_barrel, _ = _solve(p1, 0.1)
+
+    p2 = make_problem([(PLATE20, [])],
+                      rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+    p2.electrodes1 = [Electrode(rect=rect_mm((8.8, 8.8, 11.2, 11.2)),
+                                contact="all", label="pad face",
+                                polygons=[_disc(10, 10, 1.2)])]
+    r_face, _ = _solve(p2, 0.1)
+    assert r_barrel.R_ohm > r_face.R_ohm * 1.05
+
+
+def test_ring_fallback_nearest_copper():
+    """Antipad bigger than the drill: no copper at the wall ring, the
+    contact falls back to the nearest copper ring inside the pad
+    footprint (e.g. thermal-spoke tips / hole edge)."""
+    hole = [(10 + 1.2 * np.cos(a), 10 + 1.2 * np.sin(a))
+            for a in np.linspace(0, 2 * np.pi, 64, endpoint=False)]
+    p = make_problem([(PLATE20, [hole])],
+                     rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+    p.electrodes1 = [_barrel(10, 10, drill_mm=0.6, pad_mm=4.0)]
+    res, stack = _solve(p, 0.1)
+    e1, _ = raster.electrode_masks(stack, p)
+    ii, jj = np.nonzero(e1[0])
+    d = np.hypot(stack.x0_nm + (jj + 0.5) * stack.h_nm - 10 * NM,
+                 stack.y0_nm + (ii + 0.5) * stack.h_nm - 10 * NM)
+    assert len(ii) >= 8
+    assert (d >= 1.2 * NM - stack.h_nm).all()
+    assert (d <= 1.2 * NM + 2.5 * stack.h_nm).all()
+    assert np.isfinite(res.R_ohm) and res.R_ohm > 0
+
+
+def test_solder_filled_barrel_resistance():
+    """THT joints: the solder core conducts in parallel with the plating.
+    Exact parallel-area formula, and a sanity ratio for a 1 mm drill."""
+    v = ViaLink(x=0, y=0, drill_nm=1_000_000, z_top_nm=-1, z_bot_nm=1)
+    rho, sn = 1.68e-8, 1.32e-7
+    r_plain = v.barrel_resistance(1_600_000, rho, 18_000)
+    r_fill = v.barrel_resistance(1_600_000, rho, 18_000,
+                                 solder_rho_ohm_m=sn)
+    ga = math.pi * 1e-3 * 18e-6 / rho
+    ga += math.pi * (0.5e-3 - 18e-6) ** 2 / sn
+    assert r_fill == pytest.approx(1.6e-3 / ga, rel=1e-12)
+    assert 1.5 < r_plain / r_fill < 4.0
+
+
+def test_contact_solder_coat():
+    """A soldered THT contact adds an average-thickness solder buildup
+    over the pad face on its SOLDER side only (opposite the component),
+    lowering the spreading resistance vs the bare barrel contact."""
+    def prob():
+        p = make_problem([(PLATE20, [])],
+                         rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+        p.electrodes1 = [_barrel(10, 10, drill_mm=1.0, pad_mm=2.4,
+                                 solder=True, polygons=[_disc(10, 10, 1.2)])]
+        p.electrodes1[0].protrusion_side = "F.Cu"
+        return p
+
+    # solder side not among the included layers -> no coat there
+    q = prob()
+    q.electrodes1[0].protrusion_side = "B.Cu"
+    assert contact_solder_buildups(q) == []
+
+    p = prob()
+    assert contact_solder_buildups(p) == ["F.Cu"]
+    assert len(p.buildups) == 1 and p.buildups[0].layer_name == "F.Cu"
+    r_coat, stack = _solve(p, 0.1)
+    assert stack.buildup is not None and stack.buildup.any()
+
+    r_bare, _ = _solve(prob(), 0.1)          # helper not called: no coat
+    assert r_coat.R_ohm < r_bare.R_ohm
+
+
+def test_lead_fillet_profile():
+    """The protruding-lead solder cone paints thick_scale with the exact
+    per-cell formula: 1 + H*clip((rb-r)/(rb-ra), 0, 1)*(rho_cu/rho_sn)/t
+    on copper of the protrusion side; nothing elsewhere."""
+    p = make_problem([(PLATE20, [])],
+                     rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+    p.electrodes1 = [_barrel(10, 10, drill_mm=1.0, pad_mm=2.4, solder=True)]
+    p.electrodes1[0].protrusion_side = "F.Cu"
+    stack = raster.rasterize_stack(p, 0.1 * NM)
+    assert stack.thick_scale is not None
+    ny, nx = stack.shape2d
+    jj, ii = np.meshgrid(np.arange(nx), np.arange(ny))
+    r = np.hypot(stack.x0_nm + (jj + 0.5) * stack.h_nm - 10 * NM,
+                 stack.y0_nm + (ii + 0.5) * stack.h_nm - 10 * NM)
+    ra, rb, H = 0.5 * NM, 1.2 * NM, p.tht_protrusion_nm
+    t_eq = H * np.clip((rb - r) / (rb - ra), 0, 1) \
+        * (p.rho_ohm_m / p.solder_rho_ohm_m)
+    expect = np.where(stack.masks[0],
+                      1.0 + t_eq / p.layers[0].thickness_nm, 1.0)
+    assert np.allclose(stack.thick_scale[0], expect, rtol=1e-12)
+    # 1.5 mm of solder at the wall ~ 191 um copper: factor ~ 3.7 on 70 um
+    assert stack.thick_scale[0].max() > 3.0
+
+    p.electrodes1[0].protrusion_side = None    # e.g. via contact: no cone
+    s2 = raster.rasterize_stack(p, 0.1 * NM)
+    assert s2.thick_scale is None
+
+
+def test_lead_fillet_lowers_resistance(monkeypatch):
+    """The cone shorts the joint vicinity: R(with cone) < R(coat-less
+    bare barrel); the adaptive grid pins the cone cells fine and
+    matches the uniform grid."""
+    def prob(protrude=True):
+        p = make_problem([(PLATE20, [])],
+                         rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+        p.electrodes1 = [_barrel(10, 10, drill_mm=1.0, pad_mm=2.4,
+                                 solder=True)]
+        p.electrodes1[0].protrusion_side = "F.Cu"
+        if not protrude:
+            p.tht_protrusion_nm = 0
+        return p
+
+    r_cone, _ = _solve(prob(), 0.1)
+    r_bare, _ = _solve(prob(protrude=False), 0.1)
+    assert r_cone.R_ohm < r_bare.R_ohm
+
+    from fill_resistance import config
+    monkeypatch.setattr(config, "ADAPTIVE_CELLS", True)
+    r_ada, _ = _solve(prob(), 0.1)
+    assert r_ada.R_ohm == pytest.approx(r_cone.R_ohm, rel=2e-3)
+
+
+def _pad_link(populated=True):
+    return ViaLink(x=10 * NM, y=10 * NM, drill_nm=1_000_000, z_top_nm=-1,
+                   z_bot_nm=1, kind="pad", pad_nm=2_400_000,
+                   solder_filled=populated,
+                   protrusion_side="F.Cu" if populated else None)
+
+
+def test_stitching_pad_joint():
+    """A populated THT pad on the net (not a contact) gets the full
+    joint: solder-side coat, cone, and a conducting (plugged) mouth;
+    a DNP pad gets an open hole and nothing else."""
+    def prob(populated=True):
+        p = make_problem([(PLATE20, [])],
+                         rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+        p.vias = [_pad_link(populated)]
+        return p
+
+    p = prob()
+    assert tht_joint_buildups(p) == ["F.Cu"]
+    assert len(p.buildups) == 1
+    r_joint, stack = _solve(p, 0.1)
+    assert stack.thick_scale is not None and stack.thick_scale.max() > 3.0
+    assert stack.buildup is not None and stack.buildup.any()
+    assert stack.masks[0][stack.cell_of(10 * NM, 10 * NM)]   # plugged mouth
+
+    q = prob(populated=False)
+    assert tht_joint_buildups(q) == []
+    r_bare, s2 = _solve(q, 0.1)
+    assert s2.buildup is None
+    assert not s2.masks[0][s2.cell_of(10 * NM, 10 * NM)]     # DNP: open hole
+    assert r_joint.R_ohm < r_bare.R_ohm
+
+
+def test_cone_not_doubled_at_contact():
+    """A contact THT pad also appears in the net's pad list (ViaLink):
+    the cone and coat must be applied once, not squared/stacked. The
+    hole plug (this synthetic barrel spans z = -1..1, so 2 nm of lead)
+    ADDS to the cone at the mouth instead of multiplying it."""
+    p = make_problem([(PLATE20, [])],
+                     rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+    p.electrodes1 = [_barrel(10, 10, drill_mm=1.0, pad_mm=2.4, solder=True,
+                             polygons=[_disc(10, 10, 1.2)])]
+    p.electrodes1[0].protrusion_side = "F.Cu"
+    p.vias = [_pad_link()]
+    assert contact_solder_buildups(p) == ["F.Cu"]
+    assert tht_joint_buildups(p) == []       # contact center is skipped
+    stack = raster.rasterize_stack(p, 0.1 * NM)
+    t_cone = p.tht_protrusion_nm * (p.rho_ohm_m / p.solder_rho_ohm_m)
+    t_plug = 2.0 * (p.rho_ohm_m / p.tht_lead_rho_ohm_m)
+    wall = 1.0 + (t_cone + t_plug) / p.layers[0].thickness_nm
+    assert stack.thick_scale.max() == pytest.approx(wall, rel=1e-12)
+
+
+def test_lead_in_barrel_resistance():
+    """Populated hole: plating || lead cylinder || solder annulus, with
+    the lead clipped to the plating bore."""
+    v = ViaLink(x=0, y=0, drill_nm=1_000_000, z_top_nm=-1, z_bot_nm=1)
+    rho, sn = 1.68e-8, 1.32e-7
+    r_solder = v.barrel_resistance(1_600_000, rho, 18_000,
+                                   solder_rho_ohm_m=sn)
+    r_lead = v.barrel_resistance(1_600_000, rho, 18_000,
+                                 solder_rho_ohm_m=sn,
+                                 lead_nm=750_000, lead_rho_ohm_m=rho)
+    rl, rc = 0.375e-3, 0.5e-3 - 18e-6
+    ga = math.pi * 1e-3 * 18e-6 / rho
+    ga += math.pi * rl ** 2 / rho + math.pi * (rc ** 2 - rl ** 2) / sn
+    assert r_lead == pytest.approx(1.6e-3 / ga, rel=1e-12)
+    assert r_lead < r_solder
+    # a lead wider than the bore is clipped to it
+    r_big = v.barrel_resistance(1_600_000, rho, 18_000,
+                                solder_rho_ohm_m=sn,
+                                lead_nm=2_000_000, lead_rho_ohm_m=rho)
+    ga2 = math.pi * 1e-3 * 18e-6 / rho + math.pi * rc ** 2 / rho
+    assert r_big == pytest.approx(1.6e-3 / ga2, rel=1e-12)
+
+
+def test_oblong_pad_cone_uses_inscribed_dim():
+    """Oblong pads: the cone tapers to the inscribed circle (pad_min),
+    never past it, so the long pad axis is not overstated sideways."""
+    p = make_problem([(PLATE20, [])],
+                     rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+    p.vias = [_pad_link()]
+    p.vias[0].pad_min_nm = 1_600_000          # 2.4 mm max, 1.6 mm min
+    stack = raster.rasterize_stack(p, 0.1 * NM)
+    ii, jj = np.nonzero(stack.thick_scale[0] != 1.0)
+    d = np.hypot(stack.x0_nm + (jj + 0.5) * stack.h_nm - 10 * NM,
+                 stack.y0_nm + (ii + 0.5) * stack.h_nm - 10 * NM)
+    assert len(d) and d.max() < 0.8 * NM
+
+
+def test_stitching_coat_exact_shape():
+    """When KiCad supplies the exact pad polygon, the coat uses it
+    instead of the pad-diameter disc (oblong pads stay honest)."""
+    p = make_problem([(PLATE20, [])],
+                     rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+    p.vias = [_pad_link()]
+    shape = _disc(10, 10, 0.9)
+    assert tht_joint_buildups(p, {(10 * NM, 10 * NM): [shape]}) == ["F.Cu"]
+    assert p.buildups[0].polygons[0] is shape
+
+
+def test_vialink_solder_json():
+    p = make_problem([(PLATE20, [])],
+                     rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+    p.vias = [_pad_link()]
+    d = problem_to_json(p)
+    q = problem_from_json(d)
+    assert q.vias[0].solder_filled is True
+    assert q.vias[0].protrusion_side == "F.Cu"
+    # legacy dumps without the flag: THT pads counted as solder-filled,
+    # vias as plating-only
+    del d["vias"][0]["solder_filled"], d["vias"][0]["protrusion_side"]
+    q = problem_from_json(d)
+    assert q.vias[0].solder_filled is True
+    assert q.vias[0].protrusion_side is None
+    d["vias"][0]["kind"] = "via"
+    assert problem_from_json(d).vias[0].solder_filled is False
+
+
+# --- slotted (oblong) holes --------------------------------------------------
+# The lead/barrel of a slotted hole is a stadium, not a circle: modeling
+# it as a circle of the slot's LONG dimension painted contact rings,
+# mouths and cones bigger than the oblong pad itself.
+
+def _slot_dist_mm(stack, ii, jj, x_mm, y_mm, dx_nm):
+    """Distance of cells (ii, jj) to a slot axis (+-dx_nm along x)."""
+    xs = stack.x0_nm + (jj + 0.5) * stack.h_nm - x_mm * NM
+    ys = stack.y0_nm + (ii + 0.5) * stack.h_nm - y_mm * NM
+    t = np.clip(xs / dx_nm, -1.0, 1.0)
+    return np.hypot(xs - t * dx_nm, ys), xs, ys
+
+
+def test_slot_ring_hugs_slot_wall():
+    """The contact ring of a slotted THT pad follows the stadium-shaped
+    slot wall: it reaches around the end caps but never pokes past the
+    oblong pad's short side (the old circular model of the slot's long
+    dimension put cells at radius 1.5 mm straight above/below)."""
+    p = make_problem([(PLATE20, [])],
+                     rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+    e = _barrel(10, 10, drill_mm=1.0, pad_mm=3.6)   # slot 3.0 x 1.0 mm
+    e.pad_min_nm = int(1.6 * NM)                    # pad 3.6 x 1.6 mm
+    e.slot_dx_nm = 1 * NM
+    p.electrodes1 = [e]
+    stack = raster.rasterize_stack(p, 0.1 * NM)
+    e1, _ = raster.electrode_masks(stack, p)
+    ii, jj = np.nonzero(e1[0])
+    d, xs, ys = _slot_dist_mm(stack, ii, jj, 10, 10, 1 * NM)
+    assert len(ii) >= 16
+    assert (np.abs(d - 0.5 * NM) <= stack.h_nm + 1).all()
+    assert xs.max() > 1.2 * NM and xs.min() < -1.2 * NM   # rings the caps
+    assert np.abs(ys).max() < 0.8 * NM      # stays inside the 1.6 mm side
+
+
+def test_slot_mouth_is_stadium():
+    """A DNP slotted pad cuts a stadium-shaped hole: open along the whole
+    slot, copper kept just past the slot width and the end caps."""
+    p = make_problem([(PLATE20, [])],
+                     rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+    v = _pad_link(populated=False)
+    v.slot_dx_nm = 1 * NM                   # slot 3.0 x 1.0 mm along x
+    p.vias = [v]
+    stack = raster.rasterize_stack(p, 0.1 * NM)
+    m = stack.masks[0]
+    assert not m[stack.cell_of(10 * NM, 10 * NM)]
+    assert not m[stack.cell_of(int(10.9 * NM), 10 * NM)]  # slot end: open
+    assert not m[stack.cell_of(int(9.1 * NM), 10 * NM)]
+    assert m[stack.cell_of(10 * NM, int(10.8 * NM))]  # past the width: copper
+    assert m[stack.cell_of(10 * NM, int(9.2 * NM))]
+    assert m[stack.cell_of(int(11.8 * NM), 10 * NM)]  # past the cap: copper
+
+
+def test_slot_cone_follows_slot():
+    """The lead cone of a slotted oblong pad tapers from the slot WALL
+    to the pad's short dimension. The old circular-drill model (diameter
+    = the slot's long dimension) skipped the cone entirely
+    (pad_min <= drill) and, for the mouth, ate the pad's short side."""
+    p = make_problem([(PLATE20, [])],
+                     rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+    e = _barrel(10, 10, drill_mm=1.0, pad_mm=3.6, solder=True)
+    e.pad_min_nm = int(1.6 * NM)
+    e.slot_dx_nm = 1 * NM
+    e.protrusion_side = "F.Cu"
+    p.electrodes1 = [e]
+    stack = raster.rasterize_stack(p, 0.1 * NM)
+    assert stack.thick_scale is not None
+    ny, nx = stack.shape2d
+    jj, ii = np.meshgrid(np.arange(nx), np.arange(ny))
+    r, _, _ = _slot_dist_mm(stack, ii, jj, 10, 10, 1 * NM)
+    ra, rb, H = 0.5 * NM, 0.8 * NM, p.tht_protrusion_nm
+    t_eq = H * np.clip((rb - r) / (rb - ra), 0, 1) \
+        * (p.rho_ohm_m / p.solder_rho_ohm_m)
+    expect = np.where(stack.masks[0],
+                      1.0 + t_eq / p.layers[0].thickness_nm, 1.0)
+    assert np.allclose(stack.thick_scale[0], expect, rtol=1e-12)
+    assert stack.thick_scale[0].max() > 3.0
+
+
+def test_plug_conducts_on_component_side():
+    """A populated THT pad's filled hole (lead + solder plug) conducts
+    IN-PLANE across the mouth on EVERY spanned layer - the component
+    side is not bare foil. Each layer carries the FULL hole depth (the
+    pin continues beyond both mouths, so the whole plug cross-section
+    spreads current at every layer; side-to-side the only difference
+    is the solder coat + cone), converted to conduction-equivalent
+    copper: lead disc at lead resistivity, solder bore around it."""
+    p = make_multilayer([[(PLATE20, [])], [(PLATE20, [])]],
+                        rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+    p.vias = [ViaLink(x=10 * NM, y=10 * NM, drill_nm=1_000_000, z_top_nm=-1,
+                      z_bot_nm=1 * NM + 1, kind="pad", pad_nm=2_400_000,
+                      solder_filled=True, protrusion_side="L0")]
+    stack = raster.rasterize_stack(p, 0.1 * NM)
+    c = stack.cell_of(10 * NM, 10 * NM)
+    assert stack.masks[0][c] and stack.masks[1][c]     # plugged, not open
+    assert stack.plug[0][c] and stack.plug[1][c]       # drawn on both sides
+    t = p.layers[0].thickness_nm
+    # full hole depth z = -1 .. 1 mm + 1 on both layers; the mouth
+    # center lies inside the 0.75 mm lead (copper resistivity)
+    depth = 1 * NM + 2.0
+    t_cone = p.tht_protrusion_nm * (p.rho_ohm_m / p.solder_rho_ohm_m)
+    assert stack.thick_scale[0][c] == pytest.approx(
+        1.0 + (t_cone + depth * (p.rho_ohm_m / p.tht_lead_rho_ohm_m)) / t,
+        rel=1e-9)                                      # solder side: + cone
+    assert stack.thick_scale[1][c] == pytest.approx(
+        1.0 + depth * (p.rho_ohm_m / p.tht_lead_rho_ohm_m) / t, rel=1e-9)
+    # far from the joint: untouched foil
+    assert stack.thick_scale[1][stack.cell_of(14 * NM, 10 * NM)] == 1.0
+
+    # clearance swallowing the bore -> no lead, solder-only plug
+    p.tht_lead_clearance_nm = 1_000_000
+    s_sn = raster.rasterize_stack(p, 0.1 * NM)
+    assert s_sn.thick_scale[1][c] == pytest.approx(
+        1.0 + depth * (p.rho_ohm_m / p.solder_rho_ohm_m) / t, rel=1e-9)
+    p.tht_lead_clearance_nm = 250_000
+
+    # a DNP pad still cuts an open hole and gets no plug
+    p.vias[0].solder_filled = False
+    p.vias[0].protrusion_side = None
+    s2 = raster.rasterize_stack(p, 0.1 * NM)
+    assert not s2.masks[0][c] and not s2.masks[1][c]
+    assert s2.plug is None
+
+
+def test_slot_barrel_resistance():
+    """Slotted barrel: plating wall = stadium perimeter, solder core =
+    stadium bore area (both reduce to the circle for dx = dy = 0)."""
+    v = ViaLink(x=0, y=0, drill_nm=1_000_000, z_top_nm=-1, z_bot_nm=1,
+                slot_dx_nm=800_000, slot_dy_nm=600_000)   # ext = 2 mm
+    rho, sn = 1.68e-8, 1.32e-7
+    ga = (math.pi * 1e-3 + 2 * 2e-3) * 18e-6 / rho
+    r_plain = v.barrel_resistance(1_600_000, rho, 18_000)
+    assert r_plain == pytest.approx(1.6e-3 / ga, rel=1e-12)
+    rc = 0.5e-3 - 18e-6
+    ga += (math.pi * rc * rc + 2 * rc * 2e-3) / sn
+    r_fill = v.barrel_resistance(1_600_000, rho, 18_000, solder_rho_ohm_m=sn)
+    assert r_fill == pytest.approx(1.6e-3 / ga, rel=1e-12)
+
+
+def test_slot_coat_fallback_within_pad():
+    """Without an exact pad shape the stitching coat falls back to a
+    capsule along the slot (width = pad_min), not the old pad_nm disc
+    that stuck out past an oblong pad's short side."""
+    p = make_problem([(PLATE20, [])],
+                     rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+    v = _pad_link()                          # pad_nm = 2.4 mm
+    v.pad_min_nm = 1_600_000
+    v.slot_dx_nm = 1 * NM
+    p.vias = [v]
+    assert tht_joint_buildups(p) == ["F.Cu"]
+    pts = p.buildups[0].polygons[0].outline.astype(float)
+    xs, ys = pts[:, 0] - 10 * NM, pts[:, 1] - 10 * NM
+    t = np.clip(xs / (0.4 * NM), -1.0, 1.0)  # caps at +-(2.4-1.6)/2 mm
+    d = np.hypot(xs - t * 0.4 * NM, ys)
+    assert np.allclose(d, 0.8 * NM, atol=2)
+    assert np.abs(xs).max() <= 1.2 * NM + 2  # never past pad_nm / 2
+    assert np.abs(ys).max() <= 0.8 * NM + 2  # never past pad_min / 2
+
+
+def test_slot_json_roundtrip():
+    p = make_problem([(PLATE20, [])],
+                     rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+    e = _barrel(10, 10, drill_mm=1.0, pad_mm=3.6)
+    e.slot_dx_nm, e.slot_dy_nm = 700_000, -700_000
+    p.electrodes1 = [e]
+    v = _pad_link()
+    v.slot_dx_nm = 1 * NM
+    p.vias = [v]
+    q = problem_from_json(problem_to_json(p))
+    assert (q.electrodes1[0].slot_dx_nm, q.electrodes1[0].slot_dy_nm) \
+        == (700_000, -700_000)
+    assert (q.vias[0].slot_dx_nm, q.vias[0].slot_dy_nm) == (1 * NM, 0)
+    # legacy dumps: round drills
+    d = problem_to_json(p)
+    for vd in d["vias"]:
+        del vd["slot_dx_nm"], vd["slot_dy_nm"]
+    assert problem_from_json(d).vias[0].slot_dx_nm == 0
+
+
+def test_drill_info_slot_rotation():
+    """_drill_info: slot axis from the drill x/y sizes, rotated with the
+    pad (KiCad angles are CCW with y down: 90 deg sends +x to -y)."""
+    from types import SimpleNamespace as NS
+
+    from fill_resistance.board_io import _drill_info
+
+    def pad(dx_mm, dy_mm, angle_deg):
+        return NS(padstack=NS(
+            drill=NS(diameter=NS(x=int(dx_mm * NM), y=int(dy_mm * NM))),
+            angle=NS(degrees=angle_deg)))
+
+    assert _drill_info(pad(1.0, 1.0, 0.0)) == (1 * NM, 0, 0)    # round
+    assert _drill_info(pad(3.0, 1.0, 0.0)) == (1 * NM, 1 * NM, 0)
+    assert _drill_info(pad(1.0, 3.0, 0.0)) == (1 * NM, 0, 1 * NM)
+    w, dx, dy = _drill_info(pad(3.0, 1.0, 90.0))
+    assert (w, dx, dy) == (1 * NM, 0, -1 * NM)
+    w, dx, dy = _drill_info(pad(3.0, 1.0, 45.0))
+    assert w == 1 * NM
+    assert dx == pytest.approx(1 * NM / math.sqrt(2), abs=2)
+    assert dy == pytest.approx(-1 * NM / math.sqrt(2), abs=2)
+
+
+def test_barrel_electrode_json_roundtrip(tmp_path):
+    p = make_problem([(PLATE20, [])],
+                     rect1_mm=(0, 0, 1, 20), rect2_mm=(19, 0, 20, 20))
+    p.electrodes1 = [_barrel(10, 10, drill_mm=0.6, pad_mm=1.2, solder=True,
+                             polygons=[_disc(10, 10, 0.6)])]
+    p.electrodes1[0].barrel_z = (-1, 1_600_001)
+    p.electrodes1[0].protrusion_side = "B.Cu"
+    p.tht_protrusion_nm = 1_200_000
+    f = tmp_path / "d.json"
+    save_problem(p, f)
+    q = load_problem(f)
+    e = q.electrodes1[0]
+    assert e.drill_nm == 600_000 and e.pad_nm == 1_200_000
+    assert e.center == (10 * NM, 10 * NM)
+    assert e.barrel_z == (-1, 1_600_001)
+    assert e.solder is True and len(e.polygons) == 1
+    assert e.protrusion_side == "B.Cu"
+    assert q.tht_protrusion_nm == 1_200_000

tests/test_board_io_api.py

@@ -0,0 +1,201 @@
+"""board_io's kipy-facing paths, against a fake board.
+
+Real protobuf messages, a fake transport. These cover what a live KiCad
+would otherwise be needed for: the overlay push (kipy's
+Board.remove_items discards the DeleteItemsResponse, so board_io talks
+to the proto layer directly and these pin the status handling that
+depends on) and per-layer pad copper selection.
+"""
+from types import SimpleNamespace as NS
+
+import numpy as np
+import pytest
+from kipy.proto.common.commands.editor_commands_pb2 import (
+    CreateItemsResponse, DeleteItemsResponse, ItemDeletionStatus)
+from kipy.proto.common.types.base_types_pb2 import KIID
+from kipy.util.board_layer import layer_from_canonical_name
+
+from fill_resistance import board_io, config
+
+
+class _Ref:
+    """Stand-in for a reference image already on the board (kipy board
+    items carry a KIID message, not a bare id)."""
+    def __init__(self, layer_name, ident):
+        self.layer = layer_from_canonical_name(layer_name)
+        self.id = KIID(value=f"00000000-0000-0000-0000-{ident:012d}")
+
+
+class _FakeKiCad:
+    def __init__(self, delete_status=ItemDeletionStatus.IDS_OK):
+        self.delete_status = delete_status
+        self.deleted = []          # layers we were asked to clear
+        self.created = []          # ReferenceImages we were asked to add
+
+    def send(self, cmd, response_type):
+        if response_type is DeleteItemsResponse:
+            resp = DeleteItemsResponse()
+            for _ in cmd.item_ids:
+                resp.deleted_items.add().status = self.delete_status
+            self.deleted.append(len(cmd.item_ids))
+            return resp
+        if response_type is CreateItemsResponse:
+            resp = CreateItemsResponse()
+            resp.created_items.add().status.code = 1      # ISC_OK
+            self.created.append(cmd)
+            return resp
+        raise AssertionError(f"unexpected command {type(cmd).__name__}")
+
+
+class _FakeBoard:
+    def __init__(self, existing=(), delete_status=ItemDeletionStatus.IDS_OK):
+        self._kicad = _FakeKiCad(delete_status)
+        self._refs = list(existing)
+        self.commits = []
+        self.pushed = []
+        self.dropped = []
+
+    # kipy Board surface board_io actually uses
+    @property
+    def _doc(self):
+        from kipy.proto.common.types.base_types_pb2 import DocumentSpecifier
+        return DocumentSpecifier()
+
+    def get_reference_images(self):
+        return list(self._refs)
+
+    def begin_commit(self):
+        self.commits.append("open")
+        return object()
+
+    def push_commit(self, commit, message=""):
+        self.pushed.append(message)
+
+    def drop_commit(self, commit):
+        self.dropped.append(commit)
+
+
+class _Stack:
+    layer_names = ["F.Cu", "B.Cu"]
+    shape2d = (12, 16)
+    h_nm = 100_000
+    x0_nm = 0
+    y0_nm = 0
+
+
+class _Result:
+    def __init__(self, nlayers=2, ny=12, nx=16):
+        self.Jmag = np.full((nlayers, ny, nx), 1e6)
+
+
+def test_remove_overlays_counts_deleted():
+    layer = layer_from_canonical_name("User.9")
+    board = _FakeBoard(existing=[_Ref("User.9", 1), _Ref("User.9", 2),
+                                 _Ref("User.10", 3)])
+    assert board_io.remove_overlays(board, layer) == 2   # not the User.10 one
+
+
+def test_remove_overlays_no_images_is_a_noop():
+    board = _FakeBoard()
+    assert board_io.remove_overlays(
+        board, layer_from_canonical_name("User.9")) == 0
+    assert board._kicad.deleted == []      # no DeleteItems sent at all
+
+
+def test_locked_overlay_raises_instead_of_stacking():
+    """A locked image comes back IDS_IMMUTABLE while the overall request
+    still reports OK. Unchecked, the caller would add a second image on
+    top of the one it believed it had replaced."""
+    board = _FakeBoard(existing=[_Ref("User.9", 1)],
+                       delete_status=ItemDeletionStatus.IDS_IMMUTABLE)
+    with pytest.raises(RuntimeError, match="could not be removed"):
+        board_io.remove_overlays(board, layer_from_canonical_name("User.9"))
+
+
+def test_already_gone_overlay_is_not_an_error():
+    board = _FakeBoard(existing=[_Ref("User.9", 1)],
+                       delete_status=ItemDeletionStatus.IDS_NONEXISTENT)
+    assert board_io.remove_overlays(
+        board, layer_from_canonical_name("User.9")) == 1
+
+
+def test_push_clears_slots_this_run_does_not_write(monkeypatch):
+    """A 2-layer run after a 4-layer run must not leave the previous
+    solve's heatmap sitting on User.11/User.12."""
+    stale = [_Ref(n, i) for i, n in enumerate(config.OVERLAY_LAYERS)]
+    board = _FakeBoard(existing=stale)
+    board_io.push_result_overlays(board, _Stack(), _Result())
+
+    written = {c.items[0].type_url for c in board._kicad.created}
+    assert len(board._kicad.created) == 2          # F.Cu, B.Cu -> 2 slots
+    assert written                                  # images really created
+    # 2 written slots cleared + 2 unwritten slots cleared = 4 delete calls
+    assert len(board._kicad.deleted) == 4
+
+
+def test_push_is_one_undo_step():
+    board = _FakeBoard()
+    board_io.push_result_overlays(board, _Stack(), _Result())
+    assert board.commits and board.pushed and not board.dropped
+
+
+def _square(side):
+    """Minimal duck-typed PolygonWithHoles: an origin square."""
+    pts = [(0, 0), (side, 0), (side, side), (0, side)]
+    return NS(outline=NS(nodes=[NS(has_point=True, has_arc=False,
+                                   point=NS(x=x, y=y)) for x, y in pts]),
+              holes=[])
+
+
+class _PadBoard:
+    """F.Cu carries a small pad, B.Cu a deliberately larger one - KiCad
+    allows a different pad size per copper layer."""
+    def __init__(self):
+        self.f = layer_from_canonical_name("F.Cu")
+        self.b = layer_from_canonical_name("B.Cu")
+        self.asked = []
+
+    def get_pad_shapes_as_polygons(self, pad, layer):
+        self.asked.append(layer)
+        return {self.f: _square(1000), self.b: _square(5000)}.get(layer)
+
+
+def _width(polys):
+    xs = [p[0] for p in polys[0].outline]
+    return max(xs) - min(xs)
+
+
+def test_tht_pad_copper_comes_from_the_solder_side():
+    """The solder coat is sized from this shape, so a B.Cu-protruding
+    joint must not be measured with F.Cu's (here smaller) pad."""
+    board = _PadBoard()
+    polys = board_io._pad_polygons(board, pad=None, contact="all",
+                                   prefer="B.Cu")
+    assert _width(polys) == 5000
+    assert board.asked[0] == board.b          # probed before the F.Cu default
+
+
+def test_pad_copper_falls_back_when_no_side_is_known():
+    board = _PadBoard()
+    polys = board_io._pad_polygons(board, pad=None, contact="all")
+    assert _width(polys) == 1000              # F.Cu, the documented fallback
+
+
+def test_explicit_contact_layer_still_wins():
+    board = _PadBoard()
+    polys = board_io._pad_polygons(board, pad=None, contact="B.Cu",
+                                   prefer="F.Cu")
+    assert _width(polys) == 5000
+
+
+def test_push_drops_the_commit_if_it_cannot_finish(monkeypatch):
+    board = _FakeBoard()
+    monkeypatch.setattr(board_io.config, "OVERLAY_LAYERS", ("User.9",))
+
+    def boom(*a, **k):
+        raise RuntimeError("transport died")
+    monkeypatch.setattr(board, "push_commit", boom)
+
+    with pytest.raises(RuntimeError):
+        board_io.push_result_overlays(board, _Stack(), _Result())
+    assert board.dropped

tests/test_capping.py

@@ -9,10 +9,13 @@ from tests.util import NM, make_multilayer
 
 
 def _two_layer(width_mm=5.0, drill_mm=0.3, pad_mm=0.6, kind="via",
-               capped=True, cap_um=15.0, hole_mm=None):
+               capped=True, cap_um=15.0, hole_mm=None,
+               cap_max_drill_mm=10.0):
     """10 x width strip on both (outer-named) layers, e1 left on F.Cu,
     e2 right on B.Cu, one via mid-strip. Optionally a circular hole in
-    the F.Cu fill around the via (ring-bridging scenario)."""
+    the F.Cu fill around the via (ring-bridging scenario). The cap-drill
+    threshold defaults to 10 mm here (= every drill capped) so the tests
+    exercise the mouth treatment itself; the threshold has its own test."""
     y = width_mm / 2
     strip = [(0, 0), (10, 0), (10, width_mm), (0, width_mm)]
     holes = []
@@ -31,6 +34,7 @@ def _two_layer(width_mm=5.0, drill_mm=0.3, pad_mm=0.6, kind="via",
     p.vias[0].pad_nm = int(pad_mm * NM)
     p.vias_capped = capped
     p.cap_plating_nm = int(cap_um * 1000)
+    p.cap_max_drill_nm = int(cap_max_drill_mm * NM)
     return p
 
 
@@ -46,7 +50,13 @@ def test_cap_at_foil_thickness_is_identity():
     so the result equals the feature-off reference (a 'pad'-kind barrel,
     which skips rings and mouths) with the mouth fully inside copper."""
     r_cap, _ = _solve(_two_layer(capped=True, cap_um=70.0), 0.1)
-    r_ref, _ = _solve(_two_layer(kind="pad"), 0.1)
+    ref = _two_layer(kind="pad")
+    # populated pads skip rings and mouths; kill the lead + solder core
+    # so the reference barrel matches the via's plating-only resistance
+    ref.vias[0].solder_filled = True
+    ref.solder_rho_ohm_m = 1e30
+    ref.tht_lead_clearance_nm = 10 ** 9
+    r_ref, _ = _solve(ref, 0.1)
     assert r_cap.R_ohm == pytest.approx(r_ref.R_ohm, rel=1e-9)
 
 
@@ -90,14 +100,31 @@ def test_subcell_mouth_perturbs_gently():
     assert r_solid.R_ohm <= r_open.R_ohm <= 1.05 * r_solid.R_ohm
 
 
+def test_cap_drill_threshold():
+    """Drills above cap_max_drill_nm stay open even with capping on: a
+    2 mm drill over a 0.5 mm threshold behaves exactly like uncapped,
+    while a threshold above the drill restores the cap."""
+    kw = dict(drill_mm=2.0, pad_mm=2.6)
+    r_big, s_big = _solve(_two_layer(capped=True, cap_max_drill_mm=0.5,
+                                     **kw), 0.25)
+    r_open, s_open = _solve(_two_layer(capped=False, **kw), 0.25)
+    assert r_big.R_ohm == pytest.approx(r_open.R_ohm, rel=1e-12)
+    assert int(s_big.masks.sum()) == int(s_open.masks.sum())
+
+    r_cap, _ = _solve(_two_layer(capped=True, cap_max_drill_mm=2.1,
+                                 **kw), 0.25)
+    assert r_cap.R_ohm < r_open.R_ohm
+
+
 def test_capping_json_roundtrip(tmp_path):
     from fill_resistance.geometry import load_problem, save_problem
-    p = _two_layer(capped=False, cap_um=12.0)
+    p = _two_layer(capped=False, cap_um=12.0, cap_max_drill_mm=0.8)
     f = tmp_path / "d.json"
     save_problem(p, f)
     q = load_problem(f)
     assert q.vias_capped is False
     assert q.cap_plating_nm == 12_000
+    assert q.cap_max_drill_nm == 800_000
     r_p, _ = _solve(p, 0.25)
     r_q, _ = _solve(q, 0.25)
     assert r_q.R_ohm == pytest.approx(r_p.R_ohm, rel=1e-12)

tests/test_overlay.py

@@ -0,0 +1,60 @@
+"""In-KiCad overlay rendering (fill_resistance.overlay): copper-shaped
+RGBA heatmaps with a visibility floor and a soft edge bleed. The kipy
+pushing side is exercised only against a live KiCad (tools/)."""
+import io
+
+import numpy as np
+import pytest
+from PIL import Image
+
+from fill_resistance import config, overlay
+
+
+def _field(ny=20, nx=30):
+    """Two-layer |J| field: copper disc on layer 0, NaN elsewhere."""
+    data = np.full((2, ny, nx), np.nan)
+    yy, xx = np.mgrid[:ny, :nx]
+    disc = (yy - ny / 2) ** 2 + (xx - nx / 2) ** 2 <= 8 ** 2
+    data[0][disc] = 1.0 + xx[disc]              # spans the log range
+    data[1][disc] = 1e-12                       # below the global log floor
+    return data, disc
+
+
+def test_heatmap_png_shape_and_alpha():
+    data, disc = _field()
+    img = Image.open(io.BytesIO(overlay.heatmap_png(data, 0, bleed=False)))
+    assert img.size == (30, 20)
+    rgba = np.asarray(img)
+    assert (rgba[..., 3][disc] == config.OVERLAY_ALPHA).all()
+    assert (rgba[..., 3][~disc] == 0).all()
+
+
+def test_heatmap_floor_not_black():
+    """The coldest copper must stay distinguishable from a dark canvas:
+    the colormap starts FLOOR up, never at its near-black bottom."""
+    data, disc = _field()
+    rgba = np.asarray(Image.open(io.BytesIO(
+        overlay.heatmap_png(data, 1, bleed=False))))   # layer 1: all-cold
+    floor = np.array(__import__("matplotlib").colormaps[
+        config.CMAP_CURRENT](overlay.FLOOR)[:3]) * 255
+    assert np.abs(rgba[..., :3][disc] - floor).max() <= 1
+    assert rgba[..., :3][disc].sum(axis=-1).min() > 30   # not near-black
+
+
+def test_heatmap_bleed_ring():
+    """bleed=True: one pixel of half-alpha edge color outside the copper
+    (the mask stops half a cell short of the drawn outline)."""
+    from scipy import ndimage
+    data, disc = _field()
+    rgba = np.asarray(Image.open(io.BytesIO(overlay.heatmap_png(data, 0))))
+    ring = ndimage.binary_dilation(
+        disc, structure=np.ones((3, 3), dtype=bool)) & ~disc
+    assert (rgba[..., 3][ring] == config.OVERLAY_ALPHA // 2).all()
+    outside = ~disc & ~ring
+    assert (rgba[..., 3][outside] == 0).all()
+    assert (rgba[..., 3][disc] == config.OVERLAY_ALPHA).all()
+
+
+def test_heatmap_empty_field():
+    with pytest.raises(ValueError):
+        overlay.heatmap_png(np.full((1, 4, 4), np.nan), 0)

tests/test_skin.py

@@ -64,6 +64,24 @@ def test_parse_frequency():
         skin.parse_frequency("-5k")
 
 
+def test_normalize_decimal():
+    """European decimal commas parse; thousands-separator patterns are
+    rejected ('1,500' silently becoming 1.5 was a 1000x input error)."""
+    assert skin.normalize_decimal("1,5") == "1.5"
+    assert skin.normalize_decimal("0,25") == "0.25"
+    assert skin.normalize_decimal("1,5000") == "1.5000"   # 4 digits: decimal
+    assert skin.normalize_decimal("2.5") == "2.5"
+    for bad in ("1,500", "1.500,5", "1,000,000", "12,345"):
+        with pytest.raises(ValueError, match="separator"):
+            skin.normalize_decimal(bad)
+
+
+def test_parse_frequency_decimal_comma():
+    assert skin.parse_frequency("1,5k") == 1500.0
+    with pytest.raises(ValueError):
+        skin.parse_frequency("1,500")     # ambiguous, not 1.5 Hz
+
+
 def test_single_layer_ac_scales_exactly():
     """Uniform conductance scaling leaves the field shape unchanged:
     R_AC = R_DC * factor to solver precision."""

tests/test_tracks.py

@@ -259,3 +259,25 @@ def test_track_unions_with_fill():
     assert int(s_both.masks.sum()) > int(s_plate.masks.sum())
     assert r_both.R_ohm < 0.75 * r_plate.R_ohm    # bridge shortens the detour
     assert r_both.power_balance_rel < 1e-9
+
+
+def test_pad_copper_bridges_track_junction():
+    """Two traces meet ON an SMD pad, their rounded ends 0.5 mm apart:
+    the junction only exists through the pad copper (board_io stamps
+    the net's pad shapes onto their layers). Without the pad the net
+    is severed - at both track models (rasterized and 1D chain)."""
+    from fill_resistance.errors import ConnectivityError
+    tabs = [[(0, 4.5), (1, 4.5), (1, 5.5), (0, 5.5)],
+            [(19, 4.5), (20, 4.5), (20, 5.5), (19, 5.5)]]
+    pad = [(9.25, 4.4), (10.75, 4.4), (10.75, 5.6), (9.25, 5.6)]
+    segs = [_seg([(0.5, 5), (9.5, 5)], 0.5),
+            _seg([(10.5, 5), (19.5, 5)], 0.5)]
+    r1, r2 = (0, 4.5, 1, 5.5), (19, 4.5, 20, 5.5)
+
+    for h in (0.1, 0.25):          # 5 cells: outlines; 2 cells: 1D chains
+        res, _ = _solve(_seg_problem(segs, r1, r2, fills_mm=tabs + [pad]), h)
+        # ~36 squares of 0.5 mm trace + tabs/pad: sanity-band the value
+        assert 0.007 < res.R_ohm < 0.011
+
+        with pytest.raises(ConnectivityError):
+            _solve(_seg_problem(segs, r1, r2, fills_mm=tabs), h)

tools/deploy.py

@@ -18,7 +18,7 @@ from pathlib import Path
 ROOT = Path(__file__).resolve().parent.parent
 COPY_EXCLUDE = {".venv", ".git", "tests", "tools", "dist", "resources",
                 "__pycache__", ".pytest_cache", "conftest.py", "deploy.ps1",
-                ".gitignore", "metadata.json"}
+                ".gitignore", "metadata.json", "pyproject.toml", "uv.lock"}
 
 
 def plugins_dir(kicad_version: str) -> Path:

tools/gen_readme_figs.py

@@ -0,0 +1,317 @@
+"""Generate the README model figures in docs/img/.
+
+    .venv\\Scripts\\python.exe tools\\gen_readme_figs.py
+
+Real solver output wherever possible: the demo-board maps (raster map
+with the adaptive mesh, current density) and the contact-model
+comparison come straight from the plugin's own pipeline on small
+synthetic boards; only the hole-anatomy cross-section is drawn by hand.
+"""
+import sys
+from pathlib import Path
+
+import numpy as np
+
+ROOT = Path(__file__).resolve().parents[1]
+sys.path.insert(0, str(ROOT))
+
+from fill_resistance import config  # noqa: E402
+
+config.INTERACTIVE = False
+
+import matplotlib  # noqa: E402
+
+matplotlib.use("Agg", force=True)
+
+from fill_resistance import plots, raster, solver  # noqa: E402
+from fill_resistance.geometry import (Electrode, LayerFill, Polygon,  # noqa: E402
+                                      Problem, Rect, ViaLink,
+                                      contact_solder_buildups)
+
+import matplotlib.pyplot as plt  # noqa: E402
+
+plt.switch_backend("Agg")
+plots.INTERACTIVE_BACKEND = None      # save-only: no window panels
+
+NM = 1_000_000
+OUT = ROOT / "docs" / "img"
+
+COPPER = plots._COPPER
+SOLDER = plots._SOLDER
+LEAD = "#5c6570"
+CORE = "#ccd6b3"                      # FR-4
+FILLER = "#eae5dc"                    # non-conductive via fill
+INK = plots._INK
+
+
+def _poly(pts_mm, holes_mm=()) -> Polygon:
+    ring = lambda pts: np.asarray(  # noqa: E731
+        [[int(x * NM), int(y * NM)] for x, y in pts], dtype=np.int64)
+    return Polygon(outline=ring(pts_mm),
+                   holes=[ring(h) for h in holes_mm])
+
+
+def _rect(x0, y0, x1, y1, layer="User.1") -> Rect:
+    return Rect.normalized(int(x0 * NM), int(y0 * NM),
+                           int(x1 * NM), int(y1 * NM), layer)
+
+
+def _disc(x_mm, y_mm, r_mm, n=64) -> Polygon:
+    ang = np.linspace(0, 2 * np.pi, n, endpoint=False)
+    return _poly([(x_mm + r_mm * np.cos(a), y_mm + r_mm * np.sin(a))
+                  for a in ang])
+
+
+def _solve(problem, h_mm, i_test=10.0, model=None):
+    stack = raster.rasterize_stack(problem, int(h_mm * NM))
+    e1, e2 = raster.electrode_masks(stack, problem)
+    p1, p2 = raster.electrode_partition(stack, problem)
+    res = solver.run_solve(problem, stack, e1, e2, i_test,
+                           contact_model=model, parts1=p1, parts2=p2)
+    return res, stack, e1, e2
+
+
+# --- demo board: 2 layers, notched F.Cu pour, half-size B.Cu pour, ----
+# --- a soldered THT-pad contact and a stitching-via field -------------
+
+def demo_problem() -> Problem:
+    # F.Cu: full pour with a notch from the top edge down to y=6 - the
+    # current from the left contact must squeeze through the channel
+    top = _poly([(0, 0), (40, 0), (40, 22), (13, 22), (13, 6),
+                 (10, 6), (10, 22), (0, 22)])
+    # B.Cu: pour on the right half only - vias must carry the transfer
+    bot = _poly([(18, 0), (40, 0), (40, 22), (18, 22)])
+    z_bot = int(1.6 * NM)
+    vias = [ViaLink(x=int(x * NM), y=int(y * NM), drill_nm=300_000,
+                    z_top_nm=-1, z_bot_nm=z_bot + 1, pad_nm=600_000)
+            for x in (20.5, 24.5, 28.5, 32.5, 36) for y in (3, 7, 11, 15, 19)]
+    tht = Electrode(
+        rect=_rect(3.4, 9.9, 5.6, 12.1), contact="F.Cu",
+        label="THT pad (4.5, 11)", drill_nm=1_000_000,
+        pad_nm=int(2.2 * NM), pad_min_nm=int(2.2 * NM),
+        center=(int(4.5 * NM), int(11 * NM)), solder=True,
+        protrusion_side="F.Cu", polygons=[_disc(4.5, 11, 1.1)])
+    lug = Electrode(rect=_rect(37.5, 3, 39.5, 19, "User.2"),
+                    contact="B.Cu", label="lug")
+    p = Problem(
+        board_path="synthetic", net_name="DEMO",
+        rho_ohm_m=1.68e-8, plating_nm=18_000,
+        layers=[LayerFill("F.Cu", 70_000, 0, [top]),
+                LayerFill("B.Cu", 70_000, z_bot, [bot])],
+        vias=vias, electrodes1=[tht], electrodes2=[lug],
+        thickness_source="override")
+    contact_solder_buildups(p)
+    return p
+
+
+def gen_demo_maps():
+    p = demo_problem()
+    res, stack, e1, e2 = _solve(p, 0.05)
+    figs = [
+        (plots.fig_raster(stack, e1, e2, p, res), "demo-raster"),
+        (plots.fig_current(res, stack, e1, e2, p), "demo-current"),
+        (plots.fig_potential(res, stack, e1, e2, p), "demo-potential"),
+    ]
+    plots.save_and_show(figs, OUT, show=False)
+
+
+# --- contact models: equipotential vs uniform injection ---------------
+
+def gen_contact_models():
+    plate = [(0, 0), (24, 0), (24, 18), (0, 18)]
+    r_ohm, zooms = {}, {}
+    # uniform grid: the coarse adaptive leaves would pixelate the |J| zoom
+    config.ADAPTIVE_CELLS = False
+    for model in ("equipotential", "uniform"):
+        p = Problem(
+            board_path="synthetic", net_name="DEMO",
+            rho_ohm_m=1.68e-8, plating_nm=18_000,
+            layers=[LayerFill("F.Cu", 70_000, 0, [_poly(plate)])],
+            vias=[],
+            electrodes1=[Electrode(rect=_rect(4.5, 7.5, 7.5, 10.5))],
+            electrodes2=[Electrode(rect=_rect(22, 1, 23.5, 17, "User.2"))],
+            thickness_source="override")
+        res, stack, _, _ = _solve(p, 0.05, model=model)
+        h_mm = stack.h_nm / NM
+        j = res.Jmag[0] * 1e-6                       # A/mm^2
+        x0, y0 = stack.x0_nm / NM, stack.y0_nm / NM
+        c0, c1 = int((2 - x0) / h_mm), int((13 - x0) / h_mm)
+        r0, r1 = int((3 - y0) / h_mm), int((15 - y0) / h_mm)
+        zooms[model] = (j[r0:r1, c0:c1],
+                        (x0 + c0 * h_mm, x0 + c1 * h_mm,
+                         y0 + r1 * h_mm, y0 + r0 * h_mm))
+        r_ohm[model] = res.R_ohm
+    config.ADAPTIVE_CELLS = True
+
+    vmax = float(np.percentile(
+        zooms["equipotential"][0][np.isfinite(zooms["equipotential"][0])],
+        99.0))
+    fig, axes = plt.subplots(1, 2, figsize=(9.5, 4.2), sharey=True,
+                             layout="constrained")
+    titles = {"equipotential": "equipotential (ideal bonded lug):\n"
+                               "|J| crowds at the contact edges",
+              "uniform": "uniform injection (pressed conductor):\n"
+                         "|J| ramps across the contact"}
+    for ax, model in zip(axes, ("equipotential", "uniform")):
+        data, extent = zooms[model]
+        cmap = matplotlib.colormaps[config.CMAP_CURRENT].copy()
+        cmap.set_bad(plots._BG)
+        im = ax.imshow(data, cmap=cmap, vmin=0, vmax=vmax, origin="upper",
+                       extent=extent, interpolation="nearest")
+        ax.add_patch(plt.Rectangle((4.5, 7.5), 3, 3, fill=False,
+                                   ec="white", ls="--", lw=1.0))
+        ax.set_title(f"{titles[model]}\nR = {r_ohm[model] * 1e3:.3f} mΩ",
+                     fontsize=9, color=INK)
+        ax.set_xlabel("x [mm]", fontsize=8)
+        ax.tick_params(labelsize=8, colors=INK)
+    axes[0].set_ylabel("y [mm]", fontsize=8)
+    cb = fig.colorbar(im, ax=axes, shrink=0.85)
+    cb.set_label("|J| [A/mm²] @ 10 A", fontsize=9)
+    fig.suptitle("The two contact models bracket a real contact:  "
+                 "R$_{equipotential}$ ≤ R$_{real}$ ≤ R$_{uniform}$",
+                 fontsize=10, color=INK)
+    fig.savefig(OUT / "contact-models.png", dpi=config.DPI,
+                facecolor="white", bbox_inches="tight")
+    plt.close(fig)
+    print(f"saved {OUT / 'contact-models.png'}")
+
+
+# --- hole anatomy: hand-drawn cross-section of the four hole types ----
+
+CORE_T = 1.6          # substrate thickness [drawing units ~ mm]
+FOIL_T = 0.18         # foil thickness, exaggerated
+PLATE_W = 0.12        # barrel plating, exaggerated
+CAP_T = 0.07          # via cap
+COAT_T = 0.1          # pad-face solder coat
+Y_TOP = CORE_T + FOIL_T
+
+
+def _board_segment(ax, x0, x1):
+    ax.add_patch(plt.Rectangle((x0, 0), x1 - x0, CORE_T, fc=CORE, ec="none"))
+    for y in (CORE_T, -FOIL_T):
+        ax.add_patch(plt.Rectangle((x0, y), x1 - x0, FOIL_T,
+                                   fc=COPPER, ec="none"))
+
+
+def _barrel(ax, xc, drill):
+    for s in (-1, 1):
+        x = xc + s * drill / 2 - (PLATE_W if s > 0 else 0)
+        ax.add_patch(plt.Rectangle((x, -FOIL_T), PLATE_W,
+                                   CORE_T + 2 * FOIL_T, fc=COPPER,
+                                   ec="none"))
+
+
+def _label(ax, text, xy, xytext, ha="left"):
+    ax.annotate(text, xy, xytext=xytext, fontsize=7.5, color=INK, ha=ha,
+                va="center",
+                arrowprops=dict(arrowstyle="-", color=INK, lw=0.7,
+                                shrinkA=2, shrinkB=1))
+
+
+def gen_hole_anatomy():
+    fig, ax = plt.subplots(figsize=(12.5, 5.2), layout="constrained")
+    holes = [(3.0, 0.7), (10.0, 1.6), (17.5, 1.6), (25.0, 1.6)]
+    edges = [0.0]
+    for xc, d in holes:
+        edges += [xc - d / 2, xc + d / 2]
+    edges.append(28.5)
+    for x0, x1 in zip(edges[::2], edges[1::2]):
+        _board_segment(ax, x0, x1)
+    for xc, d in holes:
+        _barrel(ax, xc, d)
+
+    # 1: small via, filled + capped
+    xc, d = holes[0]
+    ax.add_patch(plt.Rectangle((xc - d / 2 + PLATE_W, -FOIL_T),
+                               d - 2 * PLATE_W, CORE_T + 2 * FOIL_T,
+                               fc=FILLER, ec="none"))
+    for y in (Y_TOP, -FOIL_T - CAP_T):
+        ax.add_patch(plt.Rectangle((xc - d / 2 - 0.12, y), d + 0.24, CAP_T,
+                                   fc=COPPER, ec="none"))
+    _label(ax, "cap, CAP_PLATING_UM (15 µm)\non both outer mouths",
+           (xc, Y_TOP + CAP_T), (xc, 3.3), ha="center")
+    _label(ax, "non-conductive fill", (xc, 0.8), (5.6, -0.9))
+
+    # 2: big via, mouth open
+    xc, d = holes[1]
+    _label(ax, "open mouth: covered cells\nremoved, sub-cell mouths\n"
+               "scale the sheet conductance",
+           (xc, Y_TOP - FOIL_T / 2), (xc, 3.2), ha="center")
+
+    # 3: populated THT pad - full solder joint
+    xc, d = holes[2]
+    lead_w = d - 0.5                      # drill - clearance, exaggerated
+    pad_r = 1.7
+    prot = 1.5
+    sn_edge = "#7d8791"                   # delineate solder sub-shapes
+    # solder fill between plating and lead
+    for s in (-1, 1):
+        x0 = xc + s * lead_w / 2 if s > 0 else xc - d / 2 + PLATE_W
+        ax.add_patch(plt.Rectangle((x0, -FOIL_T),
+                                   d / 2 - PLATE_W - lead_w / 2,
+                                   CORE_T + 2 * FOIL_T, fc=SOLDER,
+                                   ec="none"))
+    # pad-face coat, solder side only
+    ax.add_patch(plt.Rectangle((xc - pad_r, Y_TOP), 2 * pad_r, COAT_T,
+                               fc=SOLDER, ec=sn_edge, lw=0.5))
+    # solder cone: protrusion height at the wall -> 0 at the pad edge
+    for s in (-1, 1):
+        wall = xc + s * lead_w / 2
+        ax.add_patch(plt.Polygon(
+            [(wall, Y_TOP + prot), (wall, Y_TOP + COAT_T),
+             (xc + s * pad_r, Y_TOP + COAT_T)],
+            closed=True, fc=SOLDER, ec=sn_edge, lw=0.5))
+    # lead: through the hole, protruding on top, component below
+    ax.add_patch(plt.Rectangle((xc - lead_w / 2, -2.05), lead_w,
+                               2.05 + Y_TOP + prot, fc=LEAD, ec="none"))
+    ax.add_patch(plt.Rectangle((xc - 1.5, -2.75), 3.0, 0.7,
+                               fc="#8a8f96", ec="none"))
+    ax.text(xc, -2.4, "component", fontsize=7.5, color="white",
+            ha="center", va="center")
+    _label(ax, "clipped lead protrudes\nTHT_LEAD_PROTRUSION_MM (1.5 mm)",
+           (xc + lead_w / 2, Y_TOP + prot - 0.2), (xc + 3.4, 4.15))
+    _label(ax, "solder cone: full height at the\nwall, tapers to 0 at the "
+               "pad edge", (xc - (lead_w / 2 + pad_r) / 2, Y_TOP + 0.7),
+           (13.6, 4.2), ha="center")
+    _label(ax, "pad-face solder coat (50 µm),\nSOLDER side only",
+           (xc - pad_r + 0.2, Y_TOP + COAT_T / 2), (12.9, 2.35),
+           ha="center")
+    _label(ax, "solder-filled hole: lead ∥ solder ∥ plating\n"
+               "lead ⌀ = drill − THT_LEAD_CLEARANCE_MM",
+           (xc - d / 2 + PLATE_W + 0.07, 0.5), (12.3, -1.5), ha="center")
+    _label(ax, "component side:\npad face stays bare",
+           (xc + pad_r - 0.3, -FOIL_T), (xc + 4.0, -1.05))
+
+    # 4: DNP THT pad
+    xc, d = holes[3]
+    _label(ax, "open hole on every layer,\nplating-only barrel, no joint",
+           (xc, 0.8), (xc + 1.3, -2.45), ha="center")
+
+    for (xc, _), title in zip(holes, (
+            "via ≤ cap-drill\n(capped)", "via > cap-drill\n(open)",
+            "THT pad, populated\n(read from KiCad)", "THT pad, DNP")):
+        ax.text(xc, 5.6, title, fontsize=9, color=INK, ha="center",
+                va="top", fontweight="bold")
+
+    handles = [plt.Rectangle((0, 0), 1, 1, fc=c) for c in
+               (COPPER, SOLDER, LEAD, CORE, FILLER)]
+    fig.legend(handles, ("copper (foil / plating / pad)", "solder",
+                         "component lead", "FR-4", "non-conductive fill"),
+               loc="outside right center", fontsize=8, framealpha=0.95)
+    ax.set_title("How drilled holes are modeled — cross-section "
+                 "(vertical scale exaggerated)", fontsize=11, color=INK)
+    ax.set_xlim(-0.3, 29.0)
+    ax.set_ylim(-3.1, 5.7)
+    ax.set_aspect("equal")
+    ax.axis("off")
+    fig.savefig(OUT / "hole-model.png", dpi=config.DPI,
+                facecolor="white", bbox_inches="tight")
+    plt.close(fig)
+    print(f"saved {OUT / 'hole-model.png'}")
+
+
+if __name__ == "__main__":
+    OUT.mkdir(parents=True, exist_ok=True)
+    gen_hole_anatomy()
+    gen_contact_models()
+    gen_demo_maps()

tools/kicad_heatmap_overlay.py

@@ -0,0 +1,141 @@
+"""Standalone runner for the EXPERIMENTAL in-KiCad result overlays
+(also available as a dialog checkbox in the plugin): solve the open
+board headlessly and push per-layer |J| heatmaps as unlocked
+ReferenceImages, transparent outside copper.
+
+    python tools/kicad_heatmap_overlay.py --net VOUT+ --amps 45
+        -> all included copper layers onto config.OVERLAY_LAYERS
+           (User.9..User.12, stackup order, top first)
+    python tools/kicad_heatmap_overlay.py --net X --source B.Cu --dest Eco1.User
+        -> a single layer wherever you want
+
+Needs KiCad >= 10.0.1 with the board open, electrode markers or a
+selection as in a normal plugin run, and the destination layers enabled
+in Board Setup. Re-running replaces the previous overlays. Remove with
+tools/kicad_overlay_test.py --remove --layer <dest>.
+"""
+import argparse
+import sys
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
+
+from kipy.board_types import ReferenceImage
+from kipy.geometry import Vector2
+from kipy.util.board_layer import layer_from_canonical_name
+
+from fill_resistance import config, raster, solver
+from fill_resistance import board_io as bio
+from fill_resistance.overlay import heatmap_png
+
+
+def extract_problem(board, net_arg=None):
+    """Same flow as `python -m fill_resistance.board_io` (dump path)."""
+    # a clicked overlay must not switch the electrode scan into
+    # selection mode - reference images can never be contacts
+    sel = list(board.get_selection())
+    if sel and all(isinstance(s, ReferenceImage) for s in sel):
+        board.clear_selection()
+    stackup = bio.get_stackup_info(board)
+    es1, es2, net_hint = bio.get_electrodes(board, stackup)
+    if bio.any_zone_unfilled(board):
+        bio.refill(board)
+    fills = bio.gather_net_fills(board)
+    tracks = bio.gather_net_tracks(board) if config.INCLUDE_TRACKS else {}
+    copper = bio.merge_copper(fills, bio.tracks_as_polygons(tracks))
+    nets = bio.nets_overlapping(copper, es1, es2)
+    if net_arg:
+        net = net_arg
+    elif net_hint in nets:
+        net = net_hint
+    elif len(nets) == 1:
+        net = nets[0]
+    else:
+        raise SystemExit(f"candidate nets: {nets}; pass one with --net")
+    # marker rectangles may exist for SEVERAL nets (board-wide scan):
+    # keep only the parts overlapping the chosen net's copper
+    per_layer = copper.get(net, {})
+    def on_net(e):
+        return any(bio._rect_overlaps(e.rect, polys)
+                   for polys in per_layer.values())
+    es1, es2 = [e for e in es1 if on_net(e)], [e for e in es2 if on_net(e)]
+    if not es1 or not es2:
+        raise SystemExit(f"no V+/V- marker overlaps {net} copper")
+    print(f"{len(es1)} V+ / {len(es2)} V- marker(s) on {net}")
+    return bio.build_problem(board, net, list(per_layer), es1, es2,
+                             stackup, fills, tracks=tracks)
+
+
+def main() -> None:
+    ap = argparse.ArgumentParser(description=__doc__)
+    ap.add_argument("--source", default=None,
+                    help="single copper layer to overlay (default: ALL "
+                         "included layers onto config.OVERLAY_LAYERS)")
+    ap.add_argument("--dest", default=None,
+                    help="destination layer for --source (default User.9; "
+                         "must be enabled in Board Setup)")
+    ap.add_argument("--net", default=None, help="net name (default: auto)")
+    ap.add_argument("--amps", type=float, default=None,
+                    help="test current [A] (default: config)")
+    ap.add_argument("--lock", action="store_true",
+                    help="lock the overlays (default unlocked: easier to "
+                         "delete; reruns replace them either way)")
+    ap.add_argument("--alpha", type=int, default=None,
+                    help="overlay opacity over copper, 0-255 (default "
+                         "config.OVERLAY_ALPHA)")
+    args = ap.parse_args()
+    if args.alpha is not None:
+        config.OVERLAY_ALPHA = args.alpha
+
+    _, board = bio.connect()
+    problem = extract_problem(board, args.net)
+
+    h = raster.choose_cell_size(problem.copper_bbox(), len(problem.layers))
+    print(f"rasterizing at {h / 1000:.1f} um ...")
+    stack = raster.rasterize_stack(problem, h)
+    # board-wide marker scan: drop parts that land on no copper of THIS
+    # net (markers belonging to other nets' analyses)
+    for name in ("electrodes1", "electrodes2"):
+        parts = getattr(problem, name)
+        keep = [e for e in parts
+                if raster._part_mask3d(stack, problem, e).any()]
+        if len(keep) != len(parts):
+            print(f"ignoring {len(parts) - len(keep)} marker(s) off-net "
+                  f"({name[-1] == '1' and 'V+' or 'V-'})")
+        if not keep:
+            raise SystemExit(f"no {name} marker lands on this net's copper")
+        setattr(problem, name, keep)
+    e1, e2 = raster.electrode_masks(stack, problem)
+    i_test = args.amps if args.amps is not None else config.TEST_CURRENT_A
+    print(f"solving @ {i_test:g} A DC ...")
+    result = solver.run_solve(problem, stack, e1, e2, i_test)
+    print(f"R = {result.R_ohm * 1e3:.4f} mOhm, P = {result.P_total:.3f} W "
+          f"@ {i_test:g} A")
+
+    if args.source is None:
+        bio.push_result_overlays(board, stack, result, lock=args.lock)
+        return
+
+    names = stack.layer_names
+    if args.source not in names:
+        raise SystemExit(f"layer {args.source} not in solve ({names})")
+    png = heatmap_png(result.Jmag * 1e-6, names.index(args.source))
+    ny, nx = stack.shape2d
+    w_nm, h_nm = nx * stack.h_nm, ny * stack.h_nm
+    dest_name = args.dest or "User.9"
+    dest = layer_from_canonical_name(dest_name)
+    n = bio.remove_overlays(board, dest)
+    ref = ReferenceImage()
+    ref.layer = dest
+    ref.position = Vector2.from_xy(round(stack.x0_nm + w_nm / 2),
+                                   round(stack.y0_nm + h_nm / 2))
+    ref.image_scale = w_nm / (nx * bio.OVERLAY_PIX_NM)
+    ref.image_data = png
+    ref.locked = args.lock
+    bio._create_reference_image(board, ref)
+    print(f"{args.source} -> {dest_name} ({nx}x{ny} px, "
+          f"{len(png) / 1024:.0f} kB" + (f", replaced {n}" if n else "") + ")")
+
+
+if __name__ == "__main__":
+    main()

tools/kicad_overlay_test.py

@@ -0,0 +1,156 @@
+"""Route-A experiment: push a bitmap overlay into the open KiCad board as
+a locked ReferenceImage on a User layer via the IPC API.
+
+Pushes a fiducial test pattern (corner + center crosshairs, 10 mm grid,
+translucent gradient) sized to the board outline so alignment and scale
+can be verified by eye in the editor. Re-running replaces the previous
+overlay. Requires KiCad >= 10.0.1 (ReferenceImage over the API).
+
+    python tools/kicad_overlay_test.py [--layer Cmts.User] [--remove]
+    python tools/kicad_overlay_test.py --image heat.png --bbox x0,y0,x1,y1
+                                       (mm; push an arbitrary PNG instead)
+
+The overlay is editor-only: reference images never plot to gerbers.
+Delete it any time by selecting it in KiCad (it sits on the chosen
+layer) or with --remove. The layer must be enabled in Board Setup:
+User.1..User.45 usually are NOT (KiCad refuses the item with 'no
+overlapping layers with the board'); Cmts.User/Eco1.User always exist.
+"""
+import argparse
+import io
+import sys
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
+
+from kipy.board_types import ReferenceImage
+from kipy.geometry import Vector2
+from kipy.util.board_layer import canonical_name, layer_from_canonical_name
+
+from fill_resistance.board_io import (OVERLAY_PIX_NM as PIX_NM,
+                                      _create_reference_image, connect,
+                                      remove_overlays)
+
+NM = 1_000_000
+
+
+def board_bbox_nm(board):
+    """Union bbox of the Edge.Cuts shapes (fallback: all pads)."""
+    items = [s for s in board.get_shapes()
+             if canonical_name(s.layer) == "Edge.Cuts"]
+    if not items:
+        items = list(board.get_pads())
+    if not items:
+        raise SystemExit("board has no Edge.Cuts shapes and no pads")
+    x0 = y0 = None
+    x1 = y1 = None
+    for it in items:
+        box = board.get_item_bounding_box(it)
+        if box is None:
+            continue
+        lo_x, lo_y = box.pos.x, box.pos.y
+        hi_x, hi_y = lo_x + box.size.x, lo_y + box.size.y
+        x0 = lo_x if x0 is None else min(x0, lo_x)
+        y0 = lo_y if y0 is None else min(y0, lo_y)
+        x1 = hi_x if x1 is None else max(x1, hi_x)
+        y1 = hi_y if y1 is None else max(y1, hi_y)
+    return x0, y0, x1, y1
+
+
+def fiducial_png(w_nm: float, h_nm: float, px_per_mm: float = 16.0):
+    """RGBA test pattern: translucent gradient, 10 mm grid, opaque
+    crosshairs at the four corners and the center."""
+    import numpy as np
+    from PIL import Image
+
+    w_px = max(2, round(w_nm / NM * px_per_mm))
+    h_px = max(2, round(h_nm / NM * px_per_mm))
+    xx = np.linspace(0.0, 1.0, w_px)[None, :]
+    yy = np.linspace(0.0, 1.0, h_px)[:, None]
+    rgba = np.zeros((h_px, w_px, 4), dtype=np.uint8)
+    rgba[..., 0] = (255 * xx).astype(np.uint8)          # red ramp ->
+    rgba[..., 2] = (255 * yy).astype(np.uint8)          # blue ramp v
+    rgba[..., 1] = 60
+    rgba[..., 3] = 70                                   # mostly see-through
+
+    step = round(10.0 * px_per_mm)                      # 10 mm grid
+    for x in range(0, w_px, step):
+        rgba[:, x:x + 2, :3] = 255
+        rgba[:, x:x + 2, 3] = 150
+    for y in range(0, h_px, step):
+        rgba[y:y + 2, :, :3] = 255
+        rgba[y:y + 2, :, 3] = 150
+
+    def cross(cx, cy, arm=round(3 * px_per_mm)):
+        x_lo, x_hi = max(0, cx - arm), min(w_px, cx + arm + 1)
+        y_lo, y_hi = max(0, cy - arm), min(h_px, cy + arm + 1)
+        cy2 = np.clip(cy, 0, h_px - 2)
+        cx2 = np.clip(cx, 0, w_px - 2)
+        rgba[cy2:cy2 + 2, x_lo:x_hi] = (255, 0, 0, 255)
+        rgba[y_lo:y_hi, cx2:cx2 + 2] = (255, 0, 0, 255)
+
+    for cx in (0, w_px - 1):
+        for cy in (0, h_px - 1):
+            cross(cx, cy)
+    cross(w_px // 2, h_px // 2)
+
+    buf = io.BytesIO()
+    # no dpi= : without a density chunk KiCad assumes the 300 PPI default
+    Image.fromarray(rgba, "RGBA").save(buf, format="PNG")
+    return buf.getvalue(), w_px, h_px
+
+
+def main() -> None:
+    ap = argparse.ArgumentParser(description=__doc__)
+    ap.add_argument("--layer", default="Cmts.User",
+                    help="destination layer (default Cmts.User; must be "
+                         "enabled in Board Setup)")
+    ap.add_argument("--remove", action="store_true",
+                    help="only remove existing overlays on the layer")
+    ap.add_argument("--image", help="push this PNG instead of the pattern")
+    ap.add_argument("--bbox", help="x0,y0,x1,y1 [mm] for --image")
+    args = ap.parse_args()
+
+    _, board = connect()
+    layer = layer_from_canonical_name(args.layer)
+
+    n = remove_overlays(board, layer)
+    if n:
+        print(f"removed {n} previous overlay(s) on {args.layer}")
+    if args.remove:
+        return
+
+    if args.image:
+        if not args.bbox:
+            raise SystemExit("--image needs --bbox x0,y0,x1,y1 [mm]")
+        x0, y0, x1, y1 = (float(v) * NM for v in args.bbox.split(","))
+        png = Path(args.image).read_bytes()
+        from PIL import Image
+        w_px, h_px = Image.open(io.BytesIO(png)).size
+    else:
+        x0, y0, x1, y1 = board_bbox_nm(board)
+        png, w_px, h_px = fiducial_png(x1 - x0, y1 - y0)
+
+    scale = (x1 - x0) / (w_px * PIX_NM)
+
+    ref = ReferenceImage()
+    ref.layer = layer
+    ref.position = Vector2.from_xy(round((x0 + x1) / 2), round((y0 + y1) / 2))
+    ref.image_scale = scale
+    ref.image_data = png
+    ref.locked = False          # unlocked: easy to delete; reruns replace
+    _create_reference_image(board, ref)
+
+    got = [r for r in board.get_reference_images() if r.layer == layer]
+    print(f"pushed {len(png) / 1024:.0f} kB PNG ({w_px}x{h_px} px) onto "
+          f"{args.layer}: {(x1 - x0) / NM:.2f} x {(y1 - y0) / NM:.2f} mm at "
+          f"({x0 / NM:.2f}, {y0 / NM:.2f}) mm, scale {scale:.4f}")
+    for r in got:
+        print(f"readback: {r!r}")
+    print(f"-> enable layer '{args.layer}' in the Appearance panel; the "
+          f"red crosshairs must sit on the board bbox corners/center and "
+          f"the white grid must be 10 mm. Remove with --remove.")
+
+
+if __name__ == "__main__":
+    main()