Barrel contacts for vias/THT pads; configurable cap-drill threshold

- Selected vias and through-hole pads inject at the drill-wall ring on
  every spanned layer (both contact models), not the whole pad face,
  so the pad/pour spreading resistance is part of the result. Vias are
  now selectable as contacts.
- Soldered THT joints: the hole is modeled solder-filled (core in
  parallel with the plating, also for stitching THT barrels) and the
  pad face carries an average-thickness solder coat over the modeled
  copper (SOLDER_THICKNESS_UM).
- Vias with drills above a configurable threshold (dialog field,
  default CAP_MAX_DRILL_MM = 0.5) keep open mouths even with capping
  selected - the fab caps only small vias.
- Geometry dump schema v6: electrode barrel fields, cap_max_drill_nm.
- Verified against R = rho/(pi t)*acosh(d/2a) for two circular contacts
  on a sheet (+2.6% at h = 0.15 mm, a = 1 mm; uniform model above the
  equipotential one as required by the contact bracket).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

README.md

@@ -44,9 +44,10 @@ SWIG API. Requires KiCad **10.0.1+**.
    - **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,
-     SMD pad its own layer) — selected pads fill a side that has no
-     rectangles;
+   - **pads and vias** (SMD pad: real copper shape on its own layer;
+     through-hole pads and vias become **barrel contacts** — the current
+     enters at the drill wall on every spanned layer, see below) —
+     selected pads/vias fill a side that has no rectangles;
    - 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.
@@ -76,14 +77,19 @@ 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
+  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) — 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
-  remain outside the model; at f > 0 the thickness scaling is applied
-  multiplicatively to the skin-corrected sheet conductance
+  remain outside the model, but THT-pad **barrels are solder-filled**
+  (a soldered component lead): the solder core (SAC305) conducts in
+  parallel with the plating annulus. 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
@@ -110,6 +116,18 @@ 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). A
+  soldered **THT joint** additionally assumes the **hole is filled with
+  solder** (core in parallel with the plating) and the **pad face
+  carries an average-thickness solder coat** (`SOLDER_THICKNESS_UM`,
+  50 µm) over the modeled copper under the pad shape. 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

fill_resistance/board_io.py

@@ -11,7 +11,7 @@ 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 +22,8 @@ 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)
 
 MASK_TO_COPPER = {"F.Mask": "F.Cu", "B.Mask": "B.Cu"}
 
@@ -179,7 +180,8 @@ def _pad_polygons(board: Board, pad: Pad, contact: str) -> list[Polygon] | None:
     return None
 
 
-def _to_electrode(board: Board, item) -> Electrode:
+def _to_electrode(board: Board, item,
+                  stackup: StackupInfo | 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 +190,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 +215,48 @@ 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 = _pad_drill_nm(pad)
     return Electrode(rect=rect, contact=contact,
-                     polygons=_pad_polygons(board, pad, contact), label=label)
+                     polygons=_pad_polygons(board, pad, contact), label=label,
+                     # through-hole pad: current enters at the soldered
+                     # barrel; the joint is solder-filled + pad-coated
+                     drill_nm=drill, pad_nm=_padstack_pad_nm(pad),
+                     center=(pad.position.x, pad.position.y),
+                     solder=drill > 0)
 
 
-def _net_hint_of(pads: list[Pad]) -> str | None:
-    for pad in pads:
-        if pad.net is not None:
-            return pad.net.name
+def _net_hint_of(items: list) -> str | None:
+    for item in items:
+        if item.net is not None:
+            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))]
 
     if not selection:
         allr = [s for s in board.get_shapes() if isinstance(s, BoardRectangle)]
@@ -258,12 +287,12 @@ 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"layers already provide rectangles. Use pads only for a "
-                f"side that has none."
+                f"Cannot assign the {len(pads)} selected pad(s)/via(s): both "
+                f"marker layers already provide rectangles. Use pads/vias "
+                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) for p in pads]
             if not es1:
                 es1 = pad_parts
             else:
