Add vonkarman-modal submission for plate-strike-nl-v1; viewer removes rigid motion
10 changed files+424−14
solutions/plate-strike-nl-v1/vonkarman-modal.jsonadded+89−0View file
@@ -0,0 +1,89 @@
1+{
2+ "schema_version": 1,
3+ "submission_id": "vonkarman-modal",
4+ "challenge_id": "plate-strike-nl-v1",
5+ "submitter": {
6+ "name": "Jeremy Magland",
7+ "email": "jmagland@flatironinstitute.org"
8+ },
9+ "date": "2026-09-17",
10+ "description": "Von Karman thin-plate reduction (Airy stress function, in-plane inertia neglected, damping on bending only) of the challenge model, solved by a modal method: linear free-plate modes and clamped Airy modes from Bogner-Fox-Schmit finite elements, the cubic coupling tensor from their integrals, a two-step integrator exact in the linear part, and the retarded Rayleigh integral with cubic fractional delays. The chain is ordered by mesh, number of modes, and time resolution. The truncation in modes limits how far the nonlinear cascade can reach in frequency.",
11+ "code_url": "https://vault1.magland.org/challenges/plate-solutions/blob/main/solvers/vonkarman-modal/vk_modal_plate.py",
12+ "viz_url": "https://vault1.magland.org/challenges/plate-solutions/site/",
13+ "solutions": [
14+ {
15+ "solution_id": "vk-n12-w50-a50-m1",
16+ "output": "vonkarman-modal/vk-n12-w50-a50-m1.wav",
17+ "method": "von Karman plate (Airy form), modal method on Bogner-Fox-Schmit modes, Gautschi-type integrator exact in the linear part, retarded Rayleigh integral with cubic fractional delays",
18+ "parameters": {
19+ "elements_per_side": 12,
20+ "transverse_modes": 50,
21+ "airy_modes": 50,
22+ "oversampling": 1,
23+ "dt_s": 2.0833333333333333e-05,
24+ "highest_transverse_mode_hz": 1261.4
25+ },
26+ "timing_seconds": 5.742,
27+ "more_accurate_than": [],
28+ "hardware": "Intel Xeon workstation, single process, numpy/scipy"
29+ },
30+ {
31+ "solution_id": "vk-n16-w100-a100-m1",
32+ "output": "vonkarman-modal/vk-n16-w100-a100-m1.wav",
33+ "method": "von Karman plate (Airy form), modal method on Bogner-Fox-Schmit modes, Gautschi-type integrator exact in the linear part, retarded Rayleigh integral with cubic fractional delays",
34+ "parameters": {
35+ "elements_per_side": 16,
36+ "transverse_modes": 100,
37+ "airy_modes": 100,
38+ "oversampling": 1,
39+ "dt_s": 2.0833333333333333e-05,
40+ "highest_transverse_mode_hz": 2875.7
41+ },
42+ "timing_seconds": 42.933,
43+ "more_accurate_than": [
44+ "vk-n12-w50-a50-m1"
45+ ],
46+ "hardware": "Intel Xeon workstation, single process, numpy/scipy"
47+ },
48+ {
49+ "solution_id": "vk-n24-w150-a150-m2",
50+ "output": "vonkarman-modal/vk-n24-w150-a150-m2.wav",
51+ "method": "von Karman plate (Airy form), modal method on Bogner-Fox-Schmit modes, Gautschi-type integrator exact in the linear part, retarded Rayleigh integral with cubic fractional delays",
52+ "parameters": {
53+ "elements_per_side": 24,
54+ "transverse_modes": 150,
55+ "airy_modes": 150,
56+ "oversampling": 2,
57+ "dt_s": 1.0416666666666666e-05,
58+ "highest_transverse_mode_hz": 4356.9
59+ },
60+ "timing_seconds": 277.268,
61+ "more_accurate_than": [
62+ "vk-n16-w100-a100-m1"
63+ ],
64+ "hardware": "Intel Xeon workstation, single process, numpy/scipy"
65+ },
66+ {
67+ "solution_id": "vk-n24-w200-a200-m2",
68+ "output": "vonkarman-modal/vk-n24-w200-a200-m2.wav",
69+ "method": "von Karman plate (Airy form), modal method on Bogner-Fox-Schmit modes, Gautschi-type integrator exact in the linear part, retarded Rayleigh integral with cubic fractional delays",
70+ "parameters": {
71+ "elements_per_side": 24,
72+ "transverse_modes": 200,
73+ "airy_modes": 200,
74+ "oversampling": 2,
75+ "dt_s": 1.0416666666666666e-05,
76+ "highest_transverse_mode_hz": 5980.6
77+ },
78+ "timing_seconds": 623.987,
79+ "fields": {
80+ "path": "vonkarman-modal/vk-n24-w200-a200-m2.h5",
81+ "format": "hdf5"
82+ },
83+ "more_accurate_than": [
84+ "vk-n24-w150-a150-m2"
85+ ],
86+ "hardware": "Intel Xeon workstation, single process, numpy/scipy"
87+ }
88+ ]
89+}
solutions/plate-strike-nl-v1/vonkarman-modal/vk-n12-w50-a50-m1.wavadded+0−0View file
Binary file not shown.
solutions/plate-strike-nl-v1/vonkarman-modal/vk-n16-w100-a100-m1.wavadded+0−0View file
Binary file not shown.
solutions/plate-strike-nl-v1/vonkarman-modal/vk-n24-w150-a150-m2.wavadded+0−0View file
Binary file not shown.
solutions/plate-strike-nl-v1/vonkarman-modal/vk-n24-w200-a200-m2.h5added+3−0View file
@@ -0,0 +1,3 @@
1+version https://git-lfs.github.com/spec/v1
2+oid sha256:f977db9c13f54a408eea9725fabffcab02ff7e8d76db946f916052d406759e9f
3+size 37295930
solutions/plate-strike-nl-v1/vonkarman-modal/vk-n24-w200-a200-m2.wavadded+0−0View file
Binary file not shown.
solvers/kirchhoff-bfs/kl_bfs_plate.pymodified+59−8View file
@@ -119,15 +119,12 @@ class BFSPlate:
119119 def node(i, j):
120120 return i * nn + j
121121
122- K = np.zeros((ndof, ndof))
123- M = np.zeros((ndof, ndof))
124- for ex in range(n):
125- for ey in range(n):
126- nodes = [node(ex, ey), node(ex + 1, ey), node(ex, ey + 1), node(ex + 1, ey + 1)]
127- dofs = np.array([4 * nd + d for nd in nodes for d in range(4)])
128- K[np.ix_(dofs, dofs)] += Ke
129- M[np.ix_(dofs, dofs)] += Me
130122 self.node = node
123+ self._elem_quad = (XI.ravel(), ETA.ravel(), WQ)
124+ self._shape = (N, Nxx, Nyy, Nxy)
125+ K = self.assemble(Ke)
126+ M = self.assemble(Me)
127+ self.K, self.M = K, M
131128
132129 evals, evecs = eigh(K, M) # mass-normalized eigenvectors
133130 omega2 = np.maximum(evals, 0.0)
@@ -181,7 +178,61 @@ class BFSPlate:
181178 self.B2 = (self.vT - self.s1 * self.qT) / (self.s2 - self.s1)
182179 self.A2 = self.qT - self.B2
183180
181+ # ------------------------------------------------------------ assembly helpers
182+ def element_dofs(self, ex: int, ey: int) -> np.ndarray:
183+ nodes = [self.node(ex, ey), self.node(ex + 1, ey), self.node(ex, ey + 1), self.node(ex + 1, ey + 1)]
184+ return np.array([4 * nd + d for nd in nodes for d in range(4)])
185+
186+ def assemble(self, Ae: np.ndarray) -> np.ndarray:
187+ """Scatter one 16 x 16 element matrix (uniform mesh) into a dense global matrix."""
188+ A = np.zeros((self.ndof, self.ndof))
189+ for ex in range(self.n):
190+ for ey in range(self.n):
191+ dofs = self.element_dofs(ex, ey)
192+ A[np.ix_(dofs, dofs)] += Ae
193+ return A
194+
195+ def element_matrix(self, Cmat: np.ndarray) -> np.ndarray:
196+ """Element 'stiffness' for the bilinear form int [w_xx, w_yy, w_xy] Cmat [v_xx, v_yy, v_xy]^T dA (no D factor)."""
197+ _, Nxx, Nyy, Nxy = self._shape
198+ WQ = self._elem_quad[2]
199+ Ae = np.zeros((16, 16))
200+ for q in range(len(WQ)):
201+ B2 = np.stack([Nxx[q], Nyy[q], Nxy[q]])
202+ Ae += WQ[q] * (B2.T @ Cmat @ B2)
203+ return Ae
204+
205+ def boundary_dofs(self) -> np.ndarray:
206+ n = self.n
207+ b = []
208+ for i in range(n + 1):
209+ for j in range(n + 1):
210+ if i in (0, n) or j in (0, n):
211+ b.extend(4 * self.node(i, j) + d for d in range(4))
212+ return np.array(sorted(set(b)))
213+
184214 # ------------------------------------------------------------ evaluation at points
215+ def _locate(self, x: np.ndarray, y: np.ndarray):
216+ n, a, b = self.n, self.a, self.b
217+ ex = np.clip(np.floor(x / a).astype(int), 0, n - 1)
218+ ey = np.clip(np.floor(y / b).astype(int), 0, n - 1)
219+ return ex, ey, x / a - ex, y / b - ey
220+
221+ def _scatter(self, vals: np.ndarray, ex: np.ndarray, ey: np.ndarray) -> np.ndarray:
222+ out = np.zeros((len(ex), self.ndof))
223+ rows = np.arange(len(ex))
224+ for c, (cx, cy) in enumerate([(0, 0), (1, 0), (0, 1), (1, 1)]):
225+ nd = self.node(ex + cx, ey + cy)
226+ for d in range(4):
227+ out[rows, 4 * nd + d] += vals[:, 4 * c + d]
228+ return out
229+
230+ def basis_derivs_at(self, x: np.ndarray, y: np.ndarray):
231+ """Global shape-function second derivatives (xx, yy, xy) at points, each (npts, ndof)."""
232+ ex, ey, xi, eta = self._locate(x, y)
233+ _, Nxx, Nyy, Nxy = shape_functions(xi, eta, self.a, self.b)
234+ return self._scatter(Nxx, ex, ey), self._scatter(Nyy, ex, ey), self._scatter(Nxy, ex, ey)
235+
185236 def basis_at(self, x: np.ndarray, y: np.ndarray) -> np.ndarray:
186237 """Global shape-function values at points, shape (npts, ndof), dense."""
187238 n, a, b = self.n, self.a, self.b
solvers/vonkarman-modal/vk_modal_plate.pyadded+234−0View file
@@ -0,0 +1,234 @@
1+"""Von Karman (Airy stress function) thin-plate solution of plate-strike-nl-v1 by a modal method.
2+
3+The transverse displacement is expanded in the linear free-plate modes Phi_k computed with the
4+Bogner-Fox-Schmit elements of the linear solver (../kirchhoff-bfs/kl_bfs_plate.py), and the Airy
5+stress function in the clamped biharmonic eigenfunctions Psi_p on the same mesh. With
6+mass-normalized Phi and L2-normalized Psi, the equations are
7+
8+ q_k'' + 2 alpha_k q_k' + omega_k^2 q_k = f_k psi(t/T_s) + sum_{j,p} H^p_{kj} q_j eta_p,
9+ eta_p = -(E h / (2 zeta_p^4)) sum_{i,j} H^p_{ij} q_i q_j, H^p_{ij} = int Psi_p L(Phi_i, Phi_j) dA,
10+
11+where L(f, g) = f_xx g_yy + f_yy g_xx - 2 f_xy g_xy (Touze et al.). Time stepping uses a two-step
12+Gautschi-type scheme at a step of 1/(fs m): the damped linear part is propagated exactly (the
13+recurrence coefficients are those of the exact damped oscillator), and the strike and the cubic
14+coupling force are applied with a static-exact weight, which is second-order accurate and free of
15+numerical dispersion in the linear limit. The Rayleigh integral is evaluated with
16+3 x 3 Gauss points per element; each point's retarded acceleration is a cubic Lagrange
17+interpolation of the modal accelerations at neighbouring time steps, and the interpolation weights
18+are folded into a per-step table, so the pressure at a step is a short sum over the recent modal
19+accelerations. Damping acts on the bending motion only (the in-plane state is quasi-static),
20+a further approximation of the challenge model.
21+
22+Usage: python vk_modal_plate.py <challenge.json> <out-dir> [--fields tag,...] [--force F] [tag ...]
23+ where a tag is n<mesh>-w<transverse modes>-a<airy modes>-m<oversampling>, e.g. n24-w150-a150-m2
24+"""
25+from __future__ import annotations
26+
27+import json
28+import sys
29+import time
30+from pathlib import Path
31+
32+import numpy as np
33+from numpy.polynomial import legendre as leg
34+from scipy.io import wavfile
35+from scipy.linalg import eigh
36+
37+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "kirchhoff-bfs"))
38+from kl_bfs_plate import BFSPlate, load_challenge # noqa: E402
39+
40+
41+class VonKarmanModal:
42+ def __init__(self, ch: dict, n: int, n_w: int, n_airy: int, force_override: float | None = None):
43+ self.ch = ch
44+ p = dict(ch["parameters"])
45+ if force_override is not None:
46+ p["strike_peak_force"] = force_override
47+ ch2 = dict(ch)
48+ ch2["parameters"] = p
49+ self.p = p
50+ self.plate = BFSPlate(ch2, n) # linear modes, modal forces, alpha, omega
51+ pl = self.plate
52+ # transverse modes: the n_w lowest (rigid ones included; they carry no membrane coupling)
53+ self.n_w = min(n_w, pl.n_modes)
54+ self.omega = pl.omega[: self.n_w]
55+ self.alpha = pl.alpha[: self.n_w]
56+ self.f = pl.f[: self.n_w]
57+ V = pl.V[:, : self.n_w]
58+ # Airy modes: clamped biharmonic eigenfunctions, L2-normalized: K_b psi = zeta^4 M_1 psi on interior dofs
59+ Kb = pl.assemble(pl.element_matrix(np.array([[1.0, 1.0, 0.0], [1.0, 1.0, 0.0], [0.0, 0.0, 0.0]])))
60+ M1 = pl.M / pl.m # unit-density mass matrix
61+ interior = np.setdiff1d(np.arange(pl.ndof), pl.boundary_dofs())
62+ z4, psi_int = eigh(Kb[np.ix_(interior, interior)], M1[np.ix_(interior, interior)])
63+ self.n_airy = min(n_airy, len(z4))
64+ self.zeta4 = z4[: self.n_airy]
65+ Psi = np.zeros((pl.ndof, self.n_airy))
66+ Psi[interior] = psi_int[:, : self.n_airy]
67+ # coupling tensor H[p, i, j] by 3 x 3 Gauss quadrature per element
68+ g, wg = leg.leggauss(3)
69+ s = (g + 1) / 2
70+ ws = wg / 2
71+ ex, ey = np.meshgrid(np.arange(n), np.arange(n), indexing="ij")
72+ XI, ETA = np.meshgrid(s, s, indexing="ij")
73+ xq = (ex.ravel()[:, None] + XI.ravel()[None, :]).ravel() * pl.a
74+ yq = (ey.ravel()[:, None] + ETA.ravel()[None, :]).ravel() * pl.b
75+ WQ = np.tile((ws[:, None] * ws[None, :]).ravel() * pl.a * pl.b, n * n)
76+ Bxx, Byy, Bxy = pl.basis_derivs_at(xq, yq)
77+ Pxx, Pyy, Pxy = Bxx @ V, Byy @ V, Bxy @ V # (Q, n_w)
78+ PsiQ = pl.basis_at(xq, yq) @ Psi # (Q, n_airy)
79+ H = np.zeros((self.n_airy, self.n_w, self.n_w))
80+ for pi in range(self.n_airy):
81+ wp = (WQ * PsiQ[:, pi])[:, None]
82+ H[pi] = (Pxx * wp).T @ Pyy + (Pyy * wp).T @ Pxx - 2.0 * (Pxy * wp).T @ Pxy
83+ self.H = H
84+ self.Hflat = H.reshape(self.n_airy * self.n_w, self.n_w)
85+ self.coef = -(p["E"] * p["h"]) / (2.0 * self.zeta4)
86+ # Rayleigh integral: quadrature weights, delays, modal shapes at the points
87+ r = np.sqrt((p["listener_x"] - xq) ** 2 + (p["listener_y"] - yq) ** 2 + p["listener_height"] ** 2)
88+ self.tau = r / p["c_air"]
89+ self.cq = p["rho_air"] / (2 * np.pi) * WQ / r
90+ self.PhiQ = pl.basis_at(xq, yq) @ V # (Q, n_w)
91+
92+ def nonlinear_force(self, q: np.ndarray) -> np.ndarray:
93+ v = (self.Hflat @ q).reshape(self.n_airy, self.n_w) # v[p, k] = sum_j H[p,k,j] q_j
94+ eta = self.coef * (v @ q) # eta_p
95+ return v.T @ eta # F_k = sum_p v[p,k] eta_p
96+
97+ def run(self, fs: int, nsamp: int, m: int, field_grid=None, field_rate: float | None = None):
98+ p = self.p
99+ dt = 1.0 / (fs * m)
100+ T, Om = p["strike_duration"], 2 * np.pi / p["strike_duration"]
101+ nsteps = nsamp * m
102+ w2, al, f = self.omega**2, self.alpha, self.f
103+ # exact two-step recurrence of the damped linear oscillator: q+ = a1 q + a2 q- + b F
104+ rigid = self.omega == 0.0
105+ disc = np.sqrt(al**2 - w2 + 0j)
106+ E1, E2 = np.exp((-al + disc) * dt), np.exp((-al - disc) * dt)
107+ a1 = (E1 + E2).real
108+ a2 = -(E1 * E2).real
109+ with np.errstate(divide="ignore", invalid="ignore"):
110+ b = np.where(rigid, 0.0, (1.0 - a1 - a2) / np.where(rigid, 1.0, w2))
111+ if rigid.any():
112+ e2 = np.exp(-2.0 * al[rigid] * dt)
113+ a1[rigid] = 1.0 + e2
114+ a2[rigid] = -e2
115+ b[rigid] = dt * (1.0 - e2) / (2.0 * al[rigid])
116+ # fractional delays: each quadrature point's retarded acceleration is a 4-point Lagrange
117+ # interpolation of the modal accelerations at steps d0-1 .. d0+2; the interpolation weights
118+ # are folded into a per-step weight table G[step, mode]
119+ tau_steps = self.tau / dt
120+ d0 = np.floor(tau_steps).astype(int)
121+ fr = tau_steps - d0
122+ lag = np.stack([
123+ -fr * (fr - 1.0) * (fr - 2.0) / 6.0, # weight of step d0-1 (interpolating at d0+fr)
124+ (fr + 1.0) * (fr - 1.0) * (fr - 2.0) / 2.0, # d0
125+ -(fr + 1.0) * fr * (fr - 2.0) / 2.0, # d0+1
126+ (fr + 1.0) * fr * (fr - 1.0) / 6.0, # d0+2
127+ ], axis=1)
128+ dmax = int(d0.max()) + 2
129+ G = np.zeros((dmax + 1, self.n_w))
130+ for j, off in enumerate([-1, 0, 1, 2]):
131+ dd = d0 + off
132+ for u in np.unique(dd):
133+ sel = dd == u
134+ G[u] += (self.cq[sel] * lag[sel, j]) @ self.PhiQ[sel]
135+ ring = np.zeros((dmax + 1, self.n_w)) # ring[n mod (dmax+1)] = qdd at step n
136+ out = np.zeros(nsamp)
137+ q_prev = np.zeros(self.n_w)
138+ q = np.zeros(self.n_w)
139+ frames = None
140+ if field_grid is not None:
141+ fx, fy, ft = field_grid
142+ XX, YY = np.meshgrid(fx, fy, indexing="xy")
143+ PhiG = self.plate.basis_at(XX.ravel(), YY.ravel()) @ self.plate.V[:, : self.n_w]
144+ frame_every = int(round(fs * m / field_rate))
145+ frames = []
146+ for nstep in range(nsteps + 1):
147+ t = nstep * dt
148+ psi = 0.5 * (1.0 - np.cos(Om * t)) if 0.0 <= t <= T else 0.0
149+ Fext = self.nonlinear_force(q) + f * psi # forces other than the linear restoring and damping terms
150+ q_next = a1 * q + a2 * q_prev + b * Fext
151+ # acceleration from the equation of motion (velocity by centered difference)
152+ qdd = Fext - w2 * q - 2.0 * al * (q_next - q_prev) / (2.0 * dt)
153+ ring[nstep % (dmax + 1)] = qdd
154+ if nstep % m == 0: # the ring buffer starts at zero, which is the state before the strike
155+ k = nstep // m
156+ if k < nsamp:
157+ idx = (nstep - np.arange(dmax + 1)) % (dmax + 1)
158+ out[k] = np.sum(G * ring[idx])
159+ if frames is not None and nstep % frame_every == 0 and len(frames) < len(ft):
160+ frames.append((PhiG @ q).astype(np.float32))
161+ q_prev, q = q, q_next
162+ fields = None
163+ if frames is not None:
164+ fields = np.stack(frames).reshape(len(frames), len(field_grid[1]), len(field_grid[0]))
165+ return out, fields
166+
167+
168+def parse_tag(tag: str):
169+ parts = dict((kv[0], int(kv[1:])) for kv in tag.split("-"))
170+ return parts["n"], parts["w"], parts["a"], parts["m"]
171+
172+
173+def main():
174+ ch = load_challenge(Path(sys.argv[1]))
175+ out_dir = Path(sys.argv[2])
176+ args = sys.argv[3:]
177+ fields_for = set()
178+ force = None
179+ if "--fields" in args:
180+ i = args.index("--fields")
181+ fields_for = set(args[i + 1].split(","))
182+ args = args[:i] + args[i + 2:]
183+ if "--force" in args:
184+ i = args.index("--force")
185+ force = float(args[i + 1])
186+ args = args[:i] + args[i + 2:]
187+ tags = args or ["n12-w50-a50-m1", "n16-w100-a100-m1", "n24-w150-a150-m2", "n24-w200-a200-m2"]
188+ out = ch["output"]
189+ fs, nsamp = out["sample_rate_hz"], out["num_samples"]
190+ meta = []
191+ for tag in tags:
192+ n, n_w, n_a, m = parse_tag(tag)
193+ t0 = time.perf_counter()
194+ model = VonKarmanModal(ch, n, n_w, n_a, force)
195+ t1 = time.perf_counter()
196+ grid = None
197+ if tag in fields_for:
198+ fg = ch["fields"]
199+ nt, ny, nx = fg["shape"]
200+ grid = (np.linspace(0, ch["parameters"]["Lx"], nx), np.linspace(0, ch["parameters"]["Ly"], ny), np.arange(nt) / fg["frame_rate_hz"])
201+ s, fields = model.run(fs, nsamp, m, grid, ch["fields"]["frame_rate_hz"] if grid else None)
202+ dt_total = time.perf_counter() - t0
203+ sid = f"vk-{tag}"
204+ out_dir.mkdir(parents=True, exist_ok=True)
205+ wavfile.write(out_dir / f"{sid}.wav", fs, (s / out["signal_scale"]).astype(np.float32))
206+ entry = {
207+ "solution_id": sid,
208+ "output": f"{out_dir.name}/{sid}.wav",
209+ "method": "von Karman plate (Airy form), modal method on Bogner-Fox-Schmit modes, Gautschi-type integrator exact in the linear part, retarded Rayleigh integral with cubic fractional delays",
210+ "parameters": {"elements_per_side": n, "transverse_modes": model.n_w, "airy_modes": model.n_airy, "oversampling": m, "dt_s": 1.0 / (fs * m), "highest_transverse_mode_hz": round(float(model.omega.max() / (2 * np.pi)), 1)},
211+ "timing_seconds": round(dt_total, 3),
212+ "peak_si": float(np.max(np.abs(s))),
213+ }
214+ if fields is not None:
215+ import h5py
216+
217+ fg = ch["fields"]
218+ with h5py.File(out_dir / f"{sid}.h5", "w") as fh:
219+ fh.attrs["challenge_id"] = ch["challenge_id"]
220+ fh.attrs["solution_id"] = sid
221+ fh.attrs["quantity"] = fg["quantity"]
222+ fh.attrs["units"] = fg["units"]
223+ fh.create_dataset("x", data=grid[0])
224+ fh.create_dataset("y", data=grid[1])
225+ fh.create_dataset("t", data=grid[2][: fields.shape[0]])
226+ fh.create_dataset("u", data=fields, compression="gzip", compression_opts=4, chunks=(8, fields.shape[1], fields.shape[2]))
227+ entry["fields"] = {"path": f"{out_dir.name}/{sid}.h5", "format": "hdf5"}
228+ meta.append(entry)
229+ print(f"{tag}: modes w={model.n_w} a={model.n_airy} setup={t1 - t0:.1f}s total={dt_total:.1f}s peak={np.max(np.abs(s)):.4e} Pa", file=sys.stderr)
230+ print(json.dumps(meta, indent=1))
231+
232+
233+if __name__ == "__main__":
234+ main()
web/index.htmlmodified+1−0View file
@@ -38,6 +38,7 @@ button { font: inherit; padding: 4px 10px; }
3838 <button id="playBtn">Play</button>
3939 <label>Speed <input id="speed" type="range" min="-3" max="0" step="0.1" value="-2"> <span id="speedLabel"></span></label>
4040 <label>Exaggeration <input id="exag" type="range" min="0" max="3" step="0.1" value="1.5"> <span id="exagLabel"></span></label>
41+ <label><input id="derigid" type="checkbox" checked> remove rigid motion</label>
4142 <label>Time <input id="timeSlider" type="range" min="0" max="1" step="1" value="0"> <span id="timeLabel" class="muted"></span></label>
4243 </div>
4344 <div id="view"></div>
web/viewer.jsmodified+38−6View file
@@ -47,6 +47,7 @@ const speedLabel = document.getElementById("speedLabel");
4747 const exag = document.getElementById("exag");
4848 const exagLabel = document.getElementById("exagLabel");
4949 const timeSlider = document.getElementById("timeSlider");
50+const derigid = document.getElementById("derigid");
5051 const timeLabel = document.getElementById("timeLabel");
5152
5253 const renderer = new THREE.WebGLRenderer({ antialias: true });
@@ -104,16 +105,20 @@ function buildMesh() {
104105
105106 function updateSurface() {
106107 if (!field || !mesh) return;
107- const { u, nx, ny, umax } = field;
108+ const { u, nx, ny, x, y, plane, xm, ym } = field;
109+ const rigidOff = derigid.checked;
110+ const umax = rigidOff ? field.emax : field.umax;
108111 const pos = mesh.geometry.attributes.position;
109112 const col = mesh.geometry.attributes.color;
110113 const scale = umax > 0 ? (0.03 * Math.pow(10, parseFloat(exag.value)) / umax) * 0.1 : 0; // exaggerated height
111114 const base = frame * nx * ny;
115+ const c0 = plane[3 * frame], gx = plane[3 * frame + 1], gy = plane[3 * frame + 2];
112116 // PlaneGeometry vertices run row by row from +y (top) to -y; row r corresponds to y index (ny-1-r)
113117 for (let r = 0; r < ny; r++) {
114118 const j = ny - 1 - r;
115119 for (let i = 0; i < nx; i++) {
116- const v = u[base + j * nx + i];
120+ let v = u[base + j * nx + i];
121+ if (rigidOff) v -= c0 + gx * (x[i] - xm) + gy * (y[j] - ym);
117122 const k = r * nx + i;
118123 const z = v * scale;
119124 pos.setZ(k, z);
@@ -153,16 +158,42 @@ async function loadField(path) {
153158 const t = Float64Array.from(file.get("t").value);
154159 const u = Float32Array.from(file.get("u").value);
155160 file.close();
156- let umax = 0;
157- for (let i = 0; i < u.length; i++) { const a = Math.abs(u[i]); if (a > umax) umax = a; }
158- field = { x, y, t, u, nx: x.length, ny: y.length, nt: t.length, umax };
161+ // best-fit plane per frame (the free plate's rigid translation and rotations), so that the
162+ // elastic deformation can be shown on its own
163+ const nx = x.length, ny = y.length, nt = t.length;
164+ const xm = x.reduce((a, b) => a + b, 0) / nx, ym = y.reduce((a, b) => a + b, 0) / ny;
165+ let sxx = 0, syy = 0;
166+ for (let i = 0; i < nx; i++) sxx += (x[i] - xm) ** 2;
167+ for (let j = 0; j < ny; j++) syy += (y[j] - ym) ** 2;
168+ sxx *= ny; syy *= nx;
169+ const plane = new Float64Array(nt * 3); // c, gx, gy per frame
170+ for (let k = 0; k < nt; k++) {
171+ let s0 = 0, sx = 0, sy = 0;
172+ const base = k * nx * ny;
173+ for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++) {
174+ const v = u[base + j * nx + i];
175+ s0 += v; sx += v * (x[i] - xm); sy += v * (y[j] - ym);
176+ }
177+ plane[3 * k] = s0 / (nx * ny); plane[3 * k + 1] = sx / sxx; plane[3 * k + 2] = sy / syy;
178+ }
179+ let umax = 0, emax = 0;
180+ for (let k = 0; k < nt; k++) {
181+ const base = k * nx * ny;
182+ for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++) {
183+ const v = u[base + j * nx + i];
184+ const e = v - plane[3 * k] - plane[3 * k + 1] * (x[i] - xm) - plane[3 * k + 2] * (y[j] - ym);
185+ if (Math.abs(v) > umax) umax = Math.abs(v);
186+ if (Math.abs(e) > emax) emax = Math.abs(e);
187+ }
188+ }
189+ field = { x, y, t, u, nx, ny, nt, umax, emax, plane, xm, ym };
159190 timeSlider.max = String(field.nt - 1);
160191 frame = 0;
161192 simTime = 0;
162193 buildMesh();
163194 updateSurface();
164195 render();
165- status.textContent = `${path}: ${field.nx} × ${field.ny} points × ${field.nt} frames, max |u| = ${(umax * 1e3).toFixed(3)} mm`;
196+ status.textContent = `${path}: ${field.nx} × ${field.ny} points × ${field.nt} frames, max |u| = ${(umax * 1e3).toFixed(2)} mm including rigid motion, ${(emax * 1e3).toFixed(3)} mm elastic`;
166197 } catch (e) {
167198 status.textContent = `Could not load ${path}: ${e}`;
168199 }
@@ -201,6 +232,7 @@ const updSpeed = () => { speedLabel.textContent = `${Math.pow(10, parseFloat(spe
201232 const updExag = () => { exagLabel.textContent = `×${Math.pow(10, parseFloat(exag.value)).toPrecision(2)}`; updateSurface(); render(); };
202233 speed.addEventListener("input", updSpeed);
203234 exag.addEventListener("input", updExag);
235+derigid.addEventListener("change", () => { updateSurface(); render(); });
204236 updSpeed();
205237 exagLabel.textContent = `×${Math.pow(10, parseFloat(exag.value)).toPrecision(2)}`;
206238 select.addEventListener("change", () => loadField(select.value));