2"""
3crash_cymbal_fem.py — explicit nonlinear FEM simulation of a struck crash cymbal with bell, bow,
4taper, seeded imperfections, a hard-stick or soft-mallet exciter, a nonlinearity meter, and a
5spectral-onset diagnostic for the nonlinear upward energy cascade (gong-like "low → high" build-up).
7Model
8-----
9* Marguerre–von Kármán shallow shell: in-plane displacements (u, v) + transverse w on a
10 mid-surface with initial shape z0(x, y) and spatially varying thickness h(x, y). The FE lives
11 on the flat projected domain.
12* Geometry (z0_of / thickness_of): bow + bell + taper, plus
13 - hammering: N_DENTS random Gaussian dents on the bow, both signs,
14 - rim warp: low-order azimuthal harmonics growing from the shoulder to the rim,
15 - shoulder ovality (m = 2 modulation of the bell/bow junction) and a bell-centre offset.
16 All seeded via SEED, all differentiable, so ∇z0, slope check and rigid modes follow.
17* Membrane: constant-strain triangle (CST) with the full quadratic Marguerre–von Kármán strain.
18 NONLINEAR = False removes the quadratic terms (linear shallow shell, same geometry) for
19 reference runs.
20* Bending: Discrete Kirchhoff Triangle (DKT, Batoz–Bathe–Ho 1980), nodal DOFs (w, θx, θy),
21 θx = w,y θy = −w,x.
22* Internal force = −∇(strain energy) via jax.grad, gather → per-element → scatter.
23* Time: explicit leapfrog, lumped mass with scaled rotational inertia; dt from a rigorous
24 element-eigenvalue (Irons) bound. Sub-step rate is an integer multiple of the audio rate.
25* Losses: σ0 (frequency-independent) + σ1∇² (∝ ω) Rayleigh-type damping.
26* Support: free–free, 6 rigid-body modes of the curved shell projected out each step.
27* Exciter: rigid mass contacting through a patch of power-law springs F = k δ^α:
28 "stick" — Hertz (α = 1.5), hard, ~0.3 ms contact, along a line (shoulder crash);
29 "mallet" — felt (α ≈ 2.5), soft, ~3 ms contact, on a disk; injects only low frequencies,
30 so anything above a few hundred Hz must come from the nonlinear cascade.
31* Radiation at a listener point, retarded time at sub-step resolution:
32 (a) Rayleigh integral (baffled),
33 (b) unbaffled approximation: front Rayleigh − back Rayleigh routed via the shortest
34 path over the rim (Kirchhoff-style edge diffraction; NOT a BEM solve).
35* Diagnostics: membrane/bending energy split and max|w|/h per block (how nonlinear the run is);
36 per-octave onset times, spectral centroid vs time, spectrogram + band-envelope PNG.
38Deps: pip install "jax[cuda12]" gmsh scipy numpy (matplotlib optional, for the plots)
39"""
40import math, time
41import numpy as np
42import jax
43jax.config.update("jax_enable_x64", True) # set-up in f64; sim in DTYPE
44jax.config.update("jax_default_matmul_precision", "highest") # no TF32 in the stiffness matvec
45import jax.numpy as jnp
46from jax import lax
47import gmsh
48from scipy.signal import resample_poly
49from scipy.io import wavfile
51# ============================================================================ parameters
52DTYPE = jnp.float32 # simulation precision
54# material — B20 bronze
55E_MOD = 1.1e11 # Young's modulus [Pa]
56NU = 0.34
57RHO = 8700.0 # [kg/m^3]
58SIGMA0 = 0.9 # [1/s] T60 of lows ≈ 6.9/σ0
59SIGMA1 = 5e-4 # [m^2/s] ∝ω loss (T60 ≈ 0.6 s @ 5 kHz); sets the cascade ceiling
61# cymbal profile — 18" medium-thin crash, mean rim height z0 = 0
62R = 0.2286 # radius [m]
63BELL_H = 0.013 # bell height above the bow at the shoulder [m]
64BELL_R_TOP = 0.025 # radius of the (flat) bell plateau [m]
65BELL_R_SHOULDER = 0.080 # mean bell/bow junction radius [m]
66BOW_H = 0.030 # bow rise from rim to shoulder [m]
67BOW_P = 1.4 # bow power: 1 = cone, >1 = curved and flat at the shoulder
68H_RIM = 1.0e-3 # thickness at the rim [m]
69H_SHOULDER = 1.7e-3 # thickness at the bell/bow junction [m]
70H_BELL = 2.6e-3 # thickness in the bell [m]
71MAX_SLOPE = 0.4 # shallow-shell sanity bound on |∇z0|
73# symmetry breaking (all seeded)
74SEED = 1
75N_DENTS = 400 # hammer marks on the bow
76DENT_DEPTH = 0.25e-3 # typical dent depth [m], random sign (hammered from both sides)
77DENT_SIGMA = 7e-3 # dent Gaussian σ [m], ±30 % jitter (≈ 2 elements: keep ≥ H_ELEM)
78DENT_R_MIN, DENT_R_MAX = BELL_R_SHOULDER + 0.01, R - 0.01
79RIM_WARP = 0.6e-3 # m = 2 rim-warp amplitude [m]; mode m gets 2·RIM_WARP/m
80RIM_WARP_MODES = (2, 3, 4, 5)
81SHOULDER_OVALITY = 0.03 # relative m = 2 modulation of the shoulder radius
82BELL_OFFSET = 1.5e-3 # bell-centre offset from the geometric centre [m]
84# mesh
85H_ELEM = 6e-3 # target element size [m]
87# excitation
88EXCITER = "mallet" # "stick": hard shoulder crash; "mallet": soft, gong-like
89STRIKE_R, STRIKE_PHI = 0.60 * R, 0.3 # contact centre [m], [rad] (0.85 R for the stick)
90STICK_GAP = 1e-4 # initial clearance above the highest patch point [m]
91# stick (EXCITER = "stick")
92V_STRIKE, M_STICK, K_HERTZ = 7.0, 0.030, 5e8 # [m/s], [kg], total Hertz k [N/m^1.5]
93PATCH_LEN, PATCH_POINTS, PATCH_ANGLE = 0.03, 7, 0.6 # contact line [m], springs, dir. vs radial
94# mallet (EXCITER = "mallet") — power-law felt contact F = k δ^α, k set from target compression
95V_MALLET, M_MALLET = 8.0, 0.30 # [m/s], head + effective arm mass [kg] → ~10 J, tam-tam stroke
96ALPHA_MALLET = 2.5 # felt exponent (1.5 = Hertz, 2–3 = felt)
97MALLET_COMPRESSION = 12e-3 # felt compression at full stop on a rigid wall [m] → ~3 ms contact
98MALLET_RADIUS = 0.012 # head contact radius [m], spread over 7 springs
100NONLINEAR = True # False → linear shallow shell (same geometry); reference run
102# listener & air (mean rim plane is z = 0, centre at origin; drummer's ear-ish)
103LISTENER = np.array([0.40, 0.30, 0.90])
104C_AIR, RHO_AIR = 343.0, 1.2
105BACK_DIFFRACTION_GAIN = 1.0 # weight of the diffracted back-side term
107# time / numerics
108FS = 44100
109DURATION = 2.0 # [s]
110CFL_SAFETY = 0.6 # margin for nonlinearity + damping (lower to 0.5 if energy grows)
111ROT_INERTIA_FACTOR = 1 / 64 # lumped rotational inertia I = m (h²/12 + f·A_e)
112INPLANE_MASS_SCALE = 1.0 # >1 slows in-plane waves, relaxes dt (see notes)
113STEPS_PER_BLOCK = 20000
114OUT_PREFIX = "crash"
116# ============================================================================ seeded imperfections
117_rng = np.random.default_rng(SEED)
119def _uniform_annulus(n, r0, r1):
120 r = np.sqrt(_rng.uniform(r0 ** 2, r1 ** 2, n))
121 t = _rng.uniform(0, 2 * np.pi, n)
122 return r * np.cos(t), r * np.sin(t)
124_dx, _dy = _uniform_annulus(N_DENTS, DENT_R_MIN, DENT_R_MAX)
125DENT_X, DENT_Y = jnp.asarray(_dx), jnp.asarray(_dy)
126DENT_A = jnp.asarray(DENT_DEPTH * _rng.choice([-1.0, 1.0], N_DENTS) * _rng.uniform(0.5, 1.5, N_DENTS))
127DENT_S = jnp.asarray(DENT_SIGMA * _rng.uniform(0.7, 1.3, N_DENTS))
128WARP_M = jnp.asarray(RIM_WARP_MODES, dtype=float)
129WARP_A = jnp.asarray([RIM_WARP * 2.0 / m for m in RIM_WARP_MODES])
130WARP_PH = jnp.asarray(_rng.uniform(0, 2 * np.pi, len(RIM_WARP_MODES)))
131OVAL_PH = float(_rng.uniform(0, 2 * np.pi))
132_bo = _rng.uniform(0, 2 * np.pi)
133BELL_CX, BELL_CY = BELL_OFFSET * math.cos(_bo), BELL_OFFSET * math.sin(_bo)
135# ============================================================================ geometry hooks
136def _polar(x, y):
137 r = jnp.sqrt(x * x + y * y + 1e-30)
138 phi = jnp.arctan2(y, jnp.where(r < 1e-6, 1e-6, x)) # finite gradient at the centre
139 return r, phi
141def _shoulder_radius(phi):
142 return BELL_R_SHOULDER * (1.0 + SHOULDER_OVALITY * jnp.cos(2 * phi + OVAL_PH))
144def _bell_blend(r_bell, r_sh):
145 """1 on the bell plateau, 0 beyond the shoulder, C¹ cosine ramp in between."""
146 t = jnp.clip((r_bell - BELL_R_TOP) / (r_sh - BELL_R_TOP), 0.0, 1.0)
147 return 0.5 * (1.0 + jnp.cos(jnp.pi * t))
149def z0_of(x, y):
150 """Mid-surface height z0(x, y) [m]. jnp-differentiable."""
151 r, phi = _polar(x, y)
152 r_sh = _shoulder_radius(phi)
153 r_bell = jnp.sqrt((x - BELL_CX) ** 2 + (y - BELL_CY) ** 2 + 1e-30)
154 s = jnp.clip((R - r) / (R - r_sh), 0.0, 1.0) # 0 at rim → 1 at shoulder
155 u = jnp.maximum(1.0 - s, 1e-12) # 1 at rim → 0 at shoulder
156 bow = BOW_H * (1.0 - u ** BOW_P)
157 bell = BELL_H * _bell_blend(r_bell, r_sh)
158 warp = u * jnp.sum(WARP_A * jnp.cos(WARP_M * phi + WARP_PH))
159 dents = jnp.sum(DENT_A * jnp.exp(-((x - DENT_X) ** 2 + (y - DENT_Y) ** 2) / (2 * DENT_S ** 2)))
160 return bow + bell + warp + dents
162def thickness_of(x, y):
163 """Thickness h(x, y) [m]: linear taper along the bow, smooth step up into the bell."""
164 r, phi = _polar(x, y)
165 r_sh = _shoulder_radius(phi)
166 r_bell = jnp.sqrt((x - BELL_CX) ** 2 + (y - BELL_CY) ** 2 + 1e-30)
167 s = jnp.clip((R - r) / (R - r_sh), 0.0, 1.0)
168 bow = H_RIM + (H_SHOULDER - H_RIM) * s
169 return bow + (H_BELL - H_SHOULDER) * _bell_blend(r_bell, r_sh)
171_z0_grad = jax.vmap(jax.grad(lambda p: z0_of(p[0], p[1]))) # analytic ∇z0
172_z0_eval = jax.vmap(z0_of)
173_h_eval = jax.vmap(thickness_of)
175# ============================================================================ mesh (gmsh)
176def make_disk_mesh(radius, h):
177 gmsh.initialize()
178 gmsh.option.setNumber("General.Terminal", 0)
179 gmsh.model.add("cymbal")
180 gmsh.model.occ.addDisk(0, 0, 0, radius, radius)
181 gmsh.model.occ.synchronize()
182 gmsh.option.setNumber("Mesh.MeshSizeMin", 0.8 * h)
183 gmsh.option.setNumber("Mesh.MeshSizeMax", h)
184 gmsh.option.setNumber("Mesh.Algorithm", 6) # Frontal-Delaunay
185 gmsh.option.setNumber("Mesh.Smoothing", 10)
186 gmsh.model.mesh.generate(2)
187 tags, xyz, _ = gmsh.model.mesh.getNodes()
188 _, tri_tags = gmsh.model.mesh.getElementsByType(2) # 3-node triangles
189 gmsh.finalize()
190 tags = np.asarray(tags, dtype=np.int64)
191 xyz = np.asarray(xyz).reshape(-1, 3)
192 remap = np.full(tags.max() + 1, -1, dtype=np.int64)
193 remap[tags] = np.arange(len(tags))
194 tri = remap[np.asarray(tri_tags, dtype=np.int64).reshape(-1, 3)]
195 used = np.unique(tri)
196 xy = xyz[used, :2]
197 inv = np.full(len(tags), -1, dtype=np.int64); inv[used] = np.arange(len(used))
198 tri = inv[tri]
199 a = ((xy[tri[:, 1], 0] - xy[tri[:, 0], 0]) * (xy[tri[:, 2], 1] - xy[tri[:, 0], 1])
200 - (xy[tri[:, 2], 0] - xy[tri[:, 0], 0]) * (xy[tri[:, 1], 1] - xy[tri[:, 0], 1]))
201 tri[a < 0] = tri[a < 0][:, [0, 2, 1]] # counter-clockwise
202 return xy, tri
204# ============================================================================ DKT element
205def dkt_stiffness(xe, ye, D, nu):
206 """9x9 DKT bending stiffness, DOFs (w1,θx1,θy1, w2,...), θx = w,y, θy = −w,x."""
207 i, j = jnp.array([1, 2, 0]), jnp.array([2, 0, 1]) # sides k=4,5,6 ↔ ij = 23, 31, 12
208 xij, yij = xe[i] - xe[j], ye[i] - ye[j]
209 l2 = xij ** 2 + yij ** 2
210 a = -xij / l2
211 b = 0.75 * xij * yij / l2
212 c = (0.25 * xij ** 2 - 0.5 * yij ** 2) / l2
213 d = -yij / l2
214 e = (0.25 * yij ** 2 - 0.5 * xij ** 2) / l2
215 a4, a5, a6 = a; b4, b5, b6 = b; c4, c5, c6 = c; d4, d5, d6 = d; e4, e5, e6 = e
217 def H(s):
218 xi, eta = s
219 lam = 1 - xi - eta
220 N1, N2, N3 = lam * (2 * lam - 1), xi * (2 * xi - 1), eta * (2 * eta - 1)
221 N4, N5, N6 = 4 * xi * eta, 4 * eta * lam, 4 * xi * lam
222 Hx = jnp.stack([1.5 * (a6 * N6 - a5 * N5), b5 * N5 + b6 * N6, N1 - c5 * N5 - c6 * N6,
223 1.5 * (a4 * N4 - a6 * N6), b6 * N6 + b4 * N4, N2 - c4 * N4 - c6 * N6,
224 1.5 * (a5 * N5 - a4 * N4), b4 * N4 + b5 * N5, N3 - c4 * N4 - c5 * N5])
225 Hy = jnp.stack([1.5 * (d6 * N6 - d5 * N5), -N1 + e5 * N5 + e6 * N6, -b5 * N5 - b6 * N6,
226 1.5 * (d4 * N4 - d6 * N6), -N2 + e4 * N4 + e6 * N6, -b6 * N6 - b4 * N4,
227 1.5 * (d5 * N5 - d4 * N4), -N3 + e4 * N4 + e5 * N5, -b4 * N4 - b5 * N5])
228 return jnp.stack([Hx, Hy]) # (2, 9)
230 x21, x31 = xe[1] - xe[0], xe[2] - xe[0]
231 y21, y31 = ye[1] - ye[0], ye[2] - ye[0]
232 det = x21 * y31 - x31 * y21 # = 2A
233 xi_x, xi_y, eta_x, eta_y = y31 / det, -x31 / det, -y21 / det, x21 / det
235 def B_at(s):
236 J = jax.jacfwd(H)(s) # (2, 9, 2)
237 Hx_xi, Hx_eta, Hy_xi, Hy_eta = J[0, :, 0], J[0, :, 1], J[1, :, 0], J[1, :, 1]
238 Hx_x = Hx_xi * xi_x + Hx_eta * eta_x; Hx_y = Hx_xi * xi_y + Hx_eta * eta_y
239 Hy_x = Hy_xi * xi_x + Hy_eta * eta_x; Hy_y = Hy_xi * xi_y + Hy_eta * eta_y
240 return jnp.stack([Hx_x, Hy_y, Hx_y + Hy_x]) # curvature operator (3, 9)
242 Db = D * jnp.array([[1, nu, 0], [nu, 1, 0], [0, 0, 0.5 * (1 - nu)]])
243 gps = jnp.array([[1 / 6, 1 / 6], [2 / 3, 1 / 6], [1 / 6, 2 / 3]]) # exact for quadratic B
244 Bs = jax.vmap(B_at)(gps) # (3, 3, 9)
245 return (0.5 * det / 3) * jnp.einsum("gki,kl,glj->ij", Bs, Db, Bs)
247# ============================================================================ element energy
248def elem_energy_parts(qe, Gx, Gy, A, z0x, z0y, Cm, Ke):
249 """qe: (3,5) nodal [u,v,w,θx,θy]. Returns (membrane, bending) energy of one element.
250 The z0 terms are the *linear* curvature coupling; NONLINEAR gates only the quadratic terms."""
251 nl = 1.0 if NONLINEAR else 0.0
252 ux, uy = Gx @ qe[:, 0], Gy @ qe[:, 0]
253 vx, vy = Gx @ qe[:, 1], Gy @ qe[:, 1]
254 wx, wy = Gx @ qe[:, 2], Gy @ qe[:, 2]
255 exx = ux + nl * 0.5 * wx * wx + z0x * wx
256 eyy = vy + nl * 0.5 * wy * wy + z0y * wy
257 gxy = uy + vx + nl * wx * wy + z0x * wy + z0y * wx
258 C11, C12, C33 = Cm
259 Em = 0.5 * A * (C11 * (exx * exx + eyy * eyy) + 2 * C12 * exx * eyy + C33 * gxy * gxy)
260 U = qe[:, 2:].reshape(9)
261 return Em, 0.5 * (U @ Ke @ U)
263def elem_energy(*args):
264 Em, Eb = elem_energy_parts(*args)
265 return Em + Eb
267# ============================================================================ diagnostics
268def spectral_onset(p, fs, label, bands=(60, 125, 250, 500, 1000, 2000, 4000, 8000, 16000),
269 t_max=0.5):
270 """Per-band energy envelopes at ~6 ms resolution; onset = first time a band reaches −20 dB
271 of its own peak. Also plots the envelopes (log-log), which shows the cascade far better than
272 a spectrogram normalised to the low-frequency peak."""
273 from scipy.signal import stft
274 f_, t_, Z = stft(p, fs, nperseg=256, noverlap=224)
275 Pw = np.abs(Z) ** 2
276 tot_max = Pw.sum(0).max()
277 envs = []
278 print(f"\n[{label}] band onset(-20 dB) peak level vs loudest band")
279 for lo, hi in zip(bands[:-1], bands[1:]):
280 m = (f_ >= lo) & (f_ < hi)
281 env = Pw[m].sum(0)
282 if env.max() <= 0:
283 continue
284 env_db = 10 * np.log10(env / env.max() + 1e-30)
285 onset = t_[np.argmax(env_db > -20)]
286 print(f" {lo:5d}–{hi:5d} Hz {onset * 1e3:7.1f} ms {t_[np.argmax(env)] * 1e3:7.1f} ms"
287 f" {10 * np.log10(env.max() / tot_max):6.1f} dB")
288 envs.append((lo, hi, env))
289 cen = (Pw * f_[:, None]).sum(0) / (Pw.sum(0) + 1e-30)
290 for tq in (0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5):
291 if tq < t_[-1]:
292 print(f" spectral centroid at {tq * 1e3:5.0f} ms: {cen[np.searchsorted(t_, tq)]:7.0f} Hz")
293 try:
294 import matplotlib; matplotlib.use("Agg"); import matplotlib.pyplot as plt
295 fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))
296 sel = t_ <= t_max
297 ax[0].pcolormesh(t_[sel] * 1e3, f_, 10 * np.log10(Pw[:, sel] / Pw.max() + 1e-12),
298 vmin=-90, vmax=0, shading="auto")
299 ax[0].set_yscale("log"); ax[0].set_ylim(50, fs / 2)
300 ax[0].set_xlabel("t [ms]"); ax[0].set_ylabel("f [Hz]"); ax[0].set_title(label)
301 for lo, hi, env in envs:
302 ax[1].plot(t_[sel] * 1e3, 10 * np.log10(env[sel] / tot_max + 1e-30), label=f"{lo}–{hi} Hz")
303 ax[1].set_ylim(-90, 5); ax[1].set_xscale("log"); ax[1].set_xlim(2, t_max * 1e3)
304 ax[1].set_xlabel("t [ms]"); ax[1].set_ylabel("band energy [dB re loudest]")
305 ax[1].legend(fontsize=7, ncol=2); ax[1].grid(alpha=.3)
306 ax[1].set_title("band envelopes — a cascade = staggered rises, high bands last")
307 plt.tight_layout(); plt.savefig(f"{label}_cascade.png", dpi=120); plt.close()
308 print(f" wrote {label}_cascade.png")
309 except Exception:
310 pass
312# ============================================================================ main
313def main():
314 t_start = time.time()
315 tag = OUT_PREFIX if NONLINEAR else OUT_PREFIX + "_linear"
316 xy, tri = make_disk_mesh(R, H_ELEM)
317 N, E = len(xy), len(tri)
318 print(f"mesh: {N} nodes, {E} triangles, {5 * N} DOFs (NONLINEAR = {NONLINEAR})")
320 # ---- element geometry on the projected domain (numpy, f64)
321 x, y = xy[tri, 0], xy[tri, 1] # (E, 3)
322 x21, x31 = x[:, 1] - x[:, 0], x[:, 2] - x[:, 0]
323 y21, y31 = y[:, 1] - y[:, 0], y[:, 2] - y[:, 0]
324 area = 0.5 * (x21 * y31 - x31 * y21)
325 assert (area > 0).all()
326 Gx = np.stack([y[:, 1] - y[:, 2], y[:, 2] - y[:, 0], y[:, 0] - y[:, 1]], 1) / (2 * area[:, None])
327 Gy = np.stack([x[:, 2] - x[:, 1], x[:, 0] - x[:, 2], x[:, 1] - x[:, 0]], 1) / (2 * area[:, None])
328 xc, yc = x.mean(1), y.mean(1)
330 def locate(px, py):
331 """Containing element and barycentric weights of a point on the projected domain."""
332 l1 = ((y[:, 1] - y[:, 2]) * (px - x[:, 2]) + (x[:, 2] - x[:, 1]) * (py - y[:, 2])) / (2 * area)
333 l2 = ((y[:, 2] - y[:, 0]) * (px - x[:, 2]) + (x[:, 0] - x[:, 2]) * (py - y[:, 2])) / (2 * area)
334 l3 = 1 - l1 - l2
335 inside = (l1 >= -1e-9) & (l2 >= -1e-9) & (l3 >= -1e-9)
336 assert inside.any(), f"point ({px:.3f}, {py:.3f}) outside mesh"
337 es = int(np.argmax(inside))
338 return tri[es], np.array([l1[es], l2[es], l3[es]])
340 # ---- profile: thickness, initial shape and its analytic gradient
341 h_e = np.asarray(_h_eval(jnp.asarray(xc), jnp.asarray(yc)))
342 h_node = np.asarray(_h_eval(jnp.asarray(xy[:, 0]), jnp.asarray(xy[:, 1])))
343 z0_node = np.asarray(_z0_eval(jnp.asarray(xy[:, 0]), jnp.asarray(xy[:, 1])))
344 grad_z0 = np.asarray(_z0_grad(jnp.stack([jnp.asarray(xc), jnp.asarray(yc)], 1)))
345 z0x, z0y = grad_z0[:, 0], grad_z0[:, 1]
346 slope = np.hypot(z0x, z0y).max()
347 rim_node = np.hypot(xy[:, 0], xy[:, 1]) > R - 1e-6
348 print(f"profile: height {z0_node.max() * 1e3:.1f} mm, rim warp p-p "
349 f"{(z0_node[rim_node].max() - z0_node[rim_node].min()) * 1e3:.2f} mm, "
350 f"{N_DENTS} dents, h ∈ [{h_e.min() * 1e3:.2f}, {h_e.max() * 1e3:.2f}] mm, "
351 f"max |∇z0| = {slope:.3f}")
352 if slope > MAX_SLOPE:
353 print(f"WARNING: max slope {slope:.2f} exceeds shallow-shell bound {MAX_SLOPE}")
355 D_e = E_MOD * h_e ** 3 / (12 * (1 - NU ** 2))
356 Cm_e = (E_MOD * h_e / (1 - NU ** 2))[:, None] * np.array([1.0, NU, 0.5 * (1 - NU)])
358 # ---- DKT stiffness for every element
359 Ke64 = jax.vmap(dkt_stiffness, in_axes=(0, 0, 0, None))(jnp.asarray(x), jnp.asarray(y),
360 jnp.asarray(D_e), NU)
361 Ke_np = np.asarray(Ke64)
363 # element self-tests: rigid tilt → zero force; w = ½x² → energy ½·D·A
364 U_tilt = np.stack([y, np.ones_like(y), np.zeros_like(y)], 2).reshape(E, 9)
365 res = np.abs(np.einsum("eij,ej->ei", Ke_np, U_tilt)).max() / np.abs(Ke_np).max()
366 U_cur = np.stack([0.5 * x ** 2, np.zeros_like(x), -x], 2).reshape(E, 9)
367 Eb = 0.5 * np.einsum("ei,eij,ej->e", U_cur, Ke_np, U_cur)
368 err = np.abs(Eb / (0.5 * D_e * area) - 1).max()
369 print(f"DKT checks: rigid-tilt residual {res:.1e}, constant-curvature energy error {err:.1e}")
371 # ---- lumped mass (per DOF): [u, v, w, θx, θy]
372 m_e = RHO * h_e * area / 3
373 I_e = m_e * (h_e ** 2 / 12 + ROT_INERTIA_FACTOR * area)
374 M = np.zeros((N, 5))
375 for k in range(3):
376 np.add.at(M[:, 0], tri[:, k], m_e * INPLANE_MASS_SCALE)
377 np.add.at(M[:, 1], tri[:, k], m_e * INPLANE_MASS_SCALE)
378 np.add.at(M[:, 2], tri[:, k], m_e)
379 np.add.at(M[:, 3], tri[:, k], I_e)
380 np.add.at(M[:, 4], tri[:, k], I_e)
381 print(f"plate mass {M[:, 2].sum():.3f} kg")
383 # ---- stable dt from the element eigenvalue bound (Irons): λ_max(M⁻¹K) ≤ max_e λ_max(M_e⁻¹K_e)
384 Khess = jax.vmap(jax.hessian(elem_energy))(
385 jnp.zeros((E, 3, 5)), jnp.asarray(Gx), jnp.asarray(Gy), jnp.asarray(area),
386 jnp.asarray(z0x), jnp.asarray(z0y), jnp.asarray(Cm_e), Ke64).reshape(E, 15, 15)
387 Me = np.tile(np.stack([m_e * INPLANE_MASS_SCALE, m_e * INPLANE_MASS_SCALE, m_e, I_e, I_e], 1),
388 (1, 3)) # (E, 15) node-major
389 S = np.asarray(Khess) / np.sqrt(Me[:, :, None] * Me[:, None, :])
390 lam = np.linalg.eigvalsh(S)
391 n_zero = int((lam[0] < 1e-9 * lam[0].max()).sum())
392 print(f"element 0 zero-energy modes: {n_zero} (expect 6)")
393 e_worst = int(np.argmax(lam[:, -1]))
394 omega_max = math.sqrt(lam[:, -1].max())
395 dt_bound = 2.0 / omega_max
396 M_SUB = math.ceil(1.0 / (FS * CFL_SAFETY * dt_bound))
397 dt = 1.0 / (FS * M_SUB)
398 n_blocks = math.ceil(DURATION * FS * M_SUB / STEPS_PER_BLOCK)
399 n_steps = n_blocks * STEPS_PER_BLOCK
400 print(f"ω_max bound {omega_max:.3e} rad/s (limiting element at r = "
401 f"{math.hypot(xc[e_worst], yc[e_worst]):.3f} m, h = {h_e[e_worst] * 1e3:.2f} mm)")
402 print(f"dt = {dt:.3e} s, {M_SUB} sub-steps/sample, {n_steps} steps")
404 # ---- rigid-body modes of the curved shell, M-orthonormalised
405 # A rigid tilt drags the in-plane displacements along: rotation about x by φ gives
406 # v = −φ z0, w = φ y; about y gives u = +φ z0, w = −φ x (first order, shallow-shell).
407 X, Y = xy[:, 0] - xy[:, 0].mean(), xy[:, 1] - xy[:, 1].mean()
408 Rm = np.zeros((6, N, 5))
409 Rm[0, :, 0] = 1; Rm[1, :, 1] = 1 # in-plane translations
410 Rm[2, :, 0] = -Y; Rm[2, :, 1] = X # in-plane rotation
411 Rm[3, :, 2] = 1 # z translation
412 Rm[4, :, 2] = Y; Rm[4, :, 3] = 1; Rm[4, :, 1] = -z0_node # tilt about x (θx = w,y)
413 Rm[5, :, 2] = -X; Rm[5, :, 4] = 1; Rm[5, :, 0] = z0_node # tilt about y (θy = −w,x)
414 for k in range(6):
415 for l in range(k):
416 Rm[k] -= (Rm[l] * M * Rm[k]).sum() * Rm[l]
417 Rm[k] /= math.sqrt((Rm[k] * M * Rm[k]).sum())
419 # ---- contact patch: a set of power-law springs F = k δ^α on a rigid exciter of mass M_EX
420 xs, ys = STRIKE_R * math.cos(STRIKE_PHI), STRIKE_R * math.sin(STRIKE_PHI)
421 if EXCITER == "stick":
422 M_EX, V_EX, ALPHA, K_TOTAL = M_STICK, V_STRIKE, 1.5, K_HERTZ
423 ang = STRIKE_PHI + PATCH_ANGLE
424 xi = np.linspace(-0.5, 0.5, PATCH_POINTS) * PATCH_LEN
425 px, py = xs + xi * math.cos(ang), ys + xi * math.sin(ang)
426 elif EXCITER == "mallet":
427 M_EX, V_EX, ALPHA = M_MALLET, V_MALLET, ALPHA_MALLET
428 # k such that the head stops after MALLET_COMPRESSION on a rigid wall
429 K_TOTAL = (ALPHA + 1) * 0.5 * M_EX * V_EX ** 2 / MALLET_COMPRESSION ** (ALPHA + 1)
430 th = np.arange(6) * np.pi / 3
431 px = np.concatenate([[xs], xs + MALLET_RADIUS * np.cos(th)])
432 py = np.concatenate([[ys], ys + MALLET_RADIUS * np.sin(th)])
433 else:
434 raise ValueError(f"unknown EXCITER {EXCITER!r}")
435 P = len(px)
436 loc = [locate(px[k], py[k]) for k in range(P)]
437 s_idx = np.array([l[0] for l in loc]) # (P, 3)
438 s_wts = np.array([l[1] for l in loc]) # (P, 3)
439 z0_p = np.asarray(_z0_eval(jnp.asarray(px), jnp.asarray(py)))
440 # rigid exciter face: best-fit plane through the surface heights; residual → staggered contact
441 Apl = np.column_stack([px - xs, py - ys, np.ones(P)])
442 s_off = z0_p - Apl @ np.linalg.lstsq(Apl, z0_p, rcond=None)[0]
443 zs0 = s_off.max() + STICK_GAP
444 tau_est = 2.0 * MALLET_COMPRESSION / V_EX if EXCITER == "mallet" else None
445 print(f"{EXCITER}: {P}-point patch, k = {K_TOTAL:.3g} N/m^{ALPHA}, m = {M_EX} kg, v = {V_EX} m/s, "
446 f"energy {0.5 * M_EX * V_EX ** 2:.2f} J"
447 + (f", free contact ≈ {tau_est * 1e3:.1f} ms → direct excitation below ~{0.5 / tau_est:.0f} Hz"
448 if tau_est else ""))
450 # ---- radiation geometry (element centroids on the curved surface, retarded times in sub-steps)
451 cent = np.stack([xc, yc, z0_node[tri].mean(1)], 1) # (E, 3)
452 r_front = np.linalg.norm(LISTENER - cent, axis=1)
453 phi = np.linspace(0, 2 * np.pi, 720, endpoint=False)
454 rim_xy = np.stack([R * np.cos(phi), R * np.sin(phi)], 1)
455 rim_z = np.asarray(_z0_eval(jnp.asarray(rim_xy[:, 0]), jnp.asarray(rim_xy[:, 1])))
456 rim = np.column_stack([rim_xy, rim_z])
457 d1 = np.linalg.norm(cent[:, None, :] - rim[None, :, :], axis=2)
458 d2 = np.linalg.norm(rim - LISTENER, axis=1)
459 r_back = (d1 + d2[None, :]).min(1)
460 coef_f = RHO_AIR * area / (2 * np.pi * r_front)
461 coef_b = BACK_DIFFRACTION_GAIN * RHO_AIR * area / (2 * np.pi * r_back)
462 df, db = r_front / (C_AIR * dt), r_back / (C_AIR * dt)
463 i0f, frf = np.floor(df).astype(np.int32), df - np.floor(df)
464 i0b, frb = np.floor(db).astype(np.int32), db - np.floor(db)
465 buf_len = n_steps + int(max(i0f.max(), i0b.max())) + 2
467 # ---- device arrays in simulation precision
468 f = lambda a: jnp.asarray(a, dtype=DTYPE)
469 tri_j = jnp.asarray(tri, dtype=jnp.int32)
470 Gx_j, Gy_j, area_j = f(Gx), f(Gy), f(area)
471 z0x_j, z0y_j, Cm_j, Ke_j = f(z0x), f(z0y), f(Cm_e), f(Ke_np)
472 rhoh_area_j, hnode_j = f(RHO * h_e * area), f(h_node)
473 M_j, Minv_j, Rm_j = f(M), f(1.0 / M), f(Rm)
474 sidx_j, swts_j, soff_j = jnp.asarray(s_idx, dtype=jnp.int32), f(s_wts), f(s_off)
475 coef_f_j, coef_b_j, frf_j, frb_j = f(coef_f), f(coef_b), f(frf), f(frb)
476 i0f_j, i0b_j = jnp.asarray(i0f), jnp.asarray(i0b)
477 k_point = K_TOTAL / P
479 def total_energy(q):
480 return jax.vmap(elem_energy)(q[tri_j], Gx_j, Gy_j, area_j, z0x_j, z0y_j, Cm_j, Ke_j).sum()
482 f_int = jax.grad(lambda q: -total_energy(q)) # internal nodal forces
484 # shell check: tangent stiffness at rest must annihilate all six rigid modes
485 q0 = jnp.zeros((N, 5), DTYPE)
486 K_apply = lambda d: jax.jvp(f_int, (q0,), (d,))[1]
487 chk = np.random.default_rng(0)
488 worst = 0.0
489 for k in range(6):
490 rnd = chk.standard_normal((N, 5)); rnd *= np.linalg.norm(Rm[k]) / np.linalg.norm(rnd)
491 worst = max(worst, float(jnp.abs(K_apply(Rm_j[k])).max() / jnp.abs(K_apply(f(rnd))).max()))
492 print(f"rigid modes: max |K·r| / |K·random| = {worst:.1e} (should be ~1e-6 or smaller)")
494 def damp_lap(vw):
495 """−2σ1 ρh ∇²ẇ as a nodal force: 2σ1 · (ρh A_e) · Σ_e G_iᵀ (G ẇ_e)."""
496 ve = vw[tri_j]
497 vx, vy = (Gx_j * ve).sum(1), (Gy_j * ve).sum(1)
498 c = rhoh_area_j[:, None] * (Gx_j * vx[:, None] + Gy_j * vy[:, None])
499 return 2 * SIGMA1 * jax.ops.segment_sum(c.ravel(), tri_j.ravel(), num_segments=N)
501 def project_rigid(v):
502 coef = jnp.einsum("knd,nd->k", Rm_j, M_j * v)
503 return v - jnp.einsum("k,knd->nd", coef, Rm_j)
505 def step(n, st):
506 q, v, zs, vs, bf, bb, Fmax, ncont = st
507 # exciter / plate contact along the patch
508 wp = jnp.einsum("pk,pk->p", swts_j, q[sidx_j, 2]) + soff_j # surface vs exciter face
509 delta = jnp.maximum(wp - zs, 0.0) # (P,)
510 Fk = k_point * jnp.power(delta, ALPHA)
511 F = Fk.sum()
512 # forces
513 fo = f_int(q)
514 fo = fo.at[sidx_j, 2].add(-Fk[:, None] * swts_j)
515 fo = fo.at[:, 2].add(-damp_lap(v[:, 2]))
516 # leapfrog (σ0 implicit), rigid-mode projection
517 v1 = (v + dt * fo * Minv_j) / (1 + 2 * SIGMA0 * dt)
518 v1 = project_rigid(v1)
519 q1 = q + dt * v1
520 vs1 = vs + dt * F / M_EX
521 zs1 = zs + dt * vs1
522 # radiation: scatter element accelerations into retarded-time buffers
523 acc = ((v1[:, 2] - v[:, 2]) / dt)[tri_j].mean(1)
524 sf, sb = coef_f_j * acc, coef_b_j * acc
525 bf = bf.at[n + i0f_j].add(sf * (1 - frf_j)).at[n + i0f_j + 1].add(sf * frf_j)
526 bb = bb.at[n + i0b_j].add(sb * (1 - frb_j)).at[n + i0b_j + 1].add(sb * frb_j)
527 return (q1, v1, zs1, vs1, bf, bb, jnp.maximum(Fmax, F), ncont + (F > 0))
530 def run_block(n0, st):
531 return lax.fori_loop(0, STEPS_PER_BLOCK,
532 lambda k, s: step((n0 + k).astype(jnp.int32), s), st)
535 def energies(st):
536 q, v = st[0], st[1]
537 Em, Eb = jax.vmap(elem_energy_parts)(q[tri_j], Gx_j, Gy_j, area_j, z0x_j, z0y_j, Cm_j, Ke_j)
538 return (0.5 * (M_j * v * v).sum(), Em.sum(), Eb.sum(), 0.5 * M_EX * st[3] ** 2,
539 jnp.abs(q[:, 2] / hnode_j).max())
541 state = (jnp.zeros((N, 5), DTYPE), jnp.zeros((N, 5), DTYPE),
542 jnp.asarray(zs0, DTYPE), jnp.asarray(-V_EX, DTYPE),
543 jnp.zeros(buf_len, DTYPE), jnp.zeros(buf_len, DTYPE),
544 jnp.asarray(0.0, DTYPE), jnp.asarray(0, jnp.int32))
546 print(f"setup done in {time.time() - t_start:.1f} s; compiling + running on "
547 f"{jax.devices()[0].platform} ...")
548 print(" (nonlinearity meter: max|w|/h ≲ 1 → effectively linear; 3–5 → clear cascade)")
549 t_run = time.time()
550 for blk in range(n_blocks):
551 state = run_block(jnp.asarray(blk * STEPS_PER_BLOCK, jnp.int32), state)
552 if blk % max(1, n_blocks // 20) == 0 or blk == n_blocks - 1:
553 ke, em, eb, ks, woh = energies(state)
554 jax.block_until_ready(ke)
555 print(f" t = {(blk + 1) * STEPS_PER_BLOCK * dt:6.3f} s KE {float(ke):.3e} "
556 f"membrane {float(em):.3e} bending {float(eb):.3e} J exciter {float(ks):.3e} J "
557 f"max|w|/h = {float(woh):5.2f} [{time.time() - t_run:.0f} s wall]")
559 q, v, zs, vs, bf, bb, Fmax, ncont = jax.device_get(state)
560 print(f"{EXCITER}: peak force {float(Fmax):.1f} N, contact {float(ncont) * dt * 1e3:.2f} ms, "
561 f"rebound velocity {float(vs):+.2f} m/s "
562 f"(energy left in exciter {0.5 * M_EX * float(vs) ** 2:.2f} of {0.5 * M_EX * V_EX ** 2:.2f} J)")
564 # ---- decimate to audio rate, write, analyse
565 n_out = int(FS * DURATION)
566 p_ray = resample_poly(np.asarray(bf[:n_steps], np.float64), 1, M_SUB)[:n_out]
567 p_back = resample_poly(np.asarray(bb[:n_steps], np.float64), 1, M_SUB)[:n_out]
568 p_free = p_ray - p_back
569 for name, p in (("rayleigh", p_ray), ("unbaffled", p_free)):
570 peak = np.abs(p).max()
571 print(f"{name:10s}: peak {peak:.3f} Pa ({20 * np.log10(peak / 2e-5):.1f} dB SPL peak)")
572 wavfile.write(f"{tag}_{name}.wav", FS, (0.9 * p / peak * 32767).astype(np.int16))
573 print(f"wrote {tag}_rayleigh.wav, {tag}_unbaffled.wav (total {time.time() - t_start:.0f} s)")
574 spectral_onset(p_free, FS, f"{tag}_unbaffled")
576if __name__ == "__main__":
577 main()