1"""The plate as evaluated on the grid: material, thickness field, stress-free
2shape, and the derived quantities both models need (masses, the strike patch,
3the listener terms). A ``Plate`` is usually produced by ``Design.evaluate``;
4``Plate.uniform`` builds one directly for tests and experiments."""
6from __future__ import annotations
8import math
9from dataclasses import dataclass, replace
10from functools import cached_property
12import numpy as np
14from .constants import RHO_AIR, ROOM_TEMPERATURE, SPEED_OF_SOUND
15from .grid import Grid
16from .material import Material
17from .params import Listener, Strike
21class Patch:
22 """The plate nodes that a strike acts on."""
24 mask: np.ndarray
25 """``(ny, nx)`` bool: the nodes of the patch."""
26 nearest: tuple[int, int]
27 """``(j, i)`` of the plate node nearest the strike point."""
28 mass: float
29 """Total lumped mass of the patch [kg]."""
30 cx: float
31 """Mass-weighted center of the patch [m]."""
32 cy: float
36class Plate:
37 grid: Grid
38 material: Material
39 h: np.ndarray
40 """Thickness at the nodes [m], ``(ny, nx)``, zero outside the plate."""
41 w0: np.ndarray
42 """Stress-free height at the nodes [m], mean plane removed, zero outside the plate."""
43 Lx: float
44 """Bounding box [m]."""
45 Ly: float
46 temperature: float = ROOM_TEMPERATURE
47 """Temperature [K], for the thermoelastic damping."""
48 name: str = "plate"
50 def __post_init__(self) -> None:
51 if self.h.shape != self.grid.shape or self.w0.shape != self.grid.shape:
52 raise ValueError("h and w0 must have the grid's shape (ny, nx)")
53 if (self.h[self.grid.mask] <= 0).any():
54 raise ValueError("the thickness must be positive on the plate")
57 def uniform(cls, grid: Grid, material: Material, thickness: float, profile: np.ndarray | None = None, name: str = "uniform") -> "Plate":
58 """A plate of uniform thickness on a grid, optionally with a stress-free profile (its mean plane is removed)."""
59 h = np.where(grid.mask, thickness, 0.0)
60 w0 = np.zeros(grid.shape) if profile is None else np.where(grid.mask, remove_plane(np.asarray(profile, float), grid), 0.0)
61 return cls(grid, material, h, w0, (grid.nx - 1) * grid.hx, (grid.ny - 1) * grid.hy, name=name)
63 def with_thickness(self, h: np.ndarray) -> "Plate":
64 """The same plate with another thickness field (masked to the plate)."""
65 return replace(self, h=np.where(self.grid.mask, h, 0.0))
67 # ---- masses and geometry ----
70 def node_mass(self) -> np.ndarray:
71 """Lumped mass of each node, rho h dA times its weight [kg]."""
72 return self.material.rho * self.h * self.grid.node_area
75 def mass(self) -> float:
76 """Total mass [kg]."""
77 return float(self.node_mass.sum())
80 def center(self) -> tuple[float, float]:
81 """Mass-weighted center [m]."""
82 m = self.node_mass
83 return float((m * self.grid.X).sum() / self.mass), float((m * self.grid.Y).sum() / self.mass)
86 def inertia(self) -> tuple[float, float]:
87 """Moments of inertia ``(Ix, Iy)`` about the in-plane axes through the center [kg m^2]."""
88 cx, cy = self.center
89 m = self.node_mass
90 return float((m * (self.grid.Y - cy) ** 2).sum()), float((m * (self.grid.X - cx) ** 2).sum())
93 def area(self) -> float:
94 return self.grid.area
97 def mean_thickness(self) -> float:
98 """Area-weighted mean thickness [m]."""
99 w = self.grid.node_weight
100 return float((w * self.h).sum() / w.sum())
103 def thickness_range(self) -> tuple[float, float]:
104 hh = self.h[self.grid.mask]
105 return float(hh.min()), float(hh.max())
108 def cell_thickness(self) -> np.ndarray:
109 """Thickness of each cell, the mean of its four corners, ``(ny - 1, nx - 1)``."""
110 h = self.h
111 return 0.25 * (h[:-1, :-1] + h[:-1, 1:] + h[1:, 1:] + h[1:, :-1])
114 def is_flat(self) -> bool:
115 """Whether the profile vanishes, in which case bending and in-plane motion decouple in the linear model."""
116 return not np.any(self.w0)
119 def profile_stats(self) -> tuple[float, float]:
120 """Peak-to-peak height and largest slope of the profile over the plate nodes."""
121 g, w0, m = self.grid, self.w0, self.grid.mask
122 left, right, down, up = g.neighbor_masks()
123 sx = np.zeros(g.shape)
124 sy = np.zeros(g.shape)
125 wp = np.pad(w0, 1)
126 both = left & right
127 sx[both] = ((wp[1:-1, 2:] - wp[1:-1, :-2]) / (2 * g.hx))[both]
128 only = right & ~left
129 sx[only] = ((wp[1:-1, 2:] - w0) / g.hx)[only]
130 only = left & ~right
131 sx[only] = ((w0 - wp[1:-1, :-2]) / g.hx)[only]
132 both = down & up
133 sy[both] = ((wp[2:, 1:-1] - wp[:-2, 1:-1]) / (2 * g.hy))[both]
134 only = up & ~down
135 sy[only] = ((wp[2:, 1:-1] - w0) / g.hy)[only]
136 only = down & ~up
137 sy[only] = ((w0 - wp[:-2, 1:-1]) / g.hy)[only]
138 v = w0[m]
139 return float(v.max() - v.min()), float(np.hypot(sx, sy)[m].max())
141 def dispersion_frequency(self, fraction: float, h: float | None = None) -> float:
142 """Frequency [Hz] at which grid dispersion lowers a bending partial by
143 ``fraction`` (a plane wave along a grid axis has the discrete frequency
144 ratio (2 - 2 cos kh) / (kh)^2), at the thickness ``h`` (the mean by default)."""
145 lo, hi = 0.0, math.pi
146 for _ in range(60):
147 mid = (lo + hi) / 2
148 ratio = (2 - 2 * math.cos(mid)) / (mid * mid)
149 if 1 - ratio < fraction:
150 lo = mid
151 else:
152 hi = mid
153 k = lo / max(self.grid.hx, self.grid.hy)
154 return self.material.bending_constant(self.mean_thickness if h is None else h) * k * k / (2 * math.pi)
156 # ---- strike and listener ----
158 def strike_point(self, strike: Strike) -> tuple[float, float]:
159 """The strike point in plate coordinates [m]."""
160 return (strike.x - 0.5) * self.Lx, (strike.y - 0.5) * self.Ly
162 def patch_around(self, x: float, y: float, radius: float) -> Patch:
163 """The plate nodes within ``radius`` of ``(x, y)``, always including the nearest one."""
164 g = self.grid
165 d = np.hypot(g.X - x, g.Y - y)
166 d = np.where(self.node_mass > 0, d, np.inf)
167 nearest = np.unravel_index(int(np.argmin(d)), g.shape)
168 mask = d <= radius
169 mask[nearest] = True
170 m = self.node_mass[mask]
171 mass = float(m.sum())
172 return Patch(mask, (int(nearest[0]), int(nearest[1])), mass, float((m * g.X[mask]).sum() / mass), float((m * g.Y[mask]).sum() / mass))
174 def strike_patch(self, strike: Strike) -> Patch:
175 x, y = self.strike_point(strike)
176 return self.patch_around(x, y, strike.radius)
178 def listener_terms(self, listener: Listener) -> tuple[np.ndarray, np.ndarray]:
179 """Per-node terms of the Rayleigh integral: the weight rho_air dA / (2 pi r)
180 and the delay r / c [s] to the listener, ``(ny, nx)`` each, zero outside
181 the plate. The profile is included in the node positions."""
182 g = self.grid
183 lx, ly, lz = listener.position
184 r = np.sqrt((lx - g.X) ** 2 + (ly - g.Y) ** 2 + (lz - self.w0) ** 2)
185 area = g.node_area
186 weight = np.where(area > 0, RHO_AIR * area / (2 * math.pi * np.maximum(r, 1e-3)), 0.0)
187 delay = np.where(area > 0, r / SPEED_OF_SOUND, 0.0)
188 return weight, delay
190 # ---- summary ----
192 def summary(self) -> str:
193 """A few lines describing the plate: box, material, grid, thickness, mass, profile, loss constants."""
194 from .loss import LN1000, coincidence_frequency, estimated_rate, plate_geometry
196 g, m = self.grid, self.material
197 hmin, hmax = self.thickness_range
198 lines = [
199 f"{self.name}: {self.Lx * 1e3:.0f} x {self.Ly * 1e3:.0f} mm box, area {self.area * 1e4:.0f} cm^2, mass {self.mass * 1e3:.0f} g",
200 f"material: E = {m.E / 1e9:.1f} GPa, nu = {m.nu}, rho = {m.rho:.0f} kg/m^3, loss factor {m.loss_factor:g}"
201 + (", thermal constants given" if m.thermal else ", no thermal constants"),
202 f"grid: {g.nx} x {g.ny} nodes at {g.hx * 1e3:.2f} x {g.hy * 1e3:.2f} mm, {g.count} in the plate",
203 f"thickness {hmin * 1e3:.2f} to {hmax * 1e3:.2f} mm (mean {self.mean_thickness * 1e3:.2f}); "
204 f"grid dispersion -3% at {self.dispersion_frequency(0.03):.0f} Hz, -10% at {self.dispersion_frequency(0.10):.0f} Hz",
205 ]
206 if self.is_flat:
207 lines.append("flat (no profile)")
208 else:
209 p2p, slope = self.profile_stats
210 lines.append(f"profile: {p2p * 1e3:.1f} mm peak to peak, largest slope {slope:.3f}" + (" (above the shallow-shell range of about 0.1)" if slope > 0.1 else ""))
211 geom = plate_geometry(self)
212 t60 = [LN1000 / max(estimated_rate(self, geom, f).total, 1e-12) for f in (100, 1000, 5000)]
213 lines.append(
214 f"coincidence frequency {coincidence_frequency(m, self.mean_thickness):.0f} Hz; estimated T60 {t60[0]:.1f} / {t60[1]:.1f} / {t60[2]:.1f} s at 100 Hz / 1 kHz / 5 kHz"
215 )
216 return "\n".join(lines)
218 def __repr__(self) -> str:
219 return f"Plate({self.name!r}, {self.grid.nx}x{self.grid.ny} grid, {self.grid.count} nodes, mass {self.mass * 1e3:.0f} g)"
222def remove_plane(f: np.ndarray, g: Grid) -> np.ndarray:
223 """Subtracts the weighted least-squares plane a + b x + c y from ``f`` over the plate nodes."""
224 w = g.node_weight
225 sel = w > 0
226 if not sel.any():
227 return f.copy()
228 A = np.stack([np.ones(sel.sum()), g.X[sel], g.Y[sel]], axis=1)
229 sw = np.sqrt(w[sel])
230 coef, *_ = np.linalg.lstsq(A * sw[:, None], f[sel] * sw, rcond=None)
231 return f - (coef[0] + coef[1] * g.X + coef[2] * g.Y)