Concrete member and slab checks on a solved frame

/** Capacities a section needs once (both axes, both signs, the P–M diagrams). */
function sectionSetup(sec, rc) {
  const mat = { fc: rc.fc, fy: rc.fy };
  const flex = {};
  for (const axis of ['y', 'z']) for (const sign of [1, -1]) {
    const r = ACI.beamFlexure({ b: sec.b, h: sec.h, bars: sec.bars, axis, sign, ...mat });
    flex[axis + (sign > 0 ? '+' : '-')] = r.ok ? r : null;
  }
  const D = { y: ACI.interaction({ b: sec.b, h: sec.h, bars: sec.bars, axis: 'y', spiral: sec.spiral, ...mat }), z: ACI.interaction({ b: sec.b, h: sec.h, bars: sec.bars, axis: 'z', spiral: sec.spiral, ...mat }) };
  // shear along local z is carried over the depth h (the layers of bending about y), along y over b; d = the centroid of the
  // bars in the far half (the tension side of a sagging beam), 0.8 × depth when that half has none
  const shear = {};
  for (const [dir, axis] of [['z', 'y'], ['y', 'z']]) {
    const L = layersFor(sec, axis, 1), depth = L.depth;
    const t = L.layers.filter((l) => l.d > depth / 2);
    const As = t.reduce((s, l) => s + l.As, 0);
    const d = As ? t.reduce((s, l) => s + l.As * l.d, 0) / As : 0.8 * depth;
    shear[dir] = { bw: L.width, d, h: depth, As };
  }
  return { flex, D, shear, Ag: sec.b * sec.h, Tth: ACI.torsionThreshold({ b: sec.b, h: sec.h, fc: rc.fc, lambda: rc.lambda }) };
}

/** φMn at a given Pu on a φ diagram (the outline's largest moment at that axial load). */
function momentAt(D, Pu) {
  const P = D.points;
  let best = 0;
  for (let k = 0; k + 1 < P.length; k++) {
    const a = P[k], b = P[k + 1];
    if ((a.phiPn - Pu) * (b.phiPn - Pu) > 0) continue;
    const t = Math.abs(b.phiPn - a.phiPn) < 1e-12 ? 0 : (Pu - a.phiPn) / (b.phiPn - a.phiPn);
    best = Math.max(best, a.phiMn + t * (b.phiMn - a.phiMn));
  }
  return best;
}

/**
 * Non-sway moment magnification of one axis of a column for one combination (ACI 318-19 6.6.4.5). Ms = the moments along the member (N·mm, sign as the diagram's:
 * the same sign at both ends is single curvature), Pu the largest compression (N), dim the section depth in the plane of bending (mm), Ig that axis's gross inertia.
 *   → { slender: false } · { slender: true, unstable, delta, Mc, Pc, Cm, M2min, kLr, limit }
 */
export function magnifyNonSway({ Ms, Pu, dim, Ig, L, rc }) {
  const M1 = Ms[0], M2 = Ms[Ms.length - 1];
  const big = Math.abs(M1) > Math.abs(M2) ? M1 : M2, small = big === M1 ? M2 : M1;
  const ratio = Math.abs(big) > 1e-9 ? small / big : 1;
  const kLr = rc.k * L / (0.3 * dim), limit = Math.min(34 + 12 * ratio, 40);
  if (kLr <= limit) return { slender: false, kLr, limit };
  const Ec = 4700 * Math.sqrt(rc.fc), EI = 0.4 * Ec * Ig / (1 + rc.betaDns), Pc = Math.PI ** 2 * EI / (rc.k * L) ** 2;
  if (Pu >= 0.75 * Pc) return { slender: true, unstable: true, Pc, kLr, limit };
  const absBig = Math.abs(big), M2min = Math.max(Pu, 0) * (15 + 0.03 * dim);
  const mid = Ms[Math.floor(Ms.length / 2)], span = Math.max(absBig, Math.abs(mid), 1e-9);
  const transverse = Math.abs(mid - (M1 + M2) / 2) > 0.05 * span;   // the diagram is not a straight line between the ends: loads between the supports (6.6.4.5.3b)
  const Cm = transverse || absBig < M2min ? 1 : Math.max(0.4, 0.6 + 0.4 * ratio);
  const delta = Math.max(1, Cm / (1 - Pu / (0.75 * Pc)));
  return { slender: true, unstable: false, delta, Mc: delta * Math.max(absBig, M2min), M2: absBig, Pc, Cm, M2min, kLr, limit };
}

