#!/usr/bin/env python3 """ crash_cymbal_fem.py — explicit nonlinear FEM simulation of a struck crash cymbal with bell, bow, taper, seeded imperfections, a hard-stick or soft-mallet exciter, a nonlinearity meter, and a spectral-onset diagnostic for the nonlinear upward energy cascade (gong-like "low → high" build-up). Model ----- * Marguerre–von Kármán shallow shell: in-plane displacements (u, v) + transverse w on a mid-surface with initial shape z0(x, y) and spatially varying thickness h(x, y). The FE lives on the flat projected domain. * Geometry (z0_of / thickness_of): bow + bell + taper, plus - hammering: N_DENTS random Gaussian dents on the bow, both signs, - rim warp: low-order azimuthal harmonics growing from the shoulder to the rim, - shoulder ovality (m = 2 modulation of the bell/bow junction) and a bell-centre offset. All seeded via SEED, all differentiable, so ∇z0, slope check and rigid modes follow. * Membrane: constant-strain triangle (CST) with the full quadratic Marguerre–von Kármán strain. NONLINEAR = False removes the quadratic terms (linear shallow shell, same geometry) for reference runs. * Bending: Discrete Kirchhoff Triangle (DKT, Batoz–Bathe–Ho 1980), nodal DOFs (w, θx, θy), θx = w,y θy = −w,x. * Internal force = −∇(strain energy) via jax.grad, gather → per-element → scatter. * Time: explicit leapfrog, lumped mass with scaled rotational inertia; dt from a rigorous element-eigenvalue (Irons) bound. Sub-step rate is an integer multiple of the audio rate. * Losses: σ0 (frequency-independent) + σ1∇² (∝ ω) Rayleigh-type damping. * Support: free–free, 6 rigid-body modes of the curved shell projected out each step. * Exciter: rigid mass contacting through a patch of power-law springs F = k δ^α: "stick" — Hertz (α = 1.5), hard, ~0.3 ms contact, along a line (shoulder crash); "mallet" — felt (α ≈ 2.5), soft, ~3 ms contact, on a disk; injects only low frequencies, so anything above a few hundred Hz must come from the nonlinear cascade. * Radiation at a listener point, retarded time at sub-step resolution: (a) Rayleigh integral (baffled), (b) unbaffled approximation: front Rayleigh − back Rayleigh routed via the shortest path over the rim (Kirchhoff-style edge diffraction; NOT a BEM solve). * Diagnostics: membrane/bending energy split and max|w|/h per block (how nonlinear the run is); per-octave onset times, spectral centroid vs time, spectrogram + band-envelope PNG. Deps: pip install "jax[cuda12]" gmsh scipy numpy (matplotlib optional, for the plots) """ import math, time import numpy as np import jax jax.config.update("jax_enable_x64", True) # set-up in f64; sim in DTYPE jax.config.update("jax_default_matmul_precision", "highest") # no TF32 in the stiffness matvec import jax.numpy as jnp from jax import lax import gmsh from scipy.signal import resample_poly from scipy.io import wavfile # ============================================================================ parameters DTYPE = jnp.float32 # simulation precision # material — B20 bronze E_MOD = 1.1e11 # Young's modulus [Pa] NU = 0.34 RHO = 8700.0 # [kg/m^3] SIGMA0 = 0.9 # [1/s] T60 of lows ≈ 6.9/σ0 SIGMA1 = 5e-4 # [m^2/s] ∝ω loss (T60 ≈ 0.6 s @ 5 kHz); sets the cascade ceiling # cymbal profile — 18" medium-thin crash, mean rim height z0 = 0 R = 0.2286 # radius [m] BELL_H = 0.013 # bell height above the bow at the shoulder [m] BELL_R_TOP = 0.025 # radius of the (flat) bell plateau [m] BELL_R_SHOULDER = 0.080 # mean bell/bow junction radius [m] BOW_H = 0.030 # bow rise from rim to shoulder [m] BOW_P = 1.4 # bow power: 1 = cone, >1 = curved and flat at the shoulder H_RIM = 1.0e-3 # thickness at the rim [m] H_SHOULDER = 1.7e-3 # thickness at the bell/bow junction [m] H_BELL = 2.6e-3 # thickness in the bell [m] MAX_SLOPE = 0.4 # shallow-shell sanity bound on |∇z0| # symmetry breaking (all seeded) SEED = 1 N_DENTS = 400 # hammer marks on the bow DENT_DEPTH = 0.25e-3 # typical dent depth [m], random sign (hammered from both sides) DENT_SIGMA = 7e-3 # dent Gaussian σ [m], ±30 % jitter (≈ 2 elements: keep ≥ H_ELEM) DENT_R_MIN, DENT_R_MAX = BELL_R_SHOULDER + 0.01, R - 0.01 RIM_WARP = 0.6e-3 # m = 2 rim-warp amplitude [m]; mode m gets 2·RIM_WARP/m RIM_WARP_MODES = (2, 3, 4, 5) SHOULDER_OVALITY = 0.03 # relative m = 2 modulation of the shoulder radius BELL_OFFSET = 1.5e-3 # bell-centre offset from the geometric centre [m] # mesh H_ELEM = 6e-3 # target element size [m] # excitation EXCITER = "mallet" # "stick": hard shoulder crash; "mallet": soft, gong-like STRIKE_R, STRIKE_PHI = 0.60 * R, 0.3 # contact centre [m], [rad] (0.85 R for the stick) STICK_GAP = 1e-4 # initial clearance above the highest patch point [m] # stick (EXCITER = "stick") V_STRIKE, M_STICK, K_HERTZ = 7.0, 0.030, 5e8 # [m/s], [kg], total Hertz k [N/m^1.5] PATCH_LEN, PATCH_POINTS, PATCH_ANGLE = 0.03, 7, 0.6 # contact line [m], springs, dir. vs radial # mallet (EXCITER = "mallet") — power-law felt contact F = k δ^α, k set from target compression V_MALLET, M_MALLET = 8.0, 0.30 # [m/s], head + effective arm mass [kg] → ~10 J, tam-tam stroke ALPHA_MALLET = 2.5 # felt exponent (1.5 = Hertz, 2–3 = felt) MALLET_COMPRESSION = 12e-3 # felt compression at full stop on a rigid wall [m] → ~3 ms contact MALLET_RADIUS = 0.012 # head contact radius [m], spread over 7 springs NONLINEAR = True # False → linear shallow shell (same geometry); reference run # listener & air (mean rim plane is z = 0, centre at origin; drummer's ear-ish) LISTENER = np.array([0.40, 0.30, 0.90]) C_AIR, RHO_AIR = 343.0, 1.2 BACK_DIFFRACTION_GAIN = 1.0 # weight of the diffracted back-side term # time / numerics FS = 44100 DURATION = 2.0 # [s] CFL_SAFETY = 0.6 # margin for nonlinearity + damping (lower to 0.5 if energy grows) ROT_INERTIA_FACTOR = 1 / 64 # lumped rotational inertia I = m (h²/12 + f·A_e) INPLANE_MASS_SCALE = 1.0 # >1 slows in-plane waves, relaxes dt (see notes) STEPS_PER_BLOCK = 20000 OUT_PREFIX = "crash" # ============================================================================ seeded imperfections _rng = np.random.default_rng(SEED) def _uniform_annulus(n, r0, r1): r = np.sqrt(_rng.uniform(r0 ** 2, r1 ** 2, n)) t = _rng.uniform(0, 2 * np.pi, n) return r * np.cos(t), r * np.sin(t) _dx, _dy = _uniform_annulus(N_DENTS, DENT_R_MIN, DENT_R_MAX) DENT_X, DENT_Y = jnp.asarray(_dx), jnp.asarray(_dy) DENT_A = jnp.asarray(DENT_DEPTH * _rng.choice([-1.0, 1.0], N_DENTS) * _rng.uniform(0.5, 1.5, N_DENTS)) DENT_S = jnp.asarray(DENT_SIGMA * _rng.uniform(0.7, 1.3, N_DENTS)) WARP_M = jnp.asarray(RIM_WARP_MODES, dtype=float) WARP_A = jnp.asarray([RIM_WARP * 2.0 / m for m in RIM_WARP_MODES]) WARP_PH = jnp.asarray(_rng.uniform(0, 2 * np.pi, len(RIM_WARP_MODES))) OVAL_PH = float(_rng.uniform(0, 2 * np.pi)) _bo = _rng.uniform(0, 2 * np.pi) BELL_CX, BELL_CY = BELL_OFFSET * math.cos(_bo), BELL_OFFSET * math.sin(_bo) # ============================================================================ geometry hooks def _polar(x, y): r = jnp.sqrt(x * x + y * y + 1e-30) phi = jnp.arctan2(y, jnp.where(r < 1e-6, 1e-6, x)) # finite gradient at the centre return r, phi def _shoulder_radius(phi): return BELL_R_SHOULDER * (1.0 + SHOULDER_OVALITY * jnp.cos(2 * phi + OVAL_PH)) def _bell_blend(r_bell, r_sh): """1 on the bell plateau, 0 beyond the shoulder, C¹ cosine ramp in between.""" t = jnp.clip((r_bell - BELL_R_TOP) / (r_sh - BELL_R_TOP), 0.0, 1.0) return 0.5 * (1.0 + jnp.cos(jnp.pi * t)) def z0_of(x, y): """Mid-surface height z0(x, y) [m]. jnp-differentiable.""" r, phi = _polar(x, y) r_sh = _shoulder_radius(phi) r_bell = jnp.sqrt((x - BELL_CX) ** 2 + (y - BELL_CY) ** 2 + 1e-30) s = jnp.clip((R - r) / (R - r_sh), 0.0, 1.0) # 0 at rim → 1 at shoulder u = jnp.maximum(1.0 - s, 1e-12) # 1 at rim → 0 at shoulder bow = BOW_H * (1.0 - u ** BOW_P) bell = BELL_H * _bell_blend(r_bell, r_sh) warp = u * jnp.sum(WARP_A * jnp.cos(WARP_M * phi + WARP_PH)) dents = jnp.sum(DENT_A * jnp.exp(-((x - DENT_X) ** 2 + (y - DENT_Y) ** 2) / (2 * DENT_S ** 2))) return bow + bell + warp + dents def thickness_of(x, y): """Thickness h(x, y) [m]: linear taper along the bow, smooth step up into the bell.""" r, phi = _polar(x, y) r_sh = _shoulder_radius(phi) r_bell = jnp.sqrt((x - BELL_CX) ** 2 + (y - BELL_CY) ** 2 + 1e-30) s = jnp.clip((R - r) / (R - r_sh), 0.0, 1.0) bow = H_RIM + (H_SHOULDER - H_RIM) * s return bow + (H_BELL - H_SHOULDER) * _bell_blend(r_bell, r_sh) _z0_grad = jax.vmap(jax.grad(lambda p: z0_of(p[0], p[1]))) # analytic ∇z0 _z0_eval = jax.vmap(z0_of) _h_eval = jax.vmap(thickness_of) # ============================================================================ mesh (gmsh) def make_disk_mesh(radius, h): gmsh.initialize() gmsh.option.setNumber("General.Terminal", 0) gmsh.model.add("cymbal") gmsh.model.occ.addDisk(0, 0, 0, radius, radius) gmsh.model.occ.synchronize() gmsh.option.setNumber("Mesh.MeshSizeMin", 0.8 * h) gmsh.option.setNumber("Mesh.MeshSizeMax", h) gmsh.option.setNumber("Mesh.Algorithm", 6) # Frontal-Delaunay gmsh.option.setNumber("Mesh.Smoothing", 10) gmsh.model.mesh.generate(2) tags, xyz, _ = gmsh.model.mesh.getNodes() _, tri_tags = gmsh.model.mesh.getElementsByType(2) # 3-node triangles gmsh.finalize() tags = np.asarray(tags, dtype=np.int64) xyz = np.asarray(xyz).reshape(-1, 3) remap = np.full(tags.max() + 1, -1, dtype=np.int64) remap[tags] = np.arange(len(tags)) tri = remap[np.asarray(tri_tags, dtype=np.int64).reshape(-1, 3)] used = np.unique(tri) xy = xyz[used, :2] inv = np.full(len(tags), -1, dtype=np.int64); inv[used] = np.arange(len(used)) tri = inv[tri] a = ((xy[tri[:, 1], 0] - xy[tri[:, 0], 0]) * (xy[tri[:, 2], 1] - xy[tri[:, 0], 1]) - (xy[tri[:, 2], 0] - xy[tri[:, 0], 0]) * (xy[tri[:, 1], 1] - xy[tri[:, 0], 1])) tri[a < 0] = tri[a < 0][:, [0, 2, 1]] # counter-clockwise return xy, tri # ============================================================================ DKT element def dkt_stiffness(xe, ye, D, nu): """9x9 DKT bending stiffness, DOFs (w1,θx1,θy1, w2,...), θx = w,y, θy = −w,x.""" i, j = jnp.array([1, 2, 0]), jnp.array([2, 0, 1]) # sides k=4,5,6 ↔ ij = 23, 31, 12 xij, yij = xe[i] - xe[j], ye[i] - ye[j] l2 = xij ** 2 + yij ** 2 a = -xij / l2 b = 0.75 * xij * yij / l2 c = (0.25 * xij ** 2 - 0.5 * yij ** 2) / l2 d = -yij / l2 e = (0.25 * yij ** 2 - 0.5 * xij ** 2) / l2 a4, a5, a6 = a; b4, b5, b6 = b; c4, c5, c6 = c; d4, d5, d6 = d; e4, e5, e6 = e def H(s): xi, eta = s lam = 1 - xi - eta N1, N2, N3 = lam * (2 * lam - 1), xi * (2 * xi - 1), eta * (2 * eta - 1) N4, N5, N6 = 4 * xi * eta, 4 * eta * lam, 4 * xi * lam Hx = jnp.stack([1.5 * (a6 * N6 - a5 * N5), b5 * N5 + b6 * N6, N1 - c5 * N5 - c6 * N6, 1.5 * (a4 * N4 - a6 * N6), b6 * N6 + b4 * N4, N2 - c4 * N4 - c6 * N6, 1.5 * (a5 * N5 - a4 * N4), b4 * N4 + b5 * N5, N3 - c4 * N4 - c5 * N5]) Hy = jnp.stack([1.5 * (d6 * N6 - d5 * N5), -N1 + e5 * N5 + e6 * N6, -b5 * N5 - b6 * N6, 1.5 * (d4 * N4 - d6 * N6), -N2 + e4 * N4 + e6 * N6, -b6 * N6 - b4 * N4, 1.5 * (d5 * N5 - d4 * N4), -N3 + e4 * N4 + e5 * N5, -b4 * N4 - b5 * N5]) return jnp.stack([Hx, Hy]) # (2, 9) x21, x31 = xe[1] - xe[0], xe[2] - xe[0] y21, y31 = ye[1] - ye[0], ye[2] - ye[0] det = x21 * y31 - x31 * y21 # = 2A xi_x, xi_y, eta_x, eta_y = y31 / det, -x31 / det, -y21 / det, x21 / det def B_at(s): J = jax.jacfwd(H)(s) # (2, 9, 2) Hx_xi, Hx_eta, Hy_xi, Hy_eta = J[0, :, 0], J[0, :, 1], J[1, :, 0], J[1, :, 1] Hx_x = Hx_xi * xi_x + Hx_eta * eta_x; Hx_y = Hx_xi * xi_y + Hx_eta * eta_y Hy_x = Hy_xi * xi_x + Hy_eta * eta_x; Hy_y = Hy_xi * xi_y + Hy_eta * eta_y return jnp.stack([Hx_x, Hy_y, Hx_y + Hy_x]) # curvature operator (3, 9) Db = D * jnp.array([[1, nu, 0], [nu, 1, 0], [0, 0, 0.5 * (1 - nu)]]) gps = jnp.array([[1 / 6, 1 / 6], [2 / 3, 1 / 6], [1 / 6, 2 / 3]]) # exact for quadratic B Bs = jax.vmap(B_at)(gps) # (3, 3, 9) return (0.5 * det / 3) * jnp.einsum("gki,kl,glj->ij", Bs, Db, Bs) # ============================================================================ element energy def elem_energy_parts(qe, Gx, Gy, A, z0x, z0y, Cm, Ke): """qe: (3,5) nodal [u,v,w,θx,θy]. Returns (membrane, bending) energy of one element. The z0 terms are the *linear* curvature coupling; NONLINEAR gates only the quadratic terms.""" nl = 1.0 if NONLINEAR else 0.0 ux, uy = Gx @ qe[:, 0], Gy @ qe[:, 0] vx, vy = Gx @ qe[:, 1], Gy @ qe[:, 1] wx, wy = Gx @ qe[:, 2], Gy @ qe[:, 2] exx = ux + nl * 0.5 * wx * wx + z0x * wx eyy = vy + nl * 0.5 * wy * wy + z0y * wy gxy = uy + vx + nl * wx * wy + z0x * wy + z0y * wx C11, C12, C33 = Cm Em = 0.5 * A * (C11 * (exx * exx + eyy * eyy) + 2 * C12 * exx * eyy + C33 * gxy * gxy) U = qe[:, 2:].reshape(9) return Em, 0.5 * (U @ Ke @ U) def elem_energy(*args): Em, Eb = elem_energy_parts(*args) return Em + Eb # ============================================================================ diagnostics def spectral_onset(p, fs, label, bands=(60, 125, 250, 500, 1000, 2000, 4000, 8000, 16000), t_max=0.5): """Per-band energy envelopes at ~6 ms resolution; onset = first time a band reaches −20 dB of its own peak. Also plots the envelopes (log-log), which shows the cascade far better than a spectrogram normalised to the low-frequency peak.""" from scipy.signal import stft f_, t_, Z = stft(p, fs, nperseg=256, noverlap=224) Pw = np.abs(Z) ** 2 tot_max = Pw.sum(0).max() envs = [] print(f"\n[{label}] band onset(-20 dB) peak level vs loudest band") for lo, hi in zip(bands[:-1], bands[1:]): m = (f_ >= lo) & (f_ < hi) env = Pw[m].sum(0) if env.max() <= 0: continue env_db = 10 * np.log10(env / env.max() + 1e-30) onset = t_[np.argmax(env_db > -20)] print(f" {lo:5d}–{hi:5d} Hz {onset * 1e3:7.1f} ms {t_[np.argmax(env)] * 1e3:7.1f} ms" f" {10 * np.log10(env.max() / tot_max):6.1f} dB") envs.append((lo, hi, env)) cen = (Pw * f_[:, None]).sum(0) / (Pw.sum(0) + 1e-30) for tq in (0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5): if tq < t_[-1]: print(f" spectral centroid at {tq * 1e3:5.0f} ms: {cen[np.searchsorted(t_, tq)]:7.0f} Hz") try: import matplotlib; matplotlib.use("Agg"); import matplotlib.pyplot as plt fig, ax = plt.subplots(1, 2, figsize=(13, 4.5)) sel = t_ <= t_max ax[0].pcolormesh(t_[sel] * 1e3, f_, 10 * np.log10(Pw[:, sel] / Pw.max() + 1e-12), vmin=-90, vmax=0, shading="auto") ax[0].set_yscale("log"); ax[0].set_ylim(50, fs / 2) ax[0].set_xlabel("t [ms]"); ax[0].set_ylabel("f [Hz]"); ax[0].set_title(label) for lo, hi, env in envs: ax[1].plot(t_[sel] * 1e3, 10 * np.log10(env[sel] / tot_max + 1e-30), label=f"{lo}–{hi} Hz") ax[1].set_ylim(-90, 5); ax[1].set_xscale("log"); ax[1].set_xlim(2, t_max * 1e3) ax[1].set_xlabel("t [ms]"); ax[1].set_ylabel("band energy [dB re loudest]") ax[1].legend(fontsize=7, ncol=2); ax[1].grid(alpha=.3) ax[1].set_title("band envelopes — a cascade = staggered rises, high bands last") plt.tight_layout(); plt.savefig(f"{label}_cascade.png", dpi=120); plt.close() print(f" wrote {label}_cascade.png") except Exception: pass # ============================================================================ main def main(): t_start = time.time() tag = OUT_PREFIX if NONLINEAR else OUT_PREFIX + "_linear" xy, tri = make_disk_mesh(R, H_ELEM) N, E = len(xy), len(tri) print(f"mesh: {N} nodes, {E} triangles, {5 * N} DOFs (NONLINEAR = {NONLINEAR})") # ---- element geometry on the projected domain (numpy, f64) x, y = xy[tri, 0], xy[tri, 1] # (E, 3) x21, x31 = x[:, 1] - x[:, 0], x[:, 2] - x[:, 0] y21, y31 = y[:, 1] - y[:, 0], y[:, 2] - y[:, 0] area = 0.5 * (x21 * y31 - x31 * y21) assert (area > 0).all() Gx = np.stack([y[:, 1] - y[:, 2], y[:, 2] - y[:, 0], y[:, 0] - y[:, 1]], 1) / (2 * area[:, None]) Gy = np.stack([x[:, 2] - x[:, 1], x[:, 0] - x[:, 2], x[:, 1] - x[:, 0]], 1) / (2 * area[:, None]) xc, yc = x.mean(1), y.mean(1) def locate(px, py): """Containing element and barycentric weights of a point on the projected domain.""" l1 = ((y[:, 1] - y[:, 2]) * (px - x[:, 2]) + (x[:, 2] - x[:, 1]) * (py - y[:, 2])) / (2 * area) l2 = ((y[:, 2] - y[:, 0]) * (px - x[:, 2]) + (x[:, 0] - x[:, 2]) * (py - y[:, 2])) / (2 * area) l3 = 1 - l1 - l2 inside = (l1 >= -1e-9) & (l2 >= -1e-9) & (l3 >= -1e-9) assert inside.any(), f"point ({px:.3f}, {py:.3f}) outside mesh" es = int(np.argmax(inside)) return tri[es], np.array([l1[es], l2[es], l3[es]]) # ---- profile: thickness, initial shape and its analytic gradient h_e = np.asarray(_h_eval(jnp.asarray(xc), jnp.asarray(yc))) h_node = np.asarray(_h_eval(jnp.asarray(xy[:, 0]), jnp.asarray(xy[:, 1]))) z0_node = np.asarray(_z0_eval(jnp.asarray(xy[:, 0]), jnp.asarray(xy[:, 1]))) grad_z0 = np.asarray(_z0_grad(jnp.stack([jnp.asarray(xc), jnp.asarray(yc)], 1))) z0x, z0y = grad_z0[:, 0], grad_z0[:, 1] slope = np.hypot(z0x, z0y).max() rim_node = np.hypot(xy[:, 0], xy[:, 1]) > R - 1e-6 print(f"profile: height {z0_node.max() * 1e3:.1f} mm, rim warp p-p " f"{(z0_node[rim_node].max() - z0_node[rim_node].min()) * 1e3:.2f} mm, " f"{N_DENTS} dents, h ∈ [{h_e.min() * 1e3:.2f}, {h_e.max() * 1e3:.2f}] mm, " f"max |∇z0| = {slope:.3f}") if slope > MAX_SLOPE: print(f"WARNING: max slope {slope:.2f} exceeds shallow-shell bound {MAX_SLOPE}") D_e = E_MOD * h_e ** 3 / (12 * (1 - NU ** 2)) Cm_e = (E_MOD * h_e / (1 - NU ** 2))[:, None] * np.array([1.0, NU, 0.5 * (1 - NU)]) # ---- DKT stiffness for every element Ke64 = jax.vmap(dkt_stiffness, in_axes=(0, 0, 0, None))(jnp.asarray(x), jnp.asarray(y), jnp.asarray(D_e), NU) Ke_np = np.asarray(Ke64) # element self-tests: rigid tilt → zero force; w = ½x² → energy ½·D·A U_tilt = np.stack([y, np.ones_like(y), np.zeros_like(y)], 2).reshape(E, 9) res = np.abs(np.einsum("eij,ej->ei", Ke_np, U_tilt)).max() / np.abs(Ke_np).max() U_cur = np.stack([0.5 * x ** 2, np.zeros_like(x), -x], 2).reshape(E, 9) Eb = 0.5 * np.einsum("ei,eij,ej->e", U_cur, Ke_np, U_cur) err = np.abs(Eb / (0.5 * D_e * area) - 1).max() print(f"DKT checks: rigid-tilt residual {res:.1e}, constant-curvature energy error {err:.1e}") # ---- lumped mass (per DOF): [u, v, w, θx, θy] m_e = RHO * h_e * area / 3 I_e = m_e * (h_e ** 2 / 12 + ROT_INERTIA_FACTOR * area) M = np.zeros((N, 5)) for k in range(3): np.add.at(M[:, 0], tri[:, k], m_e * INPLANE_MASS_SCALE) np.add.at(M[:, 1], tri[:, k], m_e * INPLANE_MASS_SCALE) np.add.at(M[:, 2], tri[:, k], m_e) np.add.at(M[:, 3], tri[:, k], I_e) np.add.at(M[:, 4], tri[:, k], I_e) print(f"plate mass {M[:, 2].sum():.3f} kg") # ---- stable dt from the element eigenvalue bound (Irons): λ_max(M⁻¹K) ≤ max_e λ_max(M_e⁻¹K_e) Khess = jax.vmap(jax.hessian(elem_energy))( jnp.zeros((E, 3, 5)), jnp.asarray(Gx), jnp.asarray(Gy), jnp.asarray(area), jnp.asarray(z0x), jnp.asarray(z0y), jnp.asarray(Cm_e), Ke64).reshape(E, 15, 15) Me = np.tile(np.stack([m_e * INPLANE_MASS_SCALE, m_e * INPLANE_MASS_SCALE, m_e, I_e, I_e], 1), (1, 3)) # (E, 15) node-major S = np.asarray(Khess) / np.sqrt(Me[:, :, None] * Me[:, None, :]) lam = np.linalg.eigvalsh(S) n_zero = int((lam[0] < 1e-9 * lam[0].max()).sum()) print(f"element 0 zero-energy modes: {n_zero} (expect 6)") e_worst = int(np.argmax(lam[:, -1])) omega_max = math.sqrt(lam[:, -1].max()) dt_bound = 2.0 / omega_max M_SUB = math.ceil(1.0 / (FS * CFL_SAFETY * dt_bound)) dt = 1.0 / (FS * M_SUB) n_blocks = math.ceil(DURATION * FS * M_SUB / STEPS_PER_BLOCK) n_steps = n_blocks * STEPS_PER_BLOCK print(f"ω_max bound {omega_max:.3e} rad/s (limiting element at r = " f"{math.hypot(xc[e_worst], yc[e_worst]):.3f} m, h = {h_e[e_worst] * 1e3:.2f} mm)") print(f"dt = {dt:.3e} s, {M_SUB} sub-steps/sample, {n_steps} steps") # ---- rigid-body modes of the curved shell, M-orthonormalised # A rigid tilt drags the in-plane displacements along: rotation about x by φ gives # v = −φ z0, w = φ y; about y gives u = +φ z0, w = −φ x (first order, shallow-shell). X, Y = xy[:, 0] - xy[:, 0].mean(), xy[:, 1] - xy[:, 1].mean() Rm = np.zeros((6, N, 5)) Rm[0, :, 0] = 1; Rm[1, :, 1] = 1 # in-plane translations Rm[2, :, 0] = -Y; Rm[2, :, 1] = X # in-plane rotation Rm[3, :, 2] = 1 # z translation Rm[4, :, 2] = Y; Rm[4, :, 3] = 1; Rm[4, :, 1] = -z0_node # tilt about x (θx = w,y) Rm[5, :, 2] = -X; Rm[5, :, 4] = 1; Rm[5, :, 0] = z0_node # tilt about y (θy = −w,x) for k in range(6): for l in range(k): Rm[k] -= (Rm[l] * M * Rm[k]).sum() * Rm[l] Rm[k] /= math.sqrt((Rm[k] * M * Rm[k]).sum()) # ---- contact patch: a set of power-law springs F = k δ^α on a rigid exciter of mass M_EX xs, ys = STRIKE_R * math.cos(STRIKE_PHI), STRIKE_R * math.sin(STRIKE_PHI) if EXCITER == "stick": M_EX, V_EX, ALPHA, K_TOTAL = M_STICK, V_STRIKE, 1.5, K_HERTZ ang = STRIKE_PHI + PATCH_ANGLE xi = np.linspace(-0.5, 0.5, PATCH_POINTS) * PATCH_LEN px, py = xs + xi * math.cos(ang), ys + xi * math.sin(ang) elif EXCITER == "mallet": M_EX, V_EX, ALPHA = M_MALLET, V_MALLET, ALPHA_MALLET # k such that the head stops after MALLET_COMPRESSION on a rigid wall K_TOTAL = (ALPHA + 1) * 0.5 * M_EX * V_EX ** 2 / MALLET_COMPRESSION ** (ALPHA + 1) th = np.arange(6) * np.pi / 3 px = np.concatenate([[xs], xs + MALLET_RADIUS * np.cos(th)]) py = np.concatenate([[ys], ys + MALLET_RADIUS * np.sin(th)]) else: raise ValueError(f"unknown EXCITER {EXCITER!r}") P = len(px) loc = [locate(px[k], py[k]) for k in range(P)] s_idx = np.array([l[0] for l in loc]) # (P, 3) s_wts = np.array([l[1] for l in loc]) # (P, 3) z0_p = np.asarray(_z0_eval(jnp.asarray(px), jnp.asarray(py))) # rigid exciter face: best-fit plane through the surface heights; residual → staggered contact Apl = np.column_stack([px - xs, py - ys, np.ones(P)]) s_off = z0_p - Apl @ np.linalg.lstsq(Apl, z0_p, rcond=None)[0] zs0 = s_off.max() + STICK_GAP tau_est = 2.0 * MALLET_COMPRESSION / V_EX if EXCITER == "mallet" else None print(f"{EXCITER}: {P}-point patch, k = {K_TOTAL:.3g} N/m^{ALPHA}, m = {M_EX} kg, v = {V_EX} m/s, " f"energy {0.5 * M_EX * V_EX ** 2:.2f} J" + (f", free contact ≈ {tau_est * 1e3:.1f} ms → direct excitation below ~{0.5 / tau_est:.0f} Hz" if tau_est else "")) # ---- radiation geometry (element centroids on the curved surface, retarded times in sub-steps) cent = np.stack([xc, yc, z0_node[tri].mean(1)], 1) # (E, 3) r_front = np.linalg.norm(LISTENER - cent, axis=1) phi = np.linspace(0, 2 * np.pi, 720, endpoint=False) rim_xy = np.stack([R * np.cos(phi), R * np.sin(phi)], 1) rim_z = np.asarray(_z0_eval(jnp.asarray(rim_xy[:, 0]), jnp.asarray(rim_xy[:, 1]))) rim = np.column_stack([rim_xy, rim_z]) d1 = np.linalg.norm(cent[:, None, :] - rim[None, :, :], axis=2) d2 = np.linalg.norm(rim - LISTENER, axis=1) r_back = (d1 + d2[None, :]).min(1) coef_f = RHO_AIR * area / (2 * np.pi * r_front) coef_b = BACK_DIFFRACTION_GAIN * RHO_AIR * area / (2 * np.pi * r_back) df, db = r_front / (C_AIR * dt), r_back / (C_AIR * dt) i0f, frf = np.floor(df).astype(np.int32), df - np.floor(df) i0b, frb = np.floor(db).astype(np.int32), db - np.floor(db) buf_len = n_steps + int(max(i0f.max(), i0b.max())) + 2 # ---- device arrays in simulation precision f = lambda a: jnp.asarray(a, dtype=DTYPE) tri_j = jnp.asarray(tri, dtype=jnp.int32) Gx_j, Gy_j, area_j = f(Gx), f(Gy), f(area) z0x_j, z0y_j, Cm_j, Ke_j = f(z0x), f(z0y), f(Cm_e), f(Ke_np) rhoh_area_j, hnode_j = f(RHO * h_e * area), f(h_node) M_j, Minv_j, Rm_j = f(M), f(1.0 / M), f(Rm) sidx_j, swts_j, soff_j = jnp.asarray(s_idx, dtype=jnp.int32), f(s_wts), f(s_off) coef_f_j, coef_b_j, frf_j, frb_j = f(coef_f), f(coef_b), f(frf), f(frb) i0f_j, i0b_j = jnp.asarray(i0f), jnp.asarray(i0b) k_point = K_TOTAL / P def total_energy(q): return jax.vmap(elem_energy)(q[tri_j], Gx_j, Gy_j, area_j, z0x_j, z0y_j, Cm_j, Ke_j).sum() f_int = jax.grad(lambda q: -total_energy(q)) # internal nodal forces # shell check: tangent stiffness at rest must annihilate all six rigid modes q0 = jnp.zeros((N, 5), DTYPE) K_apply = lambda d: jax.jvp(f_int, (q0,), (d,))[1] chk = np.random.default_rng(0) worst = 0.0 for k in range(6): rnd = chk.standard_normal((N, 5)); rnd *= np.linalg.norm(Rm[k]) / np.linalg.norm(rnd) worst = max(worst, float(jnp.abs(K_apply(Rm_j[k])).max() / jnp.abs(K_apply(f(rnd))).max())) print(f"rigid modes: max |K·r| / |K·random| = {worst:.1e} (should be ~1e-6 or smaller)") def damp_lap(vw): """−2σ1 ρh ∇²ẇ as a nodal force: 2σ1 · (ρh A_e) · Σ_e G_iᵀ (G ẇ_e).""" ve = vw[tri_j] vx, vy = (Gx_j * ve).sum(1), (Gy_j * ve).sum(1) c = rhoh_area_j[:, None] * (Gx_j * vx[:, None] + Gy_j * vy[:, None]) return 2 * SIGMA1 * jax.ops.segment_sum(c.ravel(), tri_j.ravel(), num_segments=N) def project_rigid(v): coef = jnp.einsum("knd,nd->k", Rm_j, M_j * v) return v - jnp.einsum("k,knd->nd", coef, Rm_j) def step(n, st): q, v, zs, vs, bf, bb, Fmax, ncont = st # exciter / plate contact along the patch wp = jnp.einsum("pk,pk->p", swts_j, q[sidx_j, 2]) + soff_j # surface vs exciter face delta = jnp.maximum(wp - zs, 0.0) # (P,) Fk = k_point * jnp.power(delta, ALPHA) F = Fk.sum() # forces fo = f_int(q) fo = fo.at[sidx_j, 2].add(-Fk[:, None] * swts_j) fo = fo.at[:, 2].add(-damp_lap(v[:, 2])) # leapfrog (σ0 implicit), rigid-mode projection v1 = (v + dt * fo * Minv_j) / (1 + 2 * SIGMA0 * dt) v1 = project_rigid(v1) q1 = q + dt * v1 vs1 = vs + dt * F / M_EX zs1 = zs + dt * vs1 # radiation: scatter element accelerations into retarded-time buffers acc = ((v1[:, 2] - v[:, 2]) / dt)[tri_j].mean(1) sf, sb = coef_f_j * acc, coef_b_j * acc bf = bf.at[n + i0f_j].add(sf * (1 - frf_j)).at[n + i0f_j + 1].add(sf * frf_j) bb = bb.at[n + i0b_j].add(sb * (1 - frb_j)).at[n + i0b_j + 1].add(sb * frb_j) return (q1, v1, zs1, vs1, bf, bb, jnp.maximum(Fmax, F), ncont + (F > 0)) @jax.jit def run_block(n0, st): return lax.fori_loop(0, STEPS_PER_BLOCK, lambda k, s: step((n0 + k).astype(jnp.int32), s), st) @jax.jit def energies(st): q, v = st[0], st[1] Em, Eb = jax.vmap(elem_energy_parts)(q[tri_j], Gx_j, Gy_j, area_j, z0x_j, z0y_j, Cm_j, Ke_j) return (0.5 * (M_j * v * v).sum(), Em.sum(), Eb.sum(), 0.5 * M_EX * st[3] ** 2, jnp.abs(q[:, 2] / hnode_j).max()) state = (jnp.zeros((N, 5), DTYPE), jnp.zeros((N, 5), DTYPE), jnp.asarray(zs0, DTYPE), jnp.asarray(-V_EX, DTYPE), jnp.zeros(buf_len, DTYPE), jnp.zeros(buf_len, DTYPE), jnp.asarray(0.0, DTYPE), jnp.asarray(0, jnp.int32)) print(f"setup done in {time.time() - t_start:.1f} s; compiling + running on " f"{jax.devices()[0].platform} ...") print(" (nonlinearity meter: max|w|/h ≲ 1 → effectively linear; 3–5 → clear cascade)") t_run = time.time() for blk in range(n_blocks): state = run_block(jnp.asarray(blk * STEPS_PER_BLOCK, jnp.int32), state) if blk % max(1, n_blocks // 20) == 0 or blk == n_blocks - 1: ke, em, eb, ks, woh = energies(state) jax.block_until_ready(ke) print(f" t = {(blk + 1) * STEPS_PER_BLOCK * dt:6.3f} s KE {float(ke):.3e} " f"membrane {float(em):.3e} bending {float(eb):.3e} J exciter {float(ks):.3e} J " f"max|w|/h = {float(woh):5.2f} [{time.time() - t_run:.0f} s wall]") q, v, zs, vs, bf, bb, Fmax, ncont = jax.device_get(state) print(f"{EXCITER}: peak force {float(Fmax):.1f} N, contact {float(ncont) * dt * 1e3:.2f} ms, " f"rebound velocity {float(vs):+.2f} m/s " f"(energy left in exciter {0.5 * M_EX * float(vs) ** 2:.2f} of {0.5 * M_EX * V_EX ** 2:.2f} J)") # ---- decimate to audio rate, write, analyse n_out = int(FS * DURATION) p_ray = resample_poly(np.asarray(bf[:n_steps], np.float64), 1, M_SUB)[:n_out] p_back = resample_poly(np.asarray(bb[:n_steps], np.float64), 1, M_SUB)[:n_out] p_free = p_ray - p_back for name, p in (("rayleigh", p_ray), ("unbaffled", p_free)): peak = np.abs(p).max() print(f"{name:10s}: peak {peak:.3f} Pa ({20 * np.log10(peak / 2e-5):.1f} dB SPL peak)") wavfile.write(f"{tag}_{name}.wav", FS, (0.9 * p / peak * 32767).astype(np.int16)) print(f"wrote {tag}_rayleigh.wav, {tag}_unbaffled.wav (total {time.time() - t_start:.0f} s)") spectral_onset(p_free, FS, f"{tag}_unbaffled") if __name__ == "__main__": main()