Variable thickness: coefficient assembly of mass, bending and the Airy compatibility operator
2 changed files+72−22
solvers/vonkarman-modal-disk/bfs_disk.pymodified+38−10View file
@@ -76,18 +76,21 @@ def shape_functions(xi: np.ndarray, eta: np.ndarray, a: float):
7676
7777
7878 class BFSDisk:
79- def __init__(self, ch: dict, spacing: float, n_modes: int, symmetric_y: bool = False, max_hz: float | None = None):
79+ def __init__(self, ch: dict, spacing: float, n_modes: int, symmetric_y: bool = False, max_hz: float | None = None, thickness=None):
8080 """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."""
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"]."""
8283 p = ch["parameters"]
8384 self.p = p
84- R, h, nu = p["disk_radius"], p["h"], p["nu"]
85+ R, nu = p["disk_radius"], p["nu"]
8586 self.R = R
8687 a = float(spacing)
8788 self.a = a
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))))
8891 D = p["E"] * h**3 / (12.0 * (1.0 - nu**2))
8992 m = p["rho"] * h
90- self.D, self.m = D, m
93+ self.D, self.m = D, m # the uniform values (used only when the thickness is uniform)
9194
9295 # ---- the masked grid
9396 half = int(np.ceil(R / a)) + 1
@@ -133,12 +136,14 @@ class BFSDisk:
133136 N, Nxx, Nyy, Nxy = shape_functions(XI.ravel(), ETA.ravel(), a)
134137 self._shape = (N, Nxx, Nyy, Nxy, WQ)
135138 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)
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)
142147
143148 # ---- modes: the n_modes lowest, mass-normalized (shift-invert Lanczos around a negative shift)
144149 self.symmetric_y = symmetric_y
@@ -266,6 +271,29 @@ class BFSDisk:
266271 Ae += WQ[q] * (B2.T @ Cmat @ B2)
267272 return Ae
268273
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)
282+
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)
289+
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
296+
269297 def assemble(self, Ae: np.ndarray) -> sp.csc_matrix:
270298 rows = np.repeat(self.cell_dofs, 16, axis=1).ravel()
271299 cols = np.tile(self.cell_dofs, (1, 16)).ravel()
solvers/vonkarman-modal-disk/vk_modal_disk.pymodified+34−12View file
@@ -29,7 +29,9 @@ uses the unrestricted bases instead.
2929
3030 If the challenge has a shell-of-revolution geometry (cymbal-bow-v1: a paraboloidal bow z = Z(r);
3131 cymbal-bell-v1: the bow plus a cosine-dome bell; cymbal-hammered-v1: with raised-cosine dents, which break the
32-symmetry in y so the full modal basis is used), the same solver applies the Marguerre shallow-shell equations, which add the initial shape w0 = Z to the von Karman
32+symmetry in y so the full modal basis is used; cymbal-tapered-v1: with a thickness that varies along the
33+meridian, so the mass, the bending stiffness and the membrane compliance 1/(E h) become variable coefficients
34+in the assembled operators), the same solver applies the Marguerre shallow-shell equations, which add the initial shape w0 = Z to the von Karman
3335 system: L(w + w0, Phi) on the right of the bending equation and 2 L(w0, w) in the compatibility equation.
3436 With G^p_j = int Psi_p L(w0, Phi_j) dA this gives a linear membrane stiffening 2 sum_p |c_p| G^p (G^p)^T,
3537 which is diagonalized once in the flat-mode basis to give the shell modes (exact in the Gautschi step,
@@ -68,10 +70,23 @@ class Profile:
6870 bell of cymbal-bell-v1. Uniform thickness along the normal. The hole of the bell challenge is not represented
6971 (the solver closes it; see the module docstring)."""
7072
71- def __init__(self, R: float, H: float, h: float, bell_radius: float = 0.0, bell_height: float = 0.0, hole_radius: float = 0.0, dents=()):
73+ def __init__(self, R: float, H: float, h: float, bell_radius: float = 0.0, bell_height: float = 0.0, hole_radius: float = 0.0, dents=(), thickness_center=None, thickness_rim=None):
7274 self.R, self.H, self.h = R, H, h
7375 self.rb, self.hb, self.rh = bell_radius, bell_height, hole_radius
7476 self.dents = [(float(d["x"]), float(d["y"]), float(d["radius"]), float(d["depth"])) for d in dents]
77+ self.hc, self.he = thickness_center, thickness_rim # tapered thickness (cymbal-tapered-v1), None for uniform
78+
79+ @property
80+ def tapered(self) -> bool:
81+ return self.hc is not None
82+
83+ def h_at(self, x, y):
84+ """Thickness at (x, y): uniform, or h_c over the bell and linear to h_e at the rim."""
85+ x = np.asarray(x, dtype=float); y = np.asarray(y, dtype=float)
86+ if not self.tapered:
87+ return np.full(np.shape(x), self.h)
88+ r = np.sqrt(x * x + y * y)
89+ return np.where(r <= self.rb, self.hc, self.hc + (self.he - self.hc) * (r - self.rb) / (self.R - self.rb))
7590
7691 @property
7792 def symmetric_y(self) -> bool:
@@ -84,7 +99,8 @@ class Profile:
8499 if not g or g.get("type") != "shell-of-revolution":
85100 return None
86101 p = ch["parameters"]
87- return Profile(p["disk_radius"], p["bow_height"], p["h"], p.get("bell_radius", 0.0), p.get("bell_height", 0.0), p.get("hole_radius", 0.0), g.get("dents", ()))
102+ h = p.get("h", p.get("thickness_center"))
103+ return Profile(p["disk_radius"], p["bow_height"], h, p.get("bell_radius", 0.0), p.get("bell_height", 0.0), p.get("hole_radius", 0.0), g.get("dents", ()), p.get("thickness_center"), p.get("thickness_rim"))
88104
89105 # ---- radial profile and its derivatives
90106 def _bell(self, r):
@@ -160,9 +176,10 @@ class Profile:
160176 return -zx * nz, -zy * nz, nz
161177
162178 def top(self, x, y):
163- """The point of the top face on the normal through the mid-surface point (x, y)."""
179+ """The point of the top face on the normal through the mid-surface point (x, y), with the local thickness."""
164180 nx, ny, nz = self.normal(x, y)
165- return x + 0.5 * self.h * nx, y + 0.5 * self.h * ny, self.z(x, y) + 0.5 * self.h * nz
181+ hh = 0.5 * self.h_at(x, y)
182+ return x + hh * nx, y + hh * ny, self.z(x, y) + hh * nz
166183
167184 def max_slope(self):
168185 g = np.linspace(-self.R, self.R, 801)
@@ -203,7 +220,8 @@ class VonKarmanDisk:
203220 raise ValueError("the dented surface is not symmetric in y; the symmetric basis cannot be used")
204221 self.symmetric_y = symmetric_y
205222 t0 = time.perf_counter()
206- self.plate = BFSDisk(ch2, spacing, n_w, symmetric_y=symmetric_y) # linear modes, modal forces, alpha, omega
223+ thick = self.profile.h_at if (self.profile is not None and self.profile.tapered) else None
224+ self.plate = BFSDisk(ch2, spacing, n_w, symmetric_y=symmetric_y, thickness=thick) # linear modes, modal forces, alpha, omega
207225 pl = self.plate
208226 self.log(f" linear modes: {pl.n_modes}{' (even in y)' if symmetric_y else ''} on {pl.n_cells} cells ({pl.ndof} dof) in {time.perf_counter() - t0:.1f}s, f_max={pl.omega[-1] / (2 * np.pi):.0f} Hz, rigid={pl.n_rigid}")
209227 self.n_w = pl.n_modes
@@ -211,11 +229,15 @@ class VonKarmanDisk:
211229 self.alpha = pl.alpha
212230 self.f = pl.f
213231 V = pl.V
214- # Airy modes: clamped biharmonic eigenfunctions, L2-normalized: K_b psi = zeta^4 M_1 psi on interior dofs
232+ # Airy modes: eigenfunctions of the compatibility operator int (1/(E h)) [Psi_xx v_xx + Psi_yy v_yy
233+ # - nu (Psi_xx v_yy + Psi_yy v_xx) + 2 (1 + nu) Psi_xy v_xy] with clamped conditions, L2-normalized:
234+ # K_b psi = zeta^4 M_1 psi on interior dofs (for uniform h this is the biharmonic operator over E h)
215235 t0 = time.perf_counter()
216- Kb = pl.assemble(pl.element_matrix(np.array([[1.0, 1.0, 0.0], [1.0, 1.0, 0.0], [0.0, 0.0, 0.0]])))
217- M1 = pl.M / pl.m
218- z4, Psi = pl.eigenmodes(Kb, M1, n_airy, sigma=0.0, dofs=pl.interior_dofs)
236+ nu = p["nu"]
237+ Cnu = np.array([[1.0, -nu, 0.0], [-nu, 1.0, 0.0], [0.0, 0.0, 2.0 * (1.0 + nu)]])
238+ gx, gy = pl._gauss_xy
239+ Kb = pl.assemble_coeff(Cnu, 1.0 / (p["E"] * pl.h_at(gx, gy)))
240+ z4, Psi = pl.eigenmodes(Kb, pl.M1, n_airy, sigma=0.0, dofs=pl.interior_dofs)
219241 self.n_airy = len(z4)
220242 self.zeta4 = z4
221243 self.log(f" airy modes: {self.n_airy} in {time.perf_counter() - t0:.1f}s")
@@ -226,7 +248,7 @@ class VonKarmanDisk:
226248 Pxx, Pyy, Pxy = Bxx @ V, Byy @ V, Bxy @ V # (Q, n_w)
227249 Bn = pl.basis_at(xq, yq)
228250 PsiQ = Bn @ Psi # (Q, n_airy)
229- self.coef = -(p["E"] * p["h"]) / (2.0 * self.zeta4)
251+ self.coef = -1.0 / (2.0 * self.zeta4) # E h is inside the compatibility operator
230252 if self.torch is None:
231253 H = np.zeros((self.n_airy, self.n_w, self.n_w))
232254 for pi in range(self.n_airy):
@@ -591,7 +613,7 @@ def main():
591613 entry = {
592614 "solution_id": sid,
593615 "output": f"{out_dir.name}/{sid}.wav",
594- "method": ("Marguerre shallow shell (Airy form) on a Cartesian grid masked to the disk, modal method on shell modes built from Bogner-Fox-Schmit plate modes, Gautschi-type integrator exact in the linear part, retarded Rayleigh-type integral with cubic fractional delays" if model.profile is not None else
616+ "method": (("Marguerre shallow shell (Airy form) on a Cartesian grid masked to the disk, modal method on shell modes built from Bogner-Fox-Schmit plate modes" + (" with variable thickness" if model.profile.tapered else "") + ", Gautschi-type integrator exact in the linear part, retarded Rayleigh-type integral with cubic fractional delays") if model.profile is not None else
595617 "von Karman plate (Airy form) on a Cartesian grid masked to the disk, modal method on Bogner-Fox-Schmit modes, Gautschi-type integrator exact in the linear part, retarded Rayleigh integral with cubic fractional delays"),
596618 "parameters": {
597619 "grid_spacing_m": spacing,