Simulation and flow result plots

/**
 * Streamlines through a cell-centred field F {nx, ny, nz, h, origin, mask (1 fluid), vx, vy, vz} from seed points (world mm):
 * RK2 (midpoint) with steps of `ds` cells along the direction of flow, stopping on leaving the fluid, stalling, or after
 * `maxSteps`. → [{ points: Float32Array (xyz world), speed: Float32Array }]
 */
export function streamlines(F, seeds, { ds = 0.4, maxSteps = 3000, both = false } = {}) {
  const out = [], v = [0, 0, 0], v2 = [0, 0, 0], h = F.h, o = F.origin;
  let vmax = 0;
  for (let e = 0; e < F.mask.length; e++) if (F.mask[e] === 1) vmax = Math.max(vmax, Math.hypot(F.vx[e], F.vy[e], F.vz[e]));
  const stall = 1e-4 * (vmax || 1);
  const trace = (p0, dirSign) => {
    const pts = [], sp = [];
    let x = (p0[0] - o[0]) / h, y = (p0[1] - o[1]) / h, z = (p0[2] - o[2]) / h;
    for (let s = 0; s < maxSteps; s++) {
      if (!isFluidAt(F, x, y, z)) break;
      sampleVelocity(F, x, y, z, v);
      const m = Math.hypot(v[0], v[1], v[2]);
      pts.push(o[0] + x * h, o[1] + y * h, o[2] + z * h); sp.push(m);
      if (m < stall) break;
      const a = dirSign * ds / m;
      sampleVelocity(F, x + 0.5 * a * v[0], y + 0.5 * a * v[1], z + 0.5 * a * v[2], v2);
      const m2 = Math.hypot(v2[0], v2[1], v2[2]);
      if (m2 < stall) break;
      const b = dirSign * ds / m2;
      x += b * v2[0]; y += b * v2[1]; z += b * v2[2];
    }
    return { pts, sp };
  };
  for (const p of seeds) {
    const f = trace(p, 1);
    let pts = f.pts, sp = f.sp;
    if (both) { const r = trace(p, -1); const n = r.sp.length; const rp = []; const rs = []; for (let i = n - 1; i >= 1; i--) { rp.push(r.pts[3 * i], r.pts[3 * i + 1], r.pts[3 * i + 2]); rs.push(r.sp[i]); } pts = rp.concat(pts); sp = rs.concat(sp); }
    if (sp.length >= 2) out.push({ points: Float32Array.from(pts), speed: Float32Array.from(sp) });
  }
  return out;
}

/**
 * A slice of a cell field through the grid: axis 0/1/2, world coordinate `at` → { w, hgt, values: Float32Array (w·hgt,
 * NaN where not fluid), u0/v0 (world min corner of the slice in the in-plane axes), du, axes: [a, b], layer }.
 */
export function slice(F, values, axis, at) {
  const dims = [F.nx, F.ny, F.nz], a = axis === 0 ? 1 : 0, b = axis === 2 ? 1 : 2;
  const layer = Math.max(0, Math.min(dims[axis] - 1, Math.floor((at - F.origin[axis]) / F.h)));
  const w = dims[a], hgt = dims[b], out = new Float32Array(w * hgt);
  const st = [1, F.nx, F.nx * F.ny];
  for (let jb = 0; jb < hgt; jb++) for (let ia = 0; ia < w; ia++) {
    const e = layer * st[axis] + ia * st[a] + jb * st[b];
    out[ia + w * jb] = F.mask[e] === 1 ? values[e] : NaN;
  }
  return { w, hgt, values: out, axes: [a, b], layer, origin: [F.origin[a], F.origin[b]], h: F.h, at: F.origin[axis] + (layer + 0.5) * F.h };
}

