// The float64 reference discretization: symmetry and rigid modes of the
// linear stiffness, forces as the gradient of the energy, dissipativity,
// thickness scaling, the stand, and the free-plate frequencies.
import { test } from "node:test";
import assert from "node:assert/strict";
import { makeGrid, withMask, type Grid, type Plate, nodeMasses, standRates, insideNode, nodeX, nodeY, bendingStiffness } from "../src/params.ts";
import { uniformPlate, accelerations, addDamping, potentialEnergy, fields, maxAngularFrequency, type Fields } from "../src/model/reference.ts";
import { assemble, stiffnessEntry } from "../src/model/assemble.ts";
import { fullLoss } from "../src/model/reference.ts";

const steel = { E: 200e9, nu: 0.3, rho: 7850 };

function rnd(seed: number): () => number {
  let s = seed >>> 0;
  return () => { s = (s * 1664525 + 1013904223) >>> 0; return s / 4294967296 - 0.5; };
}
function randomFields(pl: Plate, amp: number, seed: number): Fields {
  const g = pl.grid, r = rnd(seed), f = fields(g.n);
  for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) {
    if (!insideNode(g, i, j)) continue;
    const k = j * g.nx + i;
    f.u[k] = 0.02 * amp * r(); f.v[k] = 0.02 * amp * r(); f.w[k] = amp * r();
  }
  return f;
}
function bowed(g: Grid, sag: number): Float64Array {
  const w0 = new Float64Array(g.n);
  const Lx = (g.nx - 1) * g.hx;
  for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) w0[j * g.nx + i] = sag * (1 - (2 * nodeX(g, i) / Lx) ** 2) + 0.3 * sag * Math.sin(3 * nodeY(g, j) / Lx);
  return w0;
}
function withStand(pl: Plate): Plate {
  return { ...pl, stand: { x: 0, y: 0, pad: { innerRadius: 0.005, outerRadius: 0.03, thickness: 0.01, E: 1e6, lossFactor: 0.3 }, topPad: true, preload: 1e-3 } };
}

test("the linear stiffness is symmetric with exactly three rigid modes per block (uniform, tapered, disc with hole)", () => {
  const cases: { name: string; pl: Plate; fields: 1 | 3 }[] = [];
  const g1 = makeGrid(0.12, 0.09, 0.015);
  cases.push({ name: "uniform rectangle", pl: uniformPlate(g1, steel, 1e-3), fields: 3 });
  const tapered = uniformPlate(g1, steel, 1e-3);
  for (let k = 0; k < g1.n; k++) if (tapered.h[k] > 0) tapered.h[k] = 1e-3 * (1 + 0.5 * (k % g1.nx) / g1.nx);
  cases.push({ name: "tapered rectangle", pl: tapered, fields: 3 });
  const gd = withMask(makeGrid(0.16, 0.16, 0.01), (x, y) => x * x + y * y <= 0.08 ** 2 && x * x + y * y >= 0.015 ** 2);
  cases.push({ name: "disc with hole", pl: uniformPlate(gd, steel, 1e-3), fields: 1 });
  for (const c of cases) {
    const sys = assemble(c.pl, c.fields);
    // Symmetry is built into the banded storage; check the rigid modes: K x = 0 for translations and rotations.
    const g = c.pl.grid;
    const nullVectors: Float64Array[] = [];
    const mk = (fn: (i: number, j: number, comp: number) => number) => {
      const x = new Float64Array(sys.n);
      for (let d = 0; d < sys.n; d++) {
        const k = sys.dofNode[d];
        const comp = sys.fields === 1 ? 2 : d % 3;
        x[d] = fn(k % g.nx, Math.floor(k / g.nx), comp);
      }
      return x;
    };
    nullVectors.push(mk((_i, _j, c) => (c === 2 ? 1 : 0)));
    nullVectors.push(mk((i, _j, c) => (c === 2 ? nodeX(g, i) : 0)));
    nullVectors.push(mk((_i, j, c) => (c === 2 ? nodeY(g, j) : 0)));
    if (c.fields === 3) {
      nullVectors.push(mk((_i, _j, c) => (c === 0 ? 1 : 0)));
      nullVectors.push(mk((_i, _j, c) => (c === 1 ? 1 : 0)));
      nullVectors.push(mk((i, j, c) => (c === 0 ? -nodeY(g, j) : c === 1 ? nodeX(g, i) : 0)));
    }
    const y = new Float64Array(sys.n);
    for (const x of nullVectors) {
      const { stiffnessMatvec } = require_assemble();
      stiffnessMatvec(sys, x, y);
      let ny = 0, nx = 0;
      for (let d = 0; d < sys.n; d++) { ny += y[d] * y[d]; nx += x[d] * x[d]; }
      const scale = Math.abs(stiffnessEntry(sys, 0, 0)) * Math.sqrt(nx);
      assert.ok(Math.sqrt(ny) < 1e-9 * scale, `${c.name}: rigid mode residual ${Math.sqrt(ny) / scale}`);
    }
    // A non-rigid vector must not be in the null space.
    const z = mk((i, j, c) => (c === 2 ? Math.sin(i * 0.9) * Math.cos(j * 1.1) : 0));
    const { stiffnessMatvec } = require_assemble();
    stiffnessMatvec(sys, z, y);
    assert.ok(y.some((v) => Math.abs(v) > 1e-3), `${c.name}: bending is not trivial`);
  }
});

