// The discrete shell, as in cymsim/shell.py: thickness and height sampled on the mesh, the CST
// slope of the nodal heights, element stiffnesses, lumped mass with the inflated rotational
// inertia, the six M-orthonormal rigid-body modes, the node-to-element adjacency used by the
// gather of element forces, and the Irons element bound on the largest eigenfrequency that sets
// the stable time step. Everything here is float64 and computed once per configuration.
import type { Material, Numerics } from "../config/schema.ts";
import { dktStiffness } from "./dkt.ts";
import type { CymbalProfile } from "./geometry.ts";
import { symmetricEigenvalues } from "./linalg.ts";
import type { Mesh } from "./mesh.ts";

export const N_DOF = 5;   // (u, v, w, theta_x, theta_y) per node

export interface Shell {
  mesh: Mesh; material: Material; numerics: Numerics; nl: number;
  hE: Float64Array; hNode: Float64Array; z0Node: Float64Array; z0x: Float64Array; z0y: Float64Array;
  De: Float64Array; Cm: Float64Array;      // (E,) and (E, 3) = (C11, C12, C33)
  Ke: Float64Array;                        // (E, 81)
  rhoHA: Float64Array;                     // (E,) rho h A
  M: Float64Array; Minv: Float64Array;     // (N, 5)
  Rm: Float64Array;                        // (6, N, 5), M-orthonormal
  adjOffsets: Uint32Array; adjEntries: Uint32Array;   // CSR: node -> list of (3 e + local index)
  mass: number; maxSlope: number; omegaMax: number; dtStable: number; limitingElement: number;
  zeroModes: [number, number];             // min and max count of zero-energy element modes (expect 6)
}

export function buildShell(mesh: Mesh, profile: CymbalProfile, material: Material, numerics: Numerics): Shell {
  const N = mesh.n, E = mesh.e, { xy, tri, area, Gx, Gy } = mesh;
  const nl = numerics.nonlinear ? 1 : 0;
  const hE = new Float64Array(E), hNode = new Float64Array(N), z0Node = new Float64Array(N);
  for (let e = 0; e < E; e++) hE[e] = profile.thickness(mesh.centroid[2 * e], mesh.centroid[2 * e + 1]);
  for (let n = 0; n < N; n++) { hNode[n] = profile.thickness(xy[2 * n], xy[2 * n + 1]); z0Node[n] = profile.z0(xy[2 * n], xy[2 * n + 1]); }
  const z0x = new Float64Array(E), z0y = new Float64Array(E);
  let maxSlope = 0;
  for (let e = 0; e < E; e++) {
    let sx = 0, sy = 0;
    for (let i = 0; i < 3; i++) { const z = z0Node[tri[3 * e + i]]; sx += Gx[3 * e + i] * z; sy += Gy[3 * e + i] * z; }
    z0x[e] = sx; z0y[e] = sy; maxSlope = Math.max(maxSlope, Math.hypot(sx, sy));
  }
  const { E: Ey, nu, rho } = material;
  const De = new Float64Array(E), Cm = new Float64Array(3 * E), Ke = new Float64Array(81 * E), rhoHA = new Float64Array(E);
  for (let e = 0; e < E; e++) {
    De[e] = (Ey * hE[e] ** 3) / (12 * (1 - nu * nu));
    const c = (Ey * hE[e]) / (1 - nu * nu);
    Cm[3 * e] = c; Cm[3 * e + 1] = c * nu; Cm[3 * e + 2] = c * 0.5 * (1 - nu);
    rhoHA[e] = rho * hE[e] * area[e];
    const xe = [xy[2 * tri[3 * e]], xy[2 * tri[3 * e + 1]], xy[2 * tri[3 * e + 2]]];
    const ye = [xy[2 * tri[3 * e] + 1], xy[2 * tri[3 * e + 1] + 1], xy[2 * tri[3 * e + 2] + 1]];
    Ke.set(dktStiffness(xe, ye, De[e], nu), 81 * e);
  }
  // lumped mass
  const M = new Float64Array(N_DOF * N);
  const mE = new Float64Array(E), IE = new Float64Array(E);
  for (let e = 0; e < E; e++) {
    mE[e] = (rho * hE[e] * area[e]) / 3;
    IE[e] = mE[e] * ((hE[e] * hE[e]) / 12 + numerics.rot_inertia_factor * area[e]);
    for (let k = 0; k < 3; k++) {
      const n = tri[3 * e + k];
      M[5 * n] += mE[e] * numerics.inplane_mass_scale; M[5 * n + 1] += mE[e] * numerics.inplane_mass_scale;
      M[5 * n + 2] += mE[e]; M[5 * n + 3] += IE[e]; M[5 * n + 4] += IE[e];
    }
  }
  const Minv = M.map((m) => 1 / m);
  let mass = 0;
  for (let n = 0; n < N; n++) mass += M[5 * n + 2];
  // rigid-body modes of the curved shell, M-orthonormalised by Gram-Schmidt
  let cx = 0, cy = 0;
  for (let n = 0; n < N; n++) { cx += xy[2 * n]; cy += xy[2 * n + 1]; }
  cx /= N; cy /= N;
  const Rm = new Float64Array(6 * N_DOF * N);
  const set = (k: number, n: number, d: number, val: number) => { Rm[(k * N + n) * N_DOF + d] = val; };
  for (let n = 0; n < N; n++) {
    const X = xy[2 * n] - cx, Y = xy[2 * n + 1] - cy, z = z0Node[n];
    set(0, n, 0, 1); set(1, n, 1, 1);
    set(2, n, 0, -Y); set(2, n, 1, X);
    set(3, n, 2, 1);
    set(4, n, 2, Y); set(4, n, 3, 1); set(4, n, 1, -z);
    set(5, n, 2, -X); set(5, n, 4, 1); set(5, n, 0, z);
  }
  const L = N_DOF * N;
  const mdot = (a: number, b: number) => { let s = 0; for (let i = 0; i < L; i++) s += Rm[a * L + i] * M[i] * Rm[b * L + i]; return s; };
  for (let k = 0; k < 6; k++) {
    for (let l = 0; l < k; l++) { const c = mdot(l, k); for (let i = 0; i < L; i++) Rm[k * L + i] -= c * Rm[l * L + i]; }
    const nrm = Math.sqrt(mdot(k, k));
    for (let i = 0; i < L; i++) Rm[k * L + i] /= nrm;
  }
  // adjacency (CSR)
  const deg = new Uint32Array(N);
  for (let e = 0; e < E; e++) for (let k = 0; k < 3; k++) deg[tri[3 * e + k]]++;
  const adjOffsets = new Uint32Array(N + 1);
  for (let n = 0; n < N; n++) adjOffsets[n + 1] = adjOffsets[n] + deg[n];
  const adjEntries = new Uint32Array(adjOffsets[N]), fill = new Uint32Array(N);
  for (let e = 0; e < E; e++) for (let k = 0; k < 3; k++) { const n = tri[3 * e + k]; adjEntries[adjOffsets[n] + fill[n]++] = 3 * e + k; }
  // Irons bound: lambda_max(M^-1 K) <= max_e lambda_max(M_e^-1 K_e), K_e the element tangent at rest
  let lamMax = 0, limitingElement = 0, zmin = 15, zmax = 0;
  const S = new Float64Array(225), me = new Float64Array(15);
  for (let e = 0; e < E; e++) {
    elementTangent(e, Gx, Gy, area[e], z0x[e], z0y[e], Cm, Ke, S);
    for (let i = 0; i < 3; i++) {
      me[5 * i] = mE[e] * numerics.inplane_mass_scale; me[5 * i + 1] = mE[e] * numerics.inplane_mass_scale;
      me[5 * i + 2] = mE[e]; me[5 * i + 3] = IE[e]; me[5 * i + 4] = IE[e];
    }
    for (let i = 0; i < 15; i++) for (let j = 0; j < 15; j++) S[15 * i + j] /= Math.sqrt(me[i] * me[j]);
    const ev = symmetricEigenvalues(S, 15);
    const top = ev[14];
    let nz = 0; for (const v of ev) if (v < 1e-9 * top) nz++;
    zmin = Math.min(zmin, nz); zmax = Math.max(zmax, nz);
    if (top > lamMax) { lamMax = top; limitingElement = e; }
  }
  const omegaMax = Math.sqrt(lamMax);
  return { mesh, material, numerics, nl, hE, hNode, z0Node, z0x, z0y, De, Cm, Ke, rhoHA, M, Minv, Rm, adjOffsets, adjEntries,
           mass, maxSlope, omegaMax, dtStable: 2 / omegaMax, limitingElement, zeroModes: [zmin, zmax] };
}

