Add kirchhoff-bfs submission for plate-strike-v1 (meshes 4 to 48), field export via LFS, and this repository's site with a 3D field viewer
17 changed files+927−0
.gitattributesadded+1−0View file
@@ -0,0 +1 @@
1+*.h5 filter=lfs diff=lfs merge=lfs -text
.gitignoreadded+4−0View file
@@ -0,0 +1,4 @@
1+site/
2+__pycache__/
3+challenge-repo/
4+.challenge-cache/
.mochi/workflows/results.ymladded+32−0View file
@@ -0,0 +1,32 @@
1+name: Validate and publish
2+on:
3+ push:
4+ workflow_dispatch:
5+concurrency:
6+ group: site
7+ cancel-in-progress: true
8+jobs:
9+ site:
10+ runs-on: ubuntu-latest
11+ steps:
12+ - uses: actions/checkout@v4
13+ - name: Fetch LFS objects (field exports), if git-lfs is available
14+ run: git lfs version && git lfs pull || echo "git-lfs unavailable or pull failed; field files may be pointers"
15+ - name: Get the challenge repository (specifications)
16+ run: git clone --depth 1 --branch main https://vault1.magland.org/challenges/plate challenge-repo
17+ - uses: actions/setup-python@v5
18+ with:
19+ python-version: '3.12'
20+ - run: pip install git+https://vault1.magland.org/challenges/challenge-tools
21+ - name: Validate submissions
22+ run: challenge-tools validate --challenges challenge-repo/challenges --solutions solutions
23+ - name: Build this repository's site
24+ run: python web/build.py --out site
25+ - if: github.ref == 'refs/heads/main'
26+ uses: actions/configure-pages@v4
27+ - if: github.ref == 'refs/heads/main'
28+ uses: actions/upload-pages-artifact@v3
29+ with:
30+ path: site
31+ - if: github.ref == 'refs/heads/main'
32+ uses: actions/deploy-pages@v4
README.mdadded+35−0View file
@@ -0,0 +1,35 @@
1+# Solutions to the plate challenges
2+
3+This repository holds solutions to the challenges defined in [challenges/plate](https://vault1.magland.org/challenges/plate). It is one of possibly several solution repositories registered with that challenge repository, which clones each of them, compares every solution against a common reference, and publishes the comparison at [challenges/plate/site](https://vault1.magland.org/challenges/plate/site/). This repository's own site, at [challenges/plate-solutions/site](https://vault1.magland.org/challenges/plate-solutions/site/), lists its submissions with audio players and animates the exported displacement fields as a 3D surface.
4+
5+## Layout
6+
7+```
8+solutions/<challenge-id>/<submission-id>.json one submission: metadata and one or more solutions
9+solutions/<challenge-id>/<submission-id>/*.wav the scored output files, referenced from the JSON
10+solutions/<challenge-id>/<submission-id>/*.h5 optional displacement fields for visualization (git LFS)
11+solvers/ source code of the submissions
12+web/ this repository's site (built by web/build.py)
13+```
14+
15+Field exports are about 45 MB each and are stored with git LFS (`.gitattributes` routes `*.h5` there). Cloning without git-lfs leaves pointer files in their place; the validator tolerates that and the site marks such fields as unavailable.
16+
17+## The example submission
18+
19+`kirchhoff-bfs` solves the Kirchhoff-Love thin-plate reduction of the challenge model with Bogner-Fox-Schmit C1 finite elements on uniform meshes of 4 to 48 elements per side, a closed-form modal time response, and an exact evaluation of the retarded Rayleigh integral. It is a thin-plate approximation, so its converged limit differs from the three-dimensional model that defines the challenge; the chain is ordered by mesh refinement.
20+
21+```
22+python solvers/kirchhoff-bfs/kl_bfs_plate.py <plate-repo>/challenges/plate-strike-v1/challenge.json solutions/plate-strike-v1/kirchhoff-bfs --fields 48 4 8 12 16 24 32 48
23+```
24+
25+## Adding a submission
26+
27+Follow the same steps as in [challenges/string-solutions](https://vault1.magland.org/challenges/string-solutions#adding-a-submission-to-this-repository): produce the output exactly as the [specification](https://vault1.magland.org/challenges/plate/blob/main/challenges/plate-strike-v1/README.md) says, write `solutions/<challenge-id>/<submission-id>.json` following the [schema](https://vault1.magland.org/challenges/challenge-tools/blob/main/schema.json), record known accuracy orderings in `more_accurate_than`, and open a pull request.
28+
29+## Validating locally
30+
31+```
32+pip install git+https://vault1.magland.org/challenges/challenge-tools
33+challenge-tools validate --challenges <path-to-plate-repo>/challenges --solutions solutions
34+python web/build.py --out site
35+```
solutions/plate-strike-v1/kirchhoff-bfs.jsonadded+150−0View file
@@ -0,0 +1,150 @@
1+{
2+ "schema_version": 1,
3+ "submission_id": "kirchhoff-bfs",
4+ "challenge_id": "plate-strike-v1",
5+ "submitter": {
6+ "name": "Jeremy Magland",
7+ "email": "jmagland@flatironinstitute.org"
8+ },
9+ "date": "2026-09-17",
10+ "description": "Kirchhoff-Love thin-plate reduction of the challenge model, solved with Bogner-Fox-Schmit C1 finite elements on uniform meshes, a closed-form modal time response, and an exact evaluation of the retarded Rayleigh integral. This is a thin-plate approximation, so its converged limit differs from the three-dimensional model that defines the challenge. Ordered by mesh refinement.",
11+ "code_url": "https://vault1.magland.org/challenges/plate-solutions/blob/main/solvers/kirchhoff-bfs/kl_bfs_plate.py",
12+ "viz_url": "https://vault1.magland.org/challenges/plate-solutions/site/",
13+ "solutions": [
14+ {
15+ "solution_id": "bfs-n4",
16+ "output": "kirchhoff-bfs/bfs-n4.wav",
17+ "method": "Kirchhoff-Love plate, Bogner-Fox-Schmit C1 finite elements, closed-form modal time response, exact retarded Rayleigh integral",
18+ "parameters": {
19+ "elements_per_side": 4,
20+ "element_size_m": 0.0762,
21+ "n_dof": 100,
22+ "n_modes": 100,
23+ "rigid_modes": 3,
24+ "mode_cutoff_hz": 200000.0
25+ },
26+ "timing_seconds": 0.861,
27+ "more_accurate_than": [],
28+ "hardware": "Intel Xeon workstation, single process, numpy/scipy (LAPACK eigh)",
29+ "notes": "lowest nonrigid mode frequencies (Hz): 36.345, 53.051, 67.284, 94.598, 94.598, 168.854"
30+ },
31+ {
32+ "solution_id": "bfs-n8",
33+ "output": "kirchhoff-bfs/bfs-n8.wav",
34+ "method": "Kirchhoff-Love plate, Bogner-Fox-Schmit C1 finite elements, closed-form modal time response, exact retarded Rayleigh integral",
35+ "parameters": {
36+ "elements_per_side": 8,
37+ "element_size_m": 0.0381,
38+ "n_dof": 324,
39+ "n_modes": 324,
40+ "rigid_modes": 3,
41+ "mode_cutoff_hz": 200000.0
42+ },
43+ "timing_seconds": 3.321,
44+ "more_accurate_than": [
45+ "bfs-n4"
46+ ],
47+ "hardware": "Intel Xeon workstation, single process, numpy/scipy (LAPACK eigh)",
48+ "notes": "lowest nonrigid mode frequencies (Hz): 36.322, 53.008, 67.19, 94.383, 94.383, 167.852"
49+ },
50+ {
51+ "solution_id": "bfs-n12",
52+ "output": "kirchhoff-bfs/bfs-n12.wav",
53+ "method": "Kirchhoff-Love plate, Bogner-Fox-Schmit C1 finite elements, closed-form modal time response, exact retarded Rayleigh integral",
54+ "parameters": {
55+ "elements_per_side": 12,
56+ "element_size_m": 0.0254,
57+ "n_dof": 676,
58+ "n_modes": 676,
59+ "rigid_modes": 3,
60+ "mode_cutoff_hz": 200000.0
61+ },
62+ "timing_seconds": 8.639,
63+ "more_accurate_than": [
64+ "bfs-n8"
65+ ],
66+ "hardware": "Intel Xeon workstation, single process, numpy/scipy (LAPACK eigh)",
67+ "notes": "lowest nonrigid mode frequencies (Hz): 36.32, 53.006, 67.184, 94.368, 94.368, 167.77"
68+ },
69+ {
70+ "solution_id": "bfs-n16",
71+ "output": "kirchhoff-bfs/bfs-n16.wav",
72+ "method": "Kirchhoff-Love plate, Bogner-Fox-Schmit C1 finite elements, closed-form modal time response, exact retarded Rayleigh integral",
73+ "parameters": {
74+ "elements_per_side": 16,
75+ "element_size_m": 0.01905,
76+ "n_dof": 1156,
77+ "n_modes": 1156,
78+ "rigid_modes": 3,
79+ "mode_cutoff_hz": 200000.0
80+ },
81+ "timing_seconds": 17.946,
82+ "more_accurate_than": [
83+ "bfs-n12"
84+ ],
85+ "hardware": "Intel Xeon workstation, single process, numpy/scipy (LAPACK eigh)",
86+ "notes": "lowest nonrigid mode frequencies (Hz): 36.32, 53.005, 67.183, 94.366, 94.366, 167.755"
87+ },
88+ {
89+ "solution_id": "bfs-n24",
90+ "output": "kirchhoff-bfs/bfs-n24.wav",
91+ "method": "Kirchhoff-Love plate, Bogner-Fox-Schmit C1 finite elements, closed-form modal time response, exact retarded Rayleigh integral",
92+ "parameters": {
93+ "elements_per_side": 24,
94+ "element_size_m": 0.0127,
95+ "n_dof": 2500,
96+ "n_modes": 2496,
97+ "rigid_modes": 3,
98+ "mode_cutoff_hz": 200000.0
99+ },
100+ "timing_seconds": 63.718,
101+ "more_accurate_than": [
102+ "bfs-n16"
103+ ],
104+ "hardware": "Intel Xeon workstation, single process, numpy/scipy (LAPACK eigh)",
105+ "notes": "lowest nonrigid mode frequencies (Hz): 36.32, 53.005, 67.183, 94.365, 94.365, 167.749"
106+ },
107+ {
108+ "solution_id": "bfs-n32",
109+ "output": "kirchhoff-bfs/bfs-n32.wav",
110+ "method": "Kirchhoff-Love plate, Bogner-Fox-Schmit C1 finite elements, closed-form modal time response, exact retarded Rayleigh integral",
111+ "parameters": {
112+ "elements_per_side": 32,
113+ "element_size_m": 0.009525,
114+ "n_dof": 4356,
115+ "n_modes": 4078,
116+ "rigid_modes": 3,
117+ "mode_cutoff_hz": 200000.0
118+ },
119+ "timing_seconds": 154.684,
120+ "more_accurate_than": [
121+ "bfs-n24"
122+ ],
123+ "hardware": "Intel Xeon workstation, single process, numpy/scipy (LAPACK eigh)",
124+ "notes": "lowest nonrigid mode frequencies (Hz): 36.32, 53.005, 67.183, 94.365, 94.365, 167.748"
125+ },
126+ {
127+ "solution_id": "bfs-n48",
128+ "output": "kirchhoff-bfs/bfs-n48.wav",
129+ "method": "Kirchhoff-Love plate, Bogner-Fox-Schmit C1 finite elements, closed-form modal time response, exact retarded Rayleigh integral",
130+ "parameters": {
131+ "elements_per_side": 48,
132+ "element_size_m": 0.00635,
133+ "n_dof": 9604,
134+ "n_modes": 5041,
135+ "rigid_modes": 3,
136+ "mode_cutoff_hz": 200000.0
137+ },
138+ "timing_seconds": 491.104,
139+ "fields": {
140+ "path": "kirchhoff-bfs/bfs-n48.h5",
141+ "format": "hdf5"
142+ },
143+ "more_accurate_than": [
144+ "bfs-n32"
145+ ],
146+ "hardware": "Intel Xeon workstation, single process, numpy/scipy (LAPACK eigh)",
147+ "notes": "lowest nonrigid mode frequencies (Hz): 36.32, 53.005, 67.183, 94.365, 94.365, 167.748"
148+ }
149+ ]
150+}
solutions/plate-strike-v1/kirchhoff-bfs/bfs-n12.wavadded+0−0View file
Binary file not shown.
solutions/plate-strike-v1/kirchhoff-bfs/bfs-n16.wavadded+0−0View file
Binary file not shown.
solutions/plate-strike-v1/kirchhoff-bfs/bfs-n24.wavadded+0−0View file
Binary file not shown.
solutions/plate-strike-v1/kirchhoff-bfs/bfs-n32.wavadded+0−0View file
Binary file not shown.
solutions/plate-strike-v1/kirchhoff-bfs/bfs-n4.wavadded+0−0View file
Binary file not shown.
solutions/plate-strike-v1/kirchhoff-bfs/bfs-n48.h5added+3−0View file
@@ -0,0 +1,3 @@
1+version https://git-lfs.github.com/spec/v1
2+oid sha256:9b648f1205f8c27245dc9afd5c0b1caaf890ac07fd86c6dac8ec25347cb1746e
3+size 38618270
solutions/plate-strike-v1/kirchhoff-bfs/bfs-n48.wavadded+0−0View file
Binary file not shown.
solutions/plate-strike-v1/kirchhoff-bfs/bfs-n8.wavadded+0−0View file
Binary file not shown.
solvers/kirchhoff-bfs/kl_bfs_plate.pyadded+373−0View file
@@ -0,0 +1,373 @@
1+"""Kirchhoff-Love thin-plate solution of plate-strike-v1 with Bogner-Fox-Schmit finite elements.
2+
3+Space: a uniform n x n mesh of C1-conforming bicubic Hermite rectangles (Bogner-Fox-Schmit),
4+four degrees of freedom per node (w, w_x, w_y, w_xy). Free edges are natural boundary
5+conditions of the Kirchhoff energy, so nothing is imposed. Element matrices are integrated
6+exactly with 4 x 4 Gauss-Legendre points. The generalized eigenproblem gives mass-normalized
7+modes, three of them rigid-body modes at zero frequency.
8+
9+Time: each modal equation
10+
11+ q'' + 2 alpha_n q' + omega_n^2 q = f_n psi(t / T_s), alpha_n = sigma0 + eta omega_n^2 / 2,
12+
13+is solved in closed form: a particular solution for the Hann pulse (a constant plus a cosine)
14+during the strike, and free decay afterwards.
15+
16+Radiation: the Rayleigh integral is evaluated with 3 x 3 Gauss-Legendre points per element.
17+The retardation is exact: once the strike has ended at every retarded time the pressure is a sum
18+of complex exponentials whose amplitudes absorb the delays; the few earlier samples are evaluated
19+directly at the retarded times.
20+
21+The transverse displacement w is positive in +z, so the strike enters as -p and w is reported
22+directly as u_z. Modes above 200 kHz are discarded; the pulse excites them at a negligible level
23+and their extreme damping rates would otherwise overflow the delay factors.
24+
25+Usage: python kl_bfs_plate.py <challenge.json> <out-dir> [--fields n,...] [n ...]
26+"""
27+from __future__ import annotations
28+
29+import json
30+import sys
31+import time
32+from pathlib import Path
33+
34+import numpy as np
35+from numpy.polynomial import legendre as leg
36+from scipy.io import wavfile
37+from scipy.linalg import eigh
38+
39+MODE_CUTOFF_HZ = 200_000.0
40+RIGID_OMEGA = 1.0 # rad/s; eigenvalues below this are the three rigid-body modes
41+
42+
43+def load_challenge(path: Path) -> dict:
44+ with open(path) as f:
45+ return json.load(f)
46+
47+
48+# ---------------------------------------------------------------- Hermite cubics on [0, 1]
49+def hermite(s: np.ndarray):
50+ """Values and first two derivatives (w.r.t. s) of the four cubic Hermite functions. Shapes (len(s), 4)."""
51+ s = np.asarray(s, dtype=float)
52+ H = np.stack([1 - 3 * s**2 + 2 * s**3, s - 2 * s**2 + s**3, 3 * s**2 - 2 * s**3, -(s**2) + s**3], axis=-1)
53+ dH = np.stack([-6 * s + 6 * s**2, 1 - 4 * s + 3 * s**2, 6 * s - 6 * s**2, -2 * s + 3 * s**2], axis=-1)
54+ d2H = np.stack([-6 + 12 * s, -4 + 6 * s, 6 - 12 * s, -2 + 6 * s], axis=-1)
55+ return H, dH, d2H
56+
57+
58+def shape_functions(xi: np.ndarray, eta: np.ndarray, a: float, b: float):
59+ """BFS shape functions and derivatives at local coordinates (xi, eta) in [0,1]^2 for an a x b element.
60+
61+ Returns N, N_xx, N_yy, N_xy, each of shape (npts, 16). Local dof order: for corner c in
62+ [(0,0), (1,0), (0,1), (1,1)]: w, w_x, w_y, w_xy.
63+ """
64+ Hx, dHx, d2Hx = hermite(xi)
65+ Hy, dHy, d2Hy = hermite(eta)
66+ # per-corner 1D function indices: corner at xi=0 uses (H0, H1), at xi=1 uses (H2, H3)
67+ N = np.zeros((len(xi), 16))
68+ Nxx = np.zeros_like(N)
69+ Nyy = np.zeros_like(N)
70+ Nxy = np.zeros_like(N)
71+ corners = [(0, 0), (1, 0), (0, 1), (1, 1)]
72+ for c, (cx, cy) in enumerate(corners):
73+ ix0, ix1 = (0, 1) if cx == 0 else (2, 3)
74+ iy0, iy1 = (0, 1) if cy == 0 else (2, 3)
75+ # (1D function in x, scale for derivative dof) pairs for dofs w, w_x, w_y, w_xy
76+ combos = [(ix0, iy0, 1.0), (ix1, iy0, a), (ix0, iy1, b), (ix1, iy1, a * b)]
77+ for d, (fx, fy, sc) in enumerate(combos):
78+ k = 4 * c + d
79+ N[:, k] = sc * Hx[:, fx] * Hy[:, fy]
80+ Nxx[:, k] = sc * d2Hx[:, fx] * Hy[:, fy] / a**2
81+ Nyy[:, k] = sc * Hx[:, fx] * d2Hy[:, fy] / b**2
82+ Nxy[:, k] = sc * dHx[:, fx] * dHy[:, fy] / (a * b)
83+ return N, Nxx, Nyy, Nxy
84+
85+
86+class BFSPlate:
87+ def __init__(self, ch: dict, n: int):
88+ p = ch["parameters"]
89+ self.p = p
90+ self.n = n
91+ Lx, Ly, h = p["Lx"], p["Ly"], p["h"]
92+ D = p["E"] * h**3 / (12.0 * (1.0 - p["nu"] ** 2))
93+ m = p["rho"] * h
94+ nu = p["nu"]
95+ self.D, self.m = D, m
96+ a, b = Lx / n, Ly / n
97+ self.a, self.b = a, b
98+
99+ # element matrices (uniform mesh: one of each)
100+ g, wg = leg.leggauss(4)
101+ s = (g + 1) / 2
102+ ws = wg / 2
103+ XI, ETA = np.meshgrid(s, s, indexing="ij")
104+ WQ = (ws[:, None] * ws[None, :]).ravel() * a * b
105+ N, Nxx, Nyy, Nxy = shape_functions(XI.ravel(), ETA.ravel(), a, b)
106+ C = np.array([[1.0, nu, 0.0], [nu, 1.0, 0.0], [0.0, 0.0, 2.0 * (1.0 - nu)]])
107+ Ke = np.zeros((16, 16))
108+ Me = np.zeros((16, 16))
109+ for q in range(len(WQ)):
110+ B2 = np.stack([Nxx[q], Nyy[q], Nxy[q]]) # (3, 16)
111+ Ke += WQ[q] * D * (B2.T @ C @ B2)
112+ Me += WQ[q] * m * np.outer(N[q], N[q])
113+
114+ # assembly into dense global matrices
115+ nn = n + 1
116+ ndof = 4 * nn * nn
117+ self.ndof = ndof
118+
119+ def node(i, j):
120+ return i * nn + j
121+
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
130+ self.node = node
131+
132+ evals, evecs = eigh(K, M) # mass-normalized eigenvectors
133+ omega2 = np.maximum(evals, 0.0)
134+ omega = np.sqrt(omega2)
135+ omega[omega < RIGID_OMEGA] = 0.0
136+ keep = omega <= 2 * np.pi * MODE_CUTOFF_HZ
137+ self.omega = omega[keep]
138+ self.V = evecs[:, keep]
139+ self.n_modes = int(keep.sum())
140+ self.alpha = p["sigma0"] + 0.5 * p["eta"] * self.omega**2
141+ self.rigid = self.omega == 0.0
142+ self.n_rigid = int(self.rigid.sum())
143+
144+ # modal forces: f_n = -int p_hat Phi_n dA over the strike disk
145+ R, xs, ys = p["strike_radius"], p["strike_x"], p["strike_y"]
146+ c_phi = 0.5 - 2.0 / np.pi**2
147+ gr, gw = leg.leggauss(24)
148+ r = R * (gr + 1.0) / 2.0
149+ wr = gw * R / 2.0
150+ nth = 64
151+ th = 2 * np.pi * np.arange(nth) / nth
152+ rr, tt = np.meshgrid(r, th, indexing="ij")
153+ ww = (wr[:, None] * rr * (2 * np.pi / nth)).ravel()
154+ px = (xs + rr * np.cos(tt)).ravel()
155+ py = (ys + rr * np.sin(tt)).ravel()
156+ phat = p["strike_peak_force"] / (c_phi * np.pi * R**2) * 0.5 * (1.0 + np.cos(np.pi * rr.ravel() / R))
157+ self.f = -(ww * phat) @ self.modes_at(px, py)
158+
159+ # closed-form response constants
160+ T = p["strike_duration"]
161+ Om = 2 * np.pi / T
162+ self.T, self.Om = T, Om
163+ al, w, f = self.alpha, self.omega, self.f
164+ with np.errstate(divide="ignore", invalid="ignore"):
165+ disc = np.sqrt(al**2 - w**2 + 0j)
166+ self.s1 = -al + disc
167+ self.s2 = -al - disc
168+ H = 1.0 / (w**2 - Om**2 + 2j * al * Om)
169+ self.qconst = np.where(self.rigid, 0.0, f / (2.0 * w**2))
170+ self.Hc = np.where(self.rigid, 0.0, (f / 2.0) * H)
171+ rhs1 = -(self.qconst - self.Hc.real)
172+ rhs2 = (1j * Om * self.Hc).real
173+ self.B = (rhs2 - self.s1 * rhs1) / (self.s2 - self.s1)
174+ self.A = rhs1 - self.B
175+ self.Gc = (f / 2.0) / (2.0 * al + 1j * Om)
176+ self.vconst = f / (4.0 * al)
177+ self.Ar = -(self.vconst - self.Gc.real)
178+ qT, vT, _ = self._pulse_phase(np.array([T]))
179+ self.qT, self.vT = qT[:, 0], vT[:, 0]
180+ with np.errstate(divide="ignore", invalid="ignore"):
181+ self.B2 = (self.vT - self.s1 * self.qT) / (self.s2 - self.s1)
182+ self.A2 = self.qT - self.B2
183+
184+ # ------------------------------------------------------------ evaluation at points
185+ def basis_at(self, x: np.ndarray, y: np.ndarray) -> np.ndarray:
186+ """Global shape-function values at points, shape (npts, ndof), dense."""
187+ n, a, b = self.n, self.a, self.b
188+ ex = np.clip(np.floor(x / a).astype(int), 0, n - 1)
189+ ey = np.clip(np.floor(y / b).astype(int), 0, n - 1)
190+ xi = x / a - ex
191+ eta = y / b - ey
192+ N, _, _, _ = shape_functions(xi, eta, a, b)
193+ out = np.zeros((len(x), self.ndof))
194+ rows = np.arange(len(x))
195+ for c, (cx, cy) in enumerate([(0, 0), (1, 0), (0, 1), (1, 1)]):
196+ nd = self.node(ex + cx, ey + cy)
197+ for d in range(4):
198+ out[rows, 4 * nd + d] += N[:, 4 * c + d]
199+ return out
200+
201+ def modes_at(self, x: np.ndarray, y: np.ndarray) -> np.ndarray:
202+ return self.basis_at(x, y) @ self.V
203+
204+ # ------------------------------------------------------------ modal time response
205+ def _pulse_phase(self, t: np.ndarray):
206+ t = t[None, :]
207+ Om = self.Om
208+ e = np.exp(1j * Om * t)
209+ s1, s2, A, B = self.s1[:, None], self.s2[:, None], self.A[:, None], self.B[:, None]
210+ E1, E2 = np.exp(s1 * t), np.exp(s2 * t)
211+ Hc = self.Hc[:, None]
212+ q = (A * E1 + B * E2).real + self.qconst[:, None] - (Hc * e).real
213+ v = (A * s1 * E1 + B * s2 * E2).real - (1j * Om * Hc * e).real
214+ acc = (A * s1**2 * E1 + B * s2**2 * E2).real + (Om**2 * Hc * e).real
215+ rg = self.rigid
216+ if rg.any():
217+ al = self.alpha[rg][:, None]
218+ Ar, Gc, vc = self.Ar[rg][:, None], self.Gc[rg][:, None], self.vconst[rg][:, None]
219+ Ed = np.exp(-2.0 * al * t)
220+ q[rg] = Ar * (1.0 - Ed) / (2.0 * al) + vc * t - ((e - 1.0) * Gc / (1j * Om)).real
221+ v[rg] = Ar * Ed + vc - (Gc * e).real
222+ acc[rg] = -2.0 * al * Ar * Ed - (1j * Om * Gc * e).real
223+ return q, v, acc
224+
225+ def _free_phase(self, tau: np.ndarray):
226+ tau = tau[None, :]
227+ s1, s2, A, B = self.s1[:, None], self.s2[:, None], self.A2[:, None], self.B2[:, None]
228+ E1, E2 = np.exp(s1 * tau), np.exp(s2 * tau)
229+ q = (A * E1 + B * E2).real
230+ v = (A * s1 * E1 + B * s2 * E2).real
231+ acc = (A * s1**2 * E1 + B * s2**2 * E2).real
232+ rg = self.rigid
233+ if rg.any():
234+ al = self.alpha[rg][:, None]
235+ qT, vT = self.qT[rg][:, None], self.vT[rg][:, None]
236+ Ed = np.exp(-2.0 * al * tau)
237+ q[rg] = qT + vT * (1.0 - Ed) / (2.0 * al)
238+ v[rg] = vT * Ed
239+ acc[rg] = -2.0 * al * vT * Ed
240+ return q, v, acc
241+
242+ def response(self, t: np.ndarray):
243+ q = np.zeros((self.n_modes, len(t)))
244+ v = np.zeros_like(q)
245+ acc = np.zeros_like(q)
246+ m1 = (t >= 0) & (t <= self.T)
247+ m2 = t > self.T
248+ if m1.any():
249+ q[:, m1], v[:, m1], acc[:, m1] = self._pulse_phase(t[m1])
250+ if m2.any():
251+ q[:, m2], v[:, m2], acc[:, m2] = self._free_phase(t[m2] - self.T)
252+ return q, v, acc
253+
254+ # ------------------------------------------------------------ Rayleigh integral
255+ def pressure(self, fs: int, nsamp: int) -> np.ndarray:
256+ p = self.p
257+ n, a, b = self.n, self.a, self.b
258+ g, wg = leg.leggauss(3)
259+ s = (g + 1) / 2
260+ ws = wg / 2
261+ ex, ey = np.meshgrid(np.arange(n), np.arange(n), indexing="ij")
262+ XI, ETA = np.meshgrid(s, s, indexing="ij")
263+ xq = (ex.ravel()[:, None] + XI.ravel()[None, :]).ravel() * a
264+ yq = (ey.ravel()[:, None] + ETA.ravel()[None, :]).ravel() * b
265+ WW = np.tile((ws[:, None] * ws[None, :]).ravel() * a * b, n * n)
266+ r = np.sqrt((p["listener_x"] - xq) ** 2 + (p["listener_y"] - yq) ** 2 + p["listener_height"] ** 2)
267+ tau = r / p["c_air"]
268+ cq = p["rho_air"] / (2 * np.pi) * WW / r
269+ Phi = self.modes_at(xq, yq) # (nq, n_modes)
270+ t = np.arange(nsamp) / fs
271+ out = np.zeros(nsamp)
272+ t_split = self.T + tau.max()
273+ early = (t < t_split) & (t >= tau.min())
274+ late = t >= t_split
275+ ne = int(early.sum())
276+ if ne:
277+ chunk = max(1, 4_000_000 // max(1, self.n_modes * ne))
278+ for i0 in range(0, len(xq), chunk):
279+ sl = slice(i0, i0 + chunk)
280+ tt = (t[early][None, :] - tau[sl][:, None]).ravel()
281+ _, _, acc = self.response(tt)
282+ acc = acc.reshape(self.n_modes, -1, ne)
283+ out[early] += np.einsum("q,qn,nqk->k", cq[sl], Phi[sl], acc, optimize=True)
284+ if late.any():
285+ tl = t[late] - self.T
286+ rg = self.rigid
287+ G1 = (cq[:, None] * Phi * np.exp(-np.outer(tau, self.s1))).sum(axis=0)
288+ G2 = (cq[:, None] * Phi * np.exp(-np.outer(tau, self.s2))).sum(axis=0)
289+ c1 = self.A2 * self.s1**2 * G1
290+ c2 = self.B2 * self.s2**2 * G2
291+ s1, s2 = self.s1.copy(), self.s2.copy()
292+ if rg.any():
293+ al = self.alpha[rg]
294+ Gr = (cq[:, None] * Phi[:, rg] * np.exp(np.outer(tau, 2.0 * al))).sum(axis=0)
295+ c1[rg] = -2.0 * al * self.vT[rg] * Gr
296+ c2[rg] = 0.0
297+ s1[rg] = -2.0 * al
298+ s2[rg] = 0.0
299+ coef = np.concatenate([c1, c2])
300+ rates = np.concatenate([s1, s2])
301+ res = np.zeros(len(tl))
302+ step = max(1, 20_000_000 // len(coef))
303+ for k0 in range(0, len(tl), step):
304+ tk = tl[k0:k0 + step]
305+ res[k0:k0 + step] = (np.exp(np.outer(tk, rates)) @ coef).real
306+ out[late] = res
307+ return out
308+
309+ def field(self, x: np.ndarray, y: np.ndarray, t: np.ndarray) -> np.ndarray:
310+ XX, YY = np.meshgrid(x, y, indexing="xy") # (ny, nx)
311+ Phi = self.modes_at(XX.ravel(), YY.ravel())
312+ q, _, _ = self.response(t)
313+ return (Phi @ q).T.reshape(len(t), len(y), len(x))
314+
315+
316+def main():
317+ ch = load_challenge(Path(sys.argv[1]))
318+ out_dir = Path(sys.argv[2])
319+ args = sys.argv[3:]
320+ fields_for = set()
321+ if "--fields" in args:
322+ i = args.index("--fields")
323+ fields_for = {int(v) for v in args[i + 1].split(",")}
324+ args = args[:i] + args[i + 2:]
325+ meshes = [int(a) for a in args] or [4, 6, 8, 12, 16, 24, 32]
326+ out = ch["output"]
327+ fs, nsamp = out["sample_rate_hz"], out["num_samples"]
328+ meta = []
329+ for n in meshes:
330+ t0 = time.perf_counter()
331+ plate = BFSPlate(ch, n)
332+ t1 = time.perf_counter()
333+ s = plate.pressure(fs, nsamp)
334+ dt = time.perf_counter() - t0
335+ sid = f"bfs-n{n}"
336+ out_dir.mkdir(parents=True, exist_ok=True)
337+ wavfile.write(out_dir / f"{sid}.wav", fs, (s / out["signal_scale"]).astype(np.float32))
338+ low = plate.omega[plate.omega > 0][:6] / (2 * np.pi)
339+ entry = {
340+ "solution_id": sid,
341+ "output": f"{out_dir.name}/{sid}.wav",
342+ "method": "Kirchhoff-Love plate, Bogner-Fox-Schmit C1 finite elements, closed-form modal time response, exact retarded Rayleigh integral",
343+ "parameters": {"elements_per_side": n, "element_size_m": round(plate.a, 6), "n_dof": plate.ndof, "n_modes": plate.n_modes, "rigid_modes": plate.n_rigid, "mode_cutoff_hz": MODE_CUTOFF_HZ},
344+ "timing_seconds": round(dt, 3),
345+ "peak_si": float(np.max(np.abs(s))),
346+ "lowest_frequencies_hz": [round(float(v), 3) for v in low],
347+ }
348+ if n in fields_for:
349+ import h5py
350+
351+ fg = ch["fields"]
352+ nt, ny, nx = fg["shape"]
353+ x = np.linspace(0, ch["parameters"]["Lx"], nx)
354+ y = np.linspace(0, ch["parameters"]["Ly"], ny)
355+ tt = np.arange(nt) / fg["frame_rate_hz"]
356+ u = plate.field(x, y, tt)
357+ with h5py.File(out_dir / f"{sid}.h5", "w") as f:
358+ f.attrs["challenge_id"] = ch["challenge_id"]
359+ f.attrs["solution_id"] = sid
360+ f.attrs["quantity"] = fg["quantity"]
361+ f.attrs["units"] = fg["units"]
362+ f.create_dataset("x", data=x)
363+ f.create_dataset("y", data=y)
364+ f.create_dataset("t", data=tt)
365+ f.create_dataset("u", data=u.astype(np.float32), compression="gzip", compression_opts=4, chunks=(8, ny, nx))
366+ entry["fields"] = {"path": f"{out_dir.name}/{sid}.h5", "format": "hdf5"}
367+ meta.append(entry)
368+ print(f"n={n:3d} dof={plate.ndof:5d} modes={plate.n_modes:5d} rigid={plate.n_rigid} setup={t1 - t0:6.1f}s total={dt:6.1f}s peak={np.max(np.abs(s)):.4e} Pa f_low={np.round(low, 2)}", file=sys.stderr)
369+ print(json.dumps(meta, indent=1))
370+
371+
372+if __name__ == "__main__":
373+ main()
web/build.pyadded+60−0View file
@@ -0,0 +1,60 @@
1+"""Assemble this repository's static site: a manifest of the submissions, the solution files, and the viewer page."""
2+from __future__ import annotations
3+
4+import argparse
5+import json
6+import shutil
7+from pathlib import Path
8+
9+ROOT = Path(__file__).resolve().parents[1]
10+LFS_PREFIX = b"version https://git-lfs.github.com/spec/v1"
11+
12+
13+def is_lfs_pointer(p: Path) -> bool:
14+ try:
15+ with open(p, "rb") as f:
16+ return f.read(len(LFS_PREFIX)) == LFS_PREFIX
17+ except OSError:
18+ return False
19+
20+
21+def main():
22+ ap = argparse.ArgumentParser()
23+ ap.add_argument("--out", type=Path, default=ROOT / "site")
24+ args = ap.parse_args()
25+ out: Path = args.out
26+ if out.exists():
27+ shutil.rmtree(out)
28+ out.mkdir(parents=True)
29+ manifest = {"challenges": {}}
30+ for p in sorted((ROOT / "solutions").glob("*/*.json")):
31+ sub = json.loads(p.read_text())
32+ cid = p.parent.name
33+ entry = manifest["challenges"].setdefault(cid, {"submissions": []})
34+ rel_dir = p.parent.relative_to(ROOT)
35+ for sol in sub["solutions"]:
36+ sol["output"] = f"{rel_dir.as_posix()}/{sol['output']}"
37+ if "fields" in sol and not sol["fields"]["path"].startswith("http"):
38+ local = p.parent / sol["fields"]["path"]
39+ sol["fields"]["path"] = f"{rel_dir.as_posix()}/{sol['fields']['path']}"
40+ sol["fields"]["available"] = local.exists() and not is_lfs_pointer(local)
41+ entry["submissions"].append(sub)
42+ shutil.copytree(ROOT / "solutions", out / "solutions")
43+ shutil.copyfile(ROOT / "web" / "index.html", out / "index.html")
44+ shutil.copyfile(ROOT / "web" / "viewer.js", out / "viewer.js")
45+ # the shared Web Audio player, from the installed challenge-tools package
46+ from challenge_tools.build import player_js
47+
48+ (out / "player.js").write_text(player_js())
49+ # suggested playback gains, when the challenge repository is checked out beside us (as in the workflow)
50+ gains = {}
51+ for cj in sorted((ROOT / "challenge-repo" / "challenges").glob("*/challenge.json")) if (ROOT / "challenge-repo").exists() else []:
52+ ch = json.loads(cj.read_text())
53+ gains[ch["challenge_id"]] = (ch.get("playback") or {}).get("gain_db")
54+ manifest["playback_gain_db"] = gains
55+ (out / "manifest.json").write_text(json.dumps(manifest, indent=1))
56+ print(f"site written to {out}")
57+
58+
59+if __name__ == "__main__":
60+ main()
web/index.htmladded+53−0View file
@@ -0,0 +1,53 @@
1+<!doctype html>
2+<html lang="en">
3+<head>
4+<meta charset="utf-8">
5+<meta name="viewport" content="width=device-width, initial-scale=1">
6+<title>Plate solutions</title>
7+<style>
8+:root { --fg: #1b1b1b; --bg: #fff; --muted: #666; --line: #ddd; --accent: #2a6db0; }
9+@media (prefers-color-scheme: dark) { :root { --fg: #e6e6e6; --bg: #151515; --muted: #999; --line: #333; --accent: #7ab3e8; } }
10+body { font-family: system-ui, sans-serif; color: var(--fg); background: var(--bg); margin: 0; padding: 16px; line-height: 1.45; }
11+main { max-width: 1100px; margin: 0 auto; }
12+a { color: var(--accent); }
13+table { border-collapse: collapse; width: 100%; font-size: 0.92em; }
14+th, td { border-bottom: 1px solid var(--line); padding: 6px 8px; text-align: left; vertical-align: top; }
15+td.num { text-align: right; font-variant-numeric: tabular-nums; }
16+.muted { color: var(--muted); }
17+.wrap { overflow-x: auto; }
18+#view { width: 100%; height: 480px; border: 1px solid var(--line); border-radius: 4px; touch-action: none; }
19+.controls { display: flex; flex-wrap: wrap; gap: 12px; align-items: center; margin: 8px 0; }
20+.controls label { display: flex; gap: 6px; align-items: center; }
21+input[type=range] { width: 200px; }
22+button { font: inherit; padding: 4px 10px; }
23+</style>
24+<script type="importmap">
25+{ "imports": {
26+ "three": "https://cdn.jsdelivr.net/npm/three@0.160.0/build/three.module.js",
27+ "three/addons/": "https://cdn.jsdelivr.net/npm/three@0.160.0/examples/jsm/"
28+} }
29+</script>
30+</head>
31+<body><main>
32+<h1>Solutions to the plate challenges</h1>
33+<p>This is the site of the solution repository <a href="https://vault1.magland.org/challenges/plate-solutions">challenges/plate-solutions</a>. The comparison of all registered solutions against the reference is on the <a href="https://vault1.magland.org/challenges/plate/site/">challenge repository's results site</a>; this page lists this repository's own submissions and animates the exported displacement fields as a 3D surface.</p>
34+
35+<h2>Displacement animation</h2>
36+<div class="controls">
37+ <label>Solution <select id="fieldSelect"></select></label>
38+ <button id="playBtn">Play</button>
39+ <label>Speed <input id="speed" type="range" min="-3" max="0" step="0.1" value="-2"> <span id="speedLabel"></span></label>
40+ <label>Exaggeration <input id="exag" type="range" min="0" max="3" step="0.1" value="1.5"> <span id="exagLabel"></span></label>
41+ <label>Time <input id="timeSlider" type="range" min="0" max="1" step="1" value="0"> <span id="timeLabel" class="muted"></span></label>
42+</div>
43+<div id="view"></div>
44+<p id="fieldStatus" class="muted">Loading…</p>
45+
46+<h2>Submissions</h2>
47+<div id="playback"></div>
48+<div id="submissions"></div>
49+</main>
50+<script src="player.js"></script>
51+<script type="module" src="viewer.js"></script>
52+</body>
53+</html>
web/viewer.jsadded+216−0View file
@@ -0,0 +1,216 @@
1+// Lists this repository's submissions and animates the exported HDF5 displacement fields as a 3D surface.
2+import * as THREE from "three";
3+import { OrbitControls } from "three/addons/controls/OrbitControls.js";
4+import h5wasm from "https://cdn.jsdelivr.net/npm/h5wasm@0.7/dist/esm/hdf5_hl.js";
5+
6+const manifest = await (await fetch("manifest.json")).json();
7+const esc = (s) => String(s).replace(/[&<>"']/g, (c) => ({ "&": "&", "<": "<", ">": ">", '"': """, "'": "'" }[c]));
8+
9+// ---------------------------------------------------------------- submissions table
10+const subsDiv = document.getElementById("submissions");
11+const fieldSolutions = [];
12+for (const [cid, ch] of Object.entries(manifest.challenges)) {
13+ const h3 = document.createElement("h3");
14+ h3.textContent = cid;
15+ subsDiv.appendChild(h3);
16+ for (const sub of ch.submissions) {
17+ const p = document.createElement("p");
18+ p.innerHTML = `<b>${esc(sub.submission_id)}</b> by ${esc(sub.submitter.name)}${sub.date ? " (" + esc(sub.date) + ")" : ""}. ${esc(sub.description || "")}` +
19+ (sub.code_url ? ` <a href="${esc(sub.code_url)}">code</a>` : "");
20+ subsDiv.appendChild(p);
21+ const rows = sub.solutions.map((s) => {
22+ if (s.fields && s.fields.available !== false) fieldSolutions.push({ label: `${sub.submission_id}/${s.solution_id}`, path: s.fields.path });
23+ const fieldsLink = s.fields ? (s.fields.available === false ? ' <span class="muted">(fields stored in LFS, not published here)</span>' : ` <a href="${esc(s.fields.path)}">fields</a>`) : "";
24+ return `<tr><td>${esc(s.solution_id)}</td><td>${esc(s.method || "")}<br><span class="muted"><code>${esc(JSON.stringify(s.parameters || {}))}</code></span></td>` +
25+ `<td class="num">${s.timing_seconds ?? ""}</td><td class="muted">${esc((s.more_accurate_than || []).join(", "))}</td>` +
26+ `<td><div class="cplayer" data-src="${esc(s.output)}"></div><br><a href="${esc(s.output)}">wav</a>${fieldsLink}</td></tr>`;
27+ }).join("");
28+ const wrap = document.createElement("div");
29+ wrap.className = "wrap";
30+ wrap.innerHTML = `<table><tr><th>Solution</th><th>Method</th><th>Time (s)</th><th>Claimed more accurate than</th><th>Output</th></tr>${rows}</table>`;
31+ subsDiv.appendChild(wrap);
32+ }
33+}
34+
35+const gains = manifest.playback_gain_db || {};
36+const firstGain = Object.values(gains).find((g) => typeof g === "number");
37+ChallengePlayer.controls(document.getElementById("playback"), { gainDb: typeof firstGain === "number" ? firstGain : 20, limiter: true });
38+ChallengePlayer.mount();
39+
40+// ---------------------------------------------------------------- 3D viewer
41+const view = document.getElementById("view");
42+const status = document.getElementById("fieldStatus");
43+const select = document.getElementById("fieldSelect");
44+const playBtn = document.getElementById("playBtn");
45+const speed = document.getElementById("speed");
46+const speedLabel = document.getElementById("speedLabel");
47+const exag = document.getElementById("exag");
48+const exagLabel = document.getElementById("exagLabel");
49+const timeSlider = document.getElementById("timeSlider");
50+const timeLabel = document.getElementById("timeLabel");
51+
52+const renderer = new THREE.WebGLRenderer({ antialias: true });
53+renderer.setPixelRatio(window.devicePixelRatio || 1);
54+view.appendChild(renderer.domElement);
55+const scene = new THREE.Scene();
56+const dark = matchMedia("(prefers-color-scheme: dark)").matches;
57+scene.background = new THREE.Color(dark ? 0x151515 : 0xffffff);
58+const camera = new THREE.PerspectiveCamera(40, 2, 0.01, 100);
59+camera.position.set(0.35, -0.5, 0.35);
60+camera.up.set(0, 0, 1);
61+const controls = new OrbitControls(camera, renderer.domElement);
62+controls.target.set(0, 0, 0);
63+scene.add(new THREE.AmbientLight(0xffffff, 0.6));
64+const dl = new THREE.DirectionalLight(0xffffff, 0.9);
65+dl.position.set(0.5, -0.8, 1.0);
66+scene.add(dl);
67+
68+let field = null; // {x, y, t, u, nx, ny, nt, umax}
69+let mesh = null;
70+let frame = 0;
71+let playing = false;
72+let lastTs = null;
73+let simTime = 0;
74+
75+function resize() {
76+ const w = view.clientWidth, h = view.clientHeight;
77+ renderer.setSize(w, h, false);
78+ camera.aspect = w / h;
79+ camera.updateProjectionMatrix();
80+}
81+window.addEventListener("resize", () => { resize(); render(); });
82+resize();
83+
84+function buildMesh() {
85+ if (mesh) { scene.remove(mesh); mesh.geometry.dispose(); mesh.material.dispose(); }
86+ const { x, y, nx, ny } = field;
87+ const Lx = x[nx - 1], Ly = y[ny - 1];
88+ const geo = new THREE.PlaneGeometry(Lx, Ly, nx - 1, ny - 1);
89+ geo.translate(0, 0, 0);
90+ const colors = new Float32Array(nx * ny * 3);
91+ geo.setAttribute("color", new THREE.BufferAttribute(colors, 3));
92+ const mat = new THREE.MeshStandardMaterial({ vertexColors: true, side: THREE.DoubleSide, metalness: 0.1, roughness: 0.7 });
93+ mesh = new THREE.Mesh(geo, mat);
94+ scene.add(mesh);
95+ // a thin frame at rest height for reference
96+ const edges = new THREE.EdgesGeometry(new THREE.PlaneGeometry(Lx, Ly, 1, 1));
97+ const frameLine = new THREE.LineSegments(edges, new THREE.LineBasicMaterial({ color: dark ? 0x555555 : 0xbbbbbb }));
98+ mesh.add(frameLine);
99+ // strike marker
100+ const marker = new THREE.Mesh(new THREE.SphereGeometry(0.004, 12, 12), new THREE.MeshBasicMaterial({ color: 0xd04040 }));
101+ marker.position.set(0.09144 - Lx / 2, 0.10668 - Ly / 2, 0.004);
102+ mesh.add(marker);
103+}
104+
105+function updateSurface() {
106+ if (!field || !mesh) return;
107+ const { u, nx, ny, umax } = field;
108+ const pos = mesh.geometry.attributes.position;
109+ const col = mesh.geometry.attributes.color;
110+ const scale = umax > 0 ? (0.03 * Math.pow(10, parseFloat(exag.value)) / umax) * 0.1 : 0; // exaggerated height
111+ const base = frame * nx * ny;
112+ // PlaneGeometry vertices run row by row from +y (top) to -y; row r corresponds to y index (ny-1-r)
113+ for (let r = 0; r < ny; r++) {
114+ const j = ny - 1 - r;
115+ for (let i = 0; i < nx; i++) {
116+ const v = u[base + j * nx + i];
117+ const k = r * nx + i;
118+ const z = v * scale;
119+ pos.setZ(k, z);
120+ const s = umax > 0 ? Math.max(-1, Math.min(1, v / umax)) : 0;
121+ // blue (negative) to white to red (positive)
122+ col.setXYZ(k, s > 0 ? 1 : 1 + s * 0.8, 1 - Math.abs(s) * 0.7, s < 0 ? 1 : 1 - s * 0.8);
123+ }
124+ }
125+ pos.needsUpdate = true;
126+ col.needsUpdate = true;
127+ mesh.geometry.computeVertexNormals();
128+ timeSlider.value = String(frame);
129+ timeLabel.textContent = `${(field.t[frame] * 1e3).toFixed(0)} ms`;
130+}
131+
132+function render() { renderer.render(scene, camera); }
133+controls.addEventListener("change", render);
134+
135+async function loadField(path) {
136+ status.textContent = `Loading ${path} (this may take a while, the file is large)…`;
137+ playing = false;
138+ playBtn.textContent = "Play";
139+ try {
140+ const { FS } = await h5wasm.ready;
141+ const resp = await fetch(path);
142+ const buf = new Uint8Array(await resp.arrayBuffer());
143+ if (buf.length < 200 && new TextDecoder().decode(buf).startsWith("version https://git-lfs")) {
144+ status.textContent = `${path} is a git LFS pointer; the object was not published with this site.`;
145+ return;
146+ }
147+ const name = "current.h5";
148+ try { FS.unlink(name); } catch (_) {}
149+ FS.writeFile(name, buf);
150+ const file = new h5wasm.File(name, "r");
151+ const x = Float64Array.from(file.get("x").value);
152+ const y = Float64Array.from(file.get("y").value);
153+ const t = Float64Array.from(file.get("t").value);
154+ const u = Float32Array.from(file.get("u").value);
155+ 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 };
159+ timeSlider.max = String(field.nt - 1);
160+ frame = 0;
161+ simTime = 0;
162+ buildMesh();
163+ updateSurface();
164+ render();
165+ status.textContent = `${path}: ${field.nx} × ${field.ny} points × ${field.nt} frames, max |u| = ${(umax * 1e3).toFixed(3)} mm`;
166+ } catch (e) {
167+ status.textContent = `Could not load ${path}: ${e}`;
168+ }
169+}
170+
171+function tick(ts) {
172+ if (!playing || !field) return;
173+ if (lastTs !== null) {
174+ const factor = Math.pow(10, parseFloat(speed.value));
175+ simTime += ((ts - lastTs) / 1000) * factor;
176+ const tEnd = field.t[field.nt - 1];
177+ if (simTime > tEnd) simTime = 0;
178+ const dt = field.t[1] - field.t[0];
179+ frame = Math.min(field.nt - 1, Math.round(simTime / dt));
180+ updateSurface();
181+ render();
182+ }
183+ lastTs = ts;
184+ requestAnimationFrame(tick);
185+}
186+
187+playBtn.addEventListener("click", () => {
188+ playing = !playing;
189+ playBtn.textContent = playing ? "Pause" : "Play";
190+ lastTs = null;
191+ if (playing) requestAnimationFrame(tick);
192+});
193+timeSlider.addEventListener("input", () => {
194+ if (!field) return;
195+ frame = parseInt(timeSlider.value, 10);
196+ simTime = field.t[frame];
197+ updateSurface();
198+ render();
199+});
200+const updSpeed = () => { speedLabel.textContent = `${Math.pow(10, parseFloat(speed.value)).toPrecision(2)} × real time`; };
201+const updExag = () => { exagLabel.textContent = `×${Math.pow(10, parseFloat(exag.value)).toPrecision(2)}`; updateSurface(); render(); };
202+speed.addEventListener("input", updSpeed);
203+exag.addEventListener("input", updExag);
204+updSpeed();
205+exagLabel.textContent = `×${Math.pow(10, parseFloat(exag.value)).toPrecision(2)}`;
206+select.addEventListener("change", () => loadField(select.value));
207+
208+for (const f of fieldSolutions) {
209+ const o = document.createElement("option");
210+ o.value = f.path;
211+ o.textContent = f.label;
212+ select.appendChild(o);
213+}
214+render();
215+if (fieldSolutions.length) loadField(fieldSolutions[0].path);
216+else status.textContent = "No solution in this repository exports a displacement field that is published with this site.";