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 K = np.zeros((ndof, ndof))
123 M = np.zeros((ndof, ndof))
124 for ex in range(n):
125 for ey in range(n):
126 nodes = [node(ex, ey), node(ex + 1, ey), node(ex, ey + 1), node(ex + 1, ey + 1)]
127 dofs = np.array([4 * nd + d for nd in nodes for d in range(4)])
128 K[np.ix_(dofs, dofs)] += Ke
129 M[np.ix_(dofs, dofs)] += Me
130 self.node = node
132 evals, evecs = eigh(K, M) # mass-normalized eigenvectors
133 omega2 = np.maximum(evals, 0.0)
134 omega = np.sqrt(omega2)
135 omega[omega < RIGID_OMEGA] = 0.0
136 keep = omega <= 2 * np.pi * MODE_CUTOFF_HZ
137 self.omega = omega[keep]
138 self.V = evecs[:, keep]
139 self.n_modes = int(keep.sum())
140 self.alpha = p["sigma0"] + 0.5 * p["eta"] * self.omega**2
141 self.rigid = self.omega == 0.0
142 self.n_rigid = int(self.rigid.sum())
144 # modal forces: f_n = -int p_hat Phi_n dA over the strike disk
145 R, xs, ys = p["strike_radius"], p["strike_x"], p["strike_y"]
146 c_phi = 0.5 - 2.0 / np.pi**2
147 gr, gw = leg.leggauss(24)
148 r = R * (gr + 1.0) / 2.0
149 wr = gw * R / 2.0
150 nth = 64
151 th = 2 * np.pi * np.arange(nth) / nth
152 rr, tt = np.meshgrid(r, th, indexing="ij")
153 ww = (wr[:, None] * rr * (2 * np.pi / nth)).ravel()
154 px = (xs + rr * np.cos(tt)).ravel()
155 py = (ys + rr * np.sin(tt)).ravel()
156 phat = p["strike_peak_force"] / (c_phi * np.pi * R**2) * 0.5 * (1.0 + np.cos(np.pi * rr.ravel() / R))
157 self.f = -(ww * phat) @ self.modes_at(px, py)
159 # closed-form response constants
160 T = p["strike_duration"]
161 Om = 2 * np.pi / T
162 self.T, self.Om = T, Om
163 al, w, f = self.alpha, self.omega, self.f
164 with np.errstate(divide="ignore", invalid="ignore"):
165 disc = np.sqrt(al**2 - w**2 + 0j)
166 self.s1 = -al + disc
167 self.s2 = -al - disc
168 H = 1.0 / (w**2 - Om**2 + 2j * al * Om)
169 self.qconst = np.where(self.rigid, 0.0, f / (2.0 * w**2))
170 self.Hc = np.where(self.rigid, 0.0, (f / 2.0) * H)
171 rhs1 = -(self.qconst - self.Hc.real)
172 rhs2 = (1j * Om * self.Hc).real
173 self.B = (rhs2 - self.s1 * rhs1) / (self.s2 - self.s1)
174 self.A = rhs1 - self.B
175 self.Gc = (f / 2.0) / (2.0 * al + 1j * Om)
176 self.vconst = f / (4.0 * al)
177 self.Ar = -(self.vconst - self.Gc.real)
178 qT, vT, _ = self._pulse_phase(np.array([T]))
179 self.qT, self.vT = qT[:, 0], vT[:, 0]
180 with np.errstate(divide="ignore", invalid="ignore"):
181 self.B2 = (self.vT - self.s1 * self.qT) / (self.s2 - self.s1)
182 self.A2 = self.qT - self.B2
184 # ------------------------------------------------------------ evaluation at points
185 def basis_at(self, x: np.ndarray, y: np.ndarray) -> np.ndarray:
186 """Global shape-function values at points, shape (npts, ndof), dense."""
187 n, a, b = self.n, self.a, self.b
188 ex = np.clip(np.floor(x / a).astype(int), 0, n - 1)
189 ey = np.clip(np.floor(y / b).astype(int), 0, n - 1)
190 xi = x / a - ex
191 eta = y / b - ey
192 N, _, _, _ = shape_functions(xi, eta, a, b)
193 out = np.zeros((len(x), self.ndof))
194 rows = np.arange(len(x))
195 for c, (cx, cy) in enumerate([(0, 0), (1, 0), (0, 1), (1, 1)]):
196 nd = self.node(ex + cx, ey + cy)
197 for d in range(4):
198 out[rows, 4 * nd + d] += N[:, 4 * c + d]
199 return out
201 def modes_at(self, x: np.ndarray, y: np.ndarray) -> np.ndarray:
202 return self.basis_at(x, y) @ self.V
204 # ------------------------------------------------------------ modal time response
205 def _pulse_phase(self, t: np.ndarray):
206 t = t[None, :]
207 Om = self.Om
208 e = np.exp(1j * Om * t)
209 s1, s2, A, B = self.s1[:, None], self.s2[:, None], self.A[:, None], self.B[:, None]
210 E1, E2 = np.exp(s1 * t), np.exp(s2 * t)
211 Hc = self.Hc[:, None]
212 q = (A * E1 + B * E2).real + self.qconst[:, None] - (Hc * e).real
213 v = (A * s1 * E1 + B * s2 * E2).real - (1j * Om * Hc * e).real
214 acc = (A * s1**2 * E1 + B * s2**2 * E2).real + (Om**2 * Hc * e).real
215 rg = self.rigid
216 if rg.any():
217 al = self.alpha[rg][:, None]
218 Ar, Gc, vc = self.Ar[rg][:, None], self.Gc[rg][:, None], self.vconst[rg][:, None]
219 Ed = np.exp(-2.0 * al * t)
220 q[rg] = Ar * (1.0 - Ed) / (2.0 * al) + vc * t - ((e - 1.0) * Gc / (1j * Om)).real
221 v[rg] = Ar * Ed + vc - (Gc * e).real
222 acc[rg] = -2.0 * al * Ar * Ed - (1j * Om * Gc * e).real
223 return q, v, acc
225 def _free_phase(self, tau: np.ndarray):
226 tau = tau[None, :]
227 s1, s2, A, B = self.s1[:, None], self.s2[:, None], self.A2[:, None], self.B2[:, None]
228 E1, E2 = np.exp(s1 * tau), np.exp(s2 * tau)
229 q = (A * E1 + B * E2).real
230 v = (A * s1 * E1 + B * s2 * E2).real
231 acc = (A * s1**2 * E1 + B * s2**2 * E2).real
232 rg = self.rigid
233 if rg.any():
234 al = self.alpha[rg][:, None]
235 qT, vT = self.qT[rg][:, None], self.vT[rg][:, None]
236 Ed = np.exp(-2.0 * al * tau)
237 q[rg] = qT + vT * (1.0 - Ed) / (2.0 * al)
238 v[rg] = vT * Ed
239 acc[rg] = -2.0 * al * vT * Ed
240 return q, v, acc
242 def response(self, t: np.ndarray):
243 q = np.zeros((self.n_modes, len(t)))
244 v = np.zeros_like(q)
245 acc = np.zeros_like(q)
246 m1 = (t >= 0) & (t <= self.T)
247 m2 = t > self.T
248 if m1.any():
249 q[:, m1], v[:, m1], acc[:, m1] = self._pulse_phase(t[m1])
250 if m2.any():
251 q[:, m2], v[:, m2], acc[:, m2] = self._free_phase(t[m2] - self.T)
252 return q, v, acc
254 # ------------------------------------------------------------ Rayleigh integral
255 def pressure(self, fs: int, nsamp: int) -> np.ndarray:
256 p = self.p
257 n, a, b = self.n, self.a, self.b
258 g, wg = leg.leggauss(3)
259 s = (g + 1) / 2
260 ws = wg / 2
261 ex, ey = np.meshgrid(np.arange(n), np.arange(n), indexing="ij")
262 XI, ETA = np.meshgrid(s, s, indexing="ij")
263 xq = (ex.ravel()[:, None] + XI.ravel()[None, :]).ravel() * a
264 yq = (ey.ravel()[:, None] + ETA.ravel()[None, :]).ravel() * b
265 WW = np.tile((ws[:, None] * ws[None, :]).ravel() * a * b, n * n)
266 r = np.sqrt((p["listener_x"] - xq) ** 2 + (p["listener_y"] - yq) ** 2 + p["listener_height"] ** 2)
267 tau = r / p["c_air"]
268 cq = p["rho_air"] / (2 * np.pi) * WW / r
269 Phi = self.modes_at(xq, yq) # (nq, n_modes)
270 t = np.arange(nsamp) / fs
271 out = np.zeros(nsamp)
272 t_split = self.T + tau.max()
273 early = (t < t_split) & (t >= tau.min())
274 late = t >= t_split
275 ne = int(early.sum())
276 if ne:
277 chunk = max(1, 4_000_000 // max(1, self.n_modes * ne))
278 for i0 in range(0, len(xq), chunk):
279 sl = slice(i0, i0 + chunk)
280 tt = (t[early][None, :] - tau[sl][:, None]).ravel()
281 _, _, acc = self.response(tt)
282 acc = acc.reshape(self.n_modes, -1, ne)
283 out[early] += np.einsum("q,qn,nqk->k", cq[sl], Phi[sl], acc, optimize=True)
284 if late.any():
285 tl = t[late] - self.T
286 rg = self.rigid
287 G1 = (cq[:, None] * Phi * np.exp(-np.outer(tau, self.s1))).sum(axis=0)
288 G2 = (cq[:, None] * Phi * np.exp(-np.outer(tau, self.s2))).sum(axis=0)
289 c1 = self.A2 * self.s1**2 * G1
290 c2 = self.B2 * self.s2**2 * G2
291 s1, s2 = self.s1.copy(), self.s2.copy()
292 if rg.any():
293 al = self.alpha[rg]
294 Gr = (cq[:, None] * Phi[:, rg] * np.exp(np.outer(tau, 2.0 * al))).sum(axis=0)
295 c1[rg] = -2.0 * al * self.vT[rg] * Gr
296 c2[rg] = 0.0
297 s1[rg] = -2.0 * al
298 s2[rg] = 0.0
299 coef = np.concatenate([c1, c2])
300 rates = np.concatenate([s1, s2])
301 res = np.zeros(len(tl))
302 step = max(1, 20_000_000 // len(coef))
303 for k0 in range(0, len(tl), step):
304 tk = tl[k0:k0 + step]
305 res[k0:k0 + step] = (np.exp(np.outer(tk, rates)) @ coef).real
306 out[late] = res
307 return out
309 def field(self, x: np.ndarray, y: np.ndarray, t: np.ndarray) -> np.ndarray:
310 XX, YY = np.meshgrid(x, y, indexing="xy") # (ny, nx)
311 Phi = self.modes_at(XX.ravel(), YY.ravel())
312 q, _, _ = self.response(t)
313 return (Phi @ q).T.reshape(len(t), len(y), len(x))
316def main():
317 ch = load_challenge(Path(sys.argv[1]))
318 out_dir = Path(sys.argv[2])
319 args = sys.argv[3:]
320 fields_for = set()
321 if "--fields" in args:
322 i = args.index("--fields")
323 fields_for = {int(v) for v in args[i + 1].split(",")}
324 args = args[:i] + args[i + 2:]
325 meshes = [int(a) for a in args] or [4, 6, 8, 12, 16, 24, 32]
326 out = ch["output"]
327 fs, nsamp = out["sample_rate_hz"], out["num_samples"]
328 meta = []
329 for n in meshes:
330 t0 = time.perf_counter()
331 plate = BFSPlate(ch, n)
332 t1 = time.perf_counter()
333 s = plate.pressure(fs, nsamp)
334 dt = time.perf_counter() - t0
335 sid = f"bfs-n{n}"
336 out_dir.mkdir(parents=True, exist_ok=True)
337 wavfile.write(out_dir / f"{sid}.wav", fs, (s / out["signal_scale"]).astype(np.float32))
338 low = plate.omega[plate.omega > 0][:6] / (2 * np.pi)
339 entry = {
340 "solution_id": sid,
341 "output": f"{out_dir.name}/{sid}.wav",
342 "method": "Kirchhoff-Love plate, Bogner-Fox-Schmit C1 finite elements, closed-form modal time response, exact retarded Rayleigh integral",
343 "parameters": {"elements_per_side": n, "element_size_m": round(plate.a, 6), "n_dof": plate.ndof, "n_modes": plate.n_modes, "rigid_modes": plate.n_rigid, "mode_cutoff_hz": MODE_CUTOFF_HZ},
344 "timing_seconds": round(dt, 3),
345 "peak_si": float(np.max(np.abs(s))),
346 "lowest_frequencies_hz": [round(float(v), 3) for v in low],
347 }
348 if n in fields_for:
349 import h5py
351 fg = ch["fields"]
352 nt, ny, nx = fg["shape"]
353 x = np.linspace(0, ch["parameters"]["Lx"], nx)
354 y = np.linspace(0, ch["parameters"]["Ly"], ny)
355 tt = np.arange(nt) / fg["frame_rate_hz"]
356 u = plate.field(x, y, tt)
357 with h5py.File(out_dir / f"{sid}.h5", "w") as f:
358 f.attrs["challenge_id"] = ch["challenge_id"]
359 f.attrs["solution_id"] = sid
360 f.attrs["quantity"] = fg["quantity"]
361 f.attrs["units"] = fg["units"]
362 f.create_dataset("x", data=x)
363 f.create_dataset("y", data=y)
364 f.create_dataset("t", data=tt)
365 f.create_dataset("u", data=u.astype(np.float32), compression="gzip", compression_opts=4, chunks=(8, ny, nx))
366 entry["fields"] = {"path": f"{out_dir.name}/{sid}.h5", "format": "hdf5"}
367 meta.append(entry)
368 print(f"n={n:3d} dof={plate.ndof:5d} modes={plate.n_modes:5d} rigid={plate.n_rigid} setup={t1 - t0:6.1f}s total={dt:6.1f}s peak={np.max(np.abs(s)):.4e} Pa f_low={np.round(low, 2)}", file=sys.stderr)
369 print(json.dumps(meta, indent=1))
372if __name__ == "__main__":
373 main()