diff --git a/Cargo.lock b/Cargo.lock
index 820715c..4d5a32f 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -1907,6 +1907,7 @@ name = "modelcraft-testkit"
 version = "0.1.0"
 dependencies = [
  "modelcraft-geometry",
+ "modelcraft-mesh",
  "serde_json",
  "thiserror 2.0.21",
 ]
diff --git a/ROADMAP.md b/ROADMAP.md
index 8eed8e9..319de9a 100644
--- a/ROADMAP.md
+++ b/ROADMAP.md
@@ -17,7 +17,7 @@ hex meshes today and it gives the best elements where stresses are read.
 | Milestone | Scope | Status |
 |---|---|---|
 | M0 Foundation | Workspace, quality gates, FE model, command engine with rollback and panic guard, CLI, geometry boundary on SolidCraft | done |
-| M1 Geometry in and surface preparation | Open `.solidcraft` designs with persistent face names, assemblies placed (done); face-tagged 3MF (done); test corpus, `cargo xtask corpus` (done); exact predicates (done); quality metrics and certified validity gate (done); weld, closed/manifold/orientation checks, feature edges | in progress |
+| M1 Geometry in and surface preparation | Open `.solidcraft` designs with persistent face names, assemblies placed (done); face-tagged 3MF (done); test corpus, `cargo xtask corpus` (done); exact predicates (done); quality metrics and certified validity gate (done); weld, closed/manifold/orientation and self-intersection checks, feature edges (`modelcraft-mesh::surface`, `mesh.surface.check`, done) | in progress |
 | M2 Tet meshing and size field | Size field (global, per face and edge seeds, curvature and proximity); Bowyer–Watson Delaunay (exact, symbolic perturbation, verified; `modelcraft-mesh::delaunay`, done); robust CDT boundary recovery (Diazzi et al. 2023); Delaunay refinement, ODT smoothing, sliver removal; tet4 and tet10 with mid-side nodes on the B-rep; mesh verify per set | in progress |
 | M3 First deck | Sets and surfaces bound to SolidCraft face, edge, vertex and body names, resolved after meshing; materials (isotropic and orthotropic elastic, density) with a unit-system label and consistency check; sections (solid, shell, beam profiles) assigned by set; element formulation per region mapped to each solver; BCs (fixed, displacement, symmetry); loads (force and moment on reference points, pressure, gravity); reference points with kinematic and distributing coupling; steps (linear static, frequency) with explicit per-step load and BC tables and output presets; writers for Abaqus input, Nastran SOL 101/103 and LS-DYNA implicit; the `.modelcraft` file | planned |
 | M4 Own solver and first results | Linear static and modal solver (tet4/10, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions; deformed shape, contours, mode shapes, probes, reaction sums | planned |
diff --git a/crates/engine/src/cmd/geometry.rs b/crates/engine/src/cmd/geometry.rs
index 81d4224..b78da87 100644
--- a/crates/engine/src/cmd/geometry.rs
+++ b/crates/engine/src/cmd/geometry.rs
@@ -1,7 +1,7 @@
 use std::sync::Arc;
 
 use modelcraft_geometry::tagged_mesh;
-use modelcraft_geometry::{Shape, open_design_file};
+use modelcraft_geometry::{Shape, SurfaceMesh, open_design_file};
 use serde_json::{Value, json};
 
 use crate::{CommandSpec, EngineError, Result, Session};
@@ -72,10 +72,15 @@ fn body<'a>(s: &'a Session, p: &Value) -> Result<&'a Shape> {
     s.design.as_deref().and_then(|d| d.body(name)).ok_or_else(|| EngineError::Params(format!("no body `{name}` (open a design first)")))
 }
 
-fn tessellate(s: &mut Session, p: &Value) -> Result<Value> {
+/// The body named by `body`, tessellated within `tolerance`.
+pub(crate) fn body_mesh<'a>(s: &'a Session, p: &Value) -> Result<(&'a Shape, SurfaceMesh)> {
     let shape = body(s, p)?;
     let tolerance = p.get("tolerance").and_then(Value::as_f64).ok_or_else(|| EngineError::Params("`tolerance` must be a number (mm)".into()))?;
-    let mesh = shape.tessellate(tolerance)?;
+    Ok((shape, shape.tessellate(tolerance)?))
+}
+
+fn tessellate(s: &mut Session, p: &Value) -> Result<Value> {
+    let (shape, mesh) = body_mesh(s, p)?;
     let mut per_face = vec![0usize; shape.faces().len()];
     for &f in &mesh.tri_face {
         if let Some(n) = per_face.get_mut(f as usize) {
diff --git a/crates/engine/src/cmd/mesh.rs b/crates/engine/src/cmd/mesh.rs
index e439388..9265612 100644
--- a/crates/engine/src/cmd/mesh.rs
+++ b/crates/engine/src/cmd/mesh.rs
@@ -1,9 +1,13 @@
+use std::borrow::Cow;
+
+use modelcraft_geometry::SurfaceMesh;
 use modelcraft_mesh::delaunay::Triangulation;
 use modelcraft_mesh::quality::{self, QualityReport};
+use modelcraft_mesh::surface::{BoundarySurface, Defect, FeatureKind, SurfaceOptions};
 use modelcraft_model::{ElementKind, MAX_ELEMENTS};
 use serde_json::{Value, json};
 
-use super::{label, labels, points, xyz};
+use super::{label, labels, optional_number, points, xyz};
 use crate::{CommandSpec, EngineError, Result, Session};
 
 pub const NODE_ADD: CommandSpec = CommandSpec {
@@ -57,6 +61,13 @@ fn delaunay_points(_: &mut Session, p: &Value) -> Result<Value> {
     }))
 }
 
+pub const SURFACE_CHECK: CommandSpec = CommandSpec {
+    id: "mesh.surface.check",
+    label: "Check Boundary Surface",
+    params: r#"{"body": "Body1", "tolerance": 0.05} or {"surface": 0}, optional "weld": 1e-9 (0 or 1e-15..=1e-3 of the bounding-box diagonal), "weld_distance": mm (default: the surface's own, twice SolidCraft's linear tolerance for a tessellated body, 0 for an imported mesh), "feature_angle": 45 (degrees, (0, 180]) — welds the tessellated body or imported surface and checks it can bound a volume mesh: closed, manifold, outward, no self-intersections; reports shells, genus, volume and feature edges, or the located defects"#,
+    run: surface_check,
+};
+
 fn node_add(s: &mut Session, p: &Value) -> Result<Value> {
     let l = label(p, "label")?;
     s.model.set_node(l, xyz(p, "xyz")?)?;
@@ -136,6 +147,63 @@ fn quality_report(s: &mut Session, p: &Value) -> Result<Value> {
     }))
 }
 
+/// The surface named by `body` (tessellated at `tolerance`) or `surface` (an imported index).
+fn surface_input<'a>(s: &'a Session, p: &Value) -> Result<Cow<'a, SurfaceMesh>> {
+    match (p.get("body"), p.get("surface")) {
+        (Some(_), None) => Ok(Cow::Owned(super::geometry::body_mesh(s, p)?.1)),
+        (None, Some(i)) => {
+            let n = s.surfaces.len();
+            i.as_u64()
+                .and_then(|i| usize::try_from(i).ok())
+                .and_then(|i| s.surfaces.get(i))
+                .map(|t| Cow::Borrowed(&t.mesh))
+                .ok_or_else(|| EngineError::Params(format!("`surface` must be the index of an imported surface (there are {n})")))
+        }
+        _ => Err(EngineError::Params("give either `body` (with `tolerance`) or `surface`".into())),
+    }
+}
+
+fn defect_json(d: &Defect) -> Value {
+    json!({"kind": d.kind(), "message": d.to_string(), "at": d.at()})
+}
+
+fn surface_check(s: &mut Session, p: &Value) -> Result<Value> {
+    let mesh = surface_input(s, p)?;
+    let defaults = SurfaceOptions::default();
+    let options = SurfaceOptions {
+        weld_tolerance: optional_number(p, "weld", defaults.weld_tolerance)?,
+        weld_distance: optional_number(p, "weld_distance", mesh.weld_distance)?,
+        feature_angle: optional_number(p, "feature_angle", defaults.feature_angle)?,
+    };
+    let surface = match BoundarySurface::new(&mesh.positions, &mesh.triangles, &mesh.tri_face, &options) {
+        Ok(surface) => surface,
+        Err(modelcraft_mesh::Error::Surface(d)) => {
+            return Ok(json!({"valid": false, "defects": {"count": d.total, "listed": d.listed.iter().map(defect_json).collect::<Vec<_>>()}}));
+        }
+        Err(e) => return Err(e.into()),
+    };
+    let mut faces: Vec<u32> = surface.tri_face().to_vec();
+    faces.sort_unstable();
+    faces.dedup();
+    let count = |k: FeatureKind| surface.feature_edges().iter().filter(|e| e.kind == k).count();
+    let shells: Vec<Value> = surface
+        .shells()
+        .iter()
+        .map(|sh| json!({"triangles": sh.triangles, "vertices": sh.vertices, "genus": sh.genus, "volume": sh.volume, "void": sh.void}))
+        .collect();
+    Ok(json!({
+        "valid": true,
+        "vertices": surface.points().len(),
+        "triangles": surface.triangles().len(),
+        "input": {"vertices": mesh.positions.len(), "triangles": mesh.triangles.len()},
+        "faces": faces.len(),
+        "volume": surface.volume(),
+        "shells": shells,
+        "feature_edges": {"face_boundary": count(FeatureKind::FaceBoundary), "sharp": count(FeatureKind::Sharp)},
+        "feature_corners": surface.feature_corners().len(),
+    }))
+}
+
 #[cfg(test)]
 mod tests {
     use serde_json::json;
@@ -202,4 +270,81 @@ mod tests {
         assert!(s.execute("mesh.quality", &json!({"bins": 1e9})).is_err());
         assert!(s.execute("mesh.quality", &json!({"elements": "all"})).is_err());
     }
+
+    #[test]
+    fn a_tessellated_block_is_a_valid_boundary() {
+        use modelcraft_mesh::surface::{BoundarySurface, FeatureKind, SurfaceOptions};
+        let mesh = modelcraft_geometry::block([0.0; 3], [2.0, 3.0, 4.0]).unwrap().tessellate(0.01).unwrap();
+        let options = SurfaceOptions { weld_distance: mesh.weld_distance, ..SurfaceOptions::default() };
+        let s = BoundarySurface::new(&mesh.positions, &mesh.triangles, &mesh.tri_face, &options).unwrap();
+        assert!(s.points().len() < mesh.positions.len(), "the tessellation is a soup that welding joins");
+        assert_eq!((s.shells().len(), s.shells()[0].genus), (1, 0));
+        assert!((s.volume() - 24.0).abs() < 1e-9, "{}", s.volume());
+        assert!(s.feature_edges().iter().all(|e| e.kind == FeatureKind::FaceBoundary));
+        assert_eq!(s.feature_corners().len(), 8);
+        let length: f64 = s
+            .feature_edges()
+            .iter()
+            .map(|e| {
+                let [a, b] = e.vertices.map(|v| s.points()[v as usize]);
+                ((a[0] - b[0]).powi(2) + (a[1] - b[1]).powi(2) + (a[2] - b[2]).powi(2)).sqrt()
+            })
+            .sum();
+        assert!((length - 4.0 * (2.0 + 3.0 + 4.0)).abs() < 1e-9, "the twelve edges: {length}");
+    }
+
+    #[test]
+    fn surface_check_on_an_open_design_and_on_imported_surfaces() {
+        use std::sync::Arc;
+
+        use modelcraft_geometry::sample::{Plate, plate_design};
+        use modelcraft_geometry::tagged_mesh::TaggedSurface;
+        use modelcraft_geometry::{SurfaceMesh, open_design};
+
+        let plate = Plate::default();
+        let mut s = Session { design: Some(Arc::new(open_design(&plate_design(&plate).unwrap()).unwrap())), ..Session::default() };
+        let r = s.execute("mesh.surface.check", &json!({"body": "Body1", "tolerance": 0.05})).unwrap();
+        assert_eq!(r["valid"], true, "{r}");
+        assert_eq!(r["faces"], 8);
+        assert_eq!(r["shells"].as_array().unwrap().len(), 1);
+        assert_eq!(r["shells"][0]["genus"], 1, "the hole is a handle");
+        let exact = (plate.width * plate.depth - std::f64::consts::PI * (plate.hole_diameter / 2.0).powi(2)) * plate.thickness;
+        let volume = r["volume"].as_f64().unwrap();
+        assert!(volume >= exact && volume < exact * 1.01, "a polygonal hole inscribed in the circle: {volume} vs {exact}");
+        assert!(r["feature_edges"]["face_boundary"].as_u64().unwrap() > 12);
+        assert!(r["feature_corners"].as_u64().unwrap() >= 8);
+        assert!(s.journal.is_empty(), "a check changes nothing");
+
+        // A tagged tetrahedron, and the same with a face missing.
+        let p = vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]];
+        let t = vec![[0, 2, 1], [0, 1, 3], [1, 2, 3], [0, 3, 2]];
+        let tet = SurfaceMesh { positions: p.clone(), triangles: t.clone(), tri_face: vec![0, 1, 1, 1], weld_distance: 0.0 };
+        let open = SurfaceMesh { positions: p, triangles: t[..3].to_vec(), tri_face: vec![0, 1, 1], weld_distance: 0.0 };
+        for mesh in [tet, open] {
+            s.surfaces.push(Arc::new(TaggedSurface { name: "Tet".into(), mesh, faces: 2 }));
+        }
+        let r = s.execute("mesh.surface.check", &json!({"surface": 0, "feature_angle": 180})).unwrap();
+        assert_eq!((r["valid"].clone(), r["feature_edges"].clone()), (json!(true), json!({"face_boundary": 3, "sharp": 0})));
+        let r = s.execute("mesh.surface.check", &json!({"surface": 1})).unwrap();
+        assert_eq!(r["valid"], false);
+        assert_eq!(r["defects"]["count"], 3);
+        assert_eq!(r["defects"]["listed"][0]["kind"], "open_edge");
+        assert!(r["defects"]["listed"][0]["at"].is_array());
+
+        for bad in [
+            json!({"surface": 2}),
+            json!({"surface": -1}),
+            json!({"surface": 0, "body": "Body1"}),
+            json!({}),
+            json!({"surface": 0, "weld": 1.0}),
+            json!({"surface": 0, "weld": "fine"}),
+            json!({"surface": 0, "weld_distance": -1}),
+            json!({"surface": 0, "weld_distance": 1}),
+            json!({"surface": 0, "feature_angle": 0}),
+            json!({"body": "Body1"}),
+            json!({"body": "Nope", "tolerance": 0.05}),
+        ] {
+            assert!(s.execute("mesh.surface.check", &bad).is_err(), "{bad}");
+        }
+    }
 }
diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs
index d1936c6..51ca70d 100644
--- a/crates/engine/src/cmd/mod.rs
+++ b/crates/engine/src/cmd/mod.rs
@@ -21,6 +21,7 @@ pub static ALL: &[CommandSpec] = &[
     mesh::ELEMENT_REMOVE,
     mesh::QUALITY,
     mesh::DELAUNAY_POINTS,
+    mesh::SURFACE_CHECK,
 ];
 
 /// A required label parameter (1..=u32::MAX; the model checks its own range).
@@ -63,3 +64,11 @@ pub(crate) fn labels(p: &Value, key: &str, max_len: usize) -> Result<Vec<u32>> {
         .collect::<Option<Vec<u32>>>()
         .ok_or_else(|| EngineError::Params(format!("`{key}` must be labels in 1..={MAX_LABEL}")))
 }
