2 * CPU-only checks, no GPU: the surface discretization (area and volume of the
3 * superellipsoid against the exact values) and Green's identity
4 * S_k[du/dn] - D_k[u] = u/2 at wall nodes using the brute-force reference
5 * weights. Run with: node --experimental-strip-types test/test-node.ts
6 */
7import { roundedCube, superellipsoidVolume } from '../src/surface.ts';
8import { green, greenDn } from '../src/problem.ts';
9import { refSelf, refLeaf, momentsToWeights } from '../src/reference.ts';
10import { leggauss } from '../src/gauss.ts';
12let failed = 0;
13const check = (name: string, val: number, tol: number) => {
14 const ok = val < tol;
15 if (!ok) failed++;
16 console.log(`${ok ? 'ok ' : 'FAIL'} ${name}: ${val.toExponential(2)} (tol ${tol.toExponential(0)})`);
17};
19// Gauss-Legendre against known values
20{
21 const { x, w } = leggauss(5);
22 check('leggauss(5) node', Math.abs(x[4] - 0.906179845938664), 1e-14);
23 check('leggauss(5) weight', Math.abs(w[2] - 0.568888888888889), 1e-14);
24}
26const a = 1;
27const p = 8;
28const surf = roundedCube(a, p, 3, 6);
29console.log(`surface: ${surf.npatches} patches, ${surf.npts} nodes`);
31// volume by the divergence theorem: V = (1/3) int x . n dS
32{
33 let vol = 0;
34 let area = 0;
35 const n = surf.npts;
36 for (let i = 0; i < n; i++) {
37 const xn = surf.srcvals[i] * surf.srcvals[9 * n + i] + surf.srcvals[n + i] * surf.srcvals[10 * n + i] + surf.srcvals[2 * n + i] * surf.srcvals[11 * n + i];
38 vol += (xn * surf.wts[i]) / 3;
39 area += surf.wts[i];
40 }
41 const vex = superellipsoidVolume(a, p);
42 check('volume by divergence theorem', Math.abs(vol - vex) / vex, 1e-6);
43 console.log(` area ${area.toFixed(6)}, volume ${vol.toFixed(8)} (exact ${vex.toFixed(8)})`);
44}
46// Green's identity at a few nodes of one patch, brute force over all patches
47{
48 const k = 6;
49 const xe = [1.7, 0.6, -0.5];
50 const n = surf.npts;
51 const npols = surf.npols;
52 const xs = surf.srcvals.subarray(0, 3 * n);
53 const nrm = surf.srcvals.subarray(9 * n, 12 * n);
54 const u = green(xs, xe, k);
55 const un = greenDn(xs, nrm, xe, k);
56 const P = 13;
57 for (const loc of [0, 24, 47]) {
58 const i = P * npols + loc;
59 const xt = [xs[i], xs[n + i], xs[2 * n + i]];
60 const nt = [nrm[i], nrm[n + i], nrm[2 * n + i]];
61 let Sr = 0;
62 let Si = 0;
63 let Dr = 0;
64 let Di = 0;
65 for (let j = 0; j < surf.npatches; j++) {
66 let I;
67 if (j === P) I = refSelf(surf, j, surf.uvNodes[i], surf.uvNodes[n + i], xt, nt, k, 200, 40);
68 else {
69 const d = Math.hypot(xt[0] - surf.centroids[j], xt[1] - surf.centroids[surf.npatches + j], xt[2] - surf.centroids[2 * surf.npatches + j]) / surf.rads[j];
70 I = refLeaf(surf, j, xt, nt, k, d < 3 ? 5 : 2);
71 }
72 const W = momentsToWeights(surf, I);
73 for (let l = 0; l < npols; l++) {
74 const q = j * npols + l;
75 Sr += W[4][2 * l] * un[2 * q] - W[4][2 * l + 1] * un[2 * q + 1];
76 Si += W[4][2 * l] * un[2 * q + 1] + W[4][2 * l + 1] * un[2 * q];
77 Dr += W[5][2 * l] * u[2 * q] - W[5][2 * l + 1] * u[2 * q + 1];
78 Di += W[5][2 * l] * u[2 * q + 1] + W[5][2 * l + 1] * u[2 * q];
79 }
80 }
81 const rr = Sr - Dr - u[2 * i] / 2;
82 const ri = Si - Di - u[2 * i + 1] / 2;
83 check(`Green's identity at node ${loc} of patch ${P}`, Math.hypot(rr, ri) / Math.hypot(u[2 * i], u[2 * i + 1]), 1e-5);
84 }
85}
87if (failed) {
88 console.log(`${failed} check(s) failed`);
89 process.exit(1);
90}
91console.log('all checks passed');