1"""Kirchhoff-Love plate on a Cartesian grid masked to a disk, with Bogner-Fox-Schmit elements.
3The disk of radius R is approximated by the union of the square cells of a uniform Cartesian grid
4of spacing a whose centers lie inside the disk (a staircase boundary; the grid is centered on the
5disk so the approximation has the symmetries of the square lattice). Each cell is a C1-conforming
6bicubic Hermite rectangle with four degrees of freedom per node (w, w_x, w_y, w_xy). Free edges,
7including the staircase corners, are natural boundary conditions of the Kirchhoff energy, so
8nothing is imposed. Element matrices are integrated exactly with 4 x 4 Gauss-Legendre points and
9assembled into sparse matrices; the lowest modes come from a shift-invert Lanczos eigensolver
10(ARPACK through scipy) and are mass-normalized. Three of them are rigid-body modes at zero
11frequency.
13The class also provides the closed-form linear modal response to the challenge's Hann strike
14(a particular solution during the pulse and free decay afterwards) and an exact evaluation of the
15retarded Rayleigh integral, both taken from the square-plate solver of challenges/plate-solutions,
16so that the linear limit of the nonlinear solver can be checked on the same mesh.
18Coordinates are those of the challenge: the disk is centered at the origin, w is positive in +z,
19the strike enters as -p, and w is reported directly as u_z.
20"""
21from __future__ import annotations
23import json
24from pathlib import Path
26import numpy as np
27import scipy.sparse as sp
28from numpy.polynomial import legendre as leg
29from scipy.sparse.linalg import eigsh
30from scipy.spatial import cKDTree
32RIGID_OMEGA = 1.0 # rad/s; eigenvalues below this are the three rigid-body modes
35def load_challenge(path: Path) -> dict:
36 with open(path) as f:
37 return json.load(f)
40# ---------------------------------------------------------------- Hermite cubics on [0, 1]
41def hermite(s: np.ndarray):
42 """Values and first two derivatives (w.r.t. s) of the four cubic Hermite functions. Shapes (len(s), 4)."""
43 s = np.asarray(s, dtype=float)
44 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)
45 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)
46 d2H = np.stack([-6 + 12 * s, -4 + 6 * s, 6 - 12 * s, -2 + 6 * s], axis=-1)
47 return H, dH, d2H
50CORNERS = [(0, 0), (1, 0), (0, 1), (1, 1)]
53def shape_functions(xi: np.ndarray, eta: np.ndarray, a: float):
54 """BFS shape functions and second derivatives at local coordinates (xi, eta) for square a x a elements.
56 Returns N, N_xx, N_yy, N_xy, each of shape (npts, 16). Local dof order: for corner c in
57 [(0,0), (1,0), (0,1), (1,1)]: w, w_x, w_y, w_xy.
58 """
59 Hx, dHx, d2Hx = hermite(xi)
60 Hy, dHy, d2Hy = hermite(eta)
61 N = np.zeros((len(xi), 16))
62 Nxx = np.zeros_like(N)
63 Nyy = np.zeros_like(N)
64 Nxy = np.zeros_like(N)
65 for c, (cx, cy) in enumerate(CORNERS):
66 ix0, ix1 = (0, 1) if cx == 0 else (2, 3)
67 iy0, iy1 = (0, 1) if cy == 0 else (2, 3)
68 combos = [(ix0, iy0, 1.0), (ix1, iy0, a), (ix0, iy1, a), (ix1, iy1, a * a)]
69 for d, (fx, fy, sc) in enumerate(combos):
70 k = 4 * c + d
71 N[:, k] = sc * Hx[:, fx] * Hy[:, fy]
72 Nxx[:, k] = sc * d2Hx[:, fx] * Hy[:, fy] / a**2
73 Nyy[:, k] = sc * Hx[:, fx] * d2Hy[:, fy] / a**2
74 Nxy[:, k] = sc * dHx[:, fx] * dHy[:, fy] / (a * a)
75 return N, Nxx, Nyy, Nxy
78class BFSDisk:
79 def __init__(self, ch: dict, spacing: float, n_modes: int, symmetric_y: bool = False, max_hz: float | None = None):
80 """symmetric_y restricts the basis to functions even in y (w(x, -y) = w(x, y)), which is exact for a strike
81 centered on the x axis; the grid is symmetric about that axis by construction. max_hz drops modes above it."""
82 p = ch["parameters"]
83 self.p = p
84 R, h, nu = p["disk_radius"], p["h"], p["nu"]
85 self.R = R
86 a = float(spacing)
87 self.a = a
88 D = p["E"] * h**3 / (12.0 * (1.0 - nu**2))
89 m = p["rho"] * h
90 self.D, self.m = D, m
92 # ---- the masked grid
93 half = int(np.ceil(R / a)) + 1
94 ncell = 2 * half
95 self.ncell = ncell
96 self.half = half
97 self.x0 = -half * a # coordinate of grid line 0
98 centers = self.x0 + (np.arange(ncell) + 0.5) * a
99 CX, CY = np.meshgrid(centers, centers, indexing="ij")
100 active = CX**2 + CY**2 <= R**2
101 cells = np.argwhere(active) # (nc, 2): ix, iy
102 self.cells = cells
103 self.n_cells = len(cells)
104 self.cell_id = -np.ones((ncell, ncell), dtype=int)
105 self.cell_id[cells[:, 0], cells[:, 1]] = np.arange(self.n_cells)
106 self.cell_center = np.stack([centers[cells[:, 0]], centers[cells[:, 1]]], axis=1)
107 self._tree = cKDTree(self.cell_center)
108 nn = ncell + 1
109 gnode = np.zeros((self.n_cells, 4), dtype=int)
110 for c, (cx, cy) in enumerate(CORNERS):
111 gnode[:, c] = (cells[:, 0] + cx) * nn + (cells[:, 1] + cy)
112 touch = np.bincount(gnode.ravel(), minlength=nn * nn)
113 used = touch > 0
114 self.node_map = -np.ones(nn * nn, dtype=int)
115 self.node_map[used] = np.arange(int(used.sum()))
116 self.n_nodes = int(used.sum())
117 self.ndof = 4 * self.n_nodes
118 cnode = self.node_map[gnode] # (nc, 4) compact node ids
119 gidx = np.flatnonzero(used)
120 self.node_i, self.node_j = gidx // nn, gidx % nn # grid indices of each compact node
121 self.cell_dofs = (4 * cnode[:, :, None] + np.arange(4)[None, None, :]).reshape(self.n_cells, 16)
122 bnodes = self.node_map[used & (touch < 4)]
123 self.boundary_dofs = (4 * bnodes[:, None] + np.arange(4)[None, :]).ravel()
124 self.interior_dofs = np.setdiff1d(np.arange(self.ndof), self.boundary_dofs)
125 self.area = self.n_cells * a * a
127 # ---- element matrices
128 g, wg = leg.leggauss(4)
129 s = (g + 1) / 2
130 ws = wg / 2
131 XI, ETA = np.meshgrid(s, s, indexing="ij")
132 WQ = (ws[:, None] * ws[None, :]).ravel() * a * a
133 N, Nxx, Nyy, Nxy = shape_functions(XI.ravel(), ETA.ravel(), a)
134 self._shape = (N, Nxx, Nyy, Nxy, WQ)
135 C = np.array([[1.0, nu, 0.0], [nu, 1.0, 0.0], [0.0, 0.0, 2.0 * (1.0 - nu)]])
136 Ke = D * self.element_matrix(C)
137 Me = np.zeros((16, 16))
138 for q in range(len(WQ)):
139 Me += WQ[q] * m * np.outer(N[q], N[q])
140 self.K = self.assemble(Ke)
141 self.M = self.assemble(Me)
143 # ---- modes: the n_modes lowest, mass-normalized (shift-invert Lanczos around a negative shift)
144 self.symmetric_y = symmetric_y
145 evals, evecs = self.eigenmodes(self.K, self.M, n_modes, sigma=-1.0e3)
146 omega = np.sqrt(np.maximum(evals, 0.0))
147 omega[omega < RIGID_OMEGA] = 0.0
148 if max_hz is not None:
149 keep = omega <= 2 * np.pi * max_hz
150 omega, evecs = omega[keep], evecs[:, keep]
151 self.omega = omega
152 self.V = evecs
153 self.n_modes = len(omega)
154 self.alpha = p["sigma0"] + 0.5 * p["eta"] * omega**2
155 self.rigid = omega == 0.0
156 self.n_rigid = int(self.rigid.sum())
158 # ---- modal forces: f_n = -int p_hat Phi_n dA over the strike disk (polar Gauss quadrature)
159 Rs, xs, ys = p["strike_radius"], p["strike_x"], p["strike_y"]
160 c_phi = 0.5 - 2.0 / np.pi**2
161 gr, gw = leg.leggauss(24)
162 r = Rs * (gr + 1.0) / 2.0
163 wr = gw * Rs / 2.0
164 nth = 64
165 th = 2 * np.pi * np.arange(nth) / nth
166 rr, tt = np.meshgrid(r, th, indexing="ij")
167 ww = (wr[:, None] * rr * (2 * np.pi / nth)).ravel()
168 px = (xs + rr * np.cos(tt)).ravel()
169 py = (ys + rr * np.sin(tt)).ravel()
170 phat = p["strike_peak_force"] / (c_phi * np.pi * Rs**2) * 0.5 * (1.0 + np.cos(np.pi * rr.ravel() / Rs))
171 self.f = -(ww * phat) @ self.modes_at(px, py)
173 self._listener_q = None # optional override of the listener quadrature (used by the shell solver)
174 self._setup_linear_response()
176 def _setup_linear_response(self):
177 """Closed-form constants of the damped modal response to the Hann strike (see the square-plate solver).
178 Depends on omega, alpha, f and rigid; call again after changing them."""
179 p = self.p
180 # ---- closed-form linear response constants (see the square-plate solver)
181 T = p["strike_duration"]
182 Om = 2 * np.pi / T
183 self.T, self.Om = T, Om
184 al, w, f = self.alpha, self.omega, self.f
185 with np.errstate(divide="ignore", invalid="ignore"):
186 disc = np.sqrt(al**2 - w**2 + 0j)
187 self.s1 = -al + disc
188 self.s2 = -al - disc
189 H = 1.0 / (w**2 - Om**2 + 2j * al * Om)
190 self.qconst = np.where(self.rigid, 0.0, f / (2.0 * w**2))
191 self.Hc = np.where(self.rigid, 0.0, (f / 2.0) * H)
192 rhs1 = -(self.qconst - self.Hc.real)
193 rhs2 = (1j * Om * self.Hc).real
194 self.B = (rhs2 - self.s1 * rhs1) / (self.s2 - self.s1)
195 self.A = rhs1 - self.B
196 self.Gc = (f / 2.0) / (2.0 * al + 1j * Om)
197 self.vconst = f / (4.0 * al)
198 self.Ar = -(self.vconst - self.Gc.real)
199 qT, vT, _ = self._pulse_phase(np.array([T]))
200 self.qT, self.vT = qT[:, 0], vT[:, 0]
201 with np.errstate(divide="ignore", invalid="ignore"):
202 self.B2 = (self.vT - self.s1 * self.qT) / (self.s2 - self.s1)
203 self.A2 = self.qT - self.B2
205 # ------------------------------------------------------------ symmetry and eigensolves
206 def symmetric_prolongation(self, dofs: np.ndarray | None = None) -> sp.csc_matrix:
207 """Sparse matrix whose columns span the functions even in y: each column pairs a node with its mirror image
208 (w and w_x equal, w_y and w_xy opposite), and nodes on the axis keep only w and w_x. If dofs is given, the
209 matrix is restricted to those rows (in that order) and to the columns that touch them."""
210 j_axis = self.half
211 gi, gj = self.node_i, self.node_j
212 nn = self.ncell + 1
213 mirror = self.node_map[gi * nn + (2 * j_axis - gj)]
214 rows, cols, vals = [], [], []
215 col = 0
216 signs = (1.0, 1.0, -1.0, -1.0)
217 for n in range(self.n_nodes):
218 if gj[n] > j_axis:
219 nm = mirror[n]
220 for d in range(4):
221 rows += [4 * n + d, 4 * nm + d]
222 cols += [col, col]
223 vals += [1.0, signs[d]]
224 col += 1
225 elif gj[n] == j_axis:
226 for d in range(2):
227 rows.append(4 * n + d)
228 cols.append(col)
229 vals.append(1.0)
230 col += 1
231 P = sp.csc_matrix((vals, (rows, cols)), shape=(self.ndof, col))
232 if dofs is not None:
233 P = P[dofs]
234 P = P[:, np.unique(P.nonzero()[1])]
235 return P.tocsc()
237 def eigenmodes(self, K, M, k: int, sigma: float, dofs: np.ndarray | None = None):
238 """The k lowest eigenpairs of K v = lambda M v, restricted to the dofs given (default all) and to the
239 y-symmetric subspace if the mesh was built with symmetric_y. Sorted ascending, M-normalized, and returned
240 on the full dof set (zero on excluded dofs)."""
241 Kr, Mr = (K, M) if dofs is None else (K[dofs][:, dofs], M[dofs][:, dofs])
242 P = self.symmetric_prolongation(dofs) if self.symmetric_y else None
243 if P is not None:
244 Kr, Mr = P.T @ Kr @ P, P.T @ Mr @ P
245 Kr, Mr = sp.csc_matrix(Kr), sp.csc_matrix(Mr)
246 evals, evecs = eigsh(Kr, k=min(k, Kr.shape[0] - 2), M=Mr, sigma=sigma, which="LM")
247 order = np.argsort(evals)
248 evals, evecs = evals[order], evecs[:, order]
249 if P is not None:
250 evecs = P @ evecs
251 Md = M if dofs is None else M[dofs][:, dofs]
252 evecs = evecs / np.sqrt(np.einsum("ij,ij->j", evecs, Md @ evecs))
253 if dofs is not None:
254 full = np.zeros((self.ndof, evecs.shape[1]))
255 full[dofs] = evecs
256 evecs = full
257 return evals, evecs
259 # ------------------------------------------------------------ assembly
260 def element_matrix(self, Cmat: np.ndarray) -> np.ndarray:
261 """Element matrix of int [w_xx, w_yy, w_xy] Cmat [v_xx, v_yy, v_xy]^T dA (no stiffness factor)."""
262 _, Nxx, Nyy, Nxy, WQ = self._shape
263 Ae = np.zeros((16, 16))
264 for q in range(len(WQ)):
265 B2 = np.stack([Nxx[q], Nyy[q], Nxy[q]])
266 Ae += WQ[q] * (B2.T @ Cmat @ B2)
267 return Ae
269 def assemble(self, Ae: np.ndarray) -> sp.csc_matrix:
270 rows = np.repeat(self.cell_dofs, 16, axis=1).ravel()
271 cols = np.tile(self.cell_dofs, (1, 16)).ravel()
272 vals = np.tile(Ae.ravel(), self.n_cells)
273 A = sp.coo_matrix((vals, (rows, cols)), shape=(self.ndof, self.ndof)).tocsc()
274 A.sum_duplicates()
275 return A
277 # ------------------------------------------------------------ quadrature over the masked domain
278 def quadrature(self, order: int = 3):
279 """Gauss-Legendre points and weights over all active cells: xq, yq, wq."""
280 g, wg = leg.leggauss(order)
281 s = (g + 1) / 2
282 ws = wg / 2
283 XI, ETA = np.meshgrid(s, s, indexing="ij")
284 a = self.a
285 xq = (self.x0 + (self.cells[:, 0][:, None] + XI.ravel()[None, :]) * a).ravel()
286 yq = (self.x0 + (self.cells[:, 1][:, None] + ETA.ravel()[None, :]) * a).ravel()
287 wq = np.tile((ws[:, None] * ws[None, :]).ravel() * a * a, self.n_cells)
288 return xq, yq, wq
290 # ------------------------------------------------------------ evaluation at points
291 def locate(self, x: np.ndarray, y: np.ndarray, tol: float = 1e-9):
292 """Cell id and local coordinates of points. Points inside the disk that fall in an inactive cell of the
293 staircase are assigned the nearest active cell (the element polynomial is extrapolated slightly).
294 Returns cid (-1 for points outside the disk), xi, eta."""
295 x = np.asarray(x, dtype=float)
296 y = np.asarray(y, dtype=float)
297 ix = np.clip(np.floor((x - self.x0) / self.a).astype(int), 0, self.ncell - 1)
298 iy = np.clip(np.floor((y - self.x0) / self.a).astype(int), 0, self.ncell - 1)
299 cid = self.cell_id[ix, iy].copy()
300 inside = x**2 + y**2 <= self.R**2 * (1 + tol)
301 fix = inside & (cid < 0)
302 if fix.any():
303 _, near = self._tree.query(np.stack([x[fix], y[fix]], axis=1))
304 cid[fix] = near
305 cid[~inside] = -1
306 c = np.maximum(cid, 0)
307 cx = self.x0 + self.cells[c, 0] * self.a
308 cy = self.x0 + self.cells[c, 1] * self.a
309 return cid, (x - cx) / self.a, (y - cy) / self.a
311 def _sparse_basis(self, vals: np.ndarray, cid: np.ndarray) -> sp.csr_matrix:
312 npts = len(cid)
313 ok = cid >= 0
314 rows = np.repeat(np.arange(npts)[ok], 16)
315 cols = self.cell_dofs[cid[ok]].ravel()
316 data = vals[ok].ravel()
317 return sp.csr_matrix((data, (rows, cols)), shape=(npts, self.ndof))
319 def basis_at(self, x, y) -> sp.csr_matrix:
320 """Sparse (npts x ndof) matrix of shape-function values; rows of points outside the disk are zero."""
321 cid, xi, eta = self.locate(x, y)
322 N, _, _, _ = shape_functions(xi, eta, self.a)
323 return self._sparse_basis(N, cid)
325 def basis_derivs_at(self, x, y):
326 cid, xi, eta = self.locate(x, y)
327 _, Nxx, Nyy, Nxy = shape_functions(xi, eta, self.a)
328 return self._sparse_basis(Nxx, cid), self._sparse_basis(Nyy, cid), self._sparse_basis(Nxy, cid)
330 def modes_at(self, x, y) -> np.ndarray:
331 return self.basis_at(x, y) @ self.V
333 def inside(self, x, y) -> np.ndarray:
334 return np.asarray(x) ** 2 + np.asarray(y) ** 2 <= self.R**2 * (1 + 1e-9)
336 # ------------------------------------------------------------ modal time response (linear)
337 def _pulse_phase(self, t: np.ndarray):
338 t = t[None, :]
339 Om = self.Om
340 e = np.exp(1j * Om * t)
341 s1, s2, A, B = self.s1[:, None], self.s2[:, None], self.A[:, None], self.B[:, None]
342 E1, E2 = np.exp(s1 * t), np.exp(s2 * t)
343 Hc = self.Hc[:, None]
344 q = (A * E1 + B * E2).real + self.qconst[:, None] - (Hc * e).real
345 v = (A * s1 * E1 + B * s2 * E2).real - (1j * Om * Hc * e).real
346 acc = (A * s1**2 * E1 + B * s2**2 * E2).real + (Om**2 * Hc * e).real
347 rg = self.rigid
348 if rg.any():
349 al = self.alpha[rg][:, None]
350 Ar, Gc, vc = self.Ar[rg][:, None], self.Gc[rg][:, None], self.vconst[rg][:, None]
351 Ed = np.exp(-2.0 * al * t)
352 q[rg] = Ar * (1.0 - Ed) / (2.0 * al) + vc * t - ((e - 1.0) * Gc / (1j * Om)).real
353 v[rg] = Ar * Ed + vc - (Gc * e).real
354 acc[rg] = -2.0 * al * Ar * Ed - (1j * Om * Gc * e).real
355 return q, v, acc
357 def _free_phase(self, tau: np.ndarray):
358 tau = tau[None, :]
359 s1, s2, A, B = self.s1[:, None], self.s2[:, None], self.A2[:, None], self.B2[:, None]
360 E1, E2 = np.exp(s1 * tau), np.exp(s2 * tau)
361 q = (A * E1 + B * E2).real
362 v = (A * s1 * E1 + B * s2 * E2).real
363 acc = (A * s1**2 * E1 + B * s2**2 * E2).real
364 rg = self.rigid
365 if rg.any():
366 al = self.alpha[rg][:, None]
367 qT, vT = self.qT[rg][:, None], self.vT[rg][:, None]
368 Ed = np.exp(-2.0 * al * tau)
369 q[rg] = qT + vT * (1.0 - Ed) / (2.0 * al)
370 v[rg] = vT * Ed
371 acc[rg] = -2.0 * al * vT * Ed
372 return q, v, acc
374 def response(self, t: np.ndarray):
375 q = np.zeros((self.n_modes, len(t)))
376 v = np.zeros_like(q)
377 acc = np.zeros_like(q)
378 m1 = (t >= 0) & (t <= self.T)
379 m2 = t > self.T
380 if m1.any():
381 q[:, m1], v[:, m1], acc[:, m1] = self._pulse_phase(t[m1])
382 if m2.any():
383 q[:, m2], v[:, m2], acc[:, m2] = self._free_phase(t[m2] - self.T)
384 return q, v, acc
386 # ------------------------------------------------------------ Rayleigh integral (linear, exact retardation)
387 def listener_quadrature(self, order: int = 3):
388 if self._listener_q is not None:
389 return self._listener_q
390 p = self.p
391 xq, yq, wq = self.quadrature(order)
392 r = np.sqrt((p["listener_x"] - xq) ** 2 + (p["listener_y"] - yq) ** 2 + p["listener_height"] ** 2)
393 tau = r / p["c_air"]
394 cq = p["rho_air"] / (2 * np.pi) * wq / r
395 return xq, yq, tau, cq
397 def pressure(self, fs: int, nsamp: int) -> np.ndarray:
398 xq, yq, tau, cq = self.listener_quadrature()
399 Phi = self.modes_at(xq, yq) # (nq, n_modes)
400 t = np.arange(nsamp) / fs
401 out = np.zeros(nsamp)
402 t_split = self.T + tau.max()
403 early = (t < t_split) & (t >= tau.min())
404 late = t >= t_split
405 ne = int(early.sum())
406 if ne:
407 chunk = max(1, 4_000_000 // max(1, self.n_modes * ne))
408 for i0 in range(0, len(xq), chunk):
409 sl = slice(i0, i0 + chunk)
410 tt = (t[early][None, :] - tau[sl][:, None]).ravel()
411 _, _, acc = self.response(tt)
412 acc = acc.reshape(self.n_modes, -1, ne)
413 out[early] += np.einsum("q,qn,nqk->k", cq[sl], Phi[sl], acc, optimize=True)
414 if late.any():
415 tl = t[late] - self.T
416 rg = self.rigid
417 G1 = (cq[:, None] * Phi * np.exp(-np.outer(tau, self.s1))).sum(axis=0)
418 G2 = (cq[:, None] * Phi * np.exp(-np.outer(tau, self.s2))).sum(axis=0)
419 c1 = self.A2 * self.s1**2 * G1
420 c2 = self.B2 * self.s2**2 * G2
421 s1, s2 = self.s1.copy(), self.s2.copy()
422 if rg.any():
423 al = self.alpha[rg]
424 Gr = (cq[:, None] * Phi[:, rg] * np.exp(np.outer(tau, 2.0 * al))).sum(axis=0)
425 c1[rg] = -2.0 * al * self.vT[rg] * Gr
426 c2[rg] = 0.0
427 s1[rg] = -2.0 * al
428 s2[rg] = 0.0
429 coef = np.concatenate([c1, c2])
430 rates = np.concatenate([s1, s2])
431 res = np.zeros(len(tl))
432 step = max(1, 20_000_000 // len(coef))
433 for k0 in range(0, len(tl), step):
434 tk = tl[k0:k0 + step]
435 res[k0:k0 + step] = (np.exp(np.outer(tk, rates)) @ coef).real
436 out[late] = res
437 return out