Stress, displacement and error results on tetrahedra

/** Principal stresses (descending) of a symmetric stress [xx, yy, zz, xy, yz, zx]. */
export function principal(s) {
  const [a, b, c, d, e, f] = s, m = (a + b + c) / 3, A = a - m, B = b - m, C = c - m;
  const q2 = (A * A + B * B + C * C + 2 * (d * d + e * e + f * f)) / 6, q = Math.sqrt(q2);
  if (q < 1e-14 * (Math.abs(m) + 1e-300)) return [m, m, m];
  const det = (A * (B * C - e * e) - d * (d * C - e * f) + f * (d * e - B * f)) / (q2 * q) / 2, phi = Math.acos(Math.max(-1, Math.min(1, det))) / 3;
  const s1 = m + 2 * q * Math.cos(phi), s3 = m + 2 * q * Math.cos(phi + 2 * Math.PI / 3);
  return [s1, 3 * m - s1 - s3, s3];
}

/** Von Mises equivalent strain √(2/3 e:e) of engineering strains [xx, yy, zz, γxy, γyz, γzx]. */
export const eqStrain = (e) => { const m = (e[0] + e[1] + e[2]) / 3, a = e[0] - m, b = e[1] - m, c = e[2] - m; return Math.sqrt(2 / 3 * (a * a + b * b + c * c + (e[3] * e[3] + e[4] * e[4] + e[5] * e[5]) / 2)); };

const ex = (vals, i) => { // linear field from 4 Gauss values → value at corner i (c_i)
  let S = 0; for (let k = 0; k < 4; k++) S += vals[k]; return (vals[i] - GB * S) / (GA - GB);
};

let emax = 0;
for (let a = 0; a < 10; a++) {
  const [i, j] = a < 4 ? [a, a] : EDGES[a - 4], nd = el[a];
  for (let r = 0; r < 6; r++) { sigN[6 * nd + r] += (cs[6 * i + r] + cs[6 * j + r]) / 2; epsN[6 * nd + r] += (ce[6 * i + r] + ce[6 * j + r]) / 2; }
  cnt[nd]++;
  if (a < 4) emax = Math.max(emax, vonMises(cs.subarray(6 * a, 6 * a + 6)));
}
for (let p = 0; p < 4; p++) emax = Math.max(emax, vonMises(sg[p]));
vmElem[e] = emax; if (emax > elMax) { elMax = emax; elMaxAt = e; }

const err = new Float32Array(m);
let errSum = 0, energy = 0;
const Cinv = (s) => { const [a, b, c, d, e2, f2] = s; return ((a * a + b * b + c * c) - 2 * nu * (a * b + b * c + c * a) + 2 * (1 + nu) * (d * d + e2 * e2 + f2 * f2)) / E; };
for (let e = 0; e < m; e++) {
  let s = 0, w0 = 0;
  for (let a = 0; a < 10; a++) el[a] = model.tets10[10 * e + a];
  for (let p = 0; p < 4; p++) {
    const det = jac(X, el, SH4[p].dN, g), w = GP4[p][3] * det, N = SH4[p].N, sh = elemSig.subarray((4 * e + p) * 6, (4 * e + p) * 6 + 6), d = [0, 0, 0, 0, 0, 0];
    for (let a = 0; a < 10; a++) for (let r = 0; r < 6; r++) d[r] += N[a] * sigN[6 * el[a] + r];
    for (let r = 0; r < 6; r++) d[r] -= sh[r];
    s += w * Cinv(d); w0 += w * Cinv(sh);
  }
  err[e] = Math.sqrt(Math.max(0, s)); errSum += s; energy += w0;
}

const ku = symMul(model, K, u, new Float64Array(model.N));
const reaction = [0, 0, 0], nodeR = new Map();
for (let q = 0; q < cons.nodes.length; q++) {
  const d = 3 * cons.nodes[q], r = [ku[d] - f[d], ku[d + 1] - f[d + 1], ku[d + 2] - f[d + 2]], rr = [0, 0, 0];
  for (const dir of cons.dirs[q]) { const t = r[0] * dir[0] + r[1] * dir[1] + r[2] * dir[2]; for (let a = 0; a < 3; a++) rr[a] += t * dir[a]; }
  for (let a = 0; a < 3; a++) reaction[a] += rr[a];
  nodeR.set(cons.nodes[q], rr);
}

/**
 * 0.4 P5: the error estimate of several fields at once (mode shapes, a buckling study's pre-stress and modes). Each field's
 * Zienkiewicz–Zhu indicators (tetPostprocess) are scaled to its own unit energy, so every field counts alike whatever its
 * amplitude: e_K² = Σ_f e_K,f² / U_f, Σe² = Σ_f Σe²_f / U_f, U = the number of fields, errorPct = 100·√(Σe² / (U + Σe²)) —
 * the RMS of the fields' own estimates for small errors. → { err, errSq, energy, errorPct, errorParts: [% per field] }
 */
export function combineEstimates(parts) {
  const m = parts[0].err.length, err = new Float64Array(m), pct = [];
  let errSq = 0;
  for (const q of parts) {
    const U = Math.max(q.energy, 1e-300);
    for (let e = 0; e < m; e++) err[e] += q.err[e] * q.err[e] / U;
    errSq += q.errSq / U;
    pct.push(q.errorPct);
  }
  for (let e = 0; e < m; e++) err[e] = Math.sqrt(err[e]);
  const energy = parts.length;
  return { err, errSq, energy, errorPct: 100 * Math.sqrt(errSq / (energy + errSq)), errorParts: pct };
}

/** Summary of a meshed model for the UI. */
export function meshSummary(tm, t10, model, rho) {
  let vol = 0; const m = model.tets10.length / 10, g = new Float64Array(30), el = new Int32Array(10);
  for (let e = 0; e < m; e++) { for (let a = 0; a < 10; a++) el[a] = model.tets10[10 * e + a]; for (let p = 0; p < 4; p++) vol += GP4[p][3] * jac(model.X, el, SH4[p].dN, g); }
  return { mesher: 'tet', elements: m, nodes: model.n, dofs: model.N, h: tm.stats.h, volume: vol, mass: vol * rho * 1000, quality: tm.quality, background: tm.stats.background, curvedEdges: t10.stats.curvedEdges, straightKept: t10.stats.straightKept, snapped: tm.stats.snapped, steiner: tm.stats.steiner, exactFaces: tm.stats.exactFaces, tessFaces: tm.stats.tessFaces, tessTypes: tm.stats.tessTypes };
}