Finds which part of a design space to keep for the stiffest shape at a given volume, turns the result into a solid design and a surface mesh, and checks its mass, stress and displacement.
How the precision is kept
The optimization itself runs in Float64Array: densities, filter, sensitivities and the multiplier bisection. Its result is returned as a Float32Array of densities, and the solid design (which elements are kept) is chosen from that 32-bit field. The surface mesh is smoothed and its volume measured in 64-bit, then its vertex positions are returned as a Float32Array. The method is iterative with fixed tolerances, so those set the accuracy before the 32-bit steps do.
Open item. Return the optimized densities and the surface vertex positions as Float64Array.
Calculations under this policy
15 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.
topopt.jsPartCraft 0.8.1 32-bit in part15 calculationsReviewed
Tests tools/partcraft-test/test-simulation.mjs
Density filter
64-bit
Averages each element's density with its neighbors within the minimum member size.
x̃ₘ = Σ w xₙ ÷ Σ w, w = r − distance (in elements)
Weights and their row sums are Float64Array; each filtered value is one weighted sum and one division. The radius is at least 1.0001 elements, so every element has at least its face neighbors.
topopt.js
lines 40 to 66
PartCraft 0.8.1
/** Density filter over the model's elements: apply(src) = Σ w·src / Σ w and its transpose. */
export function makeFilter(model, rElems) {
const g = model.grid, ne = model.elems.length, map = new Int32Array(g.nx * g.ny * g.nz).fill(-1);
model.elems.forEach((e, m) => { map[e] = m; });
const st = filterStencil(Math.max(rElems, 1.0001)), Hs = new Float64Array(ne);
const ijk = new Int32Array(ne * 3);
model.elems.forEach((e, m) => { const i = e % g.nx, j = Math.floor(e / g.nx) % g.ny, k = Math.floor(e / (g.nx * g.ny)); ijk[3 * m] = i; ijk[3 * m + 1] = j; ijk[3 * m + 2] = k; });
const each = (m, fn) => {
const i = ijk[3 * m], j = ijk[3 * m + 1], k = ijk[3 * m + 2];
for (const [di, dj, dk, w] of st) {
const a = i + di, b = j + dj, c = k + dk;
if (a < 0 || b < 0 || c < 0 || a >= g.nx || b >= g.ny || c >= g.nz) continue;
const n = map[a + g.nx * (b + g.ny * c)];
if (n >= 0) fn(n, w);
}
};
// neighbour lists once (CSR): the filter runs ~40 times per iteration inside the bisection
const start = new Int32Array(ne + 1), nb = [], wt = [];
for (let m = 0; m < ne; m++) { each(m, (n, w) => { nb.push(n); wt.push(w); Hs[m] += w; }); start[m + 1] = nb.length; }
const N = Int32Array.from(nb), Wt = Float64Array.from(wt);
return {
apply(src, dst = new Float64Array(ne)) { for (let m = 0; m < ne; m++) { let s = 0; for (let q = start[m]; q < start[m + 1]; q++) s += Wt[q] * src[N[q]]; dst[m] = s / Hs[m]; } return dst; },
// dst_n = Σ_m w_mn src_m / Hs_m
transpose(src, dst = new Float64Array(ne)) { dst.fill(0); for (let m = 0; m < ne; m++) { const v = src[m] / Hs[m]; for (let q = start[m]; q < start[m + 1]; q++) dst[N[q]] += Wt[q] * v; } return dst; },
stencil: st.length,
};
}
Element strain energy
64-bit
uᵀ Kₑ u for every element, the quantity the sensitivities are built from.
s = uₑᵀ Kₑ uₑ
A 24 × 24 quadratic form summed row by row in 64-bit for each element.
topopt.js
lines 68 to 79
PartCraft 0.8.1
/** uᵀ Ke u per element (unit E·h). */
function elementEnergy(model, u, out) {
const { node0, off, Ke } = model, ue = new Float64Array(24);
for (let m = 0; m < node0.length; m++) {
const n0 = node0[m];
for (let c = 0; c < 8; c++) { const d = 3 * (n0 + off[c]); ue[3 * c] = u[d]; ue[3 * c + 1] = u[d + 1]; ue[3 * c + 2] = u[d + 2]; }
let s = 0;
for (let i = 0, r = 0; i < 24; i++, r += 24) { let t = 0; for (let j = 0; j < 24; j++) t += Ke[r + j] * ue[j]; s += ue[i] * t; }
out[m] = s;
}
return out;
}
Stiffness scale, compliance and sensitivities
64-bit
The stiffness of each element from its density, the compliance of the design and how the compliance changes with each density.
E = E_min + xᵖ (1 − E_min); c = Σ E k uᵀKₑu; ∂c/∂x = −p xᵖ⁻¹ (1 − E_min) k uᵀKₑu
E_min is 1e-9 relative to the solid stiffness. The displacement behind the energies comes from a conjugate-gradient solve with a relative residual of 1e-6, so the compliance in the history carries that tolerance.
topopt.js
lines 152 to 160
PartCraft 0.8.1
for (let m = 0; m < ne; m++) scale[m] = EMIN + Math.pow(xPhys[m], p) * (1 - EMIN);
let c = 0, cg = 0; dc.fill(0);
const pre = FE.preconditioner(model, cons, scale);
cases.forEach((bc, ci) => {
const sol = FE.pcg(model, bc.cons, bc.f, { x0: us[ci], scale, ...pre, tol: 1e-6, cancelled, onProgress: (fr, msg) => onProgress((it + fr) / maxIter, t('gen.opt.iter', { n: it + 1, msg })) });
us[ci] = sol.u; cg += sol.iterations;
elementEnergy(model, sol.u, energy);
for (let m = 0; m < ne; m++) { c += model.kScale * scale[m] * energy[m]; dc[m] -= p * Math.pow(xPhys[m], p - 1) * (1 - EMIN) * model.kScale * energy[m]; }
});
Threshold projection
64-bit
Pushes filtered densities toward 0 or 1 when the projection option is on.
x = (tanh βη + tanh β(x̃ − η)) ÷ (tanh βη + tanh β(1 − η)), η = 0.5
Computed with Math.tanh in 64-bit; β doubles from 1 to 16. Preserved elements are set to exactly 1.
topopt.js
lines 140 to 141
PartCraft 0.8.1
const H = (xt, out) => { const d = Math.tanh(beta * eta) + Math.tanh(beta * (1 - eta)); for (let m = 0; m < ne; m++) out[m] = passive[m] ? 1 : proj ? (Math.tanh(beta * eta) + Math.tanh(beta * (xt[m] - eta))) / d : xt[m]; return out; };
const dH = (xt, out) => { const d = Math.tanh(beta * eta) + Math.tanh(beta * (1 - eta)); for (let m = 0; m < ne; m++) out[m] = proj ? beta * (1 - Math.tanh(beta * (xt[m] - eta)) ** 2) / d : 1; return out; };
Density update at the target volume
64-bit
Updates every density and finds the multiplier at which the design's volume meets the target.
The sensitivities are divided by their largest value first, so the bisection range 0 to 1e9 fits any units. The bisection stops at a relative width of 1e-5, so the volume of an iterate matches the target to about that tolerance, not to 15 digits.
topopt.js
lines 173 to 186
PartCraft 0.8.1
let l1 = 0, l2 = 1e9;
const target = volfrac * ne;
let xp = xPhys, xt = xTilde;
while ((l2 - l1) / (l1 + l2) > 1e-5) {
const lm = 0.5 * (l1 + l2);
for (let m = 0; m < ne; m++) {
if (passive[m]) { xnew[m] = 1; continue; }
const be = Math.sqrt(Math.max(0, -dcF[m]) / (dvF[m] * lm));
xnew[m] = Math.max(0, Math.max(x[m] - move, Math.min(1, Math.min(x[m] + move, x[m] * be))));
}
xt = tilde(xnew); xp = H(xt, new Float64Array(ne));
let vol = 0; for (let m = 0; m < ne; m++) vol += xp[m];
if (vol > target) l1 = lm; else l2 = lm;
}
Final compliance and returned densities
32-bit in part
The compliance of the final design and the density field handed to the next steps.
c = Σ f·u
The compliance is re-solved with a relative residual of 1e-8 and summed in 64-bit. The densities are returned with Float32Array.from, so every later step that reads them (the solid design, the surface) starts from 32-bit values.
The number of densest elements kept so the design holds its optimized volume, and the density of the last one kept.
n = round(Σx), at least the preserved elements; level = density of the n-th densest
The sum and the ordering use the 32-bit densities, so two elements whose 64-bit densities differed below 32-bit precision are ordered by their index. The count is rounded to a whole element. When a loaded region is cut off, the count is raised by a bisection that stops within max(1, 0.2 % of the elements).
topopt.js
lines 301 to 312
PartCraft 0.8.1
const n0 = Math.min(ne, Math.max(Math.round(total), passive.reduce((s, v) => s + v, 0)));
let n = n0, b = build(n);
if (!b.allHeld) { // grey designs: add the next densest elements until every kept (loaded) region reaches a support
let lo = n0, hi = ne;
const full = build(hi);
if (full.allHeld) {
while (hi - lo > Math.max(1, ne * 0.002)) { const mid = (lo + hi) >> 1; if (build(mid).allHeld) hi = mid; else lo = mid; }
n = hi; b = build(n);
}
}
const { solid, comp, keepComp, bridged } = b;
const level = n > 0 ? x[order[n - 1]] : 1;
Surface level that holds the optimized volume
32-bit in part
Finds the density level whose surface encloses the optimized material volume.
target = Σx × h³; 10 bisections of the level between 0.05 and 0.95
Ten bisections give the level to 0.9 ÷ 1024, about 0.001, as the function's comment says. The target volume is summed in 64-bit from the density field passed in, which is the 32-bit field when it comes from the optimization result.
topopt.js
lines 333 to 342
PartCraft 0.8.1
if (iso === 'volume') { // the level whose surface holds the optimised volume (Σ densities): bisection
let target = 0; for (let m = 0; m < x.length; m++) target += x[m];
target *= model.grid.h ** 3;
let lo = 0.05, hi = 0.95, best = null;
for (let k = 0; k < 10; k++) {
const mid = (lo + hi) / 2, s = isoSurface(model, x, { iso: mid, smooth });
best = s;
if (s.volume > target) lo = mid; else hi = mid;
}
best.target = target;
Surface vertex in a cell
64-bit
Places one surface vertex in each grid cell the surface passes through.
t = (vₐ − iso) ÷ (vₐ − v_c), kept in [0.1, 0.9]; vertex = mean of the edge crossings
The crossing parameter divides by the difference of two densities; equal densities use 0.5 and any value outside 0.1 to 0.9 is clamped, so the vertex is a bounded estimate, not the exact crossing. For a binary field the vertex is the cell's middle.
topopt.js
lines 394 to 402
PartCraft 0.8.1
let t = (va - vc) !== 0 ? (va - iso) / (va - vc) : 0.5;
if (!(t > 0.1 && t < 0.9)) t = Math.min(0.9, Math.max(0.1, Number.isFinite(t) ? t : 0.5));
const pa = P(a), pc = P(c);
s = [s[0] + pa[0] + t * (pc[0] - pa[0]), s[1] + pa[1] + t * (pc[1] - pa[1]), s[2] + pa[2] + t * (pc[2] - pa[2])]; cnt++;
}
if (!cnt) continue;
vid[i + CX * (j + CY * k)] = pos.length / 3;
if (binary) { const c = P(b); pos.push(c[0] + g.h / 2, c[1] + g.h / 2, c[2] + g.h / 2); } // the cell's middle = a grid node: the voxels' own boundary before smoothing
else pos.push(s[0] / cnt, s[1] / cnt, s[2] / cnt);
Smoothed surface mesh and its volume
32-bit in part
Smooths the surface, restores its volume, and returns the mesh with its volume.
Smoothing, the volume offset and the returned volume all use the Float64Array positions. The positions are then returned through Float32Array.from, so the mesh that leaves this function has 32-bit coordinates while its stated volume was measured on the 64-bit ones.
Smooths the surface mesh without shrinking it as much as plain averaging would.
p ← p + f (mean of neighbors − p), f alternating 0.5 and −0.53
Each pass computes all new positions from the old ones into a temporary Float64Array before copying back, so the result does not depend on vertex order.
topopt.js
lines 432 to 446
PartCraft 0.8.1
/** Taubin smoothing (λ 0.5, μ −0.53) over the mesh's edge graph, in place. */
export function taubin(pos, idx, passes) {
const nv = pos.length / 3, nb = Array.from({ length: nv }, () => new Set());
for (let t = 0; t < idx.length; t += 3) for (let c = 0; c < 3; c++) { const a = idx[t + c], b = idx[t + (c + 1) % 3]; nb[a].add(b); nb[b].add(a); }
const lists = nb.map((s) => Int32Array.from(s)), tmp = new Float64Array(pos.length);
for (let p = 0; p < passes * 2; p++) {
const f = p % 2 === 0 ? 0.5 : -0.53;
for (let v = 0; v < nv; v++) {
const l = lists[v]; if (!l.length) { tmp[3 * v] = pos[3 * v]; tmp[3 * v + 1] = pos[3 * v + 1]; tmp[3 * v + 2] = pos[3 * v + 2]; continue; }
let sx = 0, sy = 0, sz = 0; for (const w of l) { sx += pos[3 * w]; sy += pos[3 * w + 1]; sz += pos[3 * w + 2]; }
tmp[3 * v] = pos[3 * v] + f * (sx / l.length - pos[3 * v]); tmp[3 * v + 1] = pos[3 * v + 1] + f * (sy / l.length - pos[3 * v + 1]); tmp[3 * v + 2] = pos[3 * v + 2] + f * (sz / l.length - pos[3 * v + 2]);
}
pos.set(tmp);
}
}
Volume restored after smoothing
64-bit
Moves every vertex along its normal by one distance so the mesh volume is the voxel volume again.
d = (V₀ − V) ÷ area, three steps
Three fixed-point steps on d = ΔV ÷ area; a step that is not above 1e-12 in size ends the loop. The remaining volume difference is whatever three steps leave, it is not tested against a tolerance.
topopt.js
lines 448 to 465
PartCraft 0.8.1
/** Smoothing thins members one or two elements wide; move every vertex along its normal by one distance so the
* volume is v0 again (three Newton steps on d = ΔV / area). */
export function inflate(pos, idx, v0) {
const nv = pos.length / 3, nrm = new Float64Array(pos.length);
for (let step = 0; step < 3; step++) {
nrm.fill(0); let area = 0;
for (let t = 0; t < idx.length; t += 3) {
const a = 3 * idx[t], b = 3 * idx[t + 1], c = 3 * idx[t + 2];
const u = [pos[b] - pos[a], pos[b + 1] - pos[a + 1], pos[b + 2] - pos[a + 2]], w = [pos[c] - pos[a], pos[c + 1] - pos[a + 1], pos[c + 2] - pos[a + 2]];
const n = [u[1] * w[2] - u[2] * w[1], u[2] * w[0] - u[0] * w[2], u[0] * w[1] - u[1] * w[0]];
area += Math.hypot(...n) / 2;
for (const v of [a, b, c]) for (let k = 0; k < 3; k++) nrm[v + k] += n[k];
}
const d = (v0 - meshVolume(pos, idx)) / area;
if (!(Math.abs(d) > 1e-12)) break;
for (let v = 0; v < nv; v++) { const l = Math.hypot(nrm[3 * v], nrm[3 * v + 1], nrm[3 * v + 2]) || 1; for (let k = 0; k < 3; k++) pos[3 * v + k] += d * nrm[3 * v + k] / l; }
}
}
Volume of a closed mesh
64-bit
The volume enclosed by a triangle mesh.
V = Σ a · (b × c) ÷ 6
Signed tetrahedra from the origin, summed in triangle order in 64-bit. For a mesh far from the origin the terms are large with both signs, so digits are lost in the sum.
topopt.js
lines 467 to 474
PartCraft 0.8.1
export function meshVolume(P, I) {
let v = 0;
for (let t = 0; t < I.length; t += 3) {
const a = 3 * I[t], b = 3 * I[t + 1], c = 3 * I[t + 2];
v += P[a] * (P[b + 1] * P[c + 2] - P[b + 2] * P[c + 1]) - P[a + 1] * (P[b] * P[c + 2] - P[b + 2] * P[c]) + P[a + 2] * (P[b] * P[c + 1] - P[b + 1] * P[c]);
}
return v / 6;
}
Check of the solid design
64-bit
The mass, volume, largest stress and displacement and the safety factor of the kept elements, solved again with the real material.
volume = elements × h³; mass in kg = volume × ρ (t/mm³) × 1000; safety = yield ÷ max von Mises
Volume and mass are a count and three multiplications. The stress and displacement maxima are the values that postprocess in fea.js returns, taken as the worst over the load cases; their precision is recorded with that function. The solve uses a relative residual of 1e-8.
topopt.js
lines 487 to 511
PartCraft 0.8.1
/**
* The design as solid elements (solidDesign: the densest up to the optimised volume, islands dropped) solved again with the real material (h8i) and the
* same supports and loads (all cases, worst value kept) → { mass (kg), volume, maxVm, maxDisp, minSafety }.
*/
export function evaluate(job, out, onProgress = () => {}, cancelled = null, design = null) {
const mat = out.material, rho = tonnesPerMm3(mat.density), g = out.grid;
const d = design || solidDesign(out), solid = new Uint8Array(g.solid.length);
out.model.elems.forEach((e, m) => { if (d.solid[m]) solid[e] = 1; });
const g2 = { ...g, solid };
const model = FE.makeModel(g2, { E: mat.E, nu: mat.nu, type: 'h8i' });
const exactArea = job.faces ? job.faces.map((f) => f.area || 0) : null;
let worst = { maxVm: 0, maxDisp: 0, at: null, compliance: 0 };
const cases = job.cases || [job.loads || []];
cases.forEach((loads, ci) => {
const bc = FE.boundaryConditions(model, job.mesh, { constraints: job.constraints, loads }, { surfaceOnly: out.surfaceOnly, exactArea, density: rho });
const sol = FE.pcg(model, bc.cons, bc.f, { tol: 1e-8, cancelled, ...FE.preconditioner(model, bc.cons), onProgress: (f, msg) => onProgress((ci + f) / cases.length, t('gen.opt.checking', { msg })) });
const r = FE.postprocess(model, sol.u, bc.f, bc.cons, { yieldStress: mat.yield });
if (r.maxVm.value > worst.maxVm) { worst.maxVm = r.maxVm.value; worst.at = r.maxVm.at; }
worst.maxDisp = Math.max(worst.maxDisp, r.maxDisp.value);
worst.compliance += r.compliance;
});
const volume = model.elems.length * g.h ** 3;
let grey = 0; for (const v of out.x) if (v > 0.1 && v < 0.9) grey++;
return { volume, mass: volume * rho * 1000, maxVm: worst.maxVm, maxVmAt: worst.at, maxDisp: worst.maxDisp, compliance: worst.compliance, minSafety: worst.maxVm > 0 ? mat.yield / worst.maxVm : Infinity, elements: model.elems.length, grey: grey / out.x.length };
}
Compliance for another material
64-bit
Rescales the compliance of a shape when an outcome uses a different material.
c = c₀ × E₀ ÷ E
One multiplication and one division in 64-bit. It relies on the statement in the file header that the optimal shape does not depend on E for one material.
SIMP with penalty 3 and E_min = 1e-9 × E₀; the file header cites Sigmund's and Andreassen et al.'s codes. Density filter with linear hat weights; optimality-criteria update with move limit 0.2 and bisection on the multiplier; optional threshold projection (the code cites Wang, Lazarov and Sigmund 2011). Surface by surface nets at a density level, Taubin smoothing (0.5 and −0.53), then a volume-restoring offset.
Number types
Float64Array for densities, sensitivities and the surface during smoothing; Float32Array for the returned densities, the cleaned field and the returned vertex positions; Uint8Array and Int32Array for element sets.
Tolerances
Inner solves: relative residual 1e-6 per iteration, 1e-8 for the final compliance and the check. Multiplier bisection until (l₂ − l₁) ÷ (l₁ + l₂) ≤ 1e-5. Stops when no density changes by more than 0.01 or after 60 iterations. Volume-matching level: 10 bisections between 0.05 and 0.95. Volume offset: three steps, stopping when the step is not above 1e-12.
Limits
The file header lists: one material, linear static, stiffness objective only, gray elements remain at member edges, and the output surface is faceted. The check re-solves the kept elements as voxels, so its mass is the voxel volume times density.
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.