c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 1// Dihedral angle across an edge shared by two triangles, and its gradient
2// (Bridson, Marino & Fedkiw 2003). Edge (e0, e1); wings w0 and w1 on either
3// side. theta = 0 when the two triangles are coplanar; the sign tells which way
4// the hinge folds. Positions are read from X (3 per node).
6export const hingeOut = new Float64Array(13); // [theta, g_w0(3), g_w1(3), g_e0(3), g_e1(3)]
8export function hingeAngleGrad(X: Float64Array, w0: number, w1: number, e0: number, e1: number, out = hingeOut): number {
9 const a = 3 * w0;
10 const b = 3 * w1;
11 const c = 3 * e0;
12 const d = 3 * e1;
13 // E = x_e1 - x_e0
14 const Ex = X[d] - X[c];
15 const Ey = X[d + 1] - X[c + 1];
16 const Ez = X[d + 2] - X[c + 2];
17 // N1 = (x_w0 - x_e0) x (x_w0 - x_e1)
18 const p1x = X[a] - X[c], p1y = X[a + 1] - X[c + 1], p1z = X[a + 2] - X[c + 2];
19 const q1x = X[a] - X[d], q1y = X[a + 1] - X[d + 1], q1z = X[a + 2] - X[d + 2];
20 const N1x = p1y * q1z - p1z * q1y;
21 const N1y = p1z * q1x - p1x * q1z;
22 const N1z = p1x * q1y - p1y * q1x;
23 // N2 = (x_w1 - x_e1) x (x_w1 - x_e0)
24 const p2x = X[b] - X[d], p2y = X[b + 1] - X[d + 1], p2z = X[b + 2] - X[d + 2];
25 const q2x = X[b] - X[c], q2y = X[b + 1] - X[c + 1], q2z = X[b + 2] - X[c + 2];
26 const N2x = p2y * q2z - p2z * q2y;
27 const N2y = p2z * q2x - p2x * q2z;
28 const N2z = p2x * q2y - p2y * q2x;
29 const n1sq = N1x * N1x + N1y * N1y + N1z * N1z;
30 const n2sq = N2x * N2x + N2y * N2y + N2z * N2z;
31 const eLen = Math.sqrt(Ex * Ex + Ey * Ey + Ez * Ez);
32 const n1 = Math.sqrt(n1sq);
33 const n2 = Math.sqrt(n2sq);
34 // angle between unit normals, signed about the edge direction
35 const cos = (N1x * N2x + N1y * N2y + N1z * N2z) / (n1 * n2);
36 const cx = N1y * N2z - N1z * N2y;
37 const cy = N1z * N2x - N1x * N2z;
38 const cz = N1x * N2y - N1y * N2x;
39 const sin = (cx * Ex + cy * Ey + cz * Ez) / (n1 * n2 * eLen);
40 const theta = Math.atan2(sin, cos);
41 // gradients
42 const s1 = eLen / n1sq;
43 const s2 = eLen / n2sq;
44 out[0] = theta;
45 out[1] = -s1 * N1x;
46 out[2] = -s1 * N1y;
47 out[3] = -s1 * N1z;
48 out[4] = -s2 * N2x;
49 out[5] = -s2 * N2y;
50 out[6] = -s2 * N2z;
51 const t1 = (q1x * Ex + q1y * Ey + q1z * Ez) / eLen / n1sq; // (x_w0 - x_e1).E / |E| / |N1|^2
52 const t2 = (p2x * Ex + p2y * Ey + p2z * Ez) / eLen / n2sq; // (x_w1 - x_e1).E / |E| / |N2|^2
53 const r1 = (p1x * Ex + p1y * Ey + p1z * Ez) / eLen / n1sq; // (x_w0 - x_e0).E
54 const r2 = (q2x * Ex + q2y * Ey + q2z * Ez) / eLen / n2sq; // (x_w1 - x_e0).E
55 out[7] = -(t1 * N1x + t2 * N2x);
56 out[8] = -(t1 * N1y + t2 * N2y);
57 out[9] = -(t1 * N1z + t2 * N2z);
58 out[10] = r1 * N1x + r2 * N2x;
59 out[11] = r1 * N1y + r2 * N2y;
60 out[12] = r1 * N1z + r2 * N2z;
61 return theta;
62}