+
+/// An optional number parameter.
+pub(crate) fn optional_number(p: &Value, key: &str, default: f64) -> Result<f64> {
+    match p.get(key) {
+        None => Ok(default),
+        Some(v) => v.as_f64().ok_or_else(|| EngineError::Params(format!("`{key}` must be a number"))),
+    }
+}
diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs
index 1122928..8b5d2e6 100644
--- a/crates/engine/src/lib.rs
+++ b/crates/engine/src/lib.rs
@@ -184,6 +184,10 @@ mod tests {
             json!({"points": [[1e308, -1e308, 0], [0, 0, 0]]}),
             json!({"points": [[0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0]]}),
             json!({"points": vec![[1, 1, 1]; 30_000]}),
+            json!({"surface": 0}),
+            json!({"surface": -1}),
+            json!({"surface": 1e300}),
+            json!({"surface": "0", "weld": -1, "feature_angle": 1e300}),
         ];
         for spec in commands() {
             for p in &hostile {
@@ -204,6 +208,15 @@ mod tests {
             json!({"body": "Body1", "tolerance": 1e-300}),
             json!({"body": "Body1", "tolerance": 1e300}),
             json!({"body": "Body2", "tolerance": 0.1}),
+            json!({"body": "Body1", "tolerance": 0.5, "surface": 0}),
+            json!({"body": "Body1", "tolerance": 0.5, "weld": 0}),
+            json!({"body": "Body1", "tolerance": 0.5, "weld": 1e-3}),
+            json!({"body": "Body1", "tolerance": 0.5, "weld": f64::MIN_POSITIVE}),
+            json!({"body": "Body1", "tolerance": 0.5, "weld": 0, "weld_distance": 0}),
+            json!({"body": "Body1", "tolerance": 0.5, "weld_distance": 1e300}),
+            json!({"body": "Body1", "tolerance": 0.5, "weld_distance": 0.05}),
+            json!({"body": "Body1", "tolerance": 0.5, "feature_angle": 1e-300}),
+            json!({"body": "Body1", "tolerance": 0.5, "feature_angle": null}),
         ];
         for p in tess.iter().chain(&hostile) {
             for spec in commands() {
diff --git a/crates/geometry/src/lib.rs b/crates/geometry/src/lib.rs
index 99d52cd..4f4cf3c 100644
--- a/crates/geometry/src/lib.rs
+++ b/crates/geometry/src/lib.rs
@@ -47,6 +47,9 @@ impl From<KernelError> for Error {
     }
 }
 
+/// SolidCraft's linear modelling tolerance (mm): points closer than this are the same point.
+const LINEAR_TOLERANCE: f64 = solidcraft_geom::LINEAR_TOL;
+
 /// A solid body.
 #[derive(Debug, Clone)]
 pub struct Shape {
@@ -86,6 +89,11 @@ pub struct SurfaceMesh {
     pub triangles: Vec<[u32; 3]>,
     /// The B-rep face each triangle lies on.
     pub tri_face: Vec<u32>,
+    /// Copies of one point closer than this are the same point (0: only identical coordinates).
+    /// A [`Shape::tessellate`] result is a triangle soup whose faces each carry their own copy of
+    /// the boundary points they share, accurate to the kernel's linear tolerance, so two copies
+    /// can differ by up to twice it; boundary surface preparation welds at this distance.
+    pub weld_distance: f64,
 }
 
 fn xyz(v: solidcraft_geom::Vec3) -> [f64; 3] {
@@ -129,7 +137,12 @@ impl Shape {
             return Err(Error::Invalid(format!("tessellation tolerance must be a positive number, got {tolerance}")));
         }
         let m = guard("tessellate", || self.body.tessellate(tolerance))?;
-        Ok(SurfaceMesh { positions: m.positions.into_iter().map(xyz).collect(), triangles: m.triangles, tri_face: m.tri_face })
+        Ok(SurfaceMesh {
+            positions: m.positions.into_iter().map(xyz).collect(),
+            triangles: m.triangles,
+            tri_face: m.tri_face,
+            weld_distance: 2.0 * LINEAR_TOLERANCE,
+        })
     }
 }
 
diff --git a/crates/geometry/src/tagged_mesh.rs b/crates/geometry/src/tagged_mesh.rs
index e238a3f..2c91b08 100644
--- a/crates/geometry/src/tagged_mesh.rs
+++ b/crates/geometry/src/tagged_mesh.rs
@@ -341,7 +341,7 @@ fn surface(mut o: Object, k: usize, id: u32, xf: &Xf, unit: f64) -> Result<Tagge
     }
     let faces = o.faces.unwrap_or_else(|| o.tags.iter().max().map_or(1, |f| f.saturating_add(1)));
     let name = if o.name.is_empty() { format!("Mesh{}", k + 1) } else { o.name };
-    Ok(TaggedSurface { name, mesh: SurfaceMesh { positions: o.positions, triangles: o.triangles, tri_face: o.tags }, faces })
+    Ok(TaggedSurface { name, mesh: SurfaceMesh { positions: o.positions, triangles: o.triangles, tri_face: o.tags, weld_distance: 0.0 }, faces })
 }
 
 /// [`read_3mf`] of a file: a regular file of at most [`MAX_FILE_BYTES`] (SolidCraft's virtual
@@ -482,6 +482,8 @@ mod tests {
         assert_eq!(bytes, write_tagged(&[("Block", &mesh, faces), ("Again", &mesh, faces)]), "deterministic bytes");
         let got = read_3mf(&bytes).unwrap();
         assert_eq!(got.len(), 2);
+        // The file carries no weld distance.
+        let mesh = SurfaceMesh { weld_distance: 0.0, ..mesh };
         assert_eq!(got[0], TaggedSurface { name: "Block".into(), mesh: mesh.clone(), faces: 6 });
         assert_eq!(got[1].name, "Again");
         assert_eq!(got[1].mesh, mesh);
@@ -577,7 +579,11 @@ mod tests {
         }
         assert!(read_3mf(b"").is_err());
         assert!(read_3mf(b"PK\x03\x04garbage").is_err());
-        let good = write_tagged(&[("T", &SurfaceMesh { positions: vec![[0.0; 3]; 3], triangles: vec![[0, 1, 2]], tri_face: vec![0] }, 1)]);
+        let good = write_tagged(&[(
+            "T",
+            &SurfaceMesh { positions: vec![[0.0; 3]; 3], triangles: vec![[0, 1, 2]], tri_face: vec![0], weld_distance: 0.0 },
+            1,
+        )]);
         for cut in [10, good.len() / 2, good.len() - 1] {
             assert!(read_3mf(&good[..cut]).is_err(), "truncated at {cut}");
         }
diff --git a/crates/mesh/src/lib.rs b/crates/mesh/src/lib.rs
index b1c0e78..1682305 100644
--- a/crates/mesh/src/lib.rs
+++ b/crates/mesh/src/lib.rs
@@ -9,6 +9,8 @@
 
 pub mod delaunay;
 pub mod quality;
+pub mod surface;
+mod vec3;
 
 #[derive(Debug, Clone, PartialEq, thiserror::Error)]
 pub enum Error {
@@ -18,6 +20,9 @@ pub enum Error {
     Limit(String),
     #[error("meshing failed: {0}")]
     Failed(String),
+    /// The boundary surface has defects ([`surface::BoundarySurface::new`]).
+    #[error("invalid boundary surface: {0}")]
+    Surface(surface::Defects),
     #[error(transparent)]
     Model(#[from] modelcraft_model::Error),
 }
diff --git a/crates/mesh/src/predicates/indirect.rs b/crates/mesh/src/predicates/indirect.rs
index 766c882..68bdc73 100644
--- a/crates/mesh/src/predicates/indirect.rs
+++ b/crates/mesh/src/predicates/indirect.rs
@@ -275,16 +275,23 @@ fn orient2d_formula<T: Ring + Signed>([a, b, c]: [Hom<T>; 3], (u, v): (usize, us
 /// seen from the positive side of that axis (`Z` → the `xy` plane, `X` → `yz`, `Y` → `zx`).
 /// Exact. Attene's `orient2d3d` (§4.5) is this with the axis of the largest normal component.
 pub fn generic_orient2d(a: &GenericPoint, b: &GenericPoint, c: &GenericPoint, drop: Axis) -> Sign {
-    let (u, v) = drop.projected();
     if let (Some(a), Some(b), Some(c)) = (a.explicit(), b.explicit(), c.explicit()) {
-        let proj = |p: &Point3| Point2([coord(&p.0, u), coord(&p.0, v)]);
-        return orient2d(&proj(a), &proj(b), &proj(c));
+        return orient2d_projected(a, b, c, drop);
     }
+    let (u, v) = drop.projected();
     let (points, flipped) = pivot_last([a, b, c]);
     let sign = decide(points, |h| orient2d_formula(h, (u, v)), |h| orient2d_formula(h, (u, v)));
     if flipped { -sign } else { sign }
 }
 
+/// [`orient2d`] of the projections of three explicit points onto the coordinate plane that drops
+/// `drop`, seen from the positive side of that axis (the explicit case of [`generic_orient2d`]).
+pub fn orient2d_projected(a: &Point3, b: &Point3, c: &Point3, drop: Axis) -> Sign {
+    let (u, v) = drop.projected();
+    let proj = |p: &Point3| Point2([coord(&p.0, u), coord(&p.0, v)]);
+    orient2d(&proj(a), &proj(b), &proj(c))
+}
+
 /// `X_a / W_a - X_b / W_b = (W_b X_a - W_a X_b) / (W_a W_b)` along one axis.
 fn compare_formula<T: Ring + Signed>([a, b]: [Hom<T>; 2], k: usize) -> Option<Sign> {
     let diff = mulw(coord(&a.x, k), &b.w) - mulw(coord(&b.x, k), &a.w);
diff --git a/crates/mesh/src/predicates/mod.rs b/crates/mesh/src/predicates/mod.rs
index 22848c4..7ff8fb4 100644
--- a/crates/mesh/src/predicates/mod.rs
+++ b/crates/mesh/src/predicates/mod.rs
@@ -40,7 +40,7 @@ mod ring;
 mod sos;
 
 pub use adaptive::{incircle, incircle_filter, insphere, insphere_filter, orient2d, orient2d_filter, orient3d, orient3d_filter};
-pub use indirect::{Axis, GenericPoint, Lpi, Tpi, generic_compare, generic_orient2d, generic_orient3d};
+pub use indirect::{Axis, GenericPoint, Lpi, Tpi, generic_compare, generic_orient2d, generic_orient3d, orient2d_projected};
 pub use sos::{incircle_sos, insphere_sos};
 
 use crate::{Error, Result};
diff --git a/crates/mesh/src/quality/shape.rs b/crates/mesh/src/quality/shape.rs
index ad610d7..b0d94c7 100644
--- a/crates/mesh/src/quality/shape.rs
+++ b/crates/mesh/src/quality/shape.rs
@@ -14,28 +14,8 @@ use std::f64::consts::SQRT_2;
 
 use super::Metric;
 
-pub(crate) type P = [f64; 3];
-
-pub(crate) fn sub(a: P, b: P) -> P {
-    [a[0] - b[0], a[1] - b[1], a[2] - b[2]]
-}
-pub(crate) fn dot(a: P, b: P) -> f64 {
-    a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
-}
-pub(crate) fn cross(a: P, b: P) -> P {
-    [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]]
-}
-pub(crate) fn norm(a: P) -> f64 {
-    dot(a, a).sqrt()
-}
-pub(crate) fn det(a: P, b: P, c: P) -> f64 {
-    dot(a, cross(b, c))
-}
-
-/// The angle between two vectors in degrees, 0 when either is zero.
-fn angle(a: P, b: P) -> f64 {
-    norm(cross(a, b)).atan2(dot(a, b)).to_degrees()
-}
+use crate::vec3::angle_deg as angle;
+pub(crate) use crate::vec3::{P, cross, det, dot, norm, sub};
 
 /// `num / den`, or `fallback` when the denominator is zero or the quotient is not finite.
 fn ratio(num: f64, den: f64, fallback: f64) -> f64 {
diff --git a/crates/mesh/src/surface/bvh.rs b/crates/mesh/src/surface/bvh.rs
new file mode 100644
index 0000000..39edc55
--- /dev/null
+++ b/crates/mesh/src/surface/bvh.rs
@@ -0,0 +1,163 @@
+//! The self-intersection search: a bounding volume hierarchy of axis-aligned triangle boxes,
+//! built top-down by median split on the longest centroid extent, and the self-collision
+//! traversal that pairs nodes whose boxes touch (Ericson 2005, ch. 6). Boxes are compared with
+//! closed `<=` on the exact coordinates, so touching triangles are always paired; every pair is
+//! then decided exactly ([`super::intersect::surface_pair_intersects`]).
+
+use super::intersect::surface_pair_intersects;
+use super::{Box3, Defect, Defects, MAX_PAIR_TESTS, Surface};
+use crate::{Error, Result};
+
+/// Triangles per leaf.
+const LEAF: usize = 4;
+
+/// Twice the centre of a box (no rounding from halving).
+fn centre2(b: &Box3) -> [f64; 3] {
+    [0, 1, 2].map(|k| b.lo[k] + b.hi[k])
+}
+
+#[derive(Debug, Clone, Copy)]
+struct Node {
+    bounds: Box3,
+    /// The left child (the right one follows it), or `None` for a leaf over `order[start..end]`.
+    children: Option<usize>,
+    start: usize,
+    end: usize,
+}
+
+struct Bvh {
+    nodes: Vec<Node>,
+    order: Vec<u32>,
+}
+
+fn build(boxes: &[Box3]) -> Bvh {
+    let mut order: Vec<u32> = (0..boxes.len()).map(super::index).collect();
+    let empty = Node { bounds: Box3::EMPTY, children: None, start: 0, end: 0 };
+    let mut nodes = vec![empty];
+    let mut work = vec![(0usize, 0usize, boxes.len())];
+    while let Some((ni, start, end)) = work.pop() {
+        let Some(slice) = order.get_mut(start..end) else { continue };
+        let mut bounds = Box3::EMPTY;
+        let mut centres = Box3::EMPTY;
+        for &t in slice.iter() {
+            if let Some(b) = boxes.get(t as usize) {
+                bounds.add(b);
+                centres.add_point(centre2(b));
+            }
+        }
+        let mut node = Node { bounds, children: None, start, end };
+        if end - start > LEAF {
+            let extent = [0, 1, 2].map(|k| centres.hi[k] - centres.lo[k]);
+            let axis = if extent[0] >= extent[1] && extent[0] >= extent[2] {
+                0
+            } else if extent[1] >= extent[2] {
+                1
+            } else {
+                2
+            };
+            let key = |t: &u32| boxes.get(*t as usize).map_or(0.0, |b| centre2(b)[axis]);
+            let mid = (end - start) / 2;
+            slice.select_nth_unstable_by(mid, |a, b| key(a).total_cmp(&key(b)).then(a.cmp(b)));
+            let l = nodes.len();
+            nodes.push(empty);
+            nodes.push(empty);
+            node.children = Some(l);
+            work.push((l, start, start + mid));
+            work.push((l + 1, start + mid, end));
+        }
+        if let Some(n) = nodes.get_mut(ni) {
+            *n = node;
+        }
+    }
+    Bvh { nodes, order }
+}
+
+struct Search<'a> {
+    surface: &'a Surface<'a>,
+    boxes: &'a [Box3],
+    tests: u64,
+    defects: Defects,
+}
+
+impl Search<'_> {
+    /// Test one pair; `false` once the defect list is full.
+    fn pair(&mut self, s: u32, t: u32) -> Result<bool> {
+        let (Some(bs), Some(bt)) = (self.boxes.get(s as usize), self.boxes.get(t as usize)) else { return Ok(true) };
+        if s == t || !bs.touches(bt) {
+            return Ok(true);
+        }
+        self.tests += 1;
+        if self.tests > MAX_PAIR_TESTS {
+            return Err(Error::Limit(format!("self-intersection check: more than {MAX_PAIR_TESTS} triangle pairs to test")));
+        }
+        let (Some(ts), Some(tt)) = (self.surface.triangles.get(s as usize), self.surface.triangles.get(t as usize)) else { return Ok(true) };
+        if surface_pair_intersects([*ts, *tt], self.surface.corners(ts)?, self.surface.corners(tt)?) {
+            let (s, t) = (s.min(t) as usize, s.max(t) as usize);
+            self.defects
+                .push(Defect::SelfIntersection { triangles: [self.surface.source_tri(s), self.surface.source_tri(t)], at: self.surface.centroid(s) });
+        }
+        Ok(!self.defects.is_full())
+    }
+}
+
+/// Every pair of triangles that meet outside what they share, up to the listing cap.
+pub(crate) fn self_intersections(surface: &Surface) -> Result<Defects> {
+    let boxes: Vec<Box3> = surface
+        .triangles
+        .iter()
+        .map(|t| {
+            let mut b = Box3::EMPTY;
+            for v in t {
+                b.add_point(surface.xyz(*v));
+            }
+            b
+        })
+        .collect();
+    let bvh = build(&boxes);
+    let mut search = Search { surface, boxes: &boxes, tests: 0, defects: Defects::new() };
+    let leaf = |n: &Node| bvh.order.get(n.start..n.end).unwrap_or(&[]);
+    let mut stack = vec![(0usize, 0usize)];
+    'search: while let Some((a, b)) = stack.pop() {
+        let (Some(na), Some(nb)) = (bvh.nodes.get(a), bvh.nodes.get(b)) else { continue };
+        if a == b {
+            match na.children {
+                Some(l) => stack.extend([(l, l + 1), (l + 1, l + 1), (l, l)]),
+                None => {
+                    let tris = leaf(na);
+                    for (i, &s) in tris.iter().enumerate() {
+                        for &t in tris.iter().skip(i + 1) {
+                            if !search.pair(s, t)? {
+                                break 'search;
+                            }
+                        }
+                    }
+                }
+            }
+            continue;
+        }
+        if !na.bounds.touches(&nb.bounds) {
+            continue;
+        }
+        match (na.children, nb.children) {
+            (None, None) => {
+                for &s in leaf(na) {
+                    for &t in leaf(nb) {
+                        if !search.pair(s, t)? {
+                            break 'search;
+                        }
+                    }
+                }
+            }
+            (Some(l), None) => stack.extend([(l + 1, b), (l, b)]),
+            (None, Some(l)) => stack.extend([(a, l + 1), (a, l)]),
+            (Some(l), Some(_)) if na.end - na.start >= nb.end - nb.start => stack.extend([(l + 1, b), (l, b)]),
+            (Some(_), Some(l)) => stack.extend([(a, l + 1), (a, l)]),
+        }
+    }
+    let mut defects = search.defects;
+    defects.listed.sort_by_key(|d| match d {
+        Defect::SelfIntersection { triangles, .. } => *triangles,
+        _ => [u32::MAX; 2],
+    });
+    Ok(defects)
+}
diff --git a/crates/mesh/src/surface/intersect.rs b/crates/mesh/src/surface/intersect.rs
new file mode 100644
index 0000000..ad55f6a
--- /dev/null
+++ b/crates/mesh/src/surface/intersect.rs
@@ -0,0 +1,160 @@
+//! Exact triangle–triangle intersection on the orientation predicates.
+//!
+//! Every decision is a sign of [`orient3d`] or of [`orient2d_projected`] on a coordinate-plane
+//! projection, so the answers are exact for every admissible input (Shewchuk 1997; the
+//! orientation-only formulation follows Guigue and Devillers 2003). Triangles are closed sets.
+//!
+//! Two closed triangles meet exactly when an edge of one meets the other: their intersection is
+//! a convex set whose extreme points lie on the boundary of one of them. Two triangles of a
+//! surface may meet only in what they share, which [`surface_pair_intersects`] tests:
+//! - **No shared vertex:** any contact is a defect ([`triangles_meet`]).
+//! - **One shared vertex `v`:** the intersection is a convex set containing `v`; it holds another
+//!   point exactly when an edge opposite `v` meets the other triangle, or an edge from `v` runs
+//!   into the other triangle's corner at `v`.
+//! - **One shared edge:** the triangles meet beyond it only when they are coplanar and on the
+//!   same side of it (folded onto each other).
+//!
+//! Every vertex's side of the other triangle's plane is computed once per pair; a triangle
+//! strictly on one side of the other's plane rejects the pair at once.
+
+use crate::predicates::{Axis, Point3, Sign, orient2d_projected, orient3d};
+
+type Tri<'a> = [&'a Point3; 3];
+
+/// A coordinate plane on which triangle `abc` projects to a proper triangle, and the
+/// orientation of that projection; `None` when the corners are collinear.
+fn plane_axis([a, b, c]: Tri) -> Option<(Axis, Sign)> {
+    [Axis::Z, Axis::X, Axis::Y].into_iter().find_map(|axis| {
+        let s = orient2d_projected(a, b, c, axis);
+        (!s.is_zero()).then_some((axis, s))
+    })
+}
+
+/// Whether the corners are exactly collinear (or coincide).
+pub fn is_degenerate(t: Tri) -> bool {
+    plane_axis(t).is_none()
+}
+
+fn edges<'a>([a, b, c]: Tri<'a>) -> [(&'a Point3, &'a Point3); 3] {
+    [(a, b), (b, c), (c, a)]
+}
+
+/// The side of `p` against the plane of `t`.
+fn side(t: Tri, p: &Point3) -> Sign {
+    orient3d(t[0], t[1], t[2], p)
+}
+
+/// Closed point-in-triangle test in the plane of `t` (`p` must lie in it).
+fn point_in_triangle(p: &Point3, t: Tri, drop: Axis, s: Sign) -> bool {
+    edges(t).into_iter().all(|(u, v)| orient2d_projected(u, v, p, drop) != -s)
+}
+
+/// Whether `r`, collinear with `p` and `q`, lies on the closed segment `pq`.
+fn on_segment(p: &Point3, q: &Point3, r: &Point3) -> bool {
+    let (p, q, r) = (p.coords(), q.coords(), r.coords());
+    (0..3).all(|k| {
+        let (lo, hi) = if p[k] <= q[k] { (p[k], q[k]) } else { (q[k], p[k]) };
+        lo <= r[k] && r[k] <= hi
+    })
+}
+
+/// Closed segment–segment intersection of two coplanar segments, on the projection `drop`
+/// (which must be injective on their plane).
+fn segments_meet(p: &Point3, q: &Point3, a: &Point3, b: &Point3, drop: Axis) -> bool {
+    let d1 = orient2d_projected(p, q, a, drop);
+    let d2 = orient2d_projected(p, q, b, drop);
+    let d3 = orient2d_projected(a, b, p, drop);
+    let d4 = orient2d_projected(a, b, q, drop);
+    (d1 * d2 == Sign::Negative && d3 * d4 == Sign::Negative)
+        || (d1.is_zero() && on_segment(p, q, a))
+        || (d2.is_zero() && on_segment(p, q, b))
+        || (d3.is_zero() && on_segment(a, b, p))
+        || (d4.is_zero() && on_segment(a, b, q))
+}
+
+/// Whether the closed segment `pq` meets the closed, non-degenerate triangle `t`, given the
+/// sides `op`, `oq` of `p` and `q` against `t`'s plane.
+fn segment_meets(p: &Point3, q: &Point3, op: Sign, oq: Sign, t: Tri) -> bool {
+    if op == oq && !op.is_zero() {
+        return false;
+    }
+    if op.is_zero() && oq.is_zero() {
+        let Some((drop, s)) = plane_axis(t) else { return false };
+        return point_in_triangle(p, t, drop, s)
+            || point_in_triangle(q, t, drop, s)
+            || edges(t).into_iter().any(|(u, v)| segments_meet(p, q, u, v, drop));
+    }
+    // The segment crosses or touches the plane at one point, which is inside the triangle
+    // exactly when the line pq passes on no strict outside of its three edges.
+    let s = edges(t).map(|(u, v)| orient3d(p, q, u, v));
+    !(s.contains(&Sign::Positive) && s.contains(&Sign::Negative))
+}
+
+/// Whether the closed segment `pq` meets the closed, non-degenerate triangle `t`.
+pub fn segment_meets_triangle(p: &Point3, q: &Point3, t: Tri) -> bool {
+    segment_meets(p, q, side(t, p), side(t, q), t)
+}
+
+/// Whether two closed, non-degenerate triangles intersect at all.
+pub fn triangles_meet(s: Tri, t: Tri) -> bool {
+    let s_sides = s.map(|p| side(t, p));
+    let t_sides = t.map(|p| side(s, p));
+    let one_side = |x: [Sign; 3]| !x[0].is_zero() && x[0] == x[1] && x[1] == x[2];
+    if one_side(s_sides) || one_side(t_sides) {
+        return false;
+    }
+    let crosses = |a: Tri, o: [Sign; 3], b: Tri| (0..3).any(|k| segment_meets(a[k], a[(k + 1) % 3], o[k], o[(k + 1) % 3], b));
+    crosses(s, s_sides, t) || crosses(t, t_sides, s)
+}
+
+/// For triangles sharing the corner `v`: whether the edge from `v` to `a` (whose side of the
+/// plane of triangle `v c d` is `side_a`) contains a point of that triangle other than `v`: it
+/// lies in that plane, inside the triangle's corner at `v`.
+fn edge_enters_corner(v: &Point3, a: &Point3, side_a: Sign, c: &Point3, d: &Point3) -> bool {
+    if !side_a.is_zero() {
+        return false;
+    }
+    let Some((drop, s)) = plane_axis([v, c, d]) else { return false };
+    orient2d_projected(v, c, a, drop) != -s && orient2d_projected(v, a, d, drop) != -s
+}
+
+/// Whether two triangles of a surface, given by vertex ids and points, meet anywhere other than
+/// in the vertices and edge they share. Both must be non-degenerate.
+pub fn surface_pair_intersects(ids: [[u32; 3]; 2], s: Tri, t: Tri) -> bool {
+    let [si, ti] = ids;
+    // Positions in s of the corners shared with t, and their positions in t.
+    let mut shared = [(0usize, 0usize); 3];
+    let mut n = 0;
+    for (i, v) in si.iter().enumerate() {
+        if let (Some(j), Some(slot)) = (ti.iter().position(|w| w == v), shared.get_mut(n)) {
+            *slot = (i, j);
+            n += 1;
+        }
+    }
+    match shared.get(..n).unwrap_or(&[]) {
+        [] => triangles_meet(s, t),
+        &[(i, j)] => {
+            let v = s[i];
+            let (a, b) = (s[(i + 1) % 3], s[(i + 2) % 3]);
+            let (c, d) = (t[(j + 1) % 3], t[(j + 2) % 3]);
+            let (oa, ob, oc, od) = (side(t, a), side(t, b), side(s, c), side(s, d));
+            segment_meets(a, b, oa, ob, t)
+                || segment_meets(c, d, oc, od, s)
+                || edge_enters_corner(v, a, oa, c, d)
+                || edge_enters_corner(v, b, ob, c, d)
+                || edge_enters_corner(v, c, oc, a, b)
+                || edge_enters_corner(v, d, od, a, b)
+        }
+        &[(i0, j0), (i1, j1)] => {
+            // The corner of each triangle off the shared edge.
+            let (u, w, a, b) = (s[i0], s[i1], s[3 - i0 - i1], t[3 - j0 - j1]);
+            if !orient3d(u, w, a, b).is_zero() {
+                return false;
+            }
+            let Some((drop, _)) = plane_axis(s) else { return false };
+            orient2d_projected(u, w, a, drop) == orient2d_projected(u, w, b, drop)
+        }
+        // The same three vertices: a duplicate triangle.
+        _ => true,
+    }
+}
diff --git a/crates/mesh/src/surface/mod.rs b/crates/mesh/src/surface/mod.rs
new file mode 100644
index 0000000..462c911
--- /dev/null
+++ b/crates/mesh/src/surface/mod.rs
@@ -0,0 +1,535 @@
+//! Boundary surface preparation: the step between a tessellated body and volume meshing.
+//!
+//! [`BoundarySurface::new`] turns a face-tagged triangle surface (SolidCraft's tessellation, or a
+//! tagged 3MF) into one that every volume mesher can rely on, or fails with the list of
+//! [`Defect`]s, each located by input triangle or vertex index and a position:
+//!
+//! 1. **Coordinates and indices** are checked: finite, within [`modelcraft_model::MAX_COORD`],
+//!    in range.
+//! 2. **Weld.** Vertices closer than a distance (relative to the bounding-box diagonal, or
+//!    absolute: [`SurfaceOptions`]) become one, in input order (the earliest vertex of a cluster
+//!    is its representative and keeps its position), so the result does not depend on hashing.
+//!    SolidCraft's tessellation is a triangle soup (three vertices per triangle), and two faces'
+//!    copies of a shared boundary point agree only to within the kernel's tolerance, so this is
+//!    where its faces join. Triangles that collapse are dropped; coordinates below
+//!    [`crate::predicates::MIN_MAGNITUDE`] snap to zero so every point is in the exact
+//!    predicates' domain.
+//! 3. **Degenerate triangles** (exactly collinear corners) are rejected.
+//! 4. **Topology.** Every edge has exactly two triangles that traverse it in opposite directions
+//!    (closed, edge-manifold, consistently oriented), and the triangles around each vertex form a
+//!    single fan (vertex-manifold).
+//! 5. **Self-intersection.** A bounding volume hierarchy pairs triangles whose boxes touch; each
+//!    pair is decided exactly ([`intersect`]): triangles may meet only in the vertices or the edge
+//!    they share.
+//! 6. **Shells and orientation.** Connected components, their genus from the Euler
+//!    characteristic, their signed volume (divergence theorem), and their nesting from the
+//!    generalized winding number: an outer shell must enclose positive volume, a void (a shell
+//!    inside an odd number of others) negative.
+//! 7. **Feature edges**: edges between different B-rep faces, and edges whose dihedral angle
+//!    exceeds a threshold. The B-rep face of every triangle is kept, so loads and boundary
+//!    conditions on CAD faces map to mesh faces.
+//!
+//! # Sources
+//! - J. R. Shewchuk, *Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric
+//!   Predicates*, DCG 18:305–363, 1997: every intersection decision is a sign of `orient3d` or
+//!   `orient2d` ([`crate::predicates`]).
+//! - P. Guigue, O. Devillers, *Fast and Robust Triangle–Triangle Overlap Test Using Orientation
+//!   Predicates*, J. Graphics Tools 8(1):25–32, 2003: triangle overlap decided from orientation
+//!   signs alone (here through the edge–triangle reduction in [`intersect`]).
+//! - C. Ericson, *Real-Time Collision Detection*, Morgan Kaufmann, 2005, ch. 6: bounding volume
+//!   hierarchies, top-down median-split construction and the self-collision traversal.
+//! - A. Jacobson, L. Kavan, O. Sorkine-Hornung, *Robust Inside-Outside Segmentation using
+//!   Generalized Winding Numbers*, ACM TOG 32(4), 2013, with the triangle solid angle of A. van
+//!   Oosterom, J. Strackee, *The Solid Angle of a Plane Triangle*, IEEE TBME 30(2):125–126, 1983:
+//!   shell nesting.
+
+mod bvh;
+pub mod intersect;
+mod topology;
+mod weld;
+
+#[cfg(test)]
+mod tests;
+
+use std::fmt;
+
+use modelcraft_model::{MAX_COORD, MAX_ELEMENTS, MAX_NODES};
+
+use crate::predicates::Point3;
+use crate::vec3::{cross, norm, sub};
+use crate::{Error, Result};
+
+/// Defects listed in an error; the count covers the rest.
+pub const MAX_LISTED_DEFECTS: usize = 100;
+/// Cap on exact triangle-pair tests in the self-intersection check.
+pub const MAX_PAIR_TESTS: u64 = 1_000_000_000;
+/// Cap on triangle solid angles summed to nest shells (shells × triangles, after box culling).
+pub const MAX_WINDING_WORK: u64 = 500_000_000;
+
+/// The largest weld distance, as a fraction of the bounding-box diagonal.
+pub const MAX_RELATIVE_WELD: f64 = 1e-3;
+
+/// How [`BoundarySurface::new`] welds and finds features.
+///
+/// Vertices closer than the larger of `weld_tolerance` × the bounding-box diagonal and
+/// `weld_distance` are welded. Tessellations of B-rep faces agree on their shared boundary points
+/// only to within the kernel's tolerance, so a caller with CAD input sets `weld_distance` from it
+/// (twice the kernel's linear tolerance: each face's copy of a boundary point is within it of the
+/// exact point).
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct SurfaceOptions {
+    /// Relative weld distance, in `[0, 1e-3]`; nonzero values below 1e-15 are refused (finer than
+    /// the coordinates themselves). 0 with `weld_distance` 0 welds identical coordinates only.
+    pub weld_tolerance: f64,
+    /// Absolute weld distance in model units, at least 0 (the combined distance may be at most
+    /// [`MAX_RELATIVE_WELD`] of the bounding-box diagonal). A tessellation carries its own
+    /// (`modelcraft_geometry::SurfaceMesh::weld_distance`).
+    pub weld_distance: f64,
+    /// An edge inside one B-rep face is a feature when the normals of its two triangles differ
+    /// by more than this angle, in degrees, in `(0, 180]` (180: face boundaries only).
+    pub feature_angle: f64,
+}
+
+impl Default for SurfaceOptions {
+    fn default() -> Self {
+        SurfaceOptions { weld_tolerance: 1e-9, weld_distance: 0.0, feature_angle: 45.0 }
+    }
+}
+
+/// A defect that keeps a triangle surface from bounding a volume mesh. Triangle and vertex
+/// numbers are indices into the input; `at` is a position near the defect.
+#[derive(Debug, Clone, PartialEq)]
+pub enum Defect {
+    /// A coordinate is not a finite number within ±[`MAX_COORD`].
+    BadCoordinate { vertex: u32 },
+    /// A triangle refers to a vertex that does not exist.
+    BadIndex { triangle: u32 },
+    /// A triangle with exactly collinear corners.
+    Degenerate { triangle: u32, at: [f64; 3] },
+    /// An edge with only one triangle: the surface has a hole there.
+    OpenEdge { vertices: [u32; 2], triangle: u32, at: [f64; 3] },
+    /// An edge shared by more than two triangles.
+    NonManifoldEdge { vertices: [u32; 2], triangles: u32, at: [f64; 3] },
+    /// Two triangles traverse their shared edge in the same direction: one of them is flipped.
+    Flipped { vertices: [u32; 2], triangles: [u32; 2], at: [f64; 3] },
+    /// The triangles around a vertex form more than one fan (two sheets touch at a point).
+    NonManifoldVertex { vertex: u32, at: [f64; 3] },
+    /// Two triangles meet outside the vertices and edge they share.
+    SelfIntersection { triangles: [u32; 2], at: [f64; 3] },
+    /// A shell enclosing negative volume although it is not inside another shell (its triangles
+    /// face inwards), or positive volume although it is a void.
+    Inverted { shell: u32, triangle: u32, volume: f64 },
+    /// A shell whose enclosed volume is too small to tell its orientation from rounding.
+    Flat { shell: u32, triangle: u32 },
+}
+
+impl Defect {
+    /// A short, stable name for the kind of defect (`open_edge`, `self_intersection`, …).
+    pub fn kind(&self) -> &'static str {
+        match self {
+            Defect::BadCoordinate { .. } => "bad_coordinate",
+            Defect::BadIndex { .. } => "bad_index",
+            Defect::Degenerate { .. } => "degenerate",
+            Defect::OpenEdge { .. } => "open_edge",
+            Defect::NonManifoldEdge { .. } => "non_manifold_edge",
+            Defect::Flipped { .. } => "flipped",
+            Defect::NonManifoldVertex { .. } => "non_manifold_vertex",
+            Defect::SelfIntersection { .. } => "self_intersection",
+            Defect::Inverted { .. } => "inverted",
+            Defect::Flat { .. } => "flat",
+        }
+    }
+
+    /// A position near the defect, when it has one.
+    pub fn at(&self) -> Option<[f64; 3]> {
+        match self {
+            Defect::BadCoordinate { .. } | Defect::BadIndex { .. } | Defect::Inverted { .. } | Defect::Flat { .. } => None,
+            Defect::Degenerate { at, .. }
+            | Defect::OpenEdge { at, .. }
+            | Defect::NonManifoldEdge { at, .. }
+            | Defect::Flipped { at, .. }
+            | Defect::NonManifoldVertex { at, .. }
+            | Defect::SelfIntersection { at, .. } => Some(*at),
+        }
+    }
+}
+
+struct Xyz([f64; 3]);
+impl fmt::Display for Xyz {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(f, "({}, {}, {})", self.0[0], self.0[1], self.0[2])
+    }
+}
+
+impl fmt::Display for Defect {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        match *self {
+            Defect::BadCoordinate { vertex } => write!(f, "vertex {vertex}: a coordinate is not a finite number within ±{MAX_COORD}"),
+            Defect::BadIndex { triangle } => write!(f, "triangle {triangle} refers to a vertex that does not exist"),
+            Defect::Degenerate { triangle, at } => write!(f, "triangle {triangle} at {} has collinear corners", Xyz(at)),
+            Defect::OpenEdge { vertices: [a, b], triangle, at } => {
+                write!(f, "open edge {a}–{b} at {} (only triangle {triangle} uses it)", Xyz(at))
+            }
+            Defect::NonManifoldEdge { vertices: [a, b], triangles, at } => {
+                write!(f, "edge {a}–{b} at {} is shared by {triangles} triangles", Xyz(at))
+            }
+            Defect::Flipped { vertices: [a, b], triangles: [s, t], at } => {
+                write!(f, "triangles {s} and {t} are oriented inconsistently across edge {a}–{b} at {}", Xyz(at))
+            }
+            Defect::NonManifoldVertex { vertex, at } => write!(f, "vertex {vertex} at {} joins separate sheets", Xyz(at)),
+            Defect::SelfIntersection { triangles: [s, t], at } => write!(f, "triangles {s} and {t} intersect near {}", Xyz(at)),
+            Defect::Inverted { shell, triangle, volume } => {
+                write!(f, "shell {shell} (from triangle {triangle}) is inside out: it encloses {volume} where the opposite sign is expected")
+            }
+            Defect::Flat { shell, triangle } => write!(f, "shell {shell} (from triangle {triangle}) encloses no measurable volume"),
+        }
+    }
+}
+
+/// The defects found by the first check that failed: up to [`MAX_LISTED_DEFECTS`], and how many
+/// there are in all (a lower bound for self-intersections, whose search stops at the cap).
+#[derive(Debug, Clone, PartialEq)]
+pub struct Defects {
+    pub listed: Vec<Defect>,
+    pub total: usize,
+}
+
+impl Defects {
+    fn new() -> Self {
+        Defects { listed: Vec::new(), total: 0 }
+    }
+
+    fn push(&mut self, d: Defect) {
+        self.total += 1;
+        if self.listed.len() < MAX_LISTED_DEFECTS {
+            self.listed.push(d);
+        }
+    }
+
+    fn is_full(&self) -> bool {
+        self.listed.len() >= MAX_LISTED_DEFECTS
+    }
+
+    /// `Err` when there is at least one defect.
+    fn into_result(self) -> Result<()> {
+        if self.total == 0 { Ok(()) } else { Err(Error::Surface(self)) }
+    }
+}
+
+impl fmt::Display for Defects {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        let s = if self.total == 1 { "" } else { "s" };
+        write!(f, "{} defect{s}", self.total)?;
+        for d in self.listed.iter().take(3) {
+            write!(f, "; {d}")?;
+        }
+        if self.total > 3 {
+            write!(f, "; …")?;
+        }
+        Ok(())
+    }
+}
+
+/// A connected component of the surface.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct Shell {
+    /// The shell's first triangle (output index); [`BoundarySurface::tri_shell`] lists the rest.
+    pub first_triangle: u32,
+    pub triangles: u32,
+    pub vertices: u32,
+    /// Handles (through holes): 0 for a block, 1 for a plate with a hole.
+    pub genus: u32,
+    /// Signed enclosed volume: positive for an outer shell, negative for a void.
+    pub volume: f64,
+    /// Inside an odd number of other shells: a cavity in the solid.
+    pub void: bool,
+}
+
+/// Why an edge is a feature.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
+pub enum FeatureKind {
+    /// Between two B-rep faces.
+    FaceBoundary,
+    /// Inside one face, with a dihedral angle above the threshold.
+    Sharp,
+}
+
+/// An edge the volume mesh must keep: its vertices (ascending), its two triangles (the one that
+/// traverses it from `vertices[0]` to `vertices[1]` first) and the angle between their normals
+/// in degrees.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct FeatureEdge {
+    pub vertices: [u32; 2],
+    pub triangles: [u32; 2],
+    pub angle: f64,
+    pub kind: FeatureKind,
+}
+
+/// A welded, closed, 2-manifold, outward-oriented, self-intersection-free triangle surface with
+/// its B-rep face tags and feature edges: the boundary every volume mesher starts from.
+///
+/// Every point is admissible for [`crate::predicates`] ([`BoundarySurface::point3`]), triangles
+/// are counter-clockwise seen from outside the solid (voids included: their normals point into
+/// the cavity), and all output is a deterministic function of the input.
+#[derive(Debug, Clone, PartialEq)]
+pub struct BoundarySurface {
+    points: Vec<[f64; 3]>,
+    triangles: Vec<[u32; 3]>,
+    tri_face: Vec<u32>,
+    neighbors: Vec<[u32; 3]>,
+    tri_shell: Vec<u32>,
+    shells: Vec<Shell>,
+    feature_edges: Vec<FeatureEdge>,
+    source_vertex: Vec<u32>,
+    source_triangle: Vec<u32>,
+}
+
+impl BoundarySurface {
+    /// Prepare `triangles` over `positions` (counter-clockwise seen from outside) with the B-rep
+    /// face of each triangle (`tri_face`; empty means one face). Fails with
+    /// [`Error::Surface`] listing the defects of the first check that fails, with
+    /// [`Error::Invalid`] on unusable options or sizes, or with [`Error::Limit`] when a search cap
+    /// is reached.
+    pub fn new(positions: &[[f64; 3]], triangles: &[[u32; 3]], tri_face: &[u32], options: &SurfaceOptions) -> Result<Self> {
+        check_options(options)?;
+        if positions.len() > MAX_NODES || triangles.len() > MAX_ELEMENTS {
+            return Err(Error::Invalid(format!(
+                "{} vertices and {} triangles exceed the limits ({MAX_NODES}, {MAX_ELEMENTS})",
+                positions.len(),
+                triangles.len()
+            )));
+        }
+        if !tri_face.is_empty() && tri_face.len() != triangles.len() {
+            return Err(Error::Invalid(format!("{} face tags for {} triangles", tri_face.len(), triangles.len())));
+        }
+        if triangles.is_empty() {
+            return Err(Error::Invalid("the surface has no triangles".into()));
+        }
+        check_input(positions, triangles)?;
+
+        let welded = weld::weld(positions, triangles, options.weld_tolerance, options.weld_distance)?;
+        if welded.triangles.is_empty() {
+            return Err(Error::Invalid("every triangle collapses when the vertices are welded".into()));
+        }
+        let pts: Vec<Point3> = welded.points.iter().map(|&p| Point3::try_from(p)).collect::<Result<_>>()?;
+        let tri_face: Vec<u32> = welded.source_triangle.iter().map(|&t| tri_face.get(t as usize).copied().unwrap_or(0)).collect();
+        let surface =
+            Surface { pts: &pts, triangles: &welded.triangles, source_vertex: &welded.source_vertex, source_triangle: &welded.source_triangle };
+
+        let mut defects = Defects::new();
+        for (t, tri) in surface.triangles.iter().enumerate() {
+            if intersect::is_degenerate(surface.corners(tri)?) {
+                defects.push(Defect::Degenerate { triangle: surface.source_tri(t), at: surface.centroid(t) });
+            }
+        }
+        defects.into_result()?;
+
+        let topo = topology::Topology::build(&surface)?;
+        bvh::self_intersections(&surface)?.into_result()?;
+        let shells = topology::shells(&surface, &topo)?;
+        let feature_edges = topology::feature_edges(&surface, &topo, &tri_face, options.feature_angle);
+        let weld::Welded { points, triangles, source_vertex, source_triangle } = welded;
+        Ok(BoundarySurface {
+            points,
+            triangles,
+            tri_face,
+            neighbors: topo.neighbors,
+            tri_shell: shells.tri_shell,
+            shells: shells.shells,
+            feature_edges,
+            source_vertex,
+            source_triangle,
+        })
+    }
+
+    /// Welded vertex positions.
+    pub fn points(&self) -> &[[f64; 3]] {
+        &self.points
+    }
+
+    /// Vertex `v` as a predicate point.
+    pub fn point3(&self, v: u32) -> Option<Point3> {
+        self.points.get(v as usize).and_then(|&p| Point3::try_from(p).ok())
+    }
+
+    /// Triangles over [`Self::points`], counter-clockwise seen from outside the solid.
+    pub fn triangles(&self) -> &[[u32; 3]] {
+        &self.triangles
+    }
+
+    /// The B-rep face of each triangle (the input's tag, 0 when there was none).
+    pub fn tri_face(&self) -> &[u32] {
+        &self.tri_face
+    }
+
+    /// For each triangle `[a, b, c]`, the triangles across its edges `ab`, `bc` and `ca`.
+    pub fn neighbors(&self) -> &[[u32; 3]] {
+        &self.neighbors
+    }
+
+    /// The shell of each triangle (an index into [`Self::shells`]).
+    pub fn tri_shell(&self) -> &[u32] {
+        &self.tri_shell
+    }
+
+    /// Shells in order of their first triangle.
+    pub fn shells(&self) -> &[Shell] {
+        &self.shells
+    }
+
+    /// The enclosed volume: the sum of the shells' signed volumes.
+    pub fn volume(&self) -> f64 {
+        self.shells.iter().map(|s| s.volume).sum()
+    }
+
+    /// Feature edges, sorted by vertices.
+    pub fn feature_edges(&self) -> &[FeatureEdge] {
+        &self.feature_edges
+    }
+
+    /// Feature corners: vertices where the number of feature edges is neither 0 nor 2 (where
+    /// three faces meet, or where a feature curve ends), ascending.
+    pub fn feature_corners(&self) -> Vec<u32> {
+        let mut count = vec![0u32; self.points.len()];
+        for e in &self.feature_edges {
+            for v in e.vertices {
+                if let Some(c) = count.get_mut(v as usize) {
+                    *c += 1;
+                }
+            }
+        }
+        count.iter().enumerate().filter(|(_, c)| !matches!(c, 0 | 2)).filter_map(|(v, _)| u32::try_from(v).ok()).collect()
+    }
+
+    /// The input vertex each point comes from (the first of the vertices welded into it).
+    pub fn source_vertex(&self) -> &[u32] {
+        &self.source_vertex
+    }
+
+    /// The input triangle each triangle comes from (collapsed ones are dropped).
+    pub fn source_triangle(&self) -> &[u32] {
+        &self.source_triangle
+    }
+}
+
+fn check_options(o: &SurfaceOptions) -> Result<()> {
+    let w = o.weld_tolerance;
+    if !(w == 0.0 || (1e-15..=MAX_RELATIVE_WELD).contains(&w)) {
+        return Err(Error::Invalid(format!(
+            "weld tolerance must be 0 or in [1e-15, {MAX_RELATIVE_WELD}] (relative to the bounding-box diagonal), got {w}"
+        )));
+    }
+    if !(o.weld_distance >= 0.0 && o.weld_distance <= MAX_COORD) {
+        return Err(Error::Invalid(format!("weld distance must be a length of at least 0, got {}", o.weld_distance)));
+    }
+    if !(o.feature_angle > 0.0 && o.feature_angle <= 180.0) {
+        return Err(Error::Invalid(format!("feature angle must be in (0, 180] degrees, got {}", o.feature_angle)));
+    }
+    Ok(())
+}
+
+fn check_input(positions: &[[f64; 3]], triangles: &[[u32; 3]]) -> Result<()> {
+    let mut defects = Defects::new();
+    for (v, p) in positions.iter().enumerate() {
+        if p.iter().any(|c| !(c.abs() <= MAX_COORD)) {
+            defects.push(Defect::BadCoordinate { vertex: index(v) });
+        }
+    }
+    for (t, tri) in triangles.iter().enumerate() {
+        if tri.iter().any(|&v| v as usize >= positions.len()) {
+            defects.push(Defect::BadIndex { triangle: index(t) });
+        }
+    }
+    defects.into_result()
+}
+
+/// An index known to fit (sizes are capped well below `u32::MAX`).
+pub(crate) fn index(i: usize) -> u32 {
+    u32::try_from(i).unwrap_or(u32::MAX)
+}
+
+/// The welded surface the checks run on.
+pub(crate) struct Surface<'a> {
+    pts: &'a [Point3],
+    triangles: &'a [[u32; 3]],
+    source_vertex: &'a [u32],
+    source_triangle: &'a [u32],
+}
+
+impl Surface<'_> {
+    fn point(&self, v: u32) -> Result<&Point3> {
+        self.pts.get(v as usize).ok_or_else(|| Error::Failed(format!("surface: no vertex {v}")))
+    }
+
+    fn corners(&self, t: &[u32; 3]) -> Result<[&Point3; 3]> {
+        Ok([self.point(t[0])?, self.point(t[1])?, self.point(t[2])?])
+    }
+
+    fn xyz(&self, v: u32) -> [f64; 3] {
+        self.pts.get(v as usize).map_or([0.0; 3], Point3::coords)
+    }
+
+    /// The corner positions of triangle `t`.
+    fn tri_xyz(&self, t: usize) -> [[f64; 3]; 3] {
+        self.triangles.get(t).copied().unwrap_or([0; 3]).map(|v| self.xyz(v))
+    }
+
+    fn centroid(&self, t: usize) -> [f64; 3] {
+        let [a, b, c] = self.tri_xyz(t);
+        [(a[0] + b[0] + c[0]) / 3.0, (a[1] + b[1] + c[1]) / 3.0, (a[2] + b[2] + c[2]) / 3.0]
+    }
+
+    fn midpoint(&self, a: u32, b: u32) -> [f64; 3] {
+        let (a, b) = (self.xyz(a), self.xyz(b));
+        [(a[0] + b[0]) / 2.0, (a[1] + b[1]) / 2.0, (a[2] + b[2]) / 2.0]
+    }
+
+    fn source_tri(&self, t: usize) -> u32 {
+        self.source_triangle.get(t).copied().unwrap_or(u32::MAX)
+    }
+
+    fn source_vert(&self, v: u32) -> u32 {
+        self.source_vertex.get(v as usize).copied().unwrap_or(u32::MAX)
+    }
+
+    /// The unnormalized normal of triangle `t` (floating point; for angles only).
+    fn normal(&self, t: usize) -> [f64; 3] {
+        let [a, b, c] = self.tri_xyz(t);
+        cross(sub(b, a), sub(c, a))
+    }
+}
+
+/// An axis-aligned box (empty until a point is added).
+#[derive(Debug, Clone, Copy)]
+pub(crate) struct Box3 {
+    pub lo: [f64; 3],
+    pub hi: [f64; 3],
+}
+
+impl Box3 {
+    pub const EMPTY: Box3 = Box3 { lo: [f64::INFINITY; 3], hi: [f64::NEG_INFINITY; 3] };
+
+    pub fn add_point(&mut self, p: [f64; 3]) {
+        self.add(&Box3 { lo: p, hi: p });
+    }
+
+    pub fn add(&mut self, o: &Box3) {
+        for ((lo, hi), (olo, ohi)) in self.lo.iter_mut().zip(&mut self.hi).zip(o.lo.iter().zip(&o.hi)) {
+            *lo = lo.min(*olo);
+            *hi = hi.max(*ohi);
+        }
+    }
+
+    /// Closed overlap.
+    pub fn touches(&self, o: &Box3) -> bool {
+        (0..3).all(|k| self.lo[k] <= o.hi[k] && o.lo[k] <= self.hi[k])
+    }
+
+    pub fn contains(&self, p: [f64; 3]) -> bool {
+        self.touches(&Box3 { lo: p, hi: p })
+    }
+
+    pub fn diagonal(&self) -> f64 {
+        norm(sub(self.hi, self.lo))
+    }
+}
diff --git a/crates/mesh/src/surface/tests.rs b/crates/mesh/src/surface/tests.rs
new file mode 100644
index 0000000..8203a97
--- /dev/null
+++ b/crates/mesh/src/surface/tests.rs
@@ -0,0 +1,394 @@
+use std::f64::consts::PI;
+
+use super::intersect::{segment_meets_triangle, surface_pair_intersects, triangles_meet};
+use super::*;
+
+type Mesh = (Vec<[f64; 3]>, Vec<[u32; 3]>, Vec<u32>);
+
+/// An axis-aligned cube, counter-clockwise from outside, one face tag per side.
+fn cube(lo: [f64; 3], size: f64) -> Mesh {
+    let p =
+        (0..8).map(|i| [lo[0] + size * f64::from(i & 1), lo[1] + size * f64::from((i >> 1) & 1), lo[2] + size * f64::from((i >> 2) & 1)]).collect();
+    let t = vec![[0, 2, 3], [0, 3, 1], [4, 5, 7], [4, 7, 6], [0, 1, 5], [0, 5, 4], [2, 6, 7], [2, 7, 3], [0, 4, 6], [0, 6, 2], [1, 3, 7], [1, 7, 5]];
+    (p, t, vec![0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5])
+}
+
+fn flip(m: &mut Mesh) {
+    for t in &mut m.1 {
+        t.swap(1, 2);
+    }
+}
+
+fn merge(a: Mesh, b: Mesh) -> Mesh {
+    let (mut p, mut t, mut f) = a;
+    let n = p.len() as u32;
+    let faces = f.iter().max().map_or(0, |m| m + 1);
+    p.extend(b.0);
+    t.extend(b.1.iter().map(|t| t.map(|v| v + n)));
+    f.extend(b.2.iter().map(|x| x + faces));
+    (p, t, f)
+}
+
+/// Every triangle with its own three vertices, as SolidCraft's tessellator emits them.
+fn soup(m: &Mesh) -> Mesh {
+    let mut p = Vec::new();
+    let mut t = Vec::new();
+    for tri in &m.1 {
+        let base = p.len() as u32;
+        p.extend(tri.map(|v| m.0[v as usize]));
+        t.push([base, base + 1, base + 2]);
+    }
+    (p, t, m.2.clone())
+}
+
+/// A closed parametric surface on an `n × m` grid, (u, v) ↦ point, periodic in u; `poles` caps
+/// v = 0 and v = 1 with single vertices (a sphere), otherwise v is periodic too (a torus).
+fn grid(n: u32, m: u32, poles: bool, f: impl Fn(f64, f64) -> [f64; 3]) -> Mesh {
+    let mut p = Vec::new();
+    let rows = if poles { m - 1 } else { m };
+    for j in 0..rows {
+        for i in 0..n {
+            let v = if poles { f64::from(j + 1) / f64::from(m) } else { f64::from(j) / f64::from(m) };
+            p.push(f(f64::from(i) / f64::from(n), v));
+        }
+    }
+    let id = |i: u32, j: u32| (j % rows) * n + i % n;
+    let mut t = Vec::new();
+    let quads = if poles { rows - 1 } else { rows };
+    for j in 0..quads {
+        for i in 0..n {
+            let (a, b, c, d) = (id(i, j), id(i + 1, j), id(i + 1, j + 1), id(i, j + 1));
+            t.push([a, b, c]);
+            t.push([a, c, d]);
+        }
+    }
+    if poles {
+        let (south, north) = (p.len() as u32, p.len() as u32 + 1);
+        p.push(f(0.0, 0.0));
+        p.push(f(0.0, 1.0));
+        for i in 0..n {
+            t.push([south, id(i + 1, 0), id(i, 0)]);
+            t.push([north, id(i, rows - 1), id(i + 1, rows - 1)]);
+        }
+    }
+    let faces = vec![0; t.len()];
+    let mut mesh = (p, t, faces);
+    // Orient outwards.
+    let six_v: f64 = mesh.1.iter().map(|t| crate::vec3::det(mesh.0[t[0] as usize], mesh.0[t[1] as usize], mesh.0[t[2] as usize])).sum();
+    if six_v < 0.0 {
+        flip(&mut mesh);
+    }
+    mesh
+}
+
+fn torus(n: u32, m: u32) -> Mesh {
+    grid(n, m, false, |u, v| {
+        let (a, b) = (2.0 * PI * u, 2.0 * PI * v);
+        [(3.0 + b.cos()) * a.cos(), (3.0 + b.cos()) * a.sin(), b.sin()]
+    })
+}
+
+fn sphere(n: u32, m: u32) -> Mesh {
+    grid(n, m, true, |u, v| {
+        let (a, b) = (2.0 * PI * u, PI * v);
+        [b.sin() * a.cos(), b.sin() * a.sin(), -b.cos()]
+    })
+}
+
+fn moved(m: Mesh, d: [f64; 3]) -> Mesh {
+    (m.0.iter().map(|p| [p[0] + d[0], p[1] + d[1], p[2] + d[2]]).collect(), m.1, m.2)
+}
+
+fn prepare(m: &Mesh) -> Result<BoundarySurface> {
+    BoundarySurface::new(&m.0, &m.1, &m.2, &SurfaceOptions::default())
+}
+
+fn defects(m: &Mesh) -> Defects {
+    match prepare(m) {
+        Err(Error::Surface(d)) => d,
+        other => panic!("expected defects, got {other:?}"),
+    }
+}
+
+fn kinds(d: &Defects) -> Vec<&'static str> {
+    let mut k: Vec<_> = d.listed.iter().map(Defect::kind).collect();
+    k.dedup();
+    k
+}
+
+#[test]
+fn cube_is_a_valid_boundary_with_its_edges_as_features() {
+    let s = prepare(&cube([0.0; 3], 2.0)).unwrap();
+    assert_eq!((s.points().len(), s.triangles().len()), (8, 12));
+    assert_eq!(s.shells().len(), 1);
+    let sh = s.shells()[0];
+    assert_eq!((sh.genus, sh.triangles, sh.vertices, sh.void), (0, 12, 8, false));
+    assert!((s.volume() - 8.0).abs() < 1e-12);
+    assert_eq!(s.feature_edges().len(), 12);
+    assert!(s.feature_edges().iter().all(|e| e.kind == FeatureKind::FaceBoundary && (e.angle - 90.0).abs() < 1e-9));
+    assert_eq!(s.feature_corners(), (0..8).collect::<Vec<u32>>());
+    for (t, n) in s.neighbors().iter().enumerate() {
+        for (k, &m) in n.iter().enumerate() {
+            let (a, b) = (s.triangles()[t][k], s.triangles()[t][(k + 1) % 3]);
+            let other = s.triangles()[m as usize];
+            assert!((0..3).any(|j| other[j] == b && other[(j + 1) % 3] == a), "triangle {t} edge {k}");
+        }
+    }
+    // Feature edges name first the triangle that traverses them low → high.
+    for e in s.feature_edges() {
+        let t = s.triangles()[e.triangles[0] as usize];
+        assert!((0..3).any(|j| t[j] == e.vertices[0] && t[(j + 1) % 3] == e.vertices[1]));
+    }
+}
+
+#[test]
+fn a_triangle_soup_welds_into_the_solid() {
+    let c = cube([1.0, -2.0, 0.5], 3.0);
+    let s = prepare(&soup(&c)).unwrap();
+    assert_eq!((s.points().len(), s.triangles().len()), (8, 12));
+    assert_eq!(s.source_vertex().len(), 8);
+    assert_eq!(s.source_triangle(), (0..12).collect::<Vec<u32>>());
+    assert_eq!(s.tri_face(), c.2.as_slice());
+    assert!((s.volume() - 27.0).abs() < 1e-9);
+
+    // Tessellation noise below the tolerance welds; with exact welding it leaves holes.
+    let mut noisy = soup(&c);
+    for (i, p) in noisy.0.iter_mut().enumerate() {
+        p[0] += 1e-12 * (i % 3) as f64;
+    }
+    assert_eq!(prepare(&noisy).unwrap().points().len(), 8);
+    let exact = SurfaceOptions { weld_tolerance: 0.0, ..SurfaceOptions::default() };
+    assert!(matches!(BoundarySurface::new(&noisy.0, &noisy.1, &noisy.2, &exact), Err(Error::Surface(_))));
+    // Or an absolute distance (a CAD kernel's tolerance).
+    let mut noisy = soup(&c);
+    for (i, p) in noisy.0.iter_mut().enumerate() {
+        p[1] += 1.5e-6 * (i % 2) as f64;
+    }
+    assert!(matches!(prepare(&noisy), Err(Error::Surface(_))));
+    let kernel = SurfaceOptions { weld_distance: 2e-6, ..SurfaceOptions::default() };
+    assert_eq!(BoundarySurface::new(&noisy.0, &noisy.1, &noisy.2, &kernel).unwrap().points().len(), 8);
+}
+
+#[test]
+fn untagged_faces_use_the_dihedral_angle() {
+    let c = cube([0.0; 3], 1.0);
+    let s = BoundarySurface::new(&c.0, &c.1, &[], &SurfaceOptions::default()).unwrap();
+    assert!(s.tri_face().iter().all(|&f| f == 0));
+    assert_eq!(s.feature_edges().len(), 12);
+    assert!(s.feature_edges().iter().all(|e| e.kind == FeatureKind::Sharp));
+    let blunt = SurfaceOptions { feature_angle: 90.0, ..SurfaceOptions::default() };
+    assert!(BoundarySurface::new(&c.0, &c.1, &[], &blunt).unwrap().feature_edges().is_empty());
+}
+
+#[test]
+fn torus_and_sphere_have_their_genus() {
+    let t = prepare(&torus(16, 8)).unwrap();
+    assert_eq!((t.shells().len(), t.shells()[0].genus), (1, 1));
+    let s = prepare(&sphere(24, 12)).unwrap();
+    assert_eq!((s.shells().len(), s.shells()[0].genus), (1, 0));
+    // Inscribed in the unit sphere: a little less than its volume.
+    let ball = 4.0 / 3.0 * PI;
+    assert!(s.volume() < ball && s.volume() > 0.9 * ball, "{}", s.volume());
+    // A smooth surface tagged as one face has no features at a fine tessellation.
+    assert!(prepare(&torus(64, 32)).unwrap().feature_edges().is_empty());
+}
+
+#[test]
+fn a_void_is_a_negative_shell_inside_another() {
+    let mut inner = cube([1.0; 3], 1.0);
+    flip(&mut inner);
+    let s = prepare(&merge(cube([0.0; 3], 3.0), inner)).unwrap();
+    assert_eq!(s.shells().len(), 2);
+    assert!(!s.shells()[0].void && s.shells()[1].void);
+    assert!((s.volume() - 26.0).abs() < 1e-9);
+    assert!(s.tri_shell()[..12].iter().all(|&x| x == 0) && s.tri_shell()[12..].iter().all(|&x| x == 1));
+
+    // Two separate solids are two outer shells.
+    let two = prepare(&merge(cube([0.0; 3], 1.0), cube([5.0, 0.0, 0.0], 1.0))).unwrap();
+    assert!(two.shells().iter().all(|s| !s.void && s.volume > 0.0));
+
+    // A void facing outwards, and a solid facing inwards, are inverted.
+    let d = defects(&merge(cube([0.0; 3], 3.0), cube([1.0; 3], 1.0)));
+    assert_eq!(d.listed, vec![Defect::Inverted { shell: 1, triangle: 12, volume: 1.0 }]);
+    let mut c = cube([0.0; 3], 1.0);
+    flip(&mut c);
+    assert_eq!(kinds(&defects(&c)), vec!["inverted"]);
+}
+
+#[test]
+fn holes_and_bad_edges_are_located() {
+    let mut open = cube([0.0; 3], 1.0);
+    open.1.remove(3);
+    open.2.remove(3);
+    let d = defects(&open);
+    assert_eq!(d.total, 3);
+    assert!(d.listed.iter().all(|x| x.kind() == "open_edge" && x.at().is_some()));
+    assert!(d.to_string().starts_with("3 defects; open edge"), "{d}");
+
+    let mut flipped = cube([0.0; 3], 1.0);
+    flipped.1[5].swap(0, 1);
+    let d = defects(&flipped);
+    assert_eq!(kinds(&d), vec!["flipped"]);
+    assert_eq!(d.total, 3);
+    assert!(d.listed.iter().all(|x| matches!(x, Defect::Flipped { triangles, .. } if triangles.contains(&5))));
+
+    // A fin on an edge: three triangles there.
+    let mut fin = cube([0.0; 3], 1.0);
+    fin.0.push([0.5, -1.0, -1.0]);
+    fin.1.push([0, 1, 8]);
+    fin.2.push(9);
+    assert!(kinds(&defects(&fin)).contains(&"non_manifold_edge"));
+}
+
+#[test]
+fn cubes_touching_at_a_corner_are_not_a_manifold() {
+    let d = defects(&merge(cube([0.0; 3], 1.0), cube([1.0; 3], 1.0)));
+    assert_eq!(d.listed, vec![Defect::NonManifoldVertex { vertex: 7, at: [1.0; 3] }]);
+}
+
+#[test]
+fn intersecting_and_touching_solids_are_self_intersections() {
+    let d = defects(&merge(cube([0.0; 3], 2.0), cube([1.0; 3], 2.0)));
+    assert_eq!(kinds(&d), vec!["self_intersection"]);
+    // A smaller block resting on a face: coplanar contact without shared vertices.
+    let d = defects(&merge(cube([0.0; 3], 2.0), cube([2.0, 0.5, 0.5], 1.0)));
+    assert_eq!(kinds(&d), vec!["self_intersection"]);
+    // A void whose corner is pulled onto the outer surface (onto the bottom face's diagonal).
+    let mut inner = cube([1.0; 3], 1.0);
+    inner.0[0] = [1.0, 1.0, 0.0];
+    flip(&mut inner);
+    let d = defects(&merge(cube([0.0; 3], 3.0), inner));
+    assert_eq!(kinds(&d), vec!["self_intersection"]);
+}
+
+#[test]
+fn many_intersections_are_capped() {
+    // Small cubes all crossing the top of a large one.
+    let mut m = cube([0.0; 3], 100.0);
+    for i in 0..60 {
+        m = merge(m, cube([1.5 * f64::from(i) + 0.1, 0.1, 99.5], 1.0));
+    }
+    let d = defects(&m);
+    assert_eq!(d.listed.len(), MAX_LISTED_DEFECTS);
+    assert!(d.listed.windows(2).all(|w| match (&w[0], &w[1]) {
+        (Defect::SelfIntersection { triangles: a, .. }, Defect::SelfIntersection { triangles: b, .. }) => a <= b,
+        _ => false,
+    }));
+}
+
+#[test]
+fn hostile_input_is_an_error() {
+    let c = cube([0.0; 3], 1.0);
+    let mut nan = c.clone();
+    nan.0[3][1] = f64::NAN;
+    assert_eq!(defects(&nan).listed, vec![Defect::BadCoordinate { vertex: 3 }]);
+    let mut far = c.clone();
+    far.0[0][0] = 1e10;
+    assert_eq!(kinds(&defects(&far)), vec!["bad_coordinate"]);
+    let mut bad = c.clone();
+    bad.1[2][1] = 99;
+    assert_eq!(defects(&bad).listed, vec![Defect::BadIndex { triangle: 2 }]);
+    let mut degenerate = c.clone();
+    degenerate.0.push([2.0, 0.0, 0.0]);
+    degenerate.1.push([0, 1, 8]);
+    degenerate.2.push(0);
+    assert_eq!(kinds(&defects(&degenerate)), vec!["degenerate"]);
+
+    let o = SurfaceOptions::default();
+    assert!(matches!(BoundarySurface::new(&[], &[], &[], &o), Err(Error::Invalid(_))));
+    assert!(matches!(BoundarySurface::new(&c.0, &c.1, &[0, 1], &o), Err(Error::Invalid(_))));
+    assert!(matches!(BoundarySurface::new(&[[1.0; 3]; 3], &[[0, 1, 2]], &[], &o), Err(Error::Invalid(_))));
+    for w in [-1.0, f64::NAN, 1e-20, 0.5, f64::INFINITY] {
+        assert!(BoundarySurface::new(&c.0, &c.1, &c.2, &SurfaceOptions { weld_tolerance: w, ..o }).is_err(), "{w}");
+    }
+    for d in [-1.0, f64::NAN, f64::INFINITY, 0.01] {
+        assert!(BoundarySurface::new(&c.0, &c.1, &c.2, &SurfaceOptions { weld_distance: d, ..o }).is_err(), "{d}");
+    }
+    for a in [0.0, -1.0, 181.0, f64::NAN] {
+        assert!(BoundarySurface::new(&c.0, &c.1, &c.2, &SurfaceOptions { feature_angle: a, ..o }).is_err(), "{a}");
+    }
+}
+
+#[test]
+fn tessellation_noise_snaps_into_the_predicates_domain() {
+    let mut c = cube([0.0; 3], 1.0);
+    c.0[0] = [1e-40, -1e-300, 0.0];
+    let s = BoundarySurface::new(&c.0, &c.1, &c.2, &SurfaceOptions { weld_tolerance: 0.0, ..SurfaceOptions::default() }).unwrap();
+    assert_eq!(s.points()[0], [0.0; 3]);
+    assert!((0..8).all(|v| s.point3(v).is_some()));
+}
+
+#[test]
+fn output_is_deterministic() {
+    let m = merge(merge(soup(&torus(20, 10)), moved(sphere(16, 8), [0.0, 0.0, 10.0])), cube([20.0; 3], 1.0));
+    let a = prepare(&m).unwrap();
+    assert_eq!(a, prepare(&m).unwrap());
+    assert_eq!(a.shells().iter().map(|s| s.genus).collect::<Vec<_>>(), vec![1, 0, 0]);
+}
+
+#[test]
+fn random_soups_never_panic() {
+    let mut seed = 0x2545_f491_4f6c_dd1d_u64;
+    let mut next = || {
+        seed ^= seed << 13;
+        seed ^= seed >> 7;
+        seed ^= seed << 17;
+        seed
+    };
+    let specials = [0.0, -0.0, 1.0, -1.0, 1e-40, 1e9, f64::NAN, f64::INFINITY, 0.5];
+    for _ in 0..2000 {
+        let n = (next() % 9 + 1) as usize;
+        let p: Vec<[f64; 3]> =
+            (0..n).map(|_| [0, 1, 2].map(|_| if next() % 4 == 0 { specials[(next() % 9) as usize] } else { (next() % 5) as f64 })).collect();
+        let t: Vec<[u32; 3]> = (0..next() % 12).map(|_| [0, 1, 2].map(|_| (next() % (n as u64 + 1)) as u32)).collect();
+        let _ = BoundarySurface::new(&p, &t, &[], &SurfaceOptions::default());
+    }
+}
+
+fn p3(c: [f64; 3]) -> Point3 {
+    Point3::try_from(c).unwrap()
+}
+
+#[test]
+fn exact_pair_tests() {
+    let [o, x, y, z] = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]].map(p3);
+    let t = [&o, &x, &y];
+    // Segments: piercing, touching a vertex, coplanar crossing, coplanar outside, parallel.
+    assert!(segment_meets_triangle(&p3([0.2, 0.2, -1.0]), &p3([0.2, 0.2, 1.0]), t));
+    assert!(segment_meets_triangle(&p3([1.0, 0.0, 0.0]), &p3([1.0, 0.0, 5.0]), t));
+    assert!(segment_meets_triangle(&p3([-1.0, 0.3, 0.0]), &p3([2.0, 0.3, 0.0]), t));
+    assert!(!segment_meets_triangle(&p3([0.6, 0.6, 0.0]), &p3([2.0, 2.0, 0.0]), t));
+    assert!(!segment_meets_triangle(&p3([0.0, 0.0, 1.0]), &p3([1.0, 0.0, 1.0]), t));
+    // On the hypotenuse x + y = 1 of a triangle whose corners are not exact in binary.
+    let (a, b) = (p3([0.1, 0.9, 0.0]), p3([0.9, 0.1, 0.0]));
+    let mid = p3([0.5, 0.5, 0.0]);
+    assert!(
+        segment_meets_triangle(&p3([0.5, 0.5, -1.0]), &p3([0.5, 0.5, 1.0]), [&o, &a, &b])
+            == (orient2d_sign(&a, &b, &mid) * orient2d_sign(&a, &b, &o) >= 0)
+    );
+
+    assert!(triangles_meet(t, [&p3([0.1, 0.1, -1.0]), &p3([0.1, 0.1, 1.0]), &p3([0.5, 0.1, 0.0])]));
+    assert!(!triangles_meet(t, [&p3([0.1, 0.1, 0.5]), &p3([0.6, 0.1, 0.5]), &p3([0.1, 0.6, 0.5])]));
+
+    // Shared edge o–x: folded flat onto t, flat on the other side, and bent out of plane.
+    let ids = |s: [u32; 3], t: [u32; 3]| [s, t];
+    let (w, w2, w3) = (p3([0.5, 0.5, 0.0]), p3([0.5, -0.5, 0.0]), p3([0.5, 0.5, 1.0]));
+    assert!(surface_pair_intersects(ids([0, 1, 2], [1, 0, 3]), t, [&x, &o, &w]));
+    assert!(!surface_pair_intersects(ids([0, 1, 2], [1, 0, 3]), t, [&x, &o, &w2]));
+    assert!(!surface_pair_intersects(ids([0, 1, 2], [1, 0, 3]), t, [&x, &o, &w3]));
+
+    // Shared vertex o: a coplanar triangle inside t's corner, one outside it, one crossing t, one
+    // out of plane touching only at o, and one whose edge runs along t's edge o–x.
+    let at_o = ids([0, 1, 2], [0, 3, 4]);
+    assert!(surface_pair_intersects(at_o, t, [&o, &p3([0.3, 0.1, 0.0]), &p3([0.1, 0.3, 0.0])]));
+    assert!(!surface_pair_intersects(at_o, t, [&o, &p3([-0.3, -0.1, 0.0]), &p3([-0.1, -0.3, 0.0])]));
+    assert!(surface_pair_intersects(at_o, t, [&o, &p3([0.3, 0.3, -1.0]), &p3([0.3, 0.3, 1.0])]));
+    assert!(!surface_pair_intersects(at_o, t, [&o, &z, &p3([-1.0, -1.0, 1.0])]));
+    assert!(surface_pair_intersects(at_o, t, [&o, &p3([0.5, 0.0, 0.0]), &p3([0.5, -1.0, 0.0])]));
+}
+
+/// The exact orientation of `c` against the directed line `ab` seen from +z, as -1, 0, 1.
+fn orient2d_sign(a: &Point3, b: &Point3, c: &Point3) -> i32 {
+    let p = |q: &Point3| crate::predicates::Point2::new(q.coords()[0], q.coords()[1]).unwrap();
+    crate::predicates::orient2d(&p(a), &p(b), &p(c)) as i32
+}
diff --git a/crates/mesh/src/surface/topology.rs b/crates/mesh/src/surface/topology.rs
new file mode 100644
index 0000000..a8447cf
--- /dev/null
+++ b/crates/mesh/src/surface/topology.rs
@@ -0,0 +1,262 @@
+//! Edges, manifoldness, orientation consistency, shells, nesting and feature edges.
+//!
+//! Edges come from sorting half-edges (no hash order); shells are numbered by their first
+//! triangle; every loop is bounded by the input size.
+
+use std::f64::consts::PI;
+
+use super::{Box3, Defect, Defects, FeatureEdge, FeatureKind, MAX_WINDING_WORK, Shell, Surface, index};
+use crate::vec3::{angle_deg, det, dot, norm, sub};
+use crate::{Error, Result};
+
+/// No triangle (only while building).
+const NONE: u32 = u32::MAX;
+
+/// An interior edge: its vertices (ascending) and its two triangles, the one traversing it from
+/// `v[0]` to `v[1]` first.
+pub(crate) struct Edge {
+    pub v: [u32; 2],
+    pub t: [u32; 2],
+}
+
+/// One triangle's use of an edge.
+struct Half {
+    lo: u32,
+    hi: u32,
+    /// Traversed from `lo` to `hi`.
+    forward: bool,
+    tri: u32,
+    corner: usize,
+}
+
+pub(crate) struct Topology {
+    /// Per triangle `[a, b, c]`: the triangles across `ab`, `bc`, `ca`.
+    pub neighbors: Vec<[u32; 3]>,
+    pub edges: Vec<Edge>,
+}
+
+impl Topology {
+    /// Pair every edge with exactly one oppositely oriented edge and check every vertex is a
+    /// single fan.
+    pub fn build(s: &Surface) -> Result<Topology> {
+        let mut half: Vec<Half> = Vec::with_capacity(s.triangles.len() * 3);
+        for (t, tri) in s.triangles.iter().enumerate() {
+            for (k, (a, b)) in [(tri[0], tri[1]), (tri[1], tri[2]), (tri[2], tri[0])].into_iter().enumerate() {
+                half.push(Half { lo: a.min(b), hi: a.max(b), forward: a < b, tri: index(t), corner: k });
+            }
+        }
+        half.sort_unstable_by_key(|h| (h.lo, h.hi, h.tri, h.corner));
+
+        let mut neighbors = vec![[NONE; 3]; s.triangles.len()];
+        let mut edges = Vec::with_capacity(half.len() / 2);
+        let mut defects = Defects::new();
+        for run in half.chunk_by(|x, y| (x.lo, x.hi) == (y.lo, y.hi)) {
+            let Some(&Half { lo, hi, tri: t0, .. }) = run.first() else { continue };
+            let vertices = || [s.source_vert(lo), s.source_vert(hi)];
+            let at = || s.midpoint(lo, hi);
+            match run {
+                [_] => defects.push(Defect::OpenEdge { vertices: vertices(), triangle: s.source_tri(t0 as usize), at: at() }),
+                [x, y] if x.forward == y.forward => defects.push(Defect::Flipped {
+                    vertices: vertices(),
+                    triangles: [s.source_tri(x.tri as usize), s.source_tri(y.tri as usize)],
+                    at: at(),
+                }),
+                [x, y] => {
+                    for (from, to) in [(x, y), (y, x)] {
+                        if let Some(slot) = neighbors.get_mut(from.tri as usize).and_then(|n| n.get_mut(from.corner)) {
+                            *slot = to.tri;
+                        }
+                    }
+                    let t = if x.forward { [x.tri, y.tri] } else { [y.tri, x.tri] };
+                    edges.push(Edge { v: [lo, hi], t });
+                }
+                _ => defects.push(Defect::NonManifoldEdge { vertices: vertices(), triangles: index(run.len()), at: at() }),
+            }
+        }
+        defects.into_result()?;
+
+        let topo = Topology { neighbors, edges };
+        topo.check_fans(s)?;
+        Ok(topo)
+    }
+
+    /// Every vertex's triangles form one fan: walking across the edges around it from any of its
+    /// triangles visits all of them.
+    fn check_fans(&self, s: &Surface) -> Result<()> {
+        let mut degree = vec![0u32; s.pts.len()];
+        for v in s.triangles.iter().flatten() {
+            if let Some(d) = degree.get_mut(*v as usize) {
+                *d += 1;
+            }
+        }
+        let mut seen = vec![false; s.pts.len()];
+        let mut defects = Defects::new();
+        for (t0, tri) in s.triangles.iter().enumerate() {
+            for &v in tri {
+                let Some(flag) = seen.get_mut(v as usize) else { continue };
+                if *flag {
+                    continue;
+                }
+                *flag = true;
+                let deg = degree.get(v as usize).copied().unwrap_or(0);
+                // Rotate across the edge leaving `v` in each triangle until back at the start.
+                let (mut t, mut steps) = (index(t0), 0u32);
+                loop {
+                    let corner = s.triangles.get(t as usize).and_then(|tr| tr.iter().position(|&w| w == v));
+                    let next = corner.and_then(|k| self.neighbors.get(t as usize).and_then(|n| n.get(k))).copied().unwrap_or(NONE);
+                    steps += 1;
+                    if next == NONE || next == index(t0) || steps > deg {
+                        break;
+                    }
+                    t = next;
+                }
+                if steps != deg {
+                    defects.push(Defect::NonManifoldVertex { vertex: s.source_vert(v), at: s.xyz(v) });
+                }
+            }
+        }
+        defects.into_result()
+    }
+}
+
+pub(crate) struct Shells {
+    pub tri_shell: Vec<u32>,
+    pub shells: Vec<Shell>,
+}
+
+/// Connected components, their genus, signed volume and nesting; fails when a shell faces the
+/// wrong way.
+pub(crate) fn shells(s: &Surface, topo: &Topology) -> Result<Shells> {
+    let n = s.triangles.len();
+    let mut tri_shell = vec![NONE; n];
+    let mut shells = Vec::new();
+    // Per shell: its box and its reference point (first vertex of its first triangle).
+    let mut extent: Vec<(Box3, [f64; 3])> = Vec::new();
+    let mut vertex_shell = vec![NONE; s.pts.len()];
+    let mut stack = Vec::new();
+    for start in 0..n {
+        if tri_shell.get(start) != Some(&NONE) {
+            continue;
+        }
+        let id = index(shells.len());
+        let origin = s.triangles.get(start).map_or([0.0; 3], |t| s.xyz(t[0]));
+        let (mut faces, mut verts, mut six_v) = (0u32, 0u32, 0.0);
+        let mut bounds = Box3::EMPTY;
+        if let Some(x) = tri_shell.get_mut(start) {
+            *x = id;
+        }
+        stack.push(start);
+        while let Some(t) = stack.pop() {
+            faces += 1;
+            let tri = s.triangles.get(t).copied().unwrap_or([0; 3]);
+            for v in tri {
+                if let Some(vs) = vertex_shell.get_mut(v as usize)
+                    && *vs != id
+                {
+                    *vs = id;
+                    verts += 1;
+                    bounds.add_point(s.xyz(v));
+                }
+            }
+            let [a, b, c] = s.tri_xyz(t).map(|p| sub(p, origin));
+            six_v += det(a, b, c);
+            for &m in topo.neighbors.get(t).into_iter().flatten() {
+                if let Some(x) = tri_shell.get_mut(m as usize)
+                    && *x == NONE
+                {
+                    *x = id;
+                    stack.push(m as usize);
+                }
+            }
+        }
+        // Closed 2-manifold: E = 3F/2, χ = V − E + F = 2 − 2g.
+        let chi = i64::from(verts) - i64::from(faces) * 3 / 2 + i64::from(faces);
+        let genus = u32::try_from((2 - chi).max(0) / 2).unwrap_or(u32::MAX);
+        shells.push(Shell { first_triangle: index(start), triangles: faces, vertices: verts, genus, volume: six_v / 6.0, void: false });
+        extent.push((bounds, origin));
+    }
+
+    nest(s, &tri_shell, &mut shells, &extent)?;
+
+    let mut defects = Defects::new();
+    for (id, sh) in shells.iter().enumerate() {
+        let d = extent.get(id).map_or(0.0, |(b, _)| b.diagonal());
+        let triangle = s.source_tri(sh.first_triangle as usize);
+        // Each term of the sum is at most d³ and carries a few ulps of rounding.
+        if !(sh.volume.abs() > 8.0 * f64::EPSILON * f64::from(sh.triangles) * d * d * d) {
+            defects.push(Defect::Flat { shell: index(id), triangle });
+        } else if (sh.volume > 0.0) == sh.void {
+            defects.push(Defect::Inverted { shell: index(id), triangle, volume: sh.volume });
+        }
+    }
+    defects.into_result()?;
+    Ok(Shells { tri_shell, shells })
+}
+
+/// Mark the voids: shells inside an odd number of other shells, by the generalized winding
+/// number of one of their vertices (shells do not intersect, so it is ±1 inside and 0 outside,
+/// up to rounding).
+fn nest(s: &Surface, tri_shell: &[u32], shells: &mut [Shell], extent: &[(Box3, [f64; 3])]) -> Result<()> {
+    if shells.len() < 2 {
+        return Ok(());
+    }
+    // Triangles of each shell, in order.
+    let mut members: Vec<Vec<usize>> = vec![Vec::new(); shells.len()];
+    for (t, &sh) in tri_shell.iter().enumerate() {
+        if let Some(m) = members.get_mut(sh as usize) {
+            m.push(t);
+        }
+    }
+    let mut work = 0u64;
+    for (i, &(_, x)) in extent.iter().enumerate() {
+        let mut depth = 0u32;
+        for (j, (bounds, _)) in extent.iter().enumerate() {
+            if j == i || !bounds.contains(x) {
+                continue;
+            }
+            let tris = members.get(j).map_or(&[][..], Vec::as_slice);
+            work += tris.len() as u64;
+            if work > MAX_WINDING_WORK {
+                return Err(Error::Limit(format!("shell nesting: more than {MAX_WINDING_WORK} triangle solid angles")));
+            }
+            let omega: f64 = tris.iter().map(|&t| solid_angle(s, t, x)).sum();
+            if (omega / (4.0 * PI)).abs() > 0.5 {
+                depth += 1;
+            }
+        }
+        if let Some(sh) = shells.get_mut(i) {
+            sh.void = depth % 2 == 1;
+        }
+    }
+    Ok(())
+}
+
+/// The signed solid angle triangle `t` subtends at `x` (van Oosterom and Strackee 1983).
+fn solid_angle(s: &Surface, t: usize, x: [f64; 3]) -> f64 {
+    let [a, b, c] = s.tri_xyz(t).map(|p| sub(p, x));
+    let (la, lb, lc) = (norm(a), norm(b), norm(c));
+    let den = la * lb * lc + dot(a, b) * lc + dot(a, c) * lb + dot(b, c) * la;
+    2.0 * det(a, b, c).atan2(den)
+}
+
+/// Edges between different B-rep faces, and edges inside a face whose normals turn by more than
+/// `feature_angle` degrees.
+pub(crate) fn feature_edges(s: &Surface, topo: &Topology, tri_face: &[u32], feature_angle: f64) -> Vec<FeatureEdge> {
+    let face = |t: u32| tri_face.get(t as usize).copied().unwrap_or(0);
+    let normals: Vec<[f64; 3]> = (0..s.triangles.len()).map(|t| s.normal(t)).collect();
+    let normal = |t: u32| normals.get(t as usize).copied().unwrap_or([0.0; 3]);
+    topo.edges
+        .iter()
+        .filter_map(|e| {
+            let angle = angle_deg(normal(e.t[0]), normal(e.t[1]));
+            let kind = if face(e.t[0]) != face(e.t[1]) {
+                FeatureKind::FaceBoundary
+            } else if angle > feature_angle {
+                FeatureKind::Sharp
+            } else {
+                return None;
+            };
+            Some(FeatureEdge { vertices: e.v, triangles: e.t, angle, kind })
+        })
+        .collect()
+}
diff --git a/crates/mesh/src/surface/weld.rs b/crates/mesh/src/surface/weld.rs
new file mode 100644
index 0000000..b899d64
--- /dev/null
+++ b/crates/mesh/src/surface/weld.rs
@@ -0,0 +1,138 @@
+//! Vertex welding in input order: each vertex joins the earliest representative within the
+//! tolerance, or becomes one. Lookups go through a hash map, but nothing is ever iterated out of
+//! it, so the result does not depend on hashing.
+
+use std::collections::HashMap;
+
+use super::{Box3, MAX_RELATIVE_WELD, index};
+use crate::predicates::MIN_MAGNITUDE;
+use crate::vec3::{dot, sub};
+use crate::{Error, Result};
+
+pub(crate) struct Welded {
+    pub points: Vec<[f64; 3]>,
+    pub triangles: Vec<[u32; 3]>,
+    pub source_vertex: Vec<u32>,
+    pub source_triangle: Vec<u32>,
+}
+
+/// Weld the vertices `triangles` use (indices already checked) within the larger of `relative`
+/// × the bounding-box diagonal and `absolute`, snap coordinates below [`MIN_MAGNITUDE`] to zero
+/// and drop triangles that collapse. A distance above [`MAX_RELATIVE_WELD`] of the diagonal is an
+/// error; one below 1e-15 of it welds identical coordinates only.
+pub(crate) fn weld(positions: &[[f64; 3]], triangles: &[[u32; 3]], relative: f64, absolute: f64) -> Result<Welded> {
+    let mut used = vec![false; positions.len()];
+    for v in triangles.iter().flatten() {
+        if let Some(u) = used.get_mut(*v as usize) {
+            *u = true;
+        }
+    }
+    let referenced = || positions.iter().zip(&used).enumerate().filter(|(_, (_, u))| **u).map(|(v, (p, _))| (v, *p));
+    let mut bounds = Box3::EMPTY;
+    for (_, p) in referenced() {
+        bounds.add_point(p);
+    }
+    let diagonal = bounds.diagonal();
+    let tol = (relative * diagonal).max(absolute);
+    if tol > MAX_RELATIVE_WELD * diagonal {
+        return Err(Error::Invalid(format!(
+            "weld distance {tol} is more than {MAX_RELATIVE_WELD} of the surface's size ({diagonal}): it would merge distinct vertices"
+        )));
+    }
+
+    let mut lookup = if tol >= 1e-15 * diagonal && tol > 0.0 {
+        Lookup::Grid { tol, lo: bounds.lo, heads: HashMap::new(), next: Vec::new() }
+    } else {
+        Lookup::Exact(HashMap::new())
+    };
+    let mut map = vec![u32::MAX; positions.len()];
+    let mut points: Vec<[f64; 3]> = Vec::new();
+    let mut source_vertex = Vec::new();
+    for (v, p) in referenced() {
+        let w = match lookup.find(p, &points) {
+            Some(w) => w,
+            None => {
+                let w = index(points.len());
+                points.push(p);
+                source_vertex.push(index(v));
+                lookup.insert(p, w);
+                w
+            }
+        };
+        if let Some(m) = map.get_mut(v) {
+            *m = w;
+        }
+    }
+    for x in points.iter_mut().flatten() {
+        if x.abs() < MIN_MAGNITUDE {
+            *x = 0.0;
+        }
+    }
+
+    let mut out = Vec::with_capacity(triangles.len());
+    let mut source_triangle = Vec::with_capacity(triangles.len());
+    for (t, tri) in triangles.iter().enumerate() {
+        let w = tri.map(|v| map.get(v as usize).copied().unwrap_or(u32::MAX));
+        if w[0] != w[1] && w[1] != w[2] && w[2] != w[0] {
+            out.push(w);
+            source_triangle.push(index(t));
+        }
+    }
+    Ok(Welded { points, triangles: out, source_vertex, source_triangle })
+}
+
+/// Representatives found so far.
+enum Lookup {
+    /// Identical coordinates only (keyed by bits; `+ 0.0` turns -0 into +0).
+    Exact(HashMap<[u64; 3], u32>),
+    /// Cells of side `tol`: a vertex's representatives are in its cell or the 26 around it.
+    /// Coordinates over the cell size stay below 1e15 (the distance is at least 1e-15 of the
+    /// diagonal), exact in i64. Each cell keeps the head of a chain through `next`.
+    Grid { tol: f64, lo: [f64; 3], heads: HashMap<[i64; 3], u32>, next: Vec<u32> },
+}
+
+impl Lookup {
+    fn cell(tol: f64, lo: [f64; 3], p: [f64; 3]) -> [i64; 3] {
+        [0, 1, 2].map(|k| ((p[k] - lo[k]) / tol).floor() as i64)
+    }
+
+    /// The earliest representative within the distance of `p`.
+    fn find(&self, p: [f64; 3], points: &[[f64; 3]]) -> Option<u32> {
+        match self {
+            Lookup::Exact(seen) => seen.get(&p.map(|x| (x + 0.0).to_bits())).copied(),
+            Lookup::Grid { tol, lo, heads, next } => {
+                let c = Self::cell(*tol, *lo, p);
+                let mut best: Option<u32> = None;
+                for dx in -1..=1_i64 {
+                    for dy in -1..=1_i64 {
+                        for dz in -1..=1_i64 {
+                            let key = [c[0].saturating_add(dx), c[1].saturating_add(dy), c[2].saturating_add(dz)];
+                            let mut w = heads.get(&key).copied().unwrap_or(u32::MAX);
+                            while let Some(&q) = points.get(w as usize) {
+                                let d = sub(p, q);
+                                if dot(d, d) <= tol * tol && best.is_none_or(|b| w < b) {
+                                    best = Some(w);
+                                }
+                                w = next.get(w as usize).copied().unwrap_or(u32::MAX);
+                            }
+                        }
+                    }
+                }
+                best
+            }
+        }
+    }
+
+    /// Record `p` as representative `w` (the next index).
+    fn insert(&mut self, p: [f64; 3], w: u32) {
+        match self {
+            Lookup::Exact(seen) => {
+                seen.insert(p.map(|x| (x + 0.0).to_bits()), w);
+            }
+            Lookup::Grid { tol, lo, heads, next } => {
+                let head = heads.insert(Self::cell(*tol, *lo, p), w);
+                next.push(head.unwrap_or(u32::MAX));
+            }
+        }
+    }
+}
diff --git a/crates/mesh/src/vec3.rs b/crates/mesh/src/vec3.rs
new file mode 100644
index 0000000..fe47db7
--- /dev/null
+++ b/crates/mesh/src/vec3.rs
@@ -0,0 +1,25 @@
+//! Small 3-vector helpers on `[f64; 3]` (floating point: measures and display, never meshing
+//! decisions, which use [`crate::predicates`]).
+
+pub(crate) type P = [f64; 3];
+
+pub(crate) fn sub(a: P, b: P) -> P {
+    [a[0] - b[0], a[1] - b[1], a[2] - b[2]]
+}
+pub(crate) fn dot(a: P, b: P) -> f64 {
+    a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
+}
+pub(crate) fn cross(a: P, b: P) -> P {
+    [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]]
+}
+pub(crate) fn norm(a: P) -> f64 {
+    dot(a, a).sqrt()
+}
+pub(crate) fn det(a: P, b: P, c: P) -> f64 {
+    dot(a, cross(b, c))
+}
+
+/// The angle between two vectors in degrees (0 when either is zero).
+pub(crate) fn angle_deg(a: P, b: P) -> f64 {
+    norm(cross(a, b)).atan2(dot(a, b)).to_degrees()
+}
diff --git a/crates/testkit/Cargo.toml b/crates/testkit/Cargo.toml
index fd9298d..495a49e 100644
--- a/crates/testkit/Cargo.toml
+++ b/crates/testkit/Cargo.toml
@@ -11,5 +11,6 @@ workspace = true
 
 [dependencies]
 modelcraft-geometry = { workspace = true, features = ["scripts"] }