if (msec && Math.abs(msec.A - S.Ag) > 0.05 * S.Ag) partial.add('area-mismatch');
// the drawn section's inertias against the rectangle (b along local y, h along local z): gross, or the cracked 0.70 / 0.35 of ACI 6.6.3.1.1
const likeIg = (I, Ig) => [1, 0.7, 0.35].some((f) => Math.abs(I - f * Ig) <= 0.05 * f * Ig);
if (msec && !partial.has('area-mismatch') && !(likeIg(msec.Iy, sec.b * sec.h ** 3 / 12) && likeIg(msec.Iz, sec.h * sec.b ** 3 / 12))) partial.add('axes-mismatch');

if (role === 'beam') {
  const fy = S.flex['y' + (My >= 0 ? '+' : '-')], fz = S.flex['z' + (Mz >= 0 ? '+' : '-')];
  const ry = fy ? Math.abs(My) / fy.phiMn : Infinity, rz = fz ? Math.abs(Mz) / fz.phiMn : Infinity;
  keep('rc-flexure', 'ACI 318-19 22.2 + Table 21.2.2', ry + (Math.abs(Mz) > 1e-9 ? rz : 0), 1, x);

const uy = ACI.pmRatio(S.D.y, Pu, My), uz = ACI.pmRatio(S.D.z, Pu, Mz);
const both = Math.min(Math.abs(My), Math.abs(Mz)) > 0.01 * Math.max(Math.abs(My), Math.abs(Mz), 1e-9);
let r, clause;
if (!both) { r = Math.abs(My) >= Math.abs(Mz) ? uy : uz; clause = 'ACI 318-19 22.4 + Table 21.2.2'; }
else if (Pu >= 0.1 * rc.fc * S.Ag) { r = Math.max(rayUncapped(S.D.y, Pu, My) + rayUncapped(S.D.z, Pu, Mz) - Pu / (S.D.y.points[0].phi * S.D.y.Po), Pu / S.D.y.phiPnMax); clause = 'ACI 318-19 R22.4 (Bresler 1960)'; }
// i18n-ignore: clause names of the design code, the same in every language
else { r = Math.abs(My) / Math.max(momentAt(S.D.y, Pu), 1e-9) + Math.abs(Mz) / Math.max(momentAt(S.D.z, Pu), 1e-9); clause = 'ACI 318-19 R22.4 (linear, low axial load)'; }
keep('rc-pm', clause, r, 1, x);

/** Bottom and top design moments of the two layers for the plate moments Mx, My, Mxy (all per unit width, any consistent units). */
export function woodArmer(Mx, My, Mxy) {
  const t = Math.abs(Mxy);
  const layer = (a, b) => {                                    // a, b = the moments in the direction of each layer, tension on the face this layer is on
    let xa = a + t, yb = b + t;
    if (xa < 0) { yb = b + (t * t) / Math.abs(a); xa = 0; }     // here a < −t ≤ 0, so a is not zero
    else if (yb < 0) { xa = a + (t * t) / Math.abs(b); yb = 0; }
    return [Math.max(0, xa), Math.max(0, yb)];
  };
  const [bx, by] = layer(Mx, My), [tx, ty] = layer(-Mx, -My);
  return { bx, by, tx, ty };
}

export function slabDesign(model, sol, slabIn, { sets = null } = {}) {
  const slab = cleanSlab(slabIn);
  if (!slab) return { ok: false, error: 'no-slab' };
  if (!sol?.ok || !model.plates?.length) return { ok: false, error: 'no-plates' };
  const combos = (sets || model.combos.filter((c) => c.limit === 'strength').map((c) => c.id)).filter((id) => sol.results['combo:' + id]);
  if (!combos.length) return { ok: false, error: 'no-combos' };
  const prep = sol.prep;
  const fc = slab.fc, fy = slab.fy;
  const rhoMin = fy < 420 ? 0.0020 : Math.max(0.0014, 0.0018 * 420 / fy);
  const out = { ok: true, slab, combos: combos.length, plates: model.plates.length, rows: [], checks: [], max: {} };
  const layers = [['bx', 'bottomX', 0], ['by', 'bottomY', 1], ['tx', 'topX', 0], ['ty', 'topY', 1]];
  const worst = Object.fromEntries(layers.map(([k]) => [k, { M: 0, plate: null, set: null }]));
  let shear = { V: 0, plate: null, set: null, dir: '' };
  const tOf = new Map(model.plates.map((p) => [p.id, p.t]));
  for (const id of combos) {
    const R = sol.results['combo:' + id];
    model.plates.forEach((p, k) => {
      const r = plateCentre(prep, R, k);          // [Nx Ny Nxy Mx My Mxy Qx Qy] per mm of width
      const wa = woodArmer(r[3], r[4], r[5]);
      for (const [key] of layers) if (wa[key] > worst[key].M) worst[key] = { M: wa[key], plate: p.id, set: id };
      for (const [dir, V] of [['x', Math.abs(r[6])], ['y', Math.abs(r[7])]]) if (V > shear.V) shear = { V, plate: p.id, set: id, dir };
    });
  }
  const h = slab.h || Math.max(...model.plates.map((p) => p.t));
  out.h = h;
  for (const [key, field, which] of layers) {
    const b = slab[field], top = key[0] === 't';
    // x layer outermost; for the top layers too (cover to the face)
    const d = which === 0 ? h - slab.cover - b.d / 2 : h - slab.cover - slab[layers[which === 0 ? 1 : 0][1]].d - b.d / 2;
    const Mu = worst[key].M * 1000;       // N·mm per m
    const req = Mu > 0 ? requiredAs({ b: 1000, d, fc, fy, Mu, phi: PHI_FLEX }) : 0;
    const min = rhoMin * 1000 * h;
    const prov = area(b.d) * 1000 / b.s;
    const row = { key, layer: field, bar: b.d, s: b.s, d, Mu: Mu / 1e6, req: req === null ? Infinity : req, min, need: req === null ? Infinity : Math.max(req, min), prov, ratio: req === null ? Infinity : Math.max(req, min) / prov, plate: worst[key].plate, set: worst[key].set,
      clause: 'ACI 318-19 22.3 + 7.6.1.1', top };
    out.rows.push(row);
    out.max[key] = row;
  }
  // shear: per metre, at the plate centre of the worst plate (no stirrups in a slab): φVc
  const dAvg = out.rows.reduce((a, r) => a + r.d, 0) / out.rows.length;
  const AsAvg = out.rows.reduce((a, r) => a + r.prov, 0) / out.rows.length;
  const bs = beamShear({ bw: 1000, d: dAvg, h, fc, fyt: 420, Av: 0, s: 0, Vu: shear.V * 1000, As: AsAvg * 1, lambda: slab.lambda });
  out.checks.push({ id: 'slab-shear', clause: 'ACI 318-19 22.5.5.1', demand: shear.V * 1000, capacity: bs.phiVn, ratio: bs.phiVn > 0 ? shear.V * 1000 / bs.phiVn : Infinity, unit: 'N/m', plate: shear.plate, set: shear.set, dir: shear.dir, d: dAvg });
  out.util = Math.max(...out.rows.map((r) => r.ratio), out.checks[0].ratio);
  out.thin = model.plates.some((p) => tOf.get(p.id) < h - 1e-6 || tOf.get(p.id) > h + 1e-6);
  void PHI_SHEAR;
  return out;
}

export function concreteMaterial(mat) { const E = 4700 * Math.sqrt(mat.fc); return { E, G: E / 2.4, rho: 2400 }; }

export const stiffnessFactor = (mat, vertical) => (mat.cracked ? (vertical ? 0.7 : 0.35) : 1);

const area = (d) => Math.PI * d * d / 4;