@@ -277,12 +306,12 @@ def get_electrodes(board: Board
 
     items = rects + pads
     if len(items) == 2:
-        return ([_to_electrode(board, items[0])],
-                [_to_electrode(board, items[1])], _net_hint_of(pads))
+        return ([_to_electrode(board, items[0], stackup)],
+                [_to_electrode(board, items[1], stackup)], _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"contacts are needed.\n{scheme}"
+        f"layers) and {len(pads)} pad(s)/via(s); without marker layers "
+        f"exactly 2 contacts are needed.\n{scheme}"
     )
 
 
@@ -465,7 +494,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 = []
@@ -502,7 +532,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 +552,15 @@ 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),
     )
+    solder_layers = contact_solder_buildups(problem)
+    if solder_layers:
+        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)})")
+    return problem
 
 
 if __name__ == "__main__":
@@ -533,7 +571,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,13 @@ 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_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
 SKIN_SIDES = 1                  # skin-effect field config: 1 = plane facing a
                                 # return plane (conservative), 2 = isolated foil
 

fill_resistance/dialog.py

@@ -37,6 +37,7 @@ class Selection:
     extra_cu_um: float = 0.0
     include_tracks: bool = True
     vias_capped: bool = True
+    cap_max_drill_mm: float = 0.5
     adaptive: bool = False
 
 
@@ -73,6 +74,11 @@ 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)")
@@ -205,6 +211,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,6 +233,7 @@ class _Dialog(QDialog):
                          extra_cu_um=extra_cu,
                          include_tracks=self.tracks_check.isChecked(),
                          vias_capped=self.capped_check.isChecked(),
+                         cap_max_drill_mm=cap_max_drill,
                          adaptive=self.adaptive_check.isChecked())
 
     def _try_accept(self) -> None:

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,29 @@ class TrackSeg:
 
 @dataclass
 class Electrode:
-    """One PART of a current-injection terminal: a drawn rectangle or a
-    selected pad. A terminal (V+ or V-) is a LIST of parts, all merged
-    into one equipotential contact (externally bonded). `polygons`
-    (board nm) is the exact copper shape when known (pads); None means
-    the rectangle itself is the shape. `contact` = 'all' or a layer
-    name: which included layers this part touches."""
+    """One PART of a current-injection terminal: a drawn rectangle, a
+    selected pad, or a selected via. A terminal (V+ or V-) is a LIST of
+    parts, all merged into one equipotential contact (externally
+    bonded). `polygons` (board nm) is the exact copper shape when known
+    (pads); None means the rectangle itself is the shape. `contact` =
+    'all' or a layer name: which included layers this part touches.
+
+    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 carries an average-thickness solder coat
+    (Problem.solder_thickness_nm over `polygons`)."""
     rect: Rect                                # bounding box (labels/summary)
     contact: str = "all"
     polygons: list[Polygon] | None = None
     label: str = "rect"
+    drill_nm: int = 0                         # >0: barrel contact
+    pad_nm: int = 0                           # pad diameter (search bound)
+    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)
 
 
 @dataclass
