1// The whole precompute: from a cymsim configuration to the modal data the player needs.
2//
3// 1. profile, mesh (at the playback element size) and shell, as in the simulator;
4// 2. static condensation of the in-plane and rotation DOFs (assemble.ts);
5// 3. mass-normalised eigenmodes of the condensed transverse problem (eigen.ts);
6// 4. per mode: the frequency, the quadratic form of the sigma1 loss, the complex listener gains
7// (gains.ts) and the mode shape on the nodes, quantised to int16.
8//
9// What is *not* baked in, so that the design panel can change it without a new eigensolve:
10//
11// damping sigma_m = sigma0 + sigma1 q_m with q_m = phi_m^T L phi_m stored as ``dampQ``.
12// q_m is invariant under a density change, because L is proportional to rho while a
13// mass-normalised phi is proportional to 1/sqrt(rho).
14// material the condensed stiffness is proportional to E and the transverse mass to rho, both
15// exactly, so a change of either scales every frequency by sqrt((E/E0)/(rho/rho0))
16// and leaves the mode shapes unchanged apart from the normalisation factor
17// 1/sqrt(rho/rho0). Poisson's ratio does change the shapes and needs a recompute.
18// listener the gains are a sum over elements and can be redone in a fraction of a second.
19//
20// Everything else (radius, thickness, bell and bow geometry, imperfections, Poisson's ratio,
21// mesh size) changes the eigenproblem and needs the full computation again.
22import { cloneConfig } from "../config/yaml.ts";
23import { defaultMallet, defaultStick, type SimulationConfig } from "../config/schema.ts";
24import { CymbalProfile } from "../model/geometry.ts";
25import { diskMesh } from "../model/mesh.ts";
26import { buildShell } from "../model/shell.ts";
27import { FORMAT_VERSION, type ModesData, type ModesMeta } from "../modes.ts";
28import { condensedStiffness, type Progress } from "./assemble.ts";
29import { symmetricEigen } from "./eigen.ts";
30import { computeGains } from "./gains.ts";
32export interface PrecomputeOptions {
33 elementSize: number; // playback mesh element size [m]
34 fmax: number; // keep modes below this frequency [Hz]
35}
37/** The parts of a configuration that the eigenproblem depends on. Two configurations with the
38 * same structural key share a modes file; the rest is applied when the file is loaded. */
39export function structuralKey(cfg: SimulationConfig, opts: PrecomputeOptions) {
40 return {
41 profile: cfg.profile,
42 imperfections: cfg.imperfections,
43 smoothing: cfg.mesh.smoothing,
44 nu: cfg.material.nu,
45 options: opts,
46 };
47}
49export function computeModes(baseConfig: SimulationConfig, opts: PrecomputeOptions, progress: Progress): ModesData {
50 const cfg = cloneConfig(baseConfig);
51 cfg.mesh.element_size = opts.elementSize;
52 progress("mesh and shell", 0);
53 const profile = new CymbalProfile(cfg.profile, cfg.imperfections);
54 const mesh = diskMesh(cfg.profile.radius, cfg.mesh.element_size, cfg.mesh.smoothing);
55 const shell = buildShell(mesh, profile, cfg.material, cfg.numerics);
56 const N = mesh.n, E = mesh.e, tri = mesh.tri;
58 const S = condensedStiffness(mesh, shell, progress);
60 // mass normalisation with the lumped transverse mass
61 const Mw = new Float64Array(N), dinv = new Float64Array(N);
62 for (let n = 0; n < N; n++) {
63 Mw[n] = shell.M[5 * n + 2];
64 dinv[n] = 1 / Math.sqrt(Mw[n]);
65 }
66 for (let i = 0; i < N; i++) for (let j = 0; j < N; j++) S[i * N + j] *= dinv[i] * dinv[j];
67 progress("eigenmodes", 0);
68 const eig = symmetricEigen(S, N, (phase, f) => progress(`eigenmodes: ${phase}`, f));
70 // keep 5 Hz < f < fmax (drops the three rigid modes)
71 const keep: number[] = [];
72 for (let k = 0; k < N; k++) {
73 const f = Math.sqrt(Math.max(eig.values[k], 0)) / (2 * Math.PI);
74 if (f > 5 && f < opts.fmax) keep.push(k);
75 }
76 const Nm = keep.length;
77 if (Nm === 0) throw new Error("no modes below the frequency limit: check the profile and the material");
78 const freq = new Float64Array(Nm);
79 const phi = new Float64Array(N * Nm); // node-major
80 keep.forEach((k, m) => {
81 freq[m] = Math.sqrt(eig.values[k]) / (2 * Math.PI);
82 for (let n = 0; n < N; n++) phi[n * Nm + m] = eig.vectorsT[k * N + n] * dinv[n];
83 });
84 progress("damping and radiation", 0);
86 // sigma1 quadratic form q_m = phi^T L phi with L = sum_e (rho h A)_e G^T G on w
87 const dampQ = new Float64Array(Nm);
88 const gxm = new Float64Array(Nm), gym = new Float64Array(Nm);
89 for (let e = 0; e < E; e++) {
90 gxm.fill(0);
91 gym.fill(0);
92 for (let i = 0; i < 3; i++) {
93 const row = tri[3 * e + i] * Nm, gx = mesh.Gx[3 * e + i], gy = mesh.Gy[3 * e + i];
94 for (let m = 0; m < Nm; m++) {
95 gxm[m] += gx * phi[row + m];
96 gym[m] += gy * phi[row + m];
97 }
98 }
99 const c = shell.rhoHA[e];
100 for (let m = 0; m < Nm; m++) dampQ[m] += c * (gxm[m] * gxm[m] + gym[m] * gym[m]);
101 }
103 // quantise the mode shapes
104 const phiScale = new Float32Array(Nm), phiQ = new Int16Array(N * Nm);
105 for (let m = 0; m < Nm; m++) {
106 let top = 0;
107 for (let n = 0; n < N; n++) top = Math.max(top, Math.abs(phi[n * Nm + m]));
108 phiScale[m] = top / 32767 || 1;
109 }
110 for (let n = 0; n < N; n++) for (let m = 0; m < Nm; m++) phiQ[n * Nm + m] = Math.round(phi[n * Nm + m] / phiScale[m]);
112 let mass = 0;
113 for (const mw of Mw) mass += mw;
114 const meta: ModesMeta = {
115 format_version: FORMAT_VERSION,
116 n_nodes: N, n_elements: E, n_modes: Nm, fmax: opts.fmax,
117 radius: cfg.profile.radius, mass,
118 E0: cfg.material.E, rho0: cfg.material.rho,
119 config: baseConfig as unknown as Record<string, unknown>,
120 playback_element_size: opts.elementSize,
121 exciters: { stick: defaultStick(), mallet: defaultMallet() },
122 notes: "phi = phi_q * phi_scale[m] is the transverse mode shape, mass-normalised at rho0, with " +
123 "the in-plane and rotation DOFs statically condensed. Modal equation a'' + 2 sigma a' + " +
124 "(2 pi f)^2 a = phi(p) . F, with sigma_m = sigma0 + sigma1 damp_q[m] and f scaled by " +
125 "sqrt((E/E0)/(rho/rho0)). Pressure = Re(sum_m G_m * complex acceleration); unbaffled = front - back.",
126 };
128 const data: ModesData = {
129 meta, nNodes: N, nElements: E, nModes: Nm, radius: cfg.profile.radius,
130 xy: Float32Array.from(mesh.xy), tri: Int32Array.from(mesh.tri),
131 z0: Float32Array.from(shell.z0Node), h: Float32Array.from(shell.hNode),
132 freq, dampQ,
133 gainFront: new Float32Array(2 * Nm), gainBack: new Float32Array(2 * Nm),
134 phiScale, phiQ,
135 };
136 const g = computeGains(data, cfg.listener, freq, 1, undefined, (f) => progress("damping and radiation", f));
137 data.gainFront = g.front;
138 data.gainBack = g.back;
139 progress("done", 1);
140 return data;
141}