1"""Exciters, contact patches and the schedule of hits.
3Each hit is a rigid body of mass m that moves only along z, released at the hit time a small
4clearance above the surface with downward speed v, and coupled to the shell through P springs
5with the power law F = k delta^alpha, delta being the local penetration. A stick contacts along
6a short line, a mallet on a small disk. Because the surface is curved, the rigid face is the
7best-fit plane through the patch points and the residual heights give a staggered contact. A
8"force" exciter instead applies a prescribed raised-cosine force at one point, for comparison
9with models driven by a given force; it has no contact and no motion of its own. All
10hits of a run are stacked into fixed-shape arrays (padded to the largest patch) so that the time
11loop handles any number of hits with the same compiled code. See docs/model.md, section
12"Excitation".
13"""
14from __future__ import annotations
16import math
17from dataclasses import dataclass
19import jax.numpy as jnp
20import numpy as np
21from scipy.integrate import quad
23from .config import ForceExciter, Hit, MalletExciter, StickExciter
24from .geometry import CymbalProfile
25from .mesh import Mesh
27HOLD_HEIGHT = 1.0e3 # where an exciter waits before its release [m], far above any motion
31class ContactPatch:
32 """One hit's exciter, resolved on the mesh."""
33 hit: Hit
34 x: np.ndarray # (P,) patch point coordinates on the projected domain
35 y: np.ndarray
36 nodes: np.ndarray # (P, 3) containing-element nodes
37 weights: np.ndarray # (P, 3) barycentric weights
38 offsets: np.ndarray # (P,) surface height minus the exciter's best-fit face plane
39 k_total: float # total contact stiffness [N/m^alpha]
40 alpha: float
41 mass: float
42 gap: float
43 velocity: float
44 time: float
45 force_peak: float = 0.0 # prescribed pulse (force exciter only) [N]
46 force_duration: float = 0.0
49 def prescribed(self) -> bool:
50 return self.force_duration > 0
53 def n_points(self) -> int:
54 return len(self.x)
57 def energy(self) -> float:
58 return 0.5 * self.mass * self.velocity ** 2
60 def max_compression(self) -> float:
61 """Compression at full stop against a rigid wall: m v^2 / 2 = k d^(a+1) / (a+1)."""
62 a = self.alpha
63 return ((a + 1) * self.energy / self.k_total) ** (1.0 / (a + 1))
65 def free_contact_time(self) -> float:
66 """Contact duration against a rigid wall, 2 (d_max / v) int_0^1 dx / sqrt(1 - x^(a+1)).
67 For alpha = 1.5 this is the Hertz value 2.94 d_max / v. For a prescribed force, the pulse length."""
68 if self.prescribed:
69 return self.force_duration
70 a = self.alpha
71 c, _ = quad(lambda s: 1.0 / math.sqrt(1.0 - s ** (a + 1)), 0.0, 1.0)
72 return 2.0 * c * self.max_compression() / self.velocity
75def build_patch(hit: Hit, mesh: Mesh, profile: CymbalProfile) -> ContactPatch:
76 ex = hit.exciter
77 xs, ys = hit.r * math.cos(hit.phi), hit.r * math.sin(hit.phi)
78 if isinstance(ex, StickExciter):
79 ang = hit.phi + ex.patch_angle
80 xi = np.linspace(-0.5, 0.5, ex.patch_points) * ex.patch_length
81 px, py = xs + xi * math.cos(ang), ys + xi * math.sin(ang)
82 k_total, alpha = ex.k, ex.alpha
83 elif isinstance(ex, MalletExciter):
84 alpha = ex.alpha
85 # k such that the head stops after ``compression`` on a rigid wall
86 k_total = (alpha + 1) * 0.5 * ex.mass * hit.velocity ** 2 / ex.compression ** (alpha + 1)
87 th = np.arange(6) * np.pi / 3
88 px = np.concatenate([[xs], xs + ex.radius * np.cos(th)])
89 py = np.concatenate([[ys], ys + ex.radius * np.sin(th)])
90 elif isinstance(ex, ForceExciter):
91 nodes, w = mesh.locate(xs, ys)
92 return ContactPatch(hit, np.array([xs]), np.array([ys]), nodes[None], w[None], np.zeros(1), 0.0, 1.0,
93 1.0, 0.0, 0.0, float(hit.time), float(ex.peak), float(ex.duration))
94 else:
95 raise TypeError(f"unknown exciter type {type(ex).__name__}")
96 loc = [mesh.locate(px[k], py[k]) for k in range(len(px))]
97 nodes = np.array([l[0] for l in loc])
98 weights = np.array([l[1] for l in loc])
99 z0_p = profile.z0_at(np.stack([px, py], 1))
100 # rigid exciter face: best-fit plane through the surface heights of the patch points
101 A = np.column_stack([px - xs, py - ys, np.ones(len(px))])
102 offsets = z0_p - A @ np.linalg.lstsq(A, z0_p, rcond=None)[0]
103 return ContactPatch(hit, px, py, nodes, weights, offsets, float(k_total), float(alpha),
104 float(ex.mass), float(ex.gap), float(hit.velocity), float(hit.time))
107class HitSchedule:
108 """All hits of a run as stacked device arrays (K hits, padded to P points each).
110 Padded points get zero weights and a very negative offset so they never make contact.
111 """
113 def __init__(self, patches: list[ContactPatch], dt: float, dtype):
114 self.patches = patches
115 K = len(patches)
116 P = max(p.n_points for p in patches)
117 nodes = np.zeros((K, P, 3), dtype=np.int32)
118 weights = np.zeros((K, P, 3))
119 offsets = np.full((K, P), -1e6)
120 for k, p in enumerate(patches):
121 n = p.n_points
122 nodes[k, :n] = p.nodes
123 weights[k, :n] = p.weights
124 offsets[k, :n] = p.offsets
125 f = lambda a: jnp.asarray(a, dtype=dtype)
126 self.n_hits, self.n_points = K, P
127 self.nodes = jnp.asarray(nodes)
128 self.weights, self.offsets = f(weights), f(offsets)
129 self.k_point = f([p.k_total / p.n_points for p in patches])
130 self.alpha = f([p.alpha for p in patches])
131 self.mass = f([p.mass for p in patches])
132 self.gap = f([p.gap for p in patches])
133 self.velocity = f([p.velocity for p in patches])
134 self.step = jnp.asarray([int(round(p.time / dt)) for p in patches], dtype=jnp.int32)
135 self.is_force = jnp.asarray([p.prescribed for p in patches])
136 self.f_peak = f([p.force_peak for p in patches])
137 self.f_dur = f([p.force_duration if p.prescribed else 1.0 for p in patches])
138 self.first_point = f(np.arange(P)[None, :] == 0)
139 self.dt = dt
140 self.dtype = dtype
142 def initial_state(self):
143 """(face height, velocity) of every exciter before release."""
144 return (jnp.full(self.n_hits, HOLD_HEIGHT, self.dtype), jnp.zeros(self.n_hits, self.dtype))
146 def contact(self, n, q, zs, vs, f_peak=None):
147 """Contact forces at step n.
149 Returns the per-point forces Fk (K, P), the total force per hit F (K,), and the exciter
150 state (zs, vs) after applying the release of any hit whose time is step n. ``f_peak``
151 replaces the prescribed peak forces (K,), so that they can be differentiated (tangent runs).
152 """
153 f_peak = self.f_peak if f_peak is None else f_peak
154 wp = jnp.einsum("kpi,kpi->kp", self.weights, q[self.nodes, 2]) + self.offsets
155 release = n == self.step
156 active = n >= self.step
157 zs = jnp.where(release, wp.max(1) + self.gap, zs)
158 vs = jnp.where(release, -self.velocity, vs)
159 delta = jnp.maximum(wp - zs[:, None], 0.0)
160 Fk = jnp.where(active[:, None], self.k_point[:, None] * jnp.power(delta, self.alpha[:, None]), 0.0)
161 # prescribed pulses (k_point is 0 for these, so the contact term above vanishes)
162 t = (n - self.step).astype(self.dtype) * self.dt
163 on = self.is_force & (t >= 0) & (t <= self.f_dur)
164 pulse = jnp.where(on, 0.5 * f_peak * (1 - jnp.cos(2 * jnp.pi * t / self.f_dur)), 0.0)
165 Fk = Fk + pulse[:, None] * self.first_point
166 return Fk, Fk.sum(1), zs, vs
168 def apply_to_plate(self, fo, Fk):
169 """Scatter -Fk onto the w DOFs of the patch nodes."""
170 return fo.at[self.nodes, 2].add(-Fk[:, :, None] * self.weights)
172 def advance(self, zs, vs, F, dt):
173 """Symplectic Euler for the exciters (they are integrated with the plate's velocity update).
174 A prescribed force does not move its exciter."""
175 vs1 = vs + dt * jnp.where(self.is_force, 0.0, F) / self.mass
176 return zs + dt * vs1, vs1
178 def energy(self, vs):
179 return 0.5 * (self.mass * vs * vs).sum()