IFC openings, clippings and half-space cuts

/** Outward face planes [{ n, d }] of a convex closed mesh (n · p ≤ d inside), or null when it is not convex. */
export function convexPlanes(mesh, tol = 1e-6) {
  const P = mesh.pos, I = mesh.idx;
  const b = mesh.bounds();
  if (!b) return null;
  const size = Math.max(1e-9, len(sub(b.max, b.min)));
  const eps = size * tol;
  const planes = [];
  for (let t = 0; t < I.length; t += 3) {
    const a = [P[I[t] * 3], P[I[t] * 3 + 1], P[I[t] * 3 + 2]], bb = [P[I[t + 1] * 3], P[I[t + 1] * 3 + 1], P[I[t + 1] * 3 + 2]], c = [P[I[t + 2] * 3], P[I[t + 2] * 3 + 1], P[I[t + 2] * 3 + 2]];
    const cr = cross(sub(bb, a), sub(c, a));
    if (len(cr) < eps * eps) continue;
    const n = norm(cr), d = dot(n, a);
    if (planes.some((q) => dot(q.n, n) > 1 - 1e-9 && Math.abs(q.d - d) < eps)) continue;
    planes.push({ n, d });
  }
  if (planes.length < 4) return null;
  for (let i = 0; i < P.length; i += 3) for (const q of planes) if (q.n[0] * P[i] + q.n[1] * P[i + 1] + q.n[2] * P[i + 2] - q.d > eps * 10) return null;
  return planes;
}

/** Split a convex polygon by plane q: [part with n·p ≥ d, part with n·p < d]. */
function splitPoly(poly, q) {
  const front = [], back = [];
  for (let i = 0; i < poly.length; i++) {
    const a = poly[i], b = poly[(i + 1) % poly.length];
    const da = dot(q.n, a) - q.d, db = dot(q.n, b) - q.d;
    if (da >= 0) front.push(a); else back.push(a);
    if ((da >= 0) !== (db >= 0)) {
      const t = da / (da - db);
      const x = [a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t, a[2] + (b[2] - a[2]) * t];
      front.push(x); back.push(x);
    }
  }
  return [front, back];
}

/**
 * Keep the part of `mesh` on the side of the plane (point p, normal n) that n points to. Triangles are split exactly;
 * cut edges are chained into loops and capped (outer loops and the holes inside them) when they close.
 */
export function clipPlane(mesh, p, n) {
  n = norm(n);
  const out = new Mesh();
  const P = mesh.pos, N = mesh.nrm, I = mesh.idx;
  const sd = new Float64Array(P.length / 3);
  let lo = Infinity, hi = -Infinity;
  for (let i = 0; i < sd.length; i++) { const d = (P[i * 3] - p[0]) * n[0] + (P[i * 3 + 1] - p[1]) * n[1] + (P[i * 3 + 2] - p[2]) * n[2]; sd[i] = d; if (d < lo) lo = d; if (d > hi) hi = d; }
  const span = Math.max(1e-9, hi - lo);
  const eps = span * 1e-9;
  if (lo >= -eps) return out.append(mesh);
  if (hi <= eps) return out;
  const segs = [];
  const v = (i) => [P[i * 3], P[i * 3 + 1], P[i * 3 + 2]];
  const nv = (i) => [N[i * 3], N[i * 3 + 1], N[i * 3 + 2]];
  const lerp = (i, j) => { const t = sd[i] / (sd[i] - sd[j]); const a = v(i), b = v(j); return [a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t, a[2] + (b[2] - a[2]) * t]; };
  for (let t = 0; t < I.length; t += 3) {
    const ids = [I[t], I[t + 1], I[t + 2]];
    const inside = ids.map((i) => sd[i] >= -eps);
    const k = inside.filter(Boolean).length;
    const fn = nv(ids[0]);
    const emit = (pts) => { const base = out.vertexCount; for (const q of pts) out.vertex(q, fn); for (let s = 1; s + 1 < pts.length; s++) out.idx.push(base, base + s, base + s + 1); };
    if (k === 3) { emit(ids.map(v)); continue; }
    if (k === 0) continue;
    // walk the triangle's edges, keeping the inside part
    const poly = [], cut = [];
    for (let e = 0; e < 3; e++) {
      const i = ids[e], j = ids[(e + 1) % 3];
      const ii = sd[i] >= -eps, jj = sd[j] >= -eps;
      if (ii) poly.push(v(i));
      if (ii !== jj && Math.abs(sd[i]) > eps && Math.abs(sd[j]) > eps) { const q = lerp(i, j); poly.push(q); cut.push(q); }
      else if (ii && !jj && Math.abs(sd[i]) <= eps) cut.push(v(i));
      else if (!ii && jj && Math.abs(sd[j]) <= eps) cut.push(v(j));
    }
    if (poly.length >= 3) emit(poly);
    if (cut.length === 2) segs.push(cut);
  }
  capLoops(out, segs, scale(n, -1), eps * 10 + span * 1e-7);
  return out;
}

