Geometric multigrid for the voxel stress model

/** Rᵀ K R for one octant: K (24×24 fine element), w[f][C] → adds s·result into out (24×24). */
function galerkin(K, w, s, out) {
  const T = new Float64Array(576);   // T = K R  (24 × 24)
  for (let i = 0; i < 24; i++) for (let C = 0; C < 8; C++) for (let d = 0; d < 3; d++) {
    let t = 0; for (let f = 0; f < 8; f++) { const wf = w[f][C]; if (wf) t += K[i * 24 + 3 * f + d] * wf; }
    T[i * 24 + 3 * C + d] = t;
  }
  for (let C = 0; C < 8; C++) for (let a = 0; a < 3; a++) for (let j = 0; j < 24; j++) {
    let t = 0; for (let f = 0; f < 8; f++) { const wf = w[f][C]; if (wf) t += wf * T[(3 * f + a) * 24 + j]; }
    out[(3 * C + a) * 24 + j] += s * t;
  }
}

function invertBlocks(B, nn, active) {
  for (let n = 0; n < nn; n++) {
    const a = B.subarray(9 * n, 9 * n + 9);
    const det = a[0] * (a[4] * a[8] - a[5] * a[7]) - a[1] * (a[3] * a[8] - a[5] * a[6]) + a[2] * (a[3] * a[7] - a[4] * a[6]);
    if (!active[n] || !(Math.abs(det) > 0)) { a.fill(0); continue; }
    const inv = [a[4] * a[8] - a[5] * a[7], a[2] * a[7] - a[1] * a[8], a[1] * a[5] - a[2] * a[4], a[5] * a[6] - a[3] * a[8], a[0] * a[8] - a[2] * a[6], a[2] * a[3] - a[0] * a[5], a[3] * a[7] - a[4] * a[6], a[1] * a[6] - a[0] * a[7], a[0] * a[4] - a[1] * a[3]];
    for (let k = 0; k < 9; k++) a[k] = inv[k] / det;
  }
  return B;
}

const smooth = (L, b, x, r, zero) => {
  if (zero) r.fill(0); else L.apply(x, r);
  for (let n = 0; n < L.nn; n++) {
    const d = 3 * n, m = 9 * n, D = L.Dinv, e0 = b[d] - r[d], e1 = b[d + 1] - r[d + 1], e2 = b[d + 2] - r[d + 2];
    x[d] += omega * (D[m] * e0 + D[m + 1] * e1 + D[m + 2] * e2); x[d + 1] += omega * (D[m + 3] * e0 + D[m + 4] * e1 + D[m + 5] * e2); x[d + 2] += omega * (D[m + 6] * e0 + D[m + 7] * e1 + D[m + 8] * e2);
  }
};

// fine node (i,j,k) ← coarse nodes; visit(fineNode, coarseNode, weight)
function eachTransfer(F, C, fn) {
  const FX = F.nx + 1, FY = F.ny + 1, CX = C.nx + 1, CY = C.ny + 1;
  for (let k = 0; k <= F.nz; k++) for (let j = 0; j <= F.ny; j++) for (let i = 0; i <= F.nx; i++) {
    const n = i + FX * (j + FY * k);
    if (!F.active[n]) continue;
    const I = i >> 1, J = j >> 1, K = k >> 1, oi = i & 1, oj = j & 1, ok = k & 1;
    for (let c = 0; c <= ok; c++) for (let b = 0; b <= oj; b++) for (let a = 0; a <= oi; a++) {
      const w = (oi ? 0.5 : 1) * (oj ? 0.5 : 1) * (ok ? 0.5 : 1);
      fn(n, (I + a) + CX * ((J + b) + CY * (K + c)), w);
    }
  }
}