+modelcraft-mesh = { workspace = true }
 serde_json = { workspace = true }
 thiserror = { workspace = true }
diff --git a/crates/testkit/src/runner.rs b/crates/testkit/src/runner.rs
index 1647efb..9a62684 100644
--- a/crates/testkit/src/runner.rs
+++ b/crates/testkit/src/runner.rs
@@ -14,6 +14,7 @@ use std::time::Instant;
 
 use crate::corpus::{self, CATEGORIES, Part};
 use modelcraft_geometry::{Design, Measure, SurfaceMesh};
+use modelcraft_mesh::surface::{BoundarySurface, SurfaceOptions};
 
 pub const FILE: &str = "docs/meshing/corpus.md";
 /// Chordal tolerance of the boundary surface, relative to the part's bounding-box diagonal.
@@ -40,7 +41,7 @@ pub trait Stage {
 
 /// The stages, in column order. Meshers add theirs here.
 fn stages() -> Vec<Box<dyn Stage>> {
-    vec![Box::new(SurfaceVolume)]
+    vec![Box::new(SurfaceVolume), Box::new(BoundaryCheck)]
 }
 
 /// The boundary surface's enclosed volume against SolidCraft's measured volume and the closed
@@ -60,6 +61,31 @@ impl Stage for SurfaceVolume {
     }
 }
 
