Resources Numerical accuracy Calculation policy
Gathers the results of a contact study between parts and of a 2-D magnet study: maximum stress, displacement and flux density, contact force, and the per-vertex and per-element fields.
The solvers this file calls (nonlinear.js, contact.js, magnet.js) work in 64-bit, and the reported maxima are taken from their 64-bit arrays. The fields handed back to the page are then copied into Float32Array: the per-vertex stress, plastic strain and displacement of a contact study, and the flux density and potential of a magnet study. The probe readout and the field exports read those 32-bit copies.
Open item. Return the magnet study's B and A and the per-vertex result arrays as Float64Array (or keep a 64-bit copy for the probe and the exports) so that only the copy sent to the graphics card is 32-bit.
10 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.
Tests tools/partcraft-test/test-asm-contact.mjs, tools/partcraft-test/test-magnet-study.mjs
Converts a material's density from kg/m³ to tonnes/mm³, the unit the solver uses.
ρ (t/mm³) = ρ (kg/m³) × 1e-12
The factor 1e-12 is the correct power of ten but has no exact binary form, so the product is rounded once in the last place. A yield of zero or less is replaced by Infinity (never yields).
const materialOf = (m) => NL.bilinearMaterial({ E: m.E, nu: m.nu, yield: m.yield > 0 ? m.yield : Infinity, Et: m.Et || 0, density: (m.density || 0) * T_PER_KGM3 });Gives each node the mean of the values of the elements that share it.
value at a node = Σ element values ÷ number of elements at the node
Sums and counts are Float64Array; the terms are added in element order. A node with no element keeps zero, so the division is guarded.
/** Nodal averages of per-element values (each element's value goes to its ten nodes). */
function nodalFromElements(model, perElem) {
const sum = new Float64Array(model.n), cnt = new Float64Array(model.n), m = model.tets10.length / 10;
for (let e = 0; e < m; e++) for (let a = 0; a < 10; a++) { const n = model.tets10[10 * e + a]; sum[n] += perElem[e]; cnt[n]++; }
for (let i = 0; i < model.n; i++) if (cnt[i]) sum[i] /= cnt[i];
return sum;
}The maxima the study reports, and where they occur.
|u| = √(u_x² + u_y² + u_z²); maximum over the elements (stress, plastic strain) and over the nodes (displacement)
The maxima are searched in the 64-bit element and node arrays, before anything is copied to 32-bit. The stress maximum is the largest element value, not the nodal average that is plotted. An element's position is the sum of its ten nodes' coordinates, each divided by 10.
const dispN = new Float64Array(model.n);
for (let i = 0; i < model.n; i++) dispN[i] = Math.hypot(r.u[3 * i], r.u[3 * i + 1], r.u[3 * i + 2]);
const vertex = { vm: at(pm, vmN), plastic: at(pm, apN), disp: at(pm, dispN), u: at(pm, r.u, 3) };
ctx.transfer && ctx.transfer(...Object.values(vertex).map((a) => a.buffer));
const pos = (e) => { const c = [0, 0, 0]; for (let a = 0; a < 10; a++) for (let k = 0; k < 3; k++) c[k] += model.X[3 * model.tets10[10 * e + a] + k] / 10; return c; };
const ev = argmax(ex.vm), ep = argmax(ex.plasticStrain), nd = argmax(dispN);
let plasticElements = 0; for (let e = 0; e < ex.plasticStrain.length; e++) if (ex.plasticStrain[e] > 0) plasticElements++;
return {
kind: 'nonlinear', mesher: 'tet', converged: r.converged, failedAt: r.failedAt ?? null, lambda: r.lambda, summary,
steps: r.steps.map((s) => ({ lambda: s.lambda, iterations: s.iterations, plastic: s.plastic, residual: s.residual })),
maxVm: { value: ex.vm[ev], at: pos(ev) }, maxPlastic: { value: ex.plasticStrain[ep], at: pos(ep) }, plasticElements,
maxDisp: { value: dispN[nd], at: [model.X[3 * nd], model.X[3 * nd + 1], model.X[3 * nd + 2]] },Copies the nodal stress, plastic strain, displacement size and displacement vector to the vertices of the plot mesh.
The copy is a Float32Array, so each value is rounded to about 7 significant digits. The same arrays are returned as the result's vertex fields. Code that reads a result's vertex fields sees these 32-bit values: the probe in ws-simulation.js interpolates vertex.vm and vertex.disp, and bands in addons/advanced-sim/lib/results.js takes its lowest and highest value from them.
const at = (pm, values, comps = 1) => { const out = new Float32Array(pm.node.length * comps); pm.node.forEach((n, v) => { for (let c = 0; c < comps; c++) out[comps * v + c] = values[comps * n + c]; }); return out; };Moves each mesh node of a placed part to the assembly frame with the instance's 4 × 4 matrix.
x′ = M · (x, y, z, 1)
Each coordinate is three products and a translation in 64-bit (applyMat in vec.js), written back into the mesh's node array. The matrix is used as given; the code comment says rigid placements only, and nothing here checks that.
const place = jobs.map((p) => (Array.isArray(p.place) && p.place.length === 16 && p.place.every(Number.isFinite) ? p.place : null));
place.forEach((m, i) => { if (!m) return; const nd = meshes[i].t10.nodes; for (let k = 0; k < nd.length; k += 3) { const q = applyMat(m, [nd[k], nd[k + 1], nd[k + 2]]); nd[k] = q[0]; nd[k + 1] = q[1]; nd[k + 2] = q[2]; } });Chooses the stiffness of the contact springs for each pair of touching faces when none is given.
k_n = autoPenalty(min(E_slave, E_master), h_min, factor); factor = 100 (frictionless) or 1000 (bonded)
h_min is the smallest element size over all the parts, or 1 when none is finite. The penalty method lets the faces overlap slightly; the overlap is reported as maxPenetration. The formula itself is autoPenalty in contact.js.
const hMin = Math.min(...meshes.map((m) => m.tm.stats.h || Infinity), ...jobs.map((p) => p.h || Infinity));
const pairs = (job.pairs || []).map((pr) => ({ type: pr.type === 'bonded' ? 'bonded' : 'frictionless', slave: trisOf(pr.slave), master: trisOf(pr.master),
kn: pr.kn || CT.autoPenalty(Math.min(mats[pr.slave.part].E, mats[pr.master.part].E), Number.isFinite(hMin) ? hMin : 1, pr.factor > 0 ? pr.factor : pr.type === 'bonded' ? 1e3 : 100), ...(pr.kt ? { kt: pr.kt } : {}) })); // factors 1e3 / 100: the penalty's own give is ~1e-3 / 1e-2 of the element size's share (contact.js header); gap-closing loads also need small steps, belowAdds the force at every contact point and finds the highest contact pressure.
F = Σ forceᵢ; p_max = max pressureᵢ
A plain 64-bit sum over the contact report in report order, with no compensation.
let normalForce = 0, maxPressure = 0; for (const q of rep) { normalForce += q.force; if (q.pressure > maxPressure) maxPressure = q.pressure; }The largest flux density of a magnet study and where it is, the force on each requested body, and the area of each region.
B_max = max over the elements of |B|; region area = Σ element areas
The maximum is searched in the solver's own array before the 32-bit copies are made. Region areas are a plain sum of the solver's element areas in element order. The forces come from maxwell in magnet.js.
const G = r.geom, e = argmax(r.Bmag);
const forces = (job.contours || []).map((c) => {
const inside = c.box ? (x, y) => x > c.box[0] && x < c.box[2] && y > c.box[1] && y < c.box[3] : (x, y) => (c.radius !== undefined ? Math.hypot(x - c.at[0], y - c.at[1]) < c.radius : false);
return { id: c.id, ...MG.maxwell(mesh, r, inside) };
});
const area = {}; for (let i = 0; i < mesh.region.length; i++) area[mesh.region[i]] = (area[mesh.region[i]] || 0) + G.area[i];
const out = { kind: 'magnet', stats: mesh.stats, iterations: r.iterations, Bmax: { value: r.Bmag[e], at: [G.cx[e], G.cy[e]] }, areas: area, forces,The flux density of every element and the magnetic vector potential of every node, as returned to the page.
B.x, B.y, B.mag and A are converted with Float32Array.from, so each value keeps about 7 significant digits. The Magnetics add-on writes these arrays to its CSV and VTK exports with 8 significant digits and averages B.mag per region (addons/magnetics/lib/results.js).
mesh: { nodes: mesh.nodes, tris: mesh.tris, region: mesh.region }, B: { x: Float32Array.from(r.Bx), y: Float32Array.from(r.By), mag: Float32Array.from(r.Bmag) }, A: Float32Array.from(r.A), energy: r.energy ?? null, ms: Date.now() - t0 };The largest displacement and von Mises stress in a CalculiX results file read back onto the study's mesh.
σ_vm = √(½·((a − b)² + (b − c)² + (c − a)²) + 3·(d² + e² + f²)); |u| = √(u_x² + u_y² + u_z²)
64-bit throughout; the accuracy is that of the numbers printed in the file. The stress maximum is taken over the file's stress records directly. The read is refused when the file's node count differs from the mesh's.
/** The results file (.frd) the person's own CalculiX wrote for that deck, read back on the rebuilt mesh: the displacements (and the von Mises of S when it holds one). */
export function contactRead(job) {
const { model, jobs } = prepareContact(job);
const frd = CX.parseFrd(String(job.frd || ''));
if (frd.nodes.size !== model.n) throw new Error(`The results file has ${frd.nodes.size} nodes and this study's mesh has ${model.n}: the mesh is not the one the deck was written from. Write the deck again from this study and run that.`); // i18n-ignore
const u = CX.frdDisplacements(frd, model);
if (!u) throw new Error('The results file has no displacements (DISP). Ask CalculiX for U in the step\'s output.'); // i18n-ignore
let nd = 0, dmax = -1;
for (let i = 0; i < model.n; i++) { const d = Math.hypot(u[3 * i], u[3 * i + 1], u[3 * i + 2]); if (d > dmax) { dmax = d; nd = i; } }
const stress = [...frd.blocks].reverse().find((b) => b.name === 'STRESS');
let vm = null;
if (stress) { vm = 0; for (const [, v] of stress.values) if (v.length >= 6) { const [a, b, c, d, e, f] = v; vm = Math.max(vm, Math.sqrt(0.5 * ((a - b) ** 2 + (b - c) ** 2 + (c - a) ** 2) + 3 * (d * d + e * e + f * f))); } }
return { kind: 'calculix', calculix: { nodes: frd.nodes.size, stress: !!stress }, parts: jobs.length, maxDisp: { value: dmax, at: [model.X[3 * nd], model.X[3 * nd + 1], model.X[3 * nd + 2]] }, maxVm: vm === null ? null : { value: vm } };
}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.