CraftScript meshes and element kernels

export function cylinderSurface(r, h, n = 32) {
  const P = [0, 0, 0, 0, 0, h], I = [], tf = [];
  for (let i = 0; i < n; i++) { const a = 2 * Math.PI * i / n; P.push(r * Math.cos(a), r * Math.sin(a), 0, r * Math.cos(a), r * Math.sin(a), h); }
  const b = (i) => 2 + 2 * (i % n), t = (i) => 3 + 2 * (i % n);
  for (let i = 0; i < n; i++) {
    I.push(0, b(i + 1), b(i)); tf.push(0);
    I.push(1, t(i), t(i + 1)); tf.push(1);
    I.push(b(i), b(i + 1), t(i + 1), b(i), t(i + 1), t(i)); tf.push(2, 2);
  }
  return meshValue('tri', Float64Array.from(P), Int32Array.from(I), { triFace: Int32Array.from(tf) });
}

export function sphereSurface(r, n = 16) {
  const P = [0, 0, -r], I = [], rings = n, seg = 2 * n;
  for (let k = 1; k < rings; k++) { const th = Math.PI * k / rings - Math.PI / 2; for (let i = 0; i < seg; i++) { const a = 2 * Math.PI * i / seg; P.push(r * Math.cos(th) * Math.cos(a), r * Math.cos(th) * Math.sin(a), r * Math.sin(th)); } }
  P.push(0, 0, r);
  const v = (k, i) => 1 + (k - 1) * seg + (i % seg), top = P.length / 3 - 1;
  for (let i = 0; i < seg; i++) I.push(0, v(1, i + 1), v(1, i));
  for (let k = 1; k < rings - 1; k++) for (let i = 0; i < seg; i++) I.push(v(k, i), v(k, i + 1), v(k + 1, i + 1), v(k, i), v(k + 1, i + 1), v(k + 1, i));
  for (let i = 0; i < seg; i++) I.push(v(rings - 1, i), v(rings - 1, i + 1), top);
  return meshValue('tri', Float64Array.from(P), Int32Array.from(I));
}

export function hexGrid(lx, ly, lz, nx, ny, nz, o = [0, 0, 0]) {
  const N = (nx + 1) * (ny + 1) * (nz + 1), X = new Float64Array(3 * N);
  for (let k = 0, n = 0; k <= nz; k++) for (let j = 0; j <= ny; j++) for (let i = 0; i <= nx; i++, n++) { X[3 * n] = o[0] + lx * i / nx; X[3 * n + 1] = o[1] + ly * j / ny; X[3 * n + 2] = o[2] + lz * k / nz; }
  const E = new Int32Array(8 * nx * ny * nz), id = (i, j, k) => i + (nx + 1) * (j + (ny + 1) * k);
  for (let k = 0, e = 0; k < nz; k++) for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++, e++) E.set([id(i, j, k), id(i + 1, j, k), id(i + 1, j + 1, k), id(i, j + 1, k), id(i, j, k + 1), id(i + 1, j, k + 1), id(i + 1, j + 1, k + 1), id(i, j + 1, k + 1)], 8 * e);
  return meshValue('hex8', X, E, { dims: [nx, ny, nz] });
}

// ─────────────────────────────────────────────────────────────── voxelize, tet mesh, surfaces, quality
export function voxelMeshOf(grid) {
  const ids = []; for (let e = 0; e < grid.solid.length; e++) if (grid.solid[e] === 1) ids.push(e);
  const off = VX.cornerOffsets(grid), E = new Int32Array(8 * ids.length);
  ids.forEach((e, m) => { const n0 = VX.elemNode0(grid, e); for (let c = 0; c < 8; c++) E[8 * m + c] = n0 + off[c]; });
  const NX = grid.nx + 1, NY = grid.ny + 1, N = NX * NY * (grid.nz + 1), X = new Float64Array(3 * N);
  for (let n = 0; n < N; n++) { const i = n % NX, j = Math.floor(n / NX) % NY, k = Math.floor(n / (NX * NY)); X[3 * n] = grid.origin[0] + i * grid.h; X[3 * n + 1] = grid.origin[1] + j * grid.h; X[3 * n + 2] = grid.origin[2] + k * grid.h; }
  const g = { h: grid.h, nx: grid.nx, ny: grid.ny, nz: grid.nz, origin: grid.origin.slice(), solid: grid.solid, count: ids.length, oddRows: grid.oddRows || 0 };
  return meshValue('voxel', X, E, { grid: g, elemIds: Int32Array.from(ids) });
}

/** h8/h8i stiffness of a cube of side h → 24×24 (= E·h·Ke of the unit cube). */
export function hexKe(E, nu, h, type = 'h8i') {
  checkEnu(E, nu);
  const el = kHexElement(nu, hexType(type)), s = E * h;
  return mat(24, 24, el.Ke.map((v) => v * s));   // symmetric: row-major = column-major
}

