Frame analysis: loads, releases, end forces and reactions

/** Local load pieces of one load on a member (global directions split onto the local axes). */
function pieces(load, ax, out) {
  const L = ax.L;
  const comps = (dir) => {
    if (dir === 'x') return [1, 0, 0];
    if (dir === 'y') return [0, 1, 0];
    if (dir === 'z') return [0, 0, 1];
    const g = dir === 'X' ? [1, 0, 0] : dir === 'Y' ? [0, 1, 0] : [0, 0, 1];
    const d = (v) => v[0] * g[0] + v[1] * g[1] + v[2] * g[2];
    return [d(ax.ex), d(ax.ey), d(ax.ez)];
  };
  const c = comps(load.dir);
  for (let k = 0; k < 3; k++) {
    if (Math.abs(c[k]) < 1e-15) continue;
    if (load.type === 'point') {
      if (load.a > L + 1e-9) continue;
      out.push({ comp: k, kind: 'point', a: Math.min(load.a, L), P: load.P * c[k] });
    } else {
      const a = Math.max(0, Math.min(L, load.a ?? 0)), b = Math.max(0, Math.min(L, load.b ?? L));
      if (b - a < 1e-9) continue;
      const w1 = load.w, w2 = load.w2 ?? load.w;
      // a load cut short by the member's end keeps its slope
      const span = (load.b ?? L) - (load.a ?? 0);
      const at = (x) => (span > 0 ? w1 + (w2 - w1) * (x - (load.a ?? 0)) / span : w1);
      out.push({ comp: k, kind: 'dist', a, b, w1: at(a) * c[k], w2: at(b) * c[k] });
    }
  }
}

function invSmall(M, n) {        // Gauss–Jordan with partial pivoting; null when singular
  const a = Array.from({ length: n }, (_, i) => [...M.slice(i * n, i * n + n), ...Array.from({ length: n }, (_, j) => (i === j ? 1 : 0))]);
  const scale = Math.max(1e-300, ...M.map(Math.abs));
  for (let c = 0; c < n; c++) {
    let p = c;
    for (let r = c + 1; r < n; r++) if (Math.abs(a[r][c]) > Math.abs(a[p][c])) p = r;
    if (Math.abs(a[p][c]) < 1e-12 * scale) return null;
    [a[c], a[p]] = [a[p], a[c]];
    const d = a[c][c];
    for (let j = 0; j < 2 * n; j++) a[c][j] /= d;
    for (let r = 0; r < n; r++) if (r !== c && a[r][c]) { const f = a[r][c]; for (let j = 0; j < 2 * n; j++) a[r][j] -= f * a[c][j]; }
  }
  const out = new Float64Array(n * n);
  for (let i = 0; i < n; i++) for (let j = 0; j < n; j++) out[i * n + j] = a[i][n + j];
  return out;
}

/** Static condensation of released dofs: { k (12×12, zero rows/cols at released), cond (for the loads), rel } or null. */
function condense(k, relIdx) {
  if (!relIdx.length) return { k, rel: relIdx, keep: null };
  const keep = [];
  for (let i = 0; i < 12; i++) if (!relIdx.includes(i)) keep.push(i);
  const r = relIdx.length;
  const krr = new Float64Array(r * r);
  for (let a = 0; a < r; a++) for (let b = 0; b < r; b++) krr[a * r + b] = k[relIdx[a] * 12 + relIdx[b]];
  const inv = invSmall(krr, r);
  if (!inv) return null;
  // G = k_ar · k_rr⁻¹ (12 × r, rows for kept dofs)
  const Gm = new Float64Array(12 * r);
  for (const i of keep) for (let b = 0; b < r; b++) { let s = 0; for (let c = 0; c < r; c++) s += k[i * 12 + relIdx[c]] * inv[c * r + b]; Gm[i * r + b] = s; }
  const kc = new Float64Array(144);
  for (const i of keep) for (const j of keep) { let s = k[i * 12 + j]; for (let b = 0; b < r; b++) s -= Gm[i * r + b] * k[relIdx[b] * 12 + j]; kc[i * 12 + j] = s; }
  return { k: kc, rel: relIdx, keep, Gm, r, inv, kFull: k };
}

