Settlement, flexible mats, pile groups and retaining walls

/** The settlement (mm) at the corner of a flexible rectangle b × l (m) carrying q (kPa) on an elastic half-space (E MPa, ν): Timoshenko and Goodier, m = l/b ≥ 1. */
const cornerSettlement = (b, l, q, E, nu) => {
  const m = Math.max(l, b) / Math.min(l, b), sh = Math.min(l, b), r = Math.sqrt(m * m + 1);
  return q * sh * (1 - nu * nu) / (Math.PI * E * 1000) * (m * Math.log((1 + r) / m) + Math.log(m + r)) * 1000;
};

/**
 * Elastic settlement of a loaded L × B (mm) rectangle at pressure q (kPa) on a half-space with modulus Es (MPa) and Poisson's ratio nu: the centre
 * is four corners of the quarter rectangles, the corner is one corner of the whole. Flexible loading, so the centre is an upper bound for a stiff
 * footing's average; a layer of limited depth, the overburden relief and consolidation are left out (the soil is yours).
 */
export function settlement(L, B, q, Es, nu = 0.3) {
  const l = L / 1000, b = B / 1000;
  const corner = cornerSettlement(b, l, q, Es, nu), centre = 4 * cornerSettlement(b / 2, l / 2, q, Es, nu);
  return { centre, corner };
}

const smallest = Math.min(...withLoads.map((l) => Math.min(l.cx || 300, l.cy || 300)));
let size = opts.size || Math.min(Math.min(600, Math.max(150, Math.min(L, B) / 14)), Math.max(100, smallest / 3));
size = Math.max(size, Math.sqrt(L * B / 8000));
const Ec = 4700 * Math.sqrt(fc), nu = 0.2;

  const A = Math.abs((q[2].x - q[0].x) * (q[3].y - q[1].y) - (q[3].x - q[1].x) * (q[2].y - q[0].y)) / 2;
  for (const i of p.nodes) trib.set(i, (trib.get(i) || 0) + A / 4);
}
const ks = soil.ks * 1e-6;   // kN/m³ → N/mm³
json.supports = json.nodes.map((n) => ({ node: n.id, fix: [true, true, false, false, false, true], k: [0, 0, ks * trib.get(n.id), 0, 0, 0] }));

// a column's load on the nodes inside its footprint: P/n plus the moment shares, as a pile group would take them
const spread = (l, set) => {
  const hw = (l.cx || 300) / 2 + 1, hh = (l.cy || 300) / 2 + 1;
  let at = json.nodes.filter((n) => Math.abs(n.x - l.x) <= hw && Math.abs(n.y - l.y) <= hh);
  if (!at.length) { let best = null; for (const n of json.nodes) { const d = Math.hypot(n.x - l.x, n.y - l.y); if (!best || d < best.d) best = { n, d }; } at = [best.n]; }
  const sx2 = at.reduce((a, n) => a + ((n.x - l.x) / 1000) ** 2, 0), sy2 = at.reduce((a, n) => a + ((n.y - l.y) / 1000) ** 2, 0);
  return at.map((n) => {
    const dx = (n.x - l.x) / 1000, dy = (n.y - l.y) / 1000;
    const R = (set.P || 0) / at.length + (sx2 > 1e-12 ? (set.My || 0) * dx / sx2 : 0) - (sy2 > 1e-12 ? (set.Mx || 0) * dy / sy2 : 0);   // kN, down; My loads the +x side more, Mx the −y side
    return { type: 'node', node: n.id, F: [0, 0, -R * 1000, 0, 0, 0] };
  });
};

let w = 0, q = 0, wMin = Infinity;
const a = new Map();
for (const n of json.nodes) { const k = idx.get(n.id), uz = -R.service.u[6 * k + 2]; w = Math.max(w, uz); wMin = Math.min(wMin, uz); q = Math.max(q, soil.ks * uz / 1000); a.set(n.id, uz); }
out.settleMax = w; out.settleMin = wMin; out.qMax = q;
out.reaction = json.nodes.reduce((t, n) => t + ks * trib.get(n.id) * a.get(n.id), 0) / 1000;   // kN: the springs' total, equal to the service loads