/** Stress operator of a cube: rows (point·6 + component) for the 8 corners then the centre, σ = op · ue → 54×24. */
export function hexStressOp(E, nu, h, type = 'h8i') {
  checkEnu(E, nu);
  const S = kHexElement(nu, hexType(type)).S, s = E / h, out = new Float64Array(54 * 24);
  for (let r = 0; r < 54; r++) for (let j = 0; j < 24; j++) out[r + j * 54] = S[r * 24 + j] * s;
  return mat(54, 24, out);
}

/** Element-mean stress (average of the 4 Gauss points) → ne×6 [xx yy zz xy yz zx]. */
export function t10Stress(m, u, E, nu) {
  m = t10Mesh(m); checkEnu(E, nu);
  const ne = m.elems.length / 10, D = elasticD(nu), el = new Int32Array(10), ue = new Float64Array(30), sig = new Float64Array(6), eps = new Float64Array(6), out = new Float64Array(6 * ne);
  if (u.length !== m.nodes.length) throw err(t('script.mesh.t10_u', { n: m.nodes.length }));
  for (let e = 0; e < ne; e++) {
    for (let a = 0; a < 10; a++) { el[a] = m.elems[10 * e + a]; for (let i = 0; i < 3; i++) ue[3 * a + i] = u[3 * el[a] + i]; }
    for (let p = 0; p < 4; p++) { kGpStress(m.nodes, el, p, ue, E, D, sig, eps); for (let c = 0; c < 6; c++) out[e + c * ne] += sig[c] / 4; }
  }
  return mat(ne, 6, out);
}

/** Element-centre stress of a voxel mesh → ne×6, and von Mises ne×1. */
export function voxelStressOf(m, u, E, nu, type = 'h8i') {
  m = voxelArg(m); checkEnu(E, nu);
  const S = kHexElement(nu, hexType(type)).S, s = E / m.grid.h, ne = m.elems.length / 8, out = new Float64Array(6 * ne), vm = new Float64Array(ne), v = new Float64Array(6);
  if (u.length !== m.nodes.length) throw err(t('script.mesh.voxel_u', { n: m.nodes.length }));
  for (let e = 0; e < ne; e++) {
    for (let r = 0; r < 6; r++) { let q = 0; for (let a = 0; a < 8; a++) { const n = m.elems[8 * e + a]; for (let i = 0; i < 3; i++) q += S[8 * 144 + r * 24 + 3 * a + i] * u[3 * n + i]; } v[r] = q * s; out[e + r * ne] = v[r]; }
    vm[e] = kVonMises(v);
  }
  return [mat(ne, 6, out), mat(ne, 1, vm)];
}

export function voxelSolveRun(m, E, nu, fixed, F, { scale = null, tol = 1e-8, maxit = 0, type = 'h8i', x0 = null } = {}, check = () => {}) {
  m = voxelArg(m); checkEnu(E, nu);
  const t0 = Date.now(), model = kMakeModel(m.grid, { E, nu, type: hexType(type) });
  if (scale && scale.length !== model.elems.length) throw err(t('script.mesh.voxel_scale', { n: model.elems.length }));
  const cons = consFromDofs(fixed, model.nn);
  const free = kRigidCheck(model, cons);
  if (free) throw err(t('script.mesh.voxel_free', { free: free.join(', ') }));
  if (F.length !== 3 * model.nn) throw err(t('script.mesh.voxel_f', { n: 3 * model.nn }));
  const pre = kPreconditioner(model, cons, scale);
  const t1 = Date.now();
  const sol = kVoxelPcg(model, cons, F, { x0, tol, maxIter: maxit, scale, ...pre, cancelled: () => { check(); return false; } });
  return { u: sol.u, iterations: sol.iterations, residual: sol.residual, converged: sol.converged, setupMs: t1 - t0, ms: Date.now() - t0, dofs: 3 * model.active.reduce((s, a) => s + a, 0) };
}

export function lbmFieldOf(Sv) {
  const { st } = lbmOf(Sv), [nx, ny, nz] = st.dims, N = nx * ny * nz, { cells } = st.S.lat;
  const mk = () => new Float64Array(N).fill(NaN), ux = mk(), uy = mk(), uz = mk(), rho = mk();
  for (let k = 0; k < cells.length; k++) { const c = cells[k]; ux[c] = st.S.ux[k]; uy[c] = st.S.uy[k]; uz[c] = st.S.uz[k]; rho[c] = st.S.rho[k]; }
  const shape = nz === 1 ? [nx, ny] : [nx, ny, nz];
  return struct({ ux: tensor(shape, ux), uy: tensor(shape, uy), uz: tensor(shape, uz), rho: tensor(shape, rho), step: st.S.step });
}