Stress analysis

How the precision is kept

The solve is 64-bit: preconditioned conjugate gradients, stopped when the residual is under 1e-8 of the load. The stress and displacement at each node are then stored in 32-bit arrays for the plot, and the reported maximum is read back from them, so it carries about 7 digits. The voxel mesh, not the arithmetic, is the larger limit on this result.

Open item. Take the reported maximum stress and displacement from the 64-bit values.

3 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.

for (; it < cap; it++) {
  let rn = 0; for (let i = 0; i < n; i++) rn += r[i] * r[i];
  res = Math.sqrt(rn) / bn;
  if (res < tol) break;
  if (onProgress && Date.now() - t0 > 200) { t0 = Date.now(); onProgress(Math.min(0.99, Math.log(res) / Math.log(tol)), `${label}: iteration ${it}, residual ${res.toExponential(1)}`); }   // i18n-ignore: progress label, never shown (the workers pass on only the fraction)
  if (cancelled && (it & 15) === 0 && cancelled()) throw new Error('Cancelled');
  applyK(model, p, q, scale); proj(q);
  let pq = 0; for (let i = 0; i < n; i++) pq += p[i] * q[i];
  if (!(pq > 0)) break;
  const a = rz / pq;
  for (let i = 0; i < n; i++) { u[i] += a * p[i]; r[i] -= a * q[i]; }
  prec(r, z);
  let rz2 = 0; for (let i = 0; i < n; i++) rz2 += r[i] * z[i];
  const beta = rz2 / rz; rz = rz2;
  for (let i = 0; i < n; i++) p[i] = z[i] + beta * p[i];
}
return { u, iterations: it, residual: res, converged: res < tol };

export const vonMises = (s) => Math.sqrt(0.5 * ((s[0] - s[1]) ** 2 + (s[1] - s[2]) ** 2 + (s[2] - s[0]) ** 2) + 3 * (s[3] ** 2 + s[4] ** 2 + s[5] ** 2));

const vmNode = new Float32Array(nn), disp = new Float32Array(nn);
let vMax = 0, vAt = -1, dMax = 0, dAt = -1;
for (let n = 0; n < nn; n++) {
  if (!active[n]) { vmNode[n] = NaN; disp[n] = NaN; continue; }
  vmNode[n] = vmSum[n] / cnt[n];
  disp[n] = Math.hypot(u[3 * n], u[3 * n + 1], u[3 * n + 2]);
  if (vmNode[n] > vMax) { vMax = vmNode[n]; vAt = n; }
  if (disp[n] > dMax) { dMax = disp[n]; dAt = n; }
}

In detail

Method

Each solid voxel is an 8-node hexahedron, with Wilson-Taylor incompatible modes by default so that bending is right with few elements. The system is solved matrix-free by preconditioned conjugate gradients with a geometric multigrid cycle.

Number types

64-bit arrays for stiffness, loads, displacements and the solve. 32-bit arrays for the stress and displacement stored at each node and element.

Tolerances

The solver stops when the residual is under 1e-8 of the load.

Checked by

The simulation tests in the PartCraft test tools.

Limits

Small displacements, linear elastic material. Curved faces follow the voxel staircase and lose up to about half a voxel. Stress at re-entrant corners and constrained edges depends on the mesh.

Approach: 32-bit step, flagged

One step stores or exchanges 32-bit values that reach a number a person reads. It is listed here with its open item, and the step is either moved to 64-bit or the app will warn where the number is read.

Every policy that takes this approach