1/**
2 * A cymbal as a thin closed body built from quadrilateral patches. A port of
3 * fmm3dbie_gpu/cymbal.py from the CUDA proof of concept; lengths in meters.
4 *
5 * fmm3dbie's integral equations need a closed surface, so the cymbal is not a
6 * zero-thickness shell but a thin body around an axisymmetric midsurface
7 * z = F(rho), rho = x^2 + y^2, with thickness T(rho). The midsurface has a raised
8 * bell at the center and a shallow bow outside it; the thickness decreases from
9 * the center to the rim, with zero slope at the rim. The top and bottom faces are
10 * the offset surfaces M +- (T/2) N, with M the midsurface and N its upward unit
11 * normal, and a half-torus of diameter T closes the body at the rim. Every patch
12 * lies entirely on one face or on the rim, so every patch is analytic.
13 *
14 * The faces are meshed with an O-grid in the (x, y) plane: an equiangular central
15 * square of m_c x m_c patches, then rings of patches, graded toward the rim and
16 * refined across the bell. Neighboring rings need not have the same number of
17 * patches around (the mesh is not conforming). The node derivatives dx/du, dx/dv
18 * are computed with forward-mode dual numbers, exactly to rounding. The surface
19 * carries no analytic chart for the solver: the integrals use the Legendre
20 * interpolant of the node data, as in fmm3dbie.
21 */
22import { discNodes, Surface } from './surface.ts';
24// ------------------------------------------------------------ dual numbers
26/** a + b eps with eps^2 = 0: a value and its derivative along one direction. */
27class D {
28 readonly v: number;
29 readonly d: number;
30 constructor(v: number, d: number) {
31 this.v = v;
32 this.d = d;
33 }
34}
35type DN = D | number;
36const lift = (a: DN): D => (typeof a === 'number' ? new D(a, 0) : a);
37const add = (a: DN, b: DN) => {
38 const x = lift(a);
39 const y = lift(b);
40 return new D(x.v + y.v, x.d + y.d);
41};
42const sub = (a: DN, b: DN) => {
43 const x = lift(a);
44 const y = lift(b);
45 return new D(x.v - y.v, x.d - y.d);
46};
47const mul = (a: DN, b: DN) => {
48 const x = lift(a);
49 const y = lift(b);
50 return new D(x.v * y.v, x.d * y.v + x.v * y.d);
51};
52const div = (a: DN, b: DN) => {
53 const x = lift(a);
54 const y = lift(b);
55 return new D(x.v / y.v, (x.d * y.v - x.v * y.d) / (y.v * y.v));
56};
57const dsqrt = (a: D) => {
58 const s = Math.sqrt(a.v);
59 return new D(s, a.d / (2 * s));
60};
61const dexp = (a: D) => {
62 const e = Math.exp(a.v);
63 return new D(e, e * a.d);
64};
65const dcos = (a: D) => new D(Math.cos(a.v), -Math.sin(a.v) * a.d);
66const dsin = (a: D) => new D(Math.sin(a.v), Math.cos(a.v) * a.d);
67const dtan = (a: D) => {
68 const t = Math.tan(a.v);
69 return new D(t, (1 + t * t) * a.d);
70};
72// ------------------------------------------------------------------ shape
74export interface CymbalParams {
75 /** outer radius, including the rim */
76 radius: number;
77 /** thickness at the rim */
78 tRim: number;
79 /** thickness at the center / tRim */
80 taper: number;
81 bellRadius: number;
82 bellHeight: number;
83 bellWidth: number;
84 /** fall of the midsurface from the center to the rim, apart from the bell */
85 bowHeight: number;
86}
88export const CYMBAL_DEFAULTS: CymbalParams = {
89 radius: 0.2,
90 tRim: 0.01,
91 taper: 2.0,
92 bellRadius: 0.06,
93 bellHeight: 0.025,
94 bellWidth: 0.01,
95 bowHeight: 0.015,
96};
98export class CymbalShape {
99 readonly par: CymbalParams;
100 readonly rMid: number;
102 constructor(par: CymbalParams = CYMBAL_DEFAULTS) {
103 this.par = par;
104 this.rMid = par.radius - par.tRim / 2;
105 }
107 private bellArg(rho: DN): D {
108 const { bellRadius: rb, bellWidth: w } = this.par;
109 return div(sub(rho, rb * rb), 2 * rb * w);
110 }
112 /** midsurface height as a function of rho */
113 F(rho: DN): D {
114 const { bowHeight, bellHeight } = this.par;
115 const R2 = this.rMid ** 2;
116 return add(mul(bowHeight, sub(1, div(rho, R2))), div(bellHeight, add(1, dexp(this.bellArg(rho)))));
117 }
119 /** dF/drho */
120 dF(rho: DN): D {
121 const { bowHeight, bellHeight, bellRadius: rb, bellWidth: w } = this.par;
122 const e = dexp(this.bellArg(rho));
123 const ope = add(1, e);
124 return sub(-bowHeight / this.rMid ** 2, div(mul(bellHeight / (2 * rb * w), e), mul(ope, ope)));
125 }
127 /** thickness as a function of rho */
128 T(rho: DN): D {
129 const { tRim, taper } = this.par;
130 const s = sub(1, div(rho, this.rMid ** 2));
131 return mul(tRim, add(1, mul(taper - 1, mul(s, s))));
132 }
134 /** point of the top (side = 1) or bottom (side = -1) face over (x, y) of the midsurface */
135 facePoint(x: DN, y: DN, side: number): [D, D, D] {
136 const rho = add(mul(x, x), mul(y, y));
137 const g = mul(2, this.dF(rho));
138 const q = dsqrt(add(1, mul(mul(g, g), rho)));
139 const s = div(mul(side / 2, this.T(rho)), q);
140 const sg = mul(s, g);
141 return [sub(x, mul(sg, x)), sub(y, mul(sg, y)), add(this.F(rho), s)];
142 }
144 /** point of the rim at azimuth theta; phi from -pi/2 (bottom face) to pi/2 (top face) */
145 rimPoint(theta: DN, phi: DN): [D, D, D] {
146 const R = this.rMid;
147 const rho = R * R;
148 const fp = 2 * R * this.dF(rho).v;
149 const q = Math.sqrt(1 + fp * fp);
150 const h = this.T(rho).v / 2;
151 const a = mul(h / q, dcos(lift(phi)));
152 const b = mul(h / q, dsin(lift(phi)));
153 const rr = sub(add(R, a), mul(b, fp));
154 const th = lift(theta);
155 return [mul(rr, dcos(th)), mul(rr, dsin(th)), add(add(this.F(rho).v, mul(a, fp)), b)];
156 }
158 /** closed outline of the body in the (r, z) half plane, counterclockwise: [r[], z[]] */
159 crossSection(n = 400): [Float64Array, Float64Array] {
160 const r: number[] = [];
161 const z: number[] = [];
162 for (let i = 0; i < n; i++) {
163 const p = this.facePoint((this.rMid * i) / (n - 1), 0, -1);
164 r.push(p[0].v);
165 z.push(p[2].v);
166 }
167 const nr = Math.floor(n / 4);
168 for (let i = 0; i < nr; i++) {
169 const p = this.rimPoint(0, -Math.PI / 2 + (Math.PI * i) / (nr - 1));
170 r.push(p[0].v);
171 z.push(p[2].v);
172 }
173 for (let i = n - 1; i >= 0; i--) {
174 const p = this.facePoint((this.rMid * i) / (n - 1), 0, 1);
175 r.push(p[0].v);
176 z.push(p[2].v);
177 }
178 return [Float64Array.from(r), Float64Array.from(z)];
179 }
181 /** volume by Pappus' theorem on the cross section (a check) */
182 volume(n = 20000): number {
183 const [r, z] = this.crossSection(n);
184 let s = 0;
185 for (let i = 0; i < r.length; i++) {
186 const j = (i + 1) % r.length;
187 s += ((r[i] * r[i] + r[i] * r[j] + r[j] * r[j]) / 6) * (z[j] - z[i]);
188 }
189 return 2 * Math.PI * s;
190 }
191}
193/**
194 * A test for points inside the body, or within margin of its surface. Every
195 * vertical line meets the cross section in one interval, so the outline, offset
196 * outward by margin, splits at its outermost point into a lower and an upper
197 * chain, and a point is inside when z lies between them at its radius.
198 */
199export function insideTest(shape: CymbalShape, margin = 0): (x: ArrayLike<number>) => boolean {
200 const [r0, z0] = shape.crossSection(4000);
201 const n = r0.length;
202 const r = new Float64Array(n);
203 const z = new Float64Array(n);
204 for (let i = 0; i < n; i++) {
205 const a = (i + n - 1) % n;
206 const b = (i + 1) % n;
207 let tr = r0[b] - r0[a];
208 let tz = z0[b] - z0[a];
209 const l = Math.hypot(tr, tz);
210 tr /= l;
211 tz /= l;
212 r[i] = Math.max(0, r0[i] + margin * tz);
213 z[i] = z0[i] - margin * tr;
214 }
215 let im = 0;
216 for (let i = 1; i < n; i++) if (r[i] > r[im]) im = i;
217 const chain = (idx: number[]) => {
218 const cr: number[] = [];
219 const cz: number[] = [];
220 let m = -Infinity;
221 for (const i of idx) {
222 m = Math.max(m, r[i]);
223 cr.push(m);
224 cz.push(z[i]);
225 }
226 return { cr, cz };
227 };
228 const lo = chain(Array.from({ length: im + 1 }, (_, i) => i));
229 const up = chain(Array.from({ length: n - im }, (_, i) => n - 1 - i));
230 const interp = (c: { cr: number[]; cz: number[] }, x: number) => {
231 const { cr, cz } = c;
232 if (x <= cr[0]) return cz[0];
233 let a = 0;
234 let b = cr.length - 1;
235 if (x >= cr[b]) return cz[b];
236 while (b - a > 1) {
237 const m = (a + b) >> 1;
238 if (cr[m] <= x) a = m;
239 else b = m;
240 }
241 const t = cr[b] > cr[a] ? (x - cr[a]) / (cr[b] - cr[a]) : 0;
242 return cz[a] + t * (cz[b] - cz[a]);
243 };
244 const rmax = r[im];
245 return (x) => {
246 const rr = Math.hypot(x[0], x[1]);
247 return rr <= rmax && x[2] > interp(lo, rr) && x[2] < interp(up, rr);
248 };
249}
251// ----------------------------------------------------------------- the mesh
253export interface RingLayout {
254 mC: number;
255 c: number;
256 radii: number[];
257 counts: number[];
258 nrim: number;
259}
261/**
262 * Radii and patch counts of the O-grid for a target patch size h: the central
263 * square has half-width c and mC x mC patches (mC >= 2); ring i lies between
264 * radii[i] and radii[i + 1] (radii[0] = c means the square) and has counts[i]
265 * patches around. The ring around the square reaches r = 2c, so that it is 0.59c
266 * wide at the square's corners.
267 */
268export function ringLayout(shape: CymbalShape, h: number, grade = 1.5, cFrac = 0.12): RingLayout {
269 const R = shape.rMid;
270 const { tRim, bellRadius, bellWidth } = shape.par;
271 const c = cFrac * R;
272 const mC = Math.max(2, Math.ceil((2 * c) / h));
273 const first = 2 * c;
274 const b0 = bellRadius - 2.5 * bellWidth;
275 const b1 = bellRadius + 2.5 * bellWidth;
276 const wmax = (edge: number) => {
277 let w = Math.min(h, tRim + (grade - 1) * (R - edge));
278 if (edge > b0 && edge - w < b1) w = Math.min(w, 2 * bellWidth);
279 return w;
280 };
281 const out = [R];
282 for (;;) {
283 const last = out[out.length - 1];
284 const w = wmax(last);
285 if (last - w < first + 0.5 * w) break;
286 out.push(last - w);
287 }
288 const radii = [c, first, ...out.reverse()];
289 const counts: number[] = [];
290 for (let i = 0; i < radii.length - 1; i++) {
291 const width = radii[i + 1] - radii[i];
292 const arc = Math.min(h, 2.5 * width);
293 counts.push(4 * Math.max(mC, Math.ceil((2 * Math.PI * radii[i + 1]) / (4 * arc))));
294 }
295 counts[0] = 4 * mC;
296 for (let i = 1; i < counts.length; i++) counts[i] = Math.max(counts[i], counts[i - 1]);
297 return { mC, c, radii, counts, nrim: counts[counts.length - 1] };
298}
300type Chart = (u: D, v: D) => [D, D];
301type PatchMap = (u: D, v: D) => [D, D, D];
303/** planar charts (u, v) in [-1, 1]^2 -> (x, y) of the face patches, oriented so that d/du x d/dv points up */
304function faceCharts(L: RingLayout): Chart[] {
305 const charts: Chart[] = [];
306 const q = Math.PI / 4;
307 const { mC, c, radii, counts } = L;
308 for (let i = 0; i < mC; i++) {
309 for (let j = 0; j < mC; j++) {
310 const a0 = -1 + (2 * i) / mC;
311 const a1 = -1 + (2 * (i + 1)) / mC;
312 const b0 = -1 + (2 * j) / mC;
313 const b1 = -1 + (2 * (j + 1)) / mC;
314 charts.push((u, v) => {
315 const s = add(a0, mul(0.5 * (a1 - a0), add(u, 1)));
316 const t = add(b0, mul(0.5 * (b1 - b0), add(v, 1)));
317 return [mul(c, dtan(mul(q, s))), mul(c, dtan(mul(q, t)))];
318 });
319 }
320 }
321 for (let ir = 0; ir < radii.length - 1; ir++) {
322 const r0 = radii[ir];
323 const r1 = radii[ir + 1];
324 const n = counts[ir];
325 const dth = (2 * Math.PI) / n;
326 for (let k = 0; k < n; k++) {
327 const th0 = -q + k * dth;
328 if (ir === 0) {
329 charts.push((u, v) => {
330 const th = add(th0, mul(0.5 * dth, add(v, 1)));
331 const sec = Math.floor((th.v + q) / (2 * q));
332 const rin = div(c, dcos(sub(th, sec * 2 * q)));
333 const a = mul(0.5, add(u, 1));
334 const rad = add(mul(sub(1, a), rin), mul(a, r1));
335 return [mul(rad, dcos(th)), mul(rad, dsin(th))];
336 });
337 } else {
338 charts.push((u, v) => {
339 const th = add(th0, mul(0.5 * dth, add(v, 1)));
340 const rad = add(r0, mul(0.5 * (r1 - r0), add(u, 1)));
341 return [mul(rad, dcos(th)), mul(rad, dsin(th))];
342 });
343 }
344 }
345 }
346 return charts;
347}
349/** all patch maps (u, v) -> xyz with outward orientation: top faces, bottom faces, rim */
350function patchMaps(shape: CymbalShape, h: number, grade: number): { maps: PatchMap[]; layout: RingLayout } {
351 const layout = ringLayout(shape, h, grade);
352 const charts = faceCharts(layout);
353 const maps: PatchMap[] = [];
354 for (const ch of charts) maps.push((u, v) => shape.facePoint(...ch(u, v), 1));
355 for (const ch of charts) maps.push((u, v) => shape.facePoint(...ch(v, u), -1)); // swapping u and v flips the normal
356 const { nrim } = layout;
357 for (let k = 0; k < nrim; k++) {
358 const th0 = -Math.PI / 4 + (2 * Math.PI * k) / nrim;
359 maps.push((u, v) => shape.rimPoint(add(th0, mul(Math.PI / nrim, add(u, 1))), mul(0.5 * Math.PI, v)));
360 }
361 return { maps, layout };
362}
364/** Node data (12, npts) in the fmm3dbie iptype 11 convention, and the layout. */
365export function cymbalSrcvals(shape: CymbalShape, h: number, norder: number, grade = 1.5): { srcvals: Float64Array; layout: RingLayout } {
366 const { uv } = discNodes(norder);
367 const npols = (norder + 1) ** 2;
368 const { maps, layout } = patchMaps(shape, h, grade);
369 const npts = maps.length * npols;
370 const sv = new Float64Array(12 * npts);
371 for (const [j, f] of maps.entries()) {
372 for (let k = 0; k < npols; k++) {
373 const i = j * npols + k;
374 const pu = f(new D(uv[k], 1), new D(uv[npols + k], 0));
375 const pv = f(new D(uv[k], 0), new D(uv[npols + k], 1));
376 const ru = [pu[0].d, pu[1].d, pu[2].d];
377 const rv = [pv[0].d, pv[1].d, pv[2].d];
378 const nx = ru[1] * rv[2] - ru[2] * rv[1];
379 const ny = ru[2] * rv[0] - ru[0] * rv[2];
380 const nz = ru[0] * rv[1] - ru[1] * rv[0];
381 const l = Math.hypot(nx, ny, nz);
382 const vals = [pu[0].v, pu[1].v, pu[2].v, ...ru, ...rv, nx / l, ny / l, nz / l];
383 for (let c = 0; c < 12; c++) sv[c * npts + i] = vals[c];
384 }
385 }
386 return { srcvals: sv, layout };
387}
389/** The cymbal as a Surface with the interpolated geometry, for a target patch size h (m). */
390export function cymbal(shape: CymbalShape, h: number, norder: number, grade = 1.5): Surface & { layout: RingLayout } {
391 const { srcvals, layout } = cymbalSrcvals(shape, h, norder, grade);
392 const s = new Surface(norder, srcvals, null) as Surface & { layout: RingLayout };
393 s.layout = layout;
394 return s;
395}