/**
 * CPU-only checks, no GPU: the surface discretization (area and volume of the
 * superellipsoid against the exact values) and Green's identity
 * S_k[du/dn] - D_k[u] = u/2 at wall nodes using the brute-force reference
 * weights. Run with: node --experimental-strip-types test/test-node.ts
 */
import { roundedCube, superellipsoidVolume } from '../src/surface.ts';
import { green, greenDn } from '../src/problem.ts';
import { refSelf, refLeaf, momentsToWeights } from '../src/reference.ts';
import { leggauss } from '../src/gauss.ts';

let failed = 0;
const check = (name: string, val: number, tol: number) => {
  const ok = val < tol;
  if (!ok) failed++;
  console.log(`${ok ? 'ok  ' : 'FAIL'} ${name}: ${val.toExponential(2)} (tol ${tol.toExponential(0)})`);
};

// Gauss-Legendre against known values
{
  const { x, w } = leggauss(5);
  check('leggauss(5) node', Math.abs(x[4] - 0.906179845938664), 1e-14);
  check('leggauss(5) weight', Math.abs(w[2] - 0.568888888888889), 1e-14);
}

const a = 1;
const p = 8;
const surf = roundedCube(a, p, 3, 6);
console.log(`surface: ${surf.npatches} patches, ${surf.npts} nodes`);

// volume by the divergence theorem: V = (1/3) int x . n dS
{
  let vol = 0;
  let area = 0;
  const n = surf.npts;
  for (let i = 0; i < n; i++) {
    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];
    vol += (xn * surf.wts[i]) / 3;
    area += surf.wts[i];
  }
  const vex = superellipsoidVolume(a, p);
  check('volume by divergence theorem', Math.abs(vol - vex) / vex, 1e-6);
  console.log(`     area ${area.toFixed(6)}, volume ${vol.toFixed(8)} (exact ${vex.toFixed(8)})`);
}

// Green's identity at a few nodes of one patch, brute force over all patches
{
  const k = 6;
  const xe = [1.7, 0.6, -0.5];
  const n = surf.npts;
  const npols = surf.npols;
  const xs = surf.srcvals.subarray(0, 3 * n);
  const nrm = surf.srcvals.subarray(9 * n, 12 * n);
  const u = green(xs, xe, k);
  const un = greenDn(xs, nrm, xe, k);
  const P = 13;
  for (const loc of [0, 24, 47]) {
    const i = P * npols + loc;
    const xt = [xs[i], xs[n + i], xs[2 * n + i]];
    const nt = [nrm[i], nrm[n + i], nrm[2 * n + i]];
    let Sr = 0;
    let Si = 0;
    let Dr = 0;
    let Di = 0;
    for (let j = 0; j < surf.npatches; j++) {
      let I;
      if (j === P) I = refSelf(surf, j, surf.uvNodes[i], surf.uvNodes[n + i], xt, nt, k, 200, 40);
      else {
        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];
        I = refLeaf(surf, j, xt, nt, k, d < 3 ? 5 : 2);
      }
      const W = momentsToWeights(surf, I);
      for (let l = 0; l < npols; l++) {
        const q = j * npols + l;
        Sr += W[4][2 * l] * un[2 * q] - W[4][2 * l + 1] * un[2 * q + 1];
        Si += W[4][2 * l] * un[2 * q + 1] + W[4][2 * l + 1] * un[2 * q];
        Dr += W[5][2 * l] * u[2 * q] - W[5][2 * l + 1] * u[2 * q + 1];
        Di += W[5][2 * l] * u[2 * q + 1] + W[5][2 * l + 1] * u[2 * q];
      }
    }
    const rr = Sr - Dr - u[2 * i] / 2;
    const ri = Si - Di - u[2 * i + 1] / 2;
    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);
  }
}

if (failed) {
  console.log(`${failed} check(s) failed`);
  process.exit(1);
}
console.log('all checks passed');