// The assembled linear stiffness against the reference forces, on flat and
// curved plates.
import { test } from "node:test";
import assert from "node:assert/strict";
import { makeGrid, withMask, type Plate, nodeMasses, nodeX, nodeY } from "../src/params.ts";
import { uniformPlate, accelerations, fields } from "../src/model/reference.ts";
import { assemble, stiffnessMatvec, toGrid, fromGrid, isFlatPlate } from "../src/model/assemble.ts";

const bronze = { E: 110e9, nu: 0.34, rho: 8800, eta: 0, thermal: null };

function checkAgainstReference(pl: Plate, fieldsN: 1 | 3, seed = 3): number {
  const g = pl.grid;
  const sys = assemble(pl, fieldsN);
  const masses = nodeMasses(pl);
  let s = seed >>> 0;
  const rnd = () => { s = (s * 1664525 + 1013904223) >>> 0; return s / 4294967296 - 0.5; };
  const x = new Float64Array(sys.n);
  for (let i = 0; i < sys.n; i++) x[i] = rnd();
  const y = new Float64Array(sys.n);
  stiffnessMatvec(sys, x, y);
  // The reference force at a tiny amplitude, scaled back: -M a(eps x) / eps.
  const eps = 1e-9;
  const fx = toGrid(sys, x.map((v) => v * eps), g.n);
  const out = fields(g.n);
  accelerations(pl, fx, out);
  const fm = { u: out.u.map((a, k) => (-a * masses[k]) / eps), v: out.v.map((a, k) => (-a * masses[k]) / eps), w: out.w.map((a, k) => (-a * masses[k]) / eps) };
  const yref = fromGrid(sys, fm);
  let err = 0, nrm = 0;
  for (let i = 0; i < sys.n; i++) { err += (y[i] - yref[i]) ** 2; nrm += yref[i] ** 2; }
  return Math.sqrt(err / nrm);
}

test("K x equals the linearized reference force on a flat rectangle (one and three fields)", () => {
  const pl = uniformPlate(makeGrid(0.2, 0.15, 0.02), bronze, 1e-3);
  assert.ok(isFlatPlate(pl));
  assert.ok(checkAgainstReference(pl, 1) < 1e-6);
  assert.ok(checkAgainstReference(pl, 3) < 1e-6);
});

test("K x equals the linearized reference force on a bowed, tapered plate with a hole", () => {
  const g = withMask(makeGrid(0.2, 0.2, 0.0125), (x, y) => x * x + y * y <= 0.1 ** 2 && x * x + y * y >= 0.01 ** 2);
  const w0 = new Float64Array(g.n);
  for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) w0[j * g.nx + i] = 3e-3 * Math.sin(20 * nodeX(g, i)) * Math.cos(15 * nodeY(g, j));
  const pl: Plate = uniformPlate(g, bronze, 1e-3, w0);
  for (let k = 0; k < g.n; k++) if (pl.h[k] > 0) pl.h[k] = 1e-3 * (1 + 0.4 * Math.abs(nodeX(g, k % g.nx)) / 0.1);
  assert.ok(!isFlatPlate(pl));
  const sys = assemble(pl);
  assert.equal(sys.fields, 3);
  const rel = checkAgainstReference(pl, 3);
  assert.ok(rel < 1e-6, `relative error ${rel}`);
});

test("the band is narrow: the shorter grid axis is innermost", () => {
  const wide = assemble(uniformPlate(makeGrid(0.4, 0.1, 0.01), bronze, 1e-3));
  const tall = assemble(uniformPlate(makeGrid(0.1, 0.4, 0.01), bronze, 1e-3));
  assert.equal(wide.bw, tall.bw);
  assert.ok(wide.bw <= 2 * 11 + 2, `bw ${wide.bw}`);
});