1"""The bead model: beads, elements, strikes and run settings.
3Bead indices are 0-based throughout. Element tables are float matrices whose
4first columns hold bead indices, as in the numbl version of the app:
6 springs [i, j, k, L0] or [i, j, k, L0, k3]
7 bends [i, j, k, kb] or [i, j, k, kb, cb] (j the middle bead)
8 hinges [i, j, k, l, kh] or [i, j, k, l, kh, ch] (j-k the shared edge)
9 dashpots [i, j, c] or [i, j, c, ca]
10"""
12from __future__ import annotations
14import copy
15from dataclasses import dataclass, field
16from typing import List, Optional, Sequence, Union
18import numpy as np
20DEFAULT_T60 = 2.5 # s, default decay time when no damping is given
24class Hit:
25 """A strike on a set of beads at time t.
27 Either an impulse (velocity added to each bead at time t, in m/s) or a
28 raised-cosine force pulse on each bead with peak `force` (N) lasting
29 `duration` (s).
30 """
32 beads: Union[int, Sequence[int], np.ndarray]
33 t: float = 0.0
34 velocity: Optional[Sequence[float]] = None
35 force: Optional[Sequence[float]] = None
36 duration: float = 0.0
38 def to_dict(self) -> dict:
39 return {
40 "beads": [int(b) for b in np.atleast_1d(self.beads)],
41 "t": float(self.t),
42 "velocity": None if self.velocity is None else [float(c) for c in self.velocity],
43 "force": None if self.force is None else [float(c) for c in self.force],
44 "duration": float(self.duration),
45 }
48 def from_dict(d: dict) -> "Hit":
49 return Hit(beads=np.asarray(d["beads"], dtype=np.int64), t=d["t"], velocity=d.get("velocity"),
50 force=d.get("force"), duration=d.get("duration", 0.0))
54class Model:
55 """A network of beads joined by springs and other elements.
57 pos (n, 3) rest positions (m)
58 mass scalar or (n,) masses (kg)
59 pinned bead indices (or a boolean mask) held fixed
60 springs, bends, hinges, dashpots: element tables (see module docstring)
61 friction None (default damping, see below), scalar or (n,): force -c0 v on each bead (N s/m)
62 hits list of Hit
63 duration simulated time (s)
64 sample_rate audio rate of the recorded signals (Hz)
65 dt integrator step (s); 0 = automatic
66 safety fraction of the stability limit used when dt is automatic
67 method 'verlet' or 'yoshida4'
68 frame_rate rate of the stored 3D displacement frames (Hz), for animation
69 signal_dir direction along which per-bead signals are recorded (unit vector)
70 record quantities recorded for every bead at the audio rate:
71 any of 'acceleration', 'velocity', 'displacement'
72 modes compute linear modes and record their energies at each frame
73 listen bead suggested to the viewer as the initial listener ('sum' or an index)
74 name, notes free text shown by the viewer
76 When friction is None and no element carries damping, every mode decays
77 60 dB in DEFAULT_T60 seconds (friction 2 sigma m per bead).
78 """
80 pos: np.ndarray
81 mass: Union[float, np.ndarray] = 1.0
82 pinned: Union[Sequence[int], np.ndarray] = field(default_factory=list)
83 springs: Optional[np.ndarray] = None
84 bends: Optional[np.ndarray] = None
85 hinges: Optional[np.ndarray] = None
86 dashpots: Optional[np.ndarray] = None
87 friction: Union[None, float, np.ndarray] = None
88 hits: List[Hit] = field(default_factory=list)
89 duration: float = 2.0
90 sample_rate: float = 48000.0
91 dt: float = 0.0
92 safety: float = 0.5
93 method: str = "verlet"
94 frame_rate: float = 2000.0
95 signal_dir: Sequence[float] = (0.0, 0.0, 1.0)
96 record: Sequence[str] = ("acceleration", "displacement")
97 modes: bool = False
98 listen: Union[str, int] = "sum"
99 name: str = ""
100 notes: str = ""
101 decay_law: Optional[Sequence[float]] = None # [sigma0, sigma1], set by decay(); informational
104 def n(self) -> int:
105 return int(np.asarray(self.pos).shape[0])
107 def copy(self) -> "Model":
108 return copy.deepcopy(self)
110 def scaled_hits(self, factor: float) -> "Model":
111 """A copy with every strike scaled by factor (velocities and forces)."""
112 m = self.copy()
113 for h in m.hits:
114 if h.velocity is not None:
115 h.velocity = [factor * c for c in h.velocity]
116 if h.force is not None:
117 h.force = [factor * c for c in h.force]
118 return m
121class ModelError(ValueError):
122 pass
125def _table(a, name: str, cols: Sequence[int], nidx: int, n: int) -> np.ndarray:
126 """Validate an element table; return it padded to the largest column count."""
127 full = max(cols)
128 if a is None:
129 return np.zeros((0, full))
130 a = np.asarray(a, dtype=np.float64)
131 if a.size == 0:
132 return np.zeros((0, full))
133 if a.ndim != 2 or a.shape[1] not in cols:
134 raise ModelError(f"{name} must have {' or '.join(str(c) for c in cols)} columns")
135 idx = a[:, :nidx]
136 if np.any(idx != np.round(idx)) or np.any(idx < 0) or np.any(idx >= n):
137 bad = np.argwhere((idx != np.round(idx)) | (idx < 0) | (idx >= n))[0][0]
138 raise ModelError(f"{name} row {bad}: bead index not in 0..{n - 1}")
139 if a.shape[1] < full:
140 a = np.hstack([a, np.zeros((a.shape[0], full - a.shape[1]))])
141 return a
144def _indices(v, name: str, n: int) -> np.ndarray:
145 if v is None:
146 return np.zeros(0, dtype=np.int64)
147 a = np.asarray(v)
148 if a.dtype == bool:
149 if a.shape != (n,):
150 raise ModelError(f"{name}: a boolean mask must have one entry per bead ({n})")
151 return np.flatnonzero(a)
152 a = np.atleast_1d(a).ravel()
153 if a.size == 0:
154 return np.zeros(0, dtype=np.int64)
155 if np.any(a != np.round(a)) or np.any(a < 0) or np.any(a >= n):
156 raise ModelError(f"{name}: bead indices must be integers in 0..{n - 1}")
157 return a.astype(np.int64)
160def _per_bead(v, name: str, n: int) -> np.ndarray:
161 a = np.asarray(v, dtype=np.float64).ravel()
162 if a.size == 1:
163 return np.full(n, a[0])
164 if a.size != n:
165 raise ModelError(f"{name} must be a scalar or have one entry per bead ({n})")
166 return a.copy()
170class Checked:
171 """A validated model in solver form (0-based int arrays, padded tables)."""
173 model: Model
174 n: int
175 pos0: np.ndarray
176 mass: np.ndarray
177 pinned: np.ndarray # uint8 (n,)
178 springs: np.ndarray # (m, 5)
179 bends: np.ndarray # (m, 5)
180 hinges: np.ndarray # (m, 6)
181 dashpots: np.ndarray # (m, 4)
182 friction: np.ndarray # (n,)
183 hits: List[Hit]
184 signal_dir: np.ndarray
185 warnings: List[str]
188def check(m: Model) -> Checked:
189 """Validate a model and convert it to the arrays the solver uses."""
190 warnings: List[str] = []
191 pos = np.asarray(m.pos, dtype=np.float64)
192 if pos.ndim != 2 or pos.shape[1] != 3 or pos.shape[0] < 1:
193 raise ModelError("pos must be an n-by-3 array of bead positions")
194 n = pos.shape[0]
195 mass = _per_bead(m.mass, "mass", n)
196 if np.any(~(mass > 0)):
197 raise ModelError(f"mass must be positive (bead {int(np.argmin(mass > 0))})")
198 pinned = np.zeros(n, dtype=np.uint8)
199 pinned[_indices(m.pinned, "pinned", n)] = 1
201 S = _table(m.springs, "springs", (4, 5), 2, n)
202 if np.any(S[:, 0] == S[:, 1]):
203 raise ModelError(f"springs row {int(np.argmax(S[:, 0] == S[:, 1]))} joins a bead to itself")
204 if np.any(S[:, 2:5] < 0):
205 raise ModelError("springs: k, L0 and k3 must be nonnegative")
206 B = _table(m.bends, "bends", (4, 5), 3, n)
207 for r in range(B.shape[0]):
208 if len({B[r, 0], B[r, 1], B[r, 2]}) != 3:
209 raise ModelError(f"bends row {r} repeats a bead")
210 if np.any(B[:, 3:5] < 0):
211 raise ModelError("bends: kb and cb must be nonnegative")
212 H = _table(m.hinges, "hinges", (5, 6), 4, n)
213 for r in range(H.shape[0]):
214 if len({H[r, 0], H[r, 1], H[r, 2], H[r, 3]}) != 4:
215 raise ModelError(f"hinges row {r} repeats a bead")
216 if np.any(H[:, 4:6] < 0):
217 raise ModelError("hinges: kh and ch must be nonnegative")
218 D = _table(m.dashpots, "dashpots", (3, 4), 2, n)
219 if np.any(D[:, 2:4] < 0):
220 raise ModelError("dashpots: c and ca must be nonnegative")
222 damped_elements = D.shape[0] > 0 or np.any(B[:, 4] > 0) or np.any(H[:, 5] > 0)
223 if m.friction is None and not damped_elements:
224 sigma = np.log(1000.0) / DEFAULT_T60
225 friction = 2 * sigma * mass
226 warnings.append(f"No damping given: using the default, every mode falls 60 dB in {DEFAULT_T60} s. "
227 "Set friction = 0 for none.")
228 else:
229 friction = _per_bead(0.0 if m.friction is None else m.friction, "friction", n)
231 hits = []
232 for q, h in enumerate(m.hits):
233 beads = _indices(h.beads, f"hits[{q}].beads", n)
234 if beads.size == 0:
235 raise ModelError(f"hits[{q}].beads is required")
236 if h.duration > 0 and h.velocity is not None:
237 raise ModelError(f"hits[{q}]: give velocity (impulse) or force with duration (pulse), not both")
238 if h.duration <= 0 and h.force is not None:
239 raise ModelError(f"hits[{q}].force needs a positive duration")
240 for nm in ("velocity", "force"):
241 val = getattr(h, nm)
242 if val is not None and len(val) != 3:
243 raise ModelError(f"hits[{q}].{nm} must have 3 components")
244 hits.append(Hit(beads=beads, t=float(h.t), velocity=h.velocity, force=h.force, duration=float(h.duration)))
245 if not hits:
246 warnings.append("hits is empty: nothing excites the system.")
248 d = np.asarray(m.signal_dir, dtype=np.float64).ravel()
249 if d.size != 3 or np.linalg.norm(d) == 0:
250 raise ModelError("signal_dir must be a nonzero 3-vector")
251 if m.method not in ("verlet", "yoshida4"):
252 raise ModelError("method must be 'verlet' or 'yoshida4'")
253 if not m.duration > 0:
254 raise ModelError("duration must be positive")
255 if not (0 < m.safety <= 1):
256 raise ModelError("safety must be in (0, 1]")
257 for q in m.record:
258 if q not in ("acceleration", "velocity", "displacement"):
259 raise ModelError(f"record: unknown quantity {q!r}")
260 if S.shape[0] == 0 and B.shape[0] == 0 and H.shape[0] == 0:
261 warnings.append("springs, bends and hinges are empty: the beads do not interact.")
262 return Checked(model=m, n=n, pos0=pos, mass=mass, pinned=pinned, springs=S, bends=B, hinges=H,
263 dashpots=D, friction=friction, hits=hits, signal_dir=d / np.linalg.norm(d), warnings=warnings)