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:
8699301Variable thickness: coefficient assembly of mass, bending and the Airy compatibility operatorJeremy Magland 79 def __init__(self, ch: dict, spacing: float, n_modes: int, symmetric_y: bool = False, max_hz: float | None = None, thickness=None):
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 80 """symmetric_y restricts the basis to functions even in y (w(x, -y) = w(x, y)), which is exact for a strike
8699301Variable thickness: coefficient assembly of mass, bending and the Airy compatibility operatorJeremy Magland 81 centered on the x axis; the grid is symmetric about that axis by construction. max_hz drops modes above it.
82 thickness, if given, is a function h(x, y) (variable thickness); otherwise the uniform parameters["h"]."""
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 83 p = ch["parameters"]
84 self.p = p
8699301Variable thickness: coefficient assembly of mass, bending and the Airy compatibility operatorJeremy Magland 85 R, nu = p["disk_radius"], p["nu"]
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 86 self.R = R
87 a = float(spacing)
88 self.a = a
8699301Variable thickness: coefficient assembly of mass, bending and the Airy compatibility operatorJeremy Magland 89 self.h_at = thickness if thickness is not None else (lambda x, y: np.full(np.shape(x), float(p["h"])))
90 h = float(p["h"]) if "h" in p else float(np.mean(self.h_at(np.zeros(1), np.zeros(1))))
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 91 D = p["E"] * h**3 / (12.0 * (1.0 - nu**2))
92 m = p["rho"] * h
8699301Variable thickness: coefficient assembly of mass, bending and the Airy compatibility operatorJeremy Magland 93 self.D, self.m = D, m # the uniform values (used only when the thickness is uniform)
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 94
95 # ---- the masked grid
96 half = int(np.ceil(R / a)) + 1
97 ncell = 2 * half
98 self.ncell = ncell
99 self.half = half
100 self.x0 = -half * a # coordinate of grid line 0
101 centers = self.x0 + (np.arange(ncell) + 0.5) * a
102 CX, CY = np.meshgrid(centers, centers, indexing="ij")
103 active = CX**2 + CY**2 <= R**2
104 cells = np.argwhere(active) # (nc, 2): ix, iy
105 self.cells = cells
106 self.n_cells = len(cells)
107 self.cell_id = -np.ones((ncell, ncell), dtype=int)
108 self.cell_id[cells[:, 0], cells[:, 1]] = np.arange(self.n_cells)
109 self.cell_center = np.stack([centers[cells[:, 0]], centers[cells[:, 1]]], axis=1)
110 self._tree = cKDTree(self.cell_center)
111 nn = ncell + 1
112 gnode = np.zeros((self.n_cells, 4), dtype=int)
113 for c, (cx, cy) in enumerate(CORNERS):
114 gnode[:, c] = (cells[:, 0] + cx) * nn + (cells[:, 1] + cy)
115 touch = np.bincount(gnode.ravel(), minlength=nn * nn)
116 used = touch > 0
117 self.node_map = -np.ones(nn * nn, dtype=int)
118 self.node_map[used] = np.arange(int(used.sum()))
119 self.n_nodes = int(used.sum())
120 self.ndof = 4 * self.n_nodes
121 cnode = self.node_map[gnode] # (nc, 4) compact node ids
122 gidx = np.flatnonzero(used)
123 self.node_i, self.node_j = gidx // nn, gidx % nn # grid indices of each compact node
124 self.cell_dofs = (4 * cnode[:, :, None] + np.arange(4)[None, None, :]).reshape(self.n_cells, 16)
125 bnodes = self.node_map[used & (touch < 4)]
126 self.boundary_dofs = (4 * bnodes[:, None] + np.arange(4)[None, :]).ravel()
127 self.interior_dofs = np.setdiff1d(np.arange(self.ndof), self.boundary_dofs)
128 self.area = self.n_cells * a * a
130 # ---- element matrices
131 g, wg = leg.leggauss(4)
132 s = (g + 1) / 2
133 ws = wg / 2
134 XI, ETA = np.meshgrid(s, s, indexing="ij")
135 WQ = (ws[:, None] * ws[None, :]).ravel() * a * a
136 N, Nxx, Nyy, Nxy = shape_functions(XI.ravel(), ETA.ravel(), a)
137 self._shape = (N, Nxx, Nyy, Nxy, WQ)
138 C = np.array([[1.0, nu, 0.0], [nu, 1.0, 0.0], [0.0, 0.0, 2.0 * (1.0 - nu)]])
8699301Variable thickness: coefficient assembly of mass, bending and the Airy compatibility operatorJeremy Magland 139 # Gauss-point coordinates of every cell, and the thickness there
140 gx = (self.x0 + (self.cells[:, 0][:, None] + XI.ravel()[None, :]) * a)
141 gy = (self.x0 + (self.cells[:, 1][:, None] + ETA.ravel()[None, :]) * a)
142 self._gauss_xy = (gx, gy)
143 hq = self.h_at(gx, gy) # (n_cells, nq)
144 self.K = self.assemble_coeff(C, p["E"] * hq**3 / (12.0 * (1.0 - nu**2)))
145 self.M1 = self.assemble_mass(np.ones_like(hq)) # unit-density mass matrix
146 self.M = self.assemble_mass(p["rho"] * hq)
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 147
148 # ---- modes: the n_modes lowest, mass-normalized (shift-invert Lanczos around a negative shift)
149 self.symmetric_y = symmetric_y
150 evals, evecs = self.eigenmodes(self.K, self.M, n_modes, sigma=-1.0e3)
151 omega = np.sqrt(np.maximum(evals, 0.0))
152 omega[omega < RIGID_OMEGA] = 0.0
153 if max_hz is not None:
154 keep = omega <= 2 * np.pi * max_hz
155 omega, evecs = omega[keep], evecs[:, keep]
156 self.omega = omega
157 self.V = evecs
158 self.n_modes = len(omega)
159 self.alpha = p["sigma0"] + 0.5 * p["eta"] * omega**2
160 self.rigid = omega == 0.0
161 self.n_rigid = int(self.rigid.sum())
163 # ---- modal forces: f_n = -int p_hat Phi_n dA over the strike disk (polar Gauss quadrature)
164 Rs, xs, ys = p["strike_radius"], p["strike_x"], p["strike_y"]
165 c_phi = 0.5 - 2.0 / np.pi**2
166 gr, gw = leg.leggauss(24)
167 r = Rs * (gr + 1.0) / 2.0
168 wr = gw * Rs / 2.0
169 nth = 64
170 th = 2 * np.pi * np.arange(nth) / nth
171 rr, tt = np.meshgrid(r, th, indexing="ij")
172 ww = (wr[:, None] * rr * (2 * np.pi / nth)).ravel()
173 px = (xs + rr * np.cos(tt)).ravel()
174 py = (ys + rr * np.sin(tt)).ravel()
175 phat = p["strike_peak_force"] / (c_phi * np.pi * Rs**2) * 0.5 * (1.0 + np.cos(np.pi * rr.ravel() / Rs))
176 self.f = -(ww * phat) @ self.modes_at(px, py)
df0973fMarguerre shallow-shell terms for the bowed disk (cymbal-bow-v1)Jeremy Magland 178 self._listener_q = None # optional override of the listener quadrature (used by the shell solver)
179 self._setup_linear_response()
181 def _setup_linear_response(self):
182 """Closed-form constants of the damped modal response to the Hann strike (see the square-plate solver).
183 Depends on omega, alpha, f and rigid; call again after changing them."""
184 p = self.p
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 185 # ---- closed-form linear response constants (see the square-plate solver)
186 T = p["strike_duration"]
187 Om = 2 * np.pi / T
188 self.T, self.Om = T, Om
189 al, w, f = self.alpha, self.omega, self.f
190 with np.errstate(divide="ignore", invalid="ignore"):
191 disc = np.sqrt(al**2 - w**2 + 0j)
192 self.s1 = -al + disc
193 self.s2 = -al - disc
194 H = 1.0 / (w**2 - Om**2 + 2j * al * Om)
195 self.qconst = np.where(self.rigid, 0.0, f / (2.0 * w**2))
196 self.Hc = np.where(self.rigid, 0.0, (f / 2.0) * H)
197 rhs1 = -(self.qconst - self.Hc.real)
198 rhs2 = (1j * Om * self.Hc).real
199 self.B = (rhs2 - self.s1 * rhs1) / (self.s2 - self.s1)
200 self.A = rhs1 - self.B
201 self.Gc = (f / 2.0) / (2.0 * al + 1j * Om)
202 self.vconst = f / (4.0 * al)
203 self.Ar = -(self.vconst - self.Gc.real)
204 qT, vT, _ = self._pulse_phase(np.array([T]))
205 self.qT, self.vT = qT[:, 0], vT[:, 0]
206 with np.errstate(divide="ignore", invalid="ignore"):
207 self.B2 = (self.vT - self.s1 * self.qT) / (self.s2 - self.s1)
208 self.A2 = self.qT - self.B2
210 # ------------------------------------------------------------ symmetry and eigensolves
211 def symmetric_prolongation(self, dofs: np.ndarray | None = None) -> sp.csc_matrix:
212 """Sparse matrix whose columns span the functions even in y: each column pairs a node with its mirror image
213 (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
214 matrix is restricted to those rows (in that order) and to the columns that touch them."""
215 j_axis = self.half
216 gi, gj = self.node_i, self.node_j
217 nn = self.ncell + 1
218 mirror = self.node_map[gi * nn + (2 * j_axis - gj)]
219 rows, cols, vals = [], [], []
220 col = 0
221 signs = (1.0, 1.0, -1.0, -1.0)
222 for n in range(self.n_nodes):
223 if gj[n] > j_axis:
224 nm = mirror[n]
225 for d in range(4):
226 rows += [4 * n + d, 4 * nm + d]
227 cols += [col, col]
228 vals += [1.0, signs[d]]
229 col += 1
230 elif gj[n] == j_axis:
231 for d in range(2):
232 rows.append(4 * n + d)
233 cols.append(col)
234 vals.append(1.0)
235 col += 1
236 P = sp.csc_matrix((vals, (rows, cols)), shape=(self.ndof, col))
237 if dofs is not None:
238 P = P[dofs]
239 P = P[:, np.unique(P.nonzero()[1])]
240 return P.tocsc()
242 def eigenmodes(self, K, M, k: int, sigma: float, dofs: np.ndarray | None = None):
243 """The k lowest eigenpairs of K v = lambda M v, restricted to the dofs given (default all) and to the
244 y-symmetric subspace if the mesh was built with symmetric_y. Sorted ascending, M-normalized, and returned
245 on the full dof set (zero on excluded dofs)."""
246 Kr, Mr = (K, M) if dofs is None else (K[dofs][:, dofs], M[dofs][:, dofs])
247 P = self.symmetric_prolongation(dofs) if self.symmetric_y else None
248 if P is not None:
249 Kr, Mr = P.T @ Kr @ P, P.T @ Mr @ P
250 Kr, Mr = sp.csc_matrix(Kr), sp.csc_matrix(Mr)
251 evals, evecs = eigsh(Kr, k=min(k, Kr.shape[0] - 2), M=Mr, sigma=sigma, which="LM")
252 order = np.argsort(evals)
253 evals, evecs = evals[order], evecs[:, order]
254 if P is not None:
255 evecs = P @ evecs
256 Md = M if dofs is None else M[dofs][:, dofs]
257 evecs = evecs / np.sqrt(np.einsum("ij,ij->j", evecs, Md @ evecs))
258 if dofs is not None:
259 full = np.zeros((self.ndof, evecs.shape[1]))
260 full[dofs] = evecs
261 evecs = full
262 return evals, evecs
264 # ------------------------------------------------------------ assembly
265 def element_matrix(self, Cmat: np.ndarray) -> np.ndarray:
266 """Element matrix of int [w_xx, w_yy, w_xy] Cmat [v_xx, v_yy, v_xy]^T dA (no stiffness factor)."""
267 _, Nxx, Nyy, Nxy, WQ = self._shape
268 Ae = np.zeros((16, 16))
269 for q in range(len(WQ)):
270 B2 = np.stack([Nxx[q], Nyy[q], Nxy[q]])
271 Ae += WQ[q] * (B2.T @ Cmat @ B2)
272 return Ae
8699301Variable thickness: coefficient assembly of mass, bending and the Airy compatibility operatorJeremy Magland 274 def assemble_coeff(self, Cmat: np.ndarray, coef: np.ndarray) -> sp.csc_matrix:
275 """Assemble int coef(x) [w_xx, w_yy, w_xy] Cmat [v_xx, v_yy, v_xy]^T dA with a coefficient given at the
276 Gauss points of every cell (coef of shape (n_cells, nq))."""
277 _, Nxx, Nyy, Nxy, WQ = self._shape
278 B = np.stack([Nxx, Nyy, Nxy], axis=1) # (nq, 3, 16)
279 BCB = np.einsum("qai,ab,qbj->qij", B, Cmat, B) # (nq, 16, 16)
280 Ke = np.einsum("eq,qij->eij", coef * WQ[None, :], BCB)
281 return self._assemble_all(Ke)
283 def assemble_mass(self, coef: np.ndarray) -> sp.csc_matrix:
284 """Assemble int coef(x) w v dA with a Gauss-point coefficient of shape (n_cells, nq)."""
285 N, _, _, _, WQ = self._shape
286 NN = np.einsum("qi,qj->qij", N, N)
287 Me = np.einsum("eq,qij->eij", coef * WQ[None, :], NN)
288 return self._assemble_all(Me)
290 def _assemble_all(self, Ae_all: np.ndarray) -> sp.csc_matrix:
291 rows = np.repeat(self.cell_dofs, 16, axis=1).ravel()
292 cols = np.tile(self.cell_dofs, (1, 16)).ravel()
293 A = sp.coo_matrix((Ae_all.ravel(), (rows, cols)), shape=(self.ndof, self.ndof)).tocsc()
294 A.sum_duplicates()
295 return A
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 297 def assemble(self, Ae: np.ndarray) -> sp.csc_matrix:
298 rows = np.repeat(self.cell_dofs, 16, axis=1).ravel()
299 cols = np.tile(self.cell_dofs, (1, 16)).ravel()
300 vals = np.tile(Ae.ravel(), self.n_cells)
301 A = sp.coo_matrix((vals, (rows, cols)), shape=(self.ndof, self.ndof)).tocsc()
302 A.sum_duplicates()
303 return A
305 # ------------------------------------------------------------ quadrature over the masked domain
306 def quadrature(self, order: int = 3):
307 """Gauss-Legendre points and weights over all active cells: xq, yq, wq."""
308 g, wg = leg.leggauss(order)
309 s = (g + 1) / 2
310 ws = wg / 2
311 XI, ETA = np.meshgrid(s, s, indexing="ij")
312 a = self.a
313 xq = (self.x0 + (self.cells[:, 0][:, None] + XI.ravel()[None, :]) * a).ravel()
314 yq = (self.x0 + (self.cells[:, 1][:, None] + ETA.ravel()[None, :]) * a).ravel()
315 wq = np.tile((ws[:, None] * ws[None, :]).ravel() * a * a, self.n_cells)
316 return xq, yq, wq
318 # ------------------------------------------------------------ evaluation at points
319 def locate(self, x: np.ndarray, y: np.ndarray, tol: float = 1e-9):
320 """Cell id and local coordinates of points. Points inside the disk that fall in an inactive cell of the
321 staircase are assigned the nearest active cell (the element polynomial is extrapolated slightly).
322 Returns cid (-1 for points outside the disk), xi, eta."""
323 x = np.asarray(x, dtype=float)
324 y = np.asarray(y, dtype=float)
325 ix = np.clip(np.floor((x - self.x0) / self.a).astype(int), 0, self.ncell - 1)
326 iy = np.clip(np.floor((y - self.x0) / self.a).astype(int), 0, self.ncell - 1)
327 cid = self.cell_id[ix, iy].copy()
328 inside = x**2 + y**2 <= self.R**2 * (1 + tol)
329 fix = inside & (cid < 0)
330 if fix.any():
331 _, near = self._tree.query(np.stack([x[fix], y[fix]], axis=1))
332 cid[fix] = near
333 cid[~inside] = -1
334 c = np.maximum(cid, 0)
335 cx = self.x0 + self.cells[c, 0] * self.a
336 cy = self.x0 + self.cells[c, 1] * self.a
337 return cid, (x - cx) / self.a, (y - cy) / self.a
339 def _sparse_basis(self, vals: np.ndarray, cid: np.ndarray) -> sp.csr_matrix:
340 npts = len(cid)
341 ok = cid >= 0
342 rows = np.repeat(np.arange(npts)[ok], 16)
343 cols = self.cell_dofs[cid[ok]].ravel()
344 data = vals[ok].ravel()
345 return sp.csr_matrix((data, (rows, cols)), shape=(npts, self.ndof))
347 def basis_at(self, x, y) -> sp.csr_matrix:
348 """Sparse (npts x ndof) matrix of shape-function values; rows of points outside the disk are zero."""
349 cid, xi, eta = self.locate(x, y)
350 N, _, _, _ = shape_functions(xi, eta, self.a)
351 return self._sparse_basis(N, cid)
353 def basis_derivs_at(self, x, y):
354 cid, xi, eta = self.locate(x, y)
355 _, Nxx, Nyy, Nxy = shape_functions(xi, eta, self.a)
356 return self._sparse_basis(Nxx, cid), self._sparse_basis(Nyy, cid), self._sparse_basis(Nxy, cid)
358 def modes_at(self, x, y) -> np.ndarray:
359 return self.basis_at(x, y) @ self.V
361 def inside(self, x, y) -> np.ndarray:
362 return np.asarray(x) ** 2 + np.asarray(y) ** 2 <= self.R**2 * (1 + 1e-9)
364 # ------------------------------------------------------------ modal time response (linear)
365 def _pulse_phase(self, t: np.ndarray):
366 t = t[None, :]
367 Om = self.Om
368 e = np.exp(1j * Om * t)
369 s1, s2, A, B = self.s1[:, None], self.s2[:, None], self.A[:, None], self.B[:, None]
370 E1, E2 = np.exp(s1 * t), np.exp(s2 * t)
371 Hc = self.Hc[:, None]
372 q = (A * E1 + B * E2).real + self.qconst[:, None] - (Hc * e).real
373 v = (A * s1 * E1 + B * s2 * E2).real - (1j * Om * Hc * e).real
374 acc = (A * s1**2 * E1 + B * s2**2 * E2).real + (Om**2 * Hc * e).real
375 rg = self.rigid
376 if rg.any():
377 al = self.alpha[rg][:, None]
378 Ar, Gc, vc = self.Ar[rg][:, None], self.Gc[rg][:, None], self.vconst[rg][:, None]
379 Ed = np.exp(-2.0 * al * t)
380 q[rg] = Ar * (1.0 - Ed) / (2.0 * al) + vc * t - ((e - 1.0) * Gc / (1j * Om)).real
381 v[rg] = Ar * Ed + vc - (Gc * e).real
382 acc[rg] = -2.0 * al * Ar * Ed - (1j * Om * Gc * e).real
383 return q, v, acc
385 def _free_phase(self, tau: np.ndarray):
386 tau = tau[None, :]
387 s1, s2, A, B = self.s1[:, None], self.s2[:, None], self.A2[:, None], self.B2[:, None]
388 E1, E2 = np.exp(s1 * tau), np.exp(s2 * tau)
389 q = (A * E1 + B * E2).real
390 v = (A * s1 * E1 + B * s2 * E2).real
391 acc = (A * s1**2 * E1 + B * s2**2 * E2).real
392 rg = self.rigid
393 if rg.any():
394 al = self.alpha[rg][:, None]
395 qT, vT = self.qT[rg][:, None], self.vT[rg][:, None]
396 Ed = np.exp(-2.0 * al * tau)
397 q[rg] = qT + vT * (1.0 - Ed) / (2.0 * al)
398 v[rg] = vT * Ed
399 acc[rg] = -2.0 * al * vT * Ed
400 return q, v, acc
402 def response(self, t: np.ndarray):
403 q = np.zeros((self.n_modes, len(t)))
404 v = np.zeros_like(q)
405 acc = np.zeros_like(q)
406 m1 = (t >= 0) & (t <= self.T)
407 m2 = t > self.T
408 if m1.any():
409 q[:, m1], v[:, m1], acc[:, m1] = self._pulse_phase(t[m1])
410 if m2.any():
411 q[:, m2], v[:, m2], acc[:, m2] = self._free_phase(t[m2] - self.T)
412 return q, v, acc
414 # ------------------------------------------------------------ Rayleigh integral (linear, exact retardation)
415 def listener_quadrature(self, order: int = 3):
df0973fMarguerre shallow-shell terms for the bowed disk (cymbal-bow-v1)Jeremy Magland 416 if self._listener_q is not None:
417 return self._listener_q
af2a046Von Karman modal solver on a Cartesian grid masked to the disk, site, and workflowJeremy Magland 418 p = self.p
419 xq, yq, wq = self.quadrature(order)
420 r = np.sqrt((p["listener_x"] - xq) ** 2 + (p["listener_y"] - yq) ** 2 + p["listener_height"] ** 2)
421 tau = r / p["c_air"]
422 cq = p["rho_air"] / (2 * np.pi) * wq / r
423 return xq, yq, tau, cq
425 def pressure(self, fs: int, nsamp: int) -> np.ndarray:
426 xq, yq, tau, cq = self.listener_quadrature()
427 Phi = self.modes_at(xq, yq) # (nq, n_modes)
428 t = np.arange(nsamp) / fs
429 out = np.zeros(nsamp)
430 t_split = self.T + tau.max()
431 early = (t < t_split) & (t >= tau.min())
432 late = t >= t_split
433 ne = int(early.sum())
434 if ne:
435 chunk = max(1, 4_000_000 // max(1, self.n_modes * ne))
436 for i0 in range(0, len(xq), chunk):
437 sl = slice(i0, i0 + chunk)
438 tt = (t[early][None, :] - tau[sl][:, None]).ravel()
439 _, _, acc = self.response(tt)
440 acc = acc.reshape(self.n_modes, -1, ne)
441 out[early] += np.einsum("q,qn,nqk->k", cq[sl], Phi[sl], acc, optimize=True)
442 if late.any():
443 tl = t[late] - self.T
444 rg = self.rigid
445 G1 = (cq[:, None] * Phi * np.exp(-np.outer(tau, self.s1))).sum(axis=0)
446 G2 = (cq[:, None] * Phi * np.exp(-np.outer(tau, self.s2))).sum(axis=0)
447 c1 = self.A2 * self.s1**2 * G1
448 c2 = self.B2 * self.s2**2 * G2
449 s1, s2 = self.s1.copy(), self.s2.copy()
450 if rg.any():
451 al = self.alpha[rg]
452 Gr = (cq[:, None] * Phi[:, rg] * np.exp(np.outer(tau, 2.0 * al))).sum(axis=0)
453 c1[rg] = -2.0 * al * self.vT[rg] * Gr
454 c2[rg] = 0.0
455 s1[rg] = -2.0 * al
456 s2[rg] = 0.0
457 coef = np.concatenate([c1, c2])
458 rates = np.concatenate([s1, s2])
459 res = np.zeros(len(tl))
460 step = max(1, 20_000_000 // len(coef))
461 for k0 in range(0, len(tl), step):
462 tk = tl[k0:k0 + step]
463 res[k0:k0 + step] = (np.exp(np.outer(tk, rates)) @ coef).real
464 out[late] = res
465 return out