2 * The showcase problem: a point source at x0 inside a rounded cube (the
3 * superellipsoid |x|^p + |y|^p + |z|^p = a^p) with impedance walls. With time
4 * dependence exp(-i omega t) the pressure amplitude u satisfies
5 *
6 * (Delta + k^2) u = -delta(x - x0) in the room,
7 * du/dn + i k lambda u = 0 on the walls,
8 *
9 * with n the outward normal and lambda = -beta, where beta > 0 is the wall
10 * admittance. We write u = G + v with G the free-space field of the source and
11 * find v from fmm3dbie's right-preconditioned combined field equation.
12 *
13 * Complex vectors are interleaved Float64Arrays: [re0, im0, re1, im1, ...].
14 */
16export type CVec = Float64Array;
18export const ALPHA = 1.0;
20/** Free-space Green's function exp(ikr)/(4 pi r) at points x (3, n) for a source x0. */
21export function green(x: Float64Array, x0: ArrayLike<number>, k: number): CVec {
22 const n = x.length / 3;
23 const out = new Float64Array(2 * n);
24 for (let i = 0; i < n; i++) {
25 const r = Math.hypot(x[i] - x0[0], x[n + i] - x0[1], x[2 * n + i] - x0[2]);
26 const s = 1 / (4 * Math.PI * r);
27 out[2 * i] = s * Math.cos(k * r);
28 out[2 * i + 1] = s * Math.sin(k * r);
29 }
30 return out;
31}
33/** Normal derivative of the Green's function at x with unit normals nrm (3, n). */
34export function greenDn(x: Float64Array, nrm: Float64Array, x0: ArrayLike<number>, k: number): CVec {
35 const n = x.length / 3;
36 const out = new Float64Array(2 * n);
37 for (let i = 0; i < n; i++) {
38 const dx = x[i] - x0[0];
39 const dy = x[n + i] - x0[1];
40 const dz = x[2 * n + i] - x0[2];
41 const r = Math.hypot(dx, dy, dz);
42 const s = 1 / (4 * Math.PI * r);
43 const gr = s * Math.cos(k * r);
44 const gi = s * Math.sin(k * r);
45 const dn = (dx * nrm[i] + dy * nrm[n + i] + dz * nrm[2 * n + i]) / r;
46 // g * (ik - 1/r) * dn
47 const fr = -1 / r;
48 const fi = k;
49 out[2 * i] = (gr * fr - gi * fi) * dn;
50 out[2 * i + 1] = (gr * fi + gi * fr) * dn;
51 }
52 return out;
53}
55export function inside(x: ArrayLike<number>, a: number, p: number, margin = 0): boolean {
56 return Math.pow(Math.abs(x[0] / a), p) + Math.pow(Math.abs(x[1] / a), p) + Math.pow(Math.abs(x[2] / a), p) < (1 - margin) ** p;
57}
59/** Impedance data for v = u - G: -(dG/dn + i k lambda G) at the wall nodes. */
60export function roomRhs(xs: Float64Array, nrm: Float64Array, x0: ArrayLike<number>, k: number, lam: number): CVec {
61 const g = green(xs, x0, k);
62 const gn = greenDn(xs, nrm, x0, k);
63 const n = g.length / 2;
64 const out = new Float64Array(2 * n);
65 for (let i = 0; i < n; i++) {
66 // i k lam g = (-k lam g_im, k lam g_re)
67 out[2 * i] = -(gn[2 * i] - k * lam * g[2 * i + 1]);
68 out[2 * i + 1] = -(gn[2 * i + 1] + k * lam * g[2 * i]);
69 }
70 return out;
71}
73/** Impedance data of the field of a source outside the room (the self-test). */
74export function selftestRhs(xs: Float64Array, nrm: Float64Array, xe: ArrayLike<number>, k: number, lam: number): CVec {
75 const r = roomRhs(xs, nrm, xe, k, lam);
76 for (let i = 0; i < r.length; i++) r[i] = -r[i];
77 return r;
78}
80/** n deterministic pseudo-random points in the cube of half-width 0.8 a. */
81export function selftestPoints(n: number, a: number): Float64Array {
82 const out = new Float64Array(3 * n);
83 let s = 12345;
84 const rnd = () => {
85 s = (Math.imul(s, 1664525) + 1013904223) >>> 0;
86 return s / 4294967296;
87 };
88 for (let i = 0; i < 3 * n; i++) out[i] = 0.8 * a * (2 * rnd() - 1);
89 return out;
90}
92export interface Slice {
93 name: 'xy' | 'xz' | 'yz';
94 axes: [number, number];
95 fixed: number;
96 pts: Float64Array; // (3, ngrid^2), row-major over the grid with the first axis varying fastest
97 mask: Uint8Array;
98}
100/** Points on the three planes through x0; each grid runs from -a to a. */
101export function slicePoints(ngrid: number, x0: ArrayLike<number>, a: number, p: number): { grid: Float64Array; slices: Slice[] } {
102 const grid = new Float64Array(ngrid);
103 for (let i = 0; i < ngrid; i++) grid[i] = -a + (2 * a * i) / (ngrid - 1);
104 const defs: [Slice['name'], [number, number], number][] = [
105 ['xy', [0, 1], 2],
106 ['xz', [0, 2], 1],
107 ['yz', [1, 2], 0],
108 ];
109 const slices: Slice[] = [];
110 for (const [name, axes, fixed] of defs) {
111 const N = ngrid * ngrid;
112 const pts = new Float64Array(3 * N);
113 const mask = new Uint8Array(N);
114 for (let row = 0; row < ngrid; row++) {
115 for (let col = 0; col < ngrid; col++) {
116 const q = row * ngrid + col;
117 const x = [0, 0, 0];
118 x[axes[0]] = grid[col];
119 x[axes[1]] = grid[row];
120 x[fixed] = x0[fixed];
121 pts[q] = x[0];
122 pts[N + q] = x[1];
123 pts[2 * N + q] = x[2];
124 mask[q] = inside(x, a, p, 1e-3) ? 1 : 0;
125 }
126 }
127 slices.push({ name, axes, fixed, pts, mask });
128 }
129 return { grid, slices };
130}
132/** Power radiated by the source, k/(4 pi) + Im v(x0), and absorbed by the walls, -k Re(lambda) int |u|^2 dS. */
133export function energyBalance(k: number, lam: number, vAtSource: [number, number], wts: Float64Array, uWall: CVec) {
134 const pSrc = k / (4 * Math.PI) + vAtSource[1];
135 let s = 0;
136 for (let i = 0; i < wts.length; i++) s += wts[i] * (uWall[2 * i] ** 2 + uWall[2 * i + 1] ** 2);
137 const pWall = -k * lam * s;
138 return { pSrc, pWall, mismatch: Math.abs(pSrc - pWall) / Math.abs(pSrc) };
139}
141/** Relative 2-norm difference of two complex vectors. */
142export function relError(u: CVec, uex: CVec): number {
143 let num = 0;
144 let den = 0;
145 for (let i = 0; i < u.length; i++) {
146 num += (u[i] - uex[i]) ** 2;
147 den += uex[i] ** 2;
148 }
149 return Math.sqrt(num / den);
150}
152/** Select the columns of x (3, n) where mask is set. */
153export function maskedPoints(x: Float64Array, mask: Uint8Array): { pts: Float64Array; index: Int32Array } {
154 const n = mask.length;
155 let cnt = 0;
156 for (let i = 0; i < n; i++) cnt += mask[i];
157 const pts = new Float64Array(3 * cnt);
158 const index = new Int32Array(cnt);
159 let q = 0;
160 for (let i = 0; i < n; i++) {
161 if (!mask[i]) continue;
162 pts[q] = x[i];
163 pts[cnt + q] = x[n + i];
164 pts[2 * cnt + q] = x[2 * n + i];
165 index[q++] = i;
166 }
167 return { pts, index };
168}