Resources Numerical accuracy Calculation policy
Finds the factor of safety of a circular slip surface by Bishop's simplified method or the ordinary method of slices, and searches for the critical circle.
64-bit arithmetic. The factor of safety depends on the number of slices (10 to 400), on an iteration that stops at an absolute change of 1e-5, and on a search over a grid of centers and radii that returns the lowest value tried. The unit weight of water is 9.81 kN/m³.
Open item. Report when the Bishop iteration ends at its 200-pass limit without meeting the tolerance.
6 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.
Tests bin/test-civilcraft-geotech.mjs
The constant behind every pore pressure and every submerged unit weight in this file.
γw = 9.81 kN/m³
A constant rounded to three digits. With standard gravity 9.80665 m/s² and 1000 kg/m³ the value is 9.80665 kN/m³, so pore pressures here are about 0.03 % higher than that.
const G_W = 9.81; // kN/m³, unit weight of waterThe ground, a layer top or the water table at a horizontal position.
Linear interpolation between the two neighboring points, found by bisection on the index. Beyond either end the end value is held level. A vertical step returns the second point's elevation.
const interp = (poly, x) => {
if (typeof poly === 'number') return poly;
const n = poly.length;
if (x <= poly[0].x) return poly[0].z;
if (x >= poly[n - 1].x) return poly[n - 1].z;
let lo = 0, hi = n - 1;
while (hi - lo > 1) { const m = (lo + hi) >> 1; if (poly[m].x <= x) lo = m; else hi = m; }
const a = poly[lo], b = poly[hi];
return b.x === a.x ? b.z : a.z + ((x - a.x) / (b.x - a.x)) * (b.z - a.z);
};The positions where the lower arc of a trial circle crosses the ground line.
f(x) = ground(x) − (zc − √(R² − (x − xc)²))
The span is scanned in 400 equal steps for a change of sign, and each change is refined by 40 bisections. Two crossings inside one scan step cancel and are not found. The square root's argument is clamped at zero.
function circleGround(geom, c) {
// the lower arc meets the ground where ground(x) = zc − √(R² − (x − xc)²); scan, then refine by bisection
const gx = geom.ground, xa = Math.max(gx[0].x, c.xc - c.R), xb = Math.min(gx[gx.length - 1].x, c.xc + c.R), n = 400;
const f = (x) => interp(gx, x) - (c.zc - Math.sqrt(Math.max(0, c.R * c.R - (x - c.xc) ** 2))); // > 0 where the ground is above the arc (inside)
const roots = []; let px = xa, pf = f(xa);
for (let i = 1; i <= n; i++) {
const x = xa + ((xb - xa) * i) / n, v = f(x);
if ((pf > 0) !== (v > 0)) { let lo = px, hi = x; for (let k = 0; k < 40; k++) { const m = (lo + hi) / 2; if ((f(lo) > 0) === (f(m) > 0)) lo = m; else hi = m; } roots.push({ x: (lo + hi) / 2, into: v > 0 }); }
px = x; pf = v;
}
return { roots, inside0: f(xa) > 0, inside1: f(xb) > 0, xa, xb };
}Each slice's width, weight, base angle, strength and pore pressure.
W = Σ layers γ × thickness × b; sin α = −(x − xc) ÷ R; u = γw × (zw − zb), or ru × W ÷ b
Slices are of equal arc angle, so they narrow where the arc is steep. Weight uses the layer thicknesses at the slice's middle times its width, a rectangle rule. cos α is √(1 − sin²α), clamped at zero. The slice count is made whole with a bitwise OR; it is a count, not a measured value.
const n = Math.max(10, Math.min(400, slices | 0)), S = [];
const th0 = Math.asin(clamp((xl - c.xc) / c.R, -1, 1)), th1 = Math.asin(clamp((xr - c.xc) / c.R, -1, 1)), xAt = (i) => c.xc + c.R * Math.sin(th0 + ((th1 - th0) * i) / n);
const layerAt = (x, z) => {
for (let i = geom.layers.length - 1; i >= 0; i--) if (z <= interp(geom.layers[i].top, x) + 1e-12 || i === 0) return geom.layers[i];
return geom.layers[0];
};
const tops = (x) => geom.layers.map((l) => Math.min(interp(l.top, x), interp(geom.ground, x)));
for (let i = 0; i < n; i++) {
const xa = xAt(i), b = xAt(i + 1) - xa, x = xa + b / 2, zg = interp(geom.ground, x), dxc = x - c.xc, h = c.R * c.R - dxc * dxc;
if (h < 0) continue;
const zb = c.zc - Math.sqrt(h);
if (zb >= zg - 1e-9) continue;
if (zb < geom.floor - 1e-9) problems.push({ code: 'geo.slope.floor' });
// weight: the thickness of each layer between the base and the ground
const t = tops(x); let W = 0, zsum = 0;
for (let k = 0; k < geom.layers.length; k++) {
const hi = Math.min(t[k], zg), lo = k + 1 < geom.layers.length ? Math.min(t[k + 1], zg) : -Infinity, a = Math.max(lo, zb), top = hi;
if (top > a) { const w = geom.layers[k].gamma * (top - a) * b; W += w; zsum += w * (top + a) / 2; }
}
let q = 0; for (const s of geom.surcharge || []) if (x >= s.x0 && x <= s.x1) q += s.q * b;
const Wt = W + q, zcg = Wt > 0 ? (zsum + q * zg) / Wt : zb;
const lay = layerAt(x, zb + 1e-9);
const cc = lay.cu !== undefined ? lay.cu : lay.c, ph = lay.cu !== undefined ? 0 : lay.phi;
const sinA = -dxc / c.R, cosA = Math.sqrt(Math.max(0, 1 - sinA * sinA));
let u = 0;
if (geom.water !== undefined && geom.water !== null) { const zw = interp(geom.water, x); if (zw > zb) u = G_W * (zw - zb); }
else if (geom.ru) u = geom.ru * (W / b);
S.push({ x, b, zb, zg, W: Wt, sinA, cosA, c: cc, phi: ph, u, zcg });The factor of safety of one circle by the ordinary method or by Bishop's simplified method.
Bishop: FS = Σ[(c·b + (W − u·b)·tan φ) ÷ (cos α + sin α·tan φ ÷ FS)] ÷ Σ W·sin α
Bishop's equation is iterated from 1.5 until the change is below 1e-5, an absolute tolerance, or for 200 passes; there is no warning if the limit is reached. The effective weight term is clamped at zero and mα is floored at 1e-6. The ordinary method divides by cos α with no guard.
let drive = 0; for (const s of S) { drive += s.W * s.sinA; if (kh) drive += (kh * s.W * (c.zc - s.zcg)) / c.R; }
if (drive <= 1e-9) return { fs: Infinity, ok: true, slices: S, moment: drive, problems: [{ code: 'geo.slope.nodrive' }], xl, xr };
let fs;
if (method === 'ordinary') {
let res = 0; for (const s of S) res += (s.c * s.b) / s.cosA + Math.max(0, s.W * s.cosA - (s.u * s.b) / s.cosA) * Math.tan(s.phi * RAD);
fs = res / drive;
} else {
fs = 1.5;
for (let it = 0; it < 200; it++) {
let res = 0, bad = false;
for (const s of S) {
const m = s.cosA + (s.sinA * Math.tan(s.phi * RAD)) / fs;
if (m < 0.2) bad = true;
res += (s.c * s.b + Math.max(0, s.W - s.u * s.b) * Math.tan(s.phi * RAD)) / Math.max(m, 1e-6);
}
const next = res / drive;
if (bad && !problems.some((p) => p.code === 'geo.slope.ma')) problems.push({ code: 'geo.slope.ma' });
if (Math.abs(next - fs) < tol) { fs = next; break; }
fs = next;
}
}
return { fs, ok: fin(fs), slices: S, moment: drive, problems, xl, xr, mirrored: !!geom.mirrored };The circle with the lowest factor of safety among those tried.
A grid of centers (15 by 11 by default) with 14 radii each plus circles through the ground's corner points; then four rounds on a 9 by 9 grid round the best, and radius steps of 1 % up to ±6 %. The final value is recomputed with at least 80 slices. The result is the lowest sampled value, not a proven minimum.
/**
* Search the critical circle: a grid of centres above the ground, for each a fan of radii (every circle that reaches the ground and stays above
* `floor`), keeping the smallest FS; then the grid is refined around the best three times.
* opts: { x0, x1 (centre range), z0, z1, nx, nz, nr, method, slices, kh, entry: [xa, xb] (circle must enter/exit within), minDepth }
*/
export function searchSlip(geom0, opts = {}) {
const geom = faceRight(geom0) ? geom0 : mirrorGeom(geom0), g = geom.ground, xs = g.map((p) => p.x), zs = g.map((p) => p.z);
const xmin = Math.min(...xs), xmax = Math.max(...xs), zmax = Math.max(...zs), zmin = Math.min(...zs), H = zmax - zmin || 1;
const floor = geom.floor ?? zmin - H;
let bx0 = opts.x0 ?? xmin, bx1 = opts.x1 ?? xmax, bz0 = opts.z0 ?? zmax, bz1 = opts.z1 ?? zmax + 2.5 * H;
if (geom !== geom0 && opts.x0 !== undefined) { bx0 = -opts.x1; bx1 = -opts.x0; }
const nx = opts.nx ?? 14, nz = opts.nz ?? 10, nr = opts.nr ?? 14, method = opts.method || 'bishop', slices = opts.slices ?? 40;
const grid = []; let best = null, tried = 0;
const gg = { ...geom, floor };
const scan = (x0, x1, z0, z1, kx, kz, record) => {
for (let i = 0; i <= kx; i++) for (let j = 0; j <= kz; j++) {
const xc = x0 + ((x1 - x0) * i) / kx, zc = z0 + ((z1 - z0) * j) / kz;
let cellBest = null;
// radii from the one that just reaches the floor down to the smallest that cuts the ground
const rMax = zc - floor, rMin = Math.max(zc - zmax, 0.05 * H) + 1e-3;
if (rMax <= rMin) continue;
const radii = [];
for (let k = 0; k < nr; k++) radii.push(rMin + ((rMax - rMin) * (k + 0.5)) / nr);
for (const p of g) { const d = Math.hypot(p.x - xc, p.z - zc); for (const f of [0.999, 1.001]) if (d * f > rMin && d * f < rMax) radii.push(d * f); } // circles through the face's corners (the toe)
radii.push(rMax * 0.999);
for (const R of radii) {
tried++;
const circle = { xc, zc, R }, r = bishop(gg, circle, { method, slices, kh: opts.kh || 0 });
if (!r.ok || !fin(r.fs) || r.fs <= 0) continue;
if (opts.entry && (r.xl < opts.entry[0] || r.xr > opts.entry[1])) continue;
if (r.problems.some((p) => p.code === 'geo.slope.floor')) continue;
if (!cellBest || r.fs < cellBest.fs) cellBest = { xc, zc, R, fs: r.fs, xl: r.xl, xr: r.xr };
}
if (cellBest) { if (record) grid.push({ xc, zc, fs: cellBest.fs }); if (!best || cellBest.fs < best.fs) best = cellBest; }
}
};
scan(bx0, bx1, bz0, bz1, nx, nz, true);
for (let round = 0; round < 4 && best; round++) {
const sx = (bx1 - bx0) / nx, sz = (bz1 - bz0) / nz, cx = best.xc, cz = best.zc;
bx0 = cx - sx; bx1 = cx + sx; bz0 = Math.max(cz - sz, zmax + 0.01 * H); bz1 = cz + sz;
scan(bx0, bx1, bz0, bz1, 8, 8, false);
// refine the radius around the best too
for (let k = -6; k <= 6; k++) {
const R = best.R + k * 0.01 * best.R, circle = { xc: best.xc, zc: best.zc, R }, r = bishop(gg, circle, { method, slices, kh: opts.kh || 0 }); tried++;
if (r.ok && fin(r.fs) && r.fs > 0 && !r.problems.some((p) => p.code === 'geo.slope.floor') && r.fs < best.fs) best = { xc: best.xc, zc: best.zc, R, fs: r.fs, xl: r.xl, xr: r.xr };
}
}
if (!best) return { critical: null, grid, tried, problems: [{ code: 'geo.slope.nocircle' }] };
const final = bishop(gg, best, { method, slices: Math.max(slices, 80), kh: opts.kh || 0 });
const back = (cc) => (geom === geom0 ? cc : { ...cc, xc: -cc.xc, xl: -cc.xr, xr: -cc.xl });
return { critical: { circle: { xc: back(best).xc, zc: best.zc, R: best.R }, fs: final.fs, xl: back(best).xl, xr: back(best).xr, problems: final.problems }, grid: grid.map((q) => (geom === geom0 ? q : { ...q, xc: -q.xc })), tried, floor };
}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.