af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 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 # ---- closed-form linear response constants (see the square-plate solver)
174 T = p["strike_duration"]
175 Om = 2 * np.pi / T
176 self.T, self.Om = T, Om
177 al, w, f = self.alpha, self.omega, self.f
178 with np.errstate(divide="ignore", invalid="ignore"):
179 disc = np.sqrt(al**2 - w**2 + 0j)
180 self.s1 = -al + disc
181 self.s2 = -al - disc
182 H = 1.0 / (w**2 - Om**2 + 2j * al * Om)
183 self.qconst = np.where(self.rigid, 0.0, f / (2.0 * w**2))
184 self.Hc = np.where(self.rigid, 0.0, (f / 2.0) * H)
185 rhs1 = -(self.qconst - self.Hc.real)
186 rhs2 = (1j * Om * self.Hc).real
187 self.B = (rhs2 - self.s1 * rhs1) / (self.s2 - self.s1)
188 self.A = rhs1 - self.B
189 self.Gc = (f / 2.0) / (2.0 * al + 1j * Om)
190 self.vconst = f / (4.0 * al)
191 self.Ar = -(self.vconst - self.Gc.real)
192 qT, vT, _ = self._pulse_phase(np.array([T]))
193 self.qT, self.vT = qT[:, 0], vT[:, 0]
194 with np.errstate(divide="ignore", invalid="ignore"):
195 self.B2 = (self.vT - self.s1 * self.qT) / (self.s2 - self.s1)
196 self.A2 = self.qT - self.B2
198 # ------------------------------------------------------------ symmetry and eigensolves
199 def symmetric_prolongation(self, dofs: np.ndarray | None = None) -> sp.csc_matrix:
200 """Sparse matrix whose columns span the functions even in y: each column pairs a node with its mirror image
201 (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
202 matrix is restricted to those rows (in that order) and to the columns that touch them."""
203 j_axis = self.half
204 gi, gj = self.node_i, self.node_j
205 nn = self.ncell + 1
206 mirror = self.node_map[gi * nn + (2 * j_axis - gj)]
207 rows, cols, vals = [], [], []
208 col = 0
209 signs = (1.0, 1.0, -1.0, -1.0)
210 for n in range(self.n_nodes):
211 if gj[n] > j_axis:
212 nm = mirror[n]
213 for d in range(4):
214 rows += [4 * n + d, 4 * nm + d]
215 cols += [col, col]
216 vals += [1.0, signs[d]]
217 col += 1
218 elif gj[n] == j_axis:
219 for d in range(2):
220 rows.append(4 * n + d)
221 cols.append(col)
222 vals.append(1.0)
223 col += 1
224 P = sp.csc_matrix((vals, (rows, cols)), shape=(self.ndof, col))
225 if dofs is not None:
226 P = P[dofs]
227 P = P[:, np.unique(P.nonzero()[1])]
228 return P.tocsc()
230 def eigenmodes(self, K, M, k: int, sigma: float, dofs: np.ndarray | None = None):
231 """The k lowest eigenpairs of K v = lambda M v, restricted to the dofs given (default all) and to the
232 y-symmetric subspace if the mesh was built with symmetric_y. Sorted ascending, M-normalized, and returned
233 on the full dof set (zero on excluded dofs)."""
234 Kr, Mr = (K, M) if dofs is None else (K[dofs][:, dofs], M[dofs][:, dofs])
235 P = self.symmetric_prolongation(dofs) if self.symmetric_y else None
236 if P is not None:
237 Kr, Mr = P.T @ Kr @ P, P.T @ Mr @ P
238 Kr, Mr = sp.csc_matrix(Kr), sp.csc_matrix(Mr)
239 evals, evecs = eigsh(Kr, k=min(k, Kr.shape[0] - 2), M=Mr, sigma=sigma, which="LM")
240 order = np.argsort(evals)
241 evals, evecs = evals[order], evecs[:, order]
242 if P is not None:
243 evecs = P @ evecs
244 Md = M if dofs is None else M[dofs][:, dofs]
245 evecs = evecs / np.sqrt(np.einsum("ij,ij->j", evecs, Md @ evecs))
246 if dofs is not None:
247 full = np.zeros((self.ndof, evecs.shape[1]))
248 full[dofs] = evecs
249 evecs = full
250 return evals, evecs
252 # ------------------------------------------------------------ assembly
253 def element_matrix(self, Cmat: np.ndarray) -> np.ndarray:
254 """Element matrix of int [w_xx, w_yy, w_xy] Cmat [v_xx, v_yy, v_xy]^T dA (no stiffness factor)."""
255 _, Nxx, Nyy, Nxy, WQ = self._shape
256 Ae = np.zeros((16, 16))
257 for q in range(len(WQ)):
258 B2 = np.stack([Nxx[q], Nyy[q], Nxy[q]])
259 Ae += WQ[q] * (B2.T @ Cmat @ B2)
260 return Ae
262 def assemble(self, Ae: np.ndarray) -> sp.csc_matrix:
263 rows = np.repeat(self.cell_dofs, 16, axis=1).ravel()
264 cols = np.tile(self.cell_dofs, (1, 16)).ravel()
265 vals = np.tile(Ae.ravel(), self.n_cells)
266 A = sp.coo_matrix((vals, (rows, cols)), shape=(self.ndof, self.ndof)).tocsc()
267 A.sum_duplicates()
268 return A
270 # ------------------------------------------------------------ quadrature over the masked domain
271 def quadrature(self, order: int = 3):
272 """Gauss-Legendre points and weights over all active cells: xq, yq, wq."""
273 g, wg = leg.leggauss(order)
274 s = (g + 1) / 2
275 ws = wg / 2
276 XI, ETA = np.meshgrid(s, s, indexing="ij")
277 a = self.a
278 xq = (self.x0 + (self.cells[:, 0][:, None] + XI.ravel()[None, :]) * a).ravel()
279 yq = (self.x0 + (self.cells[:, 1][:, None] + ETA.ravel()[None, :]) * a).ravel()
280 wq = np.tile((ws[:, None] * ws[None, :]).ravel() * a * a, self.n_cells)
281 return xq, yq, wq
283 # ------------------------------------------------------------ evaluation at points
284 def locate(self, x: np.ndarray, y: np.ndarray, tol: float = 1e-9):
285 """Cell id and local coordinates of points. Points inside the disk that fall in an inactive cell of the
286 staircase are assigned the nearest active cell (the element polynomial is extrapolated slightly).
287 Returns cid (-1 for points outside the disk), xi, eta."""
288 x = np.asarray(x, dtype=float)
289 y = np.asarray(y, dtype=float)
290 ix = np.clip(np.floor((x - self.x0) / self.a).astype(int), 0, self.ncell - 1)
291 iy = np.clip(np.floor((y - self.x0) / self.a).astype(int), 0, self.ncell - 1)
292 cid = self.cell_id[ix, iy].copy()
293 inside = x**2 + y**2 <= self.R**2 * (1 + tol)
294 fix = inside & (cid < 0)
295 if fix.any():
296 _, near = self._tree.query(np.stack([x[fix], y[fix]], axis=1))
297 cid[fix] = near
298 cid[~inside] = -1
299 c = np.maximum(cid, 0)
300 cx = self.x0 + self.cells[c, 0] * self.a
301 cy = self.x0 + self.cells[c, 1] * self.a
302 return cid, (x - cx) / self.a, (y - cy) / self.a
304 def _sparse_basis(self, vals: np.ndarray, cid: np.ndarray) -> sp.csr_matrix:
305 npts = len(cid)
306 ok = cid >= 0
307 rows = np.repeat(np.arange(npts)[ok], 16)
308 cols = self.cell_dofs[cid[ok]].ravel()
309 data = vals[ok].ravel()
310 return sp.csr_matrix((data, (rows, cols)), shape=(npts, self.ndof))
312 def basis_at(self, x, y) -> sp.csr_matrix:
313 """Sparse (npts x ndof) matrix of shape-function values; rows of points outside the disk are zero."""
314 cid, xi, eta = self.locate(x, y)
315 N, _, _, _ = shape_functions(xi, eta, self.a)
316 return self._sparse_basis(N, cid)
318 def basis_derivs_at(self, x, y):
319 cid, xi, eta = self.locate(x, y)
320 _, Nxx, Nyy, Nxy = shape_functions(xi, eta, self.a)
321 return self._sparse_basis(Nxx, cid), self._sparse_basis(Nyy, cid), self._sparse_basis(Nxy, cid)
323 def modes_at(self, x, y) -> np.ndarray:
324 return self.basis_at(x, y) @ self.V
326 def inside(self, x, y) -> np.ndarray:
327 return np.asarray(x) ** 2 + np.asarray(y) ** 2 <= self.R**2 * (1 + 1e-9)
329 # ------------------------------------------------------------ modal time response (linear)
330 def _pulse_phase(self, t: np.ndarray):
331 t = t[None, :]
332 Om = self.Om
333 e = np.exp(1j * Om * t)
334 s1, s2, A, B = self.s1[:, None], self.s2[:, None], self.A[:, None], self.B[:, None]
335 E1, E2 = np.exp(s1 * t), np.exp(s2 * t)
336 Hc = self.Hc[:, None]
337 q = (A * E1 + B * E2).real + self.qconst[:, None] - (Hc * e).real
338 v = (A * s1 * E1 + B * s2 * E2).real - (1j * Om * Hc * e).real
339 acc = (A * s1**2 * E1 + B * s2**2 * E2).real + (Om**2 * Hc * e).real
340 rg = self.rigid
341 if rg.any():
342 al = self.alpha[rg][:, None]
343 Ar, Gc, vc = self.Ar[rg][:, None], self.Gc[rg][:, None], self.vconst[rg][:, None]
344 Ed = np.exp(-2.0 * al * t)
345 q[rg] = Ar * (1.0 - Ed) / (2.0 * al) + vc * t - ((e - 1.0) * Gc / (1j * Om)).real
346 v[rg] = Ar * Ed + vc - (Gc * e).real
347 acc[rg] = -2.0 * al * Ar * Ed - (1j * Om * Gc * e).real
348 return q, v, acc
350 def _free_phase(self, tau: np.ndarray):
351 tau = tau[None, :]
352 s1, s2, A, B = self.s1[:, None], self.s2[:, None], self.A2[:, None], self.B2[:, None]
353 E1, E2 = np.exp(s1 * tau), np.exp(s2 * tau)
354 q = (A * E1 + B * E2).real
355 v = (A * s1 * E1 + B * s2 * E2).real
356 acc = (A * s1**2 * E1 + B * s2**2 * E2).real
357 rg = self.rigid
358 if rg.any():
359 al = self.alpha[rg][:, None]
360 qT, vT = self.qT[rg][:, None], self.vT[rg][:, None]
361 Ed = np.exp(-2.0 * al * tau)
362 q[rg] = qT + vT * (1.0 - Ed) / (2.0 * al)
363 v[rg] = vT * Ed
364 acc[rg] = -2.0 * al * vT * Ed
365 return q, v, acc
367 def response(self, t: np.ndarray):
368 q = np.zeros((self.n_modes, len(t)))
369 v = np.zeros_like(q)
370 acc = np.zeros_like(q)
371 m1 = (t >= 0) & (t <= self.T)
372 m2 = t > self.T
373 if m1.any():
374 q[:, m1], v[:, m1], acc[:, m1] = self._pulse_phase(t[m1])
375 if m2.any():
376 q[:, m2], v[:, m2], acc[:, m2] = self._free_phase(t[m2] - self.T)
377 return q, v, acc
379 # ------------------------------------------------------------ Rayleigh integral (linear, exact retardation)
380 def listener_quadrature(self, order: int = 3):
381 p = self.p
382 xq, yq, wq = self.quadrature(order)
383 r = np.sqrt((p["listener_x"] - xq) ** 2 + (p["listener_y"] - yq) ** 2 + p["listener_height"] ** 2)
384 tau = r / p["c_air"]
385 cq = p["rho_air"] / (2 * np.pi) * wq / r
386 return xq, yq, tau, cq
388 def pressure(self, fs: int, nsamp: int) -> np.ndarray:
389 xq, yq, tau, cq = self.listener_quadrature()
390 Phi = self.modes_at(xq, yq) # (nq, n_modes)
391 t = np.arange(nsamp) / fs
392 out = np.zeros(nsamp)
393 t_split = self.T + tau.max()
394 early = (t < t_split) & (t >= tau.min())
395 late = t >= t_split
396 ne = int(early.sum())
397 if ne:
398 chunk = max(1, 4_000_000 // max(1, self.n_modes * ne))
399 for i0 in range(0, len(xq), chunk):
400 sl = slice(i0, i0 + chunk)
401 tt = (t[early][None, :] - tau[sl][:, None]).ravel()
402 _, _, acc = self.response(tt)
403 acc = acc.reshape(self.n_modes, -1, ne)
404 out[early] += np.einsum("q,qn,nqk->k", cq[sl], Phi[sl], acc, optimize=True)
405 if late.any():
406 tl = t[late] - self.T
407 rg = self.rigid
408 G1 = (cq[:, None] * Phi * np.exp(-np.outer(tau, self.s1))).sum(axis=0)
409 G2 = (cq[:, None] * Phi * np.exp(-np.outer(tau, self.s2))).sum(axis=0)
410 c1 = self.A2 * self.s1**2 * G1
411 c2 = self.B2 * self.s2**2 * G2
412 s1, s2 = self.s1.copy(), self.s2.copy()
413 if rg.any():
414 al = self.alpha[rg]
415 Gr = (cq[:, None] * Phi[:, rg] * np.exp(np.outer(tau, 2.0 * al))).sum(axis=0)
416 c1[rg] = -2.0 * al * self.vT[rg] * Gr
417 c2[rg] = 0.0
418 s1[rg] = -2.0 * al
419 s2[rg] = 0.0
420 coef = np.concatenate([c1, c2])
421 rates = np.concatenate([s1, s2])
422 res = np.zeros(len(tl))
423 step = max(1, 20_000_000 // len(coef))
424 for k0 in range(0, len(tl), step):
425 tk = tl[k0:k0 + step]
426 res[k0:k0 + step] = (np.exp(np.outer(tk, rates)) @ coef).real
427 out[late] = res
428 return out