import { describe, expect, it } from "vitest";
import {
  LIBRARY_FILES,
  buildCompareBundle,
  isPlainStruct,
  methodBlockers,
  getMethod,
  parseCompareResult,
  validateSpec,
} from "../src/index";
import { maxAbs, runner, solve, testSpec } from "./helpers";

/** Run a MATLAB snippet against the library and return the struct it leaves in `out`. */
async function script(source: string) {
  const run = await runner.run({
    bundle: {
      files: [{ name: "c.m", source }, ...LIBRARY_FILES],
      entry: "c.m",
      resultVariable: "out",
    },
  });
  if (!run.ok) throw new Error(`${run.error}\n${run.output.join("")}`);
  if (!isPlainStruct(run.value)) throw new Error("unexpected shape");
  return run.value.fields as Record<string, number>;
}

function ellipseSpec() {
  const s = testSpec();
  s.plate.shape = "ellipse";
  return s;
}

describe("the plate's footprint on the grid", () => {
  it("is the whole grid for a rectangle and a staircase for an ellipse, with every node carrying area", async () => {
    const out = await script(`
      r = plate_domain(struct('shape', 'rectangle', 'Lx', 0.3, 'Ly', 0.2), 19, 13);
      out.r_inside = nnz(r.inside); out.r_interior = nnz(r.interior); out.r_boundary = nnz(r.boundary);
      out.r_m_corner = r.m(1, 1); out.r_m_edge = r.m(1, 7); out.r_m_in = r.m(5, 7); out.r_full3 = nnz(r.full3);
      e = plate_domain(struct('shape', 'ellipse', 'Lx', 0.3, 'Ly', 0.2), 31, 21);
      out.e_inside = nnz(e.inside); out.e_area = nnz(e.inside) / (31 * 21);
      out.e_min_m = min(e.m(e.inside)); out.e_max_m = max(e.m(:));
      out.e_corner = e.inside(1, 1); out.e_centre = e.inside(11, 16);
      out.e_near_ok = all(e.full3(e.near(e.inside)));
      out.e_kappa_min = min(e.kappa(e.boundary)); out.e_kappa_max = max(e.kappa(e.boundary));
      nrm = sqrt(e.normal_x(e.boundary).^2 + e.normal_y(e.boundary).^2);
      out.e_normal_unit = max(abs(nrm - 1));
      % a sliver: the disk on a coarse grid still has full stencils
      d = plate_domain(struct('shape', 'ellipse', 'Lx', 1, 'Ly', 1), 9, 9);
      out.d_full3 = nnz(d.full3);
    `);
    expect(out.r_inside).toBe(19 * 13);
    expect(out.r_interior).toBe(17 * 11);
    expect(out.r_boundary).toBe(19 * 13 - 17 * 11);
    expect(out.r_m_corner).toBe(0.25);
    expect(out.r_m_edge).toBe(0.5);
    expect(out.r_m_in).toBe(1);
    expect(out.r_full3).toBe(17 * 11);
    // pi/4 of the box, to within the staircase
    expect(Math.abs(out.e_area - Math.PI / 4)).toBeLessThan(0.05);
    expect(out.e_min_m).toBeGreaterThanOrEqual(0.25);
    expect(out.e_max_m).toBe(1);
    expect(out.e_corner).toBe(false);
    expect(out.e_centre).toBe(true);
    expect(out.e_near_ok).toBe(true);
    // curvature of a 300 x 200 ellipse lies between b/a^2 and a/b^2
    expect(out.e_kappa_min).toBeGreaterThan(0.1 / 0.15 ** 2 - 1e-9);
    expect(out.e_kappa_max).toBeLessThan(0.15 / 0.1 ** 2 + 1e-9);
    expect(out.e_normal_unit).toBeLessThan(1e-12);
    expect(out.d_full3).toBeGreaterThan(0);
  });

  it("keeps the strike and the probe on the plate", () => {
    const spec = ellipseSpec();
    expect(validateSpec(spec).some((i) => i.severity === "error")).toBe(false);
    // inside the box, outside the ellipse
    spec.excitation.x = 0.01;
    spec.excitation.y = 0.01;
    const issues = validateSpec(spec);
    expect(issues.some((i) => i.path === "excitation" && /ellipse/.test(i.message))).toBe(true);
  });
});

