1// Contact patch of one hit, resolved on the modal data: patch points (a line for a stick, a
2// centre plus hexagon for a mallet), their containing triangles, the residual heights above the
3// exciter's best-fit face plane, and the mode shape values at the points. Mirrors
4// cymsim/exciter.py build_patch, with the surface height taken from the nodal z0 by linear
5// interpolation.
6import type { Exciter } from "./config/schema.ts";
7import { planeResiduals } from "./model/linalg.ts";
8import { interpolateNodal, locate, modeValuesAt, type ModesData } from "./modes.ts";
10export interface HitPatch {
11 phi: Float32Array; // (P, Nm) point-major mode values
12 offsets: Float32Array; // (P,)
13 kPoint: number;
14 alpha: number;
15 mass: number;
16 gap: number;
17}
19export function contactStiffness(ex: Exciter, velocity: number): number {
20 if (ex.kind === "stick") return ex.k;
21 return ((ex.alpha + 1) * 0.5 * ex.mass * velocity * velocity) / Math.pow(ex.compression, ex.alpha + 1);
22}
24/** Half-extent of the patch [m], used to keep hits inside the mesh. */
25export function patchExtent(ex: Exciter): number {
26 return ex.kind === "stick" ? ex.patch_length / 2 : ex.radius;
27}
29/** ``maxPoints`` limits the number of springs used in real time (the exported configuration keeps
30 * the exciter's own patch): a stick keeps the centre and the ends of its line, a mallet its centre. */
31export function buildHitPatch(data: ModesData, ex: Exciter, r: number, phiAngle: number, velocity: number,
32 maxPoints = 7, phiFactor = 1): HitPatch | null {
33 const xs = r * Math.cos(phiAngle), ys = r * Math.sin(phiAngle);
34 const px: number[] = [], py: number[] = [];
35 if (ex.kind === "stick") {
36 const ang = phiAngle + ex.patch_angle, P = Math.max(1, Math.min(Math.floor(ex.patch_points), maxPoints));
37 for (let i = 0; i < P; i++) {
38 const xi = (P === 1 ? 0 : -0.5 + i / (P - 1)) * ex.patch_length;
39 px.push(xs + xi * Math.cos(ang));
40 py.push(ys + xi * Math.sin(ang));
41 }
42 } else {
43 px.push(xs);
44 py.push(ys);
45 for (let k = 0; k < 6 && maxPoints > 1; k++) {
46 const th = (k * Math.PI) / 3;
47 px.push(xs + ex.radius * Math.cos(th));
48 py.push(ys + ex.radius * Math.sin(th));
49 }
50 }
51 const P = px.length, Nm = data.nModes;
52 const phi = new Float32Array(P * Nm), z0p = new Float64Array(P);
53 for (let p = 0; p < P; p++) {
54 const loc = locate(data, px[p], py[p]);
55 if (!loc) return null;
56 modeValuesAt(data, loc.nodes, loc.weights, phi, p * Nm, phiFactor);
57 z0p[p] = interpolateNodal(data.z0, loc.nodes, loc.weights);
58 }
59 const offsets = Float32Array.from(planeResiduals(px.map((x) => x - xs), py.map((y) => y - ys), z0p));
60 return { phi, offsets, kPoint: contactStiffness(ex, velocity) / P, alpha: ex.alpha, mass: ex.mass, gap: ex.gap };
61}