cymbal2-solutions: Hermite finite element solver (direct and modal) for the cymbal2 challenges, site builder, workflow
14 changed files+1756−0
.gitattributesadded+1−0View file
@@ -0,0 +1 @@
1+*.h5 filter=lfs diff=lfs merge=lfs -text
.gitignoreadded+5−0View file
@@ -0,0 +1,5 @@
1+__pycache__/
2+*.partial.wav
3+*.progress.json
4+*.meta.json
5+site/
.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/cymbal2 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+34−0View file
@@ -0,0 +1,34 @@
1+# Solutions to the cymbal2 challenges
2+
3+This repository holds solutions to the challenges defined in [challenges/cymbal2](https://vault1.magland.org/challenges/cymbal2). 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/cymbal2/site](https://vault1.magland.org/challenges/cymbal2/site/). This repository's own site, at [challenges/cymbal2-solutions/site](https://vault1.magland.org/challenges/cymbal2-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+## The submission hermite-fem
16+
17+`solvers/hermite-fem/fem_shell.py` solves the shallow-shell model of MODEL.md directly in the displacements, with no modal reduction. The annulus is approximated by the square cells of a uniform Cartesian grid whose centers lie inside it (a staircase rim and hole). The three fields u, v, w are expanded in bicubic Hermite (Bogner-Fox-Schmit) elements, C1 across cells, four unknowns per node and field. The weak form is integrated with 3 x 3 Gauss points per cell exactly as written: Marguerre membrane strains, bending strains, stress resultants with the Kelvin-Voigt rates, and the assembled virtual work as the internal force, which a test confirms to be the exact gradient of the discrete strain energy; free edges are the natural boundary conditions and need no treatment. The mass is lumped by diagonal scaling of the consistent element mass. The equations are stepped with the explicit leapfrog scheme at 85 percent of the stability limit, which the in-plane block sets; the lumped masses of the in-plane derivative unknowns are scaled by 16, which halves the step count (14 substeps per output sample at 8 mm, about 28 at 4 mm) and moves the in-plane resonances by 0.15 percent while leaving the bending modes unchanged. The rigid constraint of MODEL.md is enforced by applying each stick's reaction as the rigid inertial load that cancels its net force and moments, and by projecting the six mass-orthonormal rigid fields out of the state at every chunk of steps. The sticks are point masses stepped with the same scheme through the Hertz law, each started at the current surface height under it with relative speed v_0. The Rayleigh integral is accumulated at every step from the transverse acceleration at the Gauss points, spread onto the output samples at the retarded times with cubic Lagrange interpolation, so no acceleration history is stored. The step runs on a GPU through PyTorch, the constitutive and output blocks fused by torch.compile and the whole step captured as a CUDA graph: about 5100 steps per second at 8 mm (25488 unknowns, a quarter of an hour for the seven seconds).
18+
19+Checks: the linearized assembly reproduces the free annulus's lowest frequencies to within 1 percent at 8 mm (21.7, 37.6, 50.6, 85.5 Hz against 21.9, 38.0, 50.9, 85.9), the lumped mass shifts them by under 0.3 percent relative to the consistent mass, the internal force matches the energy gradient to 1e-9, and the retarded accumulation matches a direct integral to 4e-5. On the crash geometry the two tilt modes, exact zero modes of the continuum, come out at 2.95 Hz at 8 mm because the bicubic space cannot represent the in-plane displacement that accompanies a rigid tilt of the cosine bell; the constraint removes them anyway. A run of the flat challenge with the shell left free confirmed what MODEL.md says about the constraint: the shell tilted by about a radian in half a second, the partials glided down, the treble drained away, and the explicit run diverged at 0.66 s.
20+
21+```
22+bash solvers/hermite-fem/run_chain.sh <cymbal2-repo> # on a CUDA GPU: 8 mm then 4 mm for both challenges, hours
23+python solvers/hermite-fem/check_modes.py <cymbal2-repo>/challenges/cymbal2-flat-v1/challenge.json 0.008 # linear check, CPU
24+```
25+
26+While a run is going, `<tag>.progress.json` and `<tag>.partial.wav` beside the outputs are rewritten every minute with the state of the sticks and the samples computed so far.
27+
28+## Validating locally
29+
30+```
31+pip install git+https://vault1.magland.org/challenges/challenge-tools
32+challenge-tools validate --challenges <path-to-cymbal2-repo>/challenges --solutions solutions
33+python web/build.py --out site
34+```
solvers/hermite-fem/check_gradient.pyadded+49−0View file
@@ -0,0 +1,49 @@
1+"""Check that the internal force of fem_shell is the gradient of the shell's strain energy (numpy re-implementation of
2+the same Gauss-point evaluation), on the crash geometry at a random state of millimeter amplitude."""
3+import json, sys
4+import numpy as np
5+from fem_shell import ShellFEM
6+ch = json.load(open(sys.argv[1])); a = float(sys.argv[2]) if len(sys.argv) > 2 else 0.008
7+fem = ShellFEM(ch, a, log=lambda s: None)
8+nu = fem.p["nu"]; C, G = fem.n_cells, fem.G
9+B = fem.B # (81, 48)
10+def strains(q):
11+ D = (q[fem.cell_dofs] @ B.T).reshape(C, G, 9)
12+ wx, wy = D[..., 0], D[..., 1]
13+ exx = D[..., 5] + fem.Zx * wx + 0.5 * wx**2
14+ eyy = D[..., 8] + fem.Zy * wy + 0.5 * wy**2
15+ exy = 0.5 * (D[..., 6] + D[..., 7] + fem.Zx * wy + fem.Zy * wx + wx * wy)
16+ return D, exx, eyy, exy
17+def energy(q):
18+ D, exx, eyy, exy = strains(q)
19+ Um = 0.5 * fem.Cm * (exx**2 + eyy**2 + 2 * nu * exx * eyy + 2 * (1 - nu) * exy**2)
20+ Ub = 0.5 * fem.Db * (D[..., 2]**2 + D[..., 3]**2 + 2 * nu * D[..., 2] * D[..., 3] + 2 * (1 - nu) * D[..., 4]**2)
21+ return float(np.sum((Um + Ub) * fem.wq[None, :]))
22+def force(q):
23+ D, exx, eyy, exy = strains(q)
24+ wx, wy = D[..., 0], D[..., 1]; ax, ay = fem.Zx + wx, fem.Zy + wy
25+ Nxx = fem.Cm * (exx + nu * eyy); Nyy = fem.Cm * (eyy + nu * exx); Nxy = fem.Cm * (1 - nu) * exy
26+ Mxx = fem.Db * (D[..., 2] + nu * D[..., 3]); Myy = fem.Db * (D[..., 3] + nu * D[..., 2]); Mxy = fem.Db * (1 - nu) * D[..., 4]
27+ S = np.stack([Nxx * ax + Nxy * ay, Nyy * ay + Nxy * ax, Mxx, Myy, 2 * Mxy, Nxx, Nxy, Nxy, Nyy], axis=-1) * fem.wq[None, :, None]
28+ fe = S.reshape(C, G * 9) @ B
29+ f = np.zeros(fem.n_dof); np.add.at(f, fem.cell_dofs.ravel(), fe.ravel()); return f
30+rng = np.random.default_rng(1)
31+q = np.zeros(fem.n_dof)
32+# smooth random state: w of a few mm, in-plane of 0.1 mm, derivative dofs scaled accordingly
33+xy = fem.node_xy; R = fem.p["disk_radius"]
34+for blk, amp in [(0, 3e-3), (1, 1e-4), (2, 1e-4)]:
35+ for k in range(6):
36+ kx, ky, ph = rng.uniform(-40, 40, 2).tolist() + [rng.uniform(0, 2 * np.pi)]
37+ val = amp * np.cos(kx * xy[:, 0] + ky * xy[:, 1] + ph)
38+ base = blk * 4 * fem.n_nodes + 4 * np.arange(fem.n_nodes)
39+ q[base] += val; q[base + 1] += -amp * kx * np.sin(kx * xy[:, 0] + ky * xy[:, 1] + ph)
40+ q[base + 2] += -amp * ky * np.sin(kx * xy[:, 0] + ky * xy[:, 1] + ph); q[base + 3] += -amp * kx * ky * np.cos(kx * xy[:, 0] + ky * xy[:, 1] + ph)
41+f = force(q)
42+print(f"energy {energy(q):.6e} J, |f| {np.linalg.norm(f):.4e}")
43+for trial in range(4):
44+ d = rng.standard_normal(fem.n_dof) * np.abs(q).mean() * 1e-3
45+ d[4 * fem.n_nodes:] *= 0.03
46+ for eps in (1e-2, 1e-3, 1e-4):
47+ num = (energy(q + eps * d) - energy(q - eps * d)) / (2 * eps)
48+ ana = f @ d
49+ print(f" dir {trial} eps {eps:g}: dE/ds numerical {num:.10e} analytic {ana:.10e} rel diff {abs(num - ana) / abs(ana):.2e}")
solvers/hermite-fem/check_modes.pyadded+82−0View file
@@ -0,0 +1,82 @@
1+"""Linear check of the assembly: the lowest natural frequencies of the linearized shell with the lumped and the
2+consistent mass, from the same element operators the time stepping uses. For the flat annulus the free disk's
3+values are 21.9 (2,0), 38.0 (0,0), 50.9 (3,0), 85.9 (1,0), 89.6 (4,0) Hz (Leissa), the 6 mm hole changing them
4+by well under a percent. Usage: python check_modes.py <challenge.json> [spacing_m] [n_modes]"""
5+import json
6+import sys
7+
8+import numpy as np
9+import scipy.sparse as sp
10+import scipy.sparse.linalg as spla
11+
12+from fem_shell import ShellFEM, element_basis
13+
14+ch = json.load(open(sys.argv[1]))
15+a = float(sys.argv[2]) if len(sys.argv) > 2 else 0.008
16+k = int(sys.argv[3]) if len(sys.argv) > 3 else 14
17+fem = ShellFEM(ch, a, log=lambda s: print(s))
18+p = fem.p
19+C, G, nu = fem.n_cells, fem.G, p["nu"]
20+B = fem.B.reshape(G, 9, 48)
21+# L: 9 derivative quantities -> [exx, eyy, gxy, kxx, kyy, 2kxy] at the reference state
22+L = np.zeros((C, G, 6, 9))
23+L[:, :, 0, 0] = fem.Zx; L[:, :, 0, 5] = 1.0
24+L[:, :, 1, 1] = fem.Zy; L[:, :, 1, 8] = 1.0
25+L[:, :, 2, 0] = fem.Zy; L[:, :, 2, 1] = fem.Zx; L[:, :, 2, 6] = 1.0; L[:, :, 2, 7] = 1.0
26+L[:, :, 3, 2] = 1.0; L[:, :, 4, 3] = 1.0; L[:, :, 5, 4] = 2.0
27+Q = np.array([[1, nu, 0], [nu, 1, 0], [0, 0, (1 - nu) / 2]])
28+Dm = np.zeros((C, G, 6, 6))
29+Dm[:, :, :3, :3] = fem.Cm[:, :, None, None] * Q
30+Dm[:, :, 3:, 3:] = fem.Db[:, :, None, None] * Q
31+LB = np.einsum("cgij,gjk->cgik", L, B) # (C, G, 6, 48)
32+Ke = np.einsum("g,cgik,cgij,cgjl->ckl", fem.wq, LB, Dm, LB, optimize=True)
33+rows = np.repeat(fem.cell_dofs, 48, axis=1).ravel()
34+cols = np.tile(fem.cell_dofs, (1, 48)).ravel()
35+K = sp.coo_matrix((Ke.ravel(), (rows, cols)), shape=(fem.n_dof, fem.n_dof)).tocsr()
36+K = 0.5 * (K + K.T)
37+# consistent mass
38+N = element_basis(fem.sx, fem.sy, a)["N"] # (G, 16)
39+Me16 = np.einsum("cg,gi,gj->cij", p["rho"] * fem.hq * fem.wq[None, :], N, N) # (C, 16, 16)
40+Me = np.zeros((C, 48, 48))
41+for b in range(3):
42+ Me[:, 16 * b:16 * b + 16, 16 * b:16 * b + 16] = Me16
43+Mc = sp.coo_matrix((Me.ravel(), (rows, cols)), shape=(fem.n_dof, fem.n_dof)).tocsr()
44+Ml = sp.diags(fem.Mdiag)
45+print(f"total mass: consistent {Mc.sum():.5f}, lumped {fem.Mdiag[0:4 * fem.n_nodes:4].sum():.5f} (w block)")
46+for name, M in [("lumped", Ml), ("consistent", Mc)]:
47+ vals = spla.eigsh(K, k=k, M=M.tocsc(), sigma=-1e3, which="LM", return_eigenvectors=False)
48+ f = np.sqrt(np.maximum(np.sort(vals), 0)) / (2 * np.pi)
49+ print(f"{name:10s}: " + " ".join(f"{x:7.2f}" for x in f))
50+# highest frequency with the lumped mass (stability)
51+vals = spla.eigsh(K, k=1, M=Ml.tocsc(), which="LM", return_eigenvectors=False)
52+print(f"omega_max (lumped) = {np.sqrt(vals[0]):.4e} rad/s")
53+# which block sets the explicit step: bending-only and in-plane-only top frequencies
54+nw = 4 * fem.n_nodes
55+Kw = K[:nw, :nw]; Kp = K[nw:, nw:]
56+for name, Kb, Mb in [("bending block", Kw, sp.diags(fem.Mdiag[:nw])), ("in-plane block", Kp, sp.diags(fem.Mdiag[nw:]))]:
57+ v = spla.eigsh(Kb, k=1, M=Mb.tocsc(), which="LM", return_eigenvectors=False)
58+ print(f"omega_max {name}: {np.sqrt(v[0]):.4e} rad/s")
59+# and the bending block with the derivative-dof masses scaled up by 4 (rotational mass scaling)
60+Ms = fem.Mdiag[:nw].copy(); Ms[1::4] *= 4; Ms[2::4] *= 4; Ms[3::4] *= 4
61+v = spla.eigsh(Kw, k=1, M=sp.diags(Ms).tocsc(), which="LM", return_eigenvectors=False)
62+print(f"omega_max bending with derivative masses x4: {np.sqrt(v[0]):.4e} rad/s")
63+# in-plane block with its derivative-dof masses scaled, and the effect on the low modes of the whole system
64+for sc in (4, 16, 64):
65+ Ms = fem.Mdiag.copy()
66+ for b in (1, 2):
67+ base = nw * b
68+ Ms[base + 1::4][: fem.n_nodes] *= sc; Ms[base + 2::4][: fem.n_nodes] *= sc; Ms[base + 3::4][: fem.n_nodes] *= sc
69+ Msp = sp.diags(Ms)
70+ v = spla.eigsh(Kp, k=1, M=sp.diags(Ms[nw:]).tocsc(), which="LM", return_eigenvectors=False)
71+ low = spla.eigsh(K, k=14, M=Msp.tocsc(), sigma=-1e3, which="LM", return_eigenvectors=False)
72+ fl = np.sqrt(np.maximum(np.sort(low), 0)) / (2 * np.pi)
73+ print(f"in-plane derivative masses x{sc}: omega_max in-plane {np.sqrt(v[0]):.4e} rad/s; low modes " + " ".join(f"{x:6.2f}" for x in fl[6:]))
74+# lowest in-plane modes with and without the derivative-mass scaling (in-plane block alone; 6 -> 3 rigid modes there)
75+for sc in (1, 16):
76+ Ms = fem.Mdiag[nw:].copy()
77+ for b in (0, 1):
78+ base = nw * b
79+ Ms[base + 1::4][: fem.n_nodes] *= sc; Ms[base + 2::4][: fem.n_nodes] *= sc; Ms[base + 3::4][: fem.n_nodes] *= sc
80+ low = spla.eigsh(Kp, k=9, M=sp.diags(Ms).tocsc(), sigma=-1e3, which="LM", return_eigenvectors=False)
81+ fl = np.sqrt(np.maximum(np.sort(low), 0)) / (2 * np.pi)
82+ print(f"in-plane modes, derivative masses x{sc}: " + " ".join(f"{x:7.1f}" for x in fl[3:]))
solvers/hermite-fem/fem_shell.pyadded+692−0View file
@@ -0,0 +1,692 @@
1+"""Displacement-based finite element solution of the cymbal2 shallow-shell model on a Cartesian grid.
2+
3+The annulus is approximated by the square cells of a uniform grid whose centers lie inside it (a staircase
4+outer rim and a staircase hole). The three displacement fields u, v, w are each expanded in bicubic Hermite
5+(Bogner-Fox-Schmit) elements, C1 across cells, with the four nodal unknowns f, f_x, f_y, f_xy per field. The
6+weak form of MODEL.md is integrated with 3 x 3 Gauss points per cell: the membrane strains of Marguerre theory
7+and the bending strains are evaluated at the points, the stress resultants include the Kelvin-Voigt rates, and
8+the internal force is the assembled virtual work. Nothing is linearized and no modes are computed: the mass is
9+lumped (diagonal scaling of the consistent element mass), the equations are stepped explicitly with the
10+leapfrog scheme at a step below the stability limit of the in-plane waves, and the three sticks are point
11+masses integrated with the same step through the Hertz contact law. The Rayleigh integral is accumulated at
12+every step from the transverse acceleration at the Gauss points with cubic Lagrange interpolation in time at
13+the retarded times, so nothing about the acceleration history is stored. Free edges need no treatment: they
14+are the natural boundary conditions of the weak form.
15+
16+The whole time step runs on the GPU through PyTorch (one dense matmul from element unknowns to Gauss-point
17+derivatives, elementwise constitutive evaluation, one matmul back to element forces, a scatter-add), captured
18+in a CUDA graph. The CPU path is the same code and serves for tests.
19+
20+Usage:
21+ python fem_shell.py <challenge.json> <out-dir> [--spacing 0.008] [--seconds S] [--device cuda] [--fields]
22+ [--report 60] [--linear] [--speed V] [--tag NAME] [--no-graph]
23+Writes <out>/<tag>.wav, <tag>.meta.json, and while running <tag>.progress.json and <tag>.partial.wav (the samples
24+computed so far, the rest zero), rewritten every --report seconds of wall time.
25+"""
26+from __future__ import annotations
27+
28+import json
29+import math
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+
38+# ----------------------------------------------------------------------------------------------------------------
39+# geometry of a challenge
40+
41+class Profile:
42+ """Mid-surface height Z(r), its slope, and the thickness h(r) of a cymbal2 challenge."""
43+
44+ def __init__(self, p: dict):
45+ self.R = p["disk_radius"]
46+ self.rh = p["hole_radius"]
47+ self.H = p.get("bow_height", 0.0)
48+ self.rb = p.get("bell_radius", 0.0)
49+ self.hb = p.get("bell_height", 0.0)
50+ self.h0 = p.get("thickness")
51+ self.hc, self.he = p.get("thickness_center"), p.get("thickness_rim")
52+
53+ def Z(self, r):
54+ z = self.H * (1.0 - r * r / self.R**2)
55+ if self.rb > 0:
56+ z = z + np.where(r < self.rb, 0.5 * self.hb * (1.0 + np.cos(np.pi * np.minimum(r, self.rb) / self.rb)), 0.0)
57+ return z
58+
59+ def dZ(self, r):
60+ d = -2.0 * self.H * r / self.R**2
61+ if self.rb > 0:
62+ d = d + np.where(r < self.rb, -0.5 * self.hb * (np.pi / self.rb) * np.sin(np.pi * np.minimum(r, self.rb) / self.rb), 0.0)
63+ return d
64+
65+ def slope(self, x, y):
66+ r = np.hypot(x, y)
67+ d = self.dZ(r)
68+ with np.errstate(invalid="ignore", divide="ignore"):
69+ zx = np.where(r > 0, d * x / r, 0.0)
70+ zy = np.where(r > 0, d * y / r, 0.0)
71+ return zx, zy
72+
73+ def h(self, r):
74+ if self.h0 is not None:
75+ return np.full_like(np.asarray(r, dtype=float), self.h0)
76+ return np.where(r <= self.rb, self.hc, self.hc + (self.he - self.hc) * (r - self.rb) / (self.R - self.rb))
77+
78+
79+# ----------------------------------------------------------------------------------------------------------------
80+# bicubic Hermite element on a square cell of side a
81+
82+def hermite_1d(s, a):
83+ """Cubic Hermite functions on [0, a] at s = x/a: values and x-derivatives up to second order.
84+ Returns arrays (4, n) for H, H', H'' with the ordering [f(0), f'(0), f(a), f'(a)]."""
85+ s = np.asarray(s, dtype=float)
86+ H = np.stack([1 - 3 * s**2 + 2 * s**3, a * (s - 2 * s**2 + s**3), 3 * s**2 - 2 * s**3, a * (-(s**2) + s**3)])
87+ dH = np.stack([-6 * s + 6 * s**2, a * (1 - 4 * s + 3 * s**2), 6 * s - 6 * s**2, a * (-2 * s + 3 * s**2)]) / a
88+ ddH = np.stack([-6 + 12 * s, a * (-4 + 6 * s), 6 - 12 * s, a * (-2 + 6 * s)]) / a**2
89+ return H, dH, ddH
90+
91+
92+def element_basis(sx, sy, a):
93+ """Shape functions of the 16 element unknowns at points (sx, sy) in [0,1]^2: returns a dict of arrays (n, 16)
94+ for N, Nx, Ny, Nxx, Nyy, Nxy. Local unknown k = 4 * node + d with nodes (0,0), (1,0), (0,1), (1,1) and
95+ d in [f, f_x, f_y, f_xy]."""
96+ Hx, dHx, ddHx = hermite_1d(sx, a)
97+ Hy, dHy, ddHy = hermite_1d(sy, a)
98+ out = {k: np.zeros((len(np.atleast_1d(sx)), 16)) for k in ("N", "Nx", "Ny", "Nxx", "Nyy", "Nxy")}
99+ for node, (i, j) in enumerate([(0, 0), (1, 0), (0, 1), (1, 1)]):
100+ for d, (p, q) in enumerate([(0, 0), (1, 0), (0, 1), (1, 1)]): # (x-derivative flag, y-derivative flag)
101+ ix, iy = 2 * i + p, 2 * j + q
102+ k = 4 * node + d
103+ out["N"][:, k] = Hx[ix] * Hy[iy]
104+ out["Nx"][:, k] = dHx[ix] * Hy[iy]
105+ out["Ny"][:, k] = Hx[ix] * dHy[iy]
106+ out["Nxx"][:, k] = ddHx[ix] * Hy[iy]
107+ out["Nyy"][:, k] = Hx[ix] * ddHy[iy]
108+ out["Nxy"][:, k] = dHx[ix] * dHy[iy]
109+ return out
110+
111+
112+# ----------------------------------------------------------------------------------------------------------------
113+
114+class ShellFEM:
115+ NQ = 3 # Gauss points per direction
116+ QUANT = ["w_x", "w_y", "w_xx", "w_yy", "w_xy", "u_x", "u_y", "v_x", "v_y"]
117+
118+ def __init__(self, ch: dict, spacing: float, log=print, linear: bool = False, speed: float | None = None, c_air: float | None = None, no_drift: bool = False, inplane_scale: float = 1.0):
119+ self.ch = ch
120+ self.no_drift = no_drift
121+ self.inplane_scale = inplane_scale
122+ p = dict(ch["parameters"])
123+ if c_air is not None:
124+ p["c_air"] = c_air
125+ self.p = p
126+ self.a = a = spacing
127+ self.log = log
128+ self.linear = linear
129+ self.v0 = p["stick_speed"] if speed is None else speed
130+ pr = self.profile = Profile(p)
131+ R, rh = pr.R, pr.rh
132+ # ---- cells: centers inside the annulus; nodes: their corners; the origin is a node
133+ n = int(math.ceil(R / a))
134+ ii, jj = np.meshgrid(np.arange(-n, n), np.arange(-n, n), indexing="ij")
135+ cx, cy = (ii.ravel() + 0.5) * a, (jj.ravel() + 0.5) * a
136+ rc = np.hypot(cx, cy)
137+ keep = (rc < R) & (rc > rh)
138+ self.cells = np.stack([ii.ravel()[keep], jj.ravel()[keep]], axis=1) # lower-left node index (i, j)
139+ self.n_cells = len(self.cells)
140+ # node map
141+ node_ids = {}
142+ cell_nodes = np.zeros((self.n_cells, 4), dtype=np.int64)
143+ for c, (i, j) in enumerate(self.cells):
144+ for k, (di, dj) in enumerate([(0, 0), (1, 0), (0, 1), (1, 1)]):
145+ key = (int(i + di), int(j + dj))
146+ if key not in node_ids:
147+ node_ids[key] = len(node_ids)
148+ cell_nodes[c, k] = node_ids[key]
149+ self.n_nodes = len(node_ids)
150+ self.node_xy = np.zeros((self.n_nodes, 2))
151+ for (i, j), k in node_ids.items():
152+ self.node_xy[k] = (i * a, j * a)
153+ self.node_ids = node_ids
154+ self.n_dof = 12 * self.n_nodes # blocks w, u, v of 4 n_nodes each
155+ loc = (4 * cell_nodes[:, :, None] + np.arange(4)[None, None, :]).reshape(self.n_cells, 16)
156+ self.cell_dofs = np.concatenate([loc, loc + 4 * self.n_nodes, loc + 8 * self.n_nodes], axis=1) # (C, 48): w | u | v
157+ self.cell_dofs_w = loc
158+ self.area = self.n_cells * a * a
159+ # ---- quadrature and the derivative operator B: (C-independent) (NQ*NQ*9, 48)
160+ g, wg = leg.leggauss(self.NQ)
161+ s = (g + 1) / 2
162+ ws = wg / 2
163+ SX, SY = np.meshgrid(s, s, indexing="ij")
164+ self.sx, self.sy = SX.ravel(), SY.ravel()
165+ self.wq = (ws[:, None] * ws[None, :]).ravel() * a * a # (G,)
166+ G = len(self.wq)
167+ self.G = G
168+ eb = element_basis(self.sx, self.sy, a)
169+ B = np.zeros((G, 9, 48))
170+ B[:, 0, :16], B[:, 1, :16], B[:, 2, :16], B[:, 3, :16], B[:, 4, :16] = eb["Nx"], eb["Ny"], eb["Nxx"], eb["Nyy"], eb["Nxy"]
171+ B[:, 5, 16:32], B[:, 6, 16:32] = eb["Nx"], eb["Ny"]
172+ B[:, 7, 32:48], B[:, 8, 32:48] = eb["Nx"], eb["Ny"]
173+ self.B = B.reshape(G * 9, 48)
174+ self.Nw = eb["N"] # (G, 16) transverse displacement at the Gauss points
175+ # Gauss point coordinates, slopes, thickness, stiffnesses
176+ gx = self.cells[:, 0][:, None] * a + self.sx[None, :] * a # (C, G)
177+ gy = self.cells[:, 1][:, None] * a + self.sy[None, :] * a
178+ self.gx, self.gy = gx, gy
179+ gr = np.hypot(gx, gy)
180+ self.Zx, self.Zy = pr.slope(gx, gy)
181+ self.hq = pr.h(gr)
182+ E, nu = p["E"], p["nu"]
183+ self.Cm = E * self.hq / (1 - nu**2) # membrane stiffness per Gauss point (C, G)
184+ self.Db = E * self.hq**3 / (12 * (1 - nu**2))
185+ self.mass_total = float(np.sum(p["rho"] * self.hq * self.wq[None, :]))
186+ # ---- lumped mass (HRZ diagonal scaling per element, per field)
187+ Mdiag = np.zeros(self.n_dof)
188+ N2 = eb["N"] ** 2 # (G, 16)
189+ for blk in range(3):
190+ me = (p["rho"] * self.hq * self.wq[None, :]) @ N2 # (C, 16) diagonal of the consistent element mass
191+ mcell = np.sum(p["rho"] * self.hq * self.wq[None, :], axis=1) # (C,)
192+ trans = me[:, 0::4].sum(axis=1) # translational entries
193+ me = me * (mcell / trans)[:, None]
194+ np.add.at(Mdiag, self.cell_dofs[:, 16 * blk : 16 * (blk + 1)], me)
195+ # rotational mass scaling of the in-plane derivative unknowns: their lumped masses are tiny and set the explicit
196+ # time step through the in-plane block; scaling them by inplane_scale (16 by default) roughly halves the step count
197+ # and moves the in-plane resonances by about 0.15 percent, the bending modes not at all.
198+ nw = 4 * self.n_nodes
199+ for b in (1, 2):
200+ for d in (1, 2, 3):
201+ Mdiag[nw * b + d::4][: self.n_nodes] *= inplane_scale
202+ self.Mdiag = Mdiag
203+ # ---- strikes: footprint vectors phi_j = int phi N_j dA (fine quadrature on the cells they touch)
204+ self.strikes = []
205+ Rs = p["strike_radius"]
206+ c_phi = 0.5 - 2.0 / np.pi**2
207+ gf, wf = leg.leggauss(12)
208+ sf = (gf + 1) / 2
209+ wsf = wf / 2
210+ FX, FY = np.meshgrid(sf, sf, indexing="ij")
211+ ebf = element_basis(FX.ravel(), FY.ravel(), a)["N"] # (144, 16)
212+ wff = (wsf[:, None] * wsf[None, :]).ravel() * a * a
213+ for st in ch["strikes"]:
214+ xs, ys = st["x"], st["y"]
215+ phi = np.zeros(self.n_dof)
216+ near = np.flatnonzero(np.hypot(cx[keep] - xs, cy[keep] - ys) < Rs + a)
217+ for c in near:
218+ px = self.cells[c, 0] * a + FX.ravel() * a
219+ py = self.cells[c, 1] * a + FY.ravel() * a
220+ rr = np.hypot(px - xs, py - ys)
221+ bump = np.where(rr <= Rs, 0.5 * (1 + np.cos(np.pi * np.minimum(rr / Rs, 1.0))), 0.0) / (c_phi * np.pi * Rs**2)
222+ np.add.at(phi, self.cell_dofs_w[c], (wff * bump) @ ebf)
223+ self.strikes.append(dict(t=st["t"], phi=phi, integral=float(phi[0 : 4 * self.n_nodes : 4].sum())))
224+ # ---- rigid reaction: for a unit downward force at each strike point, the rigid acceleration field of the shell
225+ # (translation along z plus the two tilts) that balances its net force and moments, as a nodal force vector M g.
226+ # With --no-drift this reaction is applied together with the stick force, so the shell gains no momentum or
227+ # angular momentum and does not drift (the constraint of zero rigid motion; the elastic motion is untouched).
228+ nw = 4 * self.n_nodes
229+ mw = self.Mdiag[0:nw:4] # translational masses of the w nodes
230+ xn, yn = self.node_xy[:, 0], self.node_xy[:, 1]
231+ Mt = mw.sum()
232+ xc, yc = (mw * xn).sum() / Mt, (mw * yn).sum() / Mt
233+ Ixx, Iyy, Ixy = (mw * (yn - yc) ** 2).sum(), (mw * (xn - xc) ** 2).sum(), (mw * (xn - xc) * (yn - yc)).sum()
234+ for st in self.strikes:
235+ xs, ys = [s_ for s_ in ch["strikes"] if abs(s_["t"] - st["t"]) < 1e-12][0]["x"], [s_ for s_ in ch["strikes"] if abs(s_["t"] - st["t"]) < 1e-12][0]["y"]
236+ # rigid field a = c0 + cx (x - xc) + cy (y - yc) with net force F and moments F xs, F ys about the center (per unit F)
237+ Minv = np.linalg.inv(np.array([[Iyy, Ixy], [Ixy, Ixx]]))
238+ cx, cy = Minv @ np.array([xs - xc, ys - yc])
239+ g = np.zeros(self.n_dof)
240+ g[0:nw:4] = 1.0 / Mt + cx * (xn - xc) + cy * (yn - yc)
241+ g[1:nw:4] = cx # f_x of a linear field
242+ g[2:nw:4] = cy
243+ st["reaction"] = self.Mdiag * g # nodal force of the rigid acceleration field for unit stick force (upward)
244+ # the six rigid fields (M-orthonormalized): used with --no-drift to project the rigid components out of the state
245+ rig = np.zeros((6, self.n_dof))
246+ rig[0, 0:nw:4] = 1.0 # w translation
247+ rig[1, 0:nw:4], rig[1, 1:nw:4] = xn - xc, 1.0 # tilt about y (w = x)
248+ rig[2, 0:nw:4], rig[2, 2:nw:4] = yn - yc, 1.0 # tilt about x (w = y)
249+ rig[3, nw:2 * nw:4] = 1.0 # u translation
250+ rig[4, 2 * nw:3 * nw:4] = 1.0 # v translation
251+ rig[5, nw:2 * nw:4], rig[5, nw + 2:2 * nw:4] = -(yn - yc), -1.0 # in-plane rotation: u = -y, v = x
252+ rig[5, 2 * nw:3 * nw:4], rig[5, 2 * nw + 1:3 * nw:4] = xn - xc, 1.0
253+ for i in range(6): # Gram-Schmidt in the M inner product
254+ for j in range(i):
255+ rig[i] -= (rig[j] @ (self.Mdiag * rig[i])) * rig[j]
256+ rig[i] /= np.sqrt(rig[i] @ (self.Mdiag * rig[i]))
257+ self.rigid = rig
258+ # ---- Rayleigh integral: weights and delays per Gauss point
259+ L = np.array([p["listener_x"], p["listener_y"], p["listener_z"]])
260+ gz = pr.Z(gr)
261+ d = np.sqrt((L[0] - gx) ** 2 + (L[1] - gy) ** 2 + (L[2] - gz) ** 2) # (C, G)
262+ self.ray_c = (p["rho_air"] / (2 * np.pi)) * self.wq[None, :] / d
263+ self.ray_tau = d / p["c_air"]
264+ self.log(f" mesh: {self.n_cells} cells, {self.n_nodes} nodes, {self.n_dof} unknowns, spacing {a * 1e3:g} mm, mass {self.mass_total:.4f} kg, "
265+ f"footprint integrals {[round(s['integral'], 4) for s in self.strikes]}")
266+
267+ # ------------------------------------------------------------------------------------------------------------
268+ def assemble_linear(self):
269+ """Sparse stiffness matrix of the linearized shell (the Hessian of the strain energy at rest), from the same
270+ Gauss-point operators the time stepping uses; symmetric, with the six rigid fields in its null space."""
271+ import scipy.sparse as sp
272+
273+ p = self.p
274+ C, G, nu = self.n_cells, self.G, p["nu"]
275+ B = self.B.reshape(G, 9, 48)
276+ L = np.zeros((C, G, 6, 9)) # 9 derivative quantities -> [exx, eyy, gxy, kxx, kyy, 2kxy] at the reference state
277+ L[:, :, 0, 0] = self.Zx; L[:, :, 0, 5] = 1.0
278+ L[:, :, 1, 1] = self.Zy; L[:, :, 1, 8] = 1.0
279+ L[:, :, 2, 0] = self.Zy; L[:, :, 2, 1] = self.Zx; L[:, :, 2, 6] = 1.0; L[:, :, 2, 7] = 1.0
280+ L[:, :, 3, 2] = 1.0; L[:, :, 4, 3] = 1.0; L[:, :, 5, 4] = 2.0
281+ Q = np.array([[1, nu, 0], [nu, 1, 0], [0, 0, (1 - nu) / 2]])
282+ Dm = np.zeros((C, G, 6, 6))
283+ Dm[:, :, :3, :3] = self.Cm[:, :, None, None] * Q
284+ Dm[:, :, 3:, 3:] = self.Db[:, :, None, None] * Q
285+ LB = np.einsum("cgij,gjk->cgik", L, B)
286+ Ke = np.einsum("g,cgik,cgij,cgjl->ckl", self.wq, LB, Dm, LB, optimize=True)
287+ rows = np.repeat(self.cell_dofs, 48, axis=1).ravel()
288+ cols = np.tile(self.cell_dofs, (1, 48)).ravel()
289+ K = sp.coo_matrix((Ke.ravel(), (rows, cols)), shape=(self.n_dof, self.n_dof)).tocsr()
290+ return 0.5 * (K + K.T)
291+
292+ def build_force(self, torch, dev, compile_: bool):
293+ """The internal force f(q, v) of the assembled shell (Kelvin-Voigt rates included) as a torch function on the
294+ device, and the tensors it shares with the time stepping."""
295+ T64 = lambda x: torch.as_tensor(np.asarray(x, dtype=np.float64), device=dev)
296+ T64i = lambda x: torch.as_tensor(np.asarray(x, dtype=np.int64), device=dev)
297+ p = self.p
298+ C, G = self.n_cells, self.G
299+ cell_dofs = T64i(self.cell_dofs)
300+ BT = T64(self.B.T.copy()) # (48, 81)
301+ Bw = T64(self.B * np.repeat(self.wq, 9)[:, None]) # (81, 48) with quadrature weights folded in
302+ Zx, Zy, Cm, Db = T64(self.Zx), T64(self.Zy), T64(self.Cm), T64(self.Db)
303+ nu, eta = p["nu"], p["eta"]
304+ lin = 0.0 if self.linear else 1.0
305+ n_dof = self.n_dof
306+
307+ def constitutive(D):
308+ Dq, Dv = D[0], D[1]
309+ wx, wy = Dq[..., 0], Dq[..., 1]
310+ ax, ay = Zx + lin * wx, Zy + lin * wy # Z_x + w_x
311+ exx = Dq[..., 5] + Zx * wx + 0.5 * lin * wx * wx
312+ eyy = Dq[..., 8] + Zy * wy + 0.5 * lin * wy * wy
313+ exy = 0.5 * (Dq[..., 6] + Dq[..., 7] + Zx * wy + Zy * wx + lin * wx * wy)
314+ rxx = Dv[..., 5] + ax * Dv[..., 0]
315+ ryy = Dv[..., 8] + ay * Dv[..., 1]
316+ rxy = 0.5 * (Dv[..., 6] + Dv[..., 7] + ax * Dv[..., 1] + ay * Dv[..., 0])
317+ exx = exx + eta * rxx
318+ eyy = eyy + eta * ryy
319+ exy = exy + eta * rxy
320+ Nxx = Cm * (exx + nu * eyy)
321+ Nyy = Cm * (eyy + nu * exx)
322+ Nxy = Cm * (1 - nu) * exy
323+ kxx = Dq[..., 2] + eta * Dv[..., 2]
324+ kyy = Dq[..., 3] + eta * Dv[..., 3]
325+ kxy = Dq[..., 4] + eta * Dv[..., 4]
326+ Mxx = Db * (kxx + nu * kyy)
327+ Myy = Db * (kyy + nu * kxx)
328+ Mxy = Db * (1 - nu) * kxy
329+ return torch.stack([Nxx * ax + Nxy * ay, Nyy * ay + Nxy * ax, Mxx, Myy, 2 * Mxy, Nxx, Nxy, Nxy, Nyy], dim=-1) # (C, G, 9)
330+
331+ constitutive_c = torch.compile(constitutive, dynamic=False) if (compile_ and str(dev).startswith("cuda")) else constitutive
332+
333+ def force(q, v):
334+ qe = q[cell_dofs] # (C, 48)
335+ ve = v[cell_dofs]
336+ D = (torch.stack([qe, ve]) @ BT).view(2, C, G, 9)
337+ S = constitutive_c(D)
338+ fe = S.reshape(C, G * 9) @ Bw # (C, 48)
339+ f = torch.zeros(n_dof, dtype=torch.float64, device=dev)
340+ f.index_add_(0, cell_dofs.view(-1), fe.view(-1))
341+ return f
342+
343+ ops = dict(cell_dofs=cell_dofs, cell_dofs_w=T64i(self.cell_dofs_w), NwT=T64(self.Nw.T.copy()), invM=T64(1.0 / self.Mdiag), T64=T64, T64i=T64i)
344+ return force, ops
345+
346+ def field_interpolator(self, xs, ys):
347+ """Indices (P, 16) and weights (P, 16) giving w at grid points from the w unknowns; mask of points inside."""
348+ XX, YY = np.meshgrid(xs, ys, indexing="xy")
349+ px, py = XX.ravel(), YY.ravel()
350+ a = self.a
351+ ci = np.floor(px / a + 1e-9).astype(int)
352+ cj = np.floor(py / a + 1e-9).astype(int)
353+ lut = {(int(i), int(j)): c for c, (i, j) in enumerate(self.cells)}
354+ idx = np.zeros((len(px), 16), dtype=np.int64)
355+ wts = np.zeros((len(px), 16))
356+ inside = np.zeros(len(px), dtype=bool)
357+ for k in range(len(px)):
358+ c = lut.get((ci[k], cj[k]))
359+ if c is None:
360+ # a point on the outer edge of the last cell row belongs to the cell below/left
361+ for di, dj in [(-1, 0), (0, -1), (-1, -1)]:
362+ c = lut.get((ci[k] + di, cj[k] + dj))
363+ if c is not None:
364+ ci[k] += di
365+ cj[k] += dj
366+ break
367+ if c is None:
368+ continue
369+ inside[k] = True
370+ sx, sy = px[k] / a - ci[k], py[k] / a - cj[k]
371+ idx[k] = self.cell_dofs_w[c]
372+ wts[k] = element_basis([sx], [sy], a)["N"][0]
373+ return idx, wts, inside
374+
375+ # ------------------------------------------------------------------------------------------------------------
376+ def run(self, fs: int, nsamp: int, device: str = "cpu", m: int | None = None, field_grid=None, field_rate=None,
377+ report: float = 60.0, out_dir: Path | None = None, tag: str = "run", use_graph: bool = True, scale: float = 1.0, check: bool = False, compile_: bool = True):
378+ import torch
379+
380+ t_setup = time.perf_counter()
381+ dev = torch.device(device)
382+ force, ops = self.build_force(torch, dev, compile_)
383+ T64, T64i = ops["T64"], ops["T64i"]
384+ p = self.p
385+ C, G = self.n_cells, self.G
386+ cell_dofs, cell_dofs_w, NwT, invM = ops["cell_dofs"], ops["cell_dofs_w"], ops["NwT"], ops["invM"]
387+ sigma0 = p["sigma0"]
388+ ms, kc = p["stick_mass"], p["contact_stiffness"]
389+ v0 = self.v0
390+
391+ zero = torch.zeros(self.n_dof, dtype=torch.float64, device=dev)
392+ torch.manual_seed(0)
393+ x = torch.randn(self.n_dof, dtype=torch.float64, device=dev)
394+ lam = 0.0
395+ saved = self.linear
396+ self.linear = True
397+ force_lin, _ = self.build_force(torch, dev, False)
398+ self.linear = saved
399+ for _ in range(60):
400+ y = invM * force_lin(x, zero)
401+ lam = float(torch.dot(x, y) / torch.dot(x, x))
402+ x = y / torch.linalg.norm(y)
403+ del force_lin
404+ omega_max = math.sqrt(max(lam, 1.0)) * 1.05 # power iteration converges from below; margin
405+ dt_stab = 2.0 / omega_max
406+ if m is None:
407+ m = int(math.ceil(1.0 / (fs * 0.85 * dt_stab)))
408+ dt = 1.0 / (fs * m)
409+ self.log(f" omega_max ~ {omega_max:.3e} rad/s ({omega_max / 2 / math.pi:.0f} Hz), stable dt < {dt_stab:.3e} s; using m = {m}, dt = {dt:.3e} s")
410+ nsteps = nsamp * m
411+
412+ # ---- Rayleigh integral bookkeeping
413+ dsteps = T64(self.ray_tau / dt).view(-1) # (C*G,) delay in steps
414+ rc = T64(self.ray_c).view(-1)
415+ out = torch.zeros(nsamp, dtype=torch.float64, device=dev)
416+ # ---- strikes
417+ nst = len(self.strikes)
418+ phi = T64(np.stack([s["phi"] for s in self.strikes])) # (nst, n_dof)
419+ phi_load = T64(np.stack([s["phi"] - (s["reaction"] if self.no_drift else 0.0) for s in self.strikes])) # force on the shell per unit F, sign as phi
420+ rigM = T64(self.rigid * self.Mdiag[None, :]) # (6, n_dof): <rigid_i, .>_M
421+ rig = T64(self.rigid)
422+
423+ def project_rigid(): # remove the rigid components of the state (the constraint of MODEL.md, enforced per chunk)
424+ for x in (q, v):
425+ x.sub_(rig.T @ (rigM @ x))
426+ t_on = T64([s["t"] for s in self.strikes])
427+ n_on = T64i([int(round(s["t"] / dt)) for s in self.strikes])
428+ z = torch.zeros(nst, dtype=torch.float64, device=dev)
429+ zd = torch.zeros(nst, dtype=torch.float64, device=dev)
430+ active = torch.zeros(nst, dtype=torch.float64, device=dev)
431+ Fst = torch.zeros(nst, dtype=torch.float64, device=dev)
432+ imp = torch.zeros(nst, dtype=torch.float64, device=dev) # contact impulse per stick
433+ dur = torch.zeros(nst, dtype=torch.float64, device=dev) # time in contact
434+ fmax = torch.zeros(nst, dtype=torch.float64, device=dev)
435+ out_direct = torch.zeros(nsamp, dtype=torch.float64, device=dev) if check else None
436+ # ---- state
437+ q = torch.zeros(self.n_dof, dtype=torch.float64, device=dev)
438+ v = torch.zeros(self.n_dof, dtype=torch.float64, device=dev)
439+ acc = torch.zeros(self.n_dof, dtype=torch.float64, device=dev)
440+ nstep = torch.zeros((), dtype=torch.int64, device=dev)
441+ lag_c = 1.0 / 6.0
442+
443+ assert m >= 4, "the retarded-time interpolation assumes at least 4 steps per sample"
444+
445+ def rayleigh(n, ag):
446+ nf = n.to(torch.float64)
447+ kk = torch.ceil((nf - 2 + dsteps) / m)
448+ pos = kk * m - dsteps
449+ fl = torch.floor(pos)
450+ j = nf - fl + 1
451+ fr = pos - fl
452+ wgt = torch.where(j == 0, -fr * (fr - 1) * (fr - 2) * lag_c,
453+ torch.where(j == 1, (fr + 1) * (fr - 1) * (fr - 2) * 0.5,
454+ torch.where(j == 2, -(fr + 1) * fr * (fr - 2) * 0.5,
455+ torch.where(j == 3, (fr + 1) * fr * (fr - 1) * lag_c, torch.zeros_like(fr)))))
456+ valid = (j >= 0) & (j <= 3) & (kk >= 0) & (kk < nsamp)
457+ idx = torch.where(valid, kk, torch.zeros_like(kk)).to(torch.int64)
458+ return idx, torch.where(valid, rc * wgt * ag, torch.zeros_like(ag))
459+
460+ ray_c = torch.compile(rayleigh, dynamic=False) if (compile_ and device.startswith("cuda")) else rayleigh
461+
462+ def step():
463+ n = nstep
464+ # sticks: start the ones whose time has come (at the current surface, with relative speed v0)
465+ wbar = phi @ q
466+ wbar_dot = phi @ v
467+ starting = (n == n_on).to(torch.float64)
468+ z.copy_(z * (1 - starting) + wbar * starting)
469+ zd.copy_(zd * (1 - starting) + (wbar_dot - v0) * starting)
470+ active.copy_(torch.maximum(active, starting))
471+ pen = torch.clamp(wbar - z, min=0.0)
472+ F = active * kc * pen * torch.sqrt(pen)
473+ Fst.copy_(F)
474+ imp.add_(dt * F)
475+ dur.add_(dt * (F > 0).to(torch.float64))
476+ fmax.copy_(torch.maximum(fmax, F))
477+ f = force(q, v)
478+ f = f + F @ phi_load # the shell is pushed down (and, with --no-drift, held rigidly): M a = -f
479+ acc.copy_(-f * invM - 2 * sigma0 * v)
480+ zd.add_(dt * F / ms)
481+ z.add_(dt * zd)
482+ v.add_(dt * acc)
483+ q.add_(dt * v)
484+ # Rayleigh integral: scatter the Gauss-point accelerations onto the retarded output samples
485+ ag = (acc[cell_dofs_w] @ NwT).reshape(-1) # (C*G,)
486+ idx, val = ray_c(n, ag)
487+ out.index_add_(0, idx, val)
488+ if out_direct is not None: # the same integral without retardation, for checking the scatter above
489+ ks = torch.div(n, m, rounding_mode="floor")
490+ on_sample = ((n % m) == 0) & (ks < nsamp)
491+ out_direct.index_add_(0, torch.where(on_sample, ks, torch.zeros_like(ks)).view(1), (torch.sum(rc * ag) * on_sample.to(torch.float64)).view(1))
492+ nstep.add_(1)
493+
494+ # fields
495+ fld = None
496+ if field_grid is not None:
497+ fx, fy, ft = field_grid
498+ idx, wts, inside = self.field_interpolator(fx, fy)
499+ fld = dict(idx=T64i(idx), wts=T64(wts), inside=inside, every=int(round(fs * m / field_rate)), nt=len(ft), ny=len(fy), nx=len(fx), frames=[])
500+
501+ def take_frame():
502+ fr = (q[fld["idx"]] * fld["wts"]).sum(1).cpu().numpy().astype(np.float32)
503+ fr[~fld["inside"]] = np.nan
504+ fld["frames"].append(fr)
505+
506+ def reset_state():
507+ for tsr in (q, v, acc, out, z, zd, active, Fst, imp, dur, fmax, nstep):
508+ tsr.zero_()
509+ if out_direct is not None:
510+ out_direct.zero_()
511+
512+ graph = None
513+ if use_graph and device.startswith("cuda"):
514+ s = torch.cuda.Stream()
515+ s.wait_stream(torch.cuda.current_stream())
516+ with torch.cuda.stream(s):
517+ for _ in range(2):
518+ step()
519+ torch.cuda.current_stream().wait_stream(s)
520+ # undo the warm-up steps
521+ reset_state()
522+ graph = torch.cuda.CUDAGraph()
523+ with torch.cuda.graph(graph):
524+ step()
525+ reset_state()
526+ self.log(f" setup {time.perf_counter() - t_setup:.1f}s; stepping {nsteps} steps ({nsamp / fs:.2f} s) on {device}{' with a CUDA graph' if graph else ''}")
527+
528+ t_start = time.perf_counter()
529+ last_report = t_start
530+ chunk = fld["every"] if fld is not None else 2000
531+ if fld is not None:
532+ take_frame()
533+ stick_stats = lambda: dict(impulse_Ns=[round(float(x), 5) for x in imp.cpu().numpy()], contact_ms=[round(1e3 * float(x), 4) for x in dur.cpu().numpy()],
534+ peak_force_N=[round(float(x), 2) for x in fmax.cpu().numpy()], rebound_mps=[round(float(x), 4) for x in zd.cpu().numpy()])
535+ contact_log = []
536+ Fprev = np.zeros(nst)
537+ n_done = 0
538+ while n_done < nsteps:
539+ k = min(chunk, nsteps - n_done)
540+ for _ in range(k):
541+ if graph is not None:
542+ graph.replay()
543+ else:
544+ step()
545+ n_done += k
546+ if self.no_drift:
547+ project_rigid()
548+ if fld is not None and len(fld["frames"]) < fld["nt"]:
549+ take_frame()
550+ Fc = Fst.cpu().numpy()
551+ for i in range(nst):
552+ if (Fc[i] > 0) != (Fprev[i] > 0):
553+ contact_log.append((round(n_done * dt, 6), i + 1, "on" if Fc[i] > 0 else "off"))
554+ Fprev = Fc
555+ now = time.perf_counter()
556+ if now - last_report >= report or n_done == nsteps:
557+ last_report = now
558+ el = now - t_start
559+ qn = q.cpu().numpy()
560+ if not np.all(np.isfinite(qn)):
561+ raise RuntimeError(f"solution diverged at step {n_done} (t = {n_done * dt:.4f} s)")
562+ w = qn[0 : 4 * self.n_nodes : 4]
563+ sig = out.cpu().numpy()
564+ prog = dict(tag=tag, step=n_done, steps=nsteps, t=round(n_done * dt, 6), duration=nsamp / fs, elapsed_s=round(el, 1),
565+ eta_s=round(el * (nsteps - n_done) / max(n_done, 1), 1), steps_per_s=round(n_done / el, 1),
566+ max_abs_w_m=float(np.abs(w).max()), rms_w_m=float(np.sqrt(np.mean(w**2))),
567+ stick_z_m=[round(float(x), 6) for x in z.cpu().numpy()], stick_v_mps=[round(float(x), 4) for x in zd.cpu().numpy()],
568+ contact_force_N=[round(float(x), 2) for x in Fc], contacts=contact_log[-40:], **stick_stats(),
569+ peak_pa=float(np.abs(sig).max()), m=m, dt=dt, device=device)
570+ self.log(f" t={prog['t']:.3f}s step {n_done}/{nsteps} elapsed {el:.0f}s eta {prog['eta_s']:.0f}s ({prog['steps_per_s']:.0f} steps/s) "
571+ f"max|w|={prog['max_abs_w_m']:.3e} peak={prog['peak_pa']:.3f} Pa sticks z={prog['stick_z_m']} F={prog['contact_force_N']}")
572+ if out_dir is not None:
573+ with open(out_dir / f"{tag}.progress.json", "w") as fh:
574+ json.dump(prog, fh, indent=1)
575+ wavfile.write(out_dir / f"{tag}.partial.wav", fs, (sig / scale).astype(np.float32))
576+ elapsed = time.perf_counter() - t_start
577+ sig = out.cpu().numpy()
578+ info = dict(m=m, dt=dt, omega_max=omega_max, steps=nsteps, elapsed=elapsed, contacts=contact_log, sticks=stick_stats(),
579+ out_direct=out_direct.cpu().numpy() if out_direct is not None else None)
580+ fields = np.stack(fld["frames"]).reshape(len(fld["frames"]), fld["ny"], fld["nx"]) if fld is not None else None
581+ return sig, fields, info
582+
583+
584+# ----------------------------------------------------------------------------------------------------------------
585+
586+def hardware_string(device: str) -> str:
587+ cpu = "unknown CPU"
588+ try:
589+ for line in open("/proc/cpuinfo"):
590+ if line.startswith("model name"):
591+ cpu = line.split(":", 1)[1].strip()
592+ break
593+ except OSError:
594+ pass
595+ import os
596+
597+ s = f"{cpu} ({os.cpu_count()} threads)"
598+ if device.startswith("cuda"):
599+ import torch
600+
601+ s += f"; time stepping on {torch.cuda.get_device_name(0)} (PyTorch {torch.__version__})"
602+ return s
603+
604+
605+def main():
606+ ch = json.load(open(sys.argv[1]))
607+ out_dir = Path(sys.argv[2])
608+ args = sys.argv[3:]
609+
610+ def take(flag, conv, default=None):
611+ nonlocal args
612+ if flag in args:
613+ i = args.index(flag)
614+ val = conv(args[i + 1])
615+ args = args[:i] + args[i + 2 :]
616+ return val
617+ return default
618+
619+ spacing = take("--spacing", float, 0.008)
620+ seconds = take("--seconds", float)
621+ device = take("--device", str, "cpu")
622+ report = take("--report", float, 60.0)
623+ speed = take("--speed", float)
624+ m_over = take("--m", int)
625+ tag = take("--tag", str)
626+ c_air = take("--c-air", float)
627+ check = "--check" in args
628+ no_drift = "--no-drift" in args
629+ inplane_scale = take("--inplane-scale", float, 16.0)
630+ compile_ = "--no-compile" not in args
631+ linear = "--linear" in args
632+ fields = "--fields" in args
633+ use_graph = "--no-graph" not in args
634+ out = ch["output"]
635+ fs = out["sample_rate_hz"]
636+ nsamp = out["num_samples"] if seconds is None else int(round(fs * seconds))
637+ tag = tag or f"fem-g{spacing * 1e3:g}"
638+ log = lambda s: print(s, file=sys.stderr, flush=True)
639+ out_dir.mkdir(parents=True, exist_ok=True)
640+ t0 = time.perf_counter()
641+ log(f"{tag}: {ch['challenge_id']}, spacing {spacing * 1e3:g} mm, {nsamp / fs:.2f} s, device {device}{' (linear)' if linear else ''}")
642+ model = ShellFEM(ch, spacing, log=log, linear=linear, speed=speed, c_air=c_air, no_drift=no_drift, inplane_scale=inplane_scale)
643+ scale = out.get("signal_scale") or 1.0
644+ grid = None
645+ if fields:
646+ fg = ch["fields"]
647+ nt, ny, nx = fg["shape"]
648+ R = ch["parameters"]["disk_radius"]
649+ grid = (np.linspace(-R, R, nx), np.linspace(-R, R, ny), np.arange(nt) / fg["frame_rate_hz"])
650+ sig, fld, info = model.run(fs, nsamp, device=device, m=m_over, report=report, out_dir=out_dir, tag=tag, use_graph=use_graph, scale=scale, check=check, compile_=compile_,
651+ field_grid=grid, field_rate=ch["fields"]["frame_rate_hz"] if grid else None)
652+ total = time.perf_counter() - t0
653+ wavfile.write(out_dir / f"{tag}.wav", fs, (sig / scale).astype(np.float32))
654+ meta = {
655+ "solution_id": tag,
656+ "output": f"{out_dir.name}/{tag}.wav",
657+ "method": "bicubic Hermite (Bogner-Fox-Schmit) finite elements for u, v, w on a Cartesian grid masked to the annulus, lumped mass, explicit leapfrog, Hertz stick, Rayleigh integral accumulated at the retarded times",
658+ "parameters": {"grid_spacing_m": spacing, "cells": model.n_cells, "n_dof": model.n_dof, "gauss_points_per_cell": model.G,
659+ "oversampling": info["m"], "dt_s": info["dt"], "omega_max_rad_s": info["omega_max"], "inplane_derivative_mass_scale": inplane_scale},
660+ "timing_seconds": round(info["elapsed"], 3),
661+ "hardware": hardware_string(device),
662+ "peak_si": float(np.max(np.abs(sig))),
663+ "sticks": info["sticks"],
664+ **({"stick_speed_override": speed} if speed is not None else {}),
665+ **({"linear": True} if linear else {}),
666+ **({"no_drift": True} if no_drift else {}),
667+ }
668+ if fld is not None:
669+ import h5py
670+
671+ fg = ch["fields"]
672+ with h5py.File(out_dir / f"{tag}.h5", "w") as fh:
673+ fh.attrs["challenge_id"] = ch["challenge_id"]
674+ fh.attrs["solution_id"] = tag
675+ fh.attrs["quantity"] = fg["quantity"]
676+ fh.attrs["units"] = fg["units"]
677+ fh.create_dataset("x", data=grid[0])
678+ fh.create_dataset("y", data=grid[1])
679+ fh.create_dataset("t", data=grid[2][: fld.shape[0]])
680+ fh.create_dataset("u", data=fld, compression="gzip", compression_opts=4, chunks=(8, fld.shape[1], fld.shape[2]))
681+ meta["fields"] = {"path": f"{out_dir.name}/{tag}.h5", "format": "hdf5"}
682+ with open(out_dir / f"{tag}.meta.json", "w") as fh:
683+ json.dump(meta, fh, indent=1)
684+ if info["out_direct"] is not None:
685+ d = info["out_direct"]
686+ log(f" check: max|out - direct| = {np.abs(sig - d).max():.3e} Pa against max|direct| = {np.abs(d).max():.3e} Pa (equal only when c_air is huge)")
687+ np.save(out_dir / f"{tag}.direct.npy", d)
688+ log(f"{tag}: done in {total:.0f}s, peak {meta['peak_si']:.4f} Pa, sticks {info['sticks']}")
689+
690+
691+if __name__ == "__main__":
692+ main()
solvers/hermite-fem/fetch_progress.shadded+15−0View file
@@ -0,0 +1,15 @@
1+#!/bin/bash
2+# Copy the progress files and partial wavs of the running chain from the office workstation into this checkout,
3+# and print a one-line status per run. Usage: bash fetch_progress.sh
4+here="$(cd "$(dirname "$0")/../.." && pwd)"
5+remote=office:/home/magland/src/vault1/challenges/cymbal2-solutions
6+for c in cymbal2-flat-v1 cymbal2-crash-v1; do
7+ mkdir -p "$here/solutions/$c/hermite-fem"
8+ scp -q "$remote/solutions/$c/hermite-fem/*.progress.json" "$remote/solutions/$c/hermite-fem/*.partial.wav" "$here/solutions/$c/hermite-fem/" 2>/dev/null
9+ for j in "$here"/solutions/$c/hermite-fem/*.progress.json; do
10+ [ -f "$j" ] && python3 -c "
11+import json,sys; d=json.load(open('$j'))
12+print(f\"$c {d['tag']}: t = {d['t']:.2f} / {d['duration']:.0f} s, {d['elapsed_s']/60:.0f} min elapsed, {d['eta_s']/60:.0f} min left, peak {d['peak_pa']:.1f} Pa, max|w| {d['max_abs_w_m']*1e3:.1f} mm, sticks: impulse {d['impulse_Ns']} N s, contact {d['contact_ms']} ms, peak force {d['peak_force_N']} N, rebound {d['rebound_mps']} m/s\")"
13+ done
14+done
15+exit 0
solvers/hermite-fem/make_submission.pyadded+42−0View file
@@ -0,0 +1,42 @@
1+"""Write solutions/<challenge>/hermite-fem.json from the <tag>.meta.json files the solver leaves beside its wavs.
2+Usage: python make_submission.py <challenge-id> <submission-id> <tag-in-order-of-increasing-accuracy> ..."""
3+import json
4+import sys
5+from pathlib import Path
6+
7+ROOT = Path(__file__).resolve().parents[2]
8+DESCRIPTIONS = {
9+ "hermite-fem": "Displacement-based finite elements on a Cartesian grid masked to the annulus: bicubic Hermite (Bogner-Fox-Schmit) elements for u, v and w, lumped mass, explicit leapfrog at the in-plane stability limit, Hertz sticks as point masses stepped with the shell, Rayleigh integral accumulated at the retarded times, rigid motion held by the rigid reaction and projection. Chain ordered by grid spacing (8 mm, 5 mm).",
10+ "hermite-modal": "The same finite element system as hermite-fem, solved in its eigenmodes below a cutoff frequency: exact propagation of the damped linear part, the nonlinear force evaluated at the Gauss points of the mesh and projected on the modes, Hertz sticks, Rayleigh integral by a retarded weight table. The rigid fields are excluded from the basis. Chain ordered by grid spacing, cutoff and step.",
11+}
12+cid = sys.argv[1]
13+sid = sys.argv[2]
14+tags = sys.argv[3:]
15+out_dir = ROOT / "solutions" / cid / sid
16+sols = []
17+for k, tag in enumerate(tags):
18+ m = json.load(open(out_dir / f"{tag}.meta.json"))
19+ m.pop("peak_si", None)
20+ sticks = m.pop("sticks", None)
21+ m["notes"] = (f"Sticks: impulse {sticks['impulse_Ns']} N s, contact {sticks['contact_ms']} ms, peak force {sticks['peak_force_N']} N, "
22+ f"rebound {sticks['rebound_mps']} m/s (bell, mid, rim)." if sticks else "")
23+ m["more_accurate_than"] = tags[:k]
24+ if "output" in m:
25+ m["output"] = f"{sid}/{Path(m['output']).name}"
26+ if "fields" in m:
27+ m["fields"]["path"] = f"{sid}/{Path(m['fields']['path']).name}"
28+ sols.append(m)
29+sub = {
30+ "schema_version": 1,
31+ "submission_id": sid,
32+ "challenge_id": cid,
33+ "submitter": {"name": "Jeremy Magland", "url": "https://vault1.magland.org/challenges/cymbal2-solutions"},
34+ "date": "2026-09-18",
35+ "description": DESCRIPTIONS[sid],
36+ "code_url": "https://vault1.magland.org/challenges/cymbal2-solutions/tree/main/solvers/hermite-fem",
37+ "viz_url": "https://vault1.magland.org/challenges/cymbal2-solutions/site/",
38+ "solutions": sols,
39+}
40+with open(ROOT / "solutions" / cid / "hermite-fem.json", "w") as fh:
41+ json.dump(sub, fh, indent=1)
42+print(ROOT / "solutions" / cid / "hermite-fem.json")
solvers/hermite-fem/modal_shell.pyadded+370−0View file
@@ -0,0 +1,370 @@
1+"""Modal solution of the cymbal2 shallow-shell model on the same discretization as fem_shell.py.
2+
3+The semi-discrete system is the one fem_shell.py builds: bicubic Hermite elements for u, v, w on the Cartesian grid,
4+3 x 3 Gauss points, lumped mass (with the same in-plane derivative-mass scaling). Here the linear part of that system
5+is diagonalized once, K phi = omega^2 M phi, and only the modes below a cutoff frequency f_max are kept; the six rigid
6+fields are excluded exactly, by a penalty in the eigenproblem that lifts them out of range (the rigid constraint of
7+MODEL.md). The state is the vector of modal amplitudes a(t), and the equations are
8+
9+ a_k'' + 2 alpha_k a_k' + omega_k^2 a_k = F_k(a, a'), alpha_k = sigma0 + eta omega_k^2 / 2,
10+
11+where F is the modal projection of the nonlinear remainder of the internal force, evaluated at the Gauss points of the
12+mesh by the same routine as the finite element solver (so the nonlinear terms are not truncated to any tensor of modal
13+couplings), plus the stick loads. The time stepping is the two-step scheme of the v1 modal solver: the damped linear
14+oscillators are propagated exactly over each step and the force is applied with a static-exact weight, which is
15+second-order accurate and free of dispersion in the linear limit, so the step can be a small multiple of the output
16+sample. The sticks obey the Hertz law with the footprint-averaged displacement as in fem_shell.py. The Rayleigh
17+integral is a per-step weight table folded with a ring buffer of modal accelerations (cubic Lagrange interpolation in
18+time at the retarded time of every Gauss point).
19+
20+The two solvers therefore solve the same discrete equations and differ only in the modal truncation and the time
21+integration; their difference measures the effect of the modes above f_max and of the step.
22+
23+Usage:
24+ python modal_shell.py <challenge.json> <out-dir> [--spacing 0.008] [--fmax 20000] [--m 4] [--seconds S] [--device cuda]
25+ [--fields] [--report 60] [--tag NAME] [--inplane-scale 16]
26+"""
27+from __future__ import annotations
28+
29+import json
30+import math
31+import sys
32+import time
33+from pathlib import Path
34+
35+import numpy as np
36+import scipy.sparse as sp
37+import scipy.sparse.linalg as spla
38+from scipy.io import wavfile
39+
40+from fem_shell import ShellFEM, hardware_string
41+
42+
43+def sliced_modes(K, Mdiag, rigid, f_max: float, k_slice: int = 250, log=print, penalty: float = 1e12):
44+ """Every eigenpair of K phi = lambda M phi with lambda <= (2 pi f_max)^2, except the six rigid fields, which are
45+ lifted to the eigenvalue `penalty` by the rank-6 term M R R^T M (R the M-orthonormal rigid fields) and thereby
46+ made exactly orthogonal to every returned mode. Shift-invert ARPACK in slices of k_slice eigenvalues, each slice
47+ factorizing K - sigma M once (the low-rank term is handled with the Woodbury identity)."""
48+ n = K.shape[0]
49+ M = sp.diags(Mdiag).tocsc()
50+ U = (Mdiag[:, None] * rigid.T) * math.sqrt(penalty) # (n, 6): penalty term is U U^T
51+ Kop = spla.LinearOperator((n, n), matvec=lambda x: K @ x + U @ (U.T @ x), dtype=np.float64)
52+ target = (2 * math.pi * f_max) ** 2
53+ vals_all, vecs_all = [], []
54+ lo = -np.inf
55+ sigma = -1e3
56+ width = None
57+ tries = 0
58+ while True:
59+ t0 = time.perf_counter()
60+ lu = spla.splu((K - sigma * M).tocsc())
61+ AiU = lu.solve(U) # (n, 6)
62+ S = np.linalg.inv(np.eye(6) + U.T @ AiU)
63+
64+ def opinv(x, lu=lu, AiU=AiU, S=S):
65+ y = lu.solve(x)
66+ return y - AiU @ (S @ (U.T @ y))
67+
68+ OPinv = spla.LinearOperator((n, n), matvec=opinv, dtype=np.float64)
69+ vals, vecs = spla.eigsh(Kop, k=k_slice, M=M, sigma=sigma, which="LM", OPinv=OPinv, tol=1e-9)
70+ order = np.argsort(vals)
71+ vals, vecs = vals[order], vecs[:, order]
72+ tries += 1
73+ if np.isfinite(lo) and vals[0] > lo * (1 + 1e-9) + 1e-3:
74+ if tries > 40:
75+ raise RuntimeError("mode slicing did not close the gap")
76+ sigma = 0.5 * (sigma + lo) # the slice did not reach back to the covered range: step back
77+ log(f" slice at {math.sqrt(max(sigma, 0)) / 2 / math.pi:.0f} Hz left a gap, retrying")
78+ continue
79+ keep = vals > lo if np.isfinite(lo) else np.ones(len(vals), dtype=bool)
80+ vals_all.append(vals[keep])
81+ vecs_all.append(vecs[:, keep])
82+ lo = vals[-1]
83+ width = vals[-1] - vals[0]
84+ log(f" slice: {int(keep.sum())} new modes up to {math.sqrt(max(lo, 0)) / 2 / math.pi:.0f} Hz in {time.perf_counter() - t0:.1f}s")
85+ if lo >= target:
86+ break
87+ sigma = lo + 0.45 * width
88+ vals = np.concatenate(vals_all)
89+ vecs = np.concatenate(vecs_all, axis=1)
90+ sel = (vals <= target) & (vals < 0.5 * penalty)
91+ vals, vecs = vals[sel], vecs[:, sel]
92+ # rigid content check and clean-up of the near-null modes of the flat geometry (exact zeros the penalty removed)
93+ rc = np.abs(rigid @ (Mdiag[:, None] * vecs)).max()
94+ log(f" {len(vals)} modes below {f_max:g} Hz; largest rigid component {rc:.1e}")
95+ return np.sqrt(np.maximum(vals, 0.0)), vecs
96+
97+
98+class ModalShell:
99+ def __init__(self, fem: ShellFEM, f_max: float, log=print):
100+ self.fem = fem
101+ self.log = log
102+ t0 = time.perf_counter()
103+ K = fem.assemble_linear()
104+ self.log(f" linear stiffness assembled ({K.nnz} nonzeros) in {time.perf_counter() - t0:.1f}s")
105+ t0 = time.perf_counter()
106+ self.omega, self.Phi = sliced_modes(K, fem.Mdiag, fem.rigid, f_max, log=log)
107+ self.n_modes = len(self.omega)
108+ self.log(f" {self.n_modes} modes in {time.perf_counter() - t0:.1f}s: f = {np.round(self.omega[:5] / 2 / np.pi, 2)} ... {self.omega[-1] / 2 / np.pi:.0f} Hz")
109+ p = fem.p
110+ self.alpha = p["sigma0"] + 0.5 * p["eta"] * self.omega**2
111+
112+ def run(self, fs: int, nsamp: int, m: int, device: str = "cpu", field_grid=None, field_rate=None, report: float = 60.0,
113+ out_dir: Path | None = None, tag: str = "modal", scale: float = 1.0, compile_: bool = True, use_graph: bool = True):
114+ import torch
115+
116+ fem, p = self.fem, self.fem.p
117+ t_setup = time.perf_counter()
118+ dev = torch.device(device)
119+ force, ops = fem.build_force(torch, dev, compile_)
120+ T64, T64i = ops["T64"], ops["T64i"]
121+ N = self.n_modes
122+ dt = 1.0 / (fs * m)
123+ nsteps = nsamp * m
124+ Phi = T64(self.Phi) # (n_dof, N)
125+ PhiT = Phi.T.contiguous()
126+ w2 = T64(self.omega**2)
127+ al = T64(self.alpha)
128+ eta, sigma0 = p["eta"], p["sigma0"]
129+ ms, kc, v0 = p["stick_mass"], p["contact_stiffness"], fem.v0
130+ # exact two-step recurrence of the damped linear oscillators
131+ disc = np.sqrt(self.alpha**2 - self.omega**2 + 0j)
132+ E1, E2 = np.exp((-self.alpha + disc) * dt), np.exp((-self.alpha - disc) * dt)
133+ a1 = T64((E1 + E2).real)
134+ a2 = T64(-(E1 * E2).real)
135+ b = T64((1.0 - (E1 + E2).real + (E1 * E2).real) / self.omega**2)
136+ # Rayleigh integral: per-step weight table G[d, k] from the modes at the Gauss points
137+ tau = fem.ray_tau.reshape(-1)
138+ cq = fem.ray_c.reshape(-1)
139+ tau_steps = tau / dt
140+ d0 = np.floor(tau_steps).astype(np.int64)
141+ fr = tau_steps - d0
142+ lag = np.stack([-fr * (fr - 1.0) * (fr - 2.0) / 6.0, (fr + 1.0) * (fr - 1.0) * (fr - 2.0) / 2.0,
143+ -(fr + 1.0) * fr * (fr - 2.0) / 2.0, (fr + 1.0) * fr * (fr - 1.0) / 6.0], axis=1)
144+ dmax = int(d0.max()) + 2
145+ L = dmax + 1
146+ Nw = T64(fem.Nw) # (G, 16)
147+ Psi = torch.einsum("gi,cin->cgn", Nw, Phi[ops["cell_dofs_w"]]).reshape(-1, N) # (C*G, N) w of the modes at the Gauss points
148+ Gt = torch.zeros((L, N), dtype=torch.float64, device=dev)
149+ for j, off in enumerate([-1, 0, 1, 2]):
150+ Gt.index_add_(0, T64i(d0 + off), T64(cq * lag[:, j])[:, None] * Psi)
151+ del Psi
152+ # sticks
153+ nst = len(fem.strikes)
154+ gains = T64(np.stack([s["phi"] for s in fem.strikes])) @ Phi # (nst, N): footprint average of the modes
155+ n_on = T64i([int(round(s["t"] / dt)) for s in fem.strikes])
156+ z = torch.zeros(nst, dtype=torch.float64, device=dev)
157+ zd = torch.zeros(nst, dtype=torch.float64, device=dev)
158+ active = torch.zeros(nst, dtype=torch.float64, device=dev)
159+ Fst = torch.zeros(nst, dtype=torch.float64, device=dev)
160+ imp = torch.zeros(nst, dtype=torch.float64, device=dev)
161+ dur = torch.zeros(nst, dtype=torch.float64, device=dev)
162+ fmax_t = torch.zeros(nst, dtype=torch.float64, device=dev)
163+ # state
164+ a = torch.zeros(N, dtype=torch.float64, device=dev)
165+ a_prev = torch.zeros(N, dtype=torch.float64, device=dev)
166+ q = torch.zeros(fem.n_dof, dtype=torch.float64, device=dev)
167+ q_prev = torch.zeros(fem.n_dof, dtype=torch.float64, device=dev)
168+ ring = torch.zeros((L, N), dtype=torch.float64, device=dev)
169+ out = torch.zeros(nsamp, dtype=torch.float64, device=dev)
170+ nstep = torch.zeros((), dtype=torch.int64, device=dev)
171+ ar = torch.arange(L, device=dev)
172+ fld = None
173+ if field_grid is not None:
174+ fx, fy, ft = field_grid
175+ idx, wts, inside = fem.field_interpolator(fx, fy)
176+ Wg = (T64(wts)[:, :, None] * Phi[T64i(idx)]).sum(1) # (P, N)
177+ fld = dict(Wg=Wg, inside=inside, every=int(round(fs / field_rate)), nt=len(ft), ny=len(fy), nx=len(fx), frames=[])
178+
179+ def step():
180+ n = nstep
181+ q.copy_(Phi @ a)
182+ v = (q - q_prev) / dt
183+ adot = (a - a_prev) / dt
184+ # sticks
185+ wbar = gains @ a
186+ wbar_dot = gains @ adot
187+ starting = (n == n_on).to(torch.float64)
188+ z.copy_(z * (1 - starting) + wbar * starting)
189+ zd.copy_(zd * (1 - starting) + (wbar_dot - v0) * starting)
190+ active.copy_(torch.maximum(active, starting))
191+ pen = torch.clamp(wbar - z, min=0.0)
192+ F = active * kc * pen * torch.sqrt(pen)
193+ Fst.copy_(F)
194+ imp.add_(dt * F)
195+ dur.add_(dt * (F > 0).to(torch.float64))
196+ fmax_t.copy_(torch.maximum(fmax_t, F))
197+ # modal force: nonlinear remainder of the internal force plus the sticks
198+ f = force(q, v)
199+ Fm = -(PhiT @ f) + w2 * a + eta * w2 * adot - F @ gains
200+ a_next = a1 * a + a2 * a_prev + b * Fm
201+ add = Fm - w2 * a - al * ((a_next - a_prev) / dt)
202+ ring[n % L] = add
203+ zd.add_(dt * F / ms)
204+ z.add_(dt * zd)
205+ q_prev.copy_(q)
206+ a_prev.copy_(a)
207+ a.copy_(a_next)
208+ nstep.add_(1)
209+
210+ def sample(k):
211+ # sample k is at step k m: the first step below computes and stores the acceleration there, the output sample
212+ # folds the table with the ring (delays counted back from that step), then the remaining steps of the sample
213+ step()
214+ idx = (nstep - 1 - ar) % L
215+ out[k] = torch.sum(Gt * ring[idx])
216+ for _ in range(m - 1):
217+ step()
218+
219+ graph = None
220+ k_t = torch.zeros((), dtype=torch.int64, device=dev)
221+
222+ def sample_graph():
223+ sample(k_t)
224+
225+ def reset():
226+ for t in (z, zd, active, Fst, imp, dur, fmax_t, a, a_prev, q, q_prev, ring, out, nstep, k_t):
227+ t.zero_()
228+
229+ if use_graph and device.startswith("cuda"):
230+ s = torch.cuda.Stream()
231+ s.wait_stream(torch.cuda.current_stream())
232+ with torch.cuda.stream(s):
233+ for _ in range(2):
234+ sample_graph()
235+ torch.cuda.current_stream().wait_stream(s)
236+ reset()
237+ graph = torch.cuda.CUDAGraph()
238+ with torch.cuda.graph(graph):
239+ sample_graph()
240+ reset()
241+ self.log(f" setup {time.perf_counter() - t_setup:.1f}s; stepping {nsteps} steps ({nsamp / fs:.2f} s, m = {m}) on {device}{' with a CUDA graph' if graph else ''}")
242+
243+ def take_frame():
244+ fr_ = (fld["Wg"] @ a).cpu().numpy().astype(np.float32)
245+ fr_[~fld["inside"]] = np.nan
246+ fld["frames"].append(fr_)
247+
248+ stick_stats = lambda: dict(impulse_Ns=[round(float(x), 5) for x in imp.cpu().numpy()], contact_ms=[round(1e3 * float(x), 4) for x in dur.cpu().numpy()],
249+ peak_force_N=[round(float(x), 2) for x in fmax_t.cpu().numpy()], rebound_mps=[round(float(x), 4) for x in zd.cpu().numpy()])
250+ t_start = time.perf_counter()
251+ last_report = t_start
252+ if fld is not None:
253+ take_frame()
254+ for k in range(nsamp):
255+ if graph is not None:
256+ k_t.fill_(k)
257+ graph.replay()
258+ else:
259+ sample(k)
260+ if fld is not None and (k + 1) % fld["every"] == 0 and len(fld["frames"]) < fld["nt"]:
261+ take_frame()
262+ now = time.perf_counter()
263+ if (now - last_report >= report or k == nsamp - 1) and k > 0:
264+ last_report = now
265+ el = now - t_start
266+ an = a.cpu().numpy()
267+ if not np.all(np.isfinite(an)):
268+ raise RuntimeError(f"solution diverged at sample {k} (t = {k / fs:.4f} s); increase m")
269+ sig = out.cpu().numpy()
270+ prog = dict(tag=tag, step=(k + 1) * m, steps=nsteps, t=round((k + 1) / fs, 6), duration=nsamp / fs, elapsed_s=round(el, 1),
271+ eta_s=round(el * (nsamp - k - 1) / (k + 1), 1), steps_per_s=round((k + 1) * m / el, 1), max_abs_a=float(np.abs(an).max()),
272+ stick_z_m=[round(float(x), 6) for x in z.cpu().numpy()], stick_v_mps=[round(float(x), 4) for x in zd.cpu().numpy()],
273+ contact_force_N=[round(float(x), 2) for x in Fst.cpu().numpy()], **stick_stats(), peak_pa=float(np.abs(sig).max()), m=m, dt=dt, device=device, n_modes=N)
274+ self.log(f" t={prog['t']:.3f}s elapsed {el:.0f}s eta {prog['eta_s']:.0f}s ({prog['steps_per_s']:.0f} steps/s) max|a|={prog['max_abs_a']:.3e} "
275+ f"peak={prog['peak_pa']:.3f} Pa sticks z={prog['stick_z_m']} F={prog['contact_force_N']}")
276+ if out_dir is not None:
277+ with open(out_dir / f"{tag}.progress.json", "w") as fh:
278+ json.dump(prog, fh, indent=1)
279+ wavfile.write(out_dir / f"{tag}.partial.wav", fs, (sig / scale).astype(np.float32))
280+ elapsed = time.perf_counter() - t_start
281+ sig = out.cpu().numpy()
282+ if not np.all(np.isfinite(sig)):
283+ raise RuntimeError("output is not finite; increase m")
284+ fields = np.stack(fld["frames"]).reshape(len(fld["frames"]), fld["ny"], fld["nx"]) if fld is not None else None
285+ return sig, fields, dict(m=m, dt=dt, elapsed=elapsed, sticks=stick_stats())
286+
287+
288+def main():
289+ ch = json.load(open(sys.argv[1]))
290+ out_dir = Path(sys.argv[2])
291+ args = sys.argv[3:]
292+
293+ def take(flag, conv, default=None):
294+ nonlocal args
295+ if flag in args:
296+ i = args.index(flag)
297+ val = conv(args[i + 1])
298+ args = args[:i] + args[i + 2 :]
299+ return val
300+ return default
301+
302+ spacing = take("--spacing", float, 0.008)
303+ f_max = take("--fmax", float, 20000.0)
304+ m = take("--m", int, 4)
305+ seconds = take("--seconds", float)
306+ device = take("--device", str, "cpu")
307+ report = take("--report", float, 60.0)
308+ speed = take("--speed", float)
309+ tag = take("--tag", str)
310+ inplane_scale = take("--inplane-scale", float, 16.0)
311+ fields = "--fields" in args
312+ compile_ = "--no-compile" not in args
313+ use_graph = "--no-graph" not in args
314+ out = ch["output"]
315+ fs = out["sample_rate_hz"]
316+ nsamp = out["num_samples"] if seconds is None else int(round(fs * seconds))
317+ tag = tag or f"modal-g{spacing * 1e3:g}-f{f_max / 1e3:g}k-m{m}"
318+ log = lambda s: print(s, file=sys.stderr, flush=True)
319+ out_dir.mkdir(parents=True, exist_ok=True)
320+ t0 = time.perf_counter()
321+ log(f"{tag}: {ch['challenge_id']}, spacing {spacing * 1e3:g} mm, modes to {f_max:g} Hz, m = {m}, {nsamp / fs:.2f} s, device {device}")
322+ fem = ShellFEM(ch, spacing, log=log, speed=speed, no_drift=True, inplane_scale=inplane_scale)
323+ model = ModalShell(fem, f_max, log=log)
324+ t1 = time.perf_counter()
325+ scale = out.get("signal_scale") or 1.0
326+ grid = None
327+ if fields:
328+ fg = ch["fields"]
329+ nt, ny, nx = fg["shape"]
330+ R = ch["parameters"]["disk_radius"]
331+ grid = (np.linspace(-R, R, nx), np.linspace(-R, R, ny), np.arange(nt) / fg["frame_rate_hz"])
332+ sig, fld, info = model.run(fs, nsamp, m, device=device, field_grid=grid, field_rate=ch["fields"]["frame_rate_hz"] if grid else None,
333+ report=report, out_dir=out_dir, tag=tag, scale=scale, compile_=compile_, use_graph=use_graph)
334+ total = time.perf_counter() - t0
335+ wavfile.write(out_dir / f"{tag}.wav", fs, (sig / scale).astype(np.float32))
336+ meta = {
337+ "solution_id": tag,
338+ "output": f"{out_dir.name}/{tag}.wav",
339+ "method": "modes of the bicubic Hermite finite element system on a Cartesian grid masked to the annulus (lumped mass), truncated at a cutoff frequency, "
340+ "exact linear propagation with the nonlinear force evaluated at the Gauss points, Hertz stick, Rayleigh integral by a retarded weight table",
341+ "parameters": {"grid_spacing_m": spacing, "cells": fem.n_cells, "n_dof": fem.n_dof, "n_modes": model.n_modes, "mode_cutoff_hz": f_max,
342+ "highest_mode_hz": round(float(model.omega.max() / 2 / np.pi), 1), "oversampling": m, "dt_s": info["dt"], "inplane_derivative_mass_scale": inplane_scale},
343+ "timing_seconds": round(info["elapsed"], 3),
344+ "setup_seconds": round(t1 - t0, 3),
345+ "hardware": hardware_string(device),
346+ "peak_si": float(np.max(np.abs(sig))),
347+ "sticks": info["sticks"],
348+ **({"stick_speed_override": speed} if speed is not None else {}),
349+ }
350+ if fld is not None:
351+ import h5py
352+
353+ fg = ch["fields"]
354+ with h5py.File(out_dir / f"{tag}.h5", "w") as fh:
355+ fh.attrs["challenge_id"] = ch["challenge_id"]
356+ fh.attrs["solution_id"] = tag
357+ fh.attrs["quantity"] = fg["quantity"]
358+ fh.attrs["units"] = fg["units"]
359+ fh.create_dataset("x", data=grid[0])
360+ fh.create_dataset("y", data=grid[1])
361+ fh.create_dataset("t", data=grid[2][: fld.shape[0]])
362+ fh.create_dataset("u", data=fld, compression="gzip", compression_opts=4, chunks=(8, fld.shape[1], fld.shape[2]))
363+ meta["fields"] = {"path": f"{out_dir.name}/{tag}.h5", "format": "hdf5"}
364+ with open(out_dir / f"{tag}.meta.json", "w") as fh:
365+ json.dump(meta, fh, indent=1)
366+ log(f"{tag}: done in {total:.0f}s, peak {meta['peak_si']:.4f} Pa, sticks {info['sticks']}")
367+
368+
369+if __name__ == "__main__":
370+ main()
solvers/hermite-fem/run_chain.shadded+18−0View file
@@ -0,0 +1,18 @@
1+#!/bin/bash
2+# Produce the hermite-fem solutions of both cymbal2 challenges, crash then flat, 8 mm then 5 mm (a quarter of an hour per 8 mm run on the office GPU,
3+# under an hour at 5 mm). Run locally on the GPU workstation, never on the vault runner:
4+# nohup bash run_chain.sh <cymbal2-repo> > chain.log 2>&1 < /dev/null & disown
5+# Progress: solutions/<challenge>/hermite-fem/<tag>.progress.json and <tag>.partial.wav, rewritten every minute.
6+set -e
7+here="$(cd "$(dirname "$0")" && pwd)"
8+spec="$1"
9+py="${PYTHON:-/home/magland/src/vault1/challenges/cymbal-env/bin/python}"
10+run() { # challenge spacing_m
11+ out="$here/../../solutions/$1/hermite-fem"
12+ mkdir -p "$out"
13+ "$py" "$here/fem_shell.py" "$spec/challenges/$1/challenge.json" "$out" --spacing "$2" --device cuda --fields --report 60 --no-drift
14+}
15+run cymbal2-crash-v1 0.008
16+run cymbal2-crash-v1 0.005
17+run cymbal2-flat-v1 0.008
18+run cymbal2-flat-v1 0.005
web/build.pyadded+75−0View file
@@ -0,0 +1,75 @@
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+ # experiments: wav files with a JSON note beside them, outside the challenges (strike schedules and the like)
44+ manifest["experiments"] = []
45+ if (ROOT / "experiments").exists():
46+ shutil.copytree(ROOT / "experiments", out / "experiments")
47+ for j in sorted((ROOT / "experiments").glob("*.json")):
48+ meta = json.loads(j.read_text())
49+ wav = j.with_suffix(".wav")
50+ if wav.exists():
51+ manifest["experiments"].append({"id": j.stem, "output": f"experiments/{wav.name}", "meta": meta})
52+ from challenge_tools.build import player_version
53+
54+ html = (ROOT / "web" / "index.html").read_text().replace('src="player.js"', f'src="player.js?v={player_version()}"')
55+ (out / "index.html").write_text(html)
56+ shutil.copyfile(ROOT / "web" / "viewer.js", out / "viewer.js")
57+ # the shared Web Audio player, from the installed challenge-tools package
58+ from challenge_tools.build import player_js
59+
60+ (out / "player.js").write_text(player_js())
61+ # suggested playback gains and the challenge parameters (strike point, radius), when the challenge repository is checked out beside us (as in the workflow)
62+ gains = {}
63+ params = {}
64+ for cj in sorted((ROOT / "challenge-repo" / "challenges").glob("*/challenge.json")) if (ROOT / "challenge-repo").exists() else []:
65+ ch = json.loads(cj.read_text())
66+ gains[ch["challenge_id"]] = (ch.get("playback") or {}).get("gain_db")
67+ params[ch["challenge_id"]] = {**ch.get("parameters", {}), "strikes": ch.get("strikes", [])}
68+ manifest["playback_gain_db"] = gains
69+ manifest["parameters"] = params
70+ (out / "manifest.json").write_text(json.dumps(manifest, indent=1))
71+ print(f"site written to {out}")
72+
73+
74+if __name__ == "__main__":
75+ main()
web/index.htmladded+56−0View file
@@ -0,0 +1,56 @@
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>cymbal2 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 cymbal2 challenges</h1>
33+<p>This is the site of the solution repository <a href="https://vault1.magland.org/challenges/cymbal2-solutions">challenges/cymbal2-solutions</a>. The comparison of all registered solutions against the reference is on the <a href="https://vault1.magland.org/challenges/cymbal2/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><input id="derigid" type="checkbox" checked> remove the rigid drift (best-fit plane)</label>
42+ <label>Time <input id="timeSlider" type="range" min="0" max="1" step="1" value="0"> <span id="timeLabel" class="muted"></span></label>
43+</div>
44+<div id="view"></div>
45+<p id="fieldStatus" class="muted">Loading…</p>
46+
47+<h2>Submissions</h2>
48+<div id="playback"></div>
49+<div id="submissions"></div>
50+
51+<div id="experiments" hidden></div>
52+</main>
53+<script src="player.js"></script>
54+<script type="module" src="viewer.js"></script>
55+</body>
56+</html>
web/viewer.jsadded+285−0View file
@@ -0,0 +1,285 @@
1+// Lists this repository's submissions and animates the exported HDF5 displacement fields of the disk 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, challenge: cid });
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+// ---------------------------------------------------------------- experiments
36+const expDiv = document.getElementById("experiments");
37+for (const e of manifest.experiments || []) {
38+ const m = e.meta || {};
39+ const p = document.createElement("p");
40+ p.innerHTML = `<b>${esc(e.id)}</b>. ${esc(m.description || m.note || "")}` +
41+ (m.strikes ? ` <span class="muted">Strikes: ${esc(m.strikes.map((s) => `${s.force} N over ${(s.duration * 1e3).toFixed(1)} ms at ${s.t} s, x = ${(s.x * 1e3).toFixed(1)} mm, y = ${(s.y * 1e3).toFixed(1)} mm`).join("; "))}.</span>` : "") +
42+ (m.parameters ? ` <span class="muted"><code>${esc(JSON.stringify(m.parameters))}</code></span>` : "") +
43+ `<br><div class="cplayer" data-src="${esc(e.output)}"></div> <a href="${esc(e.output)}">wav</a>`;
44+ expDiv.appendChild(p);
45+}
46+if (!(manifest.experiments || []).length) expDiv.innerHTML = '<p class="muted">None yet.</p>';
47+
48+const gains = manifest.playback_gain_db || {};
49+const firstGain = Object.values(gains).find((g) => typeof g === "number");
50+ChallengePlayer.controls(document.getElementById("playback"), { gainDb: typeof firstGain === "number" ? firstGain : 20, limiter: true });
51+ChallengePlayer.mount();
52+
53+// ---------------------------------------------------------------- 3D viewer
54+const view = document.getElementById("view");
55+const status = document.getElementById("fieldStatus");
56+const select = document.getElementById("fieldSelect");
57+const playBtn = document.getElementById("playBtn");
58+const speed = document.getElementById("speed");
59+const speedLabel = document.getElementById("speedLabel");
60+const exag = document.getElementById("exag");
61+const exagLabel = document.getElementById("exagLabel");
62+const timeSlider = document.getElementById("timeSlider");
63+const derigid = document.getElementById("derigid");
64+const timeLabel = document.getElementById("timeLabel");
65+
66+const renderer = new THREE.WebGLRenderer({ antialias: true });
67+renderer.setPixelRatio(window.devicePixelRatio || 1);
68+view.appendChild(renderer.domElement);
69+const scene = new THREE.Scene();
70+const dark = matchMedia("(prefers-color-scheme: dark)").matches;
71+scene.background = new THREE.Color(dark ? 0x151515 : 0xffffff);
72+const camera = new THREE.PerspectiveCamera(40, 2, 0.01, 100);
73+camera.position.set(0.45, -0.6, 0.4);
74+camera.up.set(0, 0, 1);
75+const controls = new OrbitControls(camera, renderer.domElement);
76+controls.target.set(0, 0, 0);
77+scene.add(new THREE.AmbientLight(0xffffff, 0.6));
78+const dl = new THREE.DirectionalLight(0xffffff, 0.9);
79+dl.position.set(0.5, -0.8, 1.0);
80+scene.add(dl);
81+
82+let field = null; // {x, y, t, u, nx, ny, nt, umax, emax, plane, vidx, verts, params}
83+let mesh = null;
84+let frame = 0;
85+let playing = false;
86+let lastTs = null;
87+let simTime = 0;
88+
89+function resize() {
90+ const w = view.clientWidth, h = view.clientHeight;
91+ renderer.setSize(w, h, false);
92+ camera.aspect = w / h;
93+ camera.updateProjectionMatrix();
94+}
95+window.addEventListener("resize", () => { resize(); render(); });
96+resize();
97+
98+function buildMesh() {
99+ if (mesh) { scene.remove(mesh); mesh.geometry.dispose(); mesh.material.dispose(); }
100+ const { x, y, nx, ny, vidx, verts, params } = field;
101+ // vertices: the grid points inside the disk; triangles: cells whose four corners are all inside
102+ const pos = new Float32Array(verts.length * 3);
103+ for (let k = 0; k < verts.length; k++) {
104+ const [i, j] = verts[k];
105+ pos[3 * k] = x[i]; pos[3 * k + 1] = y[j]; pos[3 * k + 2] = 0;
106+ }
107+ const idx = [];
108+ for (let j = 0; j < ny - 1; j++) for (let i = 0; i < nx - 1; i++) {
109+ const a = vidx[j * nx + i], b = vidx[j * nx + i + 1], c = vidx[(j + 1) * nx + i], d = vidx[(j + 1) * nx + i + 1];
110+ if (a < 0 || b < 0 || c < 0 || d < 0) continue;
111+ idx.push(a, b, d, a, d, c);
112+ }
113+ const geo = new THREE.BufferGeometry();
114+ geo.setAttribute("position", new THREE.BufferAttribute(pos, 3));
115+ geo.setAttribute("color", new THREE.BufferAttribute(new Float32Array(verts.length * 3), 3));
116+ geo.setIndex(idx);
117+ const mat = new THREE.MeshStandardMaterial({ vertexColors: true, side: THREE.DoubleSide, metalness: 0.1, roughness: 0.7 });
118+ mesh = new THREE.Mesh(geo, mat);
119+ scene.add(mesh);
120+ // the rim at rest, for reference
121+ const R = params.disk_radius || Math.max(...x);
122+ const circ = [];
123+ for (let k = 0; k <= 128; k++) circ.push(new THREE.Vector3(R * Math.cos(2 * Math.PI * k / 128), R * Math.sin(2 * Math.PI * k / 128), 0));
124+ mesh.add(new THREE.Line(new THREE.BufferGeometry().setFromPoints(circ), new THREE.LineBasicMaterial({ color: dark ? 0x555555 : 0xbbbbbb })));
125+ // strike markers
126+ for (const st of (params.strikes || [])) {
127+ const marker = new THREE.Mesh(new THREE.SphereGeometry(0.005, 12, 12), new THREE.MeshBasicMaterial({ color: 0xd04040 }));
128+ marker.position.set(st.x, st.y, 0.004);
129+ mesh.add(marker);
130+ }
131+}
132+
133+function updateSurface() {
134+ if (!field || !mesh) return;
135+ const { u, nx, ny, x, y, plane, verts } = field;
136+ const rigidOff = derigid.checked;
137+ const umax = rigidOff ? field.emax : field.umax;
138+ const pos = mesh.geometry.attributes.position;
139+ const col = mesh.geometry.attributes.color;
140+ const scale = umax > 0 ? (0.03 * Math.pow(10, parseFloat(exag.value)) / umax) * 0.1 : 0; // exaggerated height
141+ const base = frame * nx * ny;
142+ const c0 = plane[3 * frame], gx = plane[3 * frame + 1], gy = plane[3 * frame + 2];
143+ for (let k = 0; k < verts.length; k++) {
144+ const [i, j] = verts[k];
145+ let v = u[base + j * nx + i];
146+ if (rigidOff) v -= c0 + gx * x[i] + gy * y[j];
147+ pos.setZ(k, v * scale);
148+ const s = umax > 0 ? Math.max(-1, Math.min(1, v / umax)) : 0;
149+ col.setXYZ(k, s > 0 ? 1 : 1 + s * 0.8, 1 - Math.abs(s) * 0.7, s < 0 ? 1 : 1 - s * 0.8);
150+ }
151+ pos.needsUpdate = true;
152+ col.needsUpdate = true;
153+ mesh.geometry.computeVertexNormals();
154+ timeSlider.value = String(frame);
155+ timeLabel.textContent = `${(field.t[frame] * 1e3).toFixed(0)} ms`;
156+}
157+
158+function render() { renderer.render(scene, camera); }
159+controls.addEventListener("change", render);
160+
161+function solve3(A, b) {
162+ // Gaussian elimination for the 3 x 3 normal equations of the plane fit
163+ const M = A.map((row, i) => [...row, b[i]]);
164+ for (let c = 0; c < 3; c++) {
165+ let p = c;
166+ for (let r = c + 1; r < 3; r++) if (Math.abs(M[r][c]) > Math.abs(M[p][c])) p = r;
167+ [M[c], M[p]] = [M[p], M[c]];
168+ if (Math.abs(M[c][c]) < 1e-300) return [0, 0, 0];
169+ for (let r = 0; r < 3; r++) if (r !== c) {
170+ const f = M[r][c] / M[c][c];
171+ for (let k = c; k < 4; k++) M[r][k] -= f * M[c][k];
172+ }
173+ }
174+ return [M[0][3] / M[0][0], M[1][3] / M[1][1], M[2][3] / M[2][2]];
175+}
176+
177+async function loadField(entry) {
178+ const path = entry.path;
179+ status.textContent = `Loading ${path} (this may take a while, the file is large)…`;
180+ playing = false;
181+ playBtn.textContent = "Play";
182+ try {
183+ const { FS } = await h5wasm.ready;
184+ const resp = await fetch(path);
185+ const buf = new Uint8Array(await resp.arrayBuffer());
186+ if (buf.length < 200 && new TextDecoder().decode(buf).startsWith("version https://git-lfs")) {
187+ status.textContent = `${path} is a git LFS pointer; the object was not published with this site.`;
188+ return;
189+ }
190+ const name = "current.h5";
191+ try { FS.unlink(name); } catch (_) {}
192+ FS.writeFile(name, buf);
193+ const file = new h5wasm.File(name, "r");
194+ const x = Float64Array.from(file.get("x").value);
195+ const y = Float64Array.from(file.get("y").value);
196+ const t = Float64Array.from(file.get("t").value);
197+ const u = Float32Array.from(file.get("u").value);
198+ file.close();
199+ const nx = x.length, ny = y.length, nt = t.length;
200+ // grid points inside the disk are those with finite values in the first frame
201+ const vidx = new Int32Array(nx * ny).fill(-1);
202+ const verts = [];
203+ for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++) {
204+ if (Number.isFinite(u[j * nx + i])) { vidx[j * nx + i] = verts.length; verts.push([i, j]); }
205+ }
206+ // best-fit plane per frame over the inside points (rigid translation and rotations of the free disk)
207+ let s1 = 0, sx = 0, sy = 0, sxx = 0, sxy = 0, syy = 0;
208+ for (const [i, j] of verts) { s1 += 1; sx += x[i]; sy += y[j]; sxx += x[i] * x[i]; sxy += x[i] * y[j]; syy += y[j] * y[j]; }
209+ const A = [[s1, sx, sy], [sx, sxx, sxy], [sy, sxy, syy]];
210+ const plane = new Float64Array(nt * 3);
211+ let umax = 0, emax = 0;
212+ for (let k = 0; k < nt; k++) {
213+ const base = k * nx * ny;
214+ let b0 = 0, bx = 0, by = 0;
215+ for (const [i, j] of verts) { const v = u[base + j * nx + i]; b0 += v; bx += v * x[i]; by += v * y[j]; }
216+ const [c0, gx, gy] = solve3(A, [b0, bx, by]);
217+ plane[3 * k] = c0; plane[3 * k + 1] = gx; plane[3 * k + 2] = gy;
218+ for (const [i, j] of verts) {
219+ const v = u[base + j * nx + i];
220+ const e = v - c0 - gx * x[i] - gy * y[j];
221+ if (Math.abs(v) > umax) umax = Math.abs(v);
222+ if (Math.abs(e) > emax) emax = Math.abs(e);
223+ }
224+ }
225+ const params = (manifest.parameters || {})[entry.challenge] || {};
226+ field = { x, y, t, u, nx, ny, nt, umax, emax, plane, vidx, verts, params };
227+ timeSlider.max = String(field.nt - 1);
228+ frame = 0;
229+ simTime = 0;
230+ buildMesh();
231+ updateSurface();
232+ render();
233+ status.textContent = `${path}: ${verts.length} grid points inside the disk × ${field.nt} frames, max |u| = ${(umax * 1e3).toFixed(2)} mm including rigid motion, ${(emax * 1e3).toFixed(3)} mm elastic`;
234+ } catch (e) {
235+ status.textContent = `Could not load ${path}: ${e}`;
236+ }
237+}
238+
239+function tick(ts) {
240+ if (!playing || !field) return;
241+ if (lastTs !== null) {
242+ const factor = Math.pow(10, parseFloat(speed.value));
243+ simTime += ((ts - lastTs) / 1000) * factor;
244+ const tEnd = field.t[field.nt - 1];
245+ if (simTime > tEnd) simTime = 0;
246+ const dt = field.t[1] - field.t[0];
247+ frame = Math.min(field.nt - 1, Math.round(simTime / dt));
248+ updateSurface();
249+ render();
250+ }
251+ lastTs = ts;
252+ requestAnimationFrame(tick);
253+}
254+
255+playBtn.addEventListener("click", () => {
256+ playing = !playing;
257+ playBtn.textContent = playing ? "Pause" : "Play";
258+ lastTs = null;
259+ if (playing) requestAnimationFrame(tick);
260+});
261+timeSlider.addEventListener("input", () => {
262+ if (!field) return;
263+ frame = parseInt(timeSlider.value, 10);
264+ simTime = field.t[frame];
265+ updateSurface();
266+ render();
267+});
268+const updSpeed = () => { speedLabel.textContent = `${Math.pow(10, parseFloat(speed.value)).toPrecision(2)} × real time`; };
269+const updExag = () => { exagLabel.textContent = `×${Math.pow(10, parseFloat(exag.value)).toPrecision(2)}`; updateSurface(); render(); };
270+speed.addEventListener("input", updSpeed);
271+exag.addEventListener("input", updExag);
272+derigid.addEventListener("change", () => { updateSurface(); render(); });
273+updSpeed();
274+exagLabel.textContent = `×${Math.pow(10, parseFloat(exag.value)).toPrecision(2)}`;
275+select.addEventListener("change", () => loadField(fieldSolutions[parseInt(select.value, 10)]));
276+
277+fieldSolutions.forEach((f, k) => {
278+ const o = document.createElement("option");
279+ o.value = String(k);
280+ o.textContent = f.label;
281+ select.appendChild(o);
282+});
283+render();
284+if (fieldSolutions.length) loadField(fieldSolutions[0]);
285+else status.textContent = "No solution in this repository exports a displacement field that is published with this site.";