+/// Boundary surface preparation (`modelcraft_mesh::surface`): welded, closed, manifold, outward
+/// and free of self-intersections, or the first defects found, which keep the part from being
+/// volume meshed.
+struct BoundaryCheck;
+
+impl Stage for BoundaryCheck {
+    fn headers(&self) -> &'static [&'static str] {
+        &["boundary", "genus", "feature edges", "check ms"]
+    }
+
+    fn run(&self, l: &Loaded) -> Result<Vec<String>, String> {
+        let start = Instant::now();
+        let (mut vertices, mut genus, mut features) = (0, Vec::new(), 0);
+        for (body, mesh) in l.design.bodies.iter().zip(&l.surfaces) {
+            let options = SurfaceOptions { weld_distance: mesh.weld_distance, ..SurfaceOptions::default() };
+            let s = BoundarySurface::new(&mesh.positions, &mesh.triangles, &mesh.tri_face, &options).map_err(|e| format!("{}: {e}", body.name))?;
+            vertices += s.points().len();
+            genus.extend(s.shells().iter().map(|sh| if sh.void { format!("{} (void)", sh.genus) } else { sh.genus.to_string() }));
+            features += s.feature_edges().len();
+        }
+        let ms = start.elapsed().as_secs_f64() * 1e3;
+        Ok(vec![format!("ok, {vertices} vertices"), genus.join(", "), features.to_string(), format!("{ms:.1}")])
+    }
+}
+
 fn load(part: &Part) -> Result<(Loaded<'_>, f64), String> {
     let start = Instant::now();
     let design = corpus::load(part).map_err(|e| e.to_string())?;
@@ -180,3 +206,22 @@ fn document(table: &str, parts: &[Part]) -> String {
     }
     s
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// The analytic shapes come out of boundary preparation valid, with their genus.
+    #[test]
+    fn analytic_parts_are_valid_boundaries() {
+        let expected =
+            [("analytic-sphere", "0"), ("analytic-torus", "1"), ("analytic-cylinder", "0"), ("analytic-pipe", "1"), ("analytic-plate-hole", "1")];
+        let parts = corpus::analytic_parts();
+        for (name, genus) in expected {
+            let part = parts.iter().find(|p| p.name == name).unwrap();
+            let (loaded, _) = load(part).unwrap();
+            let cells = BoundaryCheck.run(&loaded).unwrap();
+            assert_eq!(cells[1], genus, "{name}: {cells:?}");
+        }
+    }
+}
diff --git a/docs/architecture.md b/docs/architecture.md
index c72fe91..4af077f 100644
--- a/docs/architecture.md
+++ b/docs/architecture.md
@@ -51,8 +51,11 @@ Rules beyond the layer order:
    them), so a load on a named face follows it when the design is edited and opened again.
    Bodies of components are placed in the world through their occurrences (one body per
    placement), as SolidCraft displays and exports them.
