1// The Discrete Kirchhoff Triangle (Batoz, Bathe, Ho 1980): 9 x 9 bending stiffness for nodal
2// DOFs (w, theta_x, theta_y) x 3 with theta_x = dw/dy, theta_y = -dw/dx. This is the same
3// element as cymsim/elements.py, with the derivatives of the rotation shape functions written
4// out instead of obtained by automatic differentiation.
6const GAUSS = [[1 / 6, 1 / 6], [2 / 3, 1 / 6], [1 / 6, 2 / 3]];
8/** Quadratic shape functions N1..N6 and their (xi, eta) derivatives at (xi, eta). */
9function shapes(xi: number, eta: number) {
10 const lam = 1 - xi - eta;
11 const N = [lam * (2 * lam - 1), xi * (2 * xi - 1), eta * (2 * eta - 1), 4 * xi * eta, 4 * eta * lam, 4 * xi * lam];
12 const Nxi = [-(4 * lam - 1), 4 * xi - 1, 0, 4 * eta, -4 * eta, 4 * (lam - xi)];
13 const Neta = [-(4 * lam - 1), 0, 4 * eta - 1, 4 * xi, 4 * (lam - eta), -4 * xi];
14 return { N, Nxi, Neta };
15}
17/** Row vectors Hx, Hy (length 9) as linear combinations of the six shape values N. */
18function H(N: number[], a: number[], b: number[], c: number[], d: number[], e: number[]): [number[], number[]] {
19 const [N1, N2, N3, N4, N5, N6] = N;
20 const [a4, a5, a6] = a, [b4, b5, b6] = b, [c4, c5, c6] = c, [d4, d5, d6] = d, [e4, e5, e6] = e;
21 const Hx = [1.5 * (a6 * N6 - a5 * N5), b5 * N5 + b6 * N6, N1 - c5 * N5 - c6 * N6,
22 1.5 * (a4 * N4 - a6 * N6), b6 * N6 + b4 * N4, N2 - c4 * N4 - c6 * N6,
23 1.5 * (a5 * N5 - a4 * N4), b4 * N4 + b5 * N5, N3 - c4 * N4 - c5 * N5];
24 const Hy = [1.5 * (d6 * N6 - d5 * N5), -N1 + e5 * N5 + e6 * N6, -b5 * N5 - b6 * N6,
25 1.5 * (d4 * N4 - d6 * N6), -N2 + e4 * N4 + e6 * N6, -b6 * N6 - b4 * N4,
26 1.5 * (d5 * N5 - d4 * N4), -N3 + e4 * N4 + e5 * N5, -b4 * N4 - b5 * N5];
27 return [Hx, Hy];
28}
30/** DKT stiffness (row-major 9 x 9) of the triangle with nodes (xe[i], ye[i]), rigidity D, Poisson ratio nu. */
31export function dktStiffness(xe: ArrayLike<number>, ye: ArrayLike<number>, D: number, nu: number): Float64Array {
32 const I = [1, 2, 0], J = [2, 0, 1]; // sides k = 4, 5, 6 <-> ij = 23, 31, 12
33 const a: number[] = [], b: number[] = [], c: number[] = [], d: number[] = [], e: number[] = [];
34 for (let k = 0; k < 3; k++) {
35 const xij = xe[I[k]] - xe[J[k]], yij = ye[I[k]] - ye[J[k]], l2 = xij * xij + yij * yij;
36 a.push(-xij / l2); b.push((0.75 * xij * yij) / l2); c.push((0.25 * xij * xij - 0.5 * yij * yij) / l2);
37 d.push(-yij / l2); e.push((0.25 * yij * yij - 0.5 * xij * xij) / l2);
38 }
39 const x21 = xe[1] - xe[0], x31 = xe[2] - xe[0], y21 = ye[1] - ye[0], y31 = ye[2] - ye[0];
40 const det = x21 * y31 - x31 * y21;
41 const xi_x = y31 / det, xi_y = -x31 / det, eta_x = -y21 / det, eta_y = x21 / det;
42 const Db = [[D, D * nu, 0], [D * nu, D, 0], [0, 0, D * 0.5 * (1 - nu)]];
43 const Ke = new Float64Array(81);
44 for (const [xi, eta] of GAUSS) {
45 const { Nxi, Neta } = shapes(xi, eta);
46 const [Hx_xi, Hy_xi] = H(Nxi, a, b, c, d, e);
47 const [Hx_eta, Hy_eta] = H(Neta, a, b, c, d, e);
48 const B: number[][] = [[], [], []];
49 for (let i = 0; i < 9; i++) {
50 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;
51 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;
52 B[0][i] = Hx_x; B[1][i] = Hy_y; B[2][i] = Hx_y + Hy_x;
53 }
54 const w = (0.5 * det) / 3;
55 for (let i = 0; i < 9; i++) {
56 for (let j = 0; j < 9; j++) {
57 let s = 0;
58 for (let k = 0; k < 3; k++) for (let l = 0; l < 3; l++) s += B[k][i] * Db[k][l] * B[l][j];
59 Ke[9 * i + j] += w * s;
60 }
61 }
62 }
63 return Ke;
64}