1"""Kirchhoff-Love thin-plate solution of plate-strike-v1 with Bogner-Fox-Schmit finite elements.
3Space: a uniform n x n mesh of C1-conforming bicubic Hermite rectangles (Bogner-Fox-Schmit),
4four degrees of freedom per node (w, w_x, w_y, w_xy). Free edges are natural boundary
5conditions of the Kirchhoff energy, so nothing is imposed. Element matrices are integrated
6exactly with 4 x 4 Gauss-Legendre points. The generalized eigenproblem gives mass-normalized
7modes, three of them rigid-body modes at zero frequency.
9Time: each modal equation
11 q'' + 2 alpha_n q' + omega_n^2 q = f_n psi(t / T_s), alpha_n = sigma0 + eta omega_n^2 / 2,
13is solved in closed form: a particular solution for the Hann pulse (a constant plus a cosine)
14during the strike, and free decay afterwards.
16Radiation: the Rayleigh integral is evaluated with 3 x 3 Gauss-Legendre points per element.
17The retardation is exact: once the strike has ended at every retarded time the pressure is a sum
18of complex exponentials whose amplitudes absorb the delays; the few earlier samples are evaluated
19directly at the retarded times.
21The transverse displacement w is positive in +z, so the strike enters as -p and w is reported
22directly as u_z. Modes above 200 kHz are discarded; the pulse excites them at a negligible level
23and their extreme damping rates would otherwise overflow the delay factors.
25Usage: python kl_bfs_plate.py <challenge.json> <out-dir> [--fields n,...] [n ...]
26"""
27from __future__ import annotations
29import json
30import sys
31import time
32from pathlib import Path
34import numpy as np
35from numpy.polynomial import legendre as leg
36from scipy.io import wavfile
37from scipy.linalg import eigh
39MODE_CUTOFF_HZ = 200_000.0
40RIGID_OMEGA = 1.0 # rad/s; eigenvalues below this are the three rigid-body modes
43def load_challenge(path: Path) -> dict:
44 with open(path) as f:
45 return json.load(f)
48# ---------------------------------------------------------------- Hermite cubics on [0, 1]
49def 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
58def 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.
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
86class 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
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])
114 # assembly into dense global matrices
115 nn = n + 1
116 ndof = 4 * nn * nn
117 self.ndof = ndof
119 def node(i, j):
120 return i * nn + j
122 self.node = node
123 self._elem_quad = (XI.ravel(), ETA.ravel(), WQ)
124 self._shape = (N, Nxx, Nyy, Nxy)
125 K = self.assemble(Ke)
126 M = self.assemble(Me)
127 self.K, self.M = K, M
129 evals, evecs = eigh(K, M) # mass-normalized eigenvectors
130 omega2 = np.maximum(evals, 0.0)
131 omega = np.sqrt(omega2)
132 omega[omega < RIGID_OMEGA] = 0.0
133 keep = omega <= 2 * np.pi * MODE_CUTOFF_HZ
134 self.omega = omega[keep]
135 self.V = evecs[:, keep]
136 self.n_modes = int(keep.sum())
137 self.alpha = p["sigma0"] + 0.5 * p["eta"] * self.omega**2
138 self.rigid = self.omega == 0.0
139 self.n_rigid = int(self.rigid.sum())
141 # modal forces: f_n = -int p_hat Phi_n dA over the strike disk
142 R, xs, ys = p["strike_radius"], p["strike_x"], p["strike_y"]
143 c_phi = 0.5 - 2.0 / np.pi**2
144 gr, gw = leg.leggauss(24)
145 r = R * (gr + 1.0) / 2.0
146 wr = gw * R / 2.0
147 nth = 64
148 th = 2 * np.pi * np.arange(nth) / nth
149 rr, tt = np.meshgrid(r, th, indexing="ij")
150 ww = (wr[:, None] * rr * (2 * np.pi / nth)).ravel()
151 px = (xs + rr * np.cos(tt)).ravel()
152 py = (ys + rr * np.sin(tt)).ravel()
153 phat = p["strike_peak_force"] / (c_phi * np.pi * R**2) * 0.5 * (1.0 + np.cos(np.pi * rr.ravel() / R))
154 self.f = -(ww * phat) @ self.modes_at(px, py)
156 # closed-form response constants
157 T = p["strike_duration"]
158 Om = 2 * np.pi / T
159 self.T, self.Om = T, Om
160 al, w, f = self.alpha, self.omega, self.f
161 with np.errstate(divide="ignore", invalid="ignore"):
162 disc = np.sqrt(al**2 - w**2 + 0j)
163 self.s1 = -al + disc
164 self.s2 = -al - disc
165 H = 1.0 / (w**2 - Om**2 + 2j * al * Om)
166 self.qconst = np.where(self.rigid, 0.0, f / (2.0 * w**2))
167 self.Hc = np.where(self.rigid, 0.0, (f / 2.0) * H)
168 rhs1 = -(self.qconst - self.Hc.real)
169 rhs2 = (1j * Om * self.Hc).real
170 self.B = (rhs2 - self.s1 * rhs1) / (self.s2 - self.s1)
171 self.A = rhs1 - self.B
172 self.Gc = (f / 2.0) / (2.0 * al + 1j * Om)
173 self.vconst = f / (4.0 * al)
174 self.Ar = -(self.vconst - self.Gc.real)
175 qT, vT, _ = self._pulse_phase(np.array([T]))
176 self.qT, self.vT = qT[:, 0], vT[:, 0]
177 with np.errstate(divide="ignore", invalid="ignore"):
178 self.B2 = (self.vT - self.s1 * self.qT) / (self.s2 - self.s1)
179 self.A2 = self.qT - self.B2
181 # ------------------------------------------------------------ assembly helpers
182 def element_dofs(self, ex: int, ey: int) -> np.ndarray:
183 nodes = [self.node(ex, ey), self.node(ex + 1, ey), self.node(ex, ey + 1), self.node(ex + 1, ey + 1)]
184 return np.array([4 * nd + d for nd in nodes for d in range(4)])
186 def assemble(self, Ae: np.ndarray) -> np.ndarray:
187 """Scatter one 16 x 16 element matrix (uniform mesh) into a dense global matrix."""
188 A = np.zeros((self.ndof, self.ndof))
189 for ex in range(self.n):
190 for ey in range(self.n):
191 dofs = self.element_dofs(ex, ey)
192 A[np.ix_(dofs, dofs)] += Ae
193 return A
195 def element_matrix(self, Cmat: np.ndarray) -> np.ndarray:
196 """Element 'stiffness' for the bilinear form int [w_xx, w_yy, w_xy] Cmat [v_xx, v_yy, v_xy]^T dA (no D factor)."""
197 _, Nxx, Nyy, Nxy = self._shape
198 WQ = self._elem_quad[2]
199 Ae = np.zeros((16, 16))
200 for q in range(len(WQ)):
201 B2 = np.stack([Nxx[q], Nyy[q], Nxy[q]])
202 Ae += WQ[q] * (B2.T @ Cmat @ B2)
203 return Ae
205 def boundary_dofs(self) -> np.ndarray:
206 n = self.n
207 b = []
208 for i in range(n + 1):
209 for j in range(n + 1):
210 if i in (0, n) or j in (0, n):
211 b.extend(4 * self.node(i, j) + d for d in range(4))
212 return np.array(sorted(set(b)))
214 # ------------------------------------------------------------ evaluation at points
215 def _locate(self, x: np.ndarray, y: np.ndarray):
216 n, a, b = self.n, self.a, self.b
217 ex = np.clip(np.floor(x / a).astype(int), 0, n - 1)
218 ey = np.clip(np.floor(y / b).astype(int), 0, n - 1)
219 return ex, ey, x / a - ex, y / b - ey
221 def _scatter(self, vals: np.ndarray, ex: np.ndarray, ey: np.ndarray) -> np.ndarray:
222 out = np.zeros((len(ex), self.ndof))
223 rows = np.arange(len(ex))
224 for c, (cx, cy) in enumerate([(0, 0), (1, 0), (0, 1), (1, 1)]):
225 nd = self.node(ex + cx, ey + cy)
226 for d in range(4):
227 out[rows, 4 * nd + d] += vals[:, 4 * c + d]
228 return out
230 def basis_derivs_at(self, x: np.ndarray, y: np.ndarray):
231 """Global shape-function second derivatives (xx, yy, xy) at points, each (npts, ndof)."""
232 ex, ey, xi, eta = self._locate(x, y)
233 _, Nxx, Nyy, Nxy = shape_functions(xi, eta, self.a, self.b)
234 return self._scatter(Nxx, ex, ey), self._scatter(Nyy, ex, ey), self._scatter(Nxy, ex, ey)
236 def basis_at(self, x: np.ndarray, y: np.ndarray) -> np.ndarray:
237 """Global shape-function values at points, shape (npts, ndof), dense."""
238 n, a, b = self.n, self.a, self.b
239 ex = np.clip(np.floor(x / a).astype(int), 0, n - 1)
240 ey = np.clip(np.floor(y / b).astype(int), 0, n - 1)
241 xi = x / a - ex
242 eta = y / b - ey
243 N, _, _, _ = shape_functions(xi, eta, a, b)
244 out = np.zeros((len(x), self.ndof))
245 rows = np.arange(len(x))
246 for c, (cx, cy) in enumerate([(0, 0), (1, 0), (0, 1), (1, 1)]):
247 nd = self.node(ex + cx, ey + cy)
248 for d in range(4):
249 out[rows, 4 * nd + d] += N[:, 4 * c + d]
250 return out
252 def modes_at(self, x: np.ndarray, y: np.ndarray) -> np.ndarray:
253 return self.basis_at(x, y) @ self.V
255 # ------------------------------------------------------------ modal time response
256 def _pulse_phase(self, t: np.ndarray):
257 t = t[None, :]
258 Om = self.Om
259 e = np.exp(1j * Om * t)
260 s1, s2, A, B = self.s1[:, None], self.s2[:, None], self.A[:, None], self.B[:, None]
261 E1, E2 = np.exp(s1 * t), np.exp(s2 * t)
262 Hc = self.Hc[:, None]
263 q = (A * E1 + B * E2).real + self.qconst[:, None] - (Hc * e).real
264 v = (A * s1 * E1 + B * s2 * E2).real - (1j * Om * Hc * e).real
265 acc = (A * s1**2 * E1 + B * s2**2 * E2).real + (Om**2 * Hc * e).real
266 rg = self.rigid
267 if rg.any():
268 al = self.alpha[rg][:, None]
269 Ar, Gc, vc = self.Ar[rg][:, None], self.Gc[rg][:, None], self.vconst[rg][:, None]
270 Ed = np.exp(-2.0 * al * t)
271 q[rg] = Ar * (1.0 - Ed) / (2.0 * al) + vc * t - ((e - 1.0) * Gc / (1j * Om)).real
272 v[rg] = Ar * Ed + vc - (Gc * e).real
273 acc[rg] = -2.0 * al * Ar * Ed - (1j * Om * Gc * e).real
274 return q, v, acc
276 def _free_phase(self, tau: np.ndarray):
277 tau = tau[None, :]
278 s1, s2, A, B = self.s1[:, None], self.s2[:, None], self.A2[:, None], self.B2[:, None]
279 E1, E2 = np.exp(s1 * tau), np.exp(s2 * tau)
280 q = (A * E1 + B * E2).real
281 v = (A * s1 * E1 + B * s2 * E2).real
282 acc = (A * s1**2 * E1 + B * s2**2 * E2).real
283 rg = self.rigid
284 if rg.any():
285 al = self.alpha[rg][:, None]
286 qT, vT = self.qT[rg][:, None], self.vT[rg][:, None]
287 Ed = np.exp(-2.0 * al * tau)
288 q[rg] = qT + vT * (1.0 - Ed) / (2.0 * al)
289 v[rg] = vT * Ed
290 acc[rg] = -2.0 * al * vT * Ed
291 return q, v, acc
293 def response(self, t: np.ndarray):
294 q = np.zeros((self.n_modes, len(t)))
295 v = np.zeros_like(q)
296 acc = np.zeros_like(q)
297 m1 = (t >= 0) & (t <= self.T)
298 m2 = t > self.T
299 if m1.any():
300 q[:, m1], v[:, m1], acc[:, m1] = self._pulse_phase(t[m1])
301 if m2.any():
302 q[:, m2], v[:, m2], acc[:, m2] = self._free_phase(t[m2] - self.T)
303 return q, v, acc
305 # ------------------------------------------------------------ Rayleigh integral
306 def pressure(self, fs: int, nsamp: int) -> np.ndarray:
307 p = self.p
308 n, a, b = self.n, self.a, self.b
309 g, wg = leg.leggauss(3)
310 s = (g + 1) / 2
311 ws = wg / 2
312 ex, ey = np.meshgrid(np.arange(n), np.arange(n), indexing="ij")
313 XI, ETA = np.meshgrid(s, s, indexing="ij")
314 xq = (ex.ravel()[:, None] + XI.ravel()[None, :]).ravel() * a
315 yq = (ey.ravel()[:, None] + ETA.ravel()[None, :]).ravel() * b
316 WW = np.tile((ws[:, None] * ws[None, :]).ravel() * a * b, n * n)
317 r = np.sqrt((p["listener_x"] - xq) ** 2 + (p["listener_y"] - yq) ** 2 + p["listener_height"] ** 2)
318 tau = r / p["c_air"]
319 cq = p["rho_air"] / (2 * np.pi) * WW / r
320 Phi = self.modes_at(xq, yq) # (nq, n_modes)
321 t = np.arange(nsamp) / fs
322 out = np.zeros(nsamp)
323 t_split = self.T + tau.max()
324 early = (t < t_split) & (t >= tau.min())
325 late = t >= t_split
326 ne = int(early.sum())
327 if ne:
328 chunk = max(1, 4_000_000 // max(1, self.n_modes * ne))
329 for i0 in range(0, len(xq), chunk):
330 sl = slice(i0, i0 + chunk)
331 tt = (t[early][None, :] - tau[sl][:, None]).ravel()
332 _, _, acc = self.response(tt)
333 acc = acc.reshape(self.n_modes, -1, ne)
334 out[early] += np.einsum("q,qn,nqk->k", cq[sl], Phi[sl], acc, optimize=True)
335 if late.any():
336 tl = t[late] - self.T
337 rg = self.rigid
338 G1 = (cq[:, None] * Phi * np.exp(-np.outer(tau, self.s1))).sum(axis=0)
339 G2 = (cq[:, None] * Phi * np.exp(-np.outer(tau, self.s2))).sum(axis=0)
340 c1 = self.A2 * self.s1**2 * G1
341 c2 = self.B2 * self.s2**2 * G2
342 s1, s2 = self.s1.copy(), self.s2.copy()
343 if rg.any():
344 al = self.alpha[rg]
345 Gr = (cq[:, None] * Phi[:, rg] * np.exp(np.outer(tau, 2.0 * al))).sum(axis=0)
346 c1[rg] = -2.0 * al * self.vT[rg] * Gr
347 c2[rg] = 0.0
348 s1[rg] = -2.0 * al
349 s2[rg] = 0.0
350 coef = np.concatenate([c1, c2])
351 rates = np.concatenate([s1, s2])
352 res = np.zeros(len(tl))
353 step = max(1, 20_000_000 // len(coef))
354 for k0 in range(0, len(tl), step):
355 tk = tl[k0:k0 + step]
356 res[k0:k0 + step] = (np.exp(np.outer(tk, rates)) @ coef).real
357 out[late] = res
358 return out
360 def field(self, x: np.ndarray, y: np.ndarray, t: np.ndarray) -> np.ndarray:
361 XX, YY = np.meshgrid(x, y, indexing="xy") # (ny, nx)
362 Phi = self.modes_at(XX.ravel(), YY.ravel())
363 q, _, _ = self.response(t)
364 return (Phi @ q).T.reshape(len(t), len(y), len(x))
367def main():
368 ch = load_challenge(Path(sys.argv[1]))
369 out_dir = Path(sys.argv[2])
370 args = sys.argv[3:]
371 fields_for = set()
372 if "--fields" in args:
373 i = args.index("--fields")
374 fields_for = {int(v) for v in args[i + 1].split(",")}
375 args = args[:i] + args[i + 2:]
376 meshes = [int(a) for a in args] or [4, 6, 8, 12, 16, 24, 32]
377 out = ch["output"]
378 fs, nsamp = out["sample_rate_hz"], out["num_samples"]
379 meta = []
380 for n in meshes:
381 t0 = time.perf_counter()
382 plate = BFSPlate(ch, n)
383 t1 = time.perf_counter()
384 s = plate.pressure(fs, nsamp)
385 dt = time.perf_counter() - t0
386 sid = f"bfs-n{n}"
387 out_dir.mkdir(parents=True, exist_ok=True)
388 wavfile.write(out_dir / f"{sid}.wav", fs, (s / out["signal_scale"]).astype(np.float32))
389 low = plate.omega[plate.omega > 0][:6] / (2 * np.pi)
390 entry = {
391 "solution_id": sid,
392 "output": f"{out_dir.name}/{sid}.wav",
393 "method": "Kirchhoff-Love plate, Bogner-Fox-Schmit C1 finite elements, closed-form modal time response, exact retarded Rayleigh integral",
394 "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},
395 "timing_seconds": round(dt, 3),
396 "peak_si": float(np.max(np.abs(s))),
397 "lowest_frequencies_hz": [round(float(v), 3) for v in low],
398 }
399 if n in fields_for:
400 import h5py
402 fg = ch["fields"]
403 nt, ny, nx = fg["shape"]
404 x = np.linspace(0, ch["parameters"]["Lx"], nx)
405 y = np.linspace(0, ch["parameters"]["Ly"], ny)
406 tt = np.arange(nt) / fg["frame_rate_hz"]
407 u = plate.field(x, y, tt)
408 with h5py.File(out_dir / f"{sid}.h5", "w") as f:
409 f.attrs["challenge_id"] = ch["challenge_id"]
410 f.attrs["solution_id"] = sid
411 f.attrs["quantity"] = fg["quantity"]
412 f.attrs["units"] = fg["units"]
413 f.create_dataset("x", data=x)
414 f.create_dataset("y", data=y)
415 f.create_dataset("t", data=tt)
416 f.create_dataset("u", data=u.astype(np.float32), compression="gzip", compression_opts=4, chunks=(8, ny, nx))
417 entry["fields"] = {"path": f"{out_dir.name}/{sid}.h5", "format": "hdf5"}
418 meta.append(entry)
419 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)
420 print(json.dumps(meta, indent=1))
423if __name__ == "__main__":
424 main()