-2. **Surface preparation.** Weld coincident vertices, check the surface is closed, manifold and
-   consistently oriented, find feature edges (face boundaries, sharp dihedral angles).
+2. **Surface preparation** (`modelcraft-mesh::surface`, `mesh.surface.check`). Weld the tessellation's per-face
+   vertex copies (at twice SolidCraft's linear tolerance), check the surface is closed, manifold, consistently and
+   outwardly oriented (voids inwards) and free of self-intersections (exact triangle–triangle tests), and find
+   feature edges (face boundaries, sharp dihedral angles). Defects come back located by triangle, vertex and
+   position.
 3. **Tetrahedral meshing.** Bowyer–Watson Delaunay on Shewchuk's adaptive predicates (symbolic
    perturbation for degeneracies), boundary recovery by the robust constrained Delaunay method of
    Diazzi et al. 2023 (implicit Steiner points on indirect predicates, gift-wrapping fallback), so
diff --git a/docs/meshing/corpus.md b/docs/meshing/corpus.md
index f0fe9f5..57e6bd7 100644
--- a/docs/meshing/corpus.md
+++ b/docs/meshing/corpus.md
@@ -8,29 +8,29 @@ clock on the machine that ran it; every other column is deterministic.
 
 ## Results
 
-| part | categories | bodies | faces | triangles | tol. mm | build ms | enclosed mm³ | SolidCraft mm³ | rel. diff | exact mm³ |
-|---|---|---|---|---|---|---|---|---|---|---|
-| bent-tube | sweepable, thin-wall, genus | 1 | 18 | 4296 | 0.121 | 443 | 4213.1 | 4159.4 | 1.3e-2 | 4159.8 |
-| bolted-assembly | assembly, blocky | 3 | 22 | 284 | 0.117 | 53 | 33004.7 | 32999.0 | 1.7e-4 | 32998.9 |
-| cross-cylinders | boolean, sharp | 1 | 10 | 1395 | 0.087 | 7912 | 32197.4 | 32357.8 | -5.0e-3 | 32365.8 |
-| filleted-bracket | fillet, genus, sweepable | 1 | 13 | 540 | 0.088 | 387 | 32162.3 | 32145.0 | 5.4e-4 | 32144.8 |
-| i-beam | sweepable, thin-wall | 1 | 14 | 44 | 0.208 | 3 | 61600.0 | 61600.0 | 0.0e0 | 61600.0 |
-| loft-transition | organic | 1 | 6 | 292 | 0.075 | 56 | 56390.2 | 56450.6 | -1.1e-3 | – |
-| perforated-plate | genus, thin-wall, sweepable | 1 | 54 | 1644 | 0.144 | 269 | 26822.2 | 26764.6 | 2.2e-3 | 26764.2 |
-| pinhole-block | small-feature, genus | 1 | 38 | 4044 | 0.073 | 203 | 49919.4 | 49909.6 | 2.0e-4 | 49909.5 |
-| rounded-box | fillet | 1 | 26 | 2540 | 0.054 | 148 | 23325.7 | 23340.5 | -6.3e-4 | 23340.8 |
-| shallow-groove | near-tangent, boolean | 1 | 8 | 62 | 0.075 | 11 | 47795.2 | 47782.3 | 2.7e-4 | 47781.5 |
-| sharp-wedge | sharp, thin-wall, sweepable | 1 | 5 | 8 | 0.083 | 2 | 3200.0 | 3200.0 | 0.0e0 | 3200.0 |
-| sphere-on-plate | near-tangent, boolean | 1 | 12 | 2826 | 0.098 | 687 | 69353.3 | 69509.0 | -2.2e-3 | 69509.7 |
-| spline-vase | organic | 1 | 11 | 3890 | 0.106 | 218 | 94891.2 | 94948.5 | -6.0e-4 | – |
-| stepped-block | blocky, sharp | 1 | 14 | 96 | 0.078 | 32 | 47000.0 | 47000.0 | 0.0e0 | 47000.0 |
-| thin-ribs | thin-wall, blocky | 1 | 26 | 92 | 0.119 | 65 | 32000.0 | 32000.0 | -4.5e-16 | 32000.0 |
-| thin-wall-enclosure | thin-wall, blocky | 1 | 11 | 36 | 0.108 | 8 | 15564.0 | 15564.0 | 0.0e0 | 15564.0 |
-| analytic-sphere | analytic | 1 | 6 | 2784 | 0.069 | 248 | 33353.4 | 33507.9 | -4.6e-3 | 33510.3 |
-| analytic-torus | analytic, genus | 1 | 6 | 7230 | 0.109 | 338 | 37644.1 | 37893.7 | -6.6e-3 | 37899.3 |
-| analytic-cylinder | analytic, sweepable | 1 | 4 | 1004 | 0.049 | 25 | 12533.1 | 12564.7 | -2.5e-3 | 12566.4 |
-| analytic-pipe | analytic, genus, sweepable | 1 | 6 | 1448 | 0.104 | 42 | 94150.9 | 94238.0 | -9.3e-4 | 94247.8 |
-| analytic-plate-hole | analytic, genus, sweepable | 1 | 8 | 112 | 0.224 | 11 | 98449.1 | 98429.5 | 2.0e-4 | 98429.2 |
+| part | categories | bodies | faces | triangles | tol. mm | build ms | enclosed mm³ | SolidCraft mm³ | rel. diff | exact mm³ | boundary | genus | feature edges | check ms |
+|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
+| bent-tube | sweepable, thin-wall, genus | 1 | 18 | 4296 | 0.121 | 432 | 4213.1 | 4159.4 | 1.3e-2 | 4159.8 | ok, 2148 vertices | 1 | 268 | 11.4 |
+| bolted-assembly | assembly, blocky | 3 | 22 | 284 | 0.117 | 57 | 33004.7 | 32999.0 | 1.7e-4 | 32998.9 | ok, 144 vertices | 0, 1, 1 | 104 | 0.7 |
+| cross-cylinders | boolean, sharp | 1 | 10 | 1395 | 0.087 | 7548 | 32197.4 | 32357.8 | -5.0e-3 | 32365.8 | error: Cross: invalid boundary surface: 28 defects; triangles 142 and 184 are oriented inconsistently across edge 0–396 at (-15.488656021315903, 4.577998041166734, 7.853092194999726); triangles 184 and 239 are oriented inconsistently across edge 0–554 at (-18.75000000000002, 5.634596901641276, 7.632255089815496); edge 103–108 at (-0.16007784063084018, 0.03984236936862052, 9.998639270596865) is shared by 3 triangles; … |  |  |  |
+| filleted-bracket | fillet, genus, sweepable | 1 | 13 | 540 | 0.088 | 388 | 32162.3 | 32145.0 | 5.4e-4 | 32144.8 | ok, 268 vertices | 2 | 135 | 1.4 |
+| i-beam | sweepable, thin-wall | 1 | 14 | 44 | 0.208 | 4 | 61600.0 | 61600.0 | 0.0e0 | 61600.0 | ok, 24 vertices | 0 | 36 | 0.1 |
+| loft-transition | organic | 1 | 6 | 292 | 0.075 | 54 | 56390.2 | 56450.6 | -1.1e-3 | – | ok, 148 vertices | 0 | 40 | 0.7 |
+| perforated-plate | genus, thin-wall, sweepable | 1 | 54 | 1644 | 0.144 | 270 | 26822.2 | 26764.6 | 2.2e-3 | 26764.2 | ok, 776 vertices | 24 | 828 | 4.9 |
+| pinhole-block | small-feature, genus | 1 | 38 | 4044 | 0.073 | 206 | 49919.4 | 49909.6 | 2.0e-4 | 49909.5 | ok, 1992 vertices | 16 | 1548 | 12.9 |
+| rounded-box | fillet | 1 | 26 | 2540 | 0.054 | 149 | 23325.7 | 23340.5 | -6.3e-4 | 23340.8 | ok, 1272 vertices | 0 | 216 | 6.7 |
+| shallow-groove | near-tangent, boolean | 1 | 8 | 62 | 0.075 | 11 | 47795.2 | 47782.3 | 2.7e-4 | 47781.5 | ok, 33 vertices | 0 | 24 | 0.1 |
+| sharp-wedge | sharp, thin-wall, sweepable | 1 | 5 | 8 | 0.083 | 2 | 3200.0 | 3200.0 | 0.0e0 | 3200.0 | ok, 6 vertices | 0 | 9 | 0.0 |
+| sphere-on-plate | near-tangent, boolean | 1 | 12 | 2826 | 0.098 | 682 | 69353.3 | 69509.0 | -2.2e-3 | 69509.7 | ok, 1415 vertices | 0 | 124 | 6.1 |
+| spline-vase | organic | 1 | 11 | 3890 | 0.106 | 217 | 94891.2 | 94948.5 | -6.0e-4 | – | ok, 1947 vertices | 0 | 243 | 8.9 |
+| stepped-block | blocky, sharp | 1 | 14 | 96 | 0.078 | 35 | 47000.0 | 47000.0 | 0.0e0 | 47000.0 | ok, 48 vertices | 1 | 60 | 0.4 |
+| thin-ribs | thin-wall, blocky | 1 | 26 | 92 | 0.119 | 69 | 32000.0 | 32000.0 | -4.5e-16 | 32000.0 | ok, 48 vertices | 0 | 72 | 0.5 |
+| thin-wall-enclosure | thin-wall, blocky | 1 | 11 | 36 | 0.108 | 9 | 15564.0 | 15564.0 | 0.0e0 | 15564.0 | ok, 20 vertices | 0 | 28 | 0.1 |
+| analytic-sphere | analytic | 1 | 6 | 2784 | 0.069 | 245 | 33353.4 | 33507.9 | -4.6e-3 | 33510.3 | ok, 1394 vertices | 0 | 96 | 6.1 |
+| analytic-torus | analytic, genus | 1 | 6 | 7230 | 0.109 | 334 | 37644.1 | 37893.7 | -6.6e-3 | 37899.3 | ok, 3615 vertices | 1 | 180 | 16.1 |
+| analytic-cylinder | analytic, sweepable | 1 | 4 | 1004 | 0.049 | 25 | 12533.1 | 12564.7 | -2.5e-3 | 12566.4 | ok, 504 vertices | 0 | 106 | 2.7 |
+| analytic-pipe | analytic, genus, sweepable | 1 | 6 | 1448 | 0.104 | 44 | 94150.9 | 94238.0 | -9.3e-4 | 94247.8 | ok, 724 vertices | 1 | 220 | 3.9 |
+| analytic-plate-hole | analytic, genus, sweepable | 1 | 8 | 112 | 0.224 | 11 | 98449.1 | 98429.5 | 2.0e-4 | 98429.2 | ok, 56 vertices | 1 | 62 | 0.4 |
 
 ## Categories
 