function condenseLoads(c, p) {
  if (!c.rel.length) return p;
  const out = new Float64Array(12);
  for (const i of c.keep) { let s = p[i]; for (let b = 0; b < c.r; b++) s -= c.Gm[i * c.r + b] * p[c.rel[b]]; out[i] = s; }
  return out;
}

// k_global = Tᵀ k T, block by block (3×3 blocks of the 12×12)
const kg = new Float64Array(144);
for (let bi = 0; bi < 4; bi++) for (let bj = 0; bj < 4; bj++) {
  for (let r = 0; r < 3; r++) for (let q = 0; q < 3; q++) {
    let s2 = 0;
    for (let a = 0; a < 3; a++) for (let b = 0; b < 3; b++) s2 += T[a * 3 + r] * c.k[(bi * 3 + a) * 12 + bj * 3 + b] * T[b * 3 + q];
    kg[(bi * 3 + r) * 12 + bj * 3 + q] = s2;
  }
}

const idle = new Uint8Array(6 * nn);
let maxT = 0, maxR = 0;
for (let q = 0; q < 6 * nn; q++) if (q % 6 < 3) maxT = Math.max(maxT, diag[q]); else maxR = Math.max(maxR, diag[q]);
const rotSprings = [];
for (let n = 0; n < nn; n++) {
  if (!touched[n]) continue;
  for (let d = 0; d < 3; d++) {
    const q = 6 * n + d;
    if (!held[q] && diag[q] + spring[q] <= 1e-14 * maxT) { idle[q] = 1; warnings.push({ code: 'idle-dof', node: model.nodes[n].id, dof: d }); }
  }
  const blk = Float64Array.from(rotBlock[n]);
  for (let d = 0; d < 3; d++) { if (held[6 * n + 3 + d]) blk[d * 4] += maxR || 1; else blk[d * 4] += spring[6 * n + 3 + d]; }
  const { values, vectors } = eigSym3(blk);
  for (let k = 0; k < 3; k++) {
    if (values[k] > 1e-10 * (maxR || 1)) continue;
    const v = vectors[k];
    rotSprings.push({ n, v, s: maxR || 1 });
    const dof = [0, 1, 2].reduce((b, a) => (Math.abs(v[a]) > Math.abs(v[b]) ? a : b), 0);
    warnings.push({ code: 'idle-rotation', node: model.nodes[n].id, dof: 3 + dof, axis: v.map((x) => Math.round(x * 1e6) / 1e6) });
  }
}

/** Eigenvalues (ascending) and unit eigenvectors of a symmetric 3×3 (row-major), by cyclic Jacobi rotations. */
export function eigSym3(m) {
  const a = [[m[0], m[1], m[2]], [m[3], m[4], m[5]], [m[6], m[7], m[8]]];
  const V = [[1, 0, 0], [0, 1, 0], [0, 0, 1]];
  for (let sweep = 0; sweep < 50; sweep++) {
    const off = Math.abs(a[0][1]) + Math.abs(a[0][2]) + Math.abs(a[1][2]);
    if (off <= 1e-300 || off < 1e-15 * (Math.abs(a[0][0]) + Math.abs(a[1][1]) + Math.abs(a[2][2]))) break;
    for (const [p, q] of [[0, 1], [0, 2], [1, 2]]) {
      if (a[p][q] === 0) continue;
      const th = (a[q][q] - a[p][p]) / (2 * a[p][q]);
      const t = Math.sign(th || 1) / (Math.abs(th) + Math.sqrt(th * th + 1));
      const c = 1 / Math.sqrt(t * t + 1), s = t * c;
      for (let k = 0; k < 3; k++) { const x = a[k][p], y = a[k][q]; a[k][p] = c * x - s * y; a[k][q] = s * x + c * y; }
      for (let k = 0; k < 3; k++) { const x = a[p][k], y = a[q][k]; a[p][k] = c * x - s * y; a[q][k] = s * x + c * y; }
      for (let k = 0; k < 3; k++) { const x = V[k][p], y = V[k][q]; V[k][p] = c * x - s * y; V[k][q] = s * x + c * y; }
    }
  }
  const order = [0, 1, 2].sort((x, y) => a[x][x] - a[y][y]);
  return { values: order.map((k) => a[k][k]), vectors: order.map((k) => [V[0][k], V[1][k], V[2][k]]) };
}

