Adom Fields

//! Solid bodies off the board (a package, an interface pad, a finned heatsink) solved in the same steady heat system as
//! the board copper, so the heat runs from the die through the tab, the pad and the sink to the air, and down the leads
//! into the board, in one solve with one energy balance. Each body is a union of axis-aligned boxes in board millimetres
//! (x and y as KiCad, z up from the top copper face; a body under the board has z below minus the board thickness). The
//! grid keeps every box edge and subdivides to the pitch, so a 0.6 mm pad or a 1.2 mm fin is never rounded away.
//! Bodies touch through declared contacts only; exposed faces lose heat to the air through a prescribed h (not CFD, no
//! radiation); a face lying on the board's outer copper couples to the copper under it; point links (a lead to its pad)
//! join a body to one board cell.
use crate::solve::System;
use serde::{Deserialize, Serialize};
use serde_json::{json, Value};

#[derive(Clone, Deserialize, Serialize)]
#[serde(rename_all = "camelCase", deny_unknown_fields)]
pub struct Body {
    pub name: String,
    /// [x0, y0, z0, x1, y1, z1] in board mm, z up from the top copper face
    pub boxes_mm: Vec<[f64; 6]>,
    pub conductivity_w_mk: f64,
    /// heat generated in the body, spread over its volume
    #[serde(default)] pub watts: Option<f64>,
    /// the hot part (spec `hot`) whose watts this body takes; the board then gets none of that part's heat directly
    #[serde(default)] pub heat_of: Option<String>,
    /// convection on exposed faces, W/m²K, still and moving air (prescribed, e.g. calibrated to a vendor's sink rating)
    pub h_still_w_m2k: f64,
    pub h_moving_w_m2k: f64,
    /// a face lying on the board's outer copper (z = 0 on top, z = -thickness under it, within 0.05 mm) couples to it through this, m²K/W
    #[serde(default)] pub board_contact_m2k_w: Option<f64>,
    #[serde(default)] pub provenance: String,
}
#[derive(Clone, Deserialize, Serialize)]
#[serde(rename_all = "camelCase", deny_unknown_fields)]
pub struct Contact { pub bodies: [String; 2], pub resistance_m2k_w: f64, #[serde(default)] pub provenance: String }
#[derive(Clone, Deserialize, Serialize)]
#[serde(rename_all = "camelCase", deny_unknown_fields)]
pub struct BoardLink { pub body: String, pub at: [f64; 2], pub layer: String, pub resistance_k_per_w: f64, #[serde(default)] pub pad: String }
#[derive(Clone, Deserialize, Serialize)]
#[serde(rename_all = "camelCase", deny_unknown_fields)]
pub struct Solids {
    #[serde(default = "pitch")] pub pitch_mm: f64,
    pub bodies: Vec<Body>,
    #[serde(default)] pub contacts: Vec<Contact>,
    #[serde(default)] pub board_links: Vec<BoardLink>,
    #[serde(default)] pub assumptions: Vec<String>,
}
fn pitch() -> f64 { 1.0 }

/// the voxel grid: breakpoints per axis and the owning body of each cell
pub struct Grid { pub xs: Vec<f64>, pub ys: Vec<f64>, pub zs: Vec<f64>, pub owner: Vec<Option<usize>> }
impl Grid {
    pub fn dims(&self) -> (usize, usize, usize) { (self.xs.len() - 1, self.ys.len() - 1, self.zs.len() - 1) }
    pub fn idx(&self, i: usize, j: usize, k: usize) -> usize { let (nx, ny, _) = self.dims(); (k * ny + j) * nx + i }
}
fn axis(mut cuts: Vec<f64>, pitch: f64) -> Vec<f64> {
    cuts.sort_by(|a, b| a.partial_cmp(b).unwrap()); cuts.dedup_by(|a, b| (*a - *b).abs() < 1e-6);
    let mut out = vec![cuts[0]];
    for w in cuts.windows(2) { let n = ((w[1] - w[0]) / pitch).ceil().max(1.0) as usize; for s in 1..=n { out.push(w[0] + (w[1] - w[0]) * s as f64 / n as f64); } }
    out
}

impl Solids {
    pub fn validate(&self, board: [f64; 4], layers_ok: &dyn Fn(&str) -> bool, hot: &dyn Fn(&str) -> bool) -> Result<Grid, String> {
        if self.bodies.is_empty() { return Err("solids: no bodies".into()); }
        if !(self.pitch_mm.is_finite() && self.pitch_mm >= 0.1) { return Err("solids.pitchMm must be at least 0.1 mm".into()); }
        let _ = board;
        for (i, b) in self.bodies.iter().enumerate() {
            if b.name.trim().is_empty() || self.bodies[..i].iter().any(|o| o.name == b.name) { return Err(format!("solids: body {i} needs a unique name")); }
            if b.boxes_mm.is_empty() || b.boxes_mm.iter().any(|x| x.iter().any(|v| !v.is_finite()) || x[3] <= x[0] || x[4] <= x[1] || x[5] <= x[2]) { return Err(format!("solids: {}: every box is [x0, y0, z0, x1, y1, z1] with x1 > x0, y1 > y0, z1 > z0", b.name)); }
            if !(b.conductivity_w_mk.is_finite() && b.conductivity_w_mk > 0.0) || !(b.h_still_w_m2k.is_finite() && b.h_still_w_m2k >= 0.0) || !(b.h_moving_w_m2k.is_finite() && b.h_moving_w_m2k >= 0.0) { return Err(format!("solids: {}: conductivityWMk must be positive and hStillWM2k / hMovingWM2k not negative", b.name)); }
            if b.watts.is_some_and(|w| !w.is_finite() || w < 0.0) { return Err(format!("solids: {}: watts must be finite and not negative", b.name)); }
            if b.board_contact_m2k_w.is_some_and(|r| !r.is_finite() || r < 0.0) { return Err(format!("solids: {}: boardContactM2kW must be finite and not negative", b.name)); }
            if let Some(h) = &b.heat_of { if !hot(h) { return Err(format!("solids: {}: heatOf {h} is not in the spec's hot parts", b.name)); } if b.watts.is_some() { return Err(format!("solids: {}: give watts or heatOf, not both (the heat would count twice)", b.name)); } }
        }
        for c in &self.contacts {
            for n in &c.bodies { if !self.bodies.iter().any(|b| b.name == *n) { return Err(format!("solids: contact names {n}, which is not a body")); } }
            if c.bodies[0] == c.bodies[1] || !(c.resistance_m2k_w.is_finite() && c.resistance_m2k_w >= 0.0) { return Err(format!("solids: contact {:?} needs two different bodies and a resistance of zero or more", c.bodies)); }
        }
        for l in &self.board_links {
            if !self.bodies.iter().any(|b| b.name == l.body) { return Err(format!("solids: board link names {}, which is not a body", l.body)); }
            if !layers_ok(&l.layer) || !(l.resistance_k_per_w.is_finite() && l.resistance_k_per_w > 0.0) || l.at[0] < board[0] || l.at[0] > board[2] || l.at[1] < board[1] || l.at[1] > board[3] { return Err(format!("solids: board link from {} needs a copper layer of the board, a point on the board and a positive resistanceKPerW", l.body)); }
        }
        let all = || self.bodies.iter().flat_map(|b| b.boxes_mm.iter());
        let xs = axis(all().flat_map(|x| [x[0], x[3]]).collect(), self.pitch_mm);
        let ys = axis(all().flat_map(|x| [x[1], x[4]]).collect(), self.pitch_mm);
        let zs = axis(all().flat_map(|x| [x[2], x[5]]).collect(), self.pitch_mm);
        let n = (xs.len() - 1) * (ys.len() - 1) * (zs.len() - 1);
        if n > 3_000_000 { return Err(format!("solids: the grid would have {n} cells; raise pitchMm")); }
        let mut g = Grid { xs, ys, zs, owner: vec![None; n] };
        let (nx, ny, nz) = g.dims();
        for k in 0..nz { for j in 0..ny { for i in 0..nx {
            let c = [(g.xs[i] + g.xs[i + 1]) / 2.0, (g.ys[j] + g.ys[j + 1]) / 2.0, (g.zs[k] + g.zs[k + 1]) / 2.0];
            let mut who: Option<usize> = None;
            for (bi, b) in self.bodies.iter().enumerate() { if b.boxes_mm.iter().any(|x| c[0] > x[0] && c[0] < x[3] && c[1] > x[1] && c[1] < x[4] && c[2] > x[2] && c[2] < x[5]) {
                if let Some(o) = who { if o != bi { return Err(format!("solids: {} and {} overlap at ({:.2}, {:.2}, {:.2}); bodies may touch but not overlap", self.bodies[o].name, b.name, c[0], c[1], c[2])); } } who = Some(bi);
            } }
            let id = g.idx(i, j, k); g.owner[id] = who;
        } } }
        for (bi, b) in self.bodies.iter().enumerate() { if !g.owner.iter().any(|o| *o == Some(bi)) { return Err(format!("solids: {} has no cell on the grid", b.name)); } }
        Ok(g)
    }
    fn contact(&self, a: usize, b: usize) -> Option<f64> { let (na, nb) = (&self.bodies[a].name, &self.bodies[b].name); self.contacts.iter().find(|c| (c.bodies[0] == *na && c.bodies[1] == *nb) || (c.bodies[0] == *nb && c.bodies[1] == *na)).map(|c| c.resistance_m2k_w * 1e6) }
}

/// what the solve needs to read the bodies back: the node of each owned cell and every conductance by kind
pub struct Attached {
    pub cells: Vec<(usize, usize, [f64; 6])>, // (node, body, box)
    pub faces: Vec<(usize, usize, f64, usize, f64)>, // (cell a, cell b, conductance W/K, axis, face area mm²) between owned cells
    pub air: Vec<(usize, f64)>,                     // (cell, conductance to ambient)
    pub board: Vec<(usize, usize, f64)>,            // (board node, cell, conductance)
    pub sources: Vec<(usize, f64)>,                 // (node, watts)
    pub on_copper: Vec<bool>,                       // per body: a face of it couples to the board's outer copper
}

#[allow(clippy::too_many_arguments)]
pub fn attach(sys: &mut System, s: &Solids, g: &Grid, level: f64, fan_n: f64, x0: f64, y0: f64, d: f64, h_air: f64, thickness: f64, under_face: (f64, f64), hot_watts: &dyn Fn(&str) -> f64) -> Result<Attached, String> {
    let (nx, ny, nz) = g.dims(); let (w, h) = (sys.w, sys.h); let nb = sys.nl * w * h;
    let mut cell_of = vec![usize::MAX; g.owner.len()]; let mut a = Attached { cells: Vec::new(), faces: Vec::new(), air: Vec::new(), board: Vec::new(), sources: Vec::new(), on_copper: vec![false; s.bodies.len()] };
    let kmm = |b: usize| s.bodies[b].conductivity_w_mk / 1000.0;
    let hmm = |b: usize| crate::fields::fan_mix(s.bodies[b].h_still_w_m2k, s.bodies[b].h_moving_w_m2k, level, fan_n) * 1e-6;
    let span = |v: &Vec<f64>, i: usize| v[i + 1] - v[i];
    for k in 0..nz { for j in 0..ny { for i in 0..nx { let id = g.idx(i, j, k); if let Some(b) = g.owner[id] {
        let node = sys.add_node(0.0); cell_of[id] = a.cells.len(); a.cells.push((node, b, [g.xs[i], g.ys[j], g.zs[k], g.xs[i + 1], g.ys[j + 1], g.zs[k + 1]]));
    } } } }
    let mut leak = vec![0.0f64; a.cells.len()];
    let board_cells = |bx: [f64; 4]| -> Vec<(usize, f64)> {
        let mut v = Vec::new();
        let (c0, c1) = (((bx[0] - x0) / d).floor().max(0.0) as usize, (((bx[2] - x0) / d).ceil().max(0.0) as usize).min(w));
        let (r0, r1) = (((bx[1] - y0) / d).floor().max(0.0) as usize, (((bx[3] - y0) / d).ceil().max(0.0) as usize).min(h));
        for r in r0..r1 { for c in c0..c1 { let (xa, xb) = ((x0 + c as f64 * d).max(bx[0]), (x0 + (c + 1) as f64 * d).min(bx[2])); let (ya, yb) = ((y0 + r as f64 * d).max(bx[1]), (y0 + (r + 1) as f64 * d).min(bx[3])); if xb > xa && yb > ya { v.push((r * w + c, (xb - xa) * (yb - ya))); } } }
        v
    };
    for k in 0..nz { for j in 0..ny { for i in 0..nx {
        let id = g.idx(i, j, k); let Some(b) = g.owner[id] else { continue }; let ci = cell_of[id];
        let (dx, dy, dz) = (span(&g.xs, i), span(&g.ys, j), span(&g.zs, k));
        // the three positive-side neighbours: conduction or a declared contact; the exposed faces lose heat to the air
        for (ax, nb_ijk, area, half, half_n) in [
            (0usize, (i + 1 < nx).then(|| g.idx(i + 1, j, k)), dy * dz, dx / 2.0, if i + 1 < nx { span(&g.xs, i + 1) / 2.0 } else { 0.0 }),
            (1, (j + 1 < ny).then(|| g.idx(i, j + 1, k)), dx * dz, dy / 2.0, if j + 1 < ny { span(&g.ys, j + 1) / 2.0 } else { 0.0 }),
            (2, (k + 1 < nz).then(|| g.idx(i, j, k + 1)), dx * dy, dz / 2.0, if k + 1 < nz { span(&g.zs, k + 1) / 2.0 } else { 0.0 }),
        ] {
            match nb_ijk.and_then(|q| g.owner[q].map(|o| (q, o))) {
                Some((q, o)) => {
                    let rc = if o == b { 0.0 } else { s.contact(b, o).ok_or_else(|| format!("solids: {} and {} touch but have no contact; declare one (resistanceM2kW 0 means ideal)", s.bodies[b].name, s.bodies[o].name))? / area };
                    let gc = 1.0 / (half / (kmm(b) * area) + half_n / (kmm(o) * area) + rc);
                    let cj = cell_of[q]; sys.edges.push((a.cells[ci].0, a.cells[cj].0, gc)); a.faces.push((ci, cj, gc, ax, area));
                }
                None => { let hh = hmm(b); if hh > 0.0 { leak[ci] += area / (1.0 / hh + half / kmm(b)); } }
            }
        }
        for (ax, prev, area, half) in [(0usize, (i > 0).then(|| g.idx(i - 1, j, k)), dy * dz, dx / 2.0), (1, (j > 0).then(|| g.idx(i, j - 1, k)), dx * dz, dy / 2.0), (2, (k > 0).then(|| g.idx(i, j, k - 1)), dx * dy, dz / 2.0)] {
            if prev.and_then(|q| g.owner[q]).is_some() { continue; }
            let z_lo = g.zs[k]; let z_hi = g.zs[k + 1];
            // a face on the board's outer copper: the top face of a body under the board, or the bottom face of one on it
            let on_board = ax == 2 && s.bodies[b].board_contact_m2k_w.is_some() && (z_lo.abs() < 1e-6);
            if on_board {
                let r = s.bodies[b].board_contact_m2k_w.unwrap() * 1e6;
                for (bc, ar) in board_cells([g.xs[i], g.ys[j], g.xs[i + 1], g.ys[j + 1]]) { let gc = ar / (r + half / kmm(b)); sys.edges.push((bc, a.cells[ci].0, gc)); a.board.push((bc, ci, gc)); a.on_copper[b] = true; sys.leak[0][bc] = (sys.leak[0][bc] - h_air * ar).max(1e-12); }
                continue;
            }
            let _ = z_hi;
            let hh = hmm(b); if hh > 0.0 { leak[ci] += area / (1.0 / hh + half / kmm(b)); }
        }
        // the top face of a body lying under the board touches B.Cu
        let top_under = -g.zs[k + 1];
        let cover = if s.bodies[b].board_contact_m2k_w.is_some() && top_under >= under_face.0 && top_under <= under_face.1 && (k + 1 >= nz || g.owner[g.idx(i, j, k + 1)].is_none()) { board_cells([g.xs[i], g.ys[j], g.xs[i + 1], g.ys[j + 1]]) } else { Vec::new() };
        if !cover.is_empty() {
            let r = s.bodies[b].board_contact_m2k_w.unwrap() * 1e6; let bl = sys.nl - 1;
            // that face was counted as exposed above; take its air loss back (a cell past the board's edge stays in the air)
            let hh = hmm(b); if hh > 0.0 { leak[ci] -= dx * dy / (1.0 / hh + dz / 2.0 / kmm(b)); }
            a.on_copper[b] = true;
            for (bc, ar) in cover { let node = bl * w * h + bc; let gc = ar / (r + dz / 2.0 / kmm(b)); sys.edges.push((node, a.cells[ci].0, gc)); a.board.push((node, ci, gc)); sys.leak[bl][bc] = (sys.leak[bl][bc] - h_air * ar).max(1e-12); }
        }
    } } }
    for (ci, l) in leak.iter().enumerate() { if *l > 0.0 { let node = a.cells[ci].0; sys.extra_leak[node - nb] += l; a.air.push((ci, *l)); } }
    // point links: from the body's cell nearest the board point to that board cell
    for l in &s.board_links {
        let bi = s.bodies.iter().position(|b| b.name == l.body).unwrap();
        let li = if l.layer == "F.Cu" { 0 } else { sys.nl - 1 };
        let zt = if li == 0 { 0.0 } else { -thickness };
        let ci = a.cells.iter().enumerate().filter(|(_, c)| c.1 == bi).min_by(|p, q| { let dd = |c: &[f64; 6]| ((c[0] + c[3]) / 2.0 - l.at[0]).hypot((c[1] + c[4]) / 2.0 - l.at[1]).hypot((c[2] + c[5]) / 2.0 - zt); dd(&p.1 .2).partial_cmp(&dd(&q.1 .2)).unwrap() }).map(|(i, _)| i).unwrap();
        let (c, r) = (((l.at[0] - x0) / d) as usize, ((l.at[1] - y0) / d) as usize);
        let bn = li * w * h + r.min(h - 1) * w + c.min(w - 1); let gc = 1.0 / l.resistance_k_per_w;
        sys.edges.push((bn, a.cells[ci].0, gc)); a.board.push((bn, ci, gc));
    }
    // each body's heat, spread over its volume
    for (bi, b) in s.bodies.iter().enumerate() {
        let watts = b.watts.unwrap_or(0.0) + b.heat_of.as_deref().map_or(0.0, hot_watts);
        if watts <= 0.0 { continue; }
        let vol: f64 = a.cells.iter().filter(|c| c.1 == bi).map(|c| (c.2[3] - c.2[0]) * (c.2[4] - c.2[1]) * (c.2[5] - c.2[2])).sum();
        for c in a.cells.iter().filter(|c| c.1 == bi) { a.sources.push((c.0, watts * (c.2[3] - c.2[0]) * (c.2[4] - c.2[1]) * (c.2[5] - c.2[2]) / vol)); }
    }
    Ok(a)
}

fn r3(x: f64) -> f64 { (x * 1000.0).round() / 1000.0 }

impl Attached {
    /// The bodies the solve left out of the board's heat without saying so (issue 28): a declared boardContactM2kW that no
    /// face met, and a body carrying nothing, with no heat of its own and no path to the board or to a body that has one.
    /// Returns (severity, body, what, fix).
    pub fn unused(&self, s: &Solids, thickness: f64) -> Vec<(&'static str, usize, String, String)> {
        let nb = s.bodies.len();
        let mut par: Vec<usize> = (0..nb).collect();
        fn find(p: &mut [usize], mut x: usize) -> usize { while p[x] != x { p[x] = p[p[x]]; x = p[x]; } x }
        for &(ca, cb, ..) in &self.faces { let (ra, rb) = (find(&mut par, self.cells[ca].1), find(&mut par, self.cells[cb].1)); par[ra] = rb; }
        let mut live = vec![false; nb]; // a group of touching bodies with heat in it or a path to the board
        for &(_, c, _) in &self.board { let r = find(&mut par, self.cells[c].1); live[r] = true; }
        let body_of: std::collections::HashMap<usize, usize> = self.cells.iter().map(|c| (c.0, c.1)).collect();
        for (nd, _) in &self.sources { if let Some(bi) = body_of.get(nd) { let r = find(&mut par, *bi); live[r] = true; } }
        let mut out = Vec::new();
        for (bi, b) in s.bodies.iter().enumerate() {
            let (lo, hi) = (b.boxes_mm.iter().map(|x| x[2]).fold(f64::MAX, f64::min), b.boxes_mm.iter().map(|x| x[5]).fold(f64::MIN, f64::max));
            if b.board_contact_m2k_w.is_some() && !self.on_copper[bi] {
                out.push(("high", bi, format!("{}: boardContactM2kW is set but none of its faces is on the board's copper, so no heat crosses that contact: the body spans z = {:.3} to {:.3} mm, and the copper faces are z = 0 and z = {:.3} mm, within the outline", b.name, lo, hi, -thickness), format!("Put a face of {} on the copper: the top face at z = {:.3} mm for a body under the board, the bottom face at z = 0 for one on top, overlapping the outline. If the gap is real, model it as a body with its own conductivity and remove boardContactM2kW.", b.name, -thickness)));
            } else if !live[find(&mut par, bi)] {
                out.push(("medium", bi, format!("{}: touches no copper and has no heat of its own, so it carries 0 W and changes nothing: it spans z = {:.3} to {:.3} mm", b.name, lo, hi), format!("To take heat from the board, set boardContactM2kW on {} with a face at z = {:.3} mm (under) or z = 0 (on top), or add a boardLink; otherwise remove it.", b.name, -thickness)));
            }
        }
        out
    }
    /// per cell: temperature and the heat flux vector (W/mm², from the face flows), and per body: the extremes and where
    /// its heat went
    pub fn result(&self, s: &Solids, t: &[f64]) -> (Value, Value) {
        let n = self.cells.len();
        let mut q = vec![[0.0f64; 3]; n];
        for &(a, b, gc, ax, area) in &self.faces { let f = gc * (t[self.cells[a].0] - t[self.cells[b].0]) / area; q[a][ax] += f / 2.0; q[b][ax] += f / 2.0; }
        let mut bodies = Vec::new();
        for (bi, b) in s.bodies.iter().enumerate() {
            let cells: Vec<usize> = (0..n).filter(|i| self.cells[*i].1 == bi).collect();
            let tt: Vec<f64> = cells.iter().map(|i| t[self.cells[*i].0]).collect();
            let air: f64 = self.air.iter().filter(|(c, _)| self.cells[*c].1 == bi).map(|(c, g)| g * t[self.cells[*c].0]).sum();
            let board: f64 = self.board.iter().filter(|(_, c, _)| self.cells[*c].1 == bi).map(|(bn, c, g)| g * (t[self.cells[*c].0] - t[*bn])).sum();
            let src: f64 = self.sources.iter().filter(|(nd, _)| cells.iter().any(|c| self.cells[*c].0 == *nd)).map(|(_, w)| w).sum();
            let qmax = cells.iter().map(|i| (q[*i][0].powi(2) + q[*i][1].powi(2) + q[*i][2].powi(2)).sqrt()).fold(0.0, f64::max);
            bodies.push(json!({"name": b.name, "cells": cells.len(), "maxRiseC": r3(tt.iter().cloned().fold(f64::MIN, f64::max)), "minRiseC": r3(tt.iter().cloned().fold(f64::MAX, f64::min)), "sourceW": r3(src), "toAirW": r3(air), "toBoardW": r3(board), "onCopper": self.on_copper[bi], "peakFluxWPerMm2": r3(qmax)}));
        }
        let per_cell = json!({"t": (0..n).map(|i| (t[self.cells[i].0] * 100.0).round() / 100.0).collect::<Vec<_>>(), "q": q.iter().map(|v| [ (v[0] * 1e5).round() / 1e5, (v[1] * 1e5).round() / 1e5, (v[2] * 1e5).round() / 1e5 ]).collect::<Vec<_>>()});
        (json!(bodies), per_cell)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    fn body(name: &str, bx: [f64; 6], k: f64, h: f64) -> Body { Body { name: name.into(), boxes_mm: vec![bx], conductivity_w_mk: k, watts: None, heat_of: None, h_still_w_m2k: h, h_moving_w_m2k: h, board_contact_m2k_w: None, provenance: String::new() } }
    #[test]
    fn slab_over_its_top_face_matches_the_series_network() {
        // a 10 x 10 x 10 mm block, 1 W in it, only its top face convecting: T_mean - T_top = q L / 2k-ish; check the energy
        let mut b = body("slab", [0., 0., 0., 10., 10., 10.], 10.0, 0.0); b.watts = Some(1.0);
        let mut top = body("lid", [0., 0., 10., 10., 10., 11.], 200.0, 100.0); top.provenance = "test".into();
        let s = Solids { pitch_mm: 1.0, bodies: vec![b, top], contacts: vec![Contact { bodies: ["slab".into(), "lid".into()], resistance_m2k_w: 0.0, provenance: String::new() }], board_links: vec![], assumptions: vec![] };
        let g = s.validate([0., 0., 1., 1.], &|_| true, &|_| false).unwrap();
        let mut sys = System::new(1, 2, 2); sys.leak[0] = vec![1e-6; 4];
        let a = attach(&mut sys, &s, &g, 0.0, 0.5, 100.0, 100.0, 0.5, 1e-5, 1.6, (1.55, 1.65), &|_| 0.0).unwrap();
        let mut rhs = vec![0.0; sys.n()]; for (nd, w) in &a.sources { rhs[*nd] += w; }
        let (t, _, r) = sys.solve(&rhs, 1e-10, 20000); assert!(r < 1e-8);
        let (bodies, _) = a.result(&s, &t);
        let air: f64 = bodies.as_array().unwrap().iter().map(|b| b["toAirW"].as_f64().unwrap()).sum();
        assert!((air - 1.0).abs() < 2e-3, "{air}");
        // the lid convects on its top and its four 1 mm sides: 1 W over 140 mm² at 100 W/m²K is 71.4 K, plus a little conduction
        let lid = bodies[1]["maxRiseC"].as_f64().unwrap(); assert!(lid > 71.0 && lid < 73.0, "{lid}");
        assert!(bodies[0]["maxRiseC"].as_f64().unwrap() > lid);
    }
    // issue 28: a top face inside the window couples to B.Cu; a cell of the same face past the board's edge stays in the air
    #[test]
    fn under_board_face_couples_inside_the_window_and_only_over_the_board() {
        let mut bd = body("carrier", [0., 0., -3.57, 4., 2., -1.57], 50.0, 10.0); bd.board_contact_m2k_w = Some(1e-4);
        let s = Solids { pitch_mm: 2.0, bodies: vec![bd], contacts: vec![], board_links: vec![], assumptions: vec![] };
        let g = s.validate([0., 0., 2., 2.], &|_| true, &|_| false).unwrap();
        let mut sys = System::new(1, 2, 2); sys.leak[0] = vec![1e-6; 4];
        let a = attach(&mut sys, &s, &g, 0.0, 0.5, 0.0, 0.0, 1.0, 1e-5, 1.57, (1.485, 1.62), &|_| 0.0).unwrap();
        assert!(a.on_copper[0]); assert_eq!(a.board.iter().map(|x| x.1).max(), Some(0), "only the cell over the board couples");
        let top = 2.0 * 2.0 / (1.0 / 1e-5 + 1.0 / 0.05);
        assert!((a.air[1].1 - a.air[0].1 - top).abs() < 1e-9, "the off-board cell keeps its top face in the air: {:?}", a.air);
        let mut miss = System::new(1, 2, 2); miss.leak[0] = vec![1e-6; 4];
        let m = attach(&mut miss, &s, &g, 0.0, 0.5, 0.0, 0.0, 1.0, 1e-5, 1.57, (1.0, 1.5), &|_| 0.0).unwrap();
        assert!(!m.on_copper[0] && m.board.is_empty()); assert_eq!(m.unused(&s, 1.57).len(), 1);
    }
    #[test]
    fn touching_bodies_need_a_contact() {
        let s = Solids { pitch_mm: 1.0, bodies: vec![body("a", [0., 0., 0., 1., 1., 1.], 1.0, 1.0), body("b", [1., 0., 0., 2., 1., 1.], 1.0, 1.0)], contacts: vec![], board_links: vec![], assumptions: vec![] };
        let g = s.validate([0., 0., 1., 1.], &|_| true, &|_| false).unwrap();
        let mut sys = System::new(1, 2, 2);
        assert!(attach(&mut sys, &s, &g, 0.0, 0.5, 0.0, 0.0, 0.5, 1e-5, 1.6, (1.55, 1.65), &|_| 0.0).is_err());
    }
    #[test]
    fn overlapping_bodies_are_refused() {
        let s = Solids { pitch_mm: 1.0, bodies: vec![body("a", [0., 0., 0., 2., 1., 1.], 1.0, 1.0), body("b", [1., 0., 0., 3., 1., 1.], 1.0, 1.0)], contacts: vec![], board_links: vec![], assumptions: vec![] };
        assert!(s.validate([0., 0., 1., 1.], &|_| true, &|_| false).is_err());
    }
}