/** Chain cut segments into closed loops on the plane and fill them (normal `cn`, facing out of the kept part). */
function capLoops(out, segs, cn, tol) {
  if (!segs.length) return;
  const key = (q) => Math.round(q[0] / tol) + ',' + Math.round(q[1] / tol) + ',' + Math.round(q[2] / tol);
  const adj = new Map();
  for (const [a, b] of segs) {
    const ka = key(a), kb = key(b);
    if (ka === kb) continue;
    if (!adj.has(ka)) adj.set(ka, { p: a, to: [] });
    if (!adj.has(kb)) adj.set(kb, { p: b, to: [] });
    adj.get(ka).to.push(kb); adj.get(kb).to.push(ka);
  }
  const used = new Set(), loops = [];
  const ek = (a, b) => (a < b ? a + '|' + b : b + '|' + a);
  for (const [start] of adj) {
    for (const first of adj.get(start).to) {
      if (used.has(ek(start, first))) continue;
      const loop = [adj.get(start).p];
      used.add(ek(start, first));
      let prev = start, cur = first, guard = 0;
      while (cur !== start && guard++ < 1e6) {
        loop.push(adj.get(cur).p);
        const nxt = adj.get(cur).to.find((x) => x !== prev && !used.has(ek(cur, x)));
        if (!nxt) break;
        used.add(ek(cur, nxt));
        prev = cur; cur = nxt;
      }
      if (cur === start && loop.length >= 3) loops.push(loop);
    }
  }
  if (!loops.length) return;
  // 2D frame on the plane; loops CCW about cn are outers, nesting decides holes
  const u = perp(cn), w = cross(cn, u);
  const to2 = (q) => [dot(q, u), dot(q, w)];
  const L = loops.map((l) => { const p2 = l.map(to2); let a = 0; for (let i = 0; i < p2.length; i++) { const x = p2[i], y = p2[(i + 1) % p2.length]; a += x[0] * y[1] - y[0] * x[1]; } return { l, p2, area: Math.abs(a / 2), ccw: a > 0 }; });
  L.sort((x, y) => y.area - x.area);
  const inside = (pt, poly) => { let c = false; for (let i = 0, j = poly.length - 1; i < poly.length; j = i++) { const a = poly[i], b = poly[j]; if ((a[1] > pt[1]) !== (b[1] > pt[1]) && pt[0] < (b[0] - a[0]) * (pt[1] - a[1]) / (b[1] - a[1]) + a[0]) c = !c; } return c; };
  const depth = L.map((x, i) => L.slice(0, i).filter((y) => inside(x.p2[0], y.p2)).length);
  for (let i = 0; i < L.length; i++) {
    if (depth[i] % 2) continue;
    const outer = L[i].ccw ? L[i].l : L[i].l.slice().reverse();
    const holes = [];
    for (let j = i + 1; j < L.length; j++) if (depth[j] === depth[i] + 1 && inside(L[j].p2[0], L[i].p2)) holes.push(L[j].ccw ? L[j].l.slice().reverse() : L[j].l);
    const rings = [outer, ...holes];
    const tris = planeTriangulate(rings);
    const base = out.vertexCount;
    for (const r of rings) for (const q of r) out.vertex(q, cn);
    for (let k = 0; k < tris.length; k += 3) out.idx.push(base + tris[k], base + tris[k + 1], base + tris[k + 2]);
  }
}