/** The loads of one set in the solver's terms: node forces (global) + per member local pieces. */
function loadSet(prep, entries) {
  const { model, mem, idx, nn } = prep;
  const nodeF = new Float64Array(6 * nn);
  const memLoads = mem.map(() => []);
  const mIdx = new Map(model.members.map((m, k) => [m.id, k]));
  for (const { c, f } of entries) {
    if (c.selfWeight) {
      mem.forEach((e, k) => {
        const w = -c.selfWeight * f * e.rho * 1e-9 * e.A * G_ACC;     // N/mm along −Z
        if (w) pieces({ type: 'udl', dir: 'Z', w }, e.ax, memLoads[k]);
      });
    }
    for (const l of c.loads) {
      if (l.type === 'node') { const n = idx.get(l.node); for (let d = 0; d < 6; d++) nodeF[6 * n + d] += f * l.F[d]; continue; }
      const k = mIdx.get(l.member);
      const scaled = l.type === 'point' ? { ...l, P: l.P * f } : { ...l, w: l.w * f, w2: (l.w2 ?? l.w) * f };
      pieces(scaled, mem[k].ax, memLoads[k]);
    }
  }
  return { nodeF, memLoads };
}

const back = (x) => {
  const u = new Float64Array(6 * nn);
  for (let q = 0; q < 6 * nn; q++) if (eq[q] >= 0) u[q] = x[eq[q]];
  const ends = new Float64Array(12 * mem.length);
  const reactions = new Float64Array(6 * nn);
  mem.forEach((e, k) => {
    const ul = new Float64Array(12);
    for (let blk = 0; blk < 4; blk++) {
      const base = 6 * (blk < 2 ? e.i : e.j) + (blk % 2) * 3;
      for (let a = 0; a < 3; a++) ul[blk * 3 + a] = e.T[a * 3] * u[base] + e.T[a * 3 + 1] * u[base + 1] + e.T[a * 3 + 2] * u[base + 2];
    }
    const kc = e.c.k, cp = pc[k].cp;
    for (let r = 0; r < 12; r++) { let s = -cp[r]; for (let q = 0; q < 12; q++) s += kc[r * 12 + q] * ul[q]; ends[12 * k + r] = s; }
    for (const r of e.c.rel) ends[12 * k + r] = 0;
    for (let blk = 0; blk < 4; blk++) for (let r = 0; r < 3; r++) {
      let s = 0;
      for (let a = 0; a < 3; a++) s += e.T[a * 3 + r] * ends[12 * k + blk * 3 + a];
      reactions[6 * (blk < 2 ? e.i : e.j) + (blk % 2) * 3 + r] += s;
    }
  });
  for (let q = 0; q < 6 * nn; q++) reactions[q] -= nodeF[q];
  return { u, ends, reactions };
};

let x = solveBand(prep, F);
let { u, ends, reactions } = back(x);
if (!prep.rotSprings.length) {
  const r = new Float64Array(n); let big = 0, ref = 0;
  for (let q = 0; q < 6 * nn; q++) if (eq[q] >= 0 && !prep.idle[q] && !spr[q]) { r[eq[q]] = -reactions[q]; big = Math.max(big, Math.abs(r[eq[q]])); }
  for (let i = 0; i < n; i++) ref = Math.max(ref, Math.abs(F[i]));
  if (big > 1e-13 * (ref || 1)) {
    const dx = solveBand(prep, r), x2 = new Float64Array(n);
    for (let i = 0; i < n; i++) x2[i] = x[i] + dx[i];
    const t = back(x2); let big2 = 0;
    for (let q = 0; q < 6 * nn; q++) if (eq[q] >= 0 && !prep.idle[q] && !spr[q]) big2 = Math.max(big2, Math.abs(t.reactions[q]));
    if (big2 < big) { x = x2; ({ u, ends, reactions } = t); }   // kept only when it reduced the out-of-balance
  }
}