import * as assembleModule from "../src/model/assemble.ts";
function require_assemble(): typeof assembleModule { return assembleModule; }

test("forces are the exact gradient of the discrete energy (flat and bowed, tapered, free and on a stand)", () => {
  const g = makeGrid(0.1, 0.08, 0.01);
  const variants: { name: string; pl: Plate }[] = [];
  variants.push({ name: "flat free", pl: uniformPlate(g, steel, 1e-3) });
  variants.push({ name: "bowed free", pl: uniformPlate(g, steel, 1e-3, bowed(g, 3e-3)) });
  const tapered = uniformPlate(g, steel, 1e-3, bowed(g, 2e-3));
  for (let k = 0; k < g.n; k++) if (tapered.h[k] > 0) tapered.h[k] = 1e-3 * (1 + 0.4 * (k % g.nx) / g.nx);
  variants.push({ name: "bowed tapered", pl: tapered });
  variants.push({ name: "bowed on a stand", pl: withStand(uniformPlate(g, steel, 1e-3, bowed(g, 3e-3))) });
  for (const v of variants) {
    const pl = v.pl;
    const sr = standRates(pl);
    const f = randomFields(pl, 2.5e-3, 7);
    const masses = nodeMasses(pl);
    const acc = fields(g.n);
    accelerations(pl, f, acc, false, sr);
    const dir = randomFields(pl, 1, 11);
    // Directional derivative of U along dir versus -F . dir.
    const eps = 1e-7;
    const shift = (s: number): Fields => ({ u: f.u.map((x, k) => x + s * dir.u[k]), v: f.v.map((x, k) => x + s * dir.v[k]), w: f.w.map((x, k) => x + s * dir.w[k]) });
    const dU = (potentialEnergy(pl, shift(eps), sr) - potentialEnergy(pl, shift(-eps), sr)) / (2 * eps);
    let fd = 0;
    for (let k = 0; k < g.n; k++) fd -= masses[k] * (acc.u[k] * dir.u[k] + acc.v[k] * dir.v[k] + acc.w[k] * dir.w[k]);
    assert.ok(Math.abs(dU - fd) < 1e-6 * Math.abs(dU), `${v.name}: dU ${dU} vs -F.dir ${fd}`);
  }
});

test("the loss terms and the pad dashpots only remove energy", () => {
  const g = makeGrid(0.1, 0.08, 0.01);
  const pl = withStand(uniformPlate(g, steel, 1e-3));
  const sr = standRates(pl);
  const loss = fullLoss({ sigma0: 1, sigma1: 5e-4, kappa: 1 }, sr, 20);
  const masses = nodeMasses(pl);
  for (const seed of [1, 2, 3]) {
    const vel = randomFields(pl, 1, seed);
    const acc = fields(g.n);
    addDamping(pl, loss, vel, acc, sr);
    let power = 0;
    for (let k = 0; k < g.n; k++) power += masses[k] * (acc.u[k] * vel.u[k] + acc.v[k] * vel.v[k] + acc.w[k] * vel.w[k]);
    assert.ok(power < 0, `power ${power}`);
  }
});

