1"""A struck sheet of beads: the FPUT chain made two-dimensional, as a toy cymbal or gong.
3N x N unit beads on a square grid move up and down, w_ij; the boundary ring is pinned (w = 0 on it, and the
4bending is simply supported). The sheet may be built as a shallow dome of height Hd, rest shape
5z = 16 Hd x (1 - x) y (1 - y). Energy:
6 links (to the 4 neighbours, including the pinned ring): sum [ T D^2 / 2 + beta (z' D + D^2 / 2)^2 ]
7 bending: Kb / 2 sum (Lap w)^2 (5-point Laplacian)
8with D the height difference across a link and z' its rest slope. The links are the 2D version of the FPUT-beta
9springs (a height difference D stretches a link by about D^2 / 2); bending gives plate-like dispersion; the dome adds
10the quadratic (three-wave) coupling of curved shells. Without a dome the normal modes are products of sines,
11 sin(k pi i / (N + 1)) sin(l pi j / (N + 1)), omega^2 = T lam + Kb lam^2,
12 lam = 4 sin^2(k pi / (2 (N + 1))) + 4 sin^2(l pi / (2 (N + 1))),
13and the mode energies below use this basis (approximate when there is a dome). Damping: friction c0 and
14neighbour dashpots c1 (Laplacian of the velocity), applied explicitly each step. Time stepping: leapfrog.
15"""
17from __future__ import annotations
19from dataclasses import dataclass
21import numpy as np
25class Plate:
26 N: int = 24
27 T: float = 1.0
28 beta: float = 1.0
29 Kb: float = 4.0
30 Hd: float = 0.0
31 c0: float = 0.0
32 c1: float = 0.0
34 def _z(self):
35 g = np.arange(self.N + 2) / (self.N + 1)
36 return 16 * self.Hd * np.outer(g * (1 - g), g * (1 - g))
39 def zh(self):
40 """Rest slopes of the horizontal links, (N, N + 1)."""
41 return np.diff(self._z()[1:-1, :], axis=1)
44 def zv(self):
45 """Rest slopes of the vertical links, (N + 1, N)."""
46 return np.diff(self._z()[:, 1:-1], axis=0)
48 def lam(self):
49 s = 4 * np.sin(np.arange(1, self.N + 1) * np.pi / (2 * (self.N + 1))) ** 2
50 return s[:, None] + s[None, :]
52 def omega(self):
53 lam = self.lam()
54 return np.sqrt(self.T * lam + self.Kb * lam**2)
56 def sine(self):
57 N = self.N
58 j = np.arange(1, N + 1)
59 return np.sqrt(2 / (N + 1)) * np.sin(np.outer(j, j) * np.pi / (N + 1))
62 def lap(u):
63 """5-point Laplacian of (M, N, N) with zero boundary."""
64 p = np.pad(u, ((0, 0), (1, 1), (1, 1)))
65 return p[:, 2:, 1:-1] + p[:, :-2, 1:-1] + p[:, 1:-1, 2:] + p[:, 1:-1, :-2] - 4 * u
67 def force(self, w):
68 p = np.pad(w, ((0, 0), (1, 1), (1, 1)))
69 Dh = np.diff(p[:, 1:-1, :], axis=2) # (M, N, N + 1)
70 Dv = np.diff(p[:, :, 1:-1], axis=1) # (M, N + 1, N)
71 zh, zv, T, b = self.zh, self.zv, self.T, self.beta
72 eh = zh * Dh + 0.5 * Dh * Dh
73 ev = zv * Dv + 0.5 * Dv * Dv
74 fh = T * Dh + 2 * b * eh * (zh + Dh)
75 fv = T * Dv + 2 * b * ev * (zv + Dv)
76 f = (fh[:, :, 1:] - fh[:, :, :-1]) + (fv[:, 1:, :] - fv[:, :-1, :])
77 if self.Kb:
78 f -= self.Kb * self.lap(self.lap(w))
79 return f
81 def energy(self, w, v):
82 p = np.pad(w, ((0, 0), (1, 1), (1, 1)))
83 Dh = np.diff(p[:, 1:-1, :], axis=2)
84 Dv = np.diff(p[:, :, 1:-1], axis=1)
85 eh = self.zh * Dh + 0.5 * Dh * Dh
86 ev = self.zv * Dv + 0.5 * Dv * Dv
87 U = (0.5 * self.T * Dh**2 + self.beta * eh**2).sum((1, 2)) + (0.5 * self.T * Dv**2 + self.beta * ev**2).sum((1, 2))
88 U += 0.5 * self.Kb * (self.lap(w) ** 2).sum((1, 2))
89 return 0.5 * (v * v).sum((1, 2)) + U
91 def mode_energies(self, w, v):
92 S = self.sine()
93 q = np.einsum("ai,mij,jb->mab", S, w, S)
94 p = np.einsum("ai,mij,jb->mab", S, v, S)
95 return 0.5 * (p * p + (self.omega() * q) ** 2)
98def strike(pl: Plate, E, ci, cj, radius, M=1, scale=None):
99 i = np.arange(1, pl.N + 1)
100 prof = np.exp(-0.5 * (((i[:, None] - ci) ** 2 + (i[None, :] - cj) ** 2) / radius**2))
101 prof *= np.sqrt(2 * E / (prof**2).sum())
102 v = np.tile(prof, (M, 1, 1))
103 if scale is not None:
104 v *= np.asarray(scale)[:, None, None]
105 return np.zeros_like(v), v
108def run(pl: Plate, w0, v0, t_end, dt, record_every=100):
109 w, v = np.array(w0, float), np.array(v0, float)
110 n = int(round(t_end / dt))
111 rec = dict(t=[], Ek=[], E=[])
112 f = pl.force(w)
113 for i in range(n):
114 v = v + 0.5 * dt * f
115 w = w + dt * v
116 f = pl.force(w)
117 v = v + 0.5 * dt * f
118 if pl.c0 or pl.c1:
119 v = v - dt * (pl.c0 * v - pl.c1 * pl.lap(v))
120 if (i + 1) % record_every == 0:
121 rec["t"].append((i + 1) * dt)
122 rec["Ek"].append(pl.mode_energies(w, v))
123 rec["E"].append(pl.energy(w, v))
124 out = {k: np.array(x) for k, x in rec.items()}
125 out.update(w=w, v=v)
126 return out