const applied = totals(prep, nodeF, memLoads), reacted = new Float64Array(6);
for (let nd = 0; nd < nn; nd++) {
  const [x0, y0, z0] = prep.P[nd];
  const R = reactions.subarray(6 * nd, 6 * nd + 6);
  reacted[0] += R[0]; reacted[1] += R[1]; reacted[2] += R[2];
  reacted[3] += R[3] + y0 * R[2] - z0 * R[1]; reacted[4] += R[4] + z0 * R[0] - x0 * R[2]; reacted[5] += R[5] + x0 * R[1] - y0 * R[0];
}
return { u, ends, reactions, loads: memLoads, applied: [...applied], reacted: [...reacted] };

/** Resultant of every load of a set about the global origin: [Fx Fy Fz Mx My Mz]. */
function totals(prep, nodeF, memLoads) {
  const t = new Float64Array(6);
  const add = (p, F) => { t[0] += F[0]; t[1] += F[1]; t[2] += F[2]; t[3] += p[1] * F[2] - p[2] * F[1]; t[4] += p[2] * F[0] - p[0] * F[2]; t[5] += p[0] * F[1] - p[1] * F[0]; };
  for (let nd = 0; nd < prep.nn; nd++) {
    const F = nodeF.subarray(6 * nd, 6 * nd + 6);
    add(prep.P[nd], F); t[3] += F[3]; t[4] += F[4]; t[5] += F[5];
  }
  prep.mem.forEach((e, k) => {
    const o = prep.P[e.i], ax = [e.ax.ex, e.ax.ey, e.ax.ez];
    for (const pz of memLoads[k]) {
      const dir = ax[pz.comp];
      const at = (x, q) => add([o[0] + e.ax.ex[0] * x, o[1] + e.ax.ex[1] * x, o[2] + e.ax.ex[2] * x], [dir[0] * q, dir[1] * q, dir[2] * q]);
      if (pz.kind === 'point') at(pz.a, pz.P);
      else { const h = (pz.b - pz.a) / 2; for (const [g, wt] of GAUSS4) { const x = pz.a + h * (g + 1); at(x, (pz.w1 + (pz.w2 - pz.w1) * (x - pz.a) / (pz.b - pz.a)) * wt * h); } }
    }
  });
  return t;
}

/** A result as a linear combination of computed sets (results must come from the same prep). */
export function combine(prep, list, meta = {}) {
  const nn = prep.nn, M = prep.mem.length;
  const out = { ...meta, u: new Float64Array(6 * nn), reactions: new Float64Array(6 * nn), ends: new Float64Array(12 * M), loads: prep.mem.map(() => []), applied: [0, 0, 0, 0, 0, 0], reacted: [0, 0, 0, 0, 0, 0] };
  for (const { R, f } of list) {
    for (let q = 0; q < out.u.length; q++) { out.u[q] += f * R.u[q]; out.reactions[q] += f * R.reactions[q]; }
    for (let q = 0; q < out.ends.length; q++) out.ends[q] += f * R.ends[q];
    for (let k = 0; k < 6; k++) { out.applied[k] += f * R.applied[k]; out.reacted[k] += f * R.reacted[k]; }
    R.loads.forEach((ls, m) => { for (const pz of ls) out.loads[m].push(pz.kind === 'point' ? { ...pz, P: pz.P * f } : { ...pz, w1: pz.w1 * f, w2: pz.w2 * f }); });
  }
  if (prep.shell) combineShell(prep, list, out);
  return out;
}