@@ -121,11 +134,18 @@ class ViaLink:
         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) -> float:
         """Barrel segment resistance over length_nm: thin-wall annulus of
-        plating around the drill."""
-        area_m2 = math.pi * (self.drill_nm * 1e-9) * (plating_nm * 1e-9)
-        return rho_ohm_m * (length_nm * 1e-9) / area_m2
+        plating around the drill. With solder_rho_ohm_m the hole is
+        solder-filled (soldered THT joint): the solder core conducts in
+        parallel with the plating."""
+        ga = math.pi * (self.drill_nm * 1e-9) * (plating_nm * 1e-9) \
+            / rho_ohm_m                       # plating 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
+            ga += math.pi * r_core * r_core / solder_rho_ohm_m
+        return (length_nm * 1e-9) / ga
 
 
 @dataclass
@@ -147,6 +167,9 @@ 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
 
     @property
     def layer_names(self) -> list[str]:
@@ -173,6 +196,25 @@ 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 is covered in solder of
+    average thickness solder_thickness_nm. Adds one SurfaceBuildup per
+    outer layer 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}
+    outer = [n for n in ("F.Cu", "B.Cu") if n in included]
+    touched = []
+    for e in problem.electrodes1 + problem.electrodes2:
+        if not e.solder or not e.polygons:
+            continue
+        for name in outer:
+            problem.buildups.append(
+                SurfaceBuildup(layer_name=name, polygons=list(e.polygons)))
+            touched.append(name)
+    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 +363,11 @@ 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,
+        "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,
     }
 
 
@@ -331,6 +378,13 @@ 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)),
+        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)),
     )
 
 
@@ -369,6 +423,7 @@ 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,
     }
 
 
@@ -442,6 +497,7 @@ 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)),
     )
 
 

fill_resistance/main.py

@@ -33,7 +33,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,7 +95,8 @@ 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}")

fill_resistance/raster.py

@@ -268,8 +268,10 @@ 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
-    keep a thickness-scaled sheet conductance via stack.thick_scale."""
+    hole. The fab caps only small vias: drills above cap_max_drill_nm
+    stay open even with vias_capped. 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)
@@ -296,7 +298,8 @@ def _apply_via_mouths(stack: RasterStack, problem: Problem) -> None:
         if stack.thick_scale is None:
             stack.thick_scale = np.ones(stack.masks.shape)
         for li in _via_span(problem, via):
-            if problem.vias_capped and li in outer:
+            if 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:
@@ -399,26 +402,88 @@ 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),
+    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
+    out = np.zeros((ny, nx), dtype=bool)
+    j0 = max(0, math.floor((x - rw - stack.x0_nm) / h))
+    j1 = min(nx, math.floor((x + rw - stack.x0_nm) / h) + 1)
+    i0 = max(0, math.floor((y - rw - stack.y0_nm) / h))
+    i1 = min(ny, math.floor((y + rw - 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 = np.sqrt(ys[:, None] ** 2 + xs[None, :] ** 2)
+    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
-    lug / through 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."""
+    part's contact layer(s), or the barrel-wall ring for via/THT-pad
+    contacts (current enters through the soldered barrel, not the pad
+    face). contact 'all' = every included layer (bolted lug / through
+    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 = 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]
+            part = _part_mask3d(stack, problem, el)
             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 +511,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 = 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 = _part_mask3d(stack, problem, el)
             m &= ~claimed
             claimed |= m
             out.append((el.label, m))

fill_resistance/solver.py

@@ -178,8 +178,12 @@ 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,
-                                         problem.plating_nm)
+            # THT pads carry a soldered component lead: the hole is
+            # solder-filled, the core conducts in parallel with the plating
+            r_dc = via.barrel_resistance(
+                length, problem.rho_ohm_m, problem.plating_nm,
+                solder_rho_ohm_m=(problem.solder_rho_ohm_m
+                                  if via.kind == "pad" else None))
             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:

tests/test_barrel_contacts.py

@@ -0,0 +1,158 @@
+"""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,
+                                      save_problem)
+from tests.util import NM, 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 the outer layers, 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)])]
+        return p
+
+    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_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)
+    f = tmp_path / "d.json"
+    save_problem(p, f)
+    e = load_problem(f).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

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,11 @@ 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")
+    # 'pad' barrels are solder-filled; kill the core so the reference
+    # barrel matches the via's plating-only resistance exactly
+    ref.solder_rho_ohm_m = 1e30
+    r_ref, _ = _solve(ref, 0.1)
     assert r_cap.R_ohm == pytest.approx(r_ref.R_ohm, rel=1e-9)
 
 
@@ -90,14 +98,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)