2 * Quadrilateral patch surfaces in the fmm3dbie convention (iptype 11), and the
3 * rounded cube (superellipsoid) built from them. A direct port of
4 * fmm3dbie_gpu/surface.py from the CUDA proof of concept.
5 *
6 * Each patch is the image of the reference square [-1, 1]^2 and carries
7 * (norder+1)^2 tensor Gauss-Legendre nodes, with the u index varying fastest,
8 * and an unnormalized Legendre basis P_j(u) P_i(v) in the same order. Arrays
9 * with a shape (r, c) are stored row-major in one Float64Array, index r*c + i.
10 */
11import { leggauss, legendreVander } from './gauss.ts';
13/** Reference nodes (2, npols) with u in the first row, and the tensor weights. */
14export function discNodes(norder: number): { uv: Float64Array; w: Float64Array } {
15 const n1 = norder + 1;
16 const { x, w } = leggauss(n1);
17 const npols = n1 * n1;
18 const uv = new Float64Array(2 * npols);
19 const ww = new Float64Array(npols);
20 for (let i = 0; i < n1; i++) {
21 for (let j = 0; j < n1; j++) {
22 const n = i * n1 + j;
23 uv[n] = x[j];
24 uv[npols + n] = x[i];
25 ww[n] = w[i] * w[j];
26 }
27 }
28 return { uv, w: ww };
29}
31/** Basis values (npols, npts): index n = i*(norder+1) + j is P_j(u) P_i(v). */
32export function basis2d(uv: Float64Array, norder: number): Float64Array {
33 const npts = uv.length / 2;
34 const n1 = norder + 1;
35 const Pu = legendreVander(uv.subarray(0, npts), norder);
36 const Pv = legendreVander(uv.subarray(npts), norder);
37 const B = new Float64Array(n1 * n1 * npts);
38 for (let i = 0; i < n1; i++) {
39 for (let j = 0; j < n1; j++) {
40 const n = i * n1 + j;
41 for (let q = 0; q < npts; q++) B[n * npts + q] = Pv[i * npts + q] * Pu[j * npts + q];
42 }
43 }
44 return B;
45}
47/** U (npols, npols) with coefs = U @ vals for values at the tensor GL nodes. */
48export function valsToCoefsMatrix(norder: number): Float64Array {
49 const n1 = norder + 1;
50 const { x, w } = leggauss(n1);
51 const P = legendreVander(x, norder);
52 const U1 = new Float64Array(n1 * n1);
53 for (let k = 0; k < n1; k++) for (let j = 0; j < n1; j++) U1[k * n1 + j] = ((2 * k + 1) / 2) * P[k * n1 + j] * w[j];
54 const npols = n1 * n1;
55 const U = new Float64Array(npols * npols);
56 for (let i = 0; i < n1; i++) {
57 for (let ii = 0; ii < n1; ii++) {
58 for (let j = 0; j < n1; j++) {
59 for (let jj = 0; jj < n1; jj++) U[(i * n1 + ii) * npols + (j * n1 + jj)] = U1[i * n1 + j] * U1[ii * n1 + jj];
60 }
61 }
62 }
63 return U;
64}
66/** M (npols_o, npols) with values at the order-nover nodes = M @ values at the order-norder nodes. */
67export function oversampleMatrix(norder: number, nover: number): Float64Array {
68 const { uv } = discNodes(nover);
69 const npolsO = (nover + 1) ** 2;
70 const npols = (norder + 1) ** 2;
71 const B = basis2d(uv, norder); // (npols, npols_o)
72 const U = valsToCoefsMatrix(norder);
73 const M = new Float64Array(npolsO * npols);
74 for (let k = 0; k < npolsO; k++) {
75 for (let m = 0; m < npols; m++) {
76 let s = 0;
77 for (let n = 0; n < npols; n++) s += B[n * npolsO + k] * U[n * npols + m];
78 M[k * npols + m] = s;
79 }
80 }
81 return M;
82}
84/** (n, e1, e2) with e1 x e2 = n for the six faces of the cube. */
85export const FACES: [number[], number[], number[]][] = [
86 [[1, 0, 0], [0, 1, 0], [0, 0, 1]],
87 [[-1, 0, 0], [0, 0, 1], [0, 1, 0]],
88 [[0, 1, 0], [0, 0, 1], [1, 0, 0]],
89 [[0, -1, 0], [1, 0, 0], [0, 0, 1]],
90 [[0, 0, 1], [1, 0, 0], [0, 1, 0]],
91 [[0, 0, -1], [0, 1, 0], [1, 0, 0]],
92];
94/** m+1 breakpoints on [-1, 1], blending uniform and sine-clustered spacing. */
95export function gradedBreaks(m: number, grading: number): Float64Array {
96 const br = new Float64Array(m + 1);
97 for (let i = 0; i <= m; i++) {
98 const t = -1 + (2 * i) / m;
99 br[i] = (1 - grading) * t + grading * Math.sin(0.5 * Math.PI * t);
100 }
101 return br;
102}
104/**
105 * xyz, r_u, r_v (9, n) of the face chart par = (n, e1, e2, s0, s1, t0, t1) of the
106 * superellipsoid |x|^p = a^p at reference points uv (2, n).
107 */
108export function superellipsoidChart(par: ArrayLike<number>, a: number, p: number, uv: Float64Array): Float64Array {
109 const n = uv.length / 2;
110 const out = new Float64Array(9 * n);
111 const hs = 0.5 * (par[10] - par[9]);
112 const ht = 0.5 * (par[12] - par[11]);
113 const c = [0, 0, 0];
114 const gN = [0, 0, 0];
115 for (let q = 0; q < n; q++) {
116 const s = par[9] + (uv[q] + 1) * hs;
117 const t = par[11] + (uv[n + q] + 1) * ht;
118 let S = 0;
119 for (let d = 0; d < 3; d++) {
120 c[d] = par[d] + s * par[3 + d] + t * par[6 + d];
121 S += Math.pow(c[d], p);
122 }
123 const N = Math.pow(S, 1 / p);
124 let dNu = 0;
125 let dNv = 0;
126 for (let d = 0; d < 3; d++) {
127 gN[d] = Math.pow(c[d] / N, p - 1);
128 dNu += gN[d] * par[3 + d];
129 dNv += gN[d] * par[6 + d];
130 }
131 for (let d = 0; d < 3; d++) {
132 out[d * n + q] = (a * c[d]) / N;
133 out[(3 + d) * n + q] = a * hs * (par[3 + d] / N - (c[d] * dNu) / (N * N));
134 out[(6 + d) * n + q] = a * ht * (par[6 + d] / N - (c[d] * dNv) / (N * N));
135 }
136 }
137 return out;
138}
140export interface ExactGeometry {
141 params: Float64Array; // (npatches, 13)
142 a: number;
143 p: number;
144}
146export class Surface {
147 readonly npols: number;
148 readonly npts: number;
149 readonly npatches: number;
150 readonly uvs: Float64Array; // (2, npols) reference nodes
151 readonly wref: Float64Array;
152 readonly umat: Float64Array; // (npols, npols)
153 readonly srccoefs: Float64Array; // (9, npts)
154 readonly jac: Float64Array;
155 readonly wts: Float64Array;
156 readonly centroids: Float64Array; // (3, npatches)
157 readonly rads: Float64Array;
158 readonly patchId: Int32Array;
159 readonly uvNodes: Float64Array; // (2, npts)
160 readonly norder: number;
161 readonly srcvals: Float64Array; // (12, npts): xyz, r_u, r_v, unit normal
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 162 /** the analytic chart of every patch, or null for the interpolated geometry */
163 readonly exact: ExactGeometry | null;
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 165 constructor(norder: number, srcvals: Float64Array, exact: ExactGeometry | null = null) {
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 166 this.norder = norder;
167 this.srcvals = srcvals;
168 this.exact = exact;
169 this.npols = (norder + 1) ** 2;
170 this.npts = srcvals.length / 12;
171 if (this.npts % this.npols !== 0) throw new Error('srcvals is not a whole number of patches');
172 this.npatches = this.npts / this.npols;
173 const { uv, w } = discNodes(norder);
174 this.uvs = uv;
175 this.wref = w;
176 this.umat = valsToCoefsMatrix(norder);
177 const { npols, npts, npatches } = this;
178 this.srccoefs = new Float64Array(9 * npts);
179 for (let c = 0; c < 9; c++) {
180 for (let pch = 0; pch < npatches; pch++) {
181 const base = pch * npols;
182 for (let n = 0; n < npols; n++) {
183 let s = 0;
184 for (let j = 0; j < npols; j++) s += this.umat[n * npols + j] * srcvals[c * npts + base + j];
185 this.srccoefs[c * npts + base + n] = s;
186 }
187 }
188 }
189 this.jac = new Float64Array(npts);
190 this.wts = new Float64Array(npts);
191 for (let i = 0; i < npts; i++) {
192 const ru = [srcvals[3 * npts + i], srcvals[4 * npts + i], srcvals[5 * npts + i]];
193 const rv = [srcvals[6 * npts + i], srcvals[7 * npts + i], srcvals[8 * npts + i]];
194 const cx = ru[1] * rv[2] - ru[2] * rv[1];
195 const cy = ru[2] * rv[0] - ru[0] * rv[2];
196 const cz = ru[0] * rv[1] - ru[1] * rv[0];
197 this.jac[i] = Math.hypot(cx, cy, cz);
198 this.wts[i] = this.jac[i] * w[i % npols];
199 }
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 200 // centroid and radius of each patch from a 5x5 sample grid of its chart
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 201 const g5 = new Float64Array(50);
202 for (let i = 0; i < 5; i++) {
203 for (let j = 0; j < 5; j++) {
204 g5[i * 5 + j] = -1 + 0.5 * i;
205 g5[25 + i * 5 + j] = -1 + 0.5 * j;
206 }
207 }
208 this.centroids = new Float64Array(3 * npatches);
209 this.rads = new Float64Array(npatches);
210 for (let pch = 0; pch < npatches; pch++) {
211 const g = this.evalGeometry(pch, g5);
212 const cen = [0, 0, 0];
213 for (let k = 0; k < 25; k++) for (let d = 0; d < 3; d++) cen[d] += g[d * 25 + k] / 25;
214 let r2 = 0;
215 for (let k = 0; k < 25; k++) {
216 let s = 0;
217 for (let d = 0; d < 3; d++) s += (g[d * 25 + k] - cen[d]) ** 2;
218 r2 = Math.max(r2, s);
219 }
220 for (let d = 0; d < 3; d++) this.centroids[d * npatches + pch] = cen[d];
221 this.rads[pch] = Math.sqrt(r2);
222 }
223 this.patchId = new Int32Array(npts);
224 this.uvNodes = new Float64Array(2 * npts);
225 for (let i = 0; i < npts; i++) {
226 this.patchId[i] = Math.floor(i / npols);
227 this.uvNodes[i] = uv[i % npols];
228 this.uvNodes[npts + i] = uv[npols + (i % npols)];
229 }
230 }
233 * xyz, r_u, r_v (9, n) of one patch at reference points uv (2, n): the exact chart
234 * when there is one, otherwise the Legendre interpolant of the node data, as in
235 * fmm3dbie.
236 */
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 237 evalGeometry(patch: number, uv: Float64Array): Float64Array {
238 const e = this.exact;
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 239 if (e) return superellipsoidChart(e.params.subarray(13 * patch, 13 * patch + 13), e.a, e.p, uv);
240 const n = uv.length / 2;
241 const { npols, npts } = this;
242 const B = basis2d(uv, this.norder);
243 const out = new Float64Array(9 * n);
244 const base = patch * npols;
245 for (let c = 0; c < 9; c++) {
246 for (let m = 0; m < npols; m++) {
247 const cm = this.srccoefs[c * npts + base + m];
248 for (let q = 0; q < n; q++) out[c * n + q] += cm * B[m * n + q];
249 }
250 }
251 return out;
254 /** Nodes of the order-nover rule on every patch: xyz (3, N), unit normal (3, N), weights (N). */
255 oversampled(nover: number): { x: Float64Array; n: Float64Array; w: Float64Array } {
256 const { uv, w: wo } = discNodes(nover);
257 const npolsO = (nover + 1) ** 2;
258 const N = npolsO * this.npatches;
259 const x = new Float64Array(3 * N);
260 const nrm = new Float64Array(3 * N);
261 const w = new Float64Array(N);
262 for (let pch = 0; pch < this.npatches; pch++) {
263 const g = this.evalGeometry(pch, uv);
264 for (let k = 0; k < npolsO; k++) {
265 const i = pch * npolsO + k;
266 const ru = [g[3 * npolsO + k], g[4 * npolsO + k], g[5 * npolsO + k]];
267 const rv = [g[6 * npolsO + k], g[7 * npolsO + k], g[8 * npolsO + k]];
268 const cx = ru[1] * rv[2] - ru[2] * rv[1];
269 const cy = ru[2] * rv[0] - ru[0] * rv[2];
270 const cz = ru[0] * rv[1] - ru[1] * rv[0];
271 const jac = Math.hypot(cx, cy, cz);
272 x[i] = g[k];
273 x[N + i] = g[npolsO + k];
274 x[2 * N + i] = g[2 * npolsO + k];
275 nrm[i] = cx / jac;
276 nrm[N + i] = cy / jac;
277 nrm[2 * N + i] = cz / jac;
278 w[i] = jac * wo[k];
279 }
280 }
281 return { x, n: nrm, w };
282 }
283}
285/**
286 * The superellipsoid |x|^p + |y|^p + |z|^p = a^p with m x m patches per cube
287 * face, graded toward the face edges where the curvature is concentrated.
288 */
289export function roundedCube(a: number, p: number, m: number, norder: number, grading = 0.5): Surface {
290 const { uv } = discNodes(norder);
291 const npols = (norder + 1) ** 2;
292 const br = gradedBreaks(m, grading);
293 const npatches = 6 * m * m;
294 const npts = npatches * npols;
295 const params = new Float64Array(13 * npatches);
296 const srcvals = new Float64Array(12 * npts);
297 let pch = 0;
298 for (const [n, e1, e2] of FACES) {
299 for (let i = 0; i < m; i++) {
300 for (let j = 0; j < m; j++) {
301 const par = [...n, ...e1, ...e2, br[i], br[i + 1], br[j], br[j + 1]];
302 params.set(par, 13 * pch);
303 const g = superellipsoidChart(par, a, p, uv);
304 for (let k = 0; k < npols; k++) {
305 const idx = pch * npols + k;
306 for (let c = 0; c < 9; c++) srcvals[c * npts + idx] = g[c * npols + k];
307 const ru = [g[3 * npols + k], g[4 * npols + k], g[5 * npols + k]];
308 const rv = [g[6 * npols + k], g[7 * npols + k], g[8 * npols + k]];
309 const cx = ru[1] * rv[2] - ru[2] * rv[1];
310 const cy = ru[2] * rv[0] - ru[0] * rv[2];
311 const cz = ru[0] * rv[1] - ru[1] * rv[0];
312 const len = Math.hypot(cx, cy, cz);
313 srcvals[9 * npts + idx] = cx / len;
314 srcvals[10 * npts + idx] = cy / len;
315 srcvals[11 * npts + idx] = cz / len;
316 }
317 pch++;
318 }
319 }
320 }
321 return new Surface(norder, srcvals, { params, a, p });
322}
324/** Exact volume of |x|^p + |y|^p + |z|^p <= a^p. */
325export function superellipsoidVolume(a: number, p: number): number {
326 return (8 * a ** 3 * gamma(1 + 1 / p) ** 3) / gamma(1 + 3 / p);
327}
329/** Lanczos approximation of the Gamma function (used only for the volume check). */
330function gamma(z: number): number {
331 const g = 7;
332 const c = [
333 0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059,
334 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7,
335 ];
336 if (z < 0.5) return Math.PI / (Math.sin(Math.PI * z) * gamma(1 - z));
337 z -= 1;
338 let x = c[0];
339 for (let i = 1; i < g + 2; i++) x += c[i] / (z + i);
340 const t = z + g + 0.5;
341 return Math.sqrt(2 * Math.PI) * Math.pow(t, z + 0.5) * Math.exp(-t) * x;
342}