1ee087eCymbal designer: scripted design, thickness field, mount, single strikeJeremy Magland 1// Rigid-body projection: (1) mass-weighted sums of the six rigid modes for
2// displacement and velocity, kinetic energy and max |w| in one workgroup;
3// (2) subtraction of the rigid components from every node. The in-plane
4// inertia may be scaled by P.inPlaneMass; the rigid coefficients are ratios
5// and do not see the factor, the kinetic energy does. The node mass carries
6// the local thickness.
8@group(0) @binding(0) var<uniform> P: Params;
9@group(0) @binding(1) var<storage, read_write> pos: array<f32>; // [u | v | w | wLo]
10@group(0) @binding(2) var<storage, read_write> vel: array<f32>; // [du | dv | dw]
11@group(0) @binding(3) var<storage, read_write> stats: array<f32>; // 0..11 coefficients, 12 KE, 13 max|w|
13const NS = 14u;
14var<workgroup> sh: array<array<f32, NS>, 256>;
16@compute @workgroup_size(256)
17fn sums(@builtin(local_invocation_id) lid: vec3u) {
18 let t = lid.x;
19 var s: array<f32, NS>;
20 for (var q = 0u; q < NS; q++) { s[q] = 0.0; }
21 for (var k = t; k < P.n; k += 256u) {
22 let i = i32(k % P.nx); let j = i32(k / P.nx);
23 let x = nodeX(i); let y = nodeY(j);
24 let m = nodeMass(i, j);
25 let n = P.n;
26 let p = vec3f(pos[k], pos[n + k], pos[2u * n + k]);
27 let v = vec3f(vel[k], vel[n + k], vel[2u * n + k]);
28 s[0] += m * p.x; s[1] += m * p.y; s[2] += m * p.z;
29 s[3] += m * (x * p.y - y * p.x); s[4] += m * y * p.z; s[5] += m * (-x * p.z);
30 s[6] += m * v.x; s[7] += m * v.y; s[8] += m * v.z;
31 s[9] += m * (x * v.y - y * v.x); s[10] += m * y * v.z; s[11] += m * (-x * v.z);
32 s[12] += 0.5 * m * (P.inPlaneMass * (v.x * v.x + v.y * v.y) + v.z * v.z);
33 s[13] = max(s[13], abs(p.z));
34 }
35 for (var q = 0u; q < NS; q++) { sh[t][q] = s[q]; }
36 workgroupBarrier();
37 for (var stride = 128u; stride > 0u; stride >>= 1u) {
38 if (t < stride) {
39 for (var q = 0u; q < 13u; q++) { sh[t][q] += sh[t + stride][q]; }
40 sh[t][13] = max(sh[t][13], sh[t + stride][13]);
41 }
42 workgroupBarrier();
43 }
44 if (t == 0u) {
45 let inv = array<f32, 6>(1.0 / P.mass, 1.0 / P.mass, 1.0 / P.mass, 1.0 / P.Iz, 1.0 / P.Ix, 1.0 / P.Iy);
46 for (var q = 0u; q < 6u; q++) {
47 stats[q] = sh[0][q] * inv[q];
48 stats[q + 6u] = sh[0][q + 6u] * inv[q];
49 }
50 // Kinetic energy of the non-rigid motion (the six rigid modes are mass-orthogonal).
51 let ip = array<f32, 6>(P.inPlaneMass, P.inPlaneMass, 1.0, P.inPlaneMass, 1.0, 1.0);
52 var keRigid = 0.0;
53 for (var q = 6u; q < 12u; q++) { keRigid += 0.5 * ip[q - 6u] * sh[0][q] * sh[0][q] * inv[q - 6u]; }
54 stats[12] = max(sh[0][12] - keRigid, 0.0);
55 stats[13] = sh[0][13];
56 }
57}
59@compute @workgroup_size(256)
60fn apply(@builtin(global_invocation_id) gid: vec3u) {
61 let k = gid.x;
62 if (k >= P.n) { return; }
63 let i = i32(k % P.nx); let j = i32(k / P.nx);
64 if (nodeWeight(i, j) == 0.0) { return; } // outside the plate
65 let x = nodeX(i); let y = nodeY(j);
66 let n = P.n;
67 pos[k] -= stats[0] - stats[3] * y;
68 pos[n + k] -= stats[1] + stats[3] * x;
69 pos[2u * n + k] -= stats[2] + stats[4] * y - stats[5] * x;
70 vel[k] -= stats[6] - stats[9] * y;
71 vel[n + k] -= stats[7] + stats[9] * x;
72 vel[2u * n + k] -= stats[8] + stats[10] * y - stats[11] * x;
73}