/** Banded Cholesky of a level (explicit element matrices, or the implicit level 0 when it is the only one). */
function bandedCholesky(L, model0) {
  // renumber nodes with the longest axis slowest, so the band is the product of the two shortest
  const N = [L.nx + 1, L.ny + 1, L.nz + 1], ax = [0, 1, 2].sort((a, b) => N[a] - N[b]);
  const perm = new Int32Array(L.nn);
  for (let k = 0; k < N[2]; k++) for (let j = 0; j < N[1]; j++) for (let i = 0; i < N[0]; i++) { const c = [i, j, k]; perm[i + N[0] * (j + N[1] * k)] = c[ax[0]] + N[ax[0]] * (c[ax[1]] + N[ax[1]] * c[ax[2]]); }
  const n = 3 * L.nn, bw = 3 * (N[ax[0]] * N[ax[1]] + N[ax[0]] + 1) + 2;
  const A = new Float64Array(n * (bw + 1));   // row i holds columns i−bw … i at A[i·(bw+1) + (j − i + bw)]
  const pd = (d) => 3 * perm[(d / 3) | 0] + (d % 3);
  const add = (i, j, v) => { i = pd(i); j = pd(j); if (j > i) return; A[i * (bw + 1) + (j - i + bw)] += v; };
  if (L.kind === 'explicit') {
    for (let m = 0; m < L.node0.length; m++) for (let a = 0; a < 24; a++) for (let b = 0; b < 24; b++) {
      const da = 3 * (L.node0[m] + L.off[(a / 3) | 0]) + (a % 3), db = 3 * (L.node0[m] + L.off[(b / 3) | 0]) + (b % 3);
      add(da, db, L.K[576 * m + a * 24 + b]);
    }
  } else {
    const { node0, off, Ke, kScale } = model0;
    for (let m = 0; m < node0.length; m++) {
      const s = (L.scale ? L.scale[m] : 1) * kScale;
      for (let a = 0; a < 24; a++) for (let b = 0; b < 24; b++) add(3 * (node0[m] + off[(a / 3) | 0]) + (a % 3), 3 * (node0[m] + off[(b / 3) | 0]) + (b % 3), s * Ke[a * 24 + b]);
    }
    L.consNodes.forEach((nd, k) => { for (let a = 0; a < 3; a++) for (let b = 0; b < 3; b++) add(3 * nd + a, 3 * nd + b, L.pen[9 * k + a * 3 + b]); });
  }
  for (let nd = 0; nd < L.nn; nd++) if (!L.active[nd]) for (let a = 0; a < 3; a++) A[(3 * perm[nd] + a) * (bw + 1) + bw] = 1;
  for (let i = 0; i < n; i++) {
    const ri = i * (bw + 1);
    for (let j = Math.max(0, i - bw); j <= i; j++) {
      const rj = j * (bw + 1);
      let s = A[ri + (j - i + bw)];
      const k0 = Math.max(0, i - bw, j - bw);
      for (let k = k0; k < j; k++) s -= A[ri + (k - i + bw)] * A[rj + (k - j + bw)];
      if (j === i) A[ri + bw] = Math.sqrt(s > 1e-300 ? s : 1e-300);
      else A[ri + (j - i + bw)] = s / A[rj + bw];
    }
  }
  return { A, n, bw, perm, y: new Float64Array(n) };
}

function cholSolve({ A, n, bw, perm, y }, b) {
  for (let nd = 0; nd < perm.length; nd++) for (let a = 0; a < 3; a++) y[3 * perm[nd] + a] = b[3 * nd + a];
  const x = y;
  for (let i = 0; i < n; i++) { let s = x[i]; const ri = i * (bw + 1); for (let k = Math.max(0, i - bw); k < i; k++) s -= A[ri + (k - i + bw)] * x[k]; x[i] = s / A[ri + bw]; }
  for (let i = n - 1; i >= 0; i--) { x[i] /= A[i * (bw + 1) + bw]; const ri = i * (bw + 1); const xi = x[i]; for (let k = Math.max(0, i - bw); k < i; k++) x[k] -= A[ri + (k - i + bw)] * xi; }
  for (let nd = 0; nd < perm.length; nd++) for (let a = 0; a < 3; a++) b[3 * nd + a] = y[3 * perm[nd] + a];
}