Shakes a building of one mass per floor through a ground acceleration record and reports periods, mode shapes, story drifts, floor accelerations and base shear; the floor stiffness can come from the frame model.
How the precision is kept
All matrices and histories are 64-bit numbers or Float64Array. Periods come from a cyclic Jacobi eigen-solution that stops at an absolute off-diagonal size; the response comes from Newmark average-acceleration stepping, so its accuracy is set by the time step (default 0.01 s). Peak drift, floor acceleration and base shear are tracked at every step, but the peak roof displacement is taken from the thinned history that keeps every fourth step by default. A frame's floor stiffness is the inverse of a flexibility table read from the frame solver under 1000 N floor loads.
Open item. Track the peak roof displacement at every step, as the other peaks are, instead of reading it from the thinned history.
Calculations under this policy
12 calculations in 2 source files. 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 eigenvalues and eigenvectors behind the building's periods and mode shapes.
t = sign(θ) ÷ (|θ| + √(θ² + 1)), θ = (a_qq − a_pp) ÷ (2 a_pq); c = 1 ÷ √(t² + 1); s = t × c
The rotation tangent uses the form that divides by |θ| + √(θ² + 1), which does not subtract nearly equal numbers. The sweep stops when the sum of squared off-diagonal terms falls below 1e-24, an absolute value that is not scaled to the size of the matrix, or after 100 sweeps with no warning.
hazard.js
lines 75 to 92
BimCraft 0.8.1
// -------------------------------------------------------------------------------------------------------------------------------- quake
function jacobi(A, n) { // symmetric eigenproblem, cyclic Jacobi: eigenvalues ascending, eigenvectors as columns
const a = A.map((r) => r.slice()), v = Array.from({ length: n }, (_, i) => Array.from({ length: n }, (_, j) => (i === j ? 1 : 0)));
for (let sweep = 0; sweep < 100; sweep++) {
let off = 0;
for (let i = 0; i < n; i++) for (let j = i + 1; j < n; j++) off += a[i][j] * a[i][j];
if (off < 1e-24) break;
for (let p = 0; p < n; p++) for (let q = p + 1; q < n; q++) {
if (Math.abs(a[p][q]) < 1e-300) continue;
const th = (a[q][q] - a[p][p]) / (2 * a[p][q]), t = Math.sign(th || 1) / (Math.abs(th) + Math.sqrt(th * th + 1)), c = 1 / Math.sqrt(t * t + 1), s = t * c;
for (let k = 0; k < n; k++) { const kp = a[k][p], kq = a[k][q]; a[k][p] = c * kp - s * kq; a[k][q] = s * kp + c * kq; }
for (let k = 0; k < n; k++) { const pk = a[p][k], qk = a[q][k]; a[p][k] = c * pk - s * qk; a[q][k] = s * pk + c * qk; }
for (let k = 0; k < n; k++) { const kp = v[k][p], kq = v[k][q]; v[k][p] = c * kp - s * kq; v[k][q] = s * kp + c * kq; }
}
}
const idx = Array.from({ length: n }, (_, i) => i).sort((x, y) => a[x][x] - a[y][y]);
return { values: idx.map((i) => a[i][i]), vectors: idx.map((i) => v.map((r) => r[i])) };
}
Stiffness matrix of a shear building
64-bit
Builds the tridiagonal floor stiffness matrix from one stiffness per story.
Sums and negations of the story stiffnesses in 64-bit; each diagonal term is at most two additions. Masses and stiffnesses that are not greater than zero are refused.
hazard.js
lines 94 to 101
BimCraft 0.8.1
/** Floors 1…n above a fixed base: m[i] the mass of floor i, k[i] the stiffness of the storey below it (k[0] between the ground and floor 1). */
export function shearBuilding({ m, k }) {
const n = m.length;
if (!n || k.length !== n || m.some((x) => !(x > 0)) || k.some((x) => !(x > 0))) return { ok: false, reason: 'badbuilding' };
const K = Array.from({ length: n }, () => new Array(n).fill(0));
for (let i = 0; i < n; i++) { K[i][i] += k[i]; if (i + 1 < n) { K[i][i] += k[i + 1]; K[i][i + 1] -= k[i + 1]; K[i + 1][i] -= k[i + 1]; } }
return buildingFromK({ m, K, k });
}
Periods, mode shapes and effective mass
64-bit
The natural periods, the mode shapes scaled to 1 at the top floor, the participation factor and the effective mass of each mode.
S = K[i][j] ÷ √(mᵢ mⱼ); ω = √max(λ, 0); T = 2π ÷ ω; Γ = Σ m φ ÷ Σ m φ²; effective mass = (Σ m φ)² ÷ Σ m φ²
A slightly negative eigenvalue is clamped to zero before the square root, and a zero frequency gives an infinite period. Each shape is divided by its top-floor value, or by 1 when that value is zero. The story springs of a full matrix are read back as the row sum of the first row and −K[i][i−1].
hazard.js
lines 102 to 119
BimCraft 0.8.1
/**
* 0.8 EG.47 (bim-sim_PART_3): a building from a full floor stiffness matrix K (N/m, n × n, symmetric, floors bottom first) and floor masses m (kg), for example the
* frame solver's condensed stiffness (hazard-frame.js). Same shape as shearBuilding's, with `full: true` so quakeRun solves the dense system. k = the storey
* springs read from K (for the base shear and the drift check).
*/
export function buildingFromK({ m, K, k = null }) {
const n = m.length;
const S = K.map((row, i) => row.map((x, j) => x / Math.sqrt(m[i] * m[j])));
const { values, vectors } = jacobi(S, n);
const modes = values.map((l, j) => {
const shape = vectors[j].map((x, i) => x / Math.sqrt(m[i]));
const top = shape[n - 1] || 1, sh = shape.map((x) => x / top); // top floor = 1
const w = Math.sqrt(Math.max(l, 0));
const L = m.reduce((s, mi, i) => s + mi * sh[i], 0), Mn = m.reduce((s, mi, i) => s + mi * sh[i] * sh[i], 0);
return { omega: w, period: w > 0 ? 2 * Math.PI / w : Infinity, shape: sh, gamma: L / Mn, effMass: L * L / Mn };
});
return { ok: true, n, m, k: k || m.map((_, i) => (i === 0 ? K[0].reduce((s, x) => s + x, 0) : -K[i][i - 1])), K, modes, totalMass: m.reduce((s, x) => s + x, 0), full: k === null };
}
Story stiffness for a chosen first period
64-bit
The single story stiffness that gives the building the first period the person typed.
k = (2π ÷ T₁)² ÷ ω₁², ω₁ = first frequency with k = 1
Solves the unit-stiffness building once and scales by the squared frequency ratio, which is exact for the model because frequency grows with √k. No guard on T₁ = 0.
hazard.js
lines 120 to 124
BimCraft 0.8.1
/** The storey stiffness (the same for every storey) that gives these masses the first period T1. */
export function uniformFor(m, T1) {
const unit = shearBuilding({ m, k: m.map(() => 1) });
return (2 * Math.PI / T1) ** 2 / unit.modes[0].omega ** 2; // ω ∝ √k, so k = (ω wanted / ω for k = 1)²
}
Made ground acceleration record
64-bit
A Ricker pulse or filtered noise record scaled to the peak ground acceleration asked for.
Ricker: a = (1 − 2x²) × e^(−x²), x = π f (t − 1.5 ÷ f); then a × pga ÷ max|a|
The record is a Float64Array with sample times i × dt. The noise record is integrated with a semi-implicit Euler step of size dt and its random numbers come from a 32-bit integer generator; the final scaling to the peak is one 64-bit multiplication per sample.
hazard.js
lines 129 to 146
BimCraft 0.8.1
/** A ground acceleration record (m/s²) with peak `pga`: 'ricker' (a pulse of peak frequency f, centred at 1.5/f) or 'noise' (Kanai–Tajimi filtered noise, ground frequency f, with a rise–hold–fall envelope). */
export function groundMotion({ kind = 'ricker', pga = 2, duration = 20, dt = 0.01, f = 1.5, seed = 1 } = {}) {
const n = Math.floor(duration / dt) + 1, a = new Float64Array(n);
if (kind === 'ricker') {
const t0 = 1.5 / f;
for (let i = 0; i < n; i++) { const x = Math.PI * f * (i * dt - t0), e = x * x; a[i] = (1 - 2 * e) * Math.exp(-e); }
} else {
const rnd = mulberry(seed), wg = 2 * Math.PI * f, zg = 0.6, w = new Float64Array(n);
for (let i = 0; i < n; i++) w[i] = gauss(rnd);
let x = 0, v = 0; // ground filter: x'' + 2 zg wg x' + wg² x = w (the filtered output is the acceleration)
for (let i = 0; i < n; i++) { const acc = w[i] * wg * 2 - 2 * zg * wg * v - wg * wg * x; v += acc * dt; x += v * dt; a[i] = -(2 * zg * wg * v + wg * wg * x); }
for (let i = 0; i < n; i++) { const t = i * dt, rise = duration * 0.15, fall = duration * 0.6; a[i] *= t < rise ? t / rise : t < fall ? 1 : Math.max(0, 1 - (t - fall) / (duration - fall)); }
}
let pk = 0;
for (const v of a) pk = Math.max(pk, Math.abs(v));
if (pk > 0) for (let i = 0; i < n; i++) a[i] *= pga / pk;
return a;
}
Newmark time stepping
64-bit
The floor displacements, velocities and accelerations at every time step of the record.
K* = K + c₁ M + c₂ C, c₁ = 1 ÷ (β Δt²), c₂ = γ ÷ (β Δt), β = 1/4, γ = 1/2; C = a₀ M + a₁ K
State vectors are Float64Array. A shear building is solved by the Thomas algorithm without pivoting; a full matrix is factored once by Cholesky with each pivot floored at 1e-300 instead of being reported as not positive definite. The answer's accuracy is set by the time step, default 0.01 s.
hazard.js
lines 148 to 207
BimCraft 0.8.1
/** Newmark average-acceleration time history of M u'' + C u' + K u = −M 1 ag(t), Rayleigh damping `zeta` at the first and second (or highest) modes. */
export function quakeRun(b, accel, { dt = 0.01, zeta = 0.05, keep = 4 } = {}) {
if (!b?.ok) return { ok: false, reason: 'badbuilding' };
const n = b.n, w1 = b.modes[0].omega, w2 = b.modes[Math.min(1, n - 1)].omega;
let a0 = 0, a1 = 0;
if (n === 1) a0 = 2 * zeta * w1; else { a0 = 2 * zeta * w1 * w2 / (w1 + w2); a1 = 2 * zeta / (w1 + w2); }
const m = b.m, C = b.K.map((row, i) => row.map((x, j) => a1 * x + (i === j ? a0 * m[i] : 0)));
const beta = 0.25, gam = 0.5, c1 = 1 / (beta * dt * dt), c2 = gam / (beta * dt);
// effective stiffness K* = K + c1 M + c2 C. A shear building is tridiagonal (Thomas algorithm); a full matrix (buildingFromK, e.g. a frame's condensed
// stiffness) is factored once (Cholesky: K* is symmetric positive definite) and solved by substitution each step.
const diag = new Float64Array(n), up = new Float64Array(n), lo = new Float64Array(n);
let solve;
if (b.full) {
const L = Array.from({ length: n }, () => new Float64Array(n));
const Ks = b.K.map((row, i) => row.map((x, j) => x + c2 * C[i][j] + (i === j ? c1 * m[i] : 0)));
for (let i = 0; i < n; i++) for (let j = 0; j <= i; j++) {
let sum = Ks[i][j];
for (let q = 0; q < j; q++) sum -= L[i][q] * L[j][q];
if (i === j) L[i][i] = Math.sqrt(Math.max(sum, 1e-300)); else L[i][j] = sum / L[j][j];
}
solve = (rhs) => {
const y = new Float64Array(n), x = new Float64Array(n);
for (let i = 0; i < n; i++) { let sum = rhs[i]; for (let q = 0; q < i; q++) sum -= L[i][q] * y[q]; y[i] = sum / L[i][i]; }
for (let i = n - 1; i >= 0; i--) { let sum = y[i]; for (let q = i + 1; q < n; q++) sum -= L[q][i] * x[q]; x[i] = sum / L[i][i]; }
return x;
};
} else {
for (let i = 0; i < n; i++) { diag[i] = b.K[i][i] + c1 * m[i] + c2 * C[i][i]; if (i + 1 < n) { up[i] = b.K[i][i + 1] + c2 * C[i][i + 1]; lo[i + 1] = b.K[i + 1][i] + c2 * C[i + 1][i]; } }
solve = (rhs) => {
const c = new Float64Array(n), d = new Float64Array(n), x = new Float64Array(n);
c[0] = n > 1 ? up[0] / diag[0] : 0; d[0] = rhs[0] / diag[0];
for (let i = 1; i < n; i++) { const den = diag[i] - lo[i] * c[i - 1]; c[i] = i + 1 < n ? up[i] / den : 0; d[i] = (rhs[i] - lo[i] * d[i - 1]) / den; }
x[n - 1] = d[n - 1];
for (let i = n - 2; i >= 0; i--) x[i] = d[i] - c[i] * x[i + 1];
return x;
};
}
const mul = (A, x) => A.map((row) => row.reduce((s, v, j) => s + v * x[j], 0));
let u = new Float64Array(n), v = new Float64Array(n), acc = new Float64Array(n);
for (let i = 0; i < n; i++) acc[i] = -accel[0];
const steps = accel.length, hist = Array.from({ length: n }, () => []), t = [], accAbs = Array.from({ length: n }, () => 0), shearHist = [];
const peakDrift = new Float64Array(n);
let peakShear = 0;
for (let s = 0; s < steps; s++) {
if (s % keep === 0) { t.push(s * dt); for (let i = 0; i < n; i++) hist[i].push(u[i]); }
// base shear = Σ over floors of m (u'' + ag) = the storey 0 spring force
const sh = b.full ? acc.reduce((q, a1, i) => q + m[i] * (a1 + accel[s]), 0) : b.k[0] * u[0]; peakShear = Math.max(peakShear, Math.abs(sh)); if (s % keep === 0) shearHist.push(sh);
for (let i = 0; i < n; i++) { const d = i === 0 ? u[0] : u[i] - u[i - 1]; peakDrift[i] = Math.max(peakDrift[i], Math.abs(d)); accAbs[i] = Math.max(accAbs[i], Math.abs(acc[i] + accel[s])); }
if (s === steps - 1) break;
const ag = accel[s + 1];
const cu = new Float64Array(n), cv = new Float64Array(n);
for (let i = 0; i < n; i++) { cu[i] = c1 * u[i] + (1 / (beta * dt)) * v[i] + (1 / (2 * beta) - 1) * acc[i]; cv[i] = c2 * u[i] + (gam / beta - 1) * v[i] + dt * (gam / (2 * beta) - 1) * acc[i]; }
const rhs = new Float64Array(n), Mcu = new Float64Array(n), Ccv = mul(C, cv);
for (let i = 0; i < n; i++) { Mcu[i] = m[i] * cu[i]; rhs[i] = -m[i] * ag + Mcu[i] + Ccv[i]; }
const un = solve(rhs), vn = new Float64Array(n), an = new Float64Array(n);
for (let i = 0; i < n; i++) { an[i] = c1 * (un[i] - u[i]) - (1 / (beta * dt)) * v[i] - (1 / (2 * beta) - 1) * acc[i]; vn[i] = v[i] + dt * ((1 - gam) * acc[i] + gam * an[i]); }
u = un; v = vn; acc = an;
}
return { ok: true, t, u: hist, dt: dt * keep, shear: shearHist, peak: { drift: Array.from(peakDrift), acc: accAbs, shear: peakShear, top: Math.max(...hist[n - 1].map(Math.abs)) } };
}
Peak drift, floor acceleration, base shear and roof displacement
64-bit
The largest values over the run that the person reads.
drift = uᵢ − uᵢ₋₁; floor acceleration = ü + a_g; base shear = k₀ u₁ (shear building) or Σ m (ü + a_g) (full matrix)
Peak drift, acceleration and base shear are compared at every step. The peak roof displacement in the result is the maximum of the kept history, which holds every fourth step by default, so a peak between kept steps is missed. The shear-building base shear is the bottom spring force alone.
hazard.js
lines 193 to 195
BimCraft 0.8.1
// base shear = Σ over floors of m (u'' + ag) = the storey 0 spring force
const sh = b.full ? acc.reduce((q, a1, i) => q + m[i] * (a1 + accel[s]), 0) : b.k[0] * u[0]; peakShear = Math.max(peakShear, Math.abs(sh)); if (s % keep === 0) shearHist.push(sh);
for (let i = 0; i < n; i++) { const d = i === 0 ? u[0] : u[i] - u[i - 1]; peakDrift[i] = Math.max(peakDrift[i], Math.abs(d)); accAbs[i] = Math.max(accAbs[i], Math.abs(acc[i] + accel[s])); }
Sway under static floor forces
64-bit
Floor displacements and story drifts under static sideways forces, used for the wind sway.
shear building: driftᵢ = (Σ F above) ÷ kᵢ, u = running sum of drifts; full matrix: K u = F
The full matrix is solved by Gauss elimination with partial pivoting and is called singular when a pivot is below 1e-12, an absolute value in N/m that does not scale with the stiffness. The shear building needs no solve: story shears are summed from the top and drifts from the bottom.
hazard.js
lines 209 to 231
BimCraft 0.8.1
/** Floor displacements under static forces F (N) at the floors: K u = F by the same tridiagonal elimination. */
export function staticSway(b, F) {
const n = b.n, u = new Float64Array(n);
if (b.full && b.K) { // 0.8 EG.47 (bim-analysis_PART_4): a frame's condensed stiffness is a full matrix (a floor pushed alone moves its neighbours): K u = F by Gauss elimination with partial pivoting
const A = b.K.map((row, i) => [...row, F[i] || 0]);
for (let c = 0; c < n; c++) {
let p = c;
for (let r = c + 1; r < n; r++) if (Math.abs(A[r][c]) > Math.abs(A[p][c])) p = r;
if (Math.abs(A[p][c]) < 1e-12) return { u: new Array(n).fill(NaN), drift: new Array(n).fill(NaN), singular: true };
[A[c], A[p]] = [A[p], A[c]];
for (let r = c + 1; r < n; r++) { const f = A[r][c] / A[c][c]; if (f) for (let k = c; k <= n; k++) A[r][k] -= f * A[c][k]; }
}
for (let r = n - 1; r >= 0; r--) { let s = A[r][n]; for (let k = r + 1; k < n; k++) s -= A[r][k] * u[k]; u[r] = s / A[r][r]; }
return { u: Array.from(u), drift: Array.from(u, (v, i) => v - (i ? u[i - 1] : 0)) };
}
// storey shears from the top down, then drifts
let shear = 0;
const drift = new Float64Array(n);
for (let i = n - 1; i >= 0; i--) { shear += F[i]; drift[i] = shear / b.k[i]; }
let s = 0;
for (let i = 0; i < n; i++) { s += drift[i]; u[i] = s; }
return { u: Array.from(u), drift: Array.from(drift) };
}
Groups the frame's nodes into floors by height and takes the lowest group as the base.
nodes whose heights differ by ≤ tol join one group; floor height = mean of the group
Heights are sorted and chained: a node joins the group when it is within tol (150 mm by default, absolute) of the group's highest node so far. The floor height is the 64-bit mean of the group; the rounded height appears only in a note about a lone node.
hazard-frame.js
lines 29 to 48
BimCraft 0.8.1
export function frameFloors(model, { levels = null, tol = 150 } = {}) {
const nodes = model.nodes, notes = [];
if (!nodes.length) return { floors: [], base: 0, notes };
let zs;
if (levels && levels.length) zs = levels.slice().sort((a, b) => a - b);
else {
const sorted = nodes.map((n) => n.z).sort((a, b) => a - b), groups = [];
for (const z of sorted) { const g = groups[groups.length - 1]; if (g && z - g.hi <= tol) { g.hi = z; g.sum += z; g.n++; } else groups.push({ hi: z, sum: z, n: 1 }); }
zs = groups.map((g) => g.sum / g.n);
}
// the base is the lowest height (the supports' feet); a floor needs two nodes at least (a lone column top is not a floor)
const base = zs[0];
const floors = [];
for (const z of zs.slice(1)) {
const mine = []; nodes.forEach((n, i) => { if (Math.abs(n.z - z) <= tol) mine.push(i); });
if (mine.length >= 2) floors.push({ z, nodes: mine });
else notes.push({ key: 'eg47.fr.note.lone', z: Math.round(z) });
}
return { floors, base, notes };
}
Inverse of the flexibility table
64-bit
Turns the floors' flexibility table into their stiffness matrix.
Gauss–Jordan elimination with partial pivoting on plain 64-bit arrays. It returns nothing only when a pivot is below 1e-300, so a nearly singular table is inverted without a warning and no condition number is checked.
hazard-frame.js
lines 50 to 60
BimCraft 0.8.1
function solveInverse(A) { // Gauss–Jordan with partial pivoting; A is small (one row per floor)
const n = A.length, M = A.map((r, i) => [...r, ...Array.from({ length: n }, (_, j) => (i === j ? 1 : 0))]);
for (let c = 0; c < n; c++) {
let p = c; for (let r = c + 1; r < n; r++) if (Math.abs(M[r][c]) > Math.abs(M[p][c])) p = r;
if (Math.abs(M[p][c]) < 1e-300) return null;
[M[c], M[p]] = [M[p], M[c]];
const d = M[c][c]; for (let j = 0; j < 2 * n; j++) M[c][j] /= d;
for (let r = 0; r < n; r++) if (r !== c) { const f = M[r][c]; if (f) for (let j = 0; j < 2 * n; j++) M[r][j] -= f * M[c][j]; }
}
return M.map((r) => r.slice(n));
}
Floor stiffness of a frame
64-bit
The stiffness the frame gives each floor sideways, and one stiffness per story for the drift check.
F[i][j] = (mean sideways move of floor i under 1000 N on floor j) ÷ 1000 ÷ 1000; K = F⁻¹; k_story = V ÷ drift under forces ∝ height
Displacements are read from the frame solver's Float64Array in millimeters and divided by 1000 and by the 1000 N load. F and K are each made symmetric by averaging the two off-diagonal terms. A story whose drift under the triangular pattern is not greater than zero gets an infinite stiffness.
hazard-frame.js
lines 62 to 87
BimCraft 0.8.1
export function frameStiffness(model, { dir = 'x', levels = null, tol = 150 } = {}) {
const { floors, base, notes } = frameFloors(model, { levels, tol });
if (!floors.length) return { ok: false, reason: 'nofloors', notes };
const ax = AX[dir] ?? 0, P = 1000; // N in total on each floor, spread over its nodes
const idx = nodeIndex(model), cases = floors.map((f, j) => ({ id: '__eg47_' + j, name: 'EG.47 floor ' + (j + 1), kind: 'other', selfWeight: 0,
loads: f.nodes.map((i) => ({ type: 'node', node: model.nodes[i].id, F: [ax === 0 ? P / f.nodes.length : 0, ax === 1 ? P / f.nodes.length : 0, 0, 0, 0, 0] })) }));
const m2 = { ...model, cases, combos: [] };
const sol = solveFrame(m2);
if (!sol.ok) return { ok: false, reason: 'solve', error: sol.error, notes };
const n = floors.length, F = Array.from({ length: n }, () => new Array(n).fill(0));
for (let j = 0; j < n; j++) {
const u = sol.results['case:__eg47_' + j].u;
for (let i = 0; i < n; i++) {
let s = 0; for (const k of floors[i].nodes) s += u[6 * k + ax];
F[i][j] = (s / floors[i].nodes.length) / 1000 / P; // mm → m, per N
}
}
for (let i = 0; i < n; i++) for (let j = i + 1; j < n; j++) { const a = (F[i][j] + F[j][i]) / 2; F[i][j] = F[j][i] = a; } // reciprocity: equal up to round-off
const K = solveInverse(F);
if (!K) return { ok: false, reason: 'mechanism', notes };
for (let i = 0; i < n; i++) for (let j = i + 1; j < n; j++) { const a = (K[i][j] + K[j][i]) / 2; K[i][j] = K[j][i] = a; }
// inverted-triangle pattern → floor forces f_j ∝ height above the base → u = F f → storey shears and drifts
const h = floors.map((f) => (f.z - base) / 1000), f = h.map((x) => x), u = F.map((row) => row.reduce((s, x, j) => s + x * f[j], 0));
const kStorey = floors.map((_, i) => { const V = f.slice(i).reduce((s, x) => s + x, 0), d = u[i] - (i ? u[i - 1] : 0); return d > 0 ? V / d : Infinity; });
return { ok: true, floors, F, K, kStorey, base, notes };
}
Floor masses from the frame
64-bit
The mass on each floor: half of each member's weight at each end, plus a load per square meter over the floor's plan box.
member mass = ρ × A × L × 1e-9; floor area = (x extent ÷ 1000) × (y extent ÷ 1000); m = max(1, self + area × load)
Member length is Math.hypot of the node differences in millimeters and 1e-9 turns kg/m³ × mm³ into kilograms. The floor area is the bounding rectangle of the floor's nodes, not the slab outline, and every floor mass is raised to at least 1 kg.
hazard-frame.js
lines 89 to 103
BimCraft 0.8.1
export function frameMasses(model, floors, { addKgM2 = 0, tol = 150 } = {}) {
const n = floors.length, self = new Array(n).fill(0), at = (z) => { let best = -1, bd = Infinity; floors.forEach((f, i) => { const d = Math.abs(f.z - z); if (d < bd) { bd = d; best = i; } }); return bd <= tol ? best : -2; };
const idx = nodeIndex(model), sec = new Map(model.sections.map((s) => [s.id, s])), mat = new Map(model.materials.map((x) => [x.id, x]));
for (const mb of model.members) {
const a = model.nodes[idx.get(mb.i)], b = model.nodes[idx.get(mb.j)], s = sec.get(mb.section), mt = mat.get(mb.material);
if (!a || !b || !s || !mt) continue;
const L = Math.hypot(b.x - a.x, b.y - a.y, b.z - a.z), kg = mt.rho * (s.A * L) * 1e-9; // kg/m³ × mm³ → kg
for (const nd of [a, b]) { const i = at(nd.z); if (i >= 0) self[i] += kg / 2; }
}
const area = floors.map((f) => {
const xs = f.nodes.map((i) => model.nodes[i].x), ys = f.nodes.map((i) => model.nodes[i].y);
return ((Math.max(...xs) - Math.min(...xs)) / 1000) * ((Math.max(...ys) - Math.min(...ys)) / 1000);
});
return { m: self.map((s, i) => Math.max(1, s + area[i] * addKgM2)), area_m2: area, self };
}
In detail
Method
Mass-normalized stiffness S = K ÷ √(mᵢ mⱼ), cyclic Jacobi rotations for eigenvalues; Rayleigh damping on the first two modes; Newmark β = 1/4, γ = 1/2; Thomas algorithm for a shear building, Cholesky factor for a full stiffness matrix. Frame: flexibility by unit floor loads, Gauss–Jordan inverse with partial pivoting.
Number types
64-bit numbers, plain arrays and Float64Array. The made noise record draws from a 32-bit integer generator and is then scaled in 64-bit to the peak acceleration asked for.
Tolerances
Jacobi stops when the sum of squared off-diagonal terms is below 1e-24 (absolute) or after 100 sweeps; a rotation is skipped when the term is below 1e-300. Cholesky pivots are floored at 1e-300. The static solve calls the matrix singular when a pivot is below 1e-12 (absolute, N/m). The stiffness inverse gives up only when a pivot is below 1e-300. Nodes within 150 mm in height (default) are one floor.
Checked by
tools/bimcraft-test/hazard-test.mjs runs shearBuilding, quakeRun, staticSway and the hazard-frame functions in Node.
Limits
A linear shear building with one sideways degree of freedom per floor. The floor mass added per square meter is spread over the plan bounding box of the floor's nodes, which the code describes as an estimate. None of it is a design check against a code.
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.