let mxP = 0, mxN = 0, myP = 0, myN = 0, qv = 0;
const dd = (H - f.mat.cover - BARS.get(f.mat.barX).db) ;
const nearColumn = (n) => withLoads.some((l) => Math.abs(n.x - l.x) < (l.cx || 300) / 2 + dd && Math.abs(n.y - l.y) < (l.cy || 300) / 2 + dd);
for (const n of json.nodes) {
  const v = nr.get(idx.get(n.id))?.[0]?.values || nr.get(idx.get(n.id))?.values;
  if (!v) continue;
  mxP = Math.max(mxP, v[3]); mxN = Math.min(mxN, v[3]); myP = Math.max(myP, v[4]); myN = Math.min(myN, v[4]);
  if (!nearColumn(n)) qv = Math.max(qv, Math.hypot(v[6], v[7]));   // shear is read at d from the column faces, not under the load
}
// N·mm/mm → kN·m/m; N/mm → kN/m
const M = { xP: mxP / 1000, xN: mxN / 1000, yP: myP / 1000, yN: myN / 1000, V: qv };

const hog = -lay.neg;
if (hog > 0.02 * Math.max(1, lay.pos)) {
  const cap = PHI_PLAIN * 0.42 * lam0(soil) * Math.sqrt(fc) * (H * H / 6) * 1000 / 1e6;   // kN·m/m

/** The springs' pressure (kPa) at a point (mm): the bilinear field of the node settlements times ks. */
function interp(Rs, model, idx, json, x, y, ks) {
  // nodes lie on a grid-ish mesh: find the plate containing the point and interpolate by its bilinear shape in plate coordinates (a quad here is near rectangular)
  const byId = json._byId || (json._byId = new Map(json.nodes.map((n) => [n.id, n])));
  const hit = json._plates || (json._plates = json.plates.map((p) => { const q = p.nodes.map((i) => byId.get(i)); return { p, q, x0: Math.min(...q.map((n) => n.x)), x1: Math.max(...q.map((n) => n.x)), y0: Math.min(...q.map((n) => n.y)), y1: Math.max(...q.map((n) => n.y)) }; }));
  const e = hit.find((h) => x >= h.x0 - 1e-6 && x <= h.x1 + 1e-6 && y >= h.y0 - 1e-6 && y <= h.y1 + 1e-6);
  if (!e) return 0;
  const u = (n) => -Rs.u[6 * idx.get(n.id) + 2];
  const tx = e.x1 > e.x0 ? (x - e.x0) / (e.x1 - e.x0) : 0, ty = e.y1 > e.y0 ? (y - e.y0) / (e.y1 - e.y0) : 0;
  const at = (cx, cy) => e.q.find((n) => Math.abs(n.x - cx) < 1e-6 && Math.abs(n.y - cy) < 1e-6);
  const g = (cx, cy) => { const n = at(cx, cy); return n ? u(n) : 0; };
  return ks / 1000 * ((1 - tx) * (1 - ty) * g(e.x0, e.y0) + tx * (1 - ty) * g(e.x1, e.y0) + tx * ty * g(e.x1, e.y1) + (1 - tx) * ty * g(e.x0, e.y1));
}

/**
 * The elastic method: every pile takes P/n plus the moment shares, R_i = P/n + Mx·y_i/Σy² + My·x_i/Σx² (a rigid cap, equal piles, with x and y from
 * the group's centroid). `cap` = { compression, tension } in kN a pile may carry (yours: from the pile's design, not from here).
 */
export function pileGroup(piles, load, cap = {}) {
  const n = piles.length;
  const cx = piles.reduce((s, p) => s + p.x, 0) / n, cy = piles.reduce((s, p) => s + p.y, 0) / n;
  const xs = piles.map((p) => (p.x - cx) / 1000), ys = piles.map((p) => (p.y - cy) / 1000);
  const sx2 = xs.reduce((s, v) => s + v * v, 0), sy2 = ys.reduce((s, v) => s + v * v, 0);
  const P = load.P || 0, Mx = load.Mx || 0, My = load.My || 0;
  const reactions = piles.map((p, i) => P / n + (sy2 > 1e-12 ? Mx * ys[i] / sy2 : 0) + (sx2 > 1e-12 ? My * xs[i] / sx2 : 0));
  const checks = [];
  const top = Math.max(...reactions), bot = Math.min(...reactions);
  if (cap.compression > 0) checks.push(mk('pile-compression', 'pile group, elastic method', top, cap.compression, 'kN', { n, R: r1(top), cap: cap.compression, pile: reactions.indexOf(top) + 1 }));
  else checks.push(note('pile-compression', 'pile group, elastic method', { R: r1(top) }));
  if (bot < -1e-9) {
    if (cap.tension > 0) checks.push(mk('pile-tension', 'pile group, elastic method', -bot, cap.tension, 'kN', { R: r1(bot), cap: cap.tension, pile: reactions.indexOf(bot) + 1 }));
    else checks.push({ id: 'pile-tension', clause: 'pile group, elastic method', status: cap.tension === 0 ? 'fail' : 'note', unit: 'kN', demand: -bot, capacity: 0, ratio: Infinity, vars: { R: r1(bot), pile: reactions.indexOf(bot) + 1 } });
  }
  const sum = reactions.reduce((s, v) => s + v, 0);
  return { reactions, sum, centroid: [cx, cy], checks, equilibrium: Math.abs(sum - P) < 1e-6 * Math.max(1, Math.abs(P)) };
}

/**
 * A cantilever wall with a level backfill. wall = { H (retained height of soil above the base's top, mm), hb (base thickness), t (stem thickness,
 * constant), toe, heel (mm), embed (soil over the toe, mm, ignored for passive) }. soil = { gamma (kN/m³), phi (degrees), surcharge (kPa), qa,
 * mu (base friction), fsSlide, fsOt, gammaConc }. Active pressure is Rankine, Ka = tan²(45° − φ/2), on the vertical plane through the heel,
 * with the soil on the heel and the stem as weight; passive pressure is left out.
 */
export function retainingWall(wall, soil = {}) {
  const s = { gamma: 18, phi: 30, surcharge: 0, qa: 0, mu: 0, fsSlide: 1.5, fsOt: 1.5, gammaConc: 24, ...soil };
  const H = wall.H / 1000, hb = wall.hb / 1000, t = wall.t / 1000, toe = wall.toe / 1000, heel = wall.heel / 1000;
  const Bw = toe + t + heel, Ht = H + hb;
  const Ka = Math.tan((45 - s.phi / 2) * Math.PI / 180) ** 2;
  const Pa = 0.5 * Ka * s.gamma * Ht * Ht, Ps = Ka * s.surcharge * Ht;
  const MoPa = Pa * Ht / 3, MoPs = Ps * Ht / 2;
  // weights per metre of wall, arms about the toe's front edge
  const parts = [
    { name: 'base', W: s.gammaConc * Bw * hb, x: Bw / 2 },
    { name: 'stem', W: s.gammaConc * t * H, x: toe + t / 2 },
    { name: 'soil', W: s.gamma * heel * H, x: toe + t + heel / 2 },
    { name: 'surcharge', W: s.surcharge * heel, x: toe + t + heel / 2 },
  ];
  const W = parts.reduce((a, p) => a + p.W, 0), Mr = parts.reduce((a, p) => a + p.W * p.x, 0);
  const Ho = Pa + Ps, Mo = MoPa + MoPs;
  const e = Bw / 2 - (Mr - Mo) / W;   // eccentricity of the resultant from the middle of the base
  const field = pressureField(W, e * 1000, 0, Bw * 1000, 1000);
  const checks = [];
  checks.push(mk('wall-overturning', 'Rankine; overturning about the toe', s.fsOt * Mo, Mr, 'kN·m/m', { fs: r2(Mr / Mo), need: s.fsOt, Mo: r1(Mo), Mr: r1(Mr) }));
  if (s.mu > 0) checks.push(mk('wall-sliding', 'Rankine; sliding on the base', s.fsSlide * Ho, s.mu * W, 'kN/m', { fs: r2(s.mu * W / Ho), need: s.fsSlide, H: r1(Ho), W: r1(W) }));
  else checks.push(note('wall-sliding', 'Rankine; sliding on the base', { H: r1(Ho) }));
  if (!field.ok) checks.push({ id: 'wall-bearing', clause: 'base pressure', status: 'fail', unit: 'kPa', demand: Infinity, capacity: s.qa, ratio: Infinity, vars: { reason: field.reason, e: r2(e) } });
  else if (s.qa > 0) checks.push(mk('wall-bearing', 'base pressure', field.qmax, s.qa, 'kPa', { q: r1(field.qmax), qa: s.qa, e: r2(e), kern: r2(Bw / 6) }));
  else checks.push(note('wall-bearing', 'base pressure', { q: r1(field.qmax) }));
  return { Ka, Pa, Ps, W, Mr, Mo, e, Bw, field, parts, checks };
}