8354482cymmodes: linear modes of a cymsim cymbal driven in the full nonlinear model, with the energy of every mode over timeJeremy Magland 1// Triangular mesh of the projected disk, as in cymsim/mesh.py: nodes on concentric rings,
2// annuli stitched by an angular merge, optional Laplacian smoothing, and the constant
3// gradient operators of the linear (CST) interpolation on every triangle.
5export interface Mesh {
6 n: number; e: number;
7 xy: Float64Array; // (N, 2)
8 tri: Uint32Array; // (E, 3) counter-clockwise
9 boundary: Uint32Array; // rim node indices in angular order
10 area: Float64Array; // (E,)
11 Gx: Float64Array; // (E, 3): d phi/dx = sum_i Gx[e, i] phi_i
12 Gy: Float64Array;
13 centroid: Float64Array; // (E, 2)
14}
16function stitch(aIdx: Uint32Array, aAng: Float64Array, bIdx: Uint32Array, bAng: Float64Array, out: number[]) {
17 const na = aIdx.length, nb = bIdx.length, twoPi = 2 * Math.PI;
18 const bAngU = Array.from(bAng, (t) => (((t - aAng[0]) % twoPi) + twoPi) % twoPi + aAng[0]);
19 const order = Array.from(bAngU.keys()).sort((i, j) => bAngU[i] - bAngU[j]);
20 const bi = order.map((i) => bIdx[i]), ba = order.map((i) => bAngU[i]);
21 const ta = Array.from(aAng); ta.push(aAng[0] + twoPi);
22 const tb = ba.slice(); tb.push(ba[0] + twoPi);
23 let i = 0, j = 0;
24 while (i < na || j < nb) {
25 const nextA = i < na ? ta[i + 1] : Infinity, nextB = j < nb ? tb[j + 1] : Infinity;
26 if (nextA <= nextB) { out.push(aIdx[i % na], bi[j % nb], aIdx[(i + 1) % na]); i++; }
27 else { out.push(aIdx[i % na], bi[j % nb], bi[(j + 1) % nb]); j++; }
28 }
29}
31function smooth(xy: Float64Array, tri: Uint32Array, fixed: Uint32Array, nIter: number): Float64Array {
32 const N = xy.length / 2, E = tri.length / 3;
33 const src: number[] = [], dst: number[] = [];
34 for (let e = 0; e < E; e++) {
35 for (let k = 0; k < 3; k++) {
36 const a = tri[3 * e + k], b = tri[3 * e + ((k + 1) % 3)];
37 src.push(a, b); dst.push(b, a);
38 }
39 }
40 const deg = new Float64Array(N);
41 for (const s of src) deg[s]++;
42 const free = new Uint8Array(N).fill(1);
43 for (const f of fixed) free[f] = 0;
44 let cur = Float64Array.from(xy);
45 for (let it = 0; it < nIter; it++) {
46 const acc = new Float64Array(2 * N);
47 for (let k = 0; k < src.length; k++) { acc[2 * src[k]] += cur[2 * dst[k]]; acc[2 * src[k] + 1] += cur[2 * dst[k] + 1]; }
48 const next = Float64Array.from(cur);
49 for (let n = 0; n < N; n++) if (free[n]) { next[2 * n] = acc[2 * n] / deg[n]; next[2 * n + 1] = acc[2 * n + 1] / deg[n]; }
50 cur = next;
51 }
52 return cur;
53}
55export function diskMesh(radius: number, elementSize: number, smoothing = 1): Mesh {
56 const nRings = Math.max(2, Math.round(radius / elementSize));
57 const dr = radius / nRings;
58 const nodes: number[] = [0, 0];
59 let prevIdx = new Uint32Array([0]), prevAng = new Float64Array([0]);
60 const tris: number[] = [];
61 for (let k = 1; k <= nRings; k++) {
62 const r = k * dr;
63 const nk = Math.max(6, Math.round((2 * Math.PI * r) / dr));
64 const ang = new Float64Array(nk), idx = new Uint32Array(nk);
65 const base = nodes.length / 2;
66 for (let j = 0; j < nk; j++) {
67 ang[j] = (2 * Math.PI * (j + 0.5 * (k % 2))) / nk;
68 idx[j] = base + j;
69 nodes.push(r * Math.cos(ang[j]), r * Math.sin(ang[j]));
70 }
71 if (k === 1) for (let j = 0; j < nk; j++) tris.push(0, idx[j], idx[(j + 1) % nk]);
72 else stitch(prevIdx, prevAng, idx, ang, tris);
73 prevIdx = idx; prevAng = ang;
74 }
75 let xy: Float64Array = Float64Array.from(nodes);
76 const tri = Uint32Array.from(tris);
77 if (smoothing > 0) xy = smooth(xy, tri, prevIdx, smoothing);
78 const E = tri.length / 3;
79 for (let e = 0; e < E; e++) {
80 const a = tri[3 * e], b = tri[3 * e + 1], c = tri[3 * e + 2];
81 const s = (xy[2 * b] - xy[2 * a]) * (xy[2 * c + 1] - xy[2 * a + 1]) - (xy[2 * c] - xy[2 * a]) * (xy[2 * b + 1] - xy[2 * a + 1]);
82 if (s < 0) { tri[3 * e + 1] = c; tri[3 * e + 2] = b; }
83 }
84 return makeMesh(xy, tri, prevIdx);
85}
87export function makeMesh(xy: Float64Array, tri: Uint32Array, boundary: Uint32Array): Mesh {
88 const N = xy.length / 2, E = tri.length / 3;
89 const area = new Float64Array(E), Gx = new Float64Array(3 * E), Gy = new Float64Array(3 * E), centroid = new Float64Array(2 * E);
90 for (let e = 0; e < E; e++) {
91 const a = tri[3 * e], b = tri[3 * e + 1], c = tri[3 * e + 2];
92 const x0 = xy[2 * a], y0 = xy[2 * a + 1], x1 = xy[2 * b], y1 = xy[2 * b + 1], x2 = xy[2 * c], y2 = xy[2 * c + 1];
93 const A = 0.5 * ((x1 - x0) * (y2 - y0) - (x2 - x0) * (y1 - y0));
94 if (!(A > 0)) throw new Error(`element ${e} has non-positive area`);
95 area[e] = A;
96 Gx[3 * e] = (y1 - y2) / (2 * A); Gx[3 * e + 1] = (y2 - y0) / (2 * A); Gx[3 * e + 2] = (y0 - y1) / (2 * A);
97 Gy[3 * e] = (x2 - x1) / (2 * A); Gy[3 * e + 1] = (x0 - x2) / (2 * A); Gy[3 * e + 2] = (x1 - x0) / (2 * A);
98 centroid[2 * e] = (x0 + x1 + x2) / 3; centroid[2 * e + 1] = (y0 + y1 + y2) / 3;
99 }
100 return { n: N, e: E, xy, tri, boundary, area, Gx, Gy, centroid };
101}
103/** Containing element and barycentric weights of a point; throws when outside the mesh. */
104export function locate(m: Mesh, px: number, py: number): { element: number; nodes: [number, number, number]; weights: [number, number, number] } {
105 const { xy, tri, area } = m;
106 for (let e = 0; e < m.e; e++) {
107 const a = tri[3 * e], b = tri[3 * e + 1], c = tri[3 * e + 2];
108 const x0 = xy[2 * a], y0 = xy[2 * a + 1], x1 = xy[2 * b], y1 = xy[2 * b + 1], x2 = xy[2 * c], y2 = xy[2 * c + 1];
109 const a2 = 2 * area[e];
110 const l1 = ((y1 - y2) * (px - x2) + (x2 - x1) * (py - y2)) / a2;
111 const l2 = ((y2 - y0) * (px - x2) + (x0 - x2) * (py - y2)) / a2;
112 const l3 = 1 - l1 - l2;
113 if (l1 >= -1e-9 && l2 >= -1e-9 && l3 >= -1e-9) return { element: e, nodes: [a, b, c], weights: [l1, l2, l3] };
114 }
115 throw new Error(`point (${px.toFixed(4)}, ${py.toFixed(4)}) lies outside the mesh`);
116}
118export function edgeLengths(m: Mesh): Float64Array {
119 const L = new Float64Array(3 * m.e);
120 for (let e = 0; e < m.e; e++) {
121 for (let k = 0; k < 3; k++) {
122 const a = m.tri[3 * e + k], b = m.tri[3 * e + ((k + 1) % 3)];
123 L[3 * e + k] = Math.hypot(m.xy[2 * a] - m.xy[2 * b], m.xy[2 * a + 1] - m.xy[2 * b + 1]);
124 }
125 }
126 return L;
127}