2 * Double precision CPU reference for the tests: the six kernels, and brute
3 * force near weights for one (target, patch) pair (a port of
4 * tests/test_quadrature.py from the CUDA proof of concept).
5 *
6 * The six kernels, in the order used by fmm3dbie's helm_rpcomb_imp solver:
7 * 0: S'_k 1: S_kappa 2: S'_kappa 3: D'_k - D'_kappa 4: S_k 5: D_k
8 * with kappa = i|k|, S' the normal derivative at the target and D the normal
9 * derivative at the source.
10 */
11import { leggauss } from './gauss.ts';
12import { basis2d, type Surface } from './surface.ts';
14const OVER4PI = 1 / (4 * Math.PI);
16/** K (6, n) complex, interleaved per kernel: out[m][2q], out[m][2q+1]. */
17export function kernels6(xt: ArrayLike<number>, nt: ArrayLike<number>, y: Float64Array, ns: Float64Array, k: number): Float64Array[] {
18 const n = y.length / 3;
19 const K = Array.from({ length: 6 }, () => new Float64Array(2 * n));
20 for (let q = 0; q < n; q++) {
21 const dx = xt[0] - y[q];
22 const dy = xt[1] - y[n + q];
23 const dz = xt[2] - y[2 * n + q];
24 const r2 = dx * dx + dy * dy + dz * dz;
25 const r = Math.sqrt(r2);
26 const dnt = dx * nt[0] + dy * nt[1] + dz * nt[2];
27 const dns = dx * ns[q] + dy * ns[n + q] + dz * ns[2 * n + q];
28 const nn = ns[q] * nt[0] + ns[n + q] * nt[1] + ns[2 * n + q] * nt[2];
29 const e0 = OVER4PI / r;
30 const Ekr = e0 * Math.cos(k * r);
31 const Eki = e0 * Math.sin(k * r);
32 const Ec = e0 * Math.exp(-k * r);
33 // fk = 1 - i k r, fc = 1 + k r
34 const fkr = 1;
35 const fki = -k * r;
36 const fc = 1 + k * r;
37 // gk = 3 - 3 i k r - k^2 r^2, gc = 3 + 3 k r + k^2 r^2
38 const gkr = 3 - k * k * r2;
39 const gki = -3 * k * r;
40 const gc = 3 + 3 * k * r + k * k * r2;
41 const fkEkr = fkr * Ekr - fki * Eki;
42 const fkEki = fkr * Eki + fki * Ekr;
43 const gkEkr = gkr * Ekr - gki * Eki;
44 const gkEki = gkr * Eki + gki * Ekr;
45 K[4][2 * q] = Ekr;
46 K[4][2 * q + 1] = Eki;
47 K[1][2 * q] = Ec;
48 K[0][2 * q] = (-dnt * fkEkr) / r2;
49 K[0][2 * q + 1] = (-dnt * fkEki) / r2;
50 K[2][2 * q] = (-dnt * fc * Ec) / r2;
51 K[5][2 * q] = (dns * fkEkr) / r2;
52 K[5][2 * q + 1] = (dns * fkEki) / r2;
53 const c1 = nn / r2;
54 const c2 = (dns * dnt) / (r2 * r2);
55 K[3][2 * q] = c1 * fkEkr - c2 * gkEkr - (c1 * fc * Ec - c2 * gc * Ec);
56 K[3][2 * q + 1] = c1 * fkEki - c2 * gkEki;
57 }
58 return K;
59}
61/** Moments sum_nodes K_m P_n |J| w -> (6, npols) complex. */
62export function momentsAt(surf: Surface, patch: number, uv: Float64Array, w: Float64Array, xt: ArrayLike<number>, nt: ArrayLike<number>, k: number): Float64Array[] {
63 const n = w.length;
64 const g = surf.evalGeometry(patch, uv);
65 const y = g.subarray(0, 3 * n);
66 const ns = new Float64Array(3 * n);
67 const jac = new Float64Array(n);
68 for (let q = 0; q < n; q++) {
69 const ru = [g[3 * n + q], g[4 * n + q], g[5 * n + q]];
70 const rv = [g[6 * n + q], g[7 * n + q], g[8 * n + q]];
71 const cx = ru[1] * rv[2] - ru[2] * rv[1];
72 const cy = ru[2] * rv[0] - ru[0] * rv[2];
73 const cz = ru[0] * rv[1] - ru[1] * rv[0];
74 jac[q] = Math.hypot(cx, cy, cz);
75 ns[q] = cx / jac[q];
76 ns[n + q] = cy / jac[q];
77 ns[2 * n + q] = cz / jac[q];
78 }
79 const K = kernels6(xt, nt, y, ns, k);
80 const P = basis2d(uv, surf.norder);
81 const npols = surf.npols;
82 const I = Array.from({ length: 6 }, () => new Float64Array(2 * npols));
83 for (let m = 0; m < 6; m++) {
84 for (let nn = 0; nn < npols; nn++) {
85 let sr = 0;
86 let si = 0;
87 for (let q = 0; q < n; q++) {
88 const f = P[nn * n + q] * jac[q] * w[q];
89 sr += K[m][2 * q] * f;
90 si += K[m][2 * q + 1] * f;
91 }
92 I[m][2 * nn] = sr;
93 I[m][2 * nn + 1] = si;
94 }
95 }
96 return I;
97}
99function addInto(I: Float64Array[], J: Float64Array[]) {
100 for (let m = 0; m < 6; m++) for (let i = 0; i < I[m].length; i++) I[m][i] += J[m][i];
101}
103/** Self-patch moments: polar rules on the triangles (target, edge), split at the foot of the perpendicular, many nodes. */
104export function refSelf(surf: Surface, patch: number, u0: number, v0: number, xt: ArrayLike<number>, nt: ArrayLike<number>, k: number, nx = 300, nr = 60): Float64Array[] {
105 const gr = leggauss(nr);
106 const gx = leggauss(nx);
107 const I = Array.from({ length: 6 }, () => new Float64Array(2 * surf.npols));
108 for (let e = 0; e < 4; e++) {
109 const h = [1 - u0, -1 - u0, 1 - v0, -1 - v0][e];
110 if (Math.abs(h) < 1e-14) continue;
111 const [x1, x2] = e < 2 ? [-1 - v0, 1 - v0] : [-1 - u0, 1 - u0];
112 for (const xend of [x1, x2]) {
113 if (Math.abs(xend) < 1e-14) continue;
114 const xi1 = Math.asinh(Math.abs(xend) / Math.abs(h));
115 const n = nr * nx;
116 const uv = new Float64Array(2 * n);
117 const w = new Float64Array(n);
118 for (let a = 0; a < nr; a++) {
119 const rho = 0.5 * (gr.x[a] + 1);
120 const wr = 0.5 * gr.w[a];
121 for (let b = 0; b < nx; b++) {
122 const xi = 0.5 * xi1 * (gx.x[b] + 1);
123 const x = Math.sign(xend) * Math.abs(h) * Math.sinh(xi);
124 const dx = Math.abs(h) * Math.cosh(xi) * 0.5 * xi1 * gx.w[b];
125 const q = a * nx + b;
126 if (e < 2) {
127 uv[q] = u0 + rho * h;
128 uv[n + q] = v0 + rho * x;
129 } else {
130 uv[q] = u0 + rho * x;
131 uv[n + q] = v0 + rho * h;
132 }
133 w[q] = wr * dx * Math.abs(h) * rho;
134 }
135 }
136 addInto(I, momentsAt(surf, patch, uv, w, xt, nt, k));
137 }
138 }
139 return I;
140}
142/** Near-patch moments with a uniform quadtree of the given depth and an nq x nq rule per leaf. */
143export function refLeaf(surf: Surface, patch: number, xt: ArrayLike<number>, nt: ArrayLike<number>, k: number, depth = 5, nq = 8): Float64Array[] {
144 const g = leggauss(nq);
145 const nl = 2 ** depth;
146 const h = 1 / nl;
147 const n = nl * nl * nq * nq;
148 const uv = new Float64Array(2 * n);
149 const w = new Float64Array(n);
150 let q = 0;
151 for (let i = 0; i < nl; i++) {
152 for (let j = 0; j < nl; j++) {
153 const cu = -1 + h + 2 * h * i;
154 const cv = -1 + h + 2 * h * j;
155 for (let a = 0; a < nq; a++) {
156 for (let b = 0; b < nq; b++) {
157 uv[q] = cu + h * g.x[a];
158 uv[n + q] = cv + h * g.x[b];
159 w[q] = h * h * g.w[a] * g.w[b];
160 q++;
161 }
162 }
163 }
164 }
165 return momentsAt(surf, patch, uv, w, xt, nt, k);
166}
168/** Node weights w_j = sum_n U[n, j] I_n from the moments: (6, npols) complex. */
169export function momentsToWeights(surf: Surface, I: Float64Array[]): Float64Array[] {
170 const npols = surf.npols;
171 const W = Array.from({ length: 6 }, () => new Float64Array(2 * npols));
172 for (let m = 0; m < 6; m++) {
173 for (let j = 0; j < npols; j++) {
174 let sr = 0;
175 let si = 0;
176 for (let n = 0; n < npols; n++) {
177 const u = surf.umat[n * npols + j];
178 sr += u * I[m][2 * n];
179 si += u * I[m][2 * n + 1];
180 }
181 W[m][2 * j] = sr;
182 W[m][2 * j + 1] = si;
183 }
184 }
185 return W;
186}