Solves a frame that has plates and rigid floor diaphragms, and works out plate results at nodes and the drift of each story.
How the precision is kept
The stiffness is stored as a skyline in a Float64Array and factorized once by Cholesky. Rigid floors are exact constraints, not penalty springs. Displacements, reactions and plate results are 64-bit. The mesh sets the accuracy of plate moments, and the corner values shown are smoothed from the integration points.
Open item. Run the same single step of iterative refinement after the skyline solve that the frame-only solver runs.
Calculations under this policy
12 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.
The ratio the plate element uses, from the material's two moduli.
ν = E ÷ (2G) − 1
For steel the ratio is about 1.3, so taking 1 off enlarges its relative rounding error about four times. The plate is refused unless −0.999 < ν < 0.4999. The warp in the warning is rounded to 0.01 mm; the value used is not.
The constraint is applied exactly through the assembly, with no penalty stiffness, so it does not add a large number to the matrix.
frame/plate-solve.js
lines 161 to 168
CadCraft 0.8.2
const map = (q) => {
if (eq[q] >= 0) return [[eq[q], 1]];
if (!tied(q)) return [];
const s = slave[Math.floor(q / 6)], b = 3 * s.f.k, d = q % 6;
if (d === 0) return [[meq[b], 1], [meq[b + 2], -s.dy]];
if (d === 1) return [[meq[b + 1], 1], [meq[b + 2], s.dx]];
return [[meq[b + 2], 1]];
};
Skyline Cholesky factorization
64-bit
Factorizes the stiffness of the whole model once for every load set.
K = L Lᵀ; L_ii = √(K_ii − Σ L_ik²)
In place in a Float64Array, without pivoting. The run stops as a mechanism when a reduced diagonal is not above 1e-11 × the diagonal it started from, or when that diagonal is not positive.
frame/plate-solve.js
lines 205 to 223
CadCraft 0.8.2
// skyline Cholesky in place
for (let i = 0; i < n; i++) {
const fi = first[i];
for (let j = fi; j <= i; j++) {
const fj = first[j];
let s = K[at(i, j)];
const kii = s;
for (let k = Math.max(fi, fj); k < j; k++) s -= K[at(i, k)] * K[at(j, k)];
if (j < i) K[at(i, j)] = s / K[at(j, j)];
else {
if (!(s > X.PIVOT_TOL * Math.abs(kii)) || !(kii > 0)) {
const w = dofOf[i];
if (w.q !== undefined) return { error: { code: 'mechanism', node: model.nodes[Math.floor(w.q / 6)].id, dof: w.q % 6 }, warnings };
return { error: { code: 'mechanism', node: model.nodes[w.floor.nodes[0]].id, dof: [0, 1, 5][w.a], diaphragm: w.floor.id }, warnings };
}
K[at(i, i)] = Math.sqrt(s);
}
}
}
Forward and back substitution
64-bit
Solves the factorized system for one load vector.
L y = f; Lᵀ x = y
Two sweeps over the skyline in a Float64Array, each a running subtraction of products followed by one division by the diagonal.
frame/plate-solve.js
lines 228 to 234
CadCraft 0.8.2
function solveSky(prep, f) {
const { n, L } = prep, { first, at } = prep.shell;
const y = Float64Array.from(f);
for (let i = 0; i < n; i++) { let s = y[i]; for (let k = first[i]; k < i; k++) s -= L[at(i, k)] * y[k]; y[i] = s / L[at(i, i)]; }
for (let i = n - 1; i >= 0; i--) { y[i] /= L[at(i, i)]; const xi = y[i]; for (let k = first[i]; k < i; k++) y[k] -= L[at(i, k)] * xi; }
return y;
}
Loads on plates
64-bit
Turns a plate's own weight and area loads into forces at its corners.
corner force = (∫ Nᵢ dA) × p; own weight p = −factor × f × ρ × 1e-9 × g × t
The corner shares come from the element's 2 × 2 Gauss integral. Density in kg/m³ becomes kg/mm³ through 1e-9 and g is 9.80665.
frame/plate-solve.js
lines 242 to 251
CadCraft 0.8.2
const plateLoads = (k, vec) => { const p = pl[k]; p.nodes.forEach((nd, c) => { for (let a = 0; a < 3; a++) nodeF[6 * nd + a] += p.ca[c] * vec[a]; }); };
for (const { c, f } of set.entries) {
if (c.selfWeight) pl.forEach((p, k) => { const w = -c.selfWeight * f * p.rho * 1e-9 * G_ACC * p.t; if (w) plateLoads(k, [0, 0, w]); });
for (const l of c.loads) {
if (l.type !== 'area') continue;
const k = pIdx.get(l.plate), ax = pl[k].el.ax;
const g = l.dir === 'X' ? [1, 0, 0] : l.dir === 'Y' ? [0, 1, 0] : l.dir === 'Z' ? [0, 0, 1] : ax.ez;
plateLoads(k, g.map((x) => x * l.p * f));
}
}
Node displacements from the solution
64-bit
Gives every node its six displacements, including those tied to a floor's master point.
u_q = Σ c × x_e over the constraint map of q
A tied movement is the master's movement plus an offset times the master's rotation, summed in 64-bit. The solution is used as it comes from the one solve: no residual is formed and no refinement step follows, unlike the frame-only path.
frame/plate-solve.js
lines 273 to 275
CadCraft 0.8.2
const x = solveSky(prep, F);
const u = new Float64Array(6 * nn);
for (let q = 0; q < 6 * nn; q++) { let s = 0; for (const [e, c] of maps[q]) s += c * x[e]; u[q] = s; }
Plate forces into the reactions
64-bit
Adds each plate's corner forces to the support reactions and stores the plate's center results.
f = K_g u, added at the corners; center results at ξ = η = 0
The center results are stored eight per plate in a Float64Array. The expression r / 6 | 0 in this code is an array index, not a truncated value.
frame/plate-solve.js
lines 292 to 302
CadCraft 0.8.2
const plates = new Float64Array(8 * pl.length);
const plateU = [];
pl.forEach((p, k) => {
const ug = new Float64Array(24);
p.nodes.forEach((nd, c) => { for (let d = 0; d < 6; d++) ug[6 * c + d] = u[6 * nd + d]; });
const kg = p.el.kg;
for (let r = 0; r < 24; r++) { let s = 0; for (let q = 0; q < 24; q++) s += kg[r * 24 + q] * ug[q]; reactions[6 * p.nodes[r / 6 | 0] + (r % 6)] += s; }
const ul = toLocal(p.el, ug);
plateU.push(ul);
plates.set(plateResultants(p.el, ul, 0, 0), 8 * k);
});
Plate results at the nodes
64-bit
The moments and forces shown at each node of a slab or wall.
inside: mean of the corner values of the plates at the node; on an edge: a + b x + c y fitted through nearby Gauss values, read at the node
These are smoothed values, not the element's own. Plates are grouped by their normal rounded to 6 decimals; plates with the same rounded normal are averaged together at a shared node.
frame/plate-solve.js
lines 343 to 382
CadCraft 0.8.2
export function plateNodeResults(prep, R) {
const pl = prep.shell.pl;
const key = (ax) => ax.ez.map((v) => Math.round(v * 1e6)).join(',');
const groups = new Map();
pl.forEach((p, k) => { const kk = key(p.el.ax); if (!groups.has(kk)) groups.set(kk, []); groups.get(kk).push(k); });
const gauss = pl.map((p, k) => plateGaussResultants(p.el, R.plateU[k]));
const out = new Map();
for (const [kk, list] of groups) {
const at = new Map(), edges = new Map();
for (const k of list) pl[k].nodes.forEach((nd, c) => {
if (!at.has(nd)) at.set(nd, []);
at.get(nd).push([k, c]);
const nx = pl[k].nodes[(c + 1) % 4], e = Math.min(nd, nx) + ',' + Math.max(nd, nx);
edges.set(e, (edges.get(e) || 0) + 1);
});
const onEdge = new Set();
for (const [e, n] of edges) if (n === 1) e.split(',').forEach((x) => onEdge.add(Number(x)));
for (const [nd, uses] of at) {
let values;
if (!onEdge.has(nd)) {
values = new Array(8).fill(0);
for (const [k, c] of uses) { const v = extrapolate(gauss[k])[c]; for (let a = 0; a < 8; a++) values[a] += v[a] / uses.length; }
} else {
const ring = new Set();
for (const [k] of uses) for (const n2 of pl[k].nodes) for (const [k2] of at.get(n2)) ring.add(k2);
const ax = pl[uses[0][0]].el.ax, P0 = prep.P[nd];
const pts = [], vals = [];
for (const k of ring) {
const a2 = pl[k].el.ax, d = [a2.c[0] - P0[0], a2.c[1] - P0[1], a2.c[2] - P0[2]];
const ox = d[0] * ax.ex[0] + d[1] * ax.ex[1] + d[2] * ax.ex[2], oy = d[0] * ax.ey[0] + d[1] * ax.ey[1] + d[2] * ax.ey[2];
gaussPoints(a2).forEach((g, q) => { pts.push([ox + g[0], oy + g[1]]); vals.push(gauss[k][q]); });
}
values = planeFit(pts, vals);
}
if (!out.has(nd)) out.set(nd, []);
out.get(nd).push({ axes: kk, values, plates: uses.length, edge: onEdge.has(nd) });
}
}
return out;
}
Least-squares plane through Gauss values
64-bit
The value at a node on a slab's edge, from a plane fitted through the integration point values around it.
normal equations 3 × 3, first unknown by Cramer's rule
Coordinates are taken relative to the node, which keeps the sums small. The normal equations are solved by determinants; the fit falls back to the plain mean only when the determinant is below 1e-300, an absolute limit that a nearly singular fit does not reach.
frame/plate-solve.js
lines 384 to 396
CadCraft 0.8.2
/** Least-squares a + b x + c y per component through the points; returns a (the value at the origin). */
function planeFit(pts, vals) {
const S = [0, 0, 0, 0, 0, 0, 0, 0, 0];
for (const [x, y] of pts) { const r = [1, x, y]; for (let i = 0; i < 3; i++) for (let j = 0; j < 3; j++) S[i * 3 + j] += r[i] * r[j]; }
const det = (m) => m[0] * (m[4] * m[8] - m[5] * m[7]) - m[1] * (m[3] * m[8] - m[5] * m[6]) + m[2] * (m[3] * m[7] - m[4] * m[6]);
const D = det(S);
return vals[0].map((_, c) => {
const b = [0, 0, 0];
pts.forEach(([x, y], k) => { const v = vals[k][c]; b[0] += v; b[1] += x * v; b[2] += y * v; });
if (Math.abs(D) < 1e-300) return b[0] / pts.length;
return det([b[0], S[1], S[2], b[1], S[4], S[5], b[2], S[7], S[8]]) / D; // Cramer: the first unknown
});
}
Story drift
64-bit
How far each rigid floor moves relative to the floor below, at its center and at its worst node.
drift = move(floor) − move(floor below) at the same point; ratio = |drift| ÷ story height
The drift is a difference of two floors' movements, each a master movement plus offset × rotation. The largest value is searched over the floor's nodes. The ratio is left out when the story height is not positive.
frame/plate-solve.js
lines 398 to 424
CadCraft 0.8.2
/** Storey drifts of the rigid floors. */
export function storyDrifts(prep, R) {
const floors = [...prep.shell.floors].sort((a, b) => a.z - b.z);
const supZ = prep.model.supports.map((s) => prep.P[prep.idx.get(s.node)][2]);
const base = supZ.length ? Math.min(...supZ) : 0;
const move = (f, x, y) => { // a floor's rigid in-plan movement at (x, y)
if (!f) return [0, 0];
const b = 3 * f.k, Ux = R.masters[b], Uy = R.masters[b + 1], Rz = R.masters[b + 2];
return [Ux - (y - f.ym) * Rz, Uy + (x - f.xm) * Rz];
};
return floors.map((f, k) => {
const below = k ? floors[k - 1] : null;
const h = f.z - (below ? below.z : base);
const at = (x, y) => { const a = move(f, x, y), b = move(below, x, y); return [a[0] - b[0], a[1] - b[1]]; };
const c = at(f.xm, f.ym);
let maxX = { value: 0, node: null }, maxY = { value: 0, node: null };
for (const nd of f.nodes) {
const d = at(prep.P[nd][0], prep.P[nd][1]);
if (Math.abs(d[0]) > Math.abs(maxX.value)) maxX = { value: d[0], node: prep.model.nodes[nd].id };
if (Math.abs(d[1]) > Math.abs(maxY.value)) maxY = { value: d[1], node: prep.model.nodes[nd].id };
}
const b = 3 * f.k;
return { id: f.id, name: f.name, z: f.z, h, below: below ? below.id : null, U: [R.masters[b], R.masters[b + 1], R.masters[b + 2]], master: [f.xm, f.ym],
dx: c[0], dy: c[1], ratioX: h > 0 ? Math.abs(c[0]) / h : null, ratioY: h > 0 ? Math.abs(c[1]) / h : null,
maxX: { ...maxX, ratio: h > 0 ? Math.abs(maxX.value) / h : null }, maxY: { ...maxY, ratio: h > 0 ? Math.abs(maxY.value) / h : null } };
});
}
Combined plate results
64-bit
Adds plate results and floor movements of several load sets, each times its factor.
result = Σ fᵢ × resultᵢ
Sums of products in Float64Array, in the order of the list.
frame/plate-solve.js
lines 319 to 331
CadCraft 0.8.2
/** Plate results added to a combination made by frame-solve's combine(). */
export function combineShell(prep, list, out) {
const np = prep.shell.pl.length;
out.plates = new Float64Array(8 * np);
out.plateU = prep.shell.pl.map(() => new Float64Array(24));
out.masters = new Float64Array(3 * prep.shell.floors.length);
for (const { R, f } of list) {
for (let q = 0; q < out.plates.length; q++) out.plates[q] += f * R.plates[q];
R.plateU.forEach((ul, k) => { for (let q = 0; q < 24; q++) out.plateU[k][q] += f * ul[q]; });
for (let q = 0; q < out.masters.length; q++) out.masters[q] += f * R.masters[q];
}
return out;
}
In detail
Method
Master–slave constraints for rigid diaphragms (the header cites Wilson, Three-Dimensional Static and Dynamic Analysis of Structures); skyline Cholesky; corner values by averaging extrapolated Gauss values, with a least-squares plane on a slab's edge (the header cites Zienkiewicz and Zhu 1992).
Number types
64-bit numbers and Float64Array; the skyline's row pointers are whole numbers held in a Float64Array.
Tolerances
Pivot: s ≤ 1e-11 × the diagonal entry stops the run as a mechanism (relative). Poisson's ratio must lie between −0.999 and 0.4999. Warp is reported above 0.01 × √area. Plates share a result group when their normals agree to 6 decimals. The plane fit falls back to the mean when its determinant is below 1e-300 (absolute).
Limits
The plate path solves once: the refinement step of the frame-only path is not run here. Slab design reads plate center values, which change with mesh size.
Approach: 64-bit, stated limit
The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.