3D frame analysis

/** 12×12 local stiffness (row-major Float64Array). */
export function localK(E, G, A, Iy, Iz, J, L) {
  const k = new Float64Array(144);
  const s = (i, j, v) => { k[i * 12 + j] += v; if (i !== j) k[j * 12 + i] += v; };
  const ea = E * A / L, gj = G * J / L;
  s(0, 0, ea); s(0, 6, -ea); s(6, 6, ea);
  s(3, 3, gj); s(3, 9, -gj); s(9, 9, gj);
  const z1 = 12 * E * Iz / L ** 3, z2 = 6 * E * Iz / L ** 2, z3 = 4 * E * Iz / L, z4 = 2 * E * Iz / L;
  s(1, 1, z1); s(1, 5, z2); s(1, 7, -z1); s(1, 11, z2);
  s(5, 5, z3); s(5, 7, -z2); s(5, 11, z4);
  s(7, 7, z1); s(7, 11, -z2); s(11, 11, z3);
  const y1 = 12 * E * Iy / L ** 3, y2 = 6 * E * Iy / L ** 2, y3 = 4 * E * Iy / L, y4 = 2 * E * Iy / L;
  s(2, 2, y1); s(2, 4, -y2); s(2, 8, -y1); s(2, 10, -y2);
  s(4, 4, y3); s(4, 8, y2); s(4, 10, y4);
  s(8, 8, y1); s(8, 10, y2); s(10, 10, y3);
  return k;
}

/** Equivalent nodal loads (local, 12) of one load piece on a member of length L. */
export function equivalentLoads(piece, L, p = new Float64Array(12)) {
  const add = (xi, q) => {      // q = force (per length × weight, or a point force) at ξ along the member
    const N1 = 1 - 3 * xi * xi + 2 * xi ** 3, N2 = L * (xi - 2 * xi * xi + xi ** 3), N3 = 3 * xi * xi - 2 * xi ** 3, N4 = L * (-xi * xi + xi ** 3);
    if (piece.comp === 0) { p[0] += q * (1 - xi); p[6] += q * xi; }
    else if (piece.comp === 1) { p[1] += q * N1; p[5] += q * N2; p[7] += q * N3; p[11] += q * N4; }
    else { p[2] += q * N1; p[4] -= q * N2; p[8] += q * N3; p[10] -= q * N4; }
  };
  if (piece.kind === 'point') add(piece.a / L, piece.P);
  else {
    const { a, b, w1, w2 } = piece, h = (b - a) / 2;
    for (const [g, wt] of GAUSS4) { const x = a + h * (g + 1); add(x / L, (w1 + (w2 - w1) * (x - a) / (b - a)) * wt * h); }
  }
  return p;
}

// banded Cholesky, in place: L(i, j) at i·B + (i − j)
for (let j = 0; j < n; j++) {
  const kjj = K[j * B];
  let s = kjj;
  for (let k = Math.max(0, j - band); k < j; k++) { const l = K[j * B + (j - k)]; s -= l * l; }
  if (!(s > PIVOT_TOL * Math.abs(kjj)) || !(kjj > 0)) {
    const q = dofOf[j];
    return { error: { code: 'mechanism', node: model.nodes[Math.floor(q / 6)].id, dof: q % 6 }, warnings };
  }
  const d = Math.sqrt(s);
  K[j * B] = d;
  for (let i = j + 1; i <= Math.min(n - 1, j + band); i++) {
    let v = K[i * B + (i - j)];
    for (let k = Math.max(0, i - band); k < j; k++) v -= K[i * B + (i - k)] * K[j * B + (j - k)];
    K[i * B + (i - j)] = v / d;
  }
}

function solveBand(prep, f) {
  const { n, band, B, L } = prep;
  const y = new Float64Array(n);
  for (let i = 0; i < n; i++) { let s = f[i]; for (let k = Math.max(0, i - band); k < i; k++) s -= L[i * B + (i - k)] * y[k]; y[i] = s / L[i * B]; }
  for (let i = n - 1; i >= 0; i--) { let s = y[i]; for (let k = i + 1; k <= Math.min(n - 1, i + band); k++) s -= L[k * B + (k - i)] * y[k]; y[i] = s / L[i * B]; }
  return y;
}