// The Discrete Kirchhoff Triangle (Batoz, Bathe, Ho 1980): 9 x 9 bending stiffness for nodal
// DOFs (w, theta_x, theta_y) x 3 with theta_x = dw/dy, theta_y = -dw/dx. This is the same
// element as cymsim/elements.py, with the derivatives of the rotation shape functions written
// out instead of obtained by automatic differentiation.

const GAUSS = [[1 / 6, 1 / 6], [2 / 3, 1 / 6], [1 / 6, 2 / 3]];

/** Quadratic shape functions N1..N6 and their (xi, eta) derivatives at (xi, eta). */
function shapes(xi: number, eta: number) {
  const lam = 1 - xi - eta;
  const N = [lam * (2 * lam - 1), xi * (2 * xi - 1), eta * (2 * eta - 1), 4 * xi * eta, 4 * eta * lam, 4 * xi * lam];
  const Nxi = [-(4 * lam - 1), 4 * xi - 1, 0, 4 * eta, -4 * eta, 4 * (lam - xi)];
  const Neta = [-(4 * lam - 1), 0, 4 * eta - 1, 4 * xi, 4 * (lam - eta), -4 * xi];
  return { N, Nxi, Neta };
}

/** Row vectors Hx, Hy (length 9) as linear combinations of the six shape values N. */
function H(N: number[], a: number[], b: number[], c: number[], d: number[], e: number[]): [number[], number[]] {
  const [N1, N2, N3, N4, N5, N6] = N;
  const [a4, a5, a6] = a, [b4, b5, b6] = b, [c4, c5, c6] = c, [d4, d5, d6] = d, [e4, e5, e6] = e;
  const Hx = [1.5 * (a6 * N6 - a5 * N5), b5 * N5 + b6 * N6, N1 - c5 * N5 - c6 * N6,
              1.5 * (a4 * N4 - a6 * N6), b6 * N6 + b4 * N4, N2 - c4 * N4 - c6 * N6,
              1.5 * (a5 * N5 - a4 * N4), b4 * N4 + b5 * N5, N3 - c4 * N4 - c5 * N5];
  const Hy = [1.5 * (d6 * N6 - d5 * N5), -N1 + e5 * N5 + e6 * N6, -b5 * N5 - b6 * N6,
              1.5 * (d4 * N4 - d6 * N6), -N2 + e4 * N4 + e6 * N6, -b6 * N6 - b4 * N4,
              1.5 * (d5 * N5 - d4 * N4), -N3 + e4 * N4 + e5 * N5, -b4 * N4 - b5 * N5];
  return [Hx, Hy];
}

/** DKT stiffness (row-major 9 x 9) of the triangle with nodes (xe[i], ye[i]), rigidity D, Poisson ratio nu. */
export function dktStiffness(xe: ArrayLike<number>, ye: ArrayLike<number>, D: number, nu: number): Float64Array {
  const I = [1, 2, 0], J = [2, 0, 1];   // sides k = 4, 5, 6 <-> ij = 23, 31, 12
  const a: number[] = [], b: number[] = [], c: number[] = [], d: number[] = [], e: number[] = [];
  for (let k = 0; k < 3; k++) {
    const xij = xe[I[k]] - xe[J[k]], yij = ye[I[k]] - ye[J[k]], l2 = xij * xij + yij * yij;
    a.push(-xij / l2); b.push((0.75 * xij * yij) / l2); c.push((0.25 * xij * xij - 0.5 * yij * yij) / l2);
    d.push(-yij / l2); e.push((0.25 * yij * yij - 0.5 * xij * xij) / l2);
  }
  const x21 = xe[1] - xe[0], x31 = xe[2] - xe[0], y21 = ye[1] - ye[0], y31 = ye[2] - ye[0];
  const det = x21 * y31 - x31 * y21;
  const xi_x = y31 / det, xi_y = -x31 / det, eta_x = -y21 / det, eta_y = x21 / det;
  const Db = [[D, D * nu, 0], [D * nu, D, 0], [0, 0, D * 0.5 * (1 - nu)]];
  const Ke = new Float64Array(81);
  for (const [xi, eta] of GAUSS) {
    const { Nxi, Neta } = shapes(xi, eta);
    const [Hx_xi, Hy_xi] = H(Nxi, a, b, c, d, e);
    const [Hx_eta, Hy_eta] = H(Neta, a, b, c, d, e);
    const B: number[][] = [[], [], []];
    for (let i = 0; i < 9; i++) {
      const Hx_x = Hx_xi[i] * xi_x + Hx_eta[i] * eta_x, Hx_y = Hx_xi[i] * xi_y + Hx_eta[i] * eta_y;
      const Hy_x = Hy_xi[i] * xi_x + Hy_eta[i] * eta_x, Hy_y = Hy_xi[i] * xi_y + Hy_eta[i] * eta_y;
      B[0][i] = Hx_x; B[1][i] = Hy_y; B[2][i] = Hx_y + Hy_x;
    }
    const w = (0.5 * det) / 3;
    for (let i = 0; i < 9; i++) {
      for (let j = 0; j < 9; j++) {
        let s = 0;
        for (let k = 0; k < 3; k++) for (let l = 0; l < 3; l++) s += B[k][i] * Db[k][l] * B[l][j];
        Ke[9 * i + j] += w * s;
      }
    }
  }
  return Ke;
}