Modal2D — Technical Overview

Modal2D is the standard discretize-then-diagonalize modal-sound pipeline (van den Doel/Pai; O'Brien et al.) factored into a clean engine/application split. The boundary is drawn at the level of opaque modal results: libmodal forward-declares handle types ModalObject and ModalResult and exposes only free functions over pointers to them, so the viewer never sees a stiffness matrix, an eigenvector layout, or a device pointer through the public header — only modal_create_*, modal_solve, modal_strike, modal_get_displacements, modal_fill_audio, and read-only query functions.

The engine computes eigenpairs $(\omega_i,\phi_i)$ once per geometry, augments them at solve time with per-mode damping ratios and radiation coefficients, then evaluates a bank of damped oscillators

A strike projects an impulse onto the mass-orthonormal modal basis via the boundary-normal-weighted mode value at the strike point; the per-mode, per-slot excitation amplitudes feed a CUDA synthesis kernel that produces a 44.1 kHz PCM stream, post-processed by the viewer's PortAudio callback. The same modal state drives the animated mode-shape visualization and the spectral-analysis panels, so synthesis sits at the junction of audio output and visual feedback. The GPU is used in two unrelated ways: CUDA inside libmodal (a cuSOLVER dense generalized eigensolver and a per-sample audio-synthesis kernel) and OpenGL inside the viewer (mesh shaders plus a GLSL compute shader that paints a heuristic radiated-pressure field). On top of the forward pipeline sits an inverse-design outer loop: each tuning iteration remeshes, re-solves, reads back the modal spectrum, and proposes new geometric parameters until a struck partial matches a chromatic target.

The engine/app boundary as a contract

The boundary is drawn at the level of opaque modal results. libmodal forward-declares handle types ModalObject and ModalResult (modal/types.h:144) plus free functions over pointers to them; the concrete struct ModalObject / struct ModalResult are private (modal_api.cpp:34 and :152). The viewer never sees a stiffness matrix, an eigenvector layout, or a device pointer through the public header — only modal_create_*, modal_solve, modal_strike, modal_get_displacements, modal_fill_audio, and read-only query functions. This is the standard modal sound-synthesis pipeline (van den Doel/Pai; O'Brien et al.) factored into a library:

The engine computes eigenpairs $(\omega_i,\phi_i)$ once per geometry (modal_solve), then evaluates a bank of damped oscillators. A strike projects an impulse onto the modal basis via the boundary-normal-weighted mode value at the strike point: compute_strike_impulse builds $f_{\text{impulse}}$, then project_strike (oscillator.cu:15) sets the per-mode, per-slot excitation $A_{i,k}=f_{\text{impulse},i}\cdot\text{magnitude}$ (clamped to $\pm 50$, with a $|A|<0.005$ dead-zone).

Per-mode damping and radiation weighting (implemented, beyond the textbook model)

Two details make the synthesis richer than a bare modal sum. Damping is per-mode and frequency-dependent:

(modal_api.cpp:162, materialized into mode_zetas at :2838-2846 / :3299-3308 and uploaded to d_zetas), so high overtones decay faster. Each mode also carries a scalar radiation coefficient $\alpha_i$ from a discrete monopole boundary integral (radiation.h, compute_radiation_coefficients):

where $\phi_{i,n,j}$ is the outward-normal displacement of mode $i$ at boundary node $j$ and $A_j$ the node's boundary-segment half-length. $\alpha_i$ scales each mode's contribution to the output sample, so poorly-radiating modes are attenuated in the audio.

Where the GPU enters

Two numerically heavy steps map to data-parallel hardware:

1. The generalized symmetric eigenproblem $K\phi=\lambda M\phi$. For small systems the engine forms dense single-precision $K,M$ (column-major, lower-triangular) and calls cuSOLVER's $\texttt{Ssygvdx}$ with CUSOLVER_EIG_RANGE_I, index range $[1,\min(n,m{+}3)]$, returning the lowest eigenpairs of a symmetric-definite pencil with the eigenvectors overwriting $d\_K$. Frequencies are $\omega_i=\sqrt{\lambda_i}$, with a rigid-body filter dropping $\omega_i<\texttt{rigid\_threshold}$.
2. Audio synthesis, embarrassingly parallel over output samples. One CUDA thread per PCM sample $s$ computes $t=t_{\text{start}}+s\,\Delta t$ and the radiation-weighted modal sum, with a global pitch multiplier freq_shift:

where $K=\texttt{STRIKE\_RING\_SIZE}$ is the per-mode strike ring (audio_sum_kernel, synthesis.cu:24-66; $p$ is finite-guarded). For $M$ modes and a 512-sample buffer this is $\approx 512\cdot M\cdot K$ transcendental evaluations per callback in one launch.

The viewer's GPU work is graphics, not solve. The deformed mesh is drawn by passing per-node 2D displacement $\mathbf{u}_i(t)$ as a vertex attribute; the radiated field is a per-pixel heuristic Kirchhoff-style sum over (decimated) deforming boundary sources:

with $v_{n,j}$ the boundary node's outward-normal velocity, $\omega_0$ the lowest modal frequency (freqs[0], used as a single dominant tone), $\theta_j$ the source-to-pixel angle vs. the outward normal ($\cos\theta_j<0\Rightarrow$ skip, a half-space cull), source stride $\texttt{src\_step}=\max(1,n_b/60)$, and $c=343\,\text{m/s}$ hardcoded. The $1/(0.5+3\sqrt r)$ falloff is ad-hoc. This is a real-time visual cartoon, not a coupled BEM solve, and does not feed back into audio.

Per-frame data flow

App::run (app.cpp:8736) is glfwPollEvents → process_input → render. Each App::render (app.cpp:1113): (a) advances a slow-motion-scalable clock; (b) polls the background SolveWorker and swaps in a new ModalResult if ready; (c) reads animated displacements (modal_get_displacements, app.cpp:1496) and uploads them to the displacement VBO (cuda_gl_buf_.write_displacements, :1563); (d) runs the mesh GL block (wireframe/filled shaders) before ui_.begin_frame() (:1778, which calls ImGui::NewFrame); (e) enters the ImGui phase, where it draws panels/overlays AND calls draw_medium_propagation() (:4986) — this dispatches the medium compute shader and blits its texture via ImGui's background draw list, so the compute step lives in the ImGui phase, not the pre-NewFrame block; (f) glfwSwapBuffers (:5462). Audio runs on a separate PortAudio callback thread (pa_callback) that calls modal_fill_audio, which drains the strike buffer, projects strikes under a mutex, uploads the updated $(A,t_0)$ ring arrays, and launches the synthesis kernel.

The library/app split and the C-API boundary

• libmodal builds as modal_lib (STATIC): libmodal/CMakeLists.txt GLOBs src/*.cpp src/*.cu, links Eigen3::Eigen, CUDA::cudart, CUDA::cusolver, CUDA::cublas, and compiles vendored Shewchuk Triangle (extern/triangle/triangle.c, defines ANSI_DECLARATORS TRILIBRARY NO_TIMER). Built with CUDA_SEPARABLE_COMPILATION ON, CUDA_RESOLVE_DEVICE_SYMBOLS ON, POSITION_INDEPENDENT_CODE ON.
• Public surface libmodal/include/modal/modal.h: all functions in namespace modal, trafficking in opaque ModalObject*/ModalResult* plus POD structs. Creation: modal_create_from_polygon, modal_create_preset, modal_create_solid, modal_create_string, modal_create_tuning_component, modal_create_custom_instrument (note: _from_polygon, not _polygon). Configuration (modal_set_material/_medium/_curvature/_thickness/_plate_bending/...), solve (modal_solve, :130), excitation (modal_strike :133, modal_strike_uniform :228, modal_strike_solid_at :231, modal_compute_strike_amplitudes :143), readout (modal_get_displacements :153, mode-amplitude queries), audio (modal_fill_audio :208, modal_fill_audio_snapshot :215), plus read-only tuning/viz queries.
• modal_viewer builds the executable: modal_viewer/CMakeLists.txt GLOBs src/*.cpp and links modal_lib glad glfw CUDA::cudart opengl32 portaudio_static imgui_lib comdlg32. Entry: main.cpp → viewer::App app; app.run();.

GPU / CUDA usage (engine side)

• Dense eigensolve on GPU — eigensolve_dense.cu:solve_dense. Converts Eigen sparse double K,M to column-major dense float (lower triangle), then cusolverDnSsygvdx_bufferSize + cusolverDnSsygvdx with CUSOLVER_EIG_TYPE_1, CUSOLVER_EIG_MODE_VECTOR, CUSOLVER_EIG_RANGE_I, CUBLAS_FILL_MODE_LOWER, $iu=\min(n,\,\texttt{num\_modes}+3)$. Eigenvectors come back in overwritten d_K; nonzero info throws.
• Dispatch — eigensolve.cpp:solve_eigenmodes: K.rows() <= DENSE_DOF_THRESHOLD (512) goes to the GPU dense path inside a try/catch; on exception (documented as float32 Cholesky failing on FINE meshes) or for larger meshes it falls back to solve_sparse. A per-component block solve solve_eigenmodes_per_component extracts each connected component's K/M sub-blocks, solves them independently, scatters eigenvectors to full length, and merges by ascending frequency.
• CPU sparse fallback — eigensolve_sparse.cpp:solve_sparse: a hand-rolled shift-invert Lanczos for the generalized pencil (no ARPACK/Spectra). It regularizes zero-stiffness/zero-mass DOFs, factors $K-\sigma M$ with Eigen::SimplicialLDLT (from <Eigen/SparseCholesky>) at $\sigma=10^{-4}$, runs the M-orthonormal Lanczos recurrence with full reorthogonalization from a fixed-seed (mt19937(42)) start vector, then maps Ritz values via $\lambda=1/\mu+\sigma$.
• CUDA audio synthesis — synthesis.cu. SynthesisState (audio/synthesis.h) holds device arrays d_freqs/d_zetas/d_amps/d_onset/d_alphas/d_pcm, pinned h_pcm, and a dedicated cudaStream_t. audio_sum_kernel is one thread per sample (block=256) summing $\alpha_i\,A\,e^{-\zeta\omega'\,\delta}\sin(\omega_d\,\delta)$ over modes × STRIKE_RING_SIZE with freq_shift and finite-guarding. synth_fill launches on the stream, cudaMemcpyAsyncs to pinned memory, syncs, copies out (capped at 512 samples). A CPU mirror (oscillator.cu:eval_modes_cpu/eval_single_mode) serves the visual path and tests.
• GPU strike upload — modal_api.cpp:modal_fill_audio (:3839) drains strikes, calls project_strike under mode_states_mutex, then cudaMemcpyAsyncs the rebuilt $(A,t_0)$ rings into d_amps/d_onset before synth_fill. Radiation coefficients are computed at solve time (:3225 via compute_radiation_coefficients) and uploaded to d_alphas (:3324-3328).

GPU usage (viewer side)

• Mesh rendering — inline GLSL 450-core shaders in App::init_gl (app.cpp:192): a wireframe shader (reads vec2 aPos,vec2 aDisp, applies zoom/pan/aspect) and a filled shader with heat/Chladni colormaps, plus separate 3D shaders. Geometry uploaded once after solve; the per-node displacement VBO updated each frame.
• "CUDA-GL interop" is a plain CPU upload — renderer/cuda_gl_interop.cpp:CudaGLBuffer::write_displacements is just glBufferSubData of the host displacement array (2 floats/node) into the VBO. No cudaGraphicsGLRegisterBuffer/MapResources exists in shipped source; the header comment says it merely mirrors the interface for a future zero-copy upgrade. Displacements thus round-trip device→host (in libmodal) then host→VBO.
• Medium compute shader — source compiled at app.cpp:483 (GL_COMPUTE_SHADER, local_size 16x16); driven by App::draw_medium_propagation (:7183), called at :4986 inside the ImGui phase. It builds two std430 SSBOs on the CPU (deformed boundary sources {x,y,vn,out_nx,out_ny,pad...} and deformed triangles for an in-mesh test), uploads them, sets uniforms (incl. u_speed_of_sound=343, u_dom_omega=freqs[0], u_src_step), dispatches ((gw+15)/16,(gh+15)/16,1), issues GL_SHADER_IMAGE_ACCESS_BARRIER_BIT, then draws the RGBA8 result (blue=positive, red=negative pressure) as a full-viewport image via ImGui::GetBackgroundDrawList().

Data flow plumbing

• Background solving — app/solve_worker.{h,cpp}: start spawns a std::thread running modal_solve; App::render polls is_done() (acquire) and calls take_result(), with documented release/acquire ordering on done_.
• Audio device — audio/audio_engine.cpp: Pa_OpenDefaultStream (stereo, paFloat32, 44100 Hz, 512 frames). pa_callback calls modal_fill_audio into a mono scratch buffer (clamping n>512 to 512), then runs a viewer-side DSP chain (sustain resonator bank, nonlinear voice, gating, cavity/reverb proxies) before writing interleaved stereo. Lock-free double-buffering with atomic indices/sequence counters publishes parameter/strike updates to the audio thread.
• ImGui phase ordering — the mesh GL block runs before ui_.begin_frame() (ui_manager.cpp:262-264 issues ImGui_ImplOpenGL3_NewFrame/ImGui_ImplGlfw_NewFrame/ImGui::NewFrame); overlays and the medium-compute blit run in the ImGui phase (matches the project memory note that ImGui draw-list calls must not happen in the pre-NewFrame GL block).

Build system

• CMakeLists.txt (root): project(ModalSynth LANGUAGES C CXX CUDA), C++17/CUDA17, CMAKE_CUDA_ARCHITECTURES native default, find_package(CUDAToolkit REQUIRED). Eigen 3.4.0, Catch2 v3.5.2, GLFW 3.4, glad v0.1.36 (gl=4.5 core, C generator), PortAudio v19.7.0 (static), Dear ImGui (docking branch) are all FetchContent (shallow git clones; the docking branch is unpinned). imgui_lib = imgui core + GLFW/OpenGL3 backends with IMGUI_IMPL_OPENGL_LOADER_GLAD.
• Two build trees. run_cmake_configure.bat configures with -G Ninja into build/ (-DCMAKE_MAKE_PROGRAM=$NINJA, -DCMAKE_CUDA_COMPILER=$NVCC, -DCMAKE_BUILD_TYPE=Release); run_cmake_build.bat runs cmake --build build --target modal_lib (Ninja, library only). Separately, build_release.bat/build_debug.bat operate on build_vs/ (a Visual Studio generator tree) via MSBuild modal_viewer\modal_viewer.vcxproj /p:Configuration=Release|Debug /p:Platform=x64. Both build_*.bat guard on build_vs\...\modal_viewer.vcxproj and instruct the user to run run_cmake_configure_vs.bat first — but that script is absent from the repo (only run_cmake_configure.bat, run_cmake_build.bat, run_cmake_rebuild.bat exist), so the documented MSBuild flow is not reproducible from the checked-in scripts alone.
• All scripts source machine-specific tool paths from tools.local.bat (created by setup_tools.bat from tools.local.bat.example): VCVARS, MSBUILD, CMAKE (VS2022's bundled CMake), NINJA (a chocolatey install at C:\ProgramData\chocolatey\bin\ninja.exe, not bundled), NVCC.
• Release vs Debug: both configurations are provided (build_release.bat, build_debug.bat) targeting the same .vcxproj. The real-time constraints — a 512-sample @ 44.1 kHz callback (~11.6 ms) summing the full radiation-weighted mode bank, plus a per-pixel medium compute pass each frame — make the Release build the practical target, but this is an inference from the workload; nothing in code or comments documents the Debug build as non-viable.

This section is the structural spine the other sections hang on, spanning the whole pipeline: shape/mesh (Triangle + triangulator in libmodal) → K/M assembly (Eigen, CPU) → eigensolve (cuSOLVER Ssygvdx dense on GPU for ≤512 DOFs, hand-rolled Eigen shift-invert Lanczos CPU fallback) → modes + per-mode $\zeta_i$ + radiation $\alpha_i$ stored in opaque ModalResult → strike projection (CPU, mutex-guarded) → audio synthesis (CUDA kernel on a PortAudio thread) and visualization (OpenGL mesh shaders pre-NewFrame; GLSL medium compute + ImGui overlays in the ImGui phase). It defines who owns what (engine vs app), how data crosses the C-API (opaque handles + POD query buffers), which work runs on the GPU (eigensolve + audio in CUDA; mesh + medium field in OpenGL), and which libraries provide each capability.

2D triangulation: constrained, quality Delaunay

The 2D mesher solves a planar straight-line graph (PSLG) triangulation problem: given a simple polygon boundary $\partial\Omega$ (and optionally interior hole loops), produce a triangulation $\mathcal{T}$ of $\Omega$ that is conforming, respects the input segments as constrained edges, and satisfies quality bounds. Triangle implements Ruppert's Delaunay refinement: starting from the constrained Delaunay triangulation (CDT) of the PSLG, it inserts circumcenters ("Steiner points") to eliminate (a) triangles whose minimum angle falls below a bound and (b) triangles exceeding a maximum-area constraint $A_{\max}$. The empty-circumcircle (Delaunay) property

maximizes the minimum angle over all triangulations of a fixed point set, which directly benefits the conditioning of the downstream FEM mass/stiffness matrices (well-shaped CST/DKT triangles avoid ill-conditioned element Jacobians). In this code the quality+area switches are pqzea<A>Q: q = quality (Triangle's default $\sim$20$^\circ$ minimum-angle bound), a<A> = per-triangle area cap, p = PSLG, z = zero-indexing, e = emit edge list (used to recover boundary loops). The area cap $A_{\max}$ is the only resolution knob, mapped per level COARSE$\to$EXTREME as $A_{\max}\in\{0.02,\,0.005,\,0.002,\,0.001,\,3\times10^{-4},\,1\times10^{-4}\}$ (m$^2$). Since $A_{\max}\sim \ell^2$, characteristic edge length scales as $\ell\sim\sqrt{A_{\max}}$ and node count as $N\sim |\Omega|/A_{\max}$.

Crystal-lattice (non-Delaunay) interior seeding

An alternative path bypasses Ruppert refinement: interior nodes are placed on a regular grid first, then Triangle is invoked with pzYeQ (no q, no a) so it only computes the CDT of the supplied points (the Y flag preserves segment endpoints so the polygon boundary survives intact). For a square lattice of pitch $s$, points sit at $(x_{\min}+\tfrac{s}{2}+is,\ y_{\min}+\tfrac{s}{2}+js)$; for the triangular/hexagonal lattice (these two enum values share one code path), rows are offset by $s/2$ on alternate rows with row pitch $\Delta y = s\sqrt3/2$ (the equilateral-triangle grid, also the vertex set of the honeycomb). The "graphite" lattice keeps rows regular in $y$ but jitters $x$ by $\pm s/3$ (a deterministic LCG, not seeded by the porosity seed) to mimic anisotropic microstructure. Lattice pitch is $s=\sqrt{A_{\max}}$ (res_to_spacing). Points are kept only if the candidate and its four $\pm$margin axis-neighbors (margin $0.4s$) lie strictly inside $\Omega$ (ray-casting point-in-polygon), preventing slivery boundary triangles. This trades Delaunay angle-optimality for near-regular connectivity, yielding cleaner, more isotropic modal node patterns. Glass (= pure Delaunay) and Pentagon both fall back to the refined Delaunay path (Pentagon has no lattice generator).

Porosity: hole-list carving (2D)

Voids are circular disks $\{(c_x,c_y,r)\}$ scattered to approximate a target void fraction $\phi$ (clamped to $[0,0.35]$). The placement is rejection sampling against three constraints — interior containment, boundary clearance $d_{\partial\Omega}(c)\ge r+m$, and mutual non-overlap $\lVert c_i-c_j\rVert \ge r_i+r_j+m$ with margin $m=0.4\,r_0$ — accumulating area $\sum \pi r_i^2$ until it reaches $\phi\,|\Omega|$ or a 4000-attempt cap. The mean radius is $r_0=\texttt{hole\_size}\cdot\min(L_x,L_y)$ where $\min(L_x,L_y)$ is the smaller bbox extent, each radius jittered $\times U(0.75,1.25)$, drawn from a seeded std::mt19937 for reproducibility. Each accepted disk becomes a closed segment loop wound clockwise (16 segments) plus an interior hole-marker point; Triangle's flood-fill from the marker deletes the enclosed triangles, leaving a perforated-but-connected ligament network whose effective stiffness (and hence modal spectrum) drops with $\phi$.

Shell curvature: Marguerre shallow-shell height field

A flat 2D footprint is turned into a shallow shell by attaching a known mid-surface height $z=f(x,y)$ at every node (closed forms in compute_z_at). With reference radius $R=\tfrac12\min(L_x,L_y)$ and depth $h=0.3\,a\,R$ (amount $a\in[0,1]$), the kinds are:

with $r^2=x^2+y^2$ and $(x,y)$ centered at the bbox origin (the shell solve subtracts the bbox center before calling compute_z_at). This $z$-field is the geometric input to Marguerre's shallow-shell theory, which augments the in-plane (membrane) strain with curvature-coupling terms. For linearly interpolated $z$ over a CST triangle, $\partial_x z,\partial_y z$ are element-constant and the membrane strain acquires

the small-amplitude (linearized von Kármán) limit. This couples bending $w$ to stretching, so the assembled $K\phi=\omega^2 M\phi$ has membrane–bending blocks that stiffen flexural modes — geometric stiffening that distinguishes a gong/bowl from a flat plate. The code emits both the cross-coupling block $K_{mb}$ and an extra positive-definite bending block $K_{bb,\text{extra}}=G_w^{\top}(S_z^{\top}D_m S_z)G_w\,A$ that guarantees curvature raises (never lowers) flexural frequencies. At $z\equiv0$ all coupling vanishes and the shell solve reduces exactly to the flat plate.

3D structured tetrahedral meshing

Solids are meshed without a 3D Delaunay tessellator. A cube is a regular $n^3$ grid split into 6 tets per voxel by the Freudenthal/Kuhn subdivision (six tets sharing the cube main diagonal, guaranteeing face-conformity across cells; explicit tetc table). Spheres/ellipsoids reuse the cube grid and apply the $L^\infty\!\to\!L^2$ shell map $\mathbf{p}\mapsto \mathbf{p}\,\lVert\mathbf{p}\rVert_\infty/\lVert\mathbf{p}\rVert_2$ (maps cube shells to concentric spheres), then per-axis scaling (ELLIPSOID uses $(1,0.7,1.3)$). Bodies of revolution/extrusion (cylinder, cone, pyramid, capsule, torus) stack scaled copies of a 2D base cross-section (a Triangle-meshed disc or a structured square grid) into prism layers, each base-triangle prism split into 3 tets. Winding is normalized to positive signed volume $\tfrac16\det[\,\mathbf{b}-\mathbf{a},\mathbf{c}-\mathbf{a},\mathbf{d}-\mathbf{a}\,]>0$ (push_tet). The boundary surface is extracted combinatorially: a triangular face shared by exactly one tet is a boundary face, oriented outward by the sign of $\mathbf{n}\cdot(\mathbf{x}_{\text{opp}}-\mathbf{x}_A)$. 3D porosity drills straight Z-tunnels by deleting tets whose XY centroid lands in a scattered disk, then keeps only the largest face-connected component (so the free-solid 6-rigid-body-mode invariant is preserved), reverting the whole carve if it would lose $>50\%$ of the tets.

All meshing lives in libmodal/, namespace modal.

Triangle library. Shewchuk's Triangle v1.6 is vendored verbatim at extern/triangle/triangle.{c,h} and compiled with TRILIBRARY, ANSI_DECLARATORS, REAL=double. libmodal/src/mesh/triangulator.cpp includes it under extern "C" and calls ::triangulate(switches, &in, &out, nullptr).

2D triangulation (mesh/triangulator.{h,cpp})

• triangulate(verts, res, lattice) — main entry. Resolution$\to$area map is res_to_area() (lines 57-67). Glass (Delaunay) path uses switches pqzea%.6fQ (line 347); lattice path uses pzYeQ (line 343) with pre-seeded interior from generate_lattice_interior() (lines 193-266) at pitch res_to_spacing() = sqrt(area) (lines 168-170). Effective lattice resolves a global default set by set_default_lattice().
• triangulate_with_holes(verts, holes, res, lattice) (lines 380-450) — builds outer CCW loop + per-hole 16-segment CW loops (kSeg=16, line 388) + one holelist marker per hole; switches pqzea%.6fQ. The lattice argument is ignored (v1 always Delaunay interior).
• triangulate_ring(cx,cy,r_outer,r_inner,n_seg,res) (lines 452-525) — annulus via outer CCW + inner CW loop + center hole marker.
• build_ordered_boundary() (lines 23-54) walks Triangle's marked edges (edgemarkerlist>0) into an ordered boundary loop, stored in Mesh::boundary_nodes.
• Mesh struct (triangulator.h:29-38): nodes ([x,y]), triangles (CCW int3), boundary_nodes, beam_edges (1D couplings drawn by the renderer).

Porosity (same file)

• scatter_holes(polygon, porosity, hole_size_frac, seed) (lines 114-160) — seeded std::mt19937 rejection sampling; cap porosity<=0.35, kMaxAttempts=4000; mean radius r0 = hole_size_frac * min(bbox extents), margin 0.4*r0; helpers point_in_poly, dist_to_boundary (lines 84-112). Returns std::vector<HoleCircle>.
• App-level hook: set_default_porosity(porosity, hole_size, seed) stores a PorositySpec; triangulate() (lines 276-284) auto-routes through triangulate_with_holes only when porosity>0 AND at least one disk was placed. Default {0,0.07f,1} reproduces the solid mesh (regression-guarded in tests/test_porous.cpp).

Preset / arbitrary shapes (mesh/preset.{h,cpp})

• preset_vertices(shape, scale, res) (res unused) returns boundary vertices for the simple presets. Builders: make_rect (square), make_circle (32-gon), make_equilateral, make_l_shape, make_tuning_fork, make_string, make_cross, make_star, make_regular_polygon (PENTAGON/HEXAGON/OCTAGON), make_ellipse, make_semicircle, make_trapezoid, make_rhombus, plus parametric make_heart, make_teardrop, make_crescent (lune via two-circle intersection, lines 169-189). normalize_fit() (lines 15-35) scales raw point lists to fit [0,s]^2 at 0.92 margin and forces CCW via signed area. RING throws (use triangulate_ring()); the many compound/instrument/bell presets also throw and are built directly in modal_api.cpp.
• Arbitrary polygons: modal_create_from_polygon(verts, n, resolution) (modal_api.cpp:382-395) packs raw float pairs and calls triangulate(); default MeshResolution::FINE, with resolution 0/1/2 $\to$ COARSE/MEDIUM/HIGH. This is the path the viewer's free-draw / tuning mode uses.

Shell curvature (modal_api.cpp + fem/)

• compute_z_at(kind, amount, x, y, L_x, L_y) (modal_api.cpp:4109-4153) — closed-form per-kind height; public in modal.h so tests/viewer mirror it. R=0.5*min(L_x,L_y), h=amount*R*0.3; returns 0 for Flat or amount<=0.
• modal_solve shell branch (lines 2882-2907) computes the mesh bbox, subtracts the center, fills a per-node z_field, and calls assemble_KM_shell(mesh, z_field, mat, thickness). Shell is active when curvature_kind != Flat && curvature_amount > 0 (5 DOF/node).
• assemble_KM_shell (fem/assembler_shell.cpp) gathers per-triangle z_vals and emits a 5-DOF/node (u_x,u_y,w,theta_x,theta_y) sparse K/M (reserves 225 triplets/triangle); calls plate_is_thick() to pick DKT vs Mindlin bending.
• shell_element_matrices (fem/element_shell.cpp) composes CST membrane (6 DOF, cst_element_matrices) + DKT or Mindlin bending (9 DOF) into a 15x15 element via the kCstToShell/kDktToShell scatter maps, then adds the Marguerre coupling: element-constant dz_dx/dz_dy from CST gradients b_i,c_i (lines 106-121), strain-gradient map Sz (lines 149-153), cross-coupling K_mb = B_m^T D_m S_z G_w * A (lines 167-177) scattered symmetrically, plus geometric-stiffening K_bb_extra = G_w^T (S_z^T D_m S_z) G_w * A (lines 194-203) so curvature raises (never lowers) flexural frequencies. The w-gradient operator G_w carries only the w DOFs (theta contributions dropped) — a shallow-shell, CST-linear approximation acknowledged in the in-code comment.

3D solids (mesh/solid_mesh.{h,cpp})

• make_solid(kind, size, resolution) dispatches build_cube (Freudenthal 6-tet tetc table, lines 84-101), build_sphere_like (L-inf->L2 map + per-axis scale, lines 104-118; SPHERE = $(1,1,1)$, ELLIPSOID = $(1,0.7,1.3)$), build_cylinder (Triangle-meshed disc extruded into prisms, lines 120-160), and via build_stack/connect_layers/cap_to_apex: build_cone, build_pyramid, build_capsule, build_torus. Disc/square base meshes from make_disc_base/make_square_base; porosity is force-disabled (then restored) while meshing the cross-section.
• connect_layers (lines 37-46) sorts each base-triangle's indices before splitting each prism into 3 tets so adjacent prisms/layers pick the same quad diagonal (the v12.0 crack fix); push_tet flips winding to positive volume.
• extract_surface (lines 56-82) counts faces in a std::map<array<int,3>,Rec> keyed on sorted vertex tuples, keeps count-1 faces, orients outward via the opposite-vertex dot test.
• carve_tunnels(sm, porosity, hole_size, seed) (lines 266-323) reuses scatter_holes over the XY bbox rectangle, drops tets whose XY centroid is inside a disk, then runs a face-adjacency flood-fill keeping only the largest component (reverts if the largest component is <50% of the original tets), reindexes nodes, and re-extracts the surface. modal_create_solid (modal_api.cpp:2760-2773) applies it whenever the global porosity>0.

Thick/thin dispatch. plate_is_thick(nodes, thickness) (fem/element_mindlin.cpp:7-17) returns true when thickness/span >= kThickShearRatio = 0.10 (element_mindlin.h:12; span = larger bbox extent), routing both plate and shell assemblers to the Mindlin transverse-shear bending block. Default thickness=0.001 keeps presets in the thin/DKT regime.

This is the SHAPE→MESH stage, the front of the end-to-end pipeline. A preset/polygon outline (or a parametric solid kind) is converted into a concrete discretization: 2D nodes+triangles (Mesh), an optional per-node z-height field for shells, or 3D nodes+tets+surface (SolidMesh). That discretization is the direct input to the ASSEMBLE stage (CST+DKT/Mindlin plate, Marguerre shell, or tet assemblers build sparse K and M), which feeds the EIGENSOLVE ($K\phi=\omega^2 M\phi$ → modes), then strike projection / modal synthesis / audio, with the same mesh reused for the renderer's animated mode-shape and curvature visualization. Mesh resolution and porosity here set both the spectral accuracy and the qualitative timbre (perforation lowers stiffness and detunes the spectrum).

All structural elements share the small-strain linear-elastic weak form. For a domain $\Omega$ with displacement $u$, strain $\varepsilon$, stress $\sigma=D\varepsilon$ and density $\rho$, the Galerkin discretization $u\approx N u_e$ gives the element matrices

with $B=\partial N$ the strain-displacement operator. Assembling all $K_e,M_e$ and substituting $u(x,t)=\phi\,e^{i\omega t}$ into $M\ddot u+Ku=0$ yields the generalized symmetric eigenproblem $K\phi=\omega^2 M\phi$ solved downstream (solve_eigenmodes, eigensolve.cpp). The solver dispatches on size: for $\mathrm{DOF}\le 512$ it calls a dense GPU path (solve_dense, eigensolve_dense.cu, single-precision cusolverDnSsygvdx with CUSOLVER_EIG_TYPE_1, lower fill, range-I lowest indices), falling back on exception to a hand-written sparse path (solve_sparse, eigensolve_sparse.cpp — a single-vector thick-restart Lanczos with shift-invert $\sigma=10^{-4}$ factorized by Eigen SimplicialLDLT; explicitly not ARPACK/Spectra). Presets are typically free-free, so $K$ is singular with rigid-body modes that both solvers filter by $\omega<$ rigid_threshold; Dirichlet rows are eliminated only when fix_level>0 or explicit fixed_nodes are set (modal_api.cpp).

1. CST membrane / plane stress (cst_element_matrices)

Governing PDE: 2D Navier elasticity $\nabla\!\cdot\sigma=\rho\ddot u$ with plane-stress reduction. Linear shape functions $N_i=L_i$ (area coordinates) make $\varepsilon$ constant per element. With node-cyclic coefficients $b_i=y_j-y_k$, $c_i=x_k-x_j$ and area $A$,

$K_e=tA\,B^\top D B$ (closed-form, single point — $B$ is constant). The mass is the consistent CST mass $M_e=\frac{\rho t A}{12}$ with the value $2$ on each node-self $2\times2$ diagonal block and $1$ on the off-diagonal blocks (block-diagonal in $x,y$). The orthotropic variant (cst_element_matrices_ortho) builds the plane-stress $D_{\text{principal}}$ from $(E_x,E_y,G_{xy},\nu_{xy})$ with $\nu_{yx}=\nu_{xy}E_y/E_x$, $Q_{11}=E_x/(1-\nu_{xy}\nu_{yx})$, $Q_{12}=\nu_{xy}E_y/(1-\nu_{xy}\nu_{yx})$, $Q_{22}=E_y/(1-\nu_{xy}\nu_{yx})$, $Q_{33}=G_{xy}$, then rotates it to the element frame by a Voigt strain transformation $D_{\text{local}}=T_m^\top D_{\text{principal}} T_m$ (with $T_m=T(-\theta)$).

2. DKT thin-plate bending (dkt_element_matrices)

Governing PDE: Kirchhoff biharmonic $D_b\nabla^4 w=\rho h\,\ddot w$ with $D_b=Eh^3/12(1-\nu^2)$. DOFs $(w,\theta_x,\theta_y)$. The kernel evaluates the Batoz (1980) DKT rotation-interpolation functions $H^x_k(L_1,L_2,L_3)$, $H^y_k$ (the in-file Hx[9],Hy[9]) with edge constants $P_k=-6a_k/\ell_k^2$, $q_k=3a_kb_k/\ell_k^2$, $r_k=3b_k^2/\ell_k^2$, $t_k=-6b_k/\ell_k^2$ (here $a_k,b_k$ are edge $x,y$ differences and $\ell_k^2$ the edge length squared). The curvature $\kappa=(\kappa_{xx},\kappa_{yy},2\kappa_{xy})=B u$ is assembled by analytic chain-rule differentiation of $H^x,H^y$ through the area coordinates (the $H$ are linear in $L_1,L_2,L_3$, so $d/dx,d/dy$ replace $L_i\to dL_i/dx,\,dL_i/dy$), accumulated at the 3-point interior Gauss rule (weights $1/3$):

Mass is lumped: $\rho h A/3$ on each $w$ DOF and a rotary inertia $\rho h^3 A/36$ on each rotation DOF. Note $\nu$ is hardcoded to $0.3$ inside the bending $D_b$ (the scalar prefactor $D_b$ passed in does carry the true $\nu$, but the $[1,\nu;\nu,1;\dots]$ shape matrix uses $0.3$). The file's top comment misleadingly labels the element "BCIZ" / "Tocher"; the code that actually runs uses the Batoz $H^x/H^y$.

3. Mindlin-Reissner thick plate (mindlin_element_matrices)

Adds first-order transverse shear: $K_e=K_{\text{bend}}+K_{\text{shear}}$. Bending uses linear rotations $\Rightarrow$ constant curvature $B_b$, with $\kappa_{xx}=\sum b_i\theta_{y,i}$, $\kappa_{yy}=-\sum c_i\theta_{x,i}$, $\kappa_{xy}=\sum(c_i\theta_{y,i}-b_i\theta_{x,i})$ (here $b_i,c_i$ are the area-coordinate gradients $dL_i/dx,\,dL_i/dy$). The shear strain $\gamma=\nabla w+\beta$ (with $\beta_x=\theta_y$, $\beta_y=-\theta_x$) uses the linear $w$-gradient plus rotations interpolated by $N_i=L_i$, integrated at 3 Gauss points. Shear locking is suppressed by the Lyly-Stenberg-Vihinen factor:

where $h_K$ is the longest edge. This rank-sufficient linear shear (the comment notes a pure 2-DOF constant-$\gamma$ DSG-style shear would leave one spurious zero-energy mode) recovers Kirchhoff thin-plate behavior as $h\to0$. Mass is lumped identically to DKT, and Mindlin uses the true material $\nu$. Dispatch is by plate_is_thick: thickness/span $\ge$ kThickShearRatio $=0.10$.

4. Marguerre shallow shell (shell_element_matrices)

15-DOF (5/node: $u_x,u_y,w,\theta_x,\theta_y$). It composes the CST membrane (6 DOF) + DKT or Mindlin bending (9 DOF) block-diagonally (scatter maps kCstToShell,kDktToShell), then adds Marguerre bending-stretching coupling from a per-node height field $z(x,y)$. With constant element slopes $\partial z/\partial x=\sum b_i z_i$, $\partial z/\partial y=\sum c_i z_i$, the linearized von-Kármán/Marguerre membrane strain gains $\varepsilon\mathrel{+}=S_z\nabla w$, giving cross-coupling $K_{mb}=A\,B_m^\top D_m S_z G_w$ (scattered symmetrically) and an extra bending stiffness $K_{bb}^{\text{ext}}=A\,G_w^\top(S_z^\top D_m S_z)G_w$ so curvature stiffens the shell. $G_w$ carries only the $w$ DOFs (CST-linear, $\theta$ zeroed) — a deliberate shallow-shell simplification.

5. Linear tetrahedron (tet_element_matrices)

3D isotropic elasticity, $N_i$ linear $\Rightarrow$ constant strain. Shape-function gradients $b_i,c_i,d_i$ come from inverting the $[1\,x\,y\,z]$ matrix; the $6\times12$ $B$ uses strain order $(\varepsilon_{xx},\varepsilon_{yy},\varepsilon_{zz},\gamma_{xy},\gamma_{yz},\gamma_{zx})$; full 3D Lamé $D$ with $\lambda=E\nu/[(1+\nu)(1-2\nu)]$, $\mu=E/2(1+\nu)$. $K_e=V\,B^\top D B$, lumped mass $\rho V/4$ on each translational DOF.

6. Euler-Bernoulli beam (beam_element_matrices)

Axial bar + Hermite cubic bending, plus a geometric stiffness from axial tension $T$ ($\tfrac{6}{5}T/L$, $T/10$, $\tfrac{2}{15}TL$, $\dots$ on the transverse/rotational DOFs) that recovers ideal-string physics $f=\tfrac{n}{2L}\sqrt{T/\mu}$ when $TL^2\gg EI$. Consistent Hermite mass (the $156,22L,54,-13L,\dots$ pattern $\times\rho AL/420$ for bending, $\rho AL/\{3,6\}$ for axial); rotated to global by $T^\top K_l T$, $T^\top M_l T$.

7. Acoustic P1 triangle (acoustic_element_matrices)

Scalar Helmholtz/wave operator for air pressure: $K_e=\int\nabla N\!\cdot\!\nabla N=\tfrac{1}{4A}(b_ib_j+c_ic_j)$ and a consistent $M_e=\tfrac{A}{12c^2}\,[\,2\ \text{on diagonal},\,1\ \text{off}\,]$, so the generalized problem $K_a p=\omega^2 M_a p$ (with $M_a\propto1/c^2$) yields cavity angular frequencies; e.g. a closed rectangular duct returns $\omega_{10}/2\pi\approx c/2L_x$ (test_acoustic_element.cpp).

Element kernels live in libmodal/src/fem/element_*.cpp (declared in libmodal/include/fem/element_*.h), each returning fixed-size Eigen Matrix<double,n,n> $K_e,M_e$ that the matching assembler scatters into Eigen sparse triplets. Material is modal::MaterialParams (libmodal/include/modal/types.h): E=200e9, nu=0.30, rho=7800 defaults, plus orthotropy fields (anisotropic, E_y_ratio (default 0.10), G_xy_override, grain_angle_rad).

• CST: cst_element_matrices / cst_element_matrices_ortho (element_cst.cpp). Closed-form $tA\,B^\top D B$; consistent mass via the $\frac{\rho tA}{12}\{2\text{ self},1\text{ cross}\}$ block pattern (lines 49-56). Orthotropic path builds D_principal (lines 103-112) and rotates with Tm (lines 135-141). Assembled by assemble_KM / assemble_KM_multimat (assembler.cpp), 2 DOF/node, gdof = {2*tri[i], 2*tri[i]+1}; orthotropy selected by mat.anisotropic.
• DKT: dkt_element_matrices (element_dkt.cpp). Batoz Hx[9],Hy[9] (lines 149-169) differentiated analytically into dHx_dx...dHy_dy to build the $3\times9$ curvature B; 3-point Gauss accumulate Ke += B^T*DbMat*B*A*(1/3) (line 277); nu=0.3 hardcoded in the Db shape matrix (line 98). Lumped mass: rho*h*A/3 on $w$, rho*h^3/12*A/3 = rho*h^3*A/36 on each rotation (lines 283-289).
• Mindlin: mindlin_element_matrices (element_mindlin.cpp); constant-curvature Bb (lines 49-55), stabilized linear shear Bs with Ds = kappa*G*h*h^2/(h^2+alpha*hK^2), kappa=5/6, alpha=0.1 (lines 70-89). plate_is_thick + kThickShearRatio=0.10f (element_mindlin.h) drives DKT-vs-Mindlin dispatch; uses the real nu.
• Plate assembly: assemble_KM_plate (assembler_plate.cpp), 3 DOF/node, branches on plate_is_thick(mesh.nodes, thickness).
• Shell: shell_element_matrices (element_shell.cpp); kCstToShell[6]/kDktToShell[9] scatter the sub-blocks into the 15-DOF layout (lines 15-16); the bending sub-block is DKT or Mindlin depending on the thick_bending flag (lines 36-41); then Kmb[6][9] (line 167) and Kbb_ext[9][9] (line 194) add the Marguerre coupling. assemble_KM_shell (assembler_shell.cpp) takes a per-node z_field; the $z(x,y)$ closed forms per CurvatureKind (Flat/Dome/Cone/Cylinder/Saddle/Bowl) live in modal_api.cpp::compute_z_at (lines 4109+).
• Tet: tet_element_matrices (element_tet.cpp); Cinv of the $[1\,x\,y\,z]$ matrix gives gradients (lines 10-18); Lamé D (lines 32-37); lumped diagonal mass rho*V/4 (lines 42-44). assemble_KM_solid (assembler_solid.cpp) emits only the diagonal of M. Winding normalized to positive volume upstream by push_tet/signed_vol6 in mesh/solid_mesh.cpp.
• Beam: beam_element_matrices (element_beam.cpp); local Kl axial+Hermite (lines 30-39), tension geometric block (lines 47-59), Hermite consistent Ml (lines 61-78), rotation T, transform Ke=T^T*Kl*T (line 93). assemble_KM_mixed (assembler_mixed.cpp) overlays beams onto CST membrane with a shared 3-DOF/node $(u_x,u_y,\theta)$ layout — CST fills only $(u_x,u_y)$, the beam adds $\theta$; cross-section A=t*width, I=t*width^3/12.
• Acoustic: acoustic_element_matrices + assemble_acoustic_KM (element_acoustic.cpp), winding-invariant (|area|), with pressure-release Dirichlet via row/col elimination at port_nodes (lines 53-67). Cavity meshes are solved through solve_cavity_modes (acoustic subtree).
• Eigensolve: solve_eigenmodes (eigensolve.cpp) dispatches dense (eigensolve_dense.cu, cusolverDnSsygvdx) for K.rows() <= 512 (with sparse fallback on exception) else sparse (eigensolve_sparse.cpp, Lanczos+SimplicialLDLT).

This is the "shape → mesh → ASSEMBLE → eigensolve" stage. Mesh nodes/triangles (from triangulator.cpp) or tets (solid_mesh.cpp) plus a MaterialParams are fed to the appropriate assemble_KM_*, which is selected in modal_api.cpp::modal_solve by physical regime (membrane CST, plate DKT/Mindlin, Marguerre shell, tet solid, mixed string/box membrane+beam, or air cavity). The resulting sparse $(K,M)$ feed solve_eigenmodes (eigensolve.cpp), returning $\omega_i$ and mode shapes $\phi_i$; those become the modal oscillator bank (libmodal/src/modal/oscillator.cu, audio path audio/synthesis.cu) that strike events drive, and the eigenvectors drive both audio synthesis and the deformation visualization.

The material is the SI-unit struct MaterialParams (types.h:77). Isotropic defaults are steel-like: $E=200\,\mathrm{GPa}$, $\nu=0.30$, $\rho=7800\,\mathrm{kg/m^3}$, base damping $\zeta=0.008$. Geometry contributes a thickness $t$ (passed separately to the assembler) and, for 1-D beam/string elements, an axial pre-tension $T$ stored per BeamEdge together with the beam width (default $0.01\,$m).

Plane-stress membrane (CST)

The 2-D filled shapes use a constant-strain triangle. The isotropic plane-stress constitutive matrix and element stiffness are

with the consistent mass $\mathbf{M}_e=\tfrac{\rho tA}{12}\,(2\,\mathbf I$ on diagonal node blocks, $1\,\mathbf I$ off-diagonal$)$ (element_cst.cpp:38-56). Crucially $\mathbf K_e\propto E$ and $\mathbf M_e\propto\rho$ are exactly linear: $\nu$ and the geometry set the dimensionless shape of the spectrum, $E/\rho$ sets its absolute scale.

Orthotropy / grain

When anisotropic is true, the along-grain modulus is $E_x=E$ and across-grain $E_y=E\cdot\texttt{E\_y\_ratio}$ (wood default ratio $0.10$). The principal-axis plane-stress matrix uses the engineering Poisson symmetry $\nu_{yx}=\nu_{xy}E_y/E_x$ (with nu interpreted as $\nu_{xy}$):

with auto shear modulus $G_{xy}=(E_x+E_y)/\big(4(1+\nu_{xy})\big)$ (which collapses to the isotropic $E/(2(1+\nu))$ when $E_x=E_y$) unless G_xy_override$>0$. The code clamps $E_y\ge10^{-3}$ and the denominator $1-\nu_{xy}\nu_{yx}\ge10^{-6}$ to keep the inverted plane-stress form non-singular. A grain rotation $\theta=\texttt{grain\_angle\_rad}$ rotates $\mathbf D$ into the element frame via a Voigt strain-transformation $\mathbf D_{\text{loc}}=\mathbf T^{\mathsf T}\mathbf D_{\text{princ}}\mathbf T$ (element_cst.cpp:122-142; the $\theta=0$ branch skips the rotation). Mass is grain-independent (density is isotropic), so orthotropy reshapes the spectrum directionally without touching $M$.

Thin plate (BCIZ triangle / Kirchhoff)

Out-of-plane bending uses the flexural rigidity $D_b=E\,h^{3}/\big(12(1-\nu^{2})\big)$, computed in assemble_KM_plate and passed to the bending element. Note the bending element is a simplified DKT-style BCIZ (Bazeley-Cheung-Irons-Zienkiewicz) conforming thin-plate triangle, not the full Batoz DKT (see the source comment at element_dkt.cpp:6-8). Plate stiffness carries the cubic thickness dependence $K_b\propto E h^3$; the element uses a lumped mass matrix with translational entries $\propto\rho h$ and a rotary inertia term $\rho h^{3}/12$ on the $\theta_x,\theta_y$ DOFs. Thick plates auto-switch to a Mindlin (shear-deformable) element via plate_is_thick.

Ideal string via geometric tension stiffness

A 1-D beam element (Euler-Bernoulli, Hermite shape functions) has bending stiffness $\propto EI/L^3$ and axial $\propto EA/L$, with cross-section $A=t\cdot w$ and $I=t w^3/12$. Adding pre-tension $T$ contributes the geometric (stress-stiffening) matrix on the transverse/rotational DOFs, with diagonal entries $K^{g}_{vv}=\tfrac{6}{5}\tfrac{T}{L}$, $K^{g}_{\theta\theta}=\tfrac{2}{15}TL$ and couplings $\pm T/10$, $-\tfrac{1}{30}TL$ (element_beam.cpp:47-59). This is the discrete analogue of $+T\,\partial^2 w/\partial x^2$ in the string PDE, so when $TL^2\gg EI$ the bending term is negligible and the element reproduces ideal-string physics.

Because $\mathbf K=E\,\hat{\mathbf K}(\text{geom},\nu)$ and $\mathbf M=\rho\,\hat{\mathbf M}(\text{geom})$ with $\hat{\mathbf K},\hat{\mathbf M}$ purely geometric, the Rayleigh quotient gives every eigenvalue $\lambda_i=\omega_i^2=\phi^{\mathsf T}\mathbf K\phi/\phi^{\mathsf T}\mathbf M\phi\propto E/\rho$. Hence the universal law

where the dimensionless eigenvalue $c_i$ depends only on shape, boundary conditions and $\nu$. Specializing:

• Membrane / in-plane (2-D shapes): the in-plane wave speed scales as $\sqrt{E/\rho}$ (plane-stress longitudinal speed $c=\sqrt{E/(\rho(1-\nu^2))}$), so $f\propto\frac{1}{L}\sqrt{E/\rho}$ — halving the size doubles every frequency.
• Thin-plate bending: with $D_b=Eh^3/12(1-\nu^2)$ and areal mass $\rho h$, $f_i=\frac{c_i}{2\pi}\sqrt{D_b/(\rho h L^4)}\Rightarrow f\propto\frac{h}{L^{2}}\sqrt{\frac{E}{\rho\,(1-\nu^2)}}$. Bending frequencies are linear in thickness and inverse-square in size — the signature of plates/bells.
• Ideal string (tension-dominated beam): $f_n=\frac{n}{2L}\sqrt{T/\mu}$ with linear mass $\mu=\rho A$. This is recovered by the beam's geometric stiffness in the limit $TL^2\gg EI$; the residual $EI$ adds inharmonicity (slightly stretched overtones, like a real piano).
• Curvature shift (shells): the Marguerre shallow-shell solver couples membrane and bending energy, raising low modes (a dome is stiffer than a flat plate of equal $h$), so curvature shifts and reorders the spectrum without changing $E/\rho$.

Damping is applied after the eigensolve as a diagonal modal-damping ratio — there is no global Rayleigh $\alpha\mathbf M+\beta\mathbf K$ matrix; each modal oscillator is independent. With the frequency-dependent toggle on, $\zeta_i=\zeta+\text{slope}\cdot\omega_i/(2\pi\cdot 1000)$, $\zeta_i\in[10^{-5},0.95]$, where $\omega_i$ is the angular eigenfrequency in rad/s (note res->frequencies stores $\omega=\sqrt\lambda$, not Hz, so dividing by $2\pi\cdot1000$ normalizes to "$\zeta$ added per kHz"). Defaults: $\zeta=0.008$, slope $=0.0005$. When the toggle is off the effective slope is forced to $0$ and every mode shares $\zeta$ (legacy uniform damping for A/B comparison). Each mode is a damped harmonic oscillator $\ddot q_i+2\zeta_i\omega_i\dot q_i+\omega_i^2 q_i=0$, whose impulse response $q_i(t)=A_i\,e^{-\zeta_i\omega_i t}\,\sin(\omega_{d,i}\,t)$, $\omega_{d,i}=\omega_i\sqrt{1-\zeta_i^2}$ has an amplitude envelope decaying with

Because $\tau_i\propto 1/(\zeta_i\omega_i)$ and $\zeta_i$ itself grows with $\omega_i$, the high overtones decay doubly fast: a mode at $4\,$kHz has both a larger $\omega_i$ and a larger $\zeta_i$ than the fundamental. The audible consequence is a bright transient attack that mellows within tens of milliseconds toward a near-pure fundamental — the natural timbre evolution that uniform $\zeta$ cannot produce. The clamp $\zeta_i\le0.95$ keeps $\sqrt{1-\zeta_i^2}$ real (no NaN in the kernel) and the floor $\ge10^{-5}$ prevents undamped ringing.

For the 2-D (in-plane/membrane) solve, each mode's contribution to the radiated pressure is weighted by a discrete-monopole boundary-integral coefficient

i.e. medium density times angular frequency times the net normal-displacement volume flux summed over boundary nodes, where $\hat{\mathbf n}_j$ is the averaged outward in-plane boundary normal and $A_j$ is the half-length of the adjacent boundary segments. The two in-plane displacement components $\boldsymbol\phi_i(j)=(\phi_i[2j],\phi_i[2j{+}1])$ are projected onto $\hat{\mathbf n}_j$ (radiation.cpp:58-65); the result is normalized to $\max_i|\alpha_i|=1$. The $\omega_i$ factor encodes that, at equal displacement, faster modes radiate more (volume velocity $\propto\omega$); modes that breathe (net outward flux) radiate, modes whose normal motion cancels do not. The realtime audio sum is $p(t)=\sum_i \alpha_i\sum_k A_{i,k}\,e^{-\zeta_i\omega_i(t-t_{0,k})}\sin(\omega_{d,i}(t-t_{0,k}))$, evaluated one thread per sample on the GPU. Downstream the viewer applies (when post-processing is enabled) a master gain, a 1-pole low-pass, a $\tanh$ soft-clip, and an auto-normalizer that holds the peak near $0.9$, plus an optional noise gate.

The MaterialParams struct (types.h:77) carries $E=200\mathrm{e}9$, $\nu=0.30$, $\rho=7800$, $\zeta=0.008$, zeta_freq_enabled=true, zeta_freq_slope=0.0005, and orthotropy fields; MediumParams defaults speed_of_sound=343, medium_density=1.21. It is set via modal_set_material (copies into obj->mat) and the per-element palette modal_set_element_materials. The $K\propto E$, $M\propto\rho$ assembly lives in cst_element_matrices (consistent mass), cst_element_matrices_ortho (Voigt-rotated principal $\mathbf D$), assemble_KM_plate (BCIZ bending, lumped mass), and beam_element_matrices (geometric tension block). The per-mode $\zeta_i$ is computed identically in three places — the 2-D solve (modal_api.cpp:3299-3309), the solid solve (:2838-2847), and the runtime update modal_set_damping (:3784-3791) — and uploaded to GPU synth.d_zetas. The audio kernel audio_sum_kernel (synthesis.cu:44-66) builds $\omega_d$, evaluates the decaying sinusoid, and applies the freq_shift and isfinite guard. The radiation coefficients come from compute_radiation_coefficients (radiation.cpp:8), uploaded to synth.d_alphas.

Governing eigenproblem

After FE assembly the free, undamped structural vibration problem is the generalized symmetric-definite eigenproblem

with $K=K^{\!\top}\succeq 0$ the global stiffness and $M=M^{\!\top}\succ 0$ the global consistent mass matrix (assembled from the full element mass $M_e$; the CST element mass is $m=\tfrac{\rho t A}{12}$ with $2m$ on same-node $2\times2$ diagonal blocks and $m$ off-diagonal — see cst_element_matrices, not lumped). The natural frequencies are $f_i=\omega_i/2\pi=\sqrt{\lambda_i}/2\pi$ and $\{\phi_i\}$ are the mode shapes. The DOF count per node is selected per solve as shell?5 : plate?3 : mixed?3 : 2: $2$ (CST membrane $u_x,u_y$), $3$ (membrane+beam with a drilling DOF, or DKT plate $w,\theta_x,\theta_y$), $5$ (Marguerre shallow shell $u_x,u_y,w,\theta_x,\theta_y$), and $3$ for the tetrahedral solid path.

Mass orthonormalization

Both solvers return eigenvectors $M$-orthonormalized,

i.e. unit generalized (modal) mass $m_i=\phi_i^{\!\top}M\phi_i=1$ and modal stiffness $k_i=\omega_i^2$. This is the normalization cuSOLVER's ?sygvdx produces for EIG_TYPE_1, and the one the Lanczos recurrence enforces explicitly via the $M$-inner product $\langle u,v\rangle_M=u^{\!\top}Mv$ (the basis $V$ is $M$-orthonormal). Downstream amplitude scaling depends on it: the strike projection (compute_strike_impulse) reads the modal participation directly from the eigenvector at the struck node — for plate/shell the raw out-of-plane $w$-DOF value $\phi_i(\text{w-DOF})$, for in-plane CST/mixed the directional projection $\phi_i\!\cdot\!\hat n$ onto the boundary normal $\hat n$ (or $+\hat y$ for interior strikes) — which is only the correct relative modal force amplitude when modal mass is unity.

Dense path (cuSOLVER, partial, float32)

For $n=\dim K\le 512$, solve_dense densifies $K,M$ to column-major float32 (lower triangle) and calls cusolverDnSsygvdx with CUSOLVER_EIG_TYPE_1 ($Kx=\lambda Mx$), CUSOLVER_EIG_MODE_VECTOR, CUSOLVER_EIG_RANGE_I requesting indices $[1,\,\min(n,\text{num\_modes}+3)]$, and CUBLAS_FILL_MODE_LOWER. Algorithmically this is the LAPACK ?sygvd/?sygvx-style reduction to standard form: Cholesky factor $M=LL^{\!\top}$, form $C=L^{-1}KL^{-\top}$, tridiagonalize and diagonalize $C$, then back-transform $\phi=L^{-\top}y$. The "+3" over-request gives slack so the rigid-body filter does not eat into the requested elastic count. float32 Cholesky is the precision weak point: the dispatcher wraps solve_dense in try/catch and treats any thrown std::exception (typically cusolverDnSsygvdx reporting a nonzero info, i.e. failed factorization on ill-conditioned fine meshes) as a signal to fall back to the sparse path.

Sparse path (shift-invert Lanczos, double)

For $n>512$ (or after a dense throw), solve_sparse runs a shift-invert Lanczos for the smallest eigenvalues. With spectral-transform operator

an eigenpair $(\lambda,\phi)$ of the original problem maps to an eigenpair of $\mathrm{OP}$ with eigenvalue $\mu=\dfrac{1}{\lambda-\sigma}$, so the smallest $\lambda$ become the largest, well-separated $\mu$ — the regime Lanczos converges fastest in. The shift $\sigma>0$ also regularizes the near-singular $K$ of free–free (unconstrained) models. The $M$-orthogonal Lanczos recurrence

builds a tridiagonal $T=\mathrm{tridiag}(\beta_{j-1},\alpha_j,\beta_j)$ over a Krylov basis of size $\text{ncv}=\min(2\,\text{nev}+1,n)$, $\text{nev}=\min(\text{num\_modes}+5,\,n-1)$. Ritz values come from SelfAdjointEigenSolver on $T$; processing $\mu$ from largest to smallest, $\lambda=1/\mu+\sigma$, $\omega=\sqrt{\lambda}$, Ritz vectors $\phi=V\,S_{:,i}$. Stability measures: full reorthogonalization against all prior basis vectors each step (a single Gram–Schmidt sweep over $k\le j$ in the $M$-inner product), a fixed-seed (mt19937(42)) $M$-normalized random start vector, and lucky-breakdown handling that injects a fresh $M$-orthogonal random vector when $\beta_j<10^{-14}$. Despite the source comment naming it 'thick-restart', the implementation performs a single Krylov build of size ncv with no restart.

Rigid-body filtering and zero-energy regularization

Elastic modes are separated from rigid-body / spurious zero modes by the threshold $\omega<\texttt{rigid\_threshold}$ (default $1\,\mathrm{rad/s}$), counted into num_rigid (whole-system path only) and dropped. The free 3D solid is special-cased: modal_solve calls solve_eigenmodes(\dots, num\_modes+6, \text{rigid\_threshold}=0.0) — i.e. it disables threshold filtering for the solid entirely — then unconditionally erases the first 6 modes, because a free tet body has exactly 6 rigid modes and float32 leaves them at small but nonzero $\omega$ that an absolute threshold cannot catch reliably. The sparse solver additionally regularizes DOFs with simultaneously (near-)zero stiffness and mass — adding $\max(\varepsilon_{\max}\cdot10^{-10},10^{-10})$ to both diagonals — to keep the shift-invert factorization $\mathrm{SimplicialLDLT}(K-\sigma M)$ non-singular.

Files / entry points

• libmodal/src/solver/eigensolve.cpp — solve_eigenmodes(K, M, num_modes, rigid_threshold=1.0): strategy dispatch. DENSE_DOF_THRESHOLD = 512. For $K.\text{rows()} \le 512$ it tries solve_dense inside a try/catch; on any std::exception (expected float32 Cholesky/info!=0 failure on fine meshes) it silently falls through to solve_sparse. Larger systems go straight to sparse.
• libmodal/src/solver/eigensolve_dense.cu — solve_dense: builds dense column-major float32 lower-triangle h_K,h_M; cusolverDnSsygvdx_bufferSize/cusolverDnSsygvdx with CUSOLVER_EIG_TYPE_1, EIG_MODE_VECTOR, EIG_RANGE_I, FILL_MODE_LOWER, il=1, iu=min(n, num_modes+3). Eigenvalues $\to$ omega=sqrt(w_i), eigenvectors read back column-major from d_K. Manual RAII-style CUDA cleanup with a catch-all that frees and re-throws. A nonzero info from cusolverDnSsygvdx is turned into a std::runtime_error.
• libmodal/src/solver/eigensolve_sparse.cpp — solve_sparse: shift-invert Lanczos described in Theory; Eigen::SimplicialLDLT<SparseMatrix<double>> factorization of $K-\sigma M$ ($\sigma=10^{-4}$), Eigen::SelfAdjointEigenSolver on the dense tridiagonal $T$. Ritz pairs processed largest-$\mu$-first into res, then re-sorted ascending by $\omega$.
• EigenResult (libmodal/include/solver/eigensolve_dense.h): std::vector<double> frequencies ($\omega_i$ in rad/s, ascending, elastic only), std::vector<std::vector<double>> mode_shapes (mode_shapes[i][dof]), int num_rigid.

Boundary conditions (modal_solve, libmodal/src/modal_api.cpp L2928–2979): fixed-DOF selection from explicit obj->fixed_nodes, else fix_level ($\ge4$ = all boundary_nodes; else nfix=fix_level evenly spaced boundary nodes via idx=(i*nb)/nfix). Imposition is a symmetric penalty + row/col zeroing: for each fixed DOF, zero its row and column in both $K$ and $M$, then set $K_{dd}=\text{penalty}=10^6\cdot\max|K_{ij}|$ and $M_{dd}=1$. This keeps the matrices symmetric and pushes constrained modes to $\omega=\sqrt{\text{penalty}}$, far above the audio band (not explicitly filtered, just out of range).

Number of modes kept: caller-supplied. The viewer uses modal_solve(obj, 24) for the initial/instrument solve (app.cpp L150, L6736) and an adaptive num_modes = std::min(std::max(24, n_mesh_nodes/4), 120) for the async resolve (app.cpp L847/L865). The spectrogram visualizer caps its mode-power history at MAX_MODES = 32 (viz_waveform.h). The solid path requests num_modes+6 then drops 6 rigid modes by count.

Per-component solve (solve_eigenmodes_per_component, eigensolve.cpp; dispatched from modal_solve ~L3019): node components come from build_node_components (flood-fill DFS over triangle edges + the mirrored mesh.beam_edges, modal_api.cpp L3366). DOFs are grouped per component; each component's $K_c,M_c$ sub-block is extracted via InnerIterator + an unordered_map global→local index map, solved independently with solve_eigenmodes for up to num_modes modes, eigenvectors scattered back into full-length zero-padded vectors, then all components' modes merged and sorted by ascending $\omega$, truncated to num_modes; out.num_rigid is set to 0 on this path. Used only when n_components>1 and air-cavity coupling is inactive (!obj->cavity.enabled); otherwise the whole-system solve_eigenmodes runs (cavity chambers are intentionally coupled). Verified by tests/test_component_solve.cpp (two degenerate disjoint spring chains must yield localized, not hybridized, eigenvectors: $\text{lo}/\text{hi}<10^{-6}$; and n_components<=1 must match the whole-system solve).

Mass matrix: consistent CST mass in cst_element_matrices / cst_element_matrices_ortho (element_cst.cpp L47–55), $m=\rho t A/12$, diagonal blocks $2m$, off-diagonal $m$; assembled element-wise into the global $M$ (assembler.cpp, full $M_e$, no lumping).

Downstream use: per-mode damping $\zeta_i=\zeta+\text{slope}_{\text{eff}}\cdot\omega_i/(2\pi\cdot1000)$ with $\text{slope}_{\text{eff}}=\texttt{zeta\_freq\_slope}$ only if zeta_freq_enabled (else 0), clamped to $[10^{-5},0.95]$; frequencies/zetas/alphas uploaded to GPU synthesis; strike force = directional mode-shape value at the struck DOF (assumes unit modal mass). modal_compute_hf_width_derivs re-uses stored eigenvectors in a central-difference Rayleigh-quotient shape gradient: it perturbs node $x$ by $(1+\varepsilon)$, re-assembles $K_1,M_1$, and forms $\lambda=\phi^{\!\top}K\phi/\phi^{\!\top}M\phi$ at both geometries using the helper $\texttt{rayleigh\_quotient\_num}(A,\phi)=\phi^{\!\top}A\phi$ (which returns the single quadratic form; the quotient is assembled by the caller).

Validation: tests/test_eigensolve.cpp checks a free-free bar's first axial frequency vs analytical $\omega=\pi c/L$, $c=\sqrt{E/\rho}$, within 15% (CST accuracy), and that frequencies are positive, sorted, and eigenvectors have length $2\,n_\text{nodes}$.

This is the eigensolve/modal-reduction stage, sitting between FE assembly and audio synthesis: shape → mesh (triangulator / solid_mesh) → element + global assembly ($K,M$ via assembler*.cpp) → apply Dirichlet BCs + solve $K\phi=\omega^2M\phi$ (this section) → keep lowest ~24–120 elastic modes, filter rigid bodies, localize per component → per-mode damping $\zeta_i$ and radiation weights $\alpha_i$ → strike projects an impulse onto the modal basis → GPU bank of decaying modal oscillators synthesizes audio and drives the visualizers. The mode set $(\omega_i,\zeta_i,\phi_i)$ produced here is the compact state that every real-time path (audio kernel, spectrogram, strike, auto-tuning gradients) consumes.

The continuous structure is discretized to the second-order system $M\ddot{u} + C\dot{u} + K u = F(t)$, where $u\in\mathbb{R}^{n_{\text{dof}}}$ are nodal DOFs. The generalized eigenproblem $K\phi_i=\omega_i^2 M\phi_i$ yields mass-orthonormal modes $\phi_i^\top M\phi_j=\delta_{ij}$, $\phi_i^\top K\phi_j=\omega_i^2\delta_{ij}$. Writing $u=\sum_i \phi_i q_i$ and left-multiplying by $\phi_i^\top$ decouples the system (with proportional/modal damping $\zeta_i$):

A strike is an impulsive point force $F(t)=\hat F\,\delta(t-t_0)$ applied at a single node in a single direction $\hat n$. The modal force is the projection $f_i = \phi_i^\top \hat F = \phi_i(\text{struck node})\cdot\hat n \cdot \text{magnitude}$. The code treats this as the new amplitude of an underdamped free response with $q_i(0)=0$, so each mode rings the closed-form decaying sinusoid

Crucially the implementation sets $A_i=f_i\cdot\text{magnitude}$ directly (project_strike, eval_single_mode); it does not divide by $\omega_{d,i}$. A rigorous unit impulse gives $\dot q_i(0^+)=f_i$, $A_i=f_i/\omega_{d,i}$; omitting that division gives high modes a flatter excitation spectrum than a true unit impulse — a deliberate timbral choice rather than a unit-consistent impulse response. Striking on a nodal line ($\phi_i\approx 0$ there) gives $f_i\approx 0$ and that mode stays silent — the physically-correct behaviour the viewer relies on (e.g. hitting a plate at its center kills antisymmetric modes).

The projected force depends on which DOF and which direction the strike couples to. compute_strike_impulse (modal_api.cpp:3531) dispatches on dofs_per_node and plate_w_dof_local:

• In-plane (CST 2 DOF/node $u_x,u_y$; mixed 2D+beam) — the w_dof_local < 0 branch: the strike is an in-plane force projected onto a contact direction $(n_x,n_y)$ using local DOFs 0 and 1, $f_i = \phi_i^{(u_x)}n_x + \phi_i^{(u_y)}n_y$. For an interior node the direction defaults to $\hat n=(0,1)$ (a vertical tap). For a boundary node (membership tested against mesh.boundary_nodes) the direction is the outward boundary normal: the CCW boundary-loop tangent $t=(x_{\text{next}}-x_{\text{prev}})$ is rotated 90° CW, $\hat n=(t_y,-t_x)/\lVert t\rVert$ (nx=ty/tl; ny=-tx/tl), so an edge strike pushes into the material along the surface normal.
• Plate / shell (out-of-plane $w$-DOF) — the w_dof_local >= 0 branch: a tap is purely transverse, so the force is just the single $w$ component, $f_i=\phi_i^{(w)}=\texttt{mode\_shapes[i][dofs\_per\_node*best + w\_dof\_local]}$. Here $w$ lives at local index 0 for the 3-DOF DKT plate $(w,\theta_x,\theta_y)$ and at index 2 for the 5-DOF Marguerre shell $(u_x,u_y,w,\theta_x,\theta_y)$; $(n_x,n_y)$ is irrelevant. This split matters: indexing membrane DOFs on a shell would dump tap energy into in-plane eigencontent and render as broadband noise (the explicit motivation in the source comment).
• 3D solid (3 DOF/node): there is no 2D mesh to project a point strike onto. solid_node_impulse (modal_api.cpp:3684) weights each mode by the magnitude of its motion at the struck surface node, $f_i=\lVert\phi_i(\text{node})\rVert=\sqrt{\phi_x^2+\phi_y^2+\phi_z^2}$, then normalizes by the max so the most-excited mode is $\approx 1$. An off-body click (or a fallback) excites all modes equally via modal_strike_uniform (3668), which builds a flat $f_i=1$.

A cross-component-leakage filter runs only when n_components > 1 (compute_strike_impulse, modal_api.cpp:3626): modes with $|f_i|<0.02\max_j|f_j|$ are zeroed, so striking one of several disjoint shapes does not make the others "sing along" through numerical eigenvector leakage. Single-component meshes (e.g. a string mechanically coupled to a box) skip the filter, because there every weakly-driven mode is a legitimate joint mode.

Each mode keeps a ring of $K=6$ recent strike slots $(A_k,t_{0,k})$ (STRIKE_RING_SIZE=6, struct ModeState{A[6],t0[6],head} in oscillator.h). A new strike writes the next round-robin head slot rather than overwriting, and the rendered displacement sums all slots:

Without the ring, a rapid second strike would reset $t_0$ and force the previous decaying sinusoid to restart from $\sin(0)=0$, producing an audible cancellation/click. project_strike (oscillator.cu:15) skips slots with $|f_i\cdot\text{magnitude}|<0.005$ (absolute threshold) and clamps each slot to $|A|\le 50$ so a ring slot cannot blow up.

Damping is per-mode, frequency-dependent. At solve time and on every modal_set_damping / _slope / _freq_enabled call the host rebuilds a per-mode vector $\zeta_i=\zeta_{\text{base}} + s\cdot\omega_i/(2\pi\cdot 1000\,\text{Hz})$, where the base $\zeta_{\text{base}}=\texttt{mat.zeta}$ is itself first clamped to $[0,0.999]$ (critical-damping ceiling) and each $\zeta_i$ is then clamped to $[10^{-5},0.95]$. The slope $s$ is gated by mat.zeta_freq_enabled (set to 0 when off, so every mode shares $\zeta_{\text{base}}$). High overtones therefore decay faster than the fundamental. The vector mode_zetas is uploaded to GPU d_zetas: the audio kernel audio_sum_kernel reads the per-mode $\zeta_i$. The visual decay path (VisualImpulseRing::eval) and the CPU reference eval_modes_cpu are, by contrast, called with a scalar $\zeta$, so the visual envelope uses a single damping ratio while the audio uses the full per-mode vector. The radiation coefficient $\alpha_i$ multiplies each mode's contribution to the PCM sum, $p=\sum_i\alpha_i q_i$, and a global user pitch knob freq_shift scales $\omega_i$ in the kernel.

An impulse cannot sustain a tone. The ResonatorBank (header-only) replaces the impulse with a continuously-driven filter bank: one RBJ constant-0 dB-peak bandpass biquad per mode, centered at $\omega_{0,i}=2\pi f_{c,i}/f_s$ (with $f_{c,i}$ clamped to $[20\,\text{Hz},\,0.45 f_s]$) and quality factor $Q_i=1/(2\zeta_i)\in[0.5,500]$, matched to the struck synth's damping so the timbre agrees. The normalized biquad ($b_1=0$) is

run in Direct-Form-I (state $x_1,x_2,y_1,y_2$ per mode). All modes share one source: an LCG white-noise stream ($\texttt{rng}=\texttt{rng}\cdot1664525+1013904223$) low-passed by a 1-pole filter whose cutoff fraction (color) controls timbre, $\text{lp}_a=e^{-\pi c}$. Each mode is driven by $x=\text{drive}\cdot w_i$ where $w_i$ is the per-mode contact weight set by set_sustain_weights. render_add adds to the output and tracks per-block energy; has_energy() keeps it ringing until mean-square energy $<10^{-6}$ after the drive stops.

The NonlinearVoice (header-only) replaces the linear struck sound for gong/cymbal blooms. Each mode is a free-ringing 2-pole resonator with poles at radius $r=e^{-\zeta_i\omega_{0,i}/f_s}$ and angle $\omega_{d,i}/f_s$, $y[n]=-a_{1,i}q_1-a_{2,i}q_2$, $a_{1,i}=-2r\cos(\omega_d/f_s)$, $a_{2,i}=r^2$, with state advanced $q_2\leftarrow q_1,\ q_1\leftarrow y$. A strike seeds the recurrence so it rings a decaying sine of amplitude $A_i=w_i\cdot\text{mag}$ (additively: $q_1{+}{=}A\,\text{seed}_1$, $q_2{+}{=}A\,\text{seed}_2$ with $\text{seed}_1=-\sin\omega_d/r$, $\text{seed}_2=-\sin 2\omega_d/r^2$). A per-mode amplitude envelope $\text{env}_i$ (peak-tracked, release ENV_REL=0.999) is the energy proxy. At control rate (every CTL=32 samples) two energy-conserving nonlinearities act on total energy $E=\sum_i\text{env}_i^2$ through a saturating factor $G=E/(E+E_{\text{ref}})$, $E_{\text{ref}}=0.25$:

• Pitch glide: effective frequency rises with energy, all coeffs rebuilt at $\omega_i^{\text{eff}}=\omega_i(1+g\,G)$ with $g=\text{amount}\cdot\text{glide}$, so a hard strike starts sharp and glides down as it decays.
• Shimmer cascade: with $\text{sh}=\text{amount}\cdot\text{shimmer}\cdot G$, a small fraction $\text{frac}=\min(0.04\,\text{sh},0.5)$ of each mode's energy moves to its octave-up neighbour ($\arg\min_j|\omega_j-2\omega_i|$, precomputed in target_), so energy climbs the spectrum (the cymbal "bloom"). Energy is conserved by scaling the source state by $\sqrt{1-\text{frac}}$ and injecting a ring of amplitude $\sqrt{\text{move}}$ ($\text{move}=\text{frac}\,\text{env}_i^2$) into the target's seeds.

All paths funnel through one primitive, modal_strike(ModalResult*, x, y, magnitude) in libmodal/src/modal_api.cpp (line 3642). It (a) pushes a StrikeEvent{x,y,magnitude} onto a lock-free StrikeBuffer drained by the audio thread, applies apply_sympathetic_coupling (§8) to the freshly projected force, and (b) immediately adds visual impulses ($f_i\cdot\text{magnitude}$ into each mode's VisualImpulseRing) so the mesh animates without audio latency. compute_strike_impulse (modal_api.cpp:3531) does the nearest-node search, builds the contact direction, dispatches DOF-aware projection, and applies the 2% cross-component leakage filter when n_components > 1. modal_compute_strike_amplitudes (modal_api.cpp:3719) is a pure query reusing the same projection (no side effects) — used by the viewer's strike-map/preview (§14). Solid variants: modal_strike_uniform (3668), solid_node_impulse (3684), modal_strike_solid_at (3701).

The audio thread drains events in modal_fill_audio (modal_api.cpp:3839): for each event it rebuilds f_impulse, calls project_strike under mode_states_mutex, uploads the $(A_{i,k},t_{0,k})$ arrays to GPU d_amps/d_onset, and runs the synth. The GPU kernel audio_sum_kernel in libmodal/src/audio/synthesis.cu (line 24) evaluates $p=\sum_i\alpha_i\sum_k A_{i,k}e^{-\zeta_i\omega_i\,dt}\sin(\omega_{d,i}\,dt)$ with per-mode zetas[i], a global freq_shift on $\omega$, and an isfinite NaN guard. The PortAudio callback pa_callback (audio_engine.cpp:139) fills the struck synth into a mono scratch (capped at 512 samples), then mixes the sustain bank (zeroing the struck buffer when sustain_mode_on), then the nonlinear voice (which zeroes and replaces the struck buffer while it has_energy()), then noise-gate, auto-gain, cavity SVF, and reverb to stereo.

Modal reduction of the semi-discrete structural ODE

The FEM assembly produces the second-order system

with $u\in\mathbb{R}^{n_{\rm dof}}$. The generalized symmetric eigenproblem

is solved on the GPU by cusolverDnSsygvdx with CUSOLVER_EIG_TYPE_1 (the $K\phi=\lambda M\phi$ form), CUSOLVER_EIG_MODE_VECTOR, CUBLAS_FILL_MODE_LOWER, and CUSOLVER_EIG_RANGE_I requesting the lowest indices $il=1,\ iu=\min(n,\ \texttt{num\_modes}+3)$ (a small over-fetch so rigid/near-zero modes can be discarded). The returned eigenvalues are $\lambda_i$; the code takes $\omega_i=\sqrt{\lambda_i}$ (rad/s, clamped to $0$ for $\lambda_i\le0$) and discards modes with $\omega_i<$ rigid_threshold. For a TYPE_1 problem sygvd/sygvdx returns eigenvectors normalized so $\Phi^{\mathsf T} M\,\Phi = I$, i.e. mass-orthonormal (the code does not re-normalize). Expanding $u(t)=\sum_i \phi_i\,q_i(t)$ and left-multiplying by $\phi_i^{\mathsf T}$ diagonalizes the system into $N$ scalar SDOF oscillators (the defining step of modal sound synthesis):

The damping matrix $C$ is never assembled; instead the code assumes modal (proportional) damping, assigning each oscillator a ratio $\zeta_i$ directly. This is the standard diagonal-damping assumption that makes $\Phi^{\mathsf T}C\Phi$ diagonal.

Damping model (affine / generalized-Rayleigh per-mode ζ)

Rather than the classic Rayleigh form $C=\alpha M+\beta K\Rightarrow \zeta_i=\tfrac12(\alpha/\omega_i+\beta\omega_i)$, the implementation uses an affine-in-frequency per-mode ratio:

with base ratio $\zeta_0$ (mat.zeta, default $0.008$) and slope $s$ (mat.zeta_freq_slope, default $5\times10^{-4}$). The slope is gated by the zeta_freq_enabled flag (default true); setting $s=0$ recovers uniform/constant-$\zeta$ damping. $s>0$ makes high overtones decay faster, reproducing the perceptually important "bright attack, mellow tail". The $\beta$-proportional limb of Rayleigh damping makes $\zeta\propto\omega$, so the slope term is its direct analog while $\zeta_0$ plays the role of a frequency-flat floor. Note the affine form uses ordinary frequency $f_i=\omega_i/2\pi$ in the denominator, so $s$ is a per-kHz increment, not a true $\beta$.

Impulse response and synthesis sum

For a strike at $t_0$ depositing modal amplitude $A_i$, the underdamped Green's function of each oscillator is the closed-form

evaluated directly (no time-stepping / no recurrence in the main audio path — a resonator IIR is used only in the auxiliary voices). The radiated pressure proxy fed to the DAC is the radiation-weighted modal sum

where $k$ indexes a small per-mode ring of recent strikes (STRIKE_RING_SIZE = 6) so overlapping strikes superpose additively instead of overwriting each other's phase origin (which would otherwise restart the sinusoid from $\sin 0=0$ and click). Slots with $A=0$ or $t\le t_{0}$ are skipped. An optional global pitch transposition multiplies every $\omega_i$ by freq_shift $=2^{\text{semitones}/12}$ (note: $\omega_d$ is computed from the shifted $\omega$ in the kernel).

Strike = impulse projected onto mode shapes

A point force of magnitude $g$ applied at world point $(x,y)$ projects onto mode $i$ via the modal participation factor. The code snaps the click to the nearest mesh node $j^\*$ and branches on the DOF layout:

• In-plane (CST 2D / mixed 2D+beam, w_dof_local<0): $f_i = \phi_i^{(u_x)}(j^\*)\,\hat n_x + \phi_i^{(u_y)}(j^\*)\,\hat n_y$, where $\hat n$ is the strike direction. For an interior node $\hat n=(0,1)$ ("tap from above"); for a boundary node the outward normal is built from a central-difference tangent across the boundary loop, $t=(x_{\rm next}-x_{\rm prev},\,y_{\rm next}-y_{\rm prev})$, and $\hat n=(t_y,-t_x)/\lVert t\rVert$ (CCW loop → rotate tangent 90° CW).
• Plate/shell (w_dof_local≥0): a vertical tap reads the transverse $w$-DOF directly, $f_i=\phi_i^{(w)}(j^\*)$, so the strike excites bending rather than membrane modes (the $(\hat n_x,\hat n_y)$ direction is unused).

The excitation amplitude is then $A_i = g\cdot f_i$. Because $\phi_i$ is mass-normalized, $f_i$ is exactly the modal force $\phi_i^{\mathsf T}f$ for a point impulse, i.e. the impulsive RHS that integrates to a jump in $\dot q_i$. A cross-component-leakage filter (zero out $|f_i|<0.02\max_k|f_k|$) is applied only when the mesh has $>1$ disjoint connected component, to suppress numerical cross-talk between localized eigenvectors; single-component meshes skip it so legitimate weakly-coupled joint modes survive. project_strike additionally discards $|A_i|<0.005$ and clamps $A_i\in[-50,50]$ before writing the ring slot.

Acoustic radiation coupling (discrete monopole boundary integral)

Each mode is weighted by a frequency-scaled normal-velocity surface integral (a lumped Rayleigh/monopole radiation efficiency):

where $A_j$ is the boundary-segment half-length (lumped "area" in 2D) accumulated at boundary node $j$, $\hat n_j$ is the segment-averaged, re-normalized outward normal, $\phi_i(j)=(\phi_i^{(u_x)}(j),\phi_i^{(u_y)}(j))$ (this routine assumes the 2-DOF in-plane layout, indexing $\phi[2j],\phi[2j{+}1]$), and the result is normalized so $\max_i|\alpha_i|=1$. The $\rho_{\rm med}\,\omega_i$ prefactor captures the leading $\omega$-dependence of monopole radiated power and the medium-density coupling ($\rho_{\rm med}$ default $1.21\ \mathrm{kg/m^3}$).

End-to-end audio path. viewer::AudioEngine::init (modal_viewer/src/audio/audio_engine.cpp:29) opens a PortAudio stream via Pa_OpenDefaultStream(..., 0 in, 2 out ch, paFloat32, 44100, 512 frames, pa_callback, this). The callback AudioEngine::pa_callback (line 139) calls modal::modal_fill_audio(result_, mono_scratch, n) into a MONO scratch buffer, runs the post chain, and writes interleaved stereo. If the host requests n>512 it fills only 512 and zero-pads the rest (line 153-156).

Per-mode oscillator state. modal::ModeState (libmodal/include/modal/oscillator.h:21) holds A[6], t0[6], round-robin head — the STRIKE_RING_SIZE=6 strike ring. ModalResult (defined in libmodal/src/modal_api.cpp) owns mode_states, frequencies (float, rad/s), mode_zetas (per-mode ζ), alphas (radiation, double), a StrikeBuffer, and a SynthesisState synth. Audio vs. render threads share mode_states under mode_states_mutex.

The synthesis kernel. audio_sum_kernel (libmodal/src/audio/synthesis.cu:24): one CUDA thread per output sample; each thread loops all n_modes and all ring_size strike slots, accumulating A*expf(-zeta_i*omega*dt_rel)*sinf(omega_d*dt_rel) with omega = freqs[i]*freq_shift, omega_d = omega*sqrtf(1-zeta_i^2), then p += alpha*mode_q; non-finite p is zeroed. synth_fill (line 112) launches block=256, copies device→pinned host→out_pcm, and advances current_time. It throws if n_samples is outside $[1,512]$ (device buffers are sized to 512).

Strike projection. compute_strike_impulse (modal_api.cpp:3531) does the nearest-node search, central-difference boundary-tangent/outward-normal computation, the DOF-aware dot product (plate/shell w-DOF branch at 3594, in-plane branch at 3601), and the multi-component-only 2% cross-leakage threshold (3626). modal_strike (line 3642) pushes a StrikeEvent{x,y,magnitude} to the lock-free StrikeBuffer (for audio) and immediately injects into the visual vis_impulses rings under lock. modal_compute_strike_amplitudes (line 3719) is the side-effect-free query reusing the same projection (used by the auto-tuner). project_strike (oscillator.cu:15) writes A_i = f_impulse[i]*magnitude into the next ring slot, skipping |A|<0.005 and clamping to ±50. Mesh-free / solid variants: modal_strike_uniform (3668, all modes equal) and modal_strike_solid_at (3701) with solid_node_impulse (3684, per-mode weight = normalized $\lVert\phi_i$ at surface node$\rVert$).

Strike drain & GPU upload. modal_fill_audio (modal_api.cpp:3839): drains the StrikeBuffer, per event rebuilds f_impulse (point via compute_strike_impulse+apply_sympathetic_coupling, or solid_node/uniform), project_strikes under lock, flattens mode_states into amps/onsets arrays of size nm*K, cudaMemcpyAsync to d_amps/d_onset, then synth_fill. modal_fill_audio_snapshot (line 3898) skips the drain/upload for read-only spectral analysis.

Solve-time setup. In modal_solve (modal_api.cpp lines ~3224-3334): radiation via compute_radiation_coefficients, float frequencies, per-mode mode_zetas from $\zeta_0+s\,\omega/(2\pi\cdot1000)$ clamped to $[10^{-5},0.95]$, create_synthesis_state(nm,44100), and cudaMemcpy uploads of d_freqs, d_zetas, d_alphas (alphas cast double→float). modal_set_damping (line 3769) clamps the base $\zeta\in[0,0.999]$, recomputes mode_zetas with the same formula, and re-uploads d_zetas live; modal_set_damping_slope (3799) and modal_set_damping_freq_enabled (3811) route through it. The generalized eigenproblem is solve_dense (solver/eigensolve_dense.cu:12, cusolverDnSsygvdx, range-I lowest indices, rigid_threshold filter).

Radiation. compute_radiation_coefficients (libmodal/src/modal/radiation.cpp:8) implements the boundary monopole integral exactly as written above (per-node accumulated segment half-lengths Aj, segment-averaged and re-normalized outward normals, rho*omega*sum, normalize to unit max).

Post-processing (audio_engine.cpp callback). In order: optional sustained noise voice add (if sustain_enable), nonlinear voice replace (if nl_enable), an 80 Hz-high-passed noise gate (raw - LP(raw), smoothstep gate on the audible-band envelope so sub-audible modes can't hold the gate open; full close resets state and jumps to the reverb stage), a 2-pole Chamberlin SVF cavity bandpass run in parallel and mixed back (out += gain*bp, cavity_freq/q/gain), then a post chain (enable_postproc) of master gain, 1-pole LP anti-aliasing, tanh soft-clip, and peak-tracking auto-normalize toward 0.9 gated by the gate state, and a final stereo reverb stage (Reverb::process; dual-mono passthrough when disabled). A lock-free SPSC ring (scope_l_/scope_r_) captures the final stereo output for the UI oscilloscope.

Auxiliary excitation voices. modal::ResonatorBank (include/modal/resonator_bank.h): one RBJ constant-0dB-peak bandpass biquad per mode (Direct-Form-I, center $f_i$, $Q=1/(2\zeta_i)$ clamped to $[0.5,500]$), driven by a shared LCG white-noise source LP-filtered by color — sustained bowed/rubbed tone; per-mode weights $\phi_i(\text{contact})$ are published lock-free from the UI via a double buffer. modal::NonlinearVoice (include/modal/nonlinear_voice.h): a free 2-pole state-space ring per mode (a1=-2 r cos(ω_d/f_s), a2=r^2, r=exp(-ζω/f_s)) seeded to a unit decaying sine, with control-rate (every CTL=32 samples) energy-dependent pitch glide ($\omega_{\rm eff}=\omega_0(1+\text{amount}\cdot\text{glide}\cdot G)$, $G=E/(E+E_{\rm ref})$, $E_{\rm ref}=0.25$) and an octave-up energy-cascade "shimmer" to the nearest mode at $\approx2\omega_i$.

Strike triggering. Interior/boundary clicks, drags, blow (mode 6), and bow (mode 7) flow through App::strike_active_ (app.cpp:6142); the wrapper fires modal_strike on the active result (or both A/B slots when compare is active) and optionally seeds the nonlinear voice. Drag-release strikes scale force by drag length (capped 5x); bow strikes scale force by mouse (bow) speed. Solids strike directly via modal_strike_uniform/modal_strike_solid_at (app.cpp:611-613), not through strike_active_. Sustain mode (interact_mode_==9) publishes per-mode weights via set_sustain_weights and toggles sustain_active (app.cpp:1795-1796, 2249).

This section is the back half of the pipeline: shape → mesh → FEM assemble → eigensolve (produces $\omega_i$, mass-orthonormal $\phi_i$) → [modal synthesis] → audio/viz. It consumes the EigenResult (frequencies + mode shapes) from solve_dense, augments it at solve time with per-mode damping ratios (mode_zetas) and radiation coefficients (alphas), and turns user strikes into excitation amplitudes via mode-shape projection. The CUDA audio_sum_kernel is the real-time renderer that converts the modal state into a 44.1 kHz PCM stream, which the PortAudio callback post-processes and outputs. The same mode_states/impulse model (via the parallel vis_impulses rings) also drives the animated displacement visualization (modal_get_displacements) and the spectral-analysis panels (modal_fill_audio_snapshot), so synthesis sits at the junction of audio output and visual feedback.

The structure obeys $K_s\phi_s=\omega_s^2 M_s\phi_s$ and the air cavity obeys the acoustic Helmholtz P1 eigenproblem $K_a\psi_c=\omega_c^2 M_a\psi_c$, where for a linear (P1) triangle the element stiffness and mass are

with $b_i,c_i$ the shape-function gradient coefficients, $A$ the triangle area, and $c$ the speed of sound. This is the FEM discretization of $\nabla^2 p + (\omega/c)^2 p = 0$: $K_a$ is the discrete Laplacian, $M_a$ the $1/c^2$-weighted consistent mass. The assembly (assemble_acoustic_KM) supports pressure-release ($p=0$) Dirichlet BCs at port nodes via symmetric row/column zeroing with a unit diagonal, but in the shipped v1 port_nodes is always empty (port localization is deferred), so every detected cavity is solved with the default all-Neumann (rigid-wall, sealed-box) boundary condition. The fully-Neumann acoustic operator has a constant-pressure null mode at $\omega\approx 0$; solve_cavity_modes discards it (and the unit-diagonal port penalty mode at exactly $\omega=1$ rad/s, were it present) with a kCavityRigidThreshold = 1.5 rad/s strict floor.

Reduced coupled eigenproblem

Rather than assembling a monolithic fluid-structure system, the code projects onto the two uncoupled modal bases (truncated to $N_s$ structural and $N_c$ cavity modes) and solves the small dense symmetric problem $A\,v=\lambda\,v$ with

and $\lambda=\omega^2$. Because $A$ is real symmetric by construction, solve_reduced_coupled uses Eigen's SelfAdjointEigenSolver; non-positive eigenvalues (from a pathological $L_r$) are dropped and the survivors re-sorted ascending. Each coupled eigenvector splits as $v=[a_s;\,a_c]$ giving the structural participation $a_s\in\mathbb{R}^{N_s}$ and cavity participation $a_c\in\mathbb{R}^{N_c}$ (returned in CoupledModes{omega,a_struct,a_cav}). Physical fields are reconstructed as $\phi^{\text{coupled}}=\sum_s a_s\phi_s$ and $p^{\text{coupled}}=\sum_c a_c\psi_c$. This is the symmetric (potential-form) reduction of the unsymmetric $(u,p)$ FSI problem, where $C=\int_\Gamma N_s\,(\mathbf n)\,N_a\,d\Gamma$ is the interface coupling operator and the modal projection $L_r=\Phi_s^\top C\,\Psi_c$ symmetrizes the reduced form.

Mode repulsion / avoided crossing

For a single matched pair $\omega_s=\omega_c=\omega$ with scalar coupling $g=L_r(0,0)$, the $2\times2$ block $\big[\begin{smallmatrix}\omega^2 & -g\\ -g & \omega^2\end{smallmatrix}\big]$ has eigenvalues $\omega^2\pm g$, i.e. coupled frequencies $\sqrt{\omega^2-g}$ and $\sqrt{\omega^2+g}$ — the lower branch pushed down, the upper pushed up. test_cavity_couple.cpp validates exactly this with $g=0.15\,\omega^2$, checking the split midpoint against the analytic $\tfrac12(\sqrt{\omega^2-g}+\sqrt{\omega^2+g})$.

Multi-chamber coupling

Each detected chamber $k$ contributes its own $L_{r,k}\in\mathbb{R}^{N_s\times N_{c,k}}$, horizontally concatenated into a single $L_r=[L_{r,1}\,|\,L_{r,2}\,|\dots]$ with cavity frequencies stacked in detection order. The off-diagonal blocks only connect the shared structural basis to each chamber — there is no direct cavity-to-cavity term — so two sealed chambers couple only through the structure, never directly. Consequently structural per-component re-localization is skipped whenever cavity coupling is enabled (the reduced $L_r$ is built around the whole-system structural basis; re-localizing would detune it — the documented DOUBLE_BOX regression where the cavity field went silent).

Sympathetic resonance operates at strike time in the modal amplitude domain, not in $K/M$. build_node_components flood-fills the structure into connected components over triangle-edge plus beam-edge node adjacency; build_sympathetic_pairs assigns each mode $k$ to the component holding $\ge 70\%$ (LOC_THRESHOLD=0.70) of its in-plane eigenvector energy $\sum_j(\phi_{x,j}^2+\phi_{y,j}^2)$ (else it is labeled delocalized, $-1$, and excluded). If there is $\le 1$ component the pair list is empty. For every pair $(i,j)$ in different components a Lorentzian weight is computed from the relative frequency mismatch:

where $Q$ (default 50, clamped to $\ge 1$) sharpens the match window and pairs with $w_{ij}<0.02$ (WEIGHT_FLOOR) are pruned. Note $Q$ enters the denominator linearly (it is not squared), multiplying the squared relative detuning. On a strike, the projected impulse vector is augmented symmetrically, $f_j \mathrel{+}= k\,w_{ij}\,f_i$ and $f_i \mathrel{+}= k\,w_{ij}\,f_j$, with $k$ = sympathetic strength (default 0.10), accumulated into a side buffer and added in a single pass to avoid double-coupling. Energy thus flows only between near-resonant modes in mechanically disjoint bodies. The pair list is built lazily on the first strike after enabling and rebuilt when the toggle, $Q$, or strength changes. Applied identically on the visual strike path (modal_strike) and the audio strike path (modal_fill_audio).

Each beam element is a 2-node Euler-Bernoulli frame element with 6 DOFs total (order $[u_1,v_1,\theta_1,u_2,v_2,\theta_2]$: axial $u$, transverse $v$, rotation $\theta$). Axial pre-tension $T$ adds a Hermite-consistent geometric stiffness on the transverse/rotational DOFs (the axial DOFs are unaffected); the local contribution added on top of the elastic bending block is

the standard $\tfrac{T}{30L}[36,3L,-36,3L;\dots]$ Hermite geometric matrix added to the elastic bending stiffness ($12EI/L^3$, $6EI/L^2$, $4EI/L$, $2EI/L$ entries). When $TL^2\gg EI$ the geometric term dominates and the eigenfrequencies approach the ideal taut-string law $f_n=\frac{n}{2L}\sqrt{T/\mu}$, $\mu=\rho A$, so pitch is set by tension rather than bending stiffness. At $T=0$ (default) it reverts to a pure-bending beam. Element matrices are formed in local axial coordinates and rotated to global via $K_e=R^\top K_l R$.

Multi-body / jointed structures. Beams are the only mechanical link between otherwise-disjoint CST strips (beam-boxes via add_box_edge, polygonal bells, string-on-bridge instruments via add_string_to_mesh). Before the eigensolve, obj->mesh.beam_edges is cleared and repopulated with one {n0,n1} per beam so the flood-fill treats beam-joined strips as one component. When cavity coupling is off and more than one component is found, solve_eigenmodes_per_component solves each component separately, keeping mechanically-disjoint objects from blending their near-degenerate modes into shared eigenvectors.

All three mechanisms live under libmodal and are exposed through the C-API in libmodal/include/modal/modal.h. Structure-cavity: the acoustic element/assembly is element_acoustic.cpp (P1 $K_e,M_e$; port pressure-release path dormant in v1); the reduced solve is cavity_couple.cpp (solve_reduced_coupled building $A$ directly and calling SelfAdjointEigenSolver); detection/meshing is cavity_detect.cpp / cavity_mesh.cpp (rasterize, exterior flood, single-dilate close, all boundary nodes become interface_nodes). Integration in modal_solve is gated on obj->cavity.enabled && method==Modal: it builds each Lr_k via build_reduced_coupling, concatenates blocks, runs solve_reduced_coupled, and rebuilds a merged EigenResult. The cross-chamber strike-leakage fix re-derives eff_dpn = mode_shapes[0].size()/n_struct_nodes so a 3-DOF mixed solve re-expressed in 2-DOF layout no longer mis-indexes the strike (validated by test_cavity_isolation.cpp's CAVITY_GLOCKENSPIEL). Sympathetic: the SympatheticState struct, build_node_components, build_sympathetic_pairs, and apply_sympathetic_coupling all live in modal_api.cpp, invoked on both strike paths; modal_set_sympathetic forces a rebuild on toggle/$Q$/strength change. Tension: beam_element_matrices adds the geometric block; beams are mirrored into mesh.beam_edges before the per-component eigensolve (solve_eigenmodes_per_component).

1. Reverb

1.1 Algorithmic (Freeverb / Schroeder-Moorer)

The algorithmic engine is a Freeverb: per output channel, $N_c=8$ parallel lowpass-feedback comb filters summed, then $N_a=4$ all-pass filters in series. The lowpass-feedback comb embeds a one-pole lowpass in the feedback path (the Moorer damping addition to Schroeder's reverberator). The implemented per-sample recursion, with delay length $D$, feedback $g$, and damping $d_1$ ($d_2=1-d_1$), is

so $s[n]=(1-d_1)\,y[n-D]+d_1\,s[n-1]$. The comb peaks are spaced at $f_s/D$ Hz; the lowpass makes high frequencies decay faster (frequency-dependent $T_{60}$). The Schroeder all-pass with delay $D_a$ and fixed gain $a=0.5$,

has unity magnitude but dense phase dispersion. The eight comb delays $\{1116,1188,1277,1356,1422,1491,1557,1617\}$ and four all-pass delays $\{556,441,341,225\}$ are the canonical Freeverb tuning at 44.1 kHz; the right channel adds a fixed spread=23 samples per delay. Room size maps to $g=\text{roomsize}\cdot0.28+0.7\in[0.7,0.98]$, damping to $d_1=\text{damp}\cdot0.4$; input is scaled by FIXED_GAIN=0.015.

1.2 Convolution (uniform-partitioned, overlap-save)

The convolution engine realizes $y = x * h$ for a stereo IR via uniform-partitioned FFT convolution with overlap-save. The IR of length $L$ is split into $P=\lceil L/B\rceil$ partitions of $B=512$ samples; each $h_p$ is zero-padded to $N=2B=1024$ and pre-transformed $H_p=\mathrm{FFT}_N(h_p)$. Each input block forms a length-$N$ frame (previous + current block), transformed to $X_m$ and pushed into a frequency-domain delay line. The output spectrum is the partition sum

and only the last $B$ samples of $y_m$ are valid (overlap-save discard of the first $B$). The unit tests verify it matches a direct convolution to $10^{-2}$ and that a unit-impulse IR is a passthrough. IR length is capped at 2 s ($\text{MAX\_IR}=44100\cdot2$).

1.3 Synthetic IR + top-level shaping

make_synthetic_ir builds a deterministic IR: an exponentially decaying ($\sim-60$ dB across the IR), low-passed, decorrelated LCG-noise tail plus early reflections at 13 ms (L) / 19 ms (R), RMS-normalized to $\approx0.25$; kinds Room 0.45 s / Hall 1.2 s / Plate 0.8 s (brighter LP) / Cathedral 2.0 s. The Reverb wrapper adds an integer pre-delay ring, a stereo width matrix ($\text{wet}_1=w\cdot0.5+0.5$, $\text{wet}_2=(1-w)\cdot0.5$ cross-mixing wet L/R), and a linear dry/wet crossfade $y=(1-m)x+m\,\tilde y$ against the un-predelayed dry input.

2. Nonlinear Voice (gong/cymbal)

Each mode is a free-decay 2-pole resonator — the homogeneous difference equation of a damped oscillator, with no numerator/per-sample input term (the mode rings from seeded history). Coefficients: $r_i=e^{-\zeta_i\omega_i/f_s}$, $\omega_{d,i}=\omega_i\sqrt{1-\zeta_i^2}/f_s$, $a_{1,i}=-2r_i\cos\omega_{d,i}$, $a_{2,i}=r_i^2$; recurrence $y[n]=-a_{1,i}y[n-1]-a_{2,i}y[n-2]$. A strike seeds $q_1\mathrel{+}=A\,\text{seed}_1$, $q_2\mathrel{+}=A\,\text{seed}_2$ with $\text{seed}_1=-\sin(\omega_d)/r$, $\text{seed}_2=-\sin(2\omega_d)/r^2$. A per-mode envelope $\text{env}_i$ is a peak/leaky-release follower of $|y|$. Two control-rate nonlinearities run every $\text{CTL}=32$ samples. With energy proxy $E=\sum_i\text{env}_i^2$ and saturating drive $G=E/(E+E_\text{ref})$, $E_\text{ref}=0.25$:

• Pitch glide: $\omega_i^{\text{eff}}=\omega_i(1+g\,G)$, $g=\text{amount}\cdot\text{glide}$; $a_1,a_2$ recomputed (trig per mode each tick). A hard strike starts sharp and glides down.
• Shimmer cascade: each mode $i$ transfers energy to the single octave mode $t=\arg\min_j|\omega_j-2\omega_i|$. With $\text{frac}=\min(0.5,\,0.04\cdot\text{sh})$, $\text{sh}=\text{amount}\cdot\text{shimmer}\cdot G$, the source $(q_1,q_2,\text{env})$ is scaled by $\sqrt{1-\text{frac}}$ and $\text{move}=\text{frac}\cdot\text{env}_i^2$ is injected into the target by seeding $\sqrt{\text{move}}\cdot(\text{seed}_1,\text{seed}_2)$, with $\text{env}_t\leftarrow\sqrt{\text{env}_t^2+\text{move}}$. The $\sqrt{\cdot}$ bookkeeping conserves amplitude-squared energy; energy climbs the spectrum ("bloom"). In the callback this voice zeroes and replaces the linear struck sound while it has energy.

3. Sympathetic Resonance (modal-domain Lorentzian coupling)

Cross-component sympathetic resonance applied in the modal domain at strike time, not as an audio filter. After the solve, the mesh is flood-filled into connected components (build_node_components, triangle + beam-edge adjacency). Each mode $k$ is assigned to the component with the majority of its eigenvector energy, $c(k)=\arg\max_c\sum_{j\in c}\|\phi_k(j)\|^2$, accepted only if that fraction exceeds $\text{LOC\_THRESHOLD}=0.70$ (else delocalized, excluded). The eigen-frequencies in eig.frequencies are angular ($\omega_i$ in rad/s, per the solver), so coupling is computed directly in $\omega$. For each pair $(i,j)$ in different components,

a Lorentzian in fractional detuning, sharpness $Q$ (default 50, clamped $\ge1$). Pairs below $w<0.02$ ($\text{WEIGHT\_FLOOR}$) are dropped. At strike the impulse $f$ is cross-fed symmetrically through a separate accumulator: $\text{add}_j\mathrel{+}=k\,w_{ij}\,f_i$, $\text{add}_i\mathrel{+}=k\,w_{ij}\,f_j$, then $f\mathrel{+}=\text{add}$, $k=\text{strength}$ (default 0.10). The same coupling is wired into both the audio and visual strike paths.

4. Sustained / Continuous Excitation (resonator bank)

Each mode is an RBJ constant-0 dB-peak bandpass biquad at $\omega_0=2\pi f_i/f_s$ with $Q_i=1/(2\zeta_i)$ clamped to $[0.5,500]$. With $\alpha=\sin\omega_0/(2Q)$, $a_0=1+\alpha$:

run Direct-Form-I: $y[n]=b_0x[n]+b_2x[n-2]-a_1y[n-1]-a_2y[n-2]$. All modes share one filtered-noise source: an LCG white generator $r\leftarrow1664525r+1013904223$ feeding a one-pole lowpass $\text{lp}\leftarrow\alpha_n\text{lp}+(1-\alpha_n)w$, $\alpha_n=e^{-\pi c}$ ($c$ = color). Per-mode input is $\text{drive}\cdot\phi_i(\text{contact})$ from a published weight vector. While driven it sustains; on release the high-Q bandpasses ring down.

5. Analysis Transforms for Visualization

The Sound-Analysis window is fed by the analytic modal state, not a PCM FFT. Each mode is a Lorentzian line of power $p_i=|A_i|^2$ (with $A_i$ the envelope-only amplitude $\sum A\,e^{-\zeta\omega\Delta t}$, deliberately omitting $\sin(\omega_d t)$ to avoid frame-rate aliasing), center $f_i=\omega_i/2\pi$, width $W_i=\max(W_{\min},\zeta_i f_i)$, $W_{\min}=3$ Hz:

Derived views: STFT = $S_\text{dB}$ on $N_f=512$ log-spaced bins $f_k=f_\min(f_\max/f_\min)^{k/(N_f-1)}$ (default $f_\min=1$ Hz, $f_\max=22$ kHz, user-adjustable), 512 columns; cochleagram = triangular mel filterbank ($m=2595\log_{10}(1+f/700)$) on the log-bin spectrum; chroma = fold each bin to $pc=(\mathrm{round}(12\log_2(f/440)+69))\bmod12$ and sum power; "CQT" = a chromatic re-permutation of a per-mode (pc, octave) piano-roll histogram over 9 octaves (modes binned by $\text{amp}^2$, rows re-laid as $\text{oct}\cdot12+pc$) — not a constant-Q filterbank on PCM. Spectral statistics are computed over the $\le64$ tracked-mode envelope powers (not the spectrogram bins): power-weighted centroid $\sum f p/\sum p$; flatness = geometric/arithmetic mean of the per-mode powers; power-weighted RMS inharmonicity $\sqrt{\sum_i p_i\text{cents}_i^2/\sum_i p_i}$ with $\text{cents}_i=1200\log_2(f_i/(n_i f_0))$, $f_0$ = lowest tracked mode at $\ge25\%$ of max amplitude, $\pm150$-cent outlier rejection. Separately a genuine radix-2 Cooley-Tukey FFT helper (fft_helper.h: Hann window, suggested $N=4096$) exists and is unit-tested but is not in the live panel: VizAudioFFT::compute_analytical_spectrum_ instead sums analytic Lorentzians on a linear axis ($f_k=k f_s/N$), a documented pivot away from PCM FFT.

All DSP paths are header-only, deterministic, and unit-tested; wiring lives in the viewer audio engine.

Reverb — modal_viewer/src/audio/reverb.h (namespace viewer): rvb_fft (iterative radix-2). AlgoReverb (NCOMB=8, NALLPASS=4, comb {1116..1617}, allpass {556,441,341,225}, spread=23, FIXED_GAIN=0.015; set_params → feedback_=roomsize*0.28+0.7, damp1_=damp*0.4, damp2_=1-damp1_; comb() does store=out*damp2_+store*damp1_; allpass() gain 0.5). ConvolutionReverb (B=512,N=1024,MAX_IR=44100*2; set_ir builds HL_/HR_; process keeps FDL fdl_, prev prev_, valid slice Y[B..2B)). make_synthetic_ir (er 13/19 ms, RMS 0.25, kinds 0-3). Reverb (predelay ring pd_, wet1/wet2 width matrix, dry from original in[i], mix_).

Nonlinear voice — libmodal/include/modal/nonlinear_voice.h (modal::NonlinearVoice): configure builds a1_/a2_, seeds, target_[i]=argmin|w0_j-2w0_i|. strike seeds q1_/q2_. render_add(out,n,gain) runs -a1*q1-a2*q2, calls control_update() every CTL=32. control_update: glide update_coeffs(i,w0_[i]*(1+g*G)) (g=amount_*glide_), shimmer frac=min(0.5,0.04*amount_*shimmer_*G), E_REF=0.25, ENV_REL=0.999.

Sympathetic resonance — libmodal/src/modal_api.cpp: SympatheticState (line 138, Q=50,strength=0.10). build_sympathetic_pairs (3419): build_node_components (3366) flood-fill; per-mode component via energy argmax over phi[2*j+0/1]^2 (hardcoded dof_per_node=2), LOC_THRESHOLD=0.70; w=1/(1+Q*dw*dw), dw=(wi-wj)/wavg over eig.frequencies (rad/s), WEIGHT_FLOOR=0.02. apply_sympathetic_coupling (3506): symmetric add[] cross-feed, then f+=add. Wired into modal_fill_audio (line 3863, point-strike branch only — not uniform/solid) and modal_strike (line 3655). C API modal_set_sympathetic (3818) clears built on Q/toggle change. Declared modal.h:177-198. Units confirmed by eigensolve_dense.h:8 (omega in rad/s) and eigensolve_sparse.cpp:133 (omega=sqrt(lambda)).

Resonator bank — libmodal/include/modal/resonator_bank.h (modal::ResonatorBank): configure builds RBJ bandpass b0_,b2_,a1_,a2_ (b1=0), Q=1/(2*zeta) clamp [0.5,500]. render_add(out,n,weights,color,gain): LCG + one-pole LP (lp_a=exp(-pi*color)) + per-mode DF-I biquad. Owned by AudioEngine::sustain_bank_; weights via lock-free sustain_w_[2] (set_sustain_weights,configure_sustain).

Audio chain — modal_viewer/src/audio/audio_engine.cpp::pa_callback (44.1 kHz, capped to 512 frames). Verified order: modal_fill_audio (mono) → sustain bank (additive; mutes struck when sustain_mode_on) → nonlinear voice (zeroes then replaces struck while ringing) → noise gate (HP view raw-LP(raw), 80 Hz cutoff, smoothstep state) → air-cavity Chamberlin SVF → post-process (1-pole LP, tanh soft-clip, peak auto-normalize toward 0.9 x gate) → Reverb::process mono→stereo → scope ring. Helpers: configure_sustain/configure_nl/load_ir/select_synthetic_ir (lines 58/104/84/97).

Analysis — modal_viewer/src/app/sound_analysis_compute.h (viewer::sound_analysis): compute_lorentzian_column, envelope_db, freq_to_pitch_class, compute_chroma_column, hz_to_mel/mel_to_hz, build_mel_filterbank, compute_cochlea_column, ring_rotate_linear; kMinVisWidth=3. modal_viewer/src/ui/viz_sound_analysis.cpp (VizSoundAnalysis; constants in .h: kNTime=512,kNFreq=512,kModesTracked=64,kPianoRollOctaves=9,fmin_=1,fmax_=22000): rebuild_freq_caches_ builds log_bin_centers_; write_spec_column_ (line 136; envelope amps via modal_get_mode_envelope_amplitudes, width max(3,zeta*f), power a*a); write_chroma_pr_column_ (382); cqt_buf_to_rgba_ (336); write_spectral_stats_column_ (251; centroid/flatness over tracked-mode powers, inharmonicity with 25%-max f0 and +/-150 cent reject). True FFT: modal_viewer/src/app/fft_helper.h. viz_audio_fft.cpp compute_analytical_spectrum_ (line 15) sums Lorentzians on linear axis f=k*fs/N (pivot comment lines 9-14).

These subsystems sit at the tail and periphery of shape → mesh → assemble → eigensolve → modes → strike/synthesis → audio/viz. Sympathetic coupling acts at the strike→synthesis boundary inside libmodal: it rewrites the per-mode impulse vector f_impulse right after projection (which already applied a 2% cross-component leakage threshold), before mode states are seeded/uploaded to the GPU kernel, on the point-strike branch only. The nonlinear voice and resonator bank are alternative/additive synthesis voices in the audio callback that consume the same modal parameters ($\omega_i,\zeta_i$, contact weights) but replace (nonlinear voice while ringing; struck muted in sustain mode) or supplement the linear struck sound. The reverb is the terminal stereo stage and feeds the scope ring. The analysis transforms are the viz branch: they tap the modal state (read-only getters), not the audio bus, for spectrograms, chroma, the (pc,octave) "CQT" roll, and spectral statistics.

Each boundary node of the structure is treated as an elementary radiator carrying a position $x_i$, an outward unit normal $\hat n_i$, and a normal surface velocity $v_{n,i}$ driven by that node's live modal displacement. By the Kirchhoff–Huygens principle, the field at a point $x$ is the retarded superposition of their wavelets:

Here $\omega$ is the dominant ringing mode's angular frequency (a single-frequency snapshot), $c=343\ \mathrm{m/s}$, the retarded phase $\omega(t-r_i/c)$ gives an outgoing wave, $A(r)=1/(0.5+3\sqrt{r})$ is a 2D geometric-spreading amplitude falloff, and the $\max(\cos\theta_i,0)$ obliquity factor makes each source radiate only into its outward half-space (a baffled, dipole-like directivity). Pixels inside the mesh are masked by a barycentric point-in-triangle test, so the field appears only in the surrounding "air".

The field is evaluated by a single GLSL compute shader (#version 450, $16\times16$ local size, one invocation per output texel of an rgba8 image). Boundary sources live in an std430 SSBO (binding 1, a packed BSource{x,y,vn,out_nx,out_ny}) and the mesh triangles in a second SSBO (binding 2) for the interior mask; a u_src_step stride decimates sources for performance on dense boundaries. Output color is blue for $p>0$ and red for $p<0$, with $\alpha\propto|p|\cdot\texttt{exposure}$ clamped to ≈0.7, composited as an overlay and advanced by u_time. Built in App::init_gl and dispatched per frame when the medium_viz_ toggle is on.

This overlay is distinct from the acoustic cavity solve. The P1 acoustic element (§3, §2) solves an interior Helmholtz eigenproblem for the trapped-air modes that genuinely color the sound of cavity instruments and feed the audio. The medium overlay is an exterior radiation picture for intuition only — it does not contribute to synthesis.

An experimental, off-by-default coupling (apply_medium_feedback) lets the radiated field push back on the structure, so a struck object can sympathetically drive a separate nearby one "through the air". It is a cross-component displacement coupling: node→component ids are flood-filled once, then each receiver accumulates a force from sources on other disjoint components with a 2D $1/\sqrt{r}$ radiation falloff. It is subsampled (≈ every fourth frame) and disabled in slow-motion, where the visual clock barely advances and the force would degenerate into a destabilizing near-DC push. It is the quiet, airborne sibling of the louder modal-domain sympathetic coupling in §9.

All strategies optimize the same scalar objective. With detected (or mode-tracked) frequency $f$ and chromatic target $f^*$, the error in cents is

Targets are equal-tempered. In the default (absolute) mode $f^* = f_{\text{base}}\cdot 2^{s/12}$ with $s\in\{0,2,4,5,7,9,11,12\}$ (a fixed major-scale degree table et_semi[]) and octave base $f_{\text{base}}=261.63\cdot 2^{(\text{oct}-4)}$ Hz; in "flexible" mode (tuning_flexible_) $f_{\text{base}}$ is instead a captured anchor frequency tuning_flexible_base_hz_, so targets become scale-degree multiples of a measured pitch rather than absolute C-octave Hz. An optional drift compensation divides $f^*$ by the playback shift factor $2^{\text{semi}/12}$.

1. Bisection family (Fine-W / Fine-H / Fine-Size / Fine-WH / Manual-W / Filter-W)

These exploit the thin-plate/membrane scaling intuition $f \propto 1/L^2$ (larger $\Rightarrow$ lower pitch). The search variable $x$ (=$W$, $H$, or a size ratio) lives in a multiplicative bracket $[x_0/\beta,\,x_0\beta]$ with user factor $\beta=$ search_bracket_factor_ (clamped to $[1.05,30]$). The bracket update is monotone bisection: if $f>f^*$ set $x_{\text{lo}}\leftarrow x$, else $x_{\text{hi}}\leftarrow x$, then $x \leftarrow (x_{\text{lo}}+x_{\text{hi}})/2$. The relative-width collapse test is $(x_{\text{hi}}-x_{\text{lo}})/(x_{\text{hi}}+x_{\text{lo}})<10^{-3}$. Fine-WH (mode 3) is alternating-axis coordinate descent: bisect $W$ to collapse, switch to $H$, repeat (up to 12 rounds, each re-opening a window of ratio $0.4$--$2.5$ for rounds $\le6$ then $0.7$--$1.4$, with $\pm5\%$ jitter). On a non-converged collapse of mode 1/2/10 the search jitter-retries up to 8 times around the best-so-far (jitter halving every 4 retries).

2. Hellmann-Feynman Newton acceleration (Fine-W / Manual-W only)

For modes 1 and 6 the midpoint is replaced (when enabled and a derivative cache is available) by a Newton step using an analytic eigenfrequency sensitivity. The exact eigenvalue-perturbation (Hellmann-Feynman) result for $K\phi=\lambda M\phi$ ($\lambda=\omega^2$) with $\phi^TM\phi=1$ is $d\lambda_k/dp = \phi_k^{T}(dK/dp - \lambda_k\,dM/dp)\phi_k$. The code does NOT assume that normalization; instead it finite-differences the FULL Rayleigh quotient $\lambda(\phi)=\phi^TK\phi/\phi^TM\phi$ at the held eigenvector $\phi_k$, perturbing only node $x$-coordinates by $(1+\varepsilon)$ (equivalent to a horizontal stretch):

with $\omega_k$ in rad/s. The Newton step is $W_{\text{next}} = W + (f^*-f)/(df/dW)$, clamped to the live bracket. Only the width derivative is computed.

3. Fine-Grad: Gauss-Newton / minimum-norm step in log space

Working in log-frequency, model $\ln f = \ln f_0 + a\,\delta\ln W + b\,\delta\ln H$ with axis elasticities $a=\partial\ln f/\partial\ln W$, $b=\partial\ln f/\partial\ln H$ recovered from two one-sided finite-difference probes (probe fraction $\delta=$ kGradProbeDelta$=0.02$):

The desired log change $g=\ln(f^*/f_0)$ gives an underdetermined system $a\,\delta\ln W+b\,\delta\ln H=g$; the minimum-norm (pseudo-inverse) solution — equivalently the steepest-descent direction of $\tfrac12|\ln f-\ln f^*|^2$ — is

A learning rate $\eta=$ kGradLearningRate$=0.8$ damps it, each axis is capped at $|\delta\ln\cdot|\le$ kGradMaxLogStep$=\ln 1.6$, and the multiplicative update $W\!\cdot e^{\eta a k}$, $H\!\cdot e^{\eta b k}$ is clamped to the bracket. When $a^2+b^2<10^{-9}$ (both axes inert) it falls back to thin-plate uniform scaling $x\leftarrow x\sqrt{f_0/f^*}$. For a separable power law $f=CW^pH^q$ finite differencing recovers $a=p,\,b=q$, so one undamped step would land on target.

4. Bayesian optimization (Fine-WH+ / Bayes-W)

Minimize the (optionally MAC-tracked) cost $y(W,H)=|c|$ with a GP surrogate. The prior is a zero-mean GP with separable RBF kernel

The posterior at $\mathbf{x}_*$ uses $\alpha=(K+\sigma_n^2 I)^{-1}\mathbf{y}$ via an in-place float Cholesky $K+\sigma_n^2I=LL^T$:

Sampling maximizes Expected Improvement over the incumbent $y_{\text{best}}$ (minimization form, with $\Phi$ via std::erfc):

Initial design is Latin-hypercube (6 seeds for 2D, 8 for 1D); thereafter EI is scanned on a $32\times32$ grid (BO_ACQ_GRID) and the argmax is sampled. Bayes-W (mode 8) locks $H$ ($\ell_y\!\to\!10^{-3}$, $H$ bracket collapsed to a point), giving an effective 1D search.

5. Spectro power-law fit

When the Spectro strategy is engaged on a bisection mode (1/2/3/10), the per-step midpoint is overridden by OLS on the tracked-mode trajectory $\{(x_i,f_i)\}$ of $\ln f = a + b\ln x$ (natural logs):

It needs $\ge2$ tracked samples, rejects predictions outside the bracket, and arms a forced-bisection step on overshoot (cents sign flip). The shelved tune_predict.h::fit_powerlaw_predict is the same fit with a recency window and an $n{=}1$ thin-plate seed ($b=-2$, i.e. $x\sqrt{f/f^*}$).

6. Reliability layer

Spectrum verdict (classify_spectrum) picks, among audible-and-excited partials (audible window $[20,20000]$ Hz; excited = amp $\ge$ tune_excite_rel_floor_ $\times$ max-strike-amp), the one closest to the previously tracked partial within a $1.5\times$ ratio window (else the loudest); classifies Converged / KeepGoing / NoAudibleNote. A plateau counter stops when $|c|$ fails to improve by $0.5$c over tune_plateau_window_ ($=12$) steps, and a wall-clock watchdog (tune_solve_timeout_sec_ $=15$ s) aborts a stalled solve.

All control flow lives in App::render() in modal_viewer/src/app/app.cpp, driven by integer state fine_mode_ (strategy) and fine_phase_ (1=arm/wait, 2=build mesh, 3=wait-for-sound, 4=consume-sample/decide). start_fine_tune_for_target(int mode) (app.cpp:7617) computes fine_target_hz_ from the diatonic table et_semi[]={0,2,4,5,7,9,11,12} (absolute C-octave base, OR tuning_flexible_base_hz_ anchor when tuning_flexible_ is set; app.cpp:7642-7658) and sets up each strategy's bracket/seeds.

Strategy dispatch (fine_mode_): 1=Fine-W, 2=Fine-Size, 3=Fine-WH (alternating bisection), 4=Fine-W+ (adaptive), 5=Fine-WB (multi-bracket), 6=Manual-W, 7=Fine-WH+ (2D BO), 8=Bayes-W (1D BO), 9=Filter-W (decay-windowed mode pick), 10=Fine-H, 12=Fine-Grad. Mode selection + Spectro engagement gate (only modes 1/10/2/3) at app.cpp:3120-3140.

Bisection core (app.cpp:3576-3853): bracket update if (freq > fine_target_hz_) fine_lo_ = cur_val; else fine_hi_ = cur_val; (app.cpp:3684), collapse test at :3687, axis-switch logic for Fine-WH at :3689-3736, jitter-retry on non-converged collapse at :3737-3792 (8 retries kMaxCollapseRetries, jitter halving every 4). The HF Newton step is at app.cpp:3794-3827, gated on use_hf_acceleration_, (fine_mode_==1 || fine_mode_==6), and a populated df_dw_cache_.

Hellmann-Feynman derivative: df_dw_cache_ is filled once per solve at app.cpp:1221-1232 by modal_compute_hf_width_derivs (libmodal/src/modal_api.cpp:4307), but ONLY when use_hf_acceleration_ and fine_mode_ is 1 or 6. It re-assembles K/M via assemble_KM, perturbs node-x by $(1+\varepsilon)$ ($\varepsilon=10^{-3}$), and finite-differences the Rayleigh quotient rayleigh_quotient_num (modal_api.cpp:4295) at the stored eigenvector; it returns 0 (caller bisects) for plate-bending (obj->plate_bending) or beam meshes (!obj->beams.empty()). State/docs in app.h around 900-960.

Fine-Grad (mode 12): tune_grad.h::grad_wh_step implements the log-space minimum-norm / steepest-descent step; the 3-solve probe sub-state-machine (grad_probe_ 0/1/2 = base/W-probe/H-probe) is at app.cpp:3488-3574; constants kGradProbeDelta=0.02, kGradLearningRate=0.8, kGradMaxLogStep=ln(1.6) (0.470003629f), kGradMaxOuterSteps=30 in app.h:612-618. Setup at app.cpp:7732-7742. Degenerate fallback to thin-plate scale $\sqrt{f_0/f^*}$ at app.cpp:3556-3560.

Bayesian optimization (modes 7/8): gaussian_process.h provides GaussianProcess (in-place float Cholesky refit(), posterior predict()) and free function expected_improvement (uses std::erfc for $\Phi$). Driver: bo_draw_seeds() (LHS, xorshift32-permuted; app.cpp:7435), bo_pick_next() (returns next seed while seeding, else EI argmax on BO_ACQ_GRID=32 grid; app.cpp:7469), bo_reset() (app.cpp:7359). Phase-4 loop at app.cpp:3391-3478 adds samples to bo_gp_, tracks bo_best_*, refits, picks next. Kernel hyperparameters set in start_fine_tune_for_target app.cpp:7723-7728: lengthscale ~20% (1D) / 25% (2D) of bracket, set_noise(1.0), set_amp(2500). Objective optionally MAC-tracked via mac_find_best_match/select_objective_freq (app.cpp:3406-3423).

Fine-W+ adaptive (mode 4): sample collection app.cpp:3275-3318, adaptive_pick_next() (app.cpp:8211) does initial even seeding (fine_init_seeds_), a periodic GAP-explore safety net (app.cpp:8233-8244), a monotonicity-bisection shortcut, inverse-quadratic (Lagrange-form parabola) interpolation around the best sample (app.cpp:8316-8347), then a priority scorer (BRACKET/NEAR_BEST/GAP, app.cpp:8353-8439), all mode-continuity-gated. adaptive_check_terminate() (app.cpp:8009) handles convergence/budget/plateau/density-exhaustion. Objective per sample set by finalise_sample (app.cpp:7904): loudest-mode cents by default, optional multi-mode penalty $|c|+5\cdot\text{mode\_index}$ (MODE_PENALTY_PER_INDEX), optional A-weighting via a_weight_gain.

Fine-WB (mode 5): wb_pick_next() (app.cpp:8451) sweeps, finds same-mode sign-change brackets, sorts by endpoint promise, bisects each (cap wb_max_bisects_per_bracket_); wb_check_terminate() at app.cpp:8562.

Spectro power-law: spectro_track_step_ (app.cpp:8073, greedy nearest-log-freq mode tracking with a 1-semitone = 0.0833 log2 lost-guard) and spectro_pick_axis_ (app.cpp:8138, OLS log-log fit + bracket-reject + overshoot force-bisect). Engaged only on bisection modes 1/2/3/10 (gate app.cpp:3130), overrides the midpoint at app.cpp:3836-3838.

Reliability layer (toggle tune_reliability_enabled_): tune_verdict.h::classify_spectrum invoked once per new solve at app.cpp:2987-3008 with tune_excite_rel_floor_; NoAudibleNote early-stop at app.cpp:3258-3270; plateau counter at app.cpp:3619-3626/3628-3629; watchdog at app.cpp:2927-2939. Chromatic note-name display (not snapping) via nearest_note_name in modal_viewer/src/ui/ui_manager.cpp:100. Unit tests: tests/test_tune_grad.cpp, tests/test_tune_verdict.cpp (no test_tune_predict.cpp in the live tree).

This is the inverse-design / outer-loop stage that sits on top of the forward pipeline (shape → mesh → assemble K,M → generalized eigensolve → modes $f_k, \phi_k$). Each tuning iteration proposes geometric parameters ($W,H$ or size), calls build_tuning_component() to remesh and re-solve, reads back the modal spectrum (frequencies via modal_get_frequencies, strike excitation via modal_compute_strike_amplitudes, envelope amps via modal_get_mode_envelope_amplitudes), reduces it to a single tracked partial and a cents-to-target cost, and uses that to pick the next parameters. It closes the loop between the eigensolve output and the design input, ultimately producing a shape whose audible struck partial matches a chromatic target before strike/synthesis. The Hellmann-Feynman width derivative is the one place the inner FEM operators ($K, M, \phi$) are reused directly (via re-assembly + Rayleigh quotient) to compute an analytic sensitivity rather than treating the solver as a black box.

Each physical regime starts from an elastodynamic PDE — 2D Navier elasticity for membranes, the Kirchhoff/Mindlin plate equations for bending, Marguerre shallow-shell theory for curved panels, 3D elasticity for solids. Multiplying by a test function and integrating by parts gives the weak form; discretizing the domain into triangles (or tetrahedra) and interpolating the field with low-order shape functions $N_i$ — barycentric/area coordinates $L_i$ for the membrane and solid, Batoz rotation interpolants for the DKT plate — turns each element into the small dense matrices $K_e=\int_{\Omega_e} B^\top D B\,dV$ and $M_e=\int_{\Omega_e}\rho\,N^\top N\,dV$ (§3). Because a node only shares elements with its mesh neighbours, scattering these into the global $K,M$ yields sparse, banded matrices: roughly $O(\text{elements})$ nonzeros, not $O(n^2)$.

The element order is a deliberate cost choice. Audible timbre is set by the lowest few-dozen eigenpairs, and for those, constant-strain triangles and lumped-mass plates are accurate enough (the free–free bar test lands within ~15% of the analytic axial frequency) while staying cheap to assemble and giving a tight sparsity pattern. The DOF count is $\text{nodes}\times\{2\ \text{membrane},\,3\ \text{plate/mixed},\,5\ \text{shell},\,3\ \text{solid}\}$, so a typical interactive mesh sits in the low hundreds to low thousands of DOFs.

The assembled second-order system $M\ddot q + C\dot q + K q = f$ couples every DOF. Solving the generalized eigenproblem $K\phi_i=\omega_i^2 M\phi_i$ once and changing to mass-orthonormal modal coordinates $q(t)=\sum_i \phi_i\,\eta_i(t)$ diagonalizes it into $m$ independent scalar ODEs

This is the whole performance trick. The expensive work — the diagonalization, $O(n^3)$ dense or $O(n\cdot\text{nnz}\cdot\text{iters})$ sparse — happens once per geometry; thereafter every audio sample is a handful of multiply-adds per mode, an $O(m)$ recurrence completely decoupled from mesh size. A 2000-DOF plate and a 200-DOF bar cost the same per sample once each is reduced to, say, 48 modes.

The lowest modes are extracted by solve_eigenmodes (§5), dispatching on DOF count. For $n\le 512$ it densifies $K,M$ to single precision and calls cuSOLVER's cusolverDnSsygvdx — Cholesky reduction $M=LL^\top$, $C=L^{-1}KL^{-\top}$, tridiagonalize-and-diagonalize, requesting only the lowest indices (a partial spectrum, not all $n$ eigenpairs). The float32 Cholesky is the precision cliff: on ill-conditioned fine meshes it throws, which is caught and routed to the sparse path. For $n>512$ (or on that fallback), a CPU shift-invert Lanczos runs with operator $\mathrm{OP}=(K-\sigma M)^{-1}M$, $\sigma=10^{-4}$: the spectral transform sends the smallest $\lambda$ to the largest, best-separated $\mu=1/(\lambda-\sigma)$ — the regime Lanczos converges fastest in — and the shift also regularizes the singular free–free $K$. One SimplicialLDLT factorization is reused across the iteration. Both paths return only the lowest 24–120 elastic modes after filtering rigid-body modes.

• Async solve worker. Re-solves run on a background SolveWorker thread; the UI never blocks and audio keeps ringing the previous modal bank until the new one is swapped in atomically on completion.
• Solve-scoped caches. Quantities that depend only on the modal result — strike-amplitude evaluation, the Strike-Map feature fields (§14), the tuning verdict — are recomputed once per result (keyed by the result pointer), never per frame.
• Pinned tuning resolution. The inverse-design loop re-solves many times, so its meshes are pinned to a low density (segs = 12) that keeps each re-solve well under a second; the display mesh can be finer independently.

The tuning search (§11) needs the sensitivity of a partial's frequency to a geometric parameter without paying for a full re-eigensolve at every probe. Two accelerators do the heavy lifting: a Rayleigh-quotient shape derivative — perturb a coordinate, re-assemble $K,M$, and read $\lambda=\phi^\top K\phi/\phi^\top M\phi$ at the stored eigenvector (a first-order, $O(\text{nnz})$ estimate that skips diagonalization, modal_compute_hf_width_derivs); and the Fine-Grad Gauss–Newton step in log-space plus the Gaussian-process Bayesian sampler, which reach the target in far fewer full solves than naive bisection.

Free-free longitudinal bar

A $1\,\mathrm{m}\times0.05\,\mathrm{m}$ steel bar ($E=200\,\mathrm{GPa}$, $\rho=7800\,\mathrm{kg/m^3}$, the MaterialParams defaults) is meshed at MeshResolution::MEDIUM, assembled with the 2-DOF plane-stress CST, and solved for 15 modes. The first free-free axial mode is $\omega_{\text{axial}} = \pi c/L$, $c=\sqrt{E/\rho}\approx 5064\,\mathrm{m/s}\Rightarrow\omega_{\text{axial}}\approx 15894\,\mathrm{rad/s}$. The test asserts Approx(omega_axial).epsilon(0.15) — within 15%, the stated CST accuracy. The wide tolerance reflects that constant-strain triangles overestimate stiffness; FEM frequencies converge to the exact value from above.

Closed and open rectangular acoustic cavities

Air is modelled by the scalar Helmholtz/wave operator with $K_e$ independent of $c$ and $M_e\propto 1/c^2$ (asserted directly: K2==K1, M2==M1/4 under $c{:}1\to2$). For a sealed (all-Neumann) box the analytic duct modes are $f_{mn}=\frac{c}{2}\sqrt{(m/L_x)^2+(n/L_y)^2}$, so $f_{10}=c/2L_x$. A $0.60\times0.40\,\mathrm{m}$ box ($c=340$, 48×32 grid) gives the lowest mode to epsilon(0.06) (6%, "FEM on a coarse structured mesh is slightly stiff"); a pressure-release (Dirichlet) duct gives the open-closed quarter-wave $f_1=c/4L$ to 8%. Both go through solve_cavity_modes with the $1.5\,\mathrm{rad/s}$ cavity floor.

Element-level identities

The acoustic CST on a unit right triangle is checked entry-by-entry ($K_{00}{=}1$, $K_{11}{=}K_{22}{=}0.5$, $K_{12}{=}0$), zero row sums (constant-pressure null space, margin(1e-12)), symmetric, and the consistent mass $M_e=\tfrac{A}{12c^2}[2,1]$ summing to the triangle area ($\sum M_{ij}=0.5$). The membrane CST uses the consistent mass $M_e=\tfrac{\rho t A}{12}[2I,I]$; test_fem_assembly.cpp asserts $K=K^T$ (Frobenius $<10^{-3}$), a positive mass diagonal, and $2N\times2N$ sizing — it does not test full positive-definiteness. The tet and Mindlin elements use lumped mass: test_solid.cpp asserts $\sum\mathrm{diag}(M_e)/3=\rho V$ (2%) and test_thick_plate.cpp asserts the three $w$-DOF lumped masses sum to $\rho h A$ ($10^{-5}$).

Energy localization, thick-plate, shell, solid, coupled cavity

• Per-component localization: two disjoint 2-mass chains share a degenerate mode at $\omega=\sqrt{2k/m}=\sqrt 8$; solve_eigenmodes_per_component must return localized eigenvectors ($\min(e_0,e_1)/\max(e_0,e_1)<10^{-6}$). The converse — a connected body must ring across its extent — is checked by test_tri_diag.cpp (BELL_TRI cross-body ratio $>0.02$) and the glockenspiel isolation test.
• Thick plate / shear: a free plate must have exactly 3 rigid modes with a clear gap ($w[3]>1$, $w[2]<10^{-2}w[3]$), proving the stabilized Mindlin element is hourglass-free; assemble_KM_plate reproduces direct DKT below $h/\mathrm{span}=0.10$ and direct Mindlin above to $10^{-6}$.
• Shell flat-regression: at $z\equiv0$ the 15-DOF shell is block-diagonal in membrane vs bending DOFs ($<10^{-12}$) and matches assemble_KM_plate's lowest 6 bending modes within 1%; Marguerre coupling scales linearly with $z$ amplitude.
• Solid: a free 3D cube solved with rigid_threshold=0 must show exactly 6 near-zero rigid modes (counted as $f<0.01\,f_{\max}$, since float32 leaks them to small-but-nonzero $\omega$) with the 7th genuinely elastic.
• Coupled-cavity avoided crossing: equal $\omega=500$ and coupling $c=L_r(0,0)=0.15\,\omega^2$ split symmetrically with analytic midpoint $\tfrac12(\sqrt{\omega^2-c}+\sqrt{\omega^2+c})$ to $10^{-4}$; zero coupling reproduces uncoupled frequencies to $10^{-6}$.
• Audio / DSP: eval_single_mode realizes $q(t)=A\,e^{-\zeta\omega t}\sin(\omega_d t)$, checked at $t=\pi/\omega_d$; envelope_db matches $20\log_{10}(|A|e^{-\zeta\omega t})$ to $0.01$-$0.1\,\mathrm{dB}$ with an exact $-80\,\mathrm{dB}$ floor; grad_wh_step recovers the elasticities of a separable power law $f=CW^pH^q$ to 0.5% in one step.

• Discretization order. All bulk elements are lowest order: CST, linear (4-node) tetrahedra, DKT/Mindlin triangles. Spatial convergence is $O(h^2)$ in displacement but only $O(h)$ in strain/stress, so frequencies converge slowly from above; this is why the bar tolerance is 15% at MEDIUM and the duct 6-8%. MeshResolution (default MEDIUM) is the primary convergence knob.
• Mass formulation. Consistent mass (membrane CST, acoustic) gives an upper bound on frequencies; lumped mass (tet, Mindlin) gives a lower bound and removes rotary-inertia coupling. Mixing the two across element families is a deliberate accuracy/robustness trade.
• Precision split. Systems with $\le 512$ DOF go to a float32 dense path (cuSOLVER Ssygvdx, RANGE_I, lower-triangular, $iu=\min(n,\,\mathrm{num\_modes}+3)$); larger systems use a double-precision sparse shift-invert Lanczos ($\sigma=10^{-4}$, SimplicialLDLT). The dense Cholesky inside Ssygvdx can fail on stiff FINE systems — solve_eigenmodes wraps solve_dense in try/catch and falls back. Float32 also smears the theoretically-zero rigid modes into $\omega\sim10^{-2}$, which is why solids filter rigid modes by count (exactly 6), not an absolute threshold.
• Rigid-body / penalty handling. The default structural rigid floor is $\omega<1.0\,\mathrm{rad/s}$ (strict < in both paths). The sparse solver regularizes any DOF with simultaneously zero stiffness and zero mass by adding $\mathrm{reg}=\max(\varepsilon\cdot10^{-10},10^{-10})$ to both $K$ and $M$, preventing a singular shift-invert factorization. Dirichlet/port BCs use the identity-row penalty, whose $\omega=1\,\mathrm{rad/s}$ artifact forces the cavity floor to $1.5\,\mathrm{rad/s}$.
• Heuristic constants. Several thresholds are tuned, not derived: thick-plate dispatch ratio $0.10$, dense/sparse split $512$, Lanczos shift $10^{-4}$, $\mathrm{nev}=\min(\mathrm{num\_modes}+5,n-1)$, $\mathrm{ncv}=\min(2\,\mathrm{nev}+1,n)$, cavity floor $1.5\,\mathrm{rad/s}$, audible band $[20,20000]\,\mathrm{Hz}$. The tests pin the behavioral consequences of these constants, not the constants themselves.

All tests live in tests/ and are registered in tests/CMakeLists.txt; the TEST_SOURCES list (lines 6-41) holds 36 test_*.cpp files, each an add_executable linked against modal_lib + Catch2::Catch2WithMain and discovered via catch_discover_tests. Viewer-side logic headers (tune_grad.h, tune_verdict.h, strike_map.h, reverb.h, etc.) are compiled standalone by adding modal_viewer/src to the relevant test target's include path — they are intentionally dependency-free. Solver dispatch is in eigensolve.cpp (DENSE_DOF_THRESHOLD = 512, try solve_dense catch → solve_sparse); the dense path is eigensolve_dense.cu (cusolverDnSsygvdx, RANGE_I, FILL_MODE_LOWER), the sparse path eigensolve_sparse.cpp (shift-invert Lanczos, zero-K/zero-M regularization). The rigid filter is the strict if (omega < rigid_threshold) { res.num_rigid++; continue; } in both paths.

The Strike Map turns the modal result into a spatial preview of "what does it sound like if I hit it here". It is a pure header, modal_viewer/src/app/strike_map.h (unit-tested in tests/test_strike_map.cpp), that converts the per-node strike response into four perceptual scalars. For a node with modal strike-amplitudes $a_k$ (obtained from modal_compute_strike_amplitudes, i.e. the mode-shape displacement at the struck point), audible-band frequencies $f_k$, and per-mode damping $\zeta$:

• Loudness $= \lVert a\rVert_2$ — the Euclidean norm of the excited modal amplitudes.
• Brightness $= \dfrac{\sum_k f_k\,|a_k|}{\sum_k |a_k|}$ — the amplitude-weighted spectral centroid.
• Pitch $= f$ at $\arg\max_k |a_k|$ — the frequency of the most strongly excited partial.
• Decay $= \dfrac{1}{\zeta\,\omega_\text{dom}}$ — the $1/e$ time constant of the dominant partial.

The four scalars are precomputed per mesh node once per solve (cached by result pointer, $O(\text{nodes}\times\text{modes})$) and drawn as a Gouraud-shaded colormap heatmap over the 2D mesh — a selectable layer (loudness / brightness / pitch / decay) with an opacity slider — plus a hover readout that reports the predicted note, frequency in Hz, decay, and loudness for the node under the cursor. It is 2D-only, lives on the Experimental tab, and costs nothing when the layer is off. Wiring: App::recompute_strike_map_ (cache rebuild) and App::draw_strike_map_overlay in modal_viewer/src/app/app.cpp.

An OFF-by-default reliability layer, gated by tune_reliability_enabled_, makes the auto-tune state machine terminate on physically meaningful conditions rather than grinding. Its core is a pure header, modal_viewer/src/app/tune_verdict.h (unit-tested), whose classify_spectrum examines the post-solve eigen-spectrum. Among modes that are both audible and excited, it tracks the partial nearest (within a $1.5\times$ ratio window) to the previous step's tracked frequency — a continuity criterion that resists mode-swaps — and otherwise falls back to the loudest partial. It emits a NoteClass $\in$ {Converged, KeepGoing, NoAudibleNote}.

When the layer is on, that tracked partial becomes the measured frequency the search uses; a plateau "no-progress" stop ends grinding bisection once the cents error stops improving, a wall-clock watchdog aborts a state machine stuck in a non-terminating state, and a NoAudibleNote verdict stops the search immediately. When the layer is off, the classic baseline behaviour is unchanged. As an aside, a perceived "slow tuning" episode once traced to a Debug build, not the algorithm — the FEM engine must be built in Release for the real-time and search workloads to perform.

Consolidated, de-duplicated bibliography merged from all sections. Textbooks cited by several sections appear once. Academic / conference papers (SIGGRAPH, SCA, IJNME, journals) are marked with a left rule and a highlighted venue pill; the confidence pill reflects how directly the cited work matches the shipped implementation.

1. O. C. Zienkiewicz & R. L. Taylor, The Finite Element Method, Vol. 1–2.
   Textbook · Butterworth-Heinemannhigh
   Foundational for the CST plane-stress triangle, plate-bending elements, isotropic $D$ matrices, and the Galerkin $K=\int B^\top D B$, $M=\int \rho N^\top N$ formulation used throughout libmodal.
2. K.-J. Bathe, Finite Element Procedures (Prentice-Hall, 1996; 2nd ed. 2014).
   Textbookhigh
   Consistent vs lumped mass, the generalized eigenproblem $K\phi=\omega^2M\phi$, $M$-orthonormalization, the Cholesky reduction-to-standard-form that cuSOLVER's sygvdx implements, and Mindlin-Reissner transverse-shear plates.
3. T. J. R. Hughes, The Finite Element Method: Linear Static and Dynamic FEA (Dover).
   Textbookhigh
   Tetrahedral isotropic elasticity (Lamé $D$ with $\lambda,\mu$), area/volume-coordinate shape functions, and shear-locking in Mindlin plates.
4. R. D. Cook, D. S. Malkus, M. E. Plesha & R. J. Witt, Concepts and Applications of Finite Element Analysis, 4th ed. (Wiley, 2002).
   Textbook · Wileyhigh
   Cited directly in element_shell.cpp ("Cook-Malkus-Plesha s13") for the Marguerre/von-Kármán membrane-strain terms; also the reference for the Euler-Bernoulli beam Hermite stiffness/mass, the consistent CST mass, and the 15-DOF shell's CST/DKT blocks.
5. J.-L. Batoz, K.-J. Bathe & L.-W. Ho, "A study of three-node triangular plate bending elements," Int. J. Numer. Methods Eng. 15 (1980) 1771–1812.
   IJNME 1980high
   Exact source of the DKT $H^x/H^y$ rotation-interpolation functions and the $P,q,r,t$ edge constants implemented in element_dkt.cpp; the 9-DOF thin-plate block of the shell element and the flat-plate solver it reduces to at $z=0$.
6. M. Lyly, R. Stenberg & T. Vihinen, "A stable bilinear element for the Reissner-Mindlin plate model," Comput. Methods Appl. Mech. Eng. 110 (1993) 343–357.
   CMAME 1993high
   Source of the shear-stabilization factor $h^2/(h^2+\alpha h_K^2)$ implemented in mindlin_element_matrices (the code comment names it Lyly-Stenberg-Vihinen).
7. K. Marguerre, "Zur Theorie der gekrümmten Platte großer Formänderung" (Theory of shallow shells), Proc. 5th Int. Congr. Appl. Mech. (1938) 93–101.
   ICAM 1938high
   The shallow-shell theory whose linearized bending-stretching coupling ($\varepsilon\mathrel{+}=(\partial z/\partial x)(\partial w/\partial x)$ etc.) is implemented as $K_{mb}/K_{bb,\text{extra}}$ in element_shell.cpp and fed by compute_z_at.
8. K.-U. Bletzinger, M. Bischoff & E. Ramm, "A unified approach for shear-locking-free triangular and rectangular shell finite elements," Computers & Structures 75 (2000) 321–334.
   Computers & Structures 2000medium
   The DSG (Discrete Shear Gap) method named in element_mindlin.h's doc comment as the intended locking-free approach; the shipped element is NOT DSG (it is the stabilized linear-shear Lyly-Stenberg-Vihinen element), so this is cited context, not the actual method.
9. K. J. Bathe & E. N. Dvorkin, "A four-node plate bending element based on Mindlin/Reissner plate theory and a mixed interpolation," IJNME 21 (1985) 367–383.
   IJNME 1985low
   Canonical MITC/assumed-shear approach to Mindlin-plate locking; background context only — the code uses neither MITC nor DSG.
10. G. H. Golub & C. F. Van Loan, Matrix Computations, 4th ed. (Johns Hopkins Univ. Press, 2013).
    Textbookhigh
    Generalized symmetric-definite eigenproblems, Cholesky-based reduction, conditioning, and Lanczos — covers both solver paths and the float32-conditioning discussion.
11. B. N. Parlett, The Symmetric Eigenvalue Problem, SIAM Classics ed. (1998).
    SIAMhigh
    Authoritative treatment of symmetric Lanczos, loss of orthogonality and reorthogonalization, Ritz values/vectors — grounds the full-reorthogonalization sweep and lucky-breakdown handling in the sparse solver.
12. C. Lanczos, "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators," J. Research Nat. Bur. Standards 45(4):255–282, 1950.
    J. Res. NBS 1950high
    The original Lanczos tridiagonalization underlying solve_sparse's recurrence (the shift-invert spectral transformation is conceptual, not a library call).
13. T. Ericsson & A. Ruhe, "The spectral transformation Lanczos method for the numerical solution of large sparse generalized symmetric eigenvalue problems," Math. Comp. 35(152):1251–1268, 1980.
    Math. Comp. 1980high
    The shift-invert (spectral-transform) Lanczos for generalized symmetric problems — the precise algorithm of solve_sparse with $\mathrm{OP}=(K-\sigma M)^{-1}M$ mapping smallest $\lambda$ to largest $\mu$.
14. E. Anderson et al., LAPACK Users' Guide, 3rd ed. (SIAM, 1999).
    SIAMhigh
    cusolverDnSsygvdx mirrors LAPACK's ?sygvd/?sygvx (TYPE_1, RANGE_I, FILL_LOWER, $M$-orthonormal eigenvectors) — the dense GPU path's algorithmic specification.
15. NVIDIA, "cuSOLVER Library" (CUDA Toolkit documentation), dense generalized symmetric-definite eigensolver routines (sygvdx).
    NVIDIA CUDA Toolkithigh
    cusolverDnSsygvdx in eigensolve_dense.cu is the GPU symmetric-definite generalized eigensolver (RANGE_I, FILL_MODE_LOWER) producing the lowest modes for systems up to 512 DOFs.
16. G. Guennebaud, B. Jacob, et al., "Eigen v3" (2010, eigen.tuxfamily.org).
    Open-source libraryhigh
    Sparse $K/M$ storage and triplet assembly, the SimplicialLDLT factorization used by the shift-invert Lanczos fallback, and the dense tridiagonal SelfAdjointEigenSolver inside it.
17. J. F. O'Brien, C. Shen & C. M. Gatchalian, "Synthesizing sounds from rigid-body simulations," ACM SIGGRAPH/Eurographics Symposium on Computer Animation (SCA), 2002, pp. 175–181.
    ACM SCA 2002high
    The canonical discretize-then-diagonalize modal sound synthesis from FEM eigenmodes ($K\phi=\omega^2M\phi$) excited by projected impulses — the pipeline modal_solve/modal_strike implement, and the model the strike projection and sympathetic coupling modify.
18. K. van den Doel & D. K. Pai, "The sounds of physical shapes," Presence: Teleoperators and Virtual Environments 7(4):382–395, 1998.
    Presence (MIT Press) 1998high
    Foundational formulation of modal sound as a bank of decaying sinusoids $q_i(t)=a_i e^{-d_i t}\sin(2\pi f_i t)$ excited by mode-shape-projected impulses — exactly the closed-form sum in audio_sum_kernel / oscillator.cu and the projection in project_strike.
19. K. van den Doel, P. G. Kry & D. K. Pai, "FoleyAutomatic: physically-based sound effects for interactive simulation and animation," ACM SIGGRAPH 2001, pp. 537–544.
    ACM SIGGRAPH 2001high
    Real-time modal resonator bank driven by contact/impulse forces with per-mode gains from the excitation location — mirrors the per-mode (frequency, damping, gain) triplet, the strike-amplitude projection, and the resonator forms in NonlinearVoice / ResonatorBank.
20. J.-M. Adrien, "The missing link: Modal synthesis," in Representations of Musical Signals (MIT Press, 1991), pp. 269–298.
    MIT Press (chapter) 1991high
    Canonical statement of modal synthesis: a vibrating object as a parallel bank of second-order resonators with modal frequencies, dampings, and shape-dependent input/output gains — the conceptual model under both the struck synth and the ResonatorBank voice.
21. R. R. Craig Jr. & A. J. Kurdila, Fundamentals of Structural Dynamics, 2nd ed. (Wiley, 2006).
    Textbook · Wileyhigh
    The generalized eigenproblem, mass-orthonormal modal basis, modal superposition $u=\sum\phi_i q_i$, and proportional/Rayleigh damping that justify assigning a diagonal per-mode $\zeta_i$ and ground the affine $\zeta_i=\zeta_0+s\,\omega_i$ model.
22. R. Bristow-Johnson, "Cookbook formulae for audio EQ biquad filter coefficients" (Audio EQ Cookbook), ~1999.
    comp.dsp / web standardhigh
    The exact RBJ constant-0 dB-peak bandpass biquad coefficients ($\alpha=\sin\omega_0/2Q$; $b_0=\alpha/a_0$, $b_2=-\alpha/a_0$, $a_1=-2\cos\omega_0/a_0$, $a_2=(1-\alpha)/a_0$) used verbatim in resonator_bank.h.
23. H. Chamberlin, Musical Applications of Microprocessors, 2nd ed. (Hayden, 1985).
    Hayden Books · textbookmedium
    The 2-pole Chamberlin state-variable filter ($f=2\sin(\pi f_c/f_s)$, $q_{\text{inv}}=1/Q$, lp/bp recurrence) used for the cavity/body bandpass in the audio callback.
24. N. Raghuvanshi & M. C. Lin, "Interactive sound synthesis for large scale environments," ACM I3D 2006, pp. 101–108.
    ACM I3D 2006low
    Real-time mode-bank synthesis under a fixed audio-buffer budget — relevant only to the general block-based constraint here; the cited mode-pruning/quality-scaling is NOT implemented (the synth sums every mode every sample).
25. M. R. Schroeder, "Natural Sounding Artificial Reverberation," J. Audio Eng. Soc. 10(3):219–223, 1962.
    J. Audio Eng. Soc. 1962high
    AlgoReverb is the Schroeder reverberator topology: parallel combs then series all-pass; the unity-magnitude Schroeder all-pass is implemented verbatim with gain 0.5.
26. J. A. Moorer, "About This Reverberation Business," Computer Music Journal 3(2):13–28, 1979.
    Computer Music Journal 1979high
    The lowpass-feedback comb (one-pole damping in the comb feedback loop) in AlgoReverb::comb is Moorer's frequency-dependent-decay extension, giving damp1_/damp2_.
27. "Jezar at Dreampoint", Freeverb (public-domain reverb), c. 2000; analyzed in J. O. Smith, Physical Audio Signal Processing (CCRMA).
    CCRMA / public domainhigh
    The exact comb delays $\{1116..1617\}$, all-pass $\{556,441,341,225\}$, spread 23, and fixed gain 0.015 are the canonical Freeverb constants reproduced directly.
28. W. G. Gardner, "Efficient Convolution Without Input-Output Delay," J. Audio Eng. Soc. 43(3):127–136, 1995.
    J. Audio Eng. Soc. 1995high
    ConvolutionReverb implements partitioned block FFT convolution with a frequency-domain delay line — Gardner's partitioned convolution, here in its uniform-partition form.
29. A. V. Oppenheim & R. W. Schafer, Discrete-Time Signal Processing (Prentice Hall).
    Textbookhigh
    The convolution uses overlap-save (discard first $B$ of each length-$N=2B$ IFFT block); rvb_fft and fft_helper.h are standard radix-2 FFTs.
30. J. W. Cooley & J. W. Tukey, "An Algorithm for the Machine Calculation of Complex Fourier Series," Math. Comp. 19:297–301, 1965.
    Math. Comp. 1965high
    Both rvb_fft and fft_radix2_ are radix-2 Cooley-Tukey FFTs with bit-reversal.
31. J. C. Brown, "Calculation of a constant Q spectral transform," J. Acoust. Soc. Am. 89(1):425–434, 1991.
    J. Acoust. Soc. Am. 1991low
    The "CQT" tab's octave-aligned rows are inspired by the constant-Q concept, but the implementation is a modal-domain (pc, octave) histogram re-permutation, NOT a constant-Q filterbank on PCM — the connection is conceptual only.
32. D. O'Shaughnessy, Speech Communication: Human and Machine (Addison-Wesley, 1987); mel concept orig. S. S. Stevens et al., 1937.
    Textbookmedium
    compute_cochlea_column / hz_to_mel use exactly the $2595\log_{10}(1+f/700)$ form with a triangular mel filterbank.
33. L. Cohen, Time-Frequency Analysis (Prentice Hall, 1995).
    Textbookmedium
    write_spectral_stats_column_ computes power-weighted centroid and a flatness ratio (geometric/arithmetic mean) — here over the tracked-mode powers rather than STFT bins.
34. N. Bonneel, G. Drettakis, N. Tsingos, I. Viaud-Delmon & D. James, "Fast Modal Sounds with Scalable Frequency-Domain Synthesis," ACM SIGGRAPH 2008.
    ACM SIGGRAPH 2008low
    Frequency-domain modal audio; loose conceptual backdrop for the analytic Lorentzian spectrum the viz uses (the code does not implement this method).
35. J. N. Chadwick, S. S. An & D. L. James, "Harmonic Shells: A Practical Nonlinear Sound Model for Near-Rigid Thin Shells," ACM SIGGRAPH Asia 2009.
    ACM SIGGRAPH Asia 2009medium
    Nonlinear thin-shell sound (inter-mode energy transfer → pitch glides and spectral bloom) is the phenomenon NonlinearVoice's control-rate glide + octave shimmer cascade approximates heuristically, not a derived model.
36. D. L. James, J. Barbič & D. K. Pai, "Precomputed Acoustic Transfer," ACM Transactions on Graphics (SIGGRAPH), 2006.
    ACM TOG / SIGGRAPH 2006low
    Conceptual reference for modal-vibration-to-radiated-sound. The implemented radiation coefficient is a single-scalar discrete monopole boundary integral and the medium shader is a heuristic visual cartoon — neither implements PAT's output-sensitive transfer; cited as the distant target for faithful radiation.
37. J. R. Shewchuk, "Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator," Applied Computational Geometry (WACG), LNCS 1148, Springer, 1996, pp. 203–222.
    WACG / LNCS 1148, 1996high
    Triangle v1.6 is vendored verbatim in extern/triangle and compiled into modal_lib in TRILIBRARY mode; it backs every triangulate*() call (PSLG/CDT, quality+area switches, hole-list carving).
38. J. Ruppert, "A Delaunay Refinement Algorithm for Quality 2-Dimensional Mesh Generation," J. Algorithms 18(3):548–585, 1995.
    J. Algorithms 1995high
    The quality refinement (minimum-angle bound + area cap, the q/a switches on the Glass/Delaunay path) Triangle performs is Ruppert's Delaunay refinement.
39. H. W. Kuhn, "Some combinatorial lemmas in topology," IBM J. Res. Develop. 4(5):518–524, 1960 (Freudenthal/Kuhn simplicial subdivision).
    IBM J. Res. Develop. 1960medium
    build_cube splits each voxel into 6 conforming tets via the Freudenthal/Kuhn main-diagonal subdivision (named "Freudenthal" in-code); the same idea underlies the prism 3-tet split for revolved solids.
40. J. Nocedal & S. J. Wright, Numerical Optimization, 2nd ed. (Springer, 2006).
    Springer · textbookhigh
    Fine-Grad's log-space minimum-norm / steepest-descent step (with learning-rate damping and per-axis caps) and the HF Newton step are textbook instances of Gauss-Newton / Newton's method.
41. R. T. Haftka & Z. Gürdal, Elements of Structural Optimization, 3rd ed. (Kluwer, 1992).
    Kluwer · textbookhigh
    modal_compute_hf_width_derivs implements first-order eigenvalue sensitivity ($d\lambda/dp=\phi^\top(dK/dp-\lambda\,dM/dp)\phi$, here via a finite-differenced Rayleigh quotient) — the standard structural-optimization sensitivity result.
42. R. P. Feynman, "Forces in Molecules," Phys. Rev. 1939; H. Hellmann, Einführung in die Quantenchemie, 1937.
    Phys. Rev. 1939high
    The code explicitly names the eigenvalue-derivative trick (hold the eigenvector fixed while differentiating the operator) the "Hellmann-Feynman" theorem; the implementation finite-differences the full Rayleigh quotient rather than the analytic $dK/dp,dM/dp$ forms.
43. C. E. Rasmussen & C. K. I. Williams, Gaussian Processes for Machine Learning (MIT Press, 2006).
    MIT Press · textbookhigh
    gaussian_process.h is a direct implementation of GP regression with a separable RBF kernel and Cholesky-based posterior mean/variance (Algorithm 2.1).
44. D. R. Jones, M. Schonlau & W. J. Welch, "Efficient Global Optimization of Expensive Black-Box Functions" (EGO), J. Global Optim. 1998.
    J. Global Optim. 1998high
    Fine-WH+/Bayes-W are EGO: a GP surrogate over an expensive (FEM-evaluated) objective with the closed-form Expected-Improvement acquisition implemented in expected_improvement().
45. B. Shahriari, K. Swersky, Z. Wang, R. P. Adams & N. de Freitas, "Taking the Human Out of the Loop: A Review of Bayesian Optimization," Proc. IEEE, 2016.
    Proc. IEEE 2016high
    Modern survey covering the GP-surrogate + EI-acquisition + space-filling-initial-design (LHS) pipeline that modes 7/8 implement.
46. M. D. McKay, R. J. Beckman & W. J. Conover, "A Comparison of Three Methods for Selecting Values of Input Variables…" (Latin Hypercube Sampling), Technometrics, 1979.
    Technometrics 1979high
    bo_draw_seeds() is stratified Latin-hypercube sampling with a random (xorshift32) permutation, used as the BO initial design.
47. R. P. Brent, Algorithms for Minimization without Derivatives (Prentice-Hall, 1973).
    Prentice-Hall · textbookmedium
    adaptive_pick_next()'s inverse-quadratic-interpolation step (Lagrange-form parabola through three (param, cents) points, with safeguards back to bisection) is the IQI root step in the spirit of Brent's method.
48. R. J. Allemang, "The Modal Assurance Criterion (MAC): Twenty Years of Use and Abuse," Sound and Vibration, 2003.
    Sound and Vibration 2003medium
    The optional BO objective tracks the target mode across remeshes by MAC correlation (mac_find_best_match), avoiding mode-swap when frequencies cross.
49. R. M. Lewis, V. Torczon & M. W. Trosset, "Direct search methods: then and now," J. Comput. Appl. Math., 2000.
    J. Comput. Appl. Math. 2000medium
    Survey-level context for the bisection / coordinate-descent / bracket-collapse families (derivative-free direct search over 1–2 design parameters with black-box FEM evaluations); the code uses elementary bisection/coordinate descent, not the named pattern-search algorithms.
50. L. Meirovitch, Fundamentals of Vibrations (McGraw-Hill, 2001).
    Textbook · McGraw-Hillhigh
    Canonical modal-superposition and per-mode (diagonal) modal-damping framework — the generalized eigenproblem $K\phi=\omega^2M\phi$, damped frequency $\omega_d=\omega\sqrt{1-\zeta^2}$, and the oscillator $q_i=A e^{-\zeta\omega t}\sin\omega_d t$ that eval_single_mode / audio_sum_kernel implement directly (§6–§7).
51. A. W. Leissa, Vibration of Plates, NASA SP-160 (1969).
    NASA monograph SP-160high
    Source of the thin-plate dimensionless eigenvalues and the $f\propto(h/L^2)\sqrt{D/\rho}$ bending-frequency scaling with $D=Eh^3/12(1-\nu^2)$, matching assemble_KM_plate (§4).
52. N. H. Fletcher & T. D. Rossing, The Physics of Musical Instruments, 2nd ed. (Springer, 1998).
    Textbook · Springerhigh
    The ideal-string law $f_n=(n/2L)\sqrt{T/\mu}$ recovered by the beam geometric-stiffness limit $TL^2\gg EI$ (the residual $EI$ giving stretched-overtone inharmonicity), and the physical basis for sympathetic resonance between frequency-matched modes (§4–§8).
53. G. P. Bazeley, Y. K. Cheung, B. M. Irons & O. C. Zienkiewicz, "Triangular elements in plate bending — conforming and non-conforming solutions," Proc. 1st Conf. Matrix Methods in Structural Mechanics, Wright-Patterson AFB, 1965.
    Matrix Methods (WPAFB) 1965high
    The BCIZ thin-plate bending triangle the code actually implements — element_dkt.cpp states explicitly it is a simplified BCIZ element, not the full Batoz DKT (§3).
54. J. N. Reddy, Mechanics of Laminated Composite Plates and Shells, 2nd ed. (CRC Press, 2004).
    Textbook · CRC Pressmedium
    Validates the orthotropic plane-stress $\mathbf D$ matrix and the Voigt grain-angle rotation (engineering symmetry $\nu_{yx}=\nu_{xy}E_y/E_x$) in cst_element_matrices_ortho (§4).
55. H. J.-P. Morand & R. Ohayon, Fluid-Structure Interaction: Applied Numerical Methods (Wiley, 1995).
    Textbook · Wileyhigh
    Canonical reference for structure-acoustic-cavity coupling, modal projection of the interface operator, and the symmetric vs unsymmetric reduced formulations matching the $A=[[\operatorname{diag},-L_r],[-L_r^\top,\operatorname{diag}]]$ block that solve_reduced_coupled uses (§8).
56. L. Komzsik, The Lanczos Method: Evolution and Application (SIAM, 2003).
    SIAM 2003medium
    Shift-invert Lanczos recurrence, $M$-orthonormalization, reorthogonalization, and the Ritz-value mapping $\mu\to\lambda=1/\mu+\sigma$ — the exact algorithm in eigensolve_sparse.cpp (§5).
57. S. Bilbao, Numerical Sound Synthesis: Finite Difference Schemes and Simulation in Musical Acoustics (Wiley, 2009).
    Textbook · Wileymedium
    Context for the nonlinear gong/cymbal behaviour (energy-dependent pitch glide and spectral energy cascade) that NonlinearVoice approximates at control rate; the rigorous treatment is the von-Kármán nonlinear plate/shell PDEs — NonlinearVoice is a phenomenological control-rate analogue, not a PDE solve (§6).
58. catchorg/Catch2, "A modern, C++-native test framework" (github.com/catchorg/Catch2).
    Open-source libraryhigh
    The test framework (TEST_CASE / REQUIRE / Approx with .epsilon/.margin) whose semantics every asserted tolerance in §13 relies on; the 36 test_*.cpp executables link Catch2::Catch2WithMain and register via catch_discover_tests.
59. Avoided level crossing / mode repulsion (general physics of two coupled resonators with near-equal frequencies).
    General physicslow
    Explains the $\sqrt{\omega^2\pm g}$ split validated by test_cavity_couple.cpp's avoided-crossing case: a matched structural/cavity pair with coupling $g$ repels into a lower and an upper branch (§8).