diff --git a/docs/meshing/literature.md b/docs/meshing/literature.md
index 8144ed2..edd7781 100644
--- a/docs/meshing/literature.md
+++ b/docs/meshing/literature.md
@@ -110,6 +110,23 @@ Delaunay methods with Voronoi/ODT improvement, then (2) all-hex and hex-dominant
 
 ---
 
+**A0.5 Boundary surface preparation (implemented in `modelcraft-mesh::surface`).** **Verdict: BUILD ON.**
+- **Guigue and Devillers 2003.** P. Guigue, O. Devillers. *Fast and Robust Triangle–Triangle Overlap Test Using
+  Orientation Predicates.* Journal of Graphics Tools 8(1):25–32, 2003. Triangle overlap decided from `orient3d` /
+  `orient2d` signs alone, so it is exact on the A0.1 predicates. We use the equivalent edge–triangle reduction (two
+  closed triangles meet exactly when an edge of one meets the other) and add the cases of triangles that share a
+  vertex or an edge, which may meet only there.
+- **Ericson 2005.** C. Ericson. *Real-Time Collision Detection.* Morgan Kaufmann, 2005, ch. 6. Bounding volume
+  hierarchies: top-down median-split construction and the self-collision traversal that pairs the triangles to test.
+- **Jacobson, Kavan, Sorkine-Hornung 2013.** *Robust Inside-Outside Segmentation using Generalized Winding Numbers.*
+  ACM TOG 32(4):33, 2013. DOI [10.1145/2461912.2461916](https://doi.org/10.1145/2461912.2461916). The winding number
+  of a vertex of one shell with respect to another tells nesting, so voids must face inwards.
+- **van Oosterom and Strackee 1983.** *The Solid Angle of a Plane Triangle.* IEEE Trans. Biomed. Eng. 30(2):125–126,
+  1983. DOI [10.1109/TBME.1983.325207](https://doi.org/10.1109/TBME.1983.325207). The triangle solid angle the
+  winding number sums.
+
+---
+
 ### A1. Delaunay construction
 
 **A1.1 Bowyer 1981 / Watson 1981, incremental cavity insertion.** **Verdict: build on.**