/** 15 x 15 tangent stiffness of one element at rest: A B^T C B (membrane, linear part) plus DKT on the bending slots. */
export function elementTangent(e: number, Gx: Float64Array, Gy: Float64Array, A: number, z0x: number, z0y: number,
                               Cm: Float64Array, Ke: Float64Array, out: Float64Array) {
  out.fill(0);
  // B rows: exx, eyy, gxy over the 15 DOFs (node-major u, v, w, thx, thy)
  const B = [new Float64Array(15), new Float64Array(15), new Float64Array(15)];
  for (let i = 0; i < 3; i++) {
    const gx = Gx[3 * e + i], gy = Gy[3 * e + i];
    B[0][5 * i] = gx; B[0][5 * i + 2] = z0x * gx;
    B[1][5 * i + 1] = gy; B[1][5 * i + 2] = z0y * gy;
    B[2][5 * i] = gy; B[2][5 * i + 1] = gx; B[2][5 * i + 2] = z0x * gy + z0y * gx;
  }
  const C = [[Cm[3 * e], Cm[3 * e + 1], 0], [Cm[3 * e + 1], Cm[3 * e], 0], [0, 0, Cm[3 * e + 2]]];
  for (let i = 0; i < 15; i++) for (let j = 0; j < 15; j++) {
    let s = 0;
    for (let k = 0; k < 3; k++) for (let l = 0; l < 3; l++) s += B[k][i] * C[k][l] * B[l][j];
    out[15 * i + j] = A * s;
  }
  const slot = [2, 3, 4, 7, 8, 9, 12, 13, 14];   // (w, thx, thy) x 3 within the 15
  for (let i = 0; i < 9; i++) for (let j = 0; j < 9; j++) out[15 * slot[i] + slot[j]] += Ke[81 * e + 9 * i + j];
}