describe("circular plate frequencies", () => {
  // omega a^2 sqrt(rho h / D) for the disk of radius a, nu = 0.3, from the Bessel-function
  // frequency equations (clamped: J_n I_{n+1} + I_n J_{n+1} = 0; simply supported:
  // J_{n+1}/J_n + I_{n+1}/I_n = 2 lambda / (1 - nu)) and, for free edges, the tabulated values
  // (Leissa, Vibration of Plates, 1969). The staircase edge is first order, so a 25-node grid
  // is expected within a few percent.
  it("clamped: 10.216", async () => {
    const out = await script(`
      dom = plate_domain(struct('shape', 'ellipse', 'Lx', 2, 'Ly', 2), 25, 25);
      [~, K, M] = plate_matrices(dom, 'clamped', 0.3);
      ev = sort(eig(full(K)));
      out.f1 = sqrt(ev(1)); out.f2 = sqrt(ev(2)); out.f3 = sqrt(ev(3));
      out.ratio = ev(end) / (4/dom.dx^2 + 4/dom.dy^2)^2;
      out.n = size(K, 1);
    `);
    expect(Math.abs(out.f1 - 10.216) / 10.216).toBeLessThan(0.02);
    // the (1, 0) pair, degenerate on the true disk
    expect(Math.abs(out.f2 - 21.26) / 21.26).toBeLessThan(0.03);
    expect(Math.abs(out.f3 - 21.26) / 21.26).toBeLessThan(0.03);
    // the explicit step bound derived for the interior stencil still holds
    expect(out.ratio).toBeLessThan(1);
    expect(out.ratio).toBeGreaterThan(0.9);
  });

  it("simply supported: 4.935 with the curvature term, the paradox value 5.783 without it", async () => {
    const out = await script(`
      dom = plate_domain(struct('shape', 'ellipse', 'Lx', 2, 'Ly', 2), 25, 25);
      [L, K] = plate_matrices(dom, 'simply-supported', 0.3);
      ev = sort(eig(full(K)));
      out.f1 = sqrt(ev(1)); out.f2 = sqrt(ev(2));
      out.minev = ev(1) / ev(end);
      ev0 = sort(eig(full(L * L)));
      out.g1 = sqrt(ev0(1));
      out.symm = norm(full(K - K'), 1) / norm(full(K), 1);
    `);
    expect(out.symm).toBeLessThan(1e-12);
    expect(out.minev).toBeGreaterThan(0);
    expect(Math.abs(out.f1 - 4.935) / 4.935).toBeLessThan(0.03);
    expect(Math.abs(out.f2 - 13.898) / 13.898).toBeLessThan(0.03);
    // L*L alone converges to the plate with zero Laplacian on the edge, J_0's first zero squared
    expect(Math.abs(out.g1 - 5.783) / 5.783).toBeLessThan(0.03);
    expect(out.g1 / out.f1).toBeGreaterThan(1.12);
  });

  it("free: 5.358 and 9.003, with exactly three rigid-body modes", async () => {
    const out = await script(`
      dom = plate_domain(struct('shape', 'ellipse', 'Lx', 2, 'Ly', 2), 25, 25);
      [L, K, M] = plate_matrices(dom, 'free', 0.3);
      Dh = diag(1 ./ sqrt(full(diag(M)))); S = Dh * full(K) * Dh; S = (S + S') / 2;
      ev = sort(eig(S));
      out.zero3 = max(abs(ev(1:3))) / ev(end);
      out.fourth = ev(4) / ev(end);
      out.f1 = sqrt(ev(4)); out.f2 = sqrt(ev(5)); out.f3 = sqrt(ev(6));
      out.ratio = ev(end) / (4/dom.dx^2 + 4/dom.dy^2)^2;
      X = dom.X(dom.inside); Y = dom.Y(dom.inside);
      out.nullx = max(abs(K * (X - 1))) * dom.dx^4; out.nully = max(abs(K * (Y - 1))) * dom.dx^4;
      out.lnull = max(abs(L * ones(size(X))));
    `);
    expect(out.zero3).toBeLessThan(1e-9);
    expect(out.fourth).toBeGreaterThan(1e-6);
    expect(out.nullx).toBeLessThan(1e-9);
    expect(out.nully).toBeLessThan(1e-9);
    expect(out.lnull).toBeLessThan(1e-9);
    // the (2, 0) pair is split by the staircase's anisotropy; its mean is the disk's value
    expect(Math.abs((out.f1 + out.f2) / 2 - 5.358) / 5.358).toBeLessThan(0.03);
    expect(Math.abs(out.f3 - 9.003) / 9.003).toBeLessThan(0.03);
    expect(out.ratio).toBeLessThan(1);
  });
});