function fieldOf(show) {
  const r = F.result, Fd = r.fields;
  if (show === 'pressure') return { values: Fd.p, lo: r.pmin, hi: r.pmax, unit: 'Pa', title: t('study.flow.ui.field.pressure') };
  if (show === 'nut' && Fd.nut) { let hi = 0; for (let e = 0; e < Fd.nut.length; e++) if (Fd.mask[e] === 1 && Fd.nut[e] > hi) hi = Fd.nut[e]; return { values: Fd.nut, lo: 0, hi: hi || 1, unit: '× ν', title: t('study.flow.ui.field.nut') }; }
  const n = Fd.mask.length, s = new Float32Array(n);
  for (let e = 0; e < n; e++) s[e] = Fd.mask[e] === 1 ? Math.hypot(Fd.vx[e], Fd.vy[e], Fd.vz[e]) : NaN;
  return { values: s, lo: 0, hi: r.vmax || 1, unit: 'm/s', title: t('study.flow.ui.field.speed') };
}

/**
 * A display mesh fine enough to plot a field on: every edge longer than `lim` is split at its midpoint (shared
 * edges split once, so the mesh stays closed; a triangle with 1, 2 or 3 split edges becomes 2, 3 or 4), repeated
 * until none is longer, or until about maxTris triangles. → { positions: Float32Array, indices: Uint32Array, triFace }
 */
export function refineMesh(mesh, lim, maxTris = 400000) {
  let P = Array.from(mesh.positions), I = Array.from(mesh.indices), TF = mesh.triFace ? Array.from(mesh.triFace) : null;
  const L2 = lim * lim;
  for (let pass = 0; pass < 12; pass++) {
    const mid = new Map(), d2 = (a, b) => (P[3 * a] - P[3 * b]) ** 2 + (P[3 * a + 1] - P[3 * b + 1]) ** 2 + (P[3 * a + 2] - P[3 * b + 2]) ** 2;
    const key = (a, b) => (a < b ? a * 4294967296 + b : b * 4294967296 + a);
    let marked = 0;
    for (let t = 0; t < I.length; t += 3) for (let k = 0; k < 3; k++) { const a = I[t + k], b = I[t + (k + 1) % 3], q = key(a, b); if (!mid.has(q) && d2(a, b) > L2) { mid.set(q, -1); marked++; } }
    if (!marked || I.length / 3 + marked * 2 > maxTris) break;
    for (const q of mid.keys()) { const a = Math.floor(q / 4294967296), b = q % 4294967296; mid.set(q, P.length / 3); P.push((P[3 * a] + P[3 * b]) / 2, (P[3 * a + 1] + P[3 * b + 1]) / 2, (P[3 * a + 2] + P[3 * b + 2]) / 2); }
    const NI = [], NF = [];
    for (let t = 0; t < I.length; t += 3) {
      const v = [I[t], I[t + 1], I[t + 2]], m = [0, 1, 2].map((k) => mid.get(key(v[k], v[(k + 1) % 3])) ?? -1), f = TF ? TF[t / 3] : 0;
      const tri = (a, b, c) => { NI.push(a, b, c); NF.push(f); };
      const n = m.filter((x) => x >= 0).length;
      if (n === 0) tri(v[0], v[1], v[2]);
      else if (n === 3) { tri(v[0], m[0], m[2]); tri(m[0], v[1], m[1]); tri(m[2], m[1], v[2]); tri(m[0], m[1], m[2]); }
      else {
        // rotate so the first split edge is v0–v1
        let r = m.findIndex((x) => x >= 0);
        if (n === 2 && m[(r + 2) % 3] >= 0) r = (r + 2) % 3;   // the two split edges start at v[r]
        const a = v[r], b = v[(r + 1) % 3], c = v[(r + 2) % 3], ab = m[r], bc = m[(r + 1) % 3];
        if (n === 1) { tri(a, ab, c); tri(ab, b, c); }
        else { tri(a, ab, c); tri(ab, b, bc); tri(ab, bc, c); }
      }
    }
    I = NI; TF = NF;
  }
  return { positions: Float32Array.from(P), indices: Uint32Array.from(I), triFace: Int32Array.from(TF || new Array(I.length / 3).fill(0)) };
}