import test from "node:test";
import assert from "node:assert/strict";
import { dktStiffness } from "../src/model/dkt.ts";
import { symmetricEigenvalues, planeResiduals } from "../src/model/linalg.ts";
import { fixture, maxAbs, maxAbsDiff } from "./helpers.ts";

test("DKT stiffness matches the Python (autodiff) element", () => {
  const f = fixture.dkt;
  const Ke = dktStiffness(f.x, f.y, f.D, f.nu);
  assert.ok(maxAbsDiff(Ke, f.Ke) < 1e-9 * maxAbs(f.Ke));
});

test("DKT rigid modes, constant curvature and twist", () => {
  const x = [0, 0.03, 0.01], y = [0, 0.005, 0.025], D = 10, nu = 0.3;
  const Ke = dktStiffness(x, y, D, nu);
  const apply = (U: number[]) => Array.from({ length: 9 }, (_, i) => U.reduce((s, u, j) => s + Ke[9 * i + j] * u, 0));
  const energy = (U: number[]) => 0.5 * U.reduce((s, u, i) => s + u * apply(U)[i], 0);
  const A = 0.5 * ((x[1] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (y[1] - y[0]));
  const scale = maxAbs(Ke);
  for (const U of [[1, 0, 0, 1, 0, 0, 1, 0, 0], [y[0], 1, 0, y[1], 1, 0, y[2], 1, 0], [-x[0], 0, 1, -x[1], 0, 1, -x[2], 0, 1]])
    assert.ok(maxAbs(apply(U)) < 1e-12 * scale);
  const Ucur = x.flatMap((xi) => [0.5 * xi * xi, 0, -xi]);
  assert.ok(Math.abs(energy(Ucur) / (0.5 * D * A) - 1) < 1e-9);
  const Utw = x.map((xi, i) => [xi * y[i], xi, -y[i]]).flat();
  assert.ok(Math.abs(energy(Utw) / (D * (1 - nu) * A) - 1) < 1e-9);
  const ev = symmetricEigenvalues(Ke, 9);
  assert.ok(ev[2] < 1e-10 * ev[8] && ev[3] > 1e-6 * ev[8]);
});

test("Jacobi eigenvalues and plane fit", () => {
  const ev = symmetricEigenvalues([2, 1, 0, 1, 2, 0, 0, 0, 5], 3);
  assert.ok(maxAbsDiff(ev, [1, 3, 5]) < 1e-12);
  const res = planeResiduals([0, 1, 0, 1], [0, 0, 1, 1], [1, 3, 4, 6]);
  assert.ok(maxAbs(res) < 1e-12);
  // collinear points (a stick's line): the fit degenerates to a line, as numpy's lstsq gives
  const line = planeResiduals([-1, 0, 1], [-2, 0, 2], [1, 2, 3.5]);
  assert.ok(maxAbsDiff(line, [1 / 12, -1 / 6, 1 / 12]) < 1e-12);
  // a single point: zero residual
  assert.ok(maxAbs(planeResiduals([0.3], [0.1], [0.02])) < 1e-15);
});