describe("solving on an ellipse", () => {
  it.each(["simply-supported", "clamped", "free", "pinned-center"] as const)(
    "runs both finite-difference schemes with %s edges and they agree",
    async (bc) => {
      const spec = ellipseSpec();
      spec.bc = bc;
      spec.implicit = {
        nx: spec.fd.nx,
        ny: spec.fd.ny,
        duration: spec.fd.duration,
        steps_per_sample: 2,
        theta: 0.25,
      };
      const exp = await solve(spec, "fd");
      const imp = await solve(spec, "implicit");
      for (const { result } of [exp, imp]) {
        expect(result.info.stable).toBe(true);
        expect(result.info.extra.shape).toBe("ellipse");
        expect(result.probe.w.every(Number.isFinite)).toBe(true);
        expect(maxAbs(result.probe.w)).toBeGreaterThan(0);
        expect(result.info.notes.some((n) => /staircase/.test(n))).toBe(true);
        // outside the ellipse the fields are NaN, and the same nodes in every frame and field
        const { nx, ny, nf, data } = result.frames.w;
        const nanPer = Array.from(
          { length: nf },
          (_, k) =>
            data.subarray(k * nx * ny, (k + 1) * nx * ny).filter((v) => Number.isNaN(v)).length,
        );
        expect(nanPer[0]).toBeGreaterThan(0.1 * nx * ny);
        expect(nanPer[0]).toBeLessThan(0.3 * nx * ny);
        expect(new Set(nanPer).size).toBe(1);
        expect(result.frames.vm.data.filter((v) => Number.isNaN(v)).length).toBe(nanPer[0] * nf);
        expect(result.envelope.vm.data.filter((v) => Number.isNaN(v)).length).toBe(nanPer[0]);
        // the corners of the box are off the plate, the centre is on it
        expect(data[0]).toBeNaN();
        expect(data[ny - 1 + (nx - 1) * ny]).toBeNaN();
        expect(Number.isFinite(data[Math.floor(ny / 2) + Math.floor(nx / 2) * ny])).toBe(true);
        // the probe was kept where it was asked for (0.2, 0.12), on the plate
        expect(Math.abs(result.probe.x - spec.probe.x)).toBeLessThan(
          0.5 * (spec.plate.Lx / (nx - 1)),
        );
        expect(Math.abs(result.probe.y - spec.probe.y)).toBeLessThan(
          0.5 * (spec.plate.Ly / (ny - 1)),
        );
      }
      expect(Number(imp.result.info.extra.nunknowns)).toBeLessThan(spec.fd.nx * spec.fd.ny);
      expect(Number(exp.result.info.extra.nunknowns)).toBe(Number(imp.result.info.extra.nunknowns));
      const cmp = await runner.run({
        bundle: buildCompareBundle(),
        inputs: { A: exp.run.ok ? exp.run.value : 0, B: imp.run.ok ? imp.run.value : 0 },
      });
      expect(cmp.ok).toBe(true);
      if (!cmp.ok) return;
      const c = parseCompareResult(cmp.value);
      expect(c.w_rel_l2).toBeLessThan(0.03);
      expect(c.vm_rel_l2).toBeLessThan(0.1);
      expect(c.frames_rel_l2.length).toBeGreaterThan(0);
      expect(c.frames_rel_l2.every(Number.isFinite)).toBe(true);
    },
  );

  it("is stiffer than the rectangle of the same box, and diverges above the CFL limit like it", async () => {
    // A clamped plate on a smaller domain has a higher fundamental (domain monotonicity of the
    // Dirichlet problem), so the ellipse inscribed in the box must beat the box itself.
    const out = await script(`
      r = plate_domain(struct('shape', 'rectangle', 'Lx', 0.3, 'Ly', 0.2), 25, 17);
      e = plate_domain(struct('shape', 'ellipse', 'Lx', 0.3, 'Ly', 0.2), 25, 17);
      [~, Kr] = plate_matrices(r, 'clamped', 0.3); [~, Ke] = plate_matrices(e, 'clamped', 0.3);
      out.rect = min(eig(full(Kr))); out.ell = min(eig(full(Ke)));
    `);
    expect(out.ell).toBeGreaterThan(out.rect);
    // The Rayleigh quotient of the trial function (1 - x^2/a^2 - y^2/b^2)^2 is an upper bound,
    // omega^2 rho h / D = (40/3) (3/a^4 + 2/(a^2 b^2) + 3/b^4): 10.33 against 10.216 for the disk,
    // 17.39 for a/b = 1.5. The true value lies about a percent below it, and the staircase on a
    // 25-node grid a further percent or two below that (the disk gives 10.21 at this resolution).
    const a = 0.15;
    const b = 0.1;
    const bound = Math.sqrt((40 / 3) * (3 / a ** 4 + 2 / (a * b) ** 2 + 3 / b ** 4)) * a ** 2;
    expect(bound).toBeCloseTo(17.39, 1);
    const lam = Math.sqrt(out.ell) * a ** 2;
    expect(lam).toBeLessThan(bound);
    expect(lam).toBeGreaterThan(0.94 * bound);
    const ell = ellipseSpec();
    ell.fd.cfl = 1.04;
    const blown = await solve(ell, "fd");
    expect(blown.result.info.stable).toBe(false);
    expect(blown.result.probe.w.every(Number.isFinite)).toBe(true);
  });

  it("is refused by the modal method, in the registry and in the solver", async () => {
    const spec = ellipseSpec();
    expect(methodBlockers(getMethod("modal"), spec).some((b) => /ellipse/.test(b))).toBe(true);
    expect(methodBlockers(getMethod("fd"), spec)).toEqual([]);
    expect(methodBlockers(getMethod("implicit"), spec)).toEqual([]);
    await expect(solve(spec, "modal")).rejects.toThrow(/rectangle/);
  });

  it("puts the load and the probe on the plate even when the strike is near the edge", async () => {
    const spec = ellipseSpec();
    // just inside the ellipse near its right end, footprint spilling over the edge
    spec.excitation.x = 0.28;
    spec.excitation.y = 0.1;
    spec.probe.x = 0.27;
    spec.probe.y = 0.1;
    const { result } = await solve(spec, "fd");
    expect(result.info.stable).toBe(true);
    expect(maxAbs(result.probe.w)).toBeGreaterThan(0);
    expect(Number.isFinite(result.probe.x)).toBe(true);
  });
});