test("doubling the thickness doubles the bending frequencies and leaves the membrane frequencies alone", () => {
  const g = makeGrid(0.1, 0.1, 0.01);
  const w1 = maxAngularFrequency(g, steel, 1e-3), w2 = maxAngularFrequency(g, steel, 2e-3);
  // The largest frequency is set by the in-plane waves, independent of h.
  assert.ok(Math.abs(w1 / w2 - 1) < 1e-6, `in-plane omega_max ${w1} vs ${w2}`);
  const { stiffnessEntry: ke } = assembleModule;
  const s1 = assemble(uniformPlate(g, steel, 1e-3)), s2 = assemble(uniformPlate(g, steel, 2e-3));
  const c = Math.floor(s1.n / 2);
  // Bending stiffness scales as h^3 and the mass as h: omega^2 scales as h^2.
  assert.ok(Math.abs(ke(s2, c, c) / ke(s1, c, c) - 8) < 1e-9);
  assert.ok(Math.abs(s2.mass[c] / s1.mass[c] - 2) < 1e-12);
  assert.ok(Math.abs(bendingStiffness(steel, 2e-3) / bendingStiffness(steel, 1e-3) - 8) < 1e-12);
});

test("the stand resists a rigid translation with its total stiffness and a tilt with the pads' second moment", () => {
  const g = withMask(makeGrid(0.2, 0.2, 0.01), (x, y) => x * x + y * y <= 0.1 ** 2 && x * x + y * y >= 0.01 ** 2);
  const pl = withStand(uniformPlate(g, steel, 1e-3));
  const sr = standRates(pl)!;
  const masses = nodeMasses(pl);
  const f = fields(g.n), acc = fields(g.n);
  const dz = 1e-4;
  for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? dz : 0;
  accelerations(pl, f, acc, false, sr);
  let F = 0;
  for (let k = 0; k < g.n; k++) F += masses[k] * acc.w[k];
  assert.ok(Math.abs(-F / dz / sr.totalKw - 1) < 1e-9, `translation stiffness ${-F / dz} vs ${sr.totalKw}`);
  // Tilt about y: w = theta x; restoring moment = theta * sum_k kw_k x_k^2 (bending of a plane is
  // zero; the von Karman membrane strain of a rotation is theta^2 / 2, so a tiny tilt keeps it negligible).
  const th = 1e-7;
  for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? th * nodeX(g, k % g.nx) : 0;
  accelerations(pl, f, acc, false, sr);
  let M = 0, I2 = 0;
  for (let k = 0; k < g.n; k++) { M += masses[k] * acc.w[k] * nodeX(g, k % g.nx); I2 += sr.kw[k] * nodeX(g, k % g.nx) ** 2; }
  assert.ok(Math.abs(-M / th / I2 - 1) < 1e-6, `tilt stiffness ${-M / th} vs ${I2}`);
  // With lift-off, a small motion within the preload is the same linear spring; a large upward motion unloads the bottom pad.
  for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? 1e-5 : 0;
  const lin = fields(g.n), lo = fields(g.n);
  accelerations(pl, f, lin, false, sr);
  accelerations(pl, f, lo, true, sr);
  for (let k = 0; k < g.n; k++) assert.ok(Math.abs(lin.w[k] - lo.w[k]) < 1e-9 * Math.max(1, Math.abs(lin.w[k])));
  for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? 0.01 : 0;
  accelerations(pl, f, lo, true, sr);
  let Flo = 0;
  for (let k = 0; k < g.n; k++) Flo += masses[k] * lo.w[k];
  // Bottom pad unloaded: only the top pad (compressed by preload + 10 mm) pushes down, minus the rest balance.
  const kb = sr.totalKw / 2;
  const expected = -(kb * (sr.preloadTop + 0.01) - kb * sr.preloadTop) - kb * sr.preloadBottom;
  assert.ok(Math.abs(Flo / expected - 1) < 1e-6, `lift-off force ${Flo} vs ${expected}`);
});