corpus: the filleted bracket partitioned into 30 blocks, meshed all-hex

filleted-bracket-blocks splits the bracket with SolidCraft's Split Body
into three slabs across its width, the fillet corner halved along its
diagonal and a four-cell O-grid round each hole; the corpus runner meshes
partitioned parts as one conforming block mesh. A heavy test clamps the
upright, loads the foot and checks hex8 and hex20 converge (hex20 tip
deflection within 0.1 % of a fine tet10 solution). Roadmap M5:
conforming interfaces done.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

ATTRIBUTION.md

@@ -48,6 +48,7 @@ is committed and listed here with its exports.
 | assets/icons/mesh.element.add.svg, assets/icons/mesh.element.add.vectorcraft | FieldCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.element.add.json` (`cargo xtask icons`) | MIT OR Apache-2.0 |
 | assets/icons/mesh.hex.decompose.svg, assets/icons/mesh.hex.decompose.vectorcraft | FieldCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.hex.decompose.json` (`cargo xtask icons`) | MIT OR Apache-2.0 |
 | assets/icons/mesh.hex.dominant.svg, assets/icons/mesh.hex.dominant.vectorcraft | FieldCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.hex.dominant.json` (`cargo xtask icons`) | MIT OR Apache-2.0 |
+| assets/icons/mesh.hex.blocks.svg, assets/icons/mesh.hex.blocks.vectorcraft | FieldCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.hex.blocks.json` (`cargo xtask icons`) | MIT OR Apache-2.0 |
 | assets/icons/mesh.hex.mapped.svg, assets/icons/mesh.hex.mapped.vectorcraft | FieldCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.hex.mapped.json` (`cargo xtask icons`) | MIT OR Apache-2.0 |
 | assets/icons/mesh.hex.octree.svg, assets/icons/mesh.hex.octree.vectorcraft | FieldCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.hex.octree.json` (`cargo xtask icons`) | MIT OR Apache-2.0 |
 | assets/icons/mesh.hex.polycube.svg, assets/icons/mesh.hex.polycube.vectorcraft | FieldCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.hex.polycube.json` (`cargo xtask icons`) | MIT OR Apache-2.0 |

ROADMAP.md

@@ -31,7 +31,7 @@ pressure) is a degree of freedom of the existing DOF framework, not a separate s
 | M2 Tet meshing and size field | Size field (global, per face and edge seeds, curvature and proximity; the `fieldcraft_mesh::size::SizeField` trait and edge-division helpers are in); Bowyer–Watson Delaunay (exact, symbolic perturbation, verified; `fieldcraft-mesh::delaunay`, done); robust CDT boundary recovery (Diazzi et al. 2023: implicit LNC Steiner points on segments, gift-wrapping facet recovery, inside/outside labelling, B-rep face tags on boundary tet faces, few-ulp rounding repair checked by the validity gate; `fieldcraft-mesh::cdt` and `mesh.tet.generate`, done: all 20 valid corpus parts); 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; `model.check` naming every missing or bad item; `model.*` and `mesh.tag` commands (model definition: done); writers for Abaqus input (done: `deck.write`), Nastran SOL 101/103 (done: `deck.write`) and LS-DYNA implicit (done); the `.fieldcraft` file (text definition, binary mesh, geometry references reopened on load, journal; atomic saves, autosave and recovery; undo/redo on the journal: done) | in progress |
 | M4 Own solver and first results | Element library per docs/elements.md (1D and 2D first); Linear static and modal solver (tet4/10, wedge6, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions (solver: done — `fieldcraft-fem` kernels, `fieldcraft-solver` with nested dissection, sparse Cholesky and block Lanczos, `solve.static`/`solve.frequency`/`solve.probe`, patch, cantilever, Lamé, Kirsch and beam-frequency tests; reference points with kinematic and distributing couplings and moments by master–slave constraint elimination, verified by rigid motion, clamping, Saint-Venant tension, bending and tube torsion; 1D elements I: rod2/rod3 trusses, point masses, rotary inertia, springs and dashpots in both solvers and all three decks, verified by the method of joints, the three-bar truss, spring–mass, two-mass chain and torsional pendulum frequencies, and damped decay; 2D elements: plane stress, plane strain and axisymmetric tri3, tri6, quad4 (incompatible modes QM6), quad8 with two degrees of freedom per node, centrifugal loads, `mesh.planar` 2D meshes of polygons and planar faces, verified by patch tests, Kirsch, Lamé (plane strain and axisymmetric), rotating discs and cylinders, and a plane stress cantilever; deck mapping to Abaqus CPS/CPE/CAX, Nastran membranes, LS-DYNA ELFORM 12/13); shells and beams (done: quad4 MITC4 and tri3 MITC3 flat shells with a Hughes–Brezzi drilling rotation, two-node Timoshenko beams with profile properties and orientations, six degrees of freedom at shell and beam nodes shared with solids and trusses, section forces and surface stresses, verified by shell patch tests, Navier's plate, cantilever and portal-frame beams, and the Scordelis–Lo and pinched-cylinder benchmarks); results display (done: `results.*` commands, deformed shape to scale, banded contours of displacement, von Mises, principal and component stresses with a legend, mode shapes with their frequencies and animation, probing a node by clicking it, reaction force and moment sums per set, stale results marked after an edit; explicit steps: field frames by slider or animation with contours per frame, energy histories plotted) | started |
-| M5 Structured hex | 2D quad meshing of faces (`fieldcraft-mesh::quad`: transfinite mapped meshes, free all-quad meshes of polygons with holes on prescribed boundary nodes, exact planar quad validity, shared `SizeField` and even edge divisions; done for planar faces); mapped hex of six-sided regions (trivariate TFI, curved faces) and sweep meshing (multi-level extrusions with imprinted cross-sections, revolutions, one-to-one sweeps through linking faces; interval assignment, hex8/hex20 on the geometry, B-rep face tags, `mesh.hex.generate`/`sweep`/`mapped`: done, 16 of 21 corpus parts all-hex); conforming interfaces and tet–hex transition; automatic decomposition into sweepable and mappable regions, proposed as SolidCraft splits | started |
+| M5 Structured hex | 2D quad meshing of faces (`fieldcraft-mesh::quad`: transfinite mapped meshes, free all-quad meshes of polygons with holes on prescribed boundary nodes, exact planar quad validity, shared `SizeField` and even edge divisions; done for planar faces); mapped hex of six-sided regions (trivariate TFI, curved faces) and sweep meshing (multi-level extrusions with imprinted cross-sections, revolutions, one-to-one sweeps through linking faces; interval assignment, hex8/hex20 on the geometry, B-rep face tags, `mesh.hex.generate`/`sweep`/`mapped`: done, 16 of 21 corpus parts all-hex); conforming interfaces between cells (done: `mesh.hex.blocks` meshes a part split with SolidCraft's Split Body into mapped blocks as one conforming mesh — composite faces for seams and tangent faces, interval assignment over all cells, shared faces meshed once, grading off holes and fillets, free-internal-face check; the filleted bracket in 30 cells, hex20 within 0.1 % of a fine tet10 tip deflection); sweepable cells in a conforming mesh and tet–hex transition; automatic decomposition into sweepable and mappable regions, proposed as SolidCraft splits | started |
 | M6 Automatic hex | Octree/grid-dual all-hex for any part; hex-dominant by tet recombination; polycube for blocky parts; smoothing, untangling, boundary padding, layers through thin walls; mesh editing and orphan-mesh import | planned |
 | M7 Assemblies, contact and friction | Needs M15. Part instances from SolidCraft occurrences; tie between non-matching meshes (mortar); implicit contact: surface-to-surface mortar and node-to-surface, penalty and augmented Lagrange, self-contact, automatic pair detection with a bounding-volume search shared with explicit contact; friction: Coulomb with stick–slip, static/kinetic, anisotropic, pressure-, rate- and temperature-dependent coefficients; contact output (pressure, slip, status); thermal contact conductance and frictional heating; wear (Archard); bolt pretension, connectors with behaviour (elastic, damped, stops, locks); amplitudes, multiple steps, restart; nonlinear steps in all three decks. Verified by Hertz (sphere and cylinder), the Cattaneo–Mindlin partial-slip problem, the sliding block, NAFEMS contact benchmarks | planned |
 | M8 Materials, explicit, shells, composites | Plasticity, hyperelasticity, damage data per solver; explicit dynamics decks (Abaqus/Explicit, LS-DYNA); midsurface shells; composite layups and material orientations, continuum shells on swept hex; hybrid (mixed u–p) hex8/hex20/tet10 for nearly incompressible materials (done); heat transfer (solids done: steady and θ-method transient conduction with flux, convection and sources, thermal strain in static steps (solids and 2D) and explicit dynamics (with plasticity), coupled temperature–displacement, `solve.heat`; deck output, 2D/shell/beam conduction planned) | started |

assets/icons/mesh.hex.blocks.svg

@@ -0,0 +1,10 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<svg xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" width="24" height="24" viewBox="0 0 24 24">
+  <title>mesh.hex.blocks</title>
+  <g id="Layer_1" data-name="Layer 1">
+    <path d="M12 3 L20 7 L20 16 L12 20 L4 16 L4 7 Z" fill="#2d4e70" stroke="#dde2ea" stroke-width="1.5" stroke-linecap="round" stroke-linejoin="round"/>
+    <path d="M4 7 L12 11 L20 7 M12 11 L12 20" fill="none" stroke="#dde2ea" stroke-width="1.5" stroke-linecap="round" stroke-linejoin="round"/>
+    <path d="M6 8 L6 17 M10 10 L10 19 M4 9.25 L12 13.25 M4 13.75 L12 17.75 M14 10 L14 19 M18 8 L18 17 M12 13.25 L20 9.25 M12 17.75 L20 13.75 M6 6 L14 10 M10 4 L18 8 M14 4 L6 8 M18 6 L10 10" fill="none" stroke="#4f9ce8" stroke-linecap="round" stroke-linejoin="round"/>
+    <path d="M8 5 L16 9 M16 5 L8 9 M8 9 L8 18 M4 11.5 L12 15.5 M16 9 L16 18 M12 15.5 L20 11.5" fill="none" stroke="#dde2ea" stroke-width="1.5" stroke-linecap="round" stroke-linejoin="round"/>
+  </g>
+</svg>

assets/icons/src/mesh.hex.blocks.json

@@ -0,0 +1 @@
+{"shapes": [{"d": "M12 3 L20 7 L20 16 L12 20 L4 16 L4 7 Z", "fill": "soft", "stroke": "line"}, {"d": "M4 7 L12 11 L20 7 M12 11 L12 20", "stroke": "line"}, {"d": "M6 8 L6 17 M10 10 L10 19 M4 9.25 L12 13.25 M4 13.75 L12 17.75 M14 10 L14 19 M18 8 L18 17 M12 13.25 L20 9.25 M12 17.75 L20 13.75 M6 6 L14 10 M10 4 L18 8 M14 4 L6 8 M18 6 L10 10", "stroke": "accent", "width": 1}, {"d": "M8 5 L16 9 M16 5 L8 9 M8 9 L8 18 M4 11.5 L12 15.5 M16 9 L16 18 M12 15.5 L20 11.5", "stroke": "line", "width": 1.5}]}

crates/engine/Cargo.toml

@@ -27,3 +27,5 @@ log = { workspace = true }
 
 [dev-dependencies]
 zip = { workspace = true }
+# Multi-body designs from SolidCraft scripts in tests.
+fieldcraft-geometry = { workspace = true, features = ["scripts"] }

crates/engine/src/catalogue.rs

@@ -272,6 +272,13 @@ pub static CATALOGUE: &[Capability] = &[
         "Hex mesh swept along extrusions (several levels), around revolutions, or from a source face to a target",
         &["mesh.hex.generate"],
     ),
+    ready(
+        "mesh.hex.blocks",
+        "Block Hex",
+        Mesh,
+        "Structured Hex",
+        "One conforming hex mesh of a part split into six-sided blocks: shared faces meshed once",
+    ),
     soon("mesh.hex.decompose", "Decompose", Mesh, "Structured Hex", "M5", "Find sweepable and mappable regions and propose the splits"),
     soon("mesh.transition", "Transition", Mesh, "Structured Hex", "M5", "Conforming tet–hex interfaces with pyramids"),
     soon("mesh.hex.octree", "All-Hex", Mesh, "Automatic Hex", "M6", "Octree dual all-hex mesh of any part"),

crates/engine/src/cmd/hex.rs

@@ -3,7 +3,7 @@
 use std::borrow::Cow;
 
 use fieldcraft_geometry::{FaceSurface, Shape, SurfaceMesh};
-use fieldcraft_mesh::hex::{HexMesh, HexMethod, HexOptions, Surface, hex_mesh};
+use fieldcraft_mesh::hex::{Cell, HexMesh, HexMethod, HexOptions, MAX_CELLS, Surface, hex_mesh, hex_mesh_cells};
 use fieldcraft_mesh::size::UniformSize;
 use fieldcraft_mesh::surface::BoundarySurface;
 use fieldcraft_model::{GeometryRef, Model};
@@ -129,6 +129,134 @@ pub const HEX_MAPPED: CommandSpec = CommandSpec {
     run: hex_mapped,
 };
 
+pub const HEX_BLOCKS: CommandSpec = CommandSpec {
+    id: "mesh.hex.blocks",
+    label: "Block Hex Mesh",
+    params: r#"{"bodies": ["Cell1", "Cell2"] (default: every body of the open design), "size": 2.0 (element size, mm; default 1/25 of the bounding-box diagonal of all of them), "order": 1 (hex8) or 2 (hex20, mid-edge nodes on the geometry), "tolerance": 0.05 (chordal tolerance of the tessellation, mm; default 1e-3 of that diagonal)} — replaces the mesh with one conforming all-hex mesh of a part partitioned into cells (bodies that meet face to face, such as the pieces of solid.split_body in SolidCraft): every cell a mapped block (six four-sided faces; faces that meet smoothly, such as the two halves of a hole, count as one), division counts assigned over all cells together, each shared face meshed once and its nodes merged, radial lines graded off holes and fillets. Elements are tagged with their cell's body and boundary faces with its B-rep faces (`weld`, `weld_distance` and `feature_angle` as in mesh.tet.generate). Reports the cells, the shared faces, counts and quality; a cell that is not a block, or cells that touch without sharing whole faces, is an error that says which"#,
+    kind: CommandKind::Model,
+    run: hex_blocks,
+};
+
+fn hex_blocks(s: &mut Session, p: &Value) -> Result<Value> {
+    let size = number(p, "size")?;
+    let quadratic = order(p)?;
+    let design = s.design.as_deref().ok_or_else(|| EngineError::Params("open a design first".into()))?;
+    let shapes: Vec<&Shape> = match p.get("bodies") {
+        None => design.bodies.iter().collect(),
+        Some(v) => {
+            let names = v.as_array().ok_or_else(|| EngineError::Params("`bodies` must be a list of body names".into()))?;
+            // Bounded before the duplicate check, which is quadratic in the list.
+            if names.len() > MAX_CELLS {
+                return Err(EngineError::Params(format!("`bodies` lists {} bodies; at most {MAX_CELLS} are meshed together", names.len())));
+            }
+            let mut out: Vec<&Shape> = Vec::with_capacity(names.len());
+            for n in names {
+                let name = n.as_str().ok_or_else(|| EngineError::Params("`bodies` must be a list of body names".into()))?;
+                let shape = design.body(name).ok_or_else(|| EngineError::Params(format!("no body `{name}`")))?;
+                if out.iter().any(|b| b.name == shape.name) {
+                    return Err(EngineError::Params(format!("body `{name}` is listed twice")));
+                }
+                out.push(shape);
+            }
+            out
+        }
+    };
+    if shapes.is_empty() || shapes.len() > MAX_CELLS {
+        return Err(EngineError::Params(format!("{} bodies; mesh.hex.blocks meshes 1 to {MAX_CELLS}", shapes.len())));
+    }
+    // One size and one tessellation tolerance for all cells, from their joint extent, so shared
+    // curves are divided alike.
+    let (mut lo, mut hi) = ([f64::INFINITY; 3], [f64::NEG_INFINITY; 3]);
+    for b in &shapes {
+        let m = b.measure()?;
+        for k in 0..3 {
+            lo[k] = lo[k].min(m.bbox_min[k]);
+            hi[k] = hi[k].max(m.bbox_max[k]);
+        }
+    }
+    let diag = diagonal(lo, hi);
+    if !(diag > 0.0 && diag.is_finite()) {
+        return Err(EngineError::Params("the bodies have no extent".into()));
+    }
+    let tolerance = number(p, "tolerance")?.unwrap_or(DEFAULT_RELATIVE_TOLERANCE * diag);
+    if tolerance < MIN_RELATIVE_TOLERANCE * diag {
+        return Err(EngineError::Params(format!("`tolerance` must be at least {} mm", MIN_RELATIVE_TOLERANCE * diag)));
+    }
+    let h = size.unwrap_or(DEFAULT_RELATIVE_SIZE * diag);
+    if h < MIN_RELATIVE_SIZE * diag {
+        return Err(EngineError::Params(format!("`size` must be at least {} mm", MIN_RELATIVE_SIZE * diag)));
+    }
+    let inputs: Vec<Input> = shapes.iter().map(|b| shape_input(b, Some(tolerance))).collect::<Result<_>>()?;
+    let mut surfaces = Vec::with_capacity(inputs.len());
+    for b in &inputs {
+        let options = super::mesh::surface_options(p, &b.mesh)?;
+        surfaces.push(
+            BoundarySurface::new(&b.mesh.positions, &b.mesh.triangles, &b.mesh.tri_face, &options)
+                .map_err(|e| EngineError::Params(format!("body `{}`: {e}", b.name)))?,
+        );
+    }
+    let cells: Vec<Cell> = surfaces.iter().zip(&inputs).map(|(surface, b)| Cell { name: &b.name, surface, faces: &b.surfaces }).collect();
+    let options = HexOptions {
+        method: HexMethod::Mapped,
+        quadratic,
+        chordal_tolerance: tolerance,
+        max_elements: fieldcraft_model::MAX_ELEMENTS / 8,
+        ..HexOptions::default()
+    };
+    let cm = hex_mesh_cells(&cells, &UniformSize(h), &options).map_err(|e| EngineError::Params(e.to_string()))?;
+    let m = &cm.mesh;
+    // Replace the mesh: definitions stay, nodes, elements and geometry tags are new.
+    let mut model = s.model.without_mesh();
+    for (i, xyz) in m.nodes.iter().enumerate() {
+        model.set_node(label(i)?, *xyz)?;
+    }
+    let mut body_tags = Vec::with_capacity(inputs.len());
+    for b in &inputs {
+        body_tags.push(model.geometry_tag(&GeometryRef::body(b.name.as_str())?)?);
+    }
+    for e in 0..m.element_count() {
+        let cell = cm.element_cell.get(e).map_or(0, |&c| c as usize);
+        let nodes: Vec<u32> = m.element(e).unwrap_or(&[]).iter().map(|&i| i + 1).collect();
+        model.set_element(label(e)?, m.kind, nodes)?;
+        model.tag_element(label(e)?, body_tags.get(cell).copied().ok_or_else(|| EngineError::Params("element without a cell".into()))?)?;
+    }
+    tag_faces(&mut model, &inputs, |e| cm.element_cell.get(e as usize).map_or(0, |&c| c as usize), m, 1)?;
+    s.model = model;
+    let blocks: Vec<Value> =
+        inputs.iter().zip(&cm.blocks).map(|(b, d)| json!({"body": b.name, "divisions": d, "elements": d.iter().product::<usize>()})).collect();
+    Ok(json!({
+        "cells": inputs.len(),
+        "interfaces": cm.interfaces,
+        "description": m.report.description,
+        "kind": m.kind.name(),
+        "elements": m.element_count(),
+        "nodes": m.nodes.len(),
+        "boundary_faces": m.boundary.len(),
+        "quality": {"min_scaled_jacobian": m.report.min_scaled_jacobian, "mean_scaled_jacobian": m.report.mean_scaled_jacobian},
+        "volume": m.report.volume,
+        "surface_volume": m.report.surface_volume,
+        "size": h,
+        "blocks": blocks,
+    }))
+}
+
+/// A 1-based model label for 0-based index `i`.
+fn label(i: usize) -> Result<u32> {
+    u32::try_from(i).ok().and_then(|x| x.checked_add(1)).ok_or_else(|| EngineError::Params("too many nodes or elements".into()))
+}
+
+/// The `order` parameter: false for hex8, true for hex20.
+fn order(p: &Value) -> Result<bool> {
+    match p.get("order") {
+        None => Ok(false),
+        Some(v) => match v.as_u64() {
+            Some(1) => Ok(false),
+            Some(2) => Ok(true),
+            _ => Err(EngineError::Params("`order` must be 1 (hex8) or 2 (hex20)".into())),
+        },
+    }
+}
+
 fn hex_generate(s: &mut Session, p: &Value) -> Result<Value> {
     let all = [HexMethod::Auto, HexMethod::Swept, HexMethod::Extrude, HexMethod::Sweep, HexMethod::Revolve, HexMethod::Mapped];
     generate(s, p, HexMethod::Auto, &all)
@@ -145,14 +273,7 @@ fn hex_mapped(s: &mut Session, p: &Value) -> Result<Value> {
 /// Mesh with `default` or the `method` parameter, which must be one of `allowed`.
 fn generate(s: &mut Session, p: &Value, default: HexMethod, allowed: &[HexMethod]) -> Result<Value> {
     let size = number(p, "size")?;
-    let quadratic = match p.get("order") {
-        None => false,
-        Some(v) => match v.as_u64() {
-            Some(1) => false,
-            Some(2) => true,
-            _ => return Err(EngineError::Params("`order` must be 1 (hex8) or 2 (hex20)".into())),
-        },
-    };
+    let quadratic = order(p)?;
     let method = match p.get("method") {
         None => default,
         Some(v) => v.as_str().and_then(HexMethod::from_name).filter(|m| allowed.contains(m)).ok_or_else(|| {
@@ -197,7 +318,7 @@ fn generate(s: &mut Session, p: &Value, default: HexMethod, allowed: &[HexMethod
             model.tag_element(element_label, body_tag)?;
             element_label = element_label.checked_add(1).ok_or_else(|| EngineError::Params("too many elements".into()))?;
         }
-        tag_faces(&mut model, b, m, first_element)?;
+        tag_faces(&mut model, std::slice::from_ref(b), |_| 0, m, first_element)?;
         report.push(json!({
             "body": b.name,
             "method": m.report.method.name(),
@@ -215,22 +336,20 @@ fn generate(s: &mut Session, p: &Value, default: HexMethod, allowed: &[HexMethod
     Ok(json!({"bodies": report, "nodes": node_label - 1, "elements": element_label - 1}))
 }
 
-/// Tag every boundary element face with its B-rep face.
-fn tag_faces(model: &mut Model, b: &Input<'_>, m: &HexMesh, first_element: u32) -> Result<()> {
-    let mut tags: Vec<Option<u32>> = Vec::new();
+/// Tag every boundary element face with the B-rep face it lies on, of the body
+/// (`bodies[cell(element)]`) its element meshes.
+fn tag_faces(model: &mut Model, bodies: &[Input<'_>], cell: impl Fn(u32) -> usize, m: &HexMesh, first_element: u32) -> Result<()> {
+    let mut tags: std::collections::BTreeMap<(usize, u32), u32> = std::collections::BTreeMap::new();
     for f in &m.boundary {
-        let i = f.brep_face as usize;
-        if tags.len() <= i {
-            tags.resize(i + 1, None);
-        }
-        let tag = match tags.get(i).copied().flatten() {
-            Some(t) => t,
+        let c = cell(f.element);
+        let tag = match tags.get(&(c, f.brep_face)) {
+            Some(&t) => t,
             None => {
+                let b = bodies.get(c).ok_or_else(|| EngineError::Params("a boundary face without its body".into()))?;
+                let i = f.brep_face as usize;
                 let name = b.face_names.get(i).cloned().unwrap_or_else(|| i.to_string());
                 let t = model.geometry_tag(&GeometryRef::face(b.name.as_str(), name)?)?;
-                if let Some(slot) = tags.get_mut(i) {
-                    *slot = Some(t);
-                }
+                tags.insert((c, f.brep_face), t);
                 t
             }
         };
@@ -291,6 +410,89 @@ mod tests {
         assert_eq!(s.model.element_counts(), vec![(fieldcraft_model::ElementKind::Hex8, n)]);
     }
 
+    /// A 20 × 10 × 10 box split by SolidCraft into four cells (x = 10, then y = 5 on each half),
+    /// opened as a design; with `t_junction` the second half is not split, so the first half's
+    /// two cells meet it on one face.
+    fn split_box(t_junction: bool) -> (crate::tests::TempDir, Session) {
+        let mut commands = vec![
+            json!({"command": "solid.box", "params": {"length": 20, "width": 10, "height": 10, "body_name": "Box"}}),
+            json!({"command": "solid.split_body", "params": {"body": "Box", "plane": {"origin": [10, 0, 0], "normal": [1, 0, 0]}}}),
+            json!({"command": "solid.split_body", "params": {"body": "Box (1)", "plane": {"origin": [0, 5, 0], "normal": [0, 1, 0]}}}),
+        ];
+        if !t_junction {
+            commands.push(json!({"command": "solid.split_body", "params": {"body": "Box", "plane": {"origin": [0, 5, 0], "normal": [0, 1, 0]}}}));
+        }
+        let bytes = fieldcraft_geometry::script_design(&json!({"commands": commands}).to_string()).unwrap();
+        let (dir, path) = crate::tests::design_file(&bytes);
+        let mut s = Session::default();
+        s.execute("geometry.open", &json!({"path": path})).unwrap();
+        (dir, s)
+    }
+
+    #[test]
+    fn split_cells_mesh_into_one_conforming_hex_mesh() {
+        let (_dir, mut s) = split_box(false);
+        let r = s.execute("mesh.hex.blocks", &json!({"size": 2.5})).unwrap();
+        assert_eq!(r["cells"], 4, "{r}");
+        assert_eq!(r["interfaces"], 4, "{r}");
+        assert_eq!(r["elements"], 8 * 4 * 4, "{r}");
+        // One 9 × 5 × 5 grid of nodes: the cells' nodes are merged on the interfaces.
+        assert_eq!(r["nodes"], 9 * 5 * 5, "{r}");
+        assert_eq!(s.model.element_counts(), vec![(fieldcraft_model::ElementKind::Hex8, 128)]);
+        assert!((r["volume"].as_f64().unwrap() - 2000.0).abs() < 1e-6, "{r}");
+        // Each cell's elements carry its body, and the outer faces their B-rep faces: a set of
+        // all four bodies skins the box with no internal faces.
+        let bodies = s.execute("geometry.bodies", &json!({})).unwrap();
+        let names: Vec<Value> = bodies.as_array().unwrap().iter().map(|b| json!({"body": b["name"]})).collect();
+        s.execute("model.set.create", &json!({"name": "all", "geometry": names})).unwrap();
+        let all = s.execute("model.set.resolve", &json!({"name": "all"})).unwrap();
+        assert_eq!(all["elements"]["count"], 128, "{all}");
+        assert_eq!(all["faces"]["count"], r["boundary_faces"], "{all}");
+        assert_eq!(r["boundary_faces"], 2 * (8 * 4 + 8 * 4 + 4 * 4), "{r}");
+        let first = bodies[0]["name"].as_str().unwrap();
+        s.execute("model.set.create", &json!({"name": "one", "geometry": [{"body": first}]})).unwrap();
+        assert_eq!(s.execute("model.set.resolve", &json!({"name": "one"})).unwrap()["elements"]["count"], 32);
+        // Hex20 shares mid-edge nodes across the interfaces as well.
+        let q = s.execute("mesh.hex.blocks", &json!({"size": 2.5, "order": 2})).unwrap();
+        assert_eq!(q["kind"], "hex20");
+        assert_eq!(q["nodes"], 9 * 5 * 5 + 8 * 5 * 5 + 9 * 4 * 5 + 9 * 5 * 4, "{q}");
+        // A list of bodies meshes those only.
+        let two = s.execute("mesh.hex.blocks", &json!({"bodies": [bodies[0]["name"], bodies[1]["name"]], "size": 2.5})).unwrap();
+        assert_eq!(two["cells"], 2);
+    }
+
+    #[test]
+    fn cells_that_do_not_conform_and_hostile_params_are_errors() {
+        let (_dir, mut s) = split_box(true);
+        let before = s.model.clone();
+        let e = s.execute("mesh.hex.blocks", &json!({"size": 2.5})).unwrap_err().to_string();
+        assert!(e.contains("do not conform") && e.contains("`Box`"), "{e}");
+        assert_eq!(s.model, before);
+        for p in [
+            json!({"size": 0}),
+            json!({"size": -1}),
+            json!({"size": "big"}),
+            json!({"size": 1e-300}),
+            json!({"order": 3}),
+            json!({"tolerance": -1}),
+            json!({"tolerance": 1e-300}),
+            json!({"bodies": "Box"}),
+            json!({"bodies": []}),
+            json!({"bodies": ["Nope"]}),
+            json!({"bodies": ["Box", "Box"]}),
+            json!({"bodies": [7]}),
+            json!({"size": 1e-3}),
+        ] {
+            assert!(s.execute("mesh.hex.blocks", &p).is_err(), "{p}");
+            assert_eq!(s.model, before, "{p}");
+        }
+        assert!(Session::default().execute("mesh.hex.blocks", &json!({})).is_err());
+        // One cell that is no block says which.
+        let (_dir, mut plate) = opened();
+        let e = plate.execute("mesh.hex.blocks", &json!({})).unwrap_err().to_string();
+        assert!(e.contains("`Body1`") && e.contains("not a logical cube"), "{e}");
+    }
+
     #[test]
     fn hostile_params_are_errors_and_change_nothing() {
         let (_dir, mut s) = opened();

crates/engine/src/cmd/mod.rs

@@ -40,6 +40,7 @@ pub static ALL: &[CommandSpec] = &[
     hex::HEX_GENERATE,
     hex::HEX_SWEEP,
     hex::HEX_MAPPED,
+    hex::HEX_BLOCKS,
     tags::TAG,
     tags::TAG_CLEAR,
     model::UNITS,

crates/engine/src/lib.rs

@@ -530,9 +530,14 @@ mod tests {
 
     /// A design file in a fresh temporary directory (removed when the guard drops).
     pub(crate) fn sample_design_file(plate: &Plate) -> (TempDir, String) {
+        design_file(&plate_design(plate).unwrap())
+    }
+
+    /// Design file bytes written to a temporary directory.
+    pub(crate) fn design_file(bytes: &[u8]) -> (TempDir, String) {
         let dir = TempDir::new();
-        let path = dir.0.join("plate.solidcraft");
-        std::fs::write(&path, plate_design(plate).unwrap()).unwrap();
+        let path = dir.0.join("design.solidcraft");
+        std::fs::write(&path, bytes).unwrap();
         (dir, path.to_string_lossy().into_owned())
     }
 

crates/engine/src/view.rs

@@ -32,6 +32,8 @@ pub struct Palette {
     pub mesh_edge: Rgb,
     pub geometry: Rgb,
     pub geometry_edge: Rgb,
+    /// The bodies of a design with several, in turn.
+    pub bodies: [Rgb; 8],
     /// Faces of nodes without a value in a contour plot.
     pub no_result: Rgb,
     /// The contour colour scale, low to high.
@@ -52,6 +54,16 @@ impl Palette {
         mesh_edge: Rgb(0x1e, 0x23, 0x2a),
         geometry: Rgb(0xb4, 0xba, 0xc4),
         geometry_edge: Rgb(0x2a, 0x30, 0x3a),
+        bodies: [
+            Rgb(0x8f, 0xb4, 0xd9),
+            Rgb(0xe0, 0xa4, 0x6c),
+            Rgb(0x8c, 0xc4, 0x8a),
+            Rgb(0xd6, 0x8c, 0x96),
+            Rgb(0xb0, 0x9c, 0xd8),
+            Rgb(0xd8, 0xc8, 0x70),
+            Rgb(0x7c, 0xc4, 0xc0),
+            Rgb(0xc0, 0xa0, 0x84),
+        ],
         no_result: Rgb(0x9a, 0xa0, 0xaa),
         contour: [Rgb(0x2b, 0x4f, 0xd8), Rgb(0x2f, 0xb8, 0xe0), Rgb(0x3c, 0xc8, 0x64), Rgb(0xf0, 0xd2, 0x3c), Rgb(0xe0, 0x40, 0x2e)],
         glyphs: GlyphColors {
@@ -125,7 +137,7 @@ fn model_scene(d: &DisplayMesh, p: &Palette) -> Scene {
     s.tris = d.tris.iter().map(|t| SceneTri { p: t.p, normal: t.normal, color: p.element(t.kind) }).collect();
     let faces_drawn = d.tris.is_empty();
     if faces_drawn {
-        s.tris.extend(d.surface_tris.iter().map(|t| SceneTri { p: t.p, normal: t.normal, color: p.geometry }));
+        s.tris.extend(d.surface_tris.iter().map(|t| SceneTri { p: t.p, normal: t.normal, color: d.body_color(t, p.geometry, &p.bodies) }));
     }
     s.lines = d.edges.iter().map(|e| SceneLine { a: e.p[0], b: e.p[1], color: p.mesh_edge, width: EDGE_WIDTH }).collect();
     let edge = if faces_drawn { p.mesh_edge } else { p.geometry_edge };

crates/mcp/src/tools.rs

@@ -18,6 +18,7 @@ const MESHERS: &[(&str, &str, &str)] = &[
     ("hex", "mesh.hex.generate", "all-hex mesh of a body"),
     ("hex_sweep", "mesh.hex.sweep", "swept all-hex mesh"),
     ("hex_mapped", "mesh.hex.mapped", "mapped all-hex mesh"),
+    ("hex_blocks", "mesh.hex.blocks", "one conforming all-hex mesh of a design split into six-sided blocks {bodies, size, order}"),
     ("block", "mesh.block", "structured hex block"),
     ("planar", "mesh.planar", "2D elements of a polygon or a planar face"),
 ];

crates/mesh/src/hex/brep.rs

@@ -6,7 +6,16 @@
 //! The hex meshers reason on this topology the way B-rep sweepers do (Blacker, *The Cooper
 //! tool*, Proc. 5th International Meshing Roundtable, 1996): which faces are sources, targets and
 //! linking faces, and which edges bound them.
+//!
+//! **Virtual topology.** [`Brep::composite`] joins faces that meet smoothly along all of their
+//! common boundary (the two halves of a cylinder the kernel splits at a seam, a plane running
+//! tangentially into a fillet, coplanar pieces) into one composite face, so their seam is no edge
+//! of the mesher's topology: the composite surfaces of virtual topology (A. Sheffer, M. Bercovier,
+//! T. Blacker, J. Clements, *Virtual topology operators for meshing*, IJCGA 10(3):309–331, 2000).
+//! A composite face keeps an analytic surface when every member lies on the first member's;
+//! otherwise nodes go onto its triangles.
 
+use std::borrow::Cow;
 use std::collections::BTreeMap;
 
 use super::geom::{Surface, add, dist, scale, unit};
@@ -59,6 +68,9 @@ impl Edge {
 
 pub(crate) struct Brep<'a> {
     pub surface: &'a BoundarySurface,
+    /// The face tag of each triangle: the surface's own tags, or for composite faces the smallest
+    /// tag of the composite.
+    tags: Cow<'a, [u32]>,
     pub faces: Vec<Face>,
     pub edges: Vec<Edge>,
     /// Whether each surface vertex is a B-rep vertex.
@@ -72,9 +84,29 @@ pub(crate) struct Brep<'a> {
 /// Relative geometric tolerance (of the bounding-box diagonal).
 pub(crate) const RELATIVE_EPS: f64 = 1e-7;
 
+/// Largest angle (radians) between the surface normals of two faces along their common boundary
+/// for them to join into a composite face: 1°.
+const SMOOTH_ANGLE: f64 = 0.017_453_292_519_943_295;
+/// Distance from a member's vertices to the first member's analytic surface (relative to the
+/// bounding-box diagonal) within which a composite face keeps that surface.
+const RELATIVE_ON_SURFACE: f64 = 1e-6;
+
 impl<'a> Brep<'a> {
     /// Recover the topology. `surfaces[tag]` is the analytic surface of B-rep face `tag`, if any.
     pub(crate) fn new(surface: &'a BoundarySurface, surfaces: &[Option<Surface>]) -> Result<Self> {
+        Self::build(surface, surfaces, Cow::Borrowed(surface.tri_face()), false)
+    }
+
+    /// [`Brep::new`] with smoothly meeting faces joined into composite faces (virtual topology).
+    pub(crate) fn composite(surface: &'a BoundarySurface, surfaces: &[Option<Surface>]) -> Result<Self> {
+        let tags = composite_tags(surface, surfaces);
+        Self::build(surface, surfaces, Cow::Owned(tags), true)
+    }
+
+    fn build(surface: &'a BoundarySurface, surfaces: &[Option<Surface>], tags: Cow<'a, [u32]>, composite: bool) -> Result<Self> {
+        if tags.len() != surface.triangles().len() {
+            return Err(Error::Failed("a face tag per triangle is needed".into()));
+        }
         let pts = surface.points();
         let mut lo = [f64::INFINITY; 3];
         let mut hi = [f64::NEG_INFINITY; 3];
@@ -91,7 +123,7 @@ impl<'a> Brep<'a> {
         let eps = RELATIVE_EPS * diagonal;
         // Faces in tag order.
         let mut by_tag: BTreeMap<u32, Vec<u32>> = BTreeMap::new();
-        for (t, &f) in surface.tri_face().iter().enumerate() {
+        for (t, &f) in tags.iter().enumerate() {
             by_tag.entry(f).or_default().push(t as u32);
         }
         let index: BTreeMap<u32, usize> = by_tag.keys().enumerate().map(|(i, &f)| (f, i)).collect();
@@ -109,12 +141,13 @@ impl<'a> Brep<'a> {
             }
             verts.sort_unstable();
             verts.dedup();
-            let analytic = surfaces.get(id as usize).copied().flatten().and_then(|s| s.checked());
+            let analytic =
+                if composite { composite_surface(surface, surfaces, tris, RELATIVE_ON_SURFACE * diagonal) } else { own_surface(surfaces, id) };
             let normal = unit(nsum).unwrap_or([0.0; 3]);
             let plane = planarity(pts, &verts, normal, analytic, eps);
             faces.push(Face { id, triangles: tris.clone(), vertices: verts, area, surface: analytic, plane, loops: Vec::new() });
         }
-        let face_of = |t: u32| surface.tri_face().get(t as usize).and_then(|f| index.get(f)).copied();
+        let face_of = |t: u32| tags.get(t as usize).and_then(|f| index.get(f)).copied();
         // Oriented boundary segments: `from → to` with face `left` traversing it.
         struct Seg {
             from: u32,
@@ -128,6 +161,10 @@ impl<'a> Brep<'a> {
             let (Some(left), Some(right)) = (face_of(fe.triangles[0]), face_of(fe.triangles[1])) else {
                 return Err(Error::Failed("feature edge without faces".into()));
             };
+            // The seam inside a composite face is no edge.
+            if left == right {
+                continue;
+            }
             // Orient each segment as the lower-indexed face traverses it, so that the segments of
             // one edge chain head to tail.
             let [lo, hi] = fe.vertices;
@@ -236,7 +273,7 @@ impl<'a> Brep<'a> {
                 face.loops = loops;
             }
         }
-        Ok(Brep { surface, faces, edges, is_vertex, edge_triangles, diagonal })
+        Ok(Brep { surface, tags, faces, edges, is_vertex, edge_triangles, diagonal })
     }
 
     pub(crate) fn point(&self, v: u32) -> P {
@@ -276,9 +313,14 @@ impl<'a> Brep<'a> {
 
     /// The faces that use the triangle's tag (index into `faces`).
     pub(crate) fn face_of_triangle(&self, t: u32) -> Option<usize> {
-        let tag = *self.surface.tri_face().get(t as usize)?;
+        let tag = *self.tags.get(t as usize)?;
         self.faces.binary_search_by_key(&tag, |f| f.id).ok()
     }
+
+    /// The face tag of triangle `t` in this topology (a composite face's tag for its members).
+    pub(crate) fn tag(&self, t: u32) -> Option<u32> {
+        self.tags.get(t as usize).copied()
+    }
 }
 
 /// A point and outward normal of a planar face: the analytic plane, or the plane fitted through
@@ -300,3 +342,73 @@ fn planarity(pts: &[P], verts: &[u32], normal: P, analytic: Option<Surface>, eps
     // tolerance, not to ours.
     (flat || matches!(analytic, Some(Surface::Plane { .. }))).then_some((point, n))
 }
+
+/// The analytic surface of face `tag` from the geometry, checked.
+fn own_surface(surfaces: &[Option<Surface>], tag: u32) -> Option<Surface> {
+    surfaces.get(tag as usize).copied().flatten().and_then(|s| s.checked())
+}
+
+/// The analytic surface of a composite face (its triangles `tris`): the first member's, when the
+/// vertices of every member lie on it within `tol`.
+fn composite_surface(surface: &BoundarySurface, surfaces: &[Option<Surface>], tris: &[u32], tol: f64) -> Option<Surface> {
+    let tri_face = surface.tri_face();
+    let first = tris.iter().filter_map(|&t| tri_face.get(t as usize)).min()?;
+    let s = own_surface(surfaces, *first)?;
+    let pts = surface.points();
+    let on = |v: u32| pts.get(v as usize).is_some_and(|&p| s.project(p).is_some_and(|q| dist(p, q) <= tol));
+    tris.iter().all(|&t| surface.triangles().get(t as usize).is_some_and(|tri| tri.iter().all(|&v| on(v)))).then_some(s)
+}
+
+/// Face tags with smoothly meeting faces joined: two faces join when, at every boundary segment
+/// between them, their analytic surfaces' normals are within [`SMOOTH_ANGLE`] of parallel and the
+/// two triangles there face the same side. Each composite face takes its smallest member tag.
+fn composite_tags(surface: &BoundarySurface, surfaces: &[Option<Surface>]) -> Vec<u32> {
+    let tri_face = surface.tri_face();
+    let pts = surface.points();
+    let point = |v: u32| pts.get(v as usize).copied().unwrap_or([0.0; 3]);
+    let tri_normal = |t: u32| {
+        surface.triangles().get(t as usize).and_then(|tri| {
+            let [a, b, c] = tri.map(point);
+            unit(cross(sub(b, a), sub(c, a)))
+        })
+    };
+    let cos_min = SMOOTH_ANGLE.cos();
+    // Per pair of faces (low tag, high tag): whether every segment between them is smooth.
+    let mut pairs: BTreeMap<(u32, u32), bool> = BTreeMap::new();
+    for fe in surface.feature_edges().iter().filter(|e| e.kind == FeatureKind::FaceBoundary) {
+        let (Some(&a), Some(&b)) = (tri_face.get(fe.triangles[0] as usize), tri_face.get(fe.triangles[1] as usize)) else { continue };
+        if a == b {
+            continue;
+        }
+        let m = scale(add(point(fe.vertices[0]), point(fe.vertices[1])), 0.5);
+        let normals = (own_surface(surfaces, a).and_then(|s| s.normal(m)), own_surface(surfaces, b).and_then(|s| s.normal(m)));
+        let same_side = match (tri_normal(fe.triangles[0]), tri_normal(fe.triangles[1])) {
+            (Some(x), Some(y)) => dot(x, y) > 0.0,
+            _ => false,
+        };
+        let smooth = same_side && matches!(normals, (Some(x), Some(y)) if dot(x, y).abs() >= cos_min);
+        let entry = pairs.entry((a.min(b), a.max(b))).or_insert(true);
+        *entry = *entry && smooth;
+    }
+    // Union–find over tags, the root always the smaller tag.
+    let mut parent: BTreeMap<u32, u32> = BTreeMap::new();
+    fn root(parent: &BTreeMap<u32, u32>, mut x: u32) -> u32 {
+        // Each step goes to a smaller tag, so this ends.
+        while let Some(&p) = parent.get(&x) {
+            if p == x {
+                break;
+            }
+            x = p;
+        }
+        x
+    }
+    for (&(a, b), &smooth) in &pairs {
+        if smooth {
+            let (ra, rb) = (root(&parent, a), root(&parent, b));
+            if ra != rb {
+                parent.insert(ra.max(rb), ra.min(rb));
+            }
+        }
+    }
+    tri_face.iter().map(|&f| root(&parent, f)).collect()
+}

crates/mesh/src/hex/cells.rs

@@ -0,0 +1,796 @@
+//! Conforming structured hex meshes of a part partitioned into cells (multi-block meshing).
+//!
+//! An analyst meshes a part no single sweep fits by splitting it into cells that each are a
+//! logical cube (six four-sided faces, possibly composite), meshing every cell as a mapped block
+//! and joining the blocks so that they conform: the multi-block grids of J. F. Thompson, B. K.
+//! Soni, N. P. Weatherill (eds.), *Handbook of Grid Generation*, CRC Press 1999, ch. 3 and 13.
+//! [`hex_mesh_cells`] does it in five steps:
+//!
+//! 1. **Topology.** Each cell's B-rep is read with smoothly meeting faces joined
+//!    ([`super::brep::Brep::composite`]: the halves of a hole, a plane running into a fillet) and
+//!    recognized as a logical cube ([`super::mapped::recognize`]). Corners of different cells at
+//!    the same point (within the merge tolerance, the one geometric decision) are one vertex;
+//!    cell edges between the same two vertices along the same curve are one edge; cell faces
+//!    bounded by the same four edges are one face, shared by at most two cells (an interface).
+//!    Everything after this works on these identities, never on coordinates.
+//! 2. **Interval assignment** over the edges of all cells together ([`super::interval`]): the
+//!    four edges along each direction of every cell are equal, so a count is the same on both
+//!    sides of every interface. Edges that run away from a curved face (radial lines off a hole,
+//!    the diagonal off a fillet) are then graded: the spacing at that end matches the spacing
+//!    along the curved face there and grows by at most [`GROWTH`] per element towards the target
+//!    size; their counts are reassigned with the grading's goals.
+//! 3. **Edges** are divided once, on the geometry, and their nodes belong to the edge.
+//! 4. **Faces** get one transfinite grid each (Gordon and Hall), from their four edges and
+//!    projected onto the face; an interface's grid is made once and both cells use its nodes.
+//! 5. **Cells** fill their interior by trivariate transfinite interpolation (Gordon and Thiel)
+//!    of their six face grids, blended by the edges' arc-length distribution.
+//!
+//! Nodes are merged by identity (a vertex, an edge's k-th node, a face grid's node), so the
+//! cells conform exactly: no node is duplicated on an interface and none hangs. The result goes
+//! through the same finishing as a single body ([`super::finish`]): the certified validity gate,
+//! the tiling and fidelity checks, a closed 2-manifold boundary, and no boundary face lying on
+//! another cell (a free internal face, where cells touch without sharing a face).
+
+use std::collections::BTreeMap;
+
+use super::brep::Brep;
+use super::geom::{TriangleTree, add, dist, scale};
+use super::mapped::{self, Cube, face_ring, other};
+use super::{Cells, Ctx, HexMesh, HexMethod, HexOptions, RELATIVE_TOLERANCE, Raw, Surface, finish, project_to_face, project_to_faces};
+use crate::quad::{lattice_index, transfinite_grid, transfinite_volume};
+use crate::size::{SizeField, divide_polyline, ideal_divisions};
+use crate::surface::BoundarySurface;
+use crate::vec3::{P, cross, dot, norm, sub};
+use crate::{Error, Result};
+
+/// One cell: its name (for messages), its boundary surface (face-tagged) and the analytic
+/// surfaces of its faces by tag.
+#[derive(Clone, Copy)]
+pub struct Cell<'a> {
+    pub name: &'a str,
+    pub surface: &'a BoundarySurface,
+    pub faces: &'a [Option<Surface>],
+}
+
+/// A conforming mesh of several cells.
+#[derive(Debug, Clone, PartialEq)]
+pub struct CellsMesh {
+    /// The mesh; boundary faces are tagged with the B-rep faces of their element's cell.
+    pub mesh: HexMesh,
+    /// The cell of each element (index into the cells given).
+    pub element_cell: Vec<u32>,
+    /// Each cell's block: divisions along its three logical directions.
+    pub blocks: Vec<[usize; 3]>,
+    /// Faces shared by two cells.
+    pub interfaces: usize,
+}
+
+/// The most cells meshed together.
+pub const MAX_CELLS: usize = 4096;
+/// Spacing growth from one element to the next along a graded edge (1.2: each at most 20 %
+/// longer than the one before).
+pub const GROWTH: f64 = 1.2;
+/// Segments a graded edge's polyline is cut into for the integration of its grading (the size
+/// helpers take the size constant along a sample).
+const GRADED_SEGMENTS: usize = 64;
+
+/// An edge shared by cells: its end vertices (`a < b`), the cell that divides it and how.
+struct Edge {
+    a: usize,
+    b: usize,
+    /// The owning cell, its edge and whether that edge runs from `a` to `b`.
+    owner: (usize, usize, bool),
+    length: f64,
+}
+
+/// A face shared by one or two cells.
+struct Face {
+    /// The owner's ring of vertices (the order of [`face_ring`]).
+    ring: [usize; 4],
+    /// The cells, with the direction and side of the face in each; the first owns the grid.
+    uses: Vec<(usize, usize, usize)>,
+}
+
+/// Mesh `cells` (bodies that together make the part, meeting face to face) into one conforming
+/// all-hex mesh, every cell a mapped block. `opts.method` must be [`HexMethod::Auto`] or
+/// [`HexMethod::Mapped`]; a cell that is not a logical cube, or cells that meet without sharing
+/// whole faces, are [`Error::Unsupported`] with the reason.
+pub fn hex_mesh_cells(cells: &[Cell], field: &dyn SizeField, opts: &HexOptions) -> Result<CellsMesh> {
+    if !matches!(opts.method, HexMethod::Auto | HexMethod::Mapped) {
+        return Err(Error::Invalid(format!("conforming cells are meshed as mapped blocks; method `{}` does not apply", opts.method.name())));
+    }
+    if !(opts.chordal_tolerance >= 0.0 && opts.chordal_tolerance.is_finite()) {
+        return Err(Error::Invalid(format!("the chordal tolerance must be a non-negative number, got {}", opts.chordal_tolerance)));
+    }
+    if cells.is_empty() || cells.len() > MAX_CELLS {
+        return Err(Error::Invalid(format!("{} cells; 1 to {MAX_CELLS} are meshed together", cells.len())));
+    }
+    if let Some(c) = cells.iter().find(|c| c.surface.shells().len() != 1) {
+        return Err(Error::Unsupported(format!("cell `{}` has {} shells; a cell is one solid without voids", c.name, c.surface.shells().len())));
+    }
+    let surfaces: Vec<&BoundarySurface> = cells.iter().map(|c| c.surface).collect();
+    super::check_estimate(&surfaces, field, opts)?;
+    let breps: Vec<Brep> = cells.iter().map(|c| Brep::composite(c.surface, c.faces)).collect::<Result<_>>()?;
+    let trees: Vec<TriangleTree> = cells.iter().map(|c| TriangleTree::new(c.surface)).collect::<Result<_>>()?;
+    let ctxs: Vec<Ctx> =
+        breps.iter().zip(&trees).map(|(brep, tree)| Ctx { brep, tree, field, opts, tol: RELATIVE_TOLERANCE * brep.diagonal }).collect();
+    let names: Vec<&str> = cells.iter().map(|c| c.name).collect();
+    let cubes: Vec<Cube> = breps
+        .iter()
+        .zip(&names)
+        .map(|(b, name)| mapped::recognize(b).map_err(|e| Error::Unsupported(format!("cell `{name}`: {}", super::plain(&e)))))
+        .collect::<Result<_>>()?;
+    let topo = Topology::new(&breps, &cubes, &names, opts.chordal_tolerance)?;
+    let divisions = topo.intervals(&ctxs, &cubes, field)?;
+    let mut builder =
+        Builder { topo: &topo, ctxs: &ctxs, cubes: &cubes, divisions: &divisions, field, opts, nodes: Vec::new(), edge_nodes: Vec::new() };
+    let (raw, element_cell, blocks) = builder.build()?;
+    let shared = topo.interfaces();
+    let mesh = finish(&Cells { ctxs: &ctxs, of: &element_cell, names: &names }, raw, HexMethod::Mapped)?;
+    Ok(CellsMesh { mesh, element_cell, blocks, interfaces: shared })
+}
+
+/// Vertices, edges and faces shared by the cells.
+struct Topology {
+    /// Shared vertices: where they are.
+    vertices: Vec<P>,
+    edges: Vec<Edge>,
+    faces: Vec<Face>,
+    /// Per cell: the vertex at each corner, `corner[i][j][k]`.
+    corners: Vec<[[[usize; 2]; 2]; 2]>,
+    /// Per cell, per direction: the four edges (in the order of [`Cube::groups`]).
+    cell_edges: Vec<[[usize; 4]; 3]>,
+    /// Per cell: the face on each side of each direction.
+    cell_faces: Vec<[[usize; 2]; 3]>,
+}
+
+/// Points at most this far apart (relative to the diagonal of all cells) are one vertex.
+const RELATIVE_MERGE: f64 = 1e-6;
+
+impl Topology {
+    fn new(breps: &[Brep], cubes: &[Cube], names: &[&str], chordal: f64) -> Result<Topology> {
+        let name = |c: usize| names.get(c).copied().unwrap_or("?");
+        let (mut lo, mut hi) = ([f64::INFINITY; 3], [f64::NEG_INFINITY; 3]);
+        for b in breps {
+            for p in b.surface.points() {
+                for k in 0..3 {
+                    lo[k] = lo[k].min(p[k]);
+                    hi[k] = hi[k].max(p[k]);
+                }
+            }
+        }
+        let merge = RELATIVE_MERGE * dist(lo, hi);
+        if !(merge > 0.0 && merge.is_finite()) {
+            return Err(Error::Invalid("the cells have no extent".into()));
+        }
+        // Vertices: corners within `merge` of each other, found through a grid of that spacing.
+        let mut vertices: Vec<P> = Vec::new();
+        let mut grid: BTreeMap<[i64; 3], Vec<usize>> = BTreeMap::new();
+        let cellkey = |p: P| p.map(|x| (x / merge).floor() as i64);
+        let mut corners = Vec::with_capacity(cubes.len());
+        for (c, (brep, cube)) in breps.iter().zip(cubes).enumerate() {
+            let mut ids = [[[0usize; 2]; 2]; 2];
+            let mut seen: Vec<usize> = Vec::with_capacity(8);
+            for (i, plane) in ids.iter_mut().enumerate() {
+                for (j, row) in plane.iter_mut().enumerate() {
+                    for (k, slot) in row.iter_mut().enumerate() {
+                        let p = brep.point(cube.at(i, j, k));
+                        let key = cellkey(p);
+                        let mut found: Option<usize> = None;
+                        for dx in -1..=1 {
+                            for dy in -1..=1 {
+                                for dz in -1..=1 {
+                                    let near = [key[0] + dx, key[1] + dy, key[2] + dz];
+                                    for &v in grid.get(&near).map(Vec::as_slice).unwrap_or(&[]) {
+                                        if vertices.get(v).is_some_and(|&x| dist(x, p) <= merge) {
+                                            if found.is_some_and(|f| f != v) {
+                                                return Err(Error::Unsupported(format!(
+                                                    "cell `{}`: a corner at {p:?} is within the merge distance of two vertices",
+                                                    name(c)
+                                                )));
+                                            }
+                                            found = Some(v);
+                                        }
+                                    }
+                                }
+                            }
+                        }
+                        let v = match found {
+                            Some(v) => v,
+                            None => {
+                                vertices.push(p);
+                                grid.entry(key).or_default().push(vertices.len() - 1);
+                                vertices.len() - 1
+                            }
+                        };
+                        if seen.contains(&v) {
+                            return Err(Error::Unsupported(format!("cell `{}`: two corners at the same point {p:?}", name(c))));
+                        }
+                        seen.push(v);
+                        *slot = v;
+                    }
+                }
+            }
+            corners.push(ids);
+        }
+        // Edges: the same two vertices and the same curve (midpoints within the merge distance
+        // and the tessellation's chordal deviation).
+        let edge_tol = merge + 4.0 * chordal;
+        let mut edges: Vec<Edge> = Vec::new();
+        let mut by_ends: BTreeMap<(usize, usize), Vec<(usize, P)>> = BTreeMap::new();
+        let mut cell_edges = Vec::with_capacity(cubes.len());
+        for (c, (brep, cube)) in breps.iter().zip(cubes).enumerate() {
+            let ids = corners.get(c).copied().ok_or_else(|| Error::Failed("cell corners".into()))?;
+            let mut out = [[0usize; 4]; 3];
+            for (d, row) in out.iter_mut().enumerate() {
+                for (m, slot) in row.iter_mut().enumerate() {
+                    let (u, v) = edge_ends(&ids, d, m);
+                    let (e, fwd) = cube.groups.get(d).and_then(|g| g.get(m)).copied().ok_or_else(|| Error::Failed("cube groups".into()))?;
+                    let pts = brep.edge_points(e, fwd);
+                    let (mid, length) = midpoint(&pts);
+                    let (a, b) = (u.min(v), u.max(v));
+                    let list = by_ends.entry((a, b)).or_default();
+                    let id = match list.iter().find(|(_, q)| dist(*q, mid) <= edge_tol) {
+                        Some(&(id, _)) => id,
+                        None => {
+                            edges.push(Edge { a, b, owner: (c, e, fwd == (u == a)), length });
+                            list.push((edges.len() - 1, mid));
+                            edges.len() - 1
+                        }
+                    };
+                    *slot = id;
+                }
+            }
+            cell_edges.push(out);
+        }
+        // Faces: the same four edges, on the same surface, with the cells on opposite sides.
+        let mut faces: Vec<Face> = Vec::new();
+        let mut by_edges: BTreeMap<[usize; 4], usize> = BTreeMap::new();
+        let mut cell_faces = Vec::with_capacity(cubes.len());
+        for (c, (brep, cube)) in breps.iter().zip(cubes).enumerate() {
+            let ids = corners.get(c).copied().ok_or_else(|| Error::Failed("cell corners".into()))?;
+            let mut out = [[0usize; 2]; 3];
+            for (d, row) in out.iter_mut().enumerate() {
+                for (side, slot) in row.iter_mut().enumerate() {
+                    let ring = face_ring(d, side).map(|[i, j, k]| ids[i][j][k]);
+                    let mut key = [0usize; 4];
+                    for (m, k) in key.iter_mut().enumerate() {
+                        let (u, v) = (ring[m], ring[(m + 1) % 4]);
+                        let groups = cell_edges.get(c).ok_or_else(|| Error::Failed("cell edges".into()))?;
+                        *k = edge_between(&edges, groups, u, v).ok_or_else(|| Error::Failed("a face side is no cell edge".into()))?;
+                    }
+                    key.sort_unstable();
+                    let face = cube.faces.get(d).and_then(|r| r.get(side)).copied().ok_or_else(|| Error::Failed("cube faces".into()))?;
+                    let id = match by_edges.get(&key) {
+                        Some(&id) => {
+                            let f = faces.get_mut(id).ok_or_else(|| Error::Failed("face index".into()))?;
+                            if f.uses.len() >= 2 {
+                                return Err(Error::Unsupported(format!("cell `{}`: a face is shared by more than two cells", name(c))));
+                            }
+                            if let Some(&(c0, d0, s0)) = f.uses.first() {
+                                if c0 == c {
+                                    return Err(Error::Unsupported(format!("cell `{}`: two faces of the cell have the same four edges", name(c))));
+                                }
+                                let f0 = cubes.get(c0).and_then(|q| q.faces.get(d0)).and_then(|r| r.get(s0)).copied().unwrap_or(0);
+                                let b0 = breps.get(c0).ok_or_else(|| Error::Failed("cell index".into()))?;
+                                let (centre0, area0) = face_centre(b0, f0);
+                                let (centre1, area1) = face_centre(brep, face);
+                                if dist(centre0, centre1) > edge_tol + 1e-3 * norm(area0).sqrt() || dot(area0, area1) >= 0.0 {
+                                    return Err(Error::Unsupported(format!(
+                                        "cells `{}` and `{}` share the four edges of a face but not the face (different surfaces, or on the same side)",
+                                        name(c0),
+                                        name(c)
+                                    )));
+                                }
+                            }
+                            f.uses.push((c, d, side));
+                            id
+                        }
+                        None => {
+                            faces.push(Face { ring, uses: vec![(c, d, side)] });
+                            by_edges.insert(key, faces.len() - 1);
+                            faces.len() - 1
+                        }
+                    };
+                    *slot = id;
+                }
+            }
+            cell_faces.push(out);
+        }
+        Ok(Topology { vertices, edges, faces, corners, cell_edges, cell_faces })
+    }
+
+    /// Faces shared by two cells.
+    fn interfaces(&self) -> usize {
+        self.faces.iter().filter(|f| f.uses.len() == 2).count()
+    }
+
+    /// The edge between vertices `u` and `v` among cell `c`'s edges, and whether it runs from
+    /// `u` to `v`.
+    fn edge(&self, c: usize, u: usize, v: usize) -> Option<(usize, bool)> {
+        let e = edge_between(&self.edges, self.cell_edges.get(c)?, u, v)?;
+        Some((e, self.edges.get(e)?.a == u))
+    }
+
+    /// How every edge is divided, all cells together; graded edges get the counts their
+    /// grading asks for.
+    fn intervals(&self, ctxs: &[Ctx], cubes: &[Cube], field: &dyn SizeField) -> Result<Vec<Division>> {
+        let polylines: Vec<Vec<P>> = self.edges.iter().map(|e| self.polyline(ctxs, e)).collect::<Result<_>>()?;
+        let goals: Vec<f64> = polylines.iter().map(|p| ideal_divisions(field, p)).collect::<Result<_>>()?;
+        let first = self.solve(&goals)?;
+        // Ends of edges that leave a curved face: the spacing along that face there.
+        let ends = self.graded_ends(ctxs, cubes, &first);
+        let mut graded = goals;
+        let mut divisions = Vec::with_capacity(polylines.len());
+        for ((poly, ends), goal) in polylines.into_iter().zip(ends).zip(graded.iter_mut()) {
+            let poly = if ends.iter().any(Option::is_some) {
+                let dense = densify(&poly);
+                *goal = goal.max(ideal_divisions(&Grading { field, ends }, &dense)?);
+                dense
+            } else {
+                poly
+            };
+            divisions.push(Division { count: 1, ends, polyline: poly });
+        }
+        for (d, n) in divisions.iter_mut().zip(self.solve(&graded)?) {
+            d.count = n;
+        }
+        Ok(divisions)
+    }
+
+    fn solve(&self, goals: &[f64]) -> Result<Vec<usize>> {
+        let mut ia = super::interval::Intervals::default();
+        for &g in goals {
+            ia.var(g, 1);
+        }
+        for groups in &self.cell_edges {
+            for g in groups {
+                for &e in g.iter().skip(1) {
+                    ia.equal.push((vec![g[0]], vec![e]));
+                }
+            }
+        }
+        ia.solve()
+    }
+
+    /// For every edge, at each end (`a`, `b`): the point and the spacing to start from when the
+    /// edge leaves a curved face there (ends on it without lying on it) — the smallest spacing,
+    /// by `counts`, of that face's edges at the end — or `None`.
+    fn graded_ends(&self, ctxs: &[Ctx], cubes: &[Cube], counts: &[usize]) -> Vec<GradedEnds> {
+        // The edges of every curved face, by the vertices they meet.
+        let mut curved: Vec<[usize; 4]> = Vec::new();
+        let mut curved_at: BTreeMap<usize, Vec<usize>> = BTreeMap::new();
+        for f in &self.faces {
+            let Some(&(c, d, side)) = f.uses.first() else { continue };
+            let planar = cubes
+                .get(c)
+                .and_then(|q| q.faces.get(d))
+                .and_then(|r| r.get(side))
+                .and_then(|&bf| ctxs.get(c).and_then(|x| x.brep.faces.get(bf)))
+                .is_none_or(|bf| bf.plane.is_some());
+            if planar {
+                continue;
+            }
+            let mut es = [usize::MAX; 4];
+            for (m, slot) in es.iter_mut().enumerate() {
+                if let Some((e, _)) = self.edge(c, f.ring[m], f.ring[(m + 1) % 4]) {
+                    *slot = e;
+                }
+            }
+            for &v in &f.ring {
+                curved_at.entry(v).or_default().push(curved.len());
+            }
+            curved.push(es);
+        }
+        let spacing = |e: usize| -> f64 { self.edges.get(e).map_or(f64::INFINITY, |x| x.length / counts.get(e).copied().unwrap_or(1).max(1) as f64) };
+        self.edges
+            .iter()
+            .enumerate()
+            .map(|(e, edge)| {
+                [edge.a, edge.b].map(|v| {
+                    let mut best: Option<f64> = None;
+                    for face in curved_at.get(&v).map(Vec::as_slice).unwrap_or(&[]).iter().filter_map(|&k| curved.get(k)) {
+                        if face.contains(&e) {
+                            continue;
+                        }
+                        for &o in face {
+                            if self.edges.get(o).is_some_and(|x| x.a == v || x.b == v) {
+                                let s = spacing(o);
+                                best = Some(best.map_or(s, |b| b.min(s)));
+                            }
+                        }
+                    }
+                    best.and_then(|s| self.vertices.get(v).map(|&x| (x, s)))
+                })
+            })
+            .collect()
+    }
+
+    /// The owner's polyline of an edge, from `a` to `b`.
+    fn polyline(&self, ctxs: &[Ctx], e: &Edge) -> Result<Vec<P>> {
+        let (c, be, fwd) = e.owner;
+        let ctx = ctxs.get(c).ok_or_else(|| Error::Failed("edge owner".into()))?;
+        let mut pts = ctx.brep.edge_points(be, fwd);
+        // The ends are the shared vertices.
+        if let (Some(first), Some(a)) = (pts.first_mut(), self.vertices.get(e.a)) {
+            *first = *a;
+        }
+        if let (Some(last), Some(b)) = (pts.last_mut(), self.vertices.get(e.b)) {
+            *last = *b;
+        }
+        if pts.len() < 2 {
+            return Err(Error::Failed("an edge without points".into()));
+        }
+        Ok(pts)
+    }
+}
+
+/// The edge among a cell's edges (`groups`) between vertices `u` and `v`.
+fn edge_between(edges: &[Edge], groups: &[[usize; 4]; 3], u: usize, v: usize) -> Option<usize> {
+    let key = (u.min(v), u.max(v));
+    groups.iter().flatten().copied().find(|&e| edges.get(e).is_some_and(|x| (x.a, x.b) == key))
+}
+
+/// The corners (low, high) of edge `m` of direction `d` of a cell whose corners are `ids`, in
+/// the order of [`Cube::groups`].
+fn edge_ends(ids: &[[[usize; 2]; 2]; 2], d: usize, m: usize) -> (usize, usize) {
+    let (s, t) = ((m / 2) & 1, m & 1);
+    match d {
+        0 => (ids[0][s][t], ids[1][s][t]),
+        1 => (ids[s][0][t], ids[s][1][t]),
+        _ => (ids[s][t][0], ids[s][t][1]),
+    }
+}
+
+/// The point halfway along a polyline (by arc length) and its length.
+fn midpoint(pts: &[P]) -> (P, f64) {
+    let length: f64 = pts.windows(2).map(|w| dist(w[0], w[1])).sum();
+    let mut walked = 0.0;
+    for w in pts.windows(2) {
+        let l = dist(w[0], w[1]);
+        if walked + l >= 0.5 * length && l > 0.0 {
+            let t = (0.5 * length - walked) / l;
+            return (add(w[0], scale(sub(w[1], w[0]), t)), length);
+        }
+        walked += l;
+    }
+    (pts.first().copied().unwrap_or([0.0; 3]), length)
+}
+
+/// A B-rep face's area-weighted centre and its area vector (outward).
+fn face_centre(brep: &Brep, f: usize) -> (P, P) {
+    let (mut c, mut n, mut total) = ([0.0; 3], [0.0; 3], 0.0);
+    for &t in brep.faces.get(f).map(|x| x.triangles.as_slice()).unwrap_or(&[]) {
+        let Some(tri) = brep.surface.triangles().get(t as usize) else { continue };
+        let [a, b, d] = tri.map(|v| brep.point(v));
+        let area = cross(sub(b, a), sub(d, a));
+        let w = 0.5 * norm(area);
+        c = add(c, scale(add(add(a, b), d), w / 3.0));
+        n = add(n, scale(area, 0.5));
+        total += w;
+    }
+    (if total > 0.0 { scale(c, 1.0 / total) } else { c }, n)
+}
+
+/// A polyline cut into at least [`GRADED_SEGMENTS`] pieces (for the integration of a grading).
+fn densify(pts: &[P]) -> Vec<P> {
+    let length: f64 = pts.windows(2).map(|w| dist(w[0], w[1])).sum();
+    let step = length / GRADED_SEGMENTS as f64;
+    let mut out = Vec::with_capacity(pts.len() + GRADED_SEGMENTS);
+    for w in pts.windows(2) {
+        let l = dist(w[0], w[1]);
+        let k = if step > 0.0 { ((l / step).ceil() as usize).clamp(1, GRADED_SEGMENTS) } else { 1 };
+        for s in 0..k {
+            out.push(add(w[0], scale(sub(w[1], w[0]), s as f64 / k as f64)));
+        }
+    }
+    if let Some(&last) = pts.last() {
+        out.push(last);
+    }
+    out
+}
+
+/// Each end of an edge (`a`, `b`): where grading starts from and the spacing there.
+type GradedEnds = [Option<(P, f64)>; 2];
+
+/// How an edge is divided: its count, its graded ends and the polyline (from `a` to `b`; cut
+/// finer when graded) it is divided along.
+struct Division {
+    count: usize,
+    ends: GradedEnds,
+    polyline: Vec<P>,
+}
+
+/// The size along a graded edge: the field's, capped near each graded end by the spacing there
+/// growing by [`GROWTH`] per element (linearly with the distance from that end).
+struct Grading<'a> {
+    field: &'a dyn SizeField,
+    ends: GradedEnds,
+}
+
+impl SizeField for Grading<'_> {
+    fn size(&self, p: [f64; 3]) -> f64 {
+        let mut h = self.field.size(p);
+        for (q, s) in self.ends.iter().flatten() {
+            h = h.min(s + (GROWTH - 1.0) * dist(p, *q));
+        }
+        h
+    }
+}
+
+/// The nodes of a cell's lattice, or of a face's grid, by identity.
+struct Builder<'t, 'c, 'a> {
+    topo: &'t Topology,
+    ctxs: &'c [Ctx<'a>],
+    cubes: &'c [Cube],
+    divisions: &'c [Division],
+    field: &'c dyn SizeField,
+    opts: &'c HexOptions,
+    nodes: Vec<P>,
+    /// Per edge: its node ids from `a` to `b`.
+    edge_nodes: Vec<Vec<u32>>,
+}
+
+impl Builder<'_, '_, '_> {
+    fn count(&self, e: usize) -> usize {
+        self.divisions.get(e).map_or(1, |d| d.count)
+    }
+
+    fn push(&mut self, p: P) -> Result<u32> {
+        let id = u32::try_from(self.nodes.len()).map_err(|_| Error::Limit("too many nodes".into()))?;
+        self.nodes.push(p);
+        Ok(id)
+    }
+
+    fn build(&mut self) -> Result<(Raw, Vec<u32>, Vec<[usize; 3]>)> {
+        let (topo, field, max_elements) = (self.topo, self.field, self.opts.max_elements);
+        // Block sizes and the element budget before any node is made.
+        let mut blocks = Vec::with_capacity(topo.cell_edges.len());
+        let mut total = 0usize;
+        for groups in &topo.cell_edges {
+            let n = groups.map(|g| self.count(g[0]));
+            let nodes = (n[0] + 1).checked_mul(n[1] + 1).and_then(|x| x.checked_mul(n[2] + 1)).filter(|&x| x <= mapped::MAX_NODES);
+            if nodes.is_none() {
+                return Err(Error::Limit(format!("a {} × {} × {} block exceeds {} nodes", n[0], n[1], n[2], mapped::MAX_NODES)));
+            }
+            total = total.saturating_add(n[0] * n[1] * n[2]);
+            blocks.push(n);
+        }
+        if total > max_elements {
+            return Err(Error::Limit(format!("the mesh would have {total} hexahedra; the limit is {max_elements}")));
+        }
+        // Vertices.
+        for &v in &topo.vertices {
+            self.push(v)?;
+        }
+        // Edges: divided once, on the geometry.
+        for (edge, division) in topo.edges.iter().zip(self.divisions) {
+            let n = division.count;
+            let pts = if division.ends.iter().any(Option::is_some) {
+                divide_polyline(&Grading { field, ends: division.ends }, &division.polyline, n)?
+            } else {
+                divide_polyline(field, &division.polyline, n)?
+            };
+            let (c, be, _) = edge.owner;
+            let ctx = self.ctxs.get(c).ok_or_else(|| Error::Failed("edge owner".into()))?;
+            let faces = ctx.brep.edges.get(be).map(|x| x.faces.to_vec()).unwrap_or_default();
+            let mut ids = Vec::with_capacity(n + 1);
+            ids.push(edge.a as u32);
+            for q in pts.iter().take(n).skip(1) {
+                let p = project_to_faces(ctx, &faces, *q);
+                ids.push(self.push(p)?);
+            }
+            ids.push(edge.b as u32);
+            self.edge_nodes.push(ids);
+        }
+        // Faces: one grid each, from the owner's side.
+        let mut face_grids: Vec<FaceGrid> = Vec::with_capacity(topo.faces.len());
+        for face in &topo.faces {
+            face_grids.push(self.face_grid(face)?);
+        }
+        // Cells.
+        let mut hexes: Vec<[u32; 8]> = Vec::with_capacity(total);
+        let mut element_cell: Vec<u32> = Vec::with_capacity(total);
+        for (c, &n) in blocks.iter().enumerate() {
+            let ids = self.cell_lattice(c, n, &face_grids)?;
+            let [ni, nj, nk] = n;
+            let id = |i: usize, j: usize, k: usize| ids.get(lattice_index(n, [i, j, k])).copied().unwrap_or(0);
+            for k in 0..nk {
+                for j in 0..nj {
+                    for i in 0..ni {
+                        hexes.push([
+                            id(i, j, k),
+                            id(i + 1, j, k),
+                            id(i + 1, j + 1, k),
+                            id(i, j + 1, k),
+                            id(i, j, k + 1),
+                            id(i + 1, j, k + 1),
+                            id(i + 1, j + 1, k + 1),
+                            id(i, j + 1, k + 1),
+                        ]);
+                        element_cell.push(c as u32);
+                    }
+                }
+            }
+        }
+        let shared = topo.interfaces();
+        let description = format!(
+            "{} cell{} as conforming mapped blocks, {} shared face{}",
+            blocks.len(),
+            if blocks.len() == 1 { "" } else { "s" },
+            shared,
+            if shared == 1 { "" } else { "s" }
+        );
+        let nodes = std::mem::take(&mut self.nodes);
+        Ok((Raw { nodes, hexes, description }, element_cell, blocks))
+    }
+
+    /// The nodes of edge `e` from vertex `u`.
+    fn edge_from(&self, e: usize, u: usize) -> Result<Vec<u32>> {
+        let edge = self.topo.edges.get(e).ok_or_else(|| Error::Failed("edge index".into()))?;
+        let mut ids = self.edge_nodes.get(e).cloned().ok_or_else(|| Error::Failed("edge nodes".into()))?;
+        if edge.b == u {
+            ids.reverse();
+        } else if edge.a != u {
+            return Err(Error::Failed("an edge does not start at the vertex asked".into()));
+        }
+        Ok(ids)
+    }
+
+    /// A face's transfinite grid in its owner's frame.
+    fn face_grid(&mut self, face: &Face) -> Result<FaceGrid> {
+        let &(c, d, side) = face.uses.first().ok_or_else(|| Error::Failed("a face without a cell".into()))?;
+        let r = face.ring;
+        let side_nodes = |b: &Self, u: usize, v: usize| -> Result<Vec<u32>> {
+            let (e, _) = b.topo.edge(c, u, v).ok_or_else(|| Error::Failed("a face side is no edge".into()))?;
+            b.edge_from(e, u)
+        };
+        let bottom = side_nodes(self, r[0], r[1])?;
+        let right = side_nodes(self, r[1], r[2])?;
+        let top = side_nodes(self, r[3], r[2])?;
+        let left = side_nodes(self, r[0], r[3])?;
+        let (na, nb) = (bottom.len() - 1, left.len() - 1);
+        if top.len() != na + 1 || right.len() != nb + 1 {
+            return Err(Error::Failed("opposite sides of a face have different counts".into()));
+        }
+        let at = |b: &Self, ids: &[u32]| -> Vec<P> { ids.iter().map(|&i| b.nodes.get(i as usize).copied().unwrap_or([0.0; 3])).collect() };
+        let grid = transfinite_grid(&at(self, &bottom), &at(self, &right), &at(self, &top), &at(self, &left))?;
+        let ctx = self.ctxs.get(c).ok_or_else(|| Error::Failed("face owner".into()))?;
+        let brep_face = self.cubes.get(c).and_then(|q| q.faces.get(d)).and_then(|row| row.get(side)).copied();
+        let mut ids = vec![0u32; (na + 1) * (nb + 1)];
+        for b in 0..=nb {
+            for a in 0..=na {
+                let id = if b == 0 {
+                    bottom.get(a).copied()
+                } else if b == nb {
+                    top.get(a).copied()
+                } else if a == 0 {
+                    left.get(b).copied()
+                } else if a == na {
+                    right.get(b).copied()
+                } else {
+                    let q = grid.get(b * (na + 1) + a).copied().ok_or_else(|| Error::Failed("face grid".into()))?;
+                    let q = brep_face.and_then(|f| project_to_face(ctx, f, q)).unwrap_or(q);
+                    Some(self.push(q)?)
+                };
+                if let Some(slot) = ids.get_mut(b * (na + 1) + a) {
+                    *slot = id.ok_or_else(|| Error::Failed("face side node".into()))?;
+                }
+            }
+        }
+        Ok(FaceGrid { ring: r, na, nb, ids })
+    }
+
+    /// Cell `c`'s lattice of node ids: boundary from the shared face grids, interior new.
+    fn cell_lattice(&mut self, c: usize, n: [usize; 3], face_grids: &[FaceGrid]) -> Result<Vec<u32>> {
+        let topo = self.topo;
+        let corners = topo.corners.get(c).copied().ok_or_else(|| Error::Failed("cell corners".into()))?;
+        let faces = topo.cell_faces.get(c).copied().ok_or_else(|| Error::Failed("cell faces".into()))?;
+        let total = (n[0] + 1) * (n[1] + 1) * (n[2] + 1);
+        let mut ids: Vec<Option<u32>> = vec![None; total];
+        for d in 0..3 {
+            let (da, db) = other(d);
+            for side in 0..2 {
+                let gf = faces.get(d).and_then(|r| r.get(side)).copied().ok_or_else(|| Error::Failed("cell face".into()))?;
+                let grid = face_grids.get(gf).ok_or_else(|| Error::Failed("face grid".into()))?;
+                let ring = face_ring(d, side).map(|[i, j, k]| corners[i][j][k]);
+                let map = grid.map_from(ring, n[da], n[db])?;
+                for b in 0..=n[db] {
+                    for a in 0..=n[da] {
+                        let mut l = [0usize; 3];
+                        l[d] = side * n[d];
+                        l[da] = a;
+                        l[db] = b;
+                        if let Some(slot) = ids.get_mut(lattice_index(n, l)) {
+                            *slot = Some(map(a, b)?);
+                        }
+                    }
+                }
+            }
+        }
+        // Interior: trivariate interpolation, blended by the edges' arc-length distribution.
+        let mut x: Vec<Option<P>> = ids.iter().map(|i| i.and_then(|i| self.nodes.get(i as usize).copied())).collect();
+        let groups = topo.cell_edges.get(c).copied().ok_or_else(|| Error::Failed("cell edges".into()))?;
+        let mut blend: [Vec<f64>; 3] = [Vec::new(), Vec::new(), Vec::new()];
+        for (d, g) in groups.iter().enumerate() {
+            let mut sum = vec![0.0; n[d] + 1];
+            for (m, &e) in g.iter().enumerate() {
+                let pts: Vec<P> = self
+                    .edge_from(e, edge_ends(&corners, d, m).0)?
+                    .iter()
+                    .map(|&i| self.nodes.get(i as usize).copied().unwrap_or([0.0; 3]))
+                    .collect();
+                let params = crate::quad::tfi::params(&pts, "edge")?;
+                if params.len() != sum.len() {
+                    return Err(Error::Failed("an edge's count differs from its block's".into()));
+                }
+                for (acc, p) in sum.iter_mut().zip(params) {
+                    *acc += 0.25 * p;
+                }
+            }
+            if let Some(first) = sum.first_mut() {
+                *first = 0.0;
+            }
+            if let Some(last) = sum.last_mut() {
+                *last = 1.0;
+            }
+            blend[d] = sum;
+        }
+        transfinite_volume(&mut x, n, [&blend[0], &blend[1], &blend[2]])?;
+        let mut out = Vec::with_capacity(total);
+        for (slot, p) in ids.iter().zip(&x) {
+            out.push(match (slot, p) {
+                (Some(id), _) => *id,
+                (None, Some(p)) => self.push(*p)?,
+                (None, None) => return Err(Error::Failed("a block node is missing".into())),
+            });
+        }
+        Ok(out)
+    }
+}
+
+/// A face's grid of node ids in its owner's frame: `(na + 1) × (nb + 1)`, node `(a, b)` at
+/// `b·(na + 1) + a`, corners at the owner's ring.
+struct FaceGrid {
+    ring: [usize; 4],
+    na: usize,
+    nb: usize,
+    ids: Vec<u32>,
+}
+
+impl FaceGrid {
+    /// The lookup of node `(a, b)` of a cell's view of this face, whose ring (the order of
+    /// [`face_ring`]) is `ring` and counts `na × nb`.
+    fn map_from(&self, ring: [usize; 4], na: usize, nb: usize) -> Result<impl Fn(usize, usize) -> Result<u32> + '_> {
+        // The owner's lattice position of each of the cell's corners.
+        let corner = |v: usize| -> Option<[i64; 2]> {
+            let k = self.ring.iter().position(|&x| x == v)?;
+            let (a, b) = [(0, 0), (self.na, 0), (self.na, self.nb), (0, self.nb)].get(k).copied()?;
+            Some([a as i64, b as i64])
+        };
+        let (Some(o), Some(pa), Some(pb)) = (corner(ring[0]), corner(ring[1]), corner(ring[3])) else {
+            return Err(Error::Failed("a cell's face does not have the shared face's corners".into()));
+        };
+        let step = |p: [i64; 2], n: usize| -> Option<[i64; 2]> {
+            let d = [p[0] - o[0], p[1] - o[1]];
+            let n = n as i64;
+            // One axis of the owner's lattice, `n` long.
+            (n > 0 && d.iter().filter(|&&x| x != 0).count() == 1 && d.iter().all(|&x| x == 0 || x.abs() == n)).then(|| d.map(|x| x / n))
+        };
+        let (Some(du), Some(dv)) = (step(pa, na), step(pb, nb)) else {
+            return Err(Error::Failed("a cell's view of a shared face does not match its grid".into()));
+        };
+        let (na_owner, ids) = (self.na, &self.ids);
+        Ok(move |a: usize, b: usize| -> Result<u32> {
+            let (a, b) = (a as i64, b as i64);
+            let u = o[0] + a * du[0] + b * dv[0];
+            let v = o[1] + a * du[1] + b * dv[1];
+            usize::try_from(v)
+                .ok()
+                .zip(usize::try_from(u).ok())
+                .and_then(|(v, u)| ids.get(v * (na_owner + 1) + u))
+                .copied()
+                .ok_or_else(|| Error::Failed("shared face grid index".into()))
+        })
+    }
+}

crates/mesh/src/hex/geom.rs

@@ -115,6 +115,29 @@ impl Surface {
         };
         q.iter().all(|x| x.is_finite()).then_some(q)
     }
+
+    /// The unit normal of the surface at its point closest to `p`, up to sign, or `None` where
+    /// the closest point is not unique.
+    pub(crate) fn normal(&self, p: P) -> Option<P> {
+        match *self {
+            Surface::Plane { normal, .. } => unit(normal),
+            Surface::Cylinder { origin, axis, .. } => {
+                let w = sub(p, origin);
+                unit(sub(w, scale(axis, dot(w, axis))))
+            }
+            Surface::Cone { apex, axis, half_angle } => {
+                let w = sub(p, apex);
+                let radial = unit(sub(w, scale(axis, dot(w, axis))))?;
+                unit(sub(scale(radial, half_angle.cos()), scale(axis, half_angle.sin())))
+            }
+            Surface::Sphere { center, .. } => unit(sub(p, center)),
+            Surface::Torus { center, axis, major, .. } => {
+                let w = sub(p, center);
+                let ring = add(center, scale(unit(sub(w, scale(axis, dot(w, axis))))?, major));
+                unit(sub(p, ring))
+            }
+        }
+    }
 }
 
 /// The closest point to `p` on triangle `abc`.
@@ -274,8 +297,18 @@ impl<'a> TriangleTree<'a> {
 
     /// The closest triangle to `p` among those `accept` admits (by triangle index).
     pub(crate) fn closest(&self, p: P, accept: &dyn Fn(u32) -> bool) -> Option<Hit> {
+        self.search(p, f64::INFINITY, accept)
+    }
+
+    /// [`TriangleTree::closest`] among the triangles within distance `r` of `p` (the search
+    /// skips everything farther).
+    pub(crate) fn closest_within(&self, p: P, r: f64, accept: &dyn Fn(u32) -> bool) -> Option<Hit> {
+        self.search(p, r * r, accept).filter(|h| h.dist <= r)
+    }
+
+    fn search(&self, p: P, limit2: f64, accept: &dyn Fn(u32) -> bool) -> Option<Hit> {
         let mut best: Option<Hit> = None;
-        let mut best2 = f64::INFINITY;
+        let mut best2 = limit2;
         let mut stack = vec![0usize];
         let mut visits = 0usize;
         while let Some(id) = stack.pop() {

crates/mesh/src/hex/mapped.rs

@@ -12,6 +12,7 @@
 
 use std::collections::BTreeMap;
 
+use super::brep::Brep;
 use super::interval::Intervals;
 use super::{Ctx, Raw, project_to_face, project_to_faces};
 use crate::quad::{lattice_index, transfinite_faces, transfinite_volume};
@@ -20,10 +21,30 @@ use crate::vec3::{P, det, sub};
 use crate::{Error, Result};
 
 /// The most nodes of a mapped block.
-const MAX_NODES: usize = 1 << 26;
+pub(crate) const MAX_NODES: usize = 1 << 26;
 
-pub(crate) fn mesh(ctx: &Ctx) -> Result<Raw> {
-    let brep = ctx.brep;
+/// A logical cube found in a B-rep: its corners, its edges by logical direction and its faces.
+pub(crate) struct Cube {
+    /// The surface vertex at each corner, `corner[i][j][k]`, right-handed.
+    pub corner: [[[u32; 2]; 2]; 2],
+    /// The four edges along each logical direction `d`, oriented from the low corner to the high
+    /// one, at the other two coordinates `(s, t) ∈ {0, 1}²` in the order (0, 0), (0, 1), (1, 0),
+    /// (1, 1).
+    pub groups: [[(usize, bool); 4]; 3],
+    /// The B-rep face (index into `faces`) on each side of each direction: `faces[d][side]`.
+    pub faces: [[usize; 2]; 3],
+}
+
+impl Cube {
+    /// The corner at logical coordinates `(i, j, k)` (each 0 or 1).
+    pub(crate) fn at(&self, i: usize, j: usize, k: usize) -> u32 {
+        self.corner[i & 1][j & 1][k & 1]
+    }
+}
+
+/// Recognize a logical cube: six faces with one loop of four edges each, three edges at each of
+/// eight corners. Anything else is [`Error::Unsupported`] with the reason.
+pub(crate) fn recognize(brep: &Brep) -> Result<Cube> {
     if brep.faces.len() != 6 || brep.edges.len() != 12 {
         return Err(Error::Unsupported(format!("not a logical cube: {} faces and {} edges (6 and 12 needed)", brep.faces.len(), brep.edges.len())));
     }
@@ -61,7 +82,6 @@ pub(crate) fn mesh(ctx: &Ctx) -> Result<Raw> {
     };
     let c111 = common(c110, c101, a).filter(|&x| neighbours(c011).contains(&x));
     let Some(c111) = c111 else { return Err(Error::Unsupported("not a logical cube: the corner graph is not a cube".into())) };
-    // corner[i][j][k]
     let corner = [[[c000, c], [b, c011]], [[a, c101], [c110, c111]]];
     let at = |i: usize, j: usize, k: usize| -> u32 { corner[i][j][k] };
     let mut distinct: Vec<u32> = corner.iter().flatten().flatten().copied().collect();
@@ -74,20 +94,65 @@ pub(crate) fn mesh(ctx: &Ctx) -> Result<Raw> {
     let edge_between = |u: u32, v: u32| -> Option<(usize, bool)> {
         adj.get(&u)?.iter().find(|x| x.0 == v).map(|x| (x.1, brep.edges.get(x.1).is_some_and(|e| e.first() == u)))
     };
-    // The four edges along each direction, from the low corner to the high one.
-    let mut groups: [Vec<(usize, bool)>; 3] = [Vec::new(), Vec::new(), Vec::new()];
+    let mut groups = [[(0usize, true); 4]; 3];
     for (d, group) in groups.iter_mut().enumerate() {
-        for s in 0..2 {
-            for t in 0..2 {
-                let (lo, hi) = match d {
-                    0 => (at(0, s, t), at(1, s, t)),
-                    1 => (at(s, 0, t), at(s, 1, t)),
-                    _ => (at(s, t, 0), at(s, t, 1)),
-                };
-                group.push(edge_between(lo, hi).ok_or_else(|| Error::Unsupported("not a logical cube: a corner pair has no edge".into()))?);
-            }
+        for (m, slot) in group.iter_mut().enumerate() {
+            let (s, t) = (m / 2, m % 2);
+            let (lo, hi) = match d {
+                0 => (at(0, s, t), at(1, s, t)),
+                1 => (at(s, 0, t), at(s, 1, t)),
+                _ => (at(s, t, 0), at(s, t, 1)),
+            };
+            *slot = edge_between(lo, hi).ok_or_else(|| Error::Unsupported("not a logical cube: a corner pair has no edge".into()))?;
         }
     }
+    // The face bounded by the four edges of each side.
+    let mut faces = [[0usize; 2]; 3];
+    for (d, row) in faces.iter_mut().enumerate() {
+        for (side, slot) in row.iter_mut().enumerate() {
+            let ring = face_ring(d, side).map(|[i, j, k]| at(i, j, k));
+            let mut bounding: Vec<usize> = (0..4).filter_map(|m| edge_between(ring[m], ring[(m + 1) % 4]).map(|x| x.0)).collect();
+            bounding.sort_unstable();
+            *slot = brep
+                .faces
+                .iter()
+                .position(|f| {
+                    let mut es: Vec<usize> = f.loops.iter().flatten().map(|x| x.0).collect();
+                    es.sort_unstable();
+                    es == bounding
+                })
+                .ok_or_else(|| Error::Unsupported("not a logical cube: a side's four edges bound no face".into()))?;
+        }
+    }
+    Ok(Cube { corner, groups, faces })
+}
+
+/// The corners of side `side` of direction `d` in the order (0, 0), (1, 0), (1, 1), (0, 1) of its
+/// two other directions (lower first).
+pub(crate) fn face_ring(d: usize, side: usize) -> [[usize; 3]; 4] {
+    let (da, db) = other(d);
+    [(0, 0), (1, 0), (1, 1), (0, 1)].map(|(a, b)| {
+        let mut c = [0usize; 3];
+        c[d] = side;
+        c[da] = a;
+        c[db] = b;
+        c
+    })
+}
+
+/// The two directions other than `d`, lower first.
+pub(crate) fn other(d: usize) -> (usize, usize) {
+    match d {
+        0 => (1, 2),
+        1 => (0, 2),
+        _ => (0, 1),
+    }
+}
+
+pub(crate) fn mesh(ctx: &Ctx) -> Result<Raw> {
+    let brep = ctx.brep;
+    let cube = recognize(brep)?;
+    let groups = cube.groups;
     let mut ia = Intervals::default();
     for e in 0..brep.edges.len() {
         ia.var(ideal_divisions(ctx.field, &brep.edge_points(e, true))?, 1);
@@ -99,7 +164,7 @@ pub(crate) fn mesh(ctx: &Ctx) -> Result<Raw> {
         }
     }
     let counts = ia.solve()?;
-    let count = |g: &Vec<(usize, bool)>| g.first().and_then(|x| counts.get(x.0)).copied().unwrap_or(1);
+    let count = |g: &[(usize, bool); 4]| g.first().and_then(|x| counts.get(x.0)).copied().unwrap_or(1);
     let (ni, nj, nk) = (count(&groups[0]), count(&groups[1]), count(&groups[2]));
     let total = (ni + 1).checked_mul(nj + 1).and_then(|x| x.checked_mul(nk + 1)).filter(|&x| x <= MAX_NODES);
     let total = total.ok_or_else(|| Error::Limit(format!("a {ni} × {nj} × {nk} block exceeds {MAX_NODES} nodes")))?;
@@ -135,26 +200,8 @@ pub(crate) fn mesh(ctx: &Ctx) -> Result<Raw> {
     // Faces: grids from their edges, projected onto the B-rep face their four edges bound.
     let n = [ni, nj, nk];
     transfinite_faces(&mut x, n, None, |d, side, grid| {
-        let (da, db) = match d {
-            0 => (1, 2),
-            1 => (0, 2),
-            _ => (0, 1),
-        };
-        let logical = |a: usize, b: usize| {
-            let mut c = [0usize; 3];
-            c[d] = side;
-            c[da] = a;
-            c[db] = b;
-            at(c[0], c[1], c[2])
-        };
-        let ring = [logical(0, 0), logical(1, 0), logical(1, 1), logical(0, 1)];
-        let mut bounding: Vec<usize> = (0..4).filter_map(|m| edge_between(ring[m], ring[(m + 1) % 4]).map(|x| x.0)).collect();
-        bounding.sort_unstable();
-        let face = brep.faces.iter().position(|f| {
-            let mut es: Vec<usize> = f.loops.iter().flatten().map(|x| x.0).collect();
-            es.sort_unstable();
-            es == bounding
-        });
+        let (da, db) = other(d);
+        let face = cube.faces.get(d).and_then(|r| r.get(side)).copied();
         let (na, nb) = (n[da], n[db]);
         for b in 1..nb {
             for a in 1..na {

crates/mesh/src/hex/mod.rs

@@ -19,6 +19,11 @@
 //! A body none of these fits is [`Error::Unsupported`] with each method's reason: the input for
 //! decomposition into sweepable and mappable regions.
 //!
+//! **Partitioned parts.** [`hex_mesh_cells`] ([`cells`]) meshes a part split into cells that
+//! meet face to face (SolidCraft split bodies) as one conforming mesh: every cell a mapped block,
+//! interval assignment over all cells together, each shared face meshed once and its nodes
+//! merged by identity.
+//!
 //! **Common steps.** Interval assignment ([`interval`]) picks the division counts of shared curves
 //! together. Boundary nodes on analytic faces are placed on the exact surface; quadratic hexahedra
 //! (hex20) get their mid-edge nodes projected onto the faces (and onto both faces of a B-rep
@@ -31,6 +36,7 @@
 
 mod arrangement;
 mod brep;
+mod cells;
 mod geom;
 mod interval;
 mod layered;
@@ -43,6 +49,7 @@ mod tests;
 
 use std::collections::HashMap;
 
+pub use cells::{Cell, CellsMesh, GROWTH, MAX_CELLS, hex_mesh_cells};
 pub use geom::Surface;
 
 use brep::Brep;
@@ -198,16 +205,7 @@ pub fn hex_mesh(surface: &BoundarySurface, faces: &[Option<Surface>], field: &dy
             surface.shells().len()
         )));
     }
-    // A mesh far beyond the element limit is refused before any work: the solid's volume over
-    // the cube of the smallest size the field asks for at the boundary.
-    let mut h_min = f64::INFINITY;
-    for p in surface.points() {
-        h_min = h_min.min(crate::size::checked_size(field, *p)?);
-    }
-    let estimate = surface.volume().abs() / h_min.powi(3);
-    if !(estimate <= 8.0 * opts.max_elements as f64) {
-        return Err(Error::Limit(format!("about {estimate:.3e} hexahedra at this size; the limit is {}", opts.max_elements)));
-    }
+    check_estimate(&[surface], field, opts)?;
     let brep = Brep::new(surface, faces)?;
     let tree = TriangleTree::new(surface)?;
     let ctx = Ctx { brep: &brep, tree: &tree, field, opts, tol: RELATIVE_TOLERANCE * brep.diagonal };
@@ -253,7 +251,7 @@ pub fn hex_mesh(surface: &BoundarySurface, faces: &[Option<Surface>], field: &dy
                 Err(e @ Error::Limit(_)) => return Err(e),
                 Err(e) => last = Some(e),
             },
-            _ => match mapped::mesh(&ctx).and_then(|raw| finish(&ctx, raw, method)) {
+            _ => match mapped::mesh(&ctx).and_then(|raw| finish(&Cells::one(&ctx), raw, method)) {
                 Ok(mesh) => return Ok(mesh),
                 Err(e @ Error::Limit(_)) => return Err(e),
                 Err(e) => last = Some(e),
@@ -266,6 +264,23 @@ pub fn hex_mesh(surface: &BoundarySurface, faces: &[Option<Surface>], field: &dy
     Err(Error::Unsupported(format!("no structured hex method fits this body ({})", reasons.join("; "))))
 }
 
+/// Refuse a mesh far beyond the element limit before any work: the solids' volume over the cube
+/// of the smallest size the field asks for on their boundaries.
+fn check_estimate(surfaces: &[&BoundarySurface], field: &dyn SizeField, opts: &HexOptions) -> Result<()> {
+    let (mut h_min, mut volume) = (f64::INFINITY, 0.0);
+    for s in surfaces {
+        for p in s.points() {
+            h_min = h_min.min(crate::size::checked_size(field, *p)?);
+        }
+        volume += s.volume().abs();
+    }
+    let estimate = volume / h_min.powi(3);
+    if !(estimate <= 8.0 * opts.max_elements as f64) {
+        return Err(Error::Limit(format!("about {estimate:.3e} hexahedra at this size; the limit is {}", opts.max_elements)));
+    }
+    Ok(())
+}
+
 /// The frames a layered method can try, fewest levels first, or the first reason none fits.
 fn candidate_setups(ctx: &Ctx, method: HexMethod) -> std::result::Result<Vec<layered::Setup>, String> {
     let frames = if method == HexMethod::Extrude { layered::translations(ctx.brep) } else { layered::rotations(ctx.brep) };
@@ -292,7 +307,7 @@ fn candidate_setups(ctx: &Ctx, method: HexMethod) -> std::result::Result<Vec<lay
 fn with_refinement(ctx: &Ctx, method: HexMethod, run: impl Fn(u32) -> Result<Raw>) -> Result<HexMesh> {
     let mut last = Error::Failed("no attempt".into());
     for refine in 0..=MAX_REFINE {
-        match run(refine).and_then(|raw| finish(ctx, raw, method)) {
+        match run(refine).and_then(|raw| finish(&Cells::one(ctx), raw, method)) {
             Ok(mesh) => return Ok(mesh),
             Err(e @ Error::Failed(_)) => last = e,
             Err(e) => return Err(e),
@@ -314,14 +329,130 @@ fn hex_face(k: usize) -> [usize; 4] {
     ElementKind::Hex8.faces().get(k).and_then(|f| <[usize; 4]>::try_from(*f).ok()).unwrap_or([0, 1, 2, 3])
 }
 
-/// Tag, check, upgrade to hex20 when asked, gate and measure.
-fn finish(ctx: &Ctx, raw: Raw, method: HexMethod) -> Result<HexMesh> {
+/// The cells a mesh spans, for [`finish`]: one context per cell and the cell of each element.
+pub(crate) struct Cells<'c, 'a> {
+    pub ctxs: &'c [Ctx<'a>],
+    /// The cell of each hexahedron (index into `ctxs`); empty: all in the first.
+    pub of: &'c [u32],
+    /// The cells' names, for messages.
+    pub names: &'c [&'c str],
+}
+
+impl<'c, 'a> Cells<'c, 'a> {
+    pub(crate) fn one(ctx: &'c Ctx<'a>) -> Self {
+        Cells { ctxs: std::slice::from_ref(ctx), of: &[], names: &[] }
+    }
+
+    fn name(&self, c: usize) -> String {
+        self.names.get(c).map_or_else(|| c.to_string(), |n| format!("`{n}`"))
+    }
+
+    fn cell(&self, e: usize) -> usize {
+        self.of.get(e).map_or(0, |&c| c as usize)
+    }
+
+    fn ctx(&self, e: usize) -> Result<&Ctx<'a>> {
+        self.ctxs.get(self.cell(e)).ok_or_else(|| Error::Failed(format!("element {e} has no cell")))
+    }
+}
+
+/// Buckets of a uniform grid over the cells' bounding boxes: which cells may be near a point.
+struct CellGrid {
+    lo: P,
+    step: f64,
+    dims: [usize; 3],
+    buckets: Vec<Vec<usize>>,
+}
+
+/// Cosine of the largest angle between a boundary face's normal and the reversed normal of
+/// another cell's triangle for the face to lie on that triangle (about 25°): a face on another
+/// cell faces it, while faces that merely meet at an edge (a fillet beside a split) do not.
+const FACING: f64 = 0.9;
+
+/// Buckets per axis of a [`CellGrid`] at most.
+const MAX_GRID_DIM: usize = 64;
+
+impl CellGrid {
+    /// The cells' boxes grown by `gap`.
+    fn new(ctxs: &[Ctx], gap: f64) -> CellGrid {
+        let boxes: Vec<(P, P)> = ctxs
+            .iter()
+            .map(|c| {
+                c.brep.surface.points().iter().fold(([f64::INFINITY; 3], [f64::NEG_INFINITY; 3]), |(mut lo, mut hi), p| {
+                    for k in 0..3 {
+                        lo[k] = lo[k].min(p[k] - gap);
+                        hi[k] = hi[k].max(p[k] + gap);
+                    }
+                    (lo, hi)
+                })
+            })
+            .collect();
+        let (lo, hi) = boxes.iter().fold(([f64::INFINITY; 3], [f64::NEG_INFINITY; 3]), |(mut lo, mut hi), (a, b)| {
+            for k in 0..3 {
+                lo[k] = lo[k].min(a[k]);
+                hi[k] = hi[k].max(b[k]);
+            }
+            (lo, hi)
+        });
+        let per_axis = ((ctxs.len() as f64).cbrt().ceil() as usize * 2).clamp(1, MAX_GRID_DIM);
+        let extent = (0..3).map(|k| hi[k] - lo[k]).fold(0.0, f64::max);
+        let step = if extent > 0.0 && extent.is_finite() { extent / per_axis as f64 } else { 1.0 };
+        let dim = |k: usize| (((hi[k] - lo[k]) / step).ceil() as usize).clamp(1, MAX_GRID_DIM);
+        let dims = [dim(0), dim(1), dim(2)];
+        let mut grid = CellGrid { lo, step, dims, buckets: vec![Vec::new(); dims[0] * dims[1] * dims[2]] };
+        for (c, (a, b)) in boxes.iter().enumerate() {
+            let (i0, i1) = (grid.index(*a), grid.index(*b));
+            for i in i0[0]..=i1[0] {
+                for j in i0[1]..=i1[1] {
+                    for k in i0[2]..=i1[2] {
+                        if let Some(bucket) = grid.buckets.get_mut((k * dims[1] + j) * dims[0] + i) {
+                            bucket.push(c);
+                        }
+                    }
+                }
+            }
+        }
+        grid
+    }
+
+    fn index(&self, p: P) -> [usize; 3] {
+        [0, 1, 2].map(|k| {
+            let x = ((p[k] - self.lo[k]) / self.step).floor();
+            if x.is_finite() && x > 0.0 { (x as usize).min(self.dims[k] - 1) } else { 0 }
+        })
+    }
+
+    fn near(&self, p: P) -> &[usize] {
+        let [i, j, k] = self.index(p);
+        self.buckets.get((k * self.dims[1] + j) * self.dims[0] + i).map(Vec::as_slice).unwrap_or(&[])
+    }
+}
+
+/// A boundary quadrilateral of a mesh being finished: the element face, its cell and the B-rep
+/// face (index into the cell's faces) it lies on.
+struct Skin {
+    element: u32,
+    face: u8,
+    cell: usize,
+    brep: usize,
+}
+
+/// Tag, check, upgrade to hex20 when asked, gate and measure. Besides the gate, the tiling and
+/// the fidelity checks, no boundary face of a mesh of several cells may lie on another cell's
+/// boundary (a free internal face: cells that do not conform there), and the mesh boundary must
+/// be a closed 2-manifold of quadrilaterals (every edge on exactly two).
+fn finish(cells: &Cells, raw: Raw, method: HexMethod) -> Result<HexMesh> {
     let Raw { mut nodes, hexes, description } = raw;
+    let first = cells.ctxs.first().ok_or_else(|| Error::Failed("no cell to finish".into()))?;
+    let opts = first.opts;
     if hexes.is_empty() {
         return Err(Error::Failed("the method produced no elements".into()));
     }
-    if hexes.len() > ctx.opts.max_elements {
-        return Err(Error::Limit(format!("the mesh has {} hexahedra; the limit is {}", hexes.len(), ctx.opts.max_elements)));
+    if hexes.len() > opts.max_elements {
+        return Err(Error::Limit(format!("the mesh has {} hexahedra; the limit is {}", hexes.len(), opts.max_elements)));
+    }
+    if !cells.of.is_empty() && cells.of.len() != hexes.len() {
+        return Err(Error::Failed("a cell is needed for every element".into()));
     }
     fieldcraft_model::check_coordinates(nodes.as_flattened())?;
     let at = |nodes: &[P], i: u32| nodes.get(i as usize).copied().ok_or_else(|| Error::Failed(format!("node {i} out of range")));
@@ -337,8 +468,14 @@ fn finish(ctx: &Ctx, raw: Raw, method: HexMethod) -> Result<HexMesh> {
     if let Some(n) = uses.values().copied().max().filter(|&n| n > 2) {
         return Err(Error::Failed(format!("a hexahedron face is shared by {n} elements (the mesh is not a manifold)")));
     }
+    // A boundary face this close to another cell's boundary lies on it.
+    let gap = 2.0 * opts.chordal_tolerance + cells.ctxs.iter().map(|c| c.tol).fold(0.0, f64::max);
+    let others = (cells.ctxs.len() > 1).then(|| CellGrid::new(cells.ctxs, gap));
+    let mut skin: Vec<Skin> = Vec::new();
     let mut boundary = Vec::new();
     let mut fidelity = 0.0;
+    // Boundary edges of the quadrilaterals: each must be on exactly two of them.
+    let mut rim: HashMap<(u32, u32), u32> = HashMap::new();
     for (e, h) in hexes.iter().enumerate() {
         for k in 0..6 {
             let corners = hex_face(k).map(|i| h.get(i).copied().unwrap_or(0));
@@ -347,6 +484,11 @@ fn finish(ctx: &Ctx, raw: Raw, method: HexMethod) -> Result<HexMesh> {
             if uses.get(&key) != Some(&1) {
                 continue;
             }
+            for m in 0..4 {
+                let (a, b) = (corners[m], corners[(m + 1) % 4]);
+                *rim.entry((a.min(b), a.max(b))).or_insert(0) += 1;
+            }
+            let ctx = cells.ctx(e)?;
             let [a, b, c, d] = [at(&nodes, corners[0])?, at(&nodes, corners[1])?, at(&nodes, corners[2])?, at(&nodes, corners[3])?];
             let centre = scale(add(add(a, b), add(c, d)), 0.25);
             let n = cross(sub(c, a), sub(d, b));
@@ -357,14 +499,35 @@ fn finish(ctx: &Ctx, raw: Raw, method: HexMethod) -> Result<HexMesh> {
                 .closest(centre, &|t| dot(ctx.tree.normal(t), n) > 0.0)
                 .ok_or_else(|| Error::Failed(format!("no boundary face faces the same way as element {e}'s face {}", k + 1)))?;
             let f = ctx.brep.face_of_triangle(hit.triangle).ok_or_else(|| Error::Failed("triangle without a face".into()))?;
-            fidelity += area * (hit.dist + ctx.opts.chordal_tolerance);
-            boundary.push(BoundaryFace { element: e as u32, face: k as u8 + 1, brep_face: ctx.brep.face_tag(f) });
+            let cell = cells.cell(e);
+            if let Some(grid) = &others {
+                for &o in grid.near(centre) {
+                    let Some(other) = cells.ctxs.get(o).filter(|_| o != cell) else { continue };
+                    if other.tree.closest_within(centre, gap, &|t| dot(other.tree.normal(t), n) < -FACING).is_some() {
+                        return Err(Error::Failed(format!(
+                            "element {e}'s face {} lies on cell {}'s boundary without a neighbour there: cells {} and {} do not conform (split them so that they meet face to face)",
+                            k + 1,
+                            cells.name(o),
+                            cells.name(cell),
+                            cells.name(o)
+                        )));
+                    }
+                }
+            }
+            // The geometry's own face tag (a member of a composite face keeps its own).
+            let tag = ctx.brep.surface.tri_face().get(hit.triangle as usize).copied().unwrap_or_else(|| ctx.brep.face_tag(f));
+            fidelity += area * (hit.dist + opts.chordal_tolerance);
+            boundary.push(BoundaryFace { element: e as u32, face: k as u8 + 1, brep_face: tag });
+            skin.push(Skin { element: e as u32, face: k as u8 + 1, cell, brep: f });
         }
     }
-    let kind = if ctx.opts.quadratic { ElementKind::Hex20 } else { ElementKind::Hex8 };
+    if let Some((&(a, b), &n)) = rim.iter().filter(|(_, n)| **n != 2).min_by_key(|(k, _)| **k) {
+        return Err(Error::Failed(format!("the mesh boundary is not closed: the edge between nodes {a} and {b} is on {n} boundary faces")));
+    }
+    let kind = if opts.quadratic { ElementKind::Hex20 } else { ElementKind::Hex8 };
     let mut connectivity: Vec<u32> = Vec::with_capacity(hexes.len() * kind.node_count());
-    if ctx.opts.quadratic {
-        let mids = mid_edge_nodes(ctx, &mut nodes, &hexes, &boundary)?;
+    if opts.quadratic {
+        let mids = mid_edge_nodes(cells, &mut nodes, &hexes, &skin)?;
         for (h, m) in hexes.iter().zip(&mids) {
             connectivity.extend_from_slice(h);
             connectivity.extend_from_slice(m);
@@ -398,7 +561,7 @@ fn finish(ctx: &Ctx, raw: Raw, method: HexMethod) -> Result<HexMesh> {
     let gate = validity_gate_checks(checks.iter().map(|(l, c)| (*l, c)), &GateOptions { expected_volume: Some(enclosed), volume_tolerance: 1e-9 })?
         .into_result()
         .map_err(|e| Error::Failed(format!("{} ({description})", plain(&e))))?;
-    let surface_volume = ctx.brep.surface.volume();
+    let surface_volume: f64 = cells.ctxs.iter().map(|c| c.brep.surface.volume()).sum();
     let slack = 2.0 * fidelity + 1e-9 * surface_volume.abs();
     if (gate.volume - surface_volume).abs() > slack {
         return Err(Error::Failed(format!(
@@ -478,46 +641,52 @@ pub(crate) fn project_to_face(ctx: &Ctx, f: usize, p: P) -> Option<P> {
         return Some(q);
     }
     let tag = face.id;
-    let tri_face = ctx.brep.surface.tri_face();
-    ctx.tree.closest(p, &|t| tri_face.get(t as usize) == Some(&tag)).map(|h| h.point)
+    ctx.tree.closest(p, &|t| ctx.brep.tag(t) == Some(tag)).map(|h| h.point)
 }
 
 /// The point on every listed face (their common curve or corner) nearest `p`, by alternating
 /// projections.
 pub(crate) fn project_to_faces(ctx: &Ctx, faces: &[usize], p: P) -> P {
+    let pairs: Vec<(&Ctx, usize)> = faces.iter().map(|&f| (ctx, f)).collect();
+    project_alternating(&pairs, ctx.tol, p)
+}
+
+/// The point on every listed face of the listed cells nearest `p`, by alternating projections
+/// until a round moves it less than a thousandth of `tol`.
+fn project_alternating(faces: &[(&Ctx, usize)], tol: f64, p: P) -> P {
     let mut x = p;
     for _ in 0..PROJECTION_ROUNDS {
         let before = x;
-        for &f in faces {
+        for &(ctx, f) in faces {
             if let Some(q) = project_to_face(ctx, f, x) {
                 x = q;
             }
         }
-        if faces.len() < 2 || dist(before, x) <= ctx.tol * 1e-3 {
+        if faces.len() < 2 || dist(before, x) <= tol * 1e-3 {
             break;
         }
     }
     x
 }
 
-/// Mid-edge nodes of every hexahedron (Hex20 order: positions 8..20), shared between elements;
-/// boundary edges get theirs on the B-rep faces of the boundary quads around them.
-fn mid_edge_nodes(ctx: &Ctx, nodes: &mut Vec<P>, hexes: &[[u32; 8]], boundary: &[BoundaryFace]) -> Result<Vec<[u32; 12]>> {
-    // B-rep faces at each boundary edge.
-    let mut edge_faces: HashMap<(u32, u32), Vec<usize>> = HashMap::new();
-    for b in boundary {
+/// Mid-edge nodes of every hexahedron (Hex20 order: positions 8..20), shared between elements
+/// (of any cell); boundary edges get theirs on the B-rep faces of the boundary quads around them.
+fn mid_edge_nodes(cells: &Cells, nodes: &mut Vec<P>, hexes: &[[u32; 8]], skin: &[Skin]) -> Result<Vec<[u32; 12]>> {
+    // B-rep faces (cell, face) at each boundary edge.
+    let mut edge_faces: HashMap<(u32, u32), Vec<(usize, usize)>> = HashMap::new();
+    for b in skin {
         let Some(h) = hexes.get(b.element as usize) else { continue };
         let c = hex_face(b.face as usize - 1).map(|i| h.get(i).copied().unwrap_or(0));
-        let f = ctx.brep.faces.binary_search_by_key(&b.brep_face, |x| x.id).map_err(|_| Error::Failed("boundary tag".into()))?;
         for k in 0..4 {
             let (a, z) = (c[k], c[(k + 1) % 4]);
             let list = edge_faces.entry((a.min(z), a.max(z))).or_default();
-            if !list.contains(&f) {
-                list.push(f);
+            if !list.contains(&(b.cell, b.brep)) {
+                list.push((b.cell, b.brep));
             }
         }
     }
     let edges = ElementKind::Hex20.edges();
+    let tol = cells.ctxs.iter().map(|c| c.tol).fold(f64::INFINITY, f64::min);
     let mut mid: HashMap<(u32, u32), u32> = HashMap::new();
     let mut out = Vec::with_capacity(hexes.len());
     for h in hexes {
@@ -533,7 +702,10 @@ fn mid_edge_nodes(ctx: &Ctx, nodes: &mut Vec<P>, hexes: &[[u32; 8]], boundary: &
                     let (a, b) = (nodes.get(i as usize).copied().unwrap_or([0.0; 3]), nodes.get(j as usize).copied().unwrap_or([0.0; 3]));
                     let p = scale(add(a, b), 0.5);
                     let p = match edge_faces.get(&key) {
-                        Some(fs) => project_to_faces(ctx, fs, p),
+                        Some(fs) => {
+                            let pairs: Vec<(&Ctx, usize)> = fs.iter().filter_map(|&(c, f)| cells.ctxs.get(c).map(|ctx| (ctx, f))).collect();
+                            project_alternating(&pairs, tol, p)
+                        }
                         None => p,
                     };
                     nodes.push(p);

crates/mesh/src/hex/tests.rs

@@ -292,3 +292,191 @@ fn unsweepable_bodies_say_why_and_hostile_input_is_an_error() {
     let m = hex_mesh(&s, &garbage, &UniformSize(1.0), &HexOptions::default()).unwrap();
     assert_eq!(m.element_count(), 24);
 }
+
+// ---- Conforming cells ----
+
+/// A prism over a counter-clockwise polygon from z0 to z1 whose side from vertex k to k + 1 gets
+/// tag `side_tags[k]` (bottom 0, top 1); quads of one tag make one face.
+fn tagged_prism(poly: &[[f64; 2]], side_tags: &[u32], z0: f64, z1: f64) -> BoundarySurface {
+    let n = poly.len();
+    let mut pts: Vec<[f64; 3]> = poly.iter().map(|p| [p[0], p[1], z0]).collect();
+    pts.extend(poly.iter().map(|p| [p[0], p[1], z1]));
+    let mut tris = Vec::new();
+    let mut tags = Vec::new();
+    for t in ears(&pts, &(0..n as u32).rev().collect::<Vec<_>>()) {
+        tris.push(t);
+        tags.push(0);
+    }
+    for t in ears(&pts, &(n as u32..2 * n as u32).collect::<Vec<_>>()) {
+        tris.push(t);
+        tags.push(1);
+    }
+    for (k, &tag) in side_tags.iter().enumerate().take(n) {
+        let (a, b) = (k as u32, ((k + 1) % n) as u32);
+        let (c, d) = (b + n as u32, a + n as u32);
+        tris.push([a, b, c]);
+        tris.push([a, c, d]);
+        tags.extend([tag, tag]);
+    }
+    BoundarySurface::new(&pts, &tris, &tags, &SurfaceOptions::default()).unwrap()
+}
+
+fn cells_of<'a>(surfaces: &'a [BoundarySurface], faces: &'a [Vec<Option<Surface>>]) -> Vec<Cell<'a>> {
+    surfaces.iter().zip(faces).map(|(surface, faces)| Cell { name: "cell", surface, faces }).collect()
+}
+
+/// Boxes `[x0, x1] × [y0, y1] × [0, h]` as prisms (faces 0 bottom, 1 top, 2..6 the sides).
+fn boxes(spans: &[([f64; 2], [f64; 2])], h: f64) -> Vec<BoundarySurface> {
+    spans.iter().map(|(x, y)| prism(&[[x[0], y[0]], [x[1], y[0]], [x[1], y[1]], [x[0], y[1]]], 0.0, h)).collect()
+}
+
+/// Every node of the mesh at a distinct position: the cells' nodes were merged.
+fn distinct_positions(m: &HexMesh) -> bool {
+    let mut keys: Vec<[u64; 3]> = m.nodes.iter().map(|p| p.map(|x| (x + 0.0).to_bits())).collect();
+    keys.sort_unstable();
+    keys.windows(2).all(|w| w[0] != w[1])
+}
+
+#[test]
+fn boxes_side_by_side_conform() {
+    // Two boxes sharing x = 2, then a 2 × 2 arrangement whose middle edge four cells share.
+    let pairs = boxes(&[([0.0, 2.0], [0.0, 1.0]), ([2.0, 3.0], [0.0, 1.0])], 1.0);
+    let none = vec![Vec::new(); 2];
+    let m = hex_mesh_cells(&cells_of(&pairs, &none), &UniformSize(0.5), &HexOptions::default()).unwrap();
+    assert_eq!(m.interfaces, 1);
+    assert_eq!(m.blocks.iter().map(|b| b[0] * b[1] * b[2]).sum::<usize>(), m.mesh.element_count());
+    assert_eq!(m.mesh.element_count(), 6 * 2 * 2);
+    // One structured grid: 7 × 3 × 3 nodes, none duplicated on the interface.
+    assert_eq!(m.mesh.nodes.len(), 7 * 3 * 3);
+    assert!(distinct_positions(&m.mesh));
+    assert!((m.mesh.report.volume - 3.0).abs() < 1e-12);
+    assert!(m.mesh.report.gate.passed());
+    // The boundary is the outside only: 2·(6·2 + 6·2 + 2·2) quads.
+    assert_eq!(m.mesh.boundary.len(), 56);
+    assert_eq!(m.element_cell.iter().filter(|&&c| c == 0).count(), 16);
+
+    let grid = boxes(&[([0.0, 1.0], [0.0, 1.0]), ([1.0, 2.0], [0.0, 1.0]), ([0.0, 1.0], [1.0, 2.0]), ([1.0, 2.0], [1.0, 2.0])], 1.0);
+    let none = vec![Vec::new(); 4];
+    let a = hex_mesh_cells(&cells_of(&grid, &none), &UniformSize(0.25), &HexOptions::default()).unwrap();
+    assert_eq!(a.interfaces, 4);
+    assert_eq!(a.mesh.nodes.len(), 9 * 9 * 5);
+    assert!(distinct_positions(&a.mesh));
+    // Deterministic, and hex20 shares its mid-edge nodes across the interfaces too.
+    let b = hex_mesh_cells(&cells_of(&grid, &none), &UniformSize(0.25), &HexOptions::default()).unwrap();
+    assert_eq!(a, b);
+    let q = hex_mesh_cells(&cells_of(&grid, &none), &UniformSize(0.25), &HexOptions { quadratic: true, ..HexOptions::default() }).unwrap();
+    assert_eq!(q.mesh.kind, ElementKind::Hex20);
+    // Corner nodes plus one node per lattice edge: 9·9·5 + (8·9·5)·2 + 9·9·4.
+    assert_eq!(q.mesh.nodes.len(), 405 + 720 + 324);
+    assert!(distinct_positions(&q.mesh));
+}
+
+#[test]
+fn cells_that_do_not_meet_face_to_face_are_refused() {
+    // One box against two half-height boxes: its face at x = 1 meets two faces (a T-junction).
+    let t = vec![
+        prism(&[[0.0, 0.0], [1.0, 0.0], [1.0, 2.0], [0.0, 2.0]], 0.0, 1.0),
+        prism(&[[1.0, 0.0], [2.0, 0.0], [2.0, 1.0], [1.0, 1.0]], 0.0, 1.0),
+        prism(&[[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0]], 0.0, 1.0),
+    ];
+    let none = vec![Vec::new(); 3];
+    match hex_mesh_cells(&cells_of(&t, &none), &UniformSize(0.5), &HexOptions::default()) {
+        Err(Error::Failed(why)) => assert!(why.contains("interface") || why.contains("conform"), "{why}"),
+        other => panic!("{other:?}"),
+    }
+    // A cell that is no logical cube says which.
+    let l = vec![prism(&[[0.0, 0.0], [6.0, 0.0], [6.0, 2.0], [2.0, 2.0], [2.0, 6.0], [0.0, 6.0]], 0.0, 2.0)];
+    match hex_mesh_cells(&cells_of(&l, &[Vec::new()]), &UniformSize(1.0), &HexOptions::default()) {
+        Err(Error::Unsupported(why)) => assert!(why.starts_with("cell `cell`: not a logical cube"), "{why}"),
+        other => panic!("{other:?}"),
+    }
+}
+
+#[test]
+fn hostile_cell_input_is_an_error() {
+    let two = boxes(&[([0.0, 1.0], [0.0, 1.0]), ([1.0, 2.0], [0.0, 1.0])], 1.0);
+    let none = vec![Vec::new(); 2];
+    let cells = cells_of(&two, &none);
+    assert!(hex_mesh_cells(&[], &UniformSize(1.0), &HexOptions::default()).is_err());
+    for method in [HexMethod::Sweep, HexMethod::Extrude, HexMethod::Revolve, HexMethod::Swept] {
+        assert!(matches!(hex_mesh_cells(&cells, &UniformSize(1.0), &HexOptions { method, ..HexOptions::default() }), Err(Error::Invalid(_))));
+    }
+    let bad_sizes: [&dyn crate::size::SizeField; 3] = [&UniformSize(f64::NAN), &UniformSize(0.0), &UniformSize(-1.0)];
+    for f in bad_sizes {
+        assert!(hex_mesh_cells(&cells, f, &HexOptions::default()).is_err());
+    }
+    assert!(matches!(hex_mesh_cells(&cells, &UniformSize(1e-4), &HexOptions::default()), Err(Error::Limit(_))));
+    assert!(matches!(hex_mesh_cells(&cells, &UniformSize(0.1), &HexOptions { max_elements: 100, ..HexOptions::default() }), Err(Error::Limit(_))));
+    assert!(hex_mesh_cells(&cells, &UniformSize(1.0), &HexOptions { chordal_tolerance: -1.0, ..HexOptions::default() }).is_err());
+    // The same cell twice: two cells on the same side of every face.
+    let same = [two[0].clone(), two[0].clone()];
+    assert!(hex_mesh_cells(&cells_of(&same, &none), &UniformSize(0.5), &HexOptions::default()).is_err());
+    // Garbage analytic surfaces are not trusted.
+    let garbage = vec![vec![Some(Surface::Cylinder { origin: [f64::NAN; 3], axis: [0.0; 3], radius: -1.0 }); 6]; 2];
+    let m = hex_mesh_cells(&cells_of(&two, &garbage), &UniformSize(0.5), &HexOptions::default()).unwrap();
+    assert_eq!(m.mesh.element_count(), 16);
+}
+
+/// The four cells of an O-grid round a hole of radius 1 in the square [−2, 2]², height `h`: each
+/// cell the hole's quarter between two diagonals and a side of the square. The hole's quarter is
+/// one face (tag 5) or two halves of the cylinder (tags 5 and 6: virtual topology joins them).
+fn o_grid(h: f64, arc_segments: usize, split_arcs: bool) -> (Vec<BoundarySurface>, Vec<Vec<Option<Surface>>>) {
+    use std::f64::consts::FRAC_PI_4;
+    let mut surfaces = Vec::new();
+    let mut faces = Vec::new();
+    for q in 0..4 {
+        let turn = q as f64 * 2.0 * FRAC_PI_4;
+        let rot = |p: [f64; 2]| [p[0] * turn.cos() - p[1] * turn.sin(), p[0] * turn.sin() + p[1] * turn.cos()];
+        let mut poly = vec![rot([FRAC_PI_4.cos(), -FRAC_PI_4.sin()]), rot([2.0, -2.0]), rot([2.0, 2.0])];
+        let mut tags = vec![2, 3, 4];
+        for k in 0..arc_segments {
+            let t = FRAC_PI_4 - 2.0 * FRAC_PI_4 * k as f64 / arc_segments as f64;
+            poly.push(rot([t.cos(), t.sin()]));
+            tags.push(if split_arcs && 2 * k >= arc_segments { 6 } else { 5 });
+        }
+        // Exact corners where the cells meet.
+        let snap = |p: [f64; 2]| p.map(|x| if (x.abs() - 2.0).abs() < 1e-12 { 2.0 * x.signum() } else { x });
+        let poly: Vec<[f64; 2]> = poly.into_iter().map(snap).collect();
+        surfaces.push(tagged_prism(&poly, &tags, 0.0, h));
+        let cyl = Some(Surface::Cylinder { origin: [0.0; 3], axis: [0.0, 0.0, 1.0], radius: 1.0 });
+        faces.push(vec![None, None, None, None, None, cyl, cyl]);
+    }
+    (surfaces, faces)
+}
+
+#[test]
+fn an_o_grid_round_a_hole_conforms_with_graded_radial_lines() {
+    for split_arcs in [false, true] {
+        let (surfaces, faces) = o_grid(0.5, 24, split_arcs);
+        let opts = HexOptions { chordal_tolerance: 1.0 - (std::f64::consts::PI / 48.0).cos(), ..HexOptions::default() };
+        let m = hex_mesh_cells(&cells_of(&surfaces, &faces), &UniformSize(0.5), &opts).unwrap();
+        assert_eq!(m.interfaces, 4, "{split_arcs}");
+        assert!(m.mesh.report.gate.passed());
+        assert!(m.mesh.report.min_scaled_jacobian > 0.3, "{}", m.mesh.report.min_scaled_jacobian);
+        assert!(distinct_positions(&m.mesh));
+        // A ring of quads round the hole: every block is (around) × (radial) × (through).
+        let sorted = |b: [usize; 3]| {
+            let mut b = b;
+            b.sort_unstable();
+            b
+        };
+        assert!(m.blocks.iter().all(|&b| sorted(b) == sorted(m.blocks[0])), "{:?}", m.blocks);
+        let per_block = m.blocks[0].iter().product::<usize>();
+        assert_eq!(m.mesh.element_count(), 4 * per_block);
+        // Nodes on the hole lie on it.
+        let hole: Vec<u32> = if split_arcs { vec![5, 6] } else { vec![5] };
+        assert!(on_radius(&m.mesh, &hole, 1.0), "{split_arcs}");
+        let on_hole = m.mesh.boundary.iter().filter(|b| hole.contains(&b.brep_face)).count();
+        assert!(on_hole > 0);
+        // Radial grading: the first spacing off the hole is smaller than the last at the square.
+        let radii: std::collections::BTreeSet<u64> =
+            m.mesh.nodes.iter().filter(|p| p[2] == 0.0 && p[1].abs() < 1e-12 && p[0] > 0.0).map(|p| p[0].to_bits()).collect();
+        let r: Vec<f64> = radii.into_iter().map(f64::from_bits).collect();
+        assert!(r.len() >= 3, "{r:?}");
+        let first = r[1] - r[0];
+        let last = r[r.len() - 1] - r[r.len() - 2];
+        assert!(first < 0.75 * last, "{r:?}");
+        // Spacing grows by at most the growth rate plus the rounding of the count.
+        assert!(r.windows(3).all(|w| (w[2] - w[1]) <= 1.35 * (w[1] - w[0])), "{r:?}");
+    }
+}

crates/mesh/src/quad/tfi.rs

@@ -21,7 +21,7 @@ fn dist<const D: usize>(a: &[f64; D], b: &[f64; D]) -> f64 {
 }
 
 /// Normalized cumulative arc length of a side, from 0 to 1.
-fn params<const D: usize>(side: &[[f64; D]], name: &str) -> Result<Vec<f64>> {
+pub(crate) fn params<const D: usize>(side: &[[f64; D]], name: &str) -> Result<Vec<f64>> {
     let mut acc = Vec::with_capacity(side.len());
     let mut s = 0.0;
     acc.push(0.0);

crates/render/src/display.rs

@@ -55,6 +55,8 @@ pub struct DisplayNode {
 pub struct SurfaceTri {
     pub p: [Vec3; 3],
     pub normal: Vec3,
+    /// Which surface (body) it belongs to: the order [`DisplayMesh::add_surface`] got them in.
+    pub body: u32,
 }
 
 /// Display data for a model and its geometry.
@@ -72,6 +74,8 @@ pub struct DisplayMesh {
     pub surface_tris: Vec<SurfaceTri>,
     /// Geometry edges: where two faces meet, and open boundaries.
     pub surface_edges: Vec<[Vec3; 2]>,
+    /// Geometry surfaces (bodies) added.
+    pub bodies: u32,
     pub bounds: Aabb3,
 }
 
@@ -225,7 +229,7 @@ impl DisplayMesh {
             for p in [a, b, c] {
                 self.bounds.add(p);
             }
-            self.surface_tris.push(SurfaceTri { p: [a, b, c], normal });
+            self.surface_tris.push(SurfaceTri { p: [a, b, c], normal, body: self.bodies });
             let face = tri_face.get(t).copied().unwrap_or(u32::MAX);
             let [v0, v1, v2] = [vertex(i0), vertex(i1), vertex(i2)];
             for (u, v) in [(v0, v1), (v1, v2), (v2, v0)] {
@@ -241,6 +245,17 @@ impl DisplayMesh {
                 self.surface_edges.push([a, b]);
             }
         }
+        self.bodies = self.bodies.saturating_add(1);
+    }
+
+    /// The colour of a geometry triangle: `single` when there is one body, else body `k`'s
+    /// colour from `palette` (cycling), so the bodies of a design, such as the cells of a
+    /// partitioned part, tell apart.
+    pub fn body_color<C: Copy>(&self, tri: &SurfaceTri, single: C, palette: &[C]) -> C {
+        if self.bodies <= 1 || palette.is_empty() {
+            return single;
+        }
+        palette.get(tri.body as usize % palette.len()).copied().unwrap_or(single)
     }
 
     /// Where node `label` is drawn, if it is (nodes are kept in label order).
@@ -401,6 +416,13 @@ mod tests {
         welded.add_surface(&soup, &[[0, 1, 2], [3, 4, 5]], &[0, 0]);
         assert_eq!(welded.surface_edges.len(), 4);
         let mut hostile = DisplayMesh::default();
+        let mut two = DisplayMesh::default();
+        two.add_surface(&p, &[[0, 1, 2]], &[0]);
+        two.add_surface(&p, &[[0, 2, 3]], &[0]);
+        assert_eq!(two.bodies, 2);
+        let colors: Vec<u8> = two.surface_tris.iter().map(|t| two.body_color(t, 0u8, &[7, 8])).collect();
+        assert_eq!(colors, vec![7, 8]);
+        assert_eq!(one_face.body_color(&one_face.surface_tris[0], 0u8, &[7, 8]), 0);
         hostile.add_surface(&[[f64::NAN, 0., 0.], [1., 0., 0.], [0., 1., 0.]], &[[0, 1, 2], [u32::MAX, 0, 1]], &[]);
         assert!(hostile.surface_tris.is_empty() && hostile.surface_edges.is_empty());
     }

crates/testkit/src/corpus.rs

@@ -40,6 +40,7 @@ pub const CATEGORIES: &[(&str, &str)] = &[
     ("blocky", "axis-aligned blocks (polycube-friendly)"),
     ("organic", "free-form, lofted or revolved spline surfaces"),
     ("assembly", "several bodies placed by component transforms"),
+    ("partitioned", "split into cells that each mesh as a mapped block (conforming multi-block hex)"),
     ("analytic", "closed-form geometry for solver verification"),
 ];
 

crates/testkit/src/runner.rs

@@ -15,7 +15,7 @@ use std::time::Instant;
 use crate::corpus::{self, CATEGORIES, Part};
 use fieldcraft_geometry::{Design, FaceSurface, Measure, SurfaceMesh};
 use fieldcraft_mesh::cdt::{CdtOptions, tetrahedralize};
-use fieldcraft_mesh::hex::{HexOptions, Surface, hex_mesh};
+use fieldcraft_mesh::hex::{Cell, HexOptions, Surface, hex_mesh, hex_mesh_cells};
 use fieldcraft_mesh::size::UniformSize;
 use fieldcraft_mesh::surface::{BoundarySurface, SurfaceOptions};
 
@@ -123,7 +123,8 @@ const REASON_CHARS: usize = 240;
 
 /// Structured all-hex meshing (`fieldcraft_mesh::hex`): the method that meshed each body, the
 /// element count and the smallest scaled Jacobian of the hex8 mesh and of the hex20 mesh (mid-edge
-/// nodes on the geometry), or why no structured method fits.
+/// nodes on the geometry), or why no structured method fits. A `partitioned` part's bodies are
+/// the cells of one conforming mesh (`blocks`).
 struct StructuredHex;
 
 impl Stage for StructuredHex {
@@ -144,6 +145,12 @@ impl Stage for StructuredHex {
         let size = UniformSize(HEX_RELATIVE_SIZE * diagonal);
         let linear = HexOptions { chordal_tolerance: l.tolerance, ..HexOptions::default() };
         let quadratic = HexOptions { quadratic: true, ..linear.clone() };
+        if l.part.categories.iter().any(|c| c == "partitioned") {
+            return blocks(l, &size, &linear, &quadratic).map(|mut cells| {
+                cells.push(format!("{:.0}", start.elapsed().as_secs_f64() * 1e3));
+                cells
+            });
+        }
         let (mut methods, mut elements, mut sj8, mut sj20) = (Vec::new(), 0usize, f64::INFINITY, Ok(f64::INFINITY));
         for (body, mesh) in l.design.bodies.iter().zip(&l.surfaces) {
             let surface_options = SurfaceOptions { weld_distance: mesh.weld_distance, ..SurfaceOptions::default() };
@@ -180,6 +187,39 @@ impl Stage for StructuredHex {
     }
 }
 
+/// The bodies of a partitioned part as the cells of one conforming mesh: method, elements and
+/// the hex8 and hex20 smallest scaled Jacobians.
+fn blocks(l: &Loaded, size: &UniformSize, linear: &HexOptions, quadratic: &HexOptions) -> Result<Vec<String>, String> {
+    let mut surfaces = Vec::with_capacity(l.surfaces.len());
+    let mut faces: Vec<Vec<Option<Surface>>> = Vec::with_capacity(l.surfaces.len());
+    for (body, mesh) in l.design.bodies.iter().zip(&l.surfaces) {
+        let options = SurfaceOptions { weld_distance: mesh.weld_distance, ..SurfaceOptions::default() };
+        surfaces.push(BoundarySurface::new(&mesh.positions, &mesh.triangles, &mesh.tri_face, &options).map_err(|e| format!("{}: {e}", body.name))?);
+        faces.push(body.face_surfaces().map_err(|e| format!("{}: {e}", body.name))?.into_iter().map(|f| f.map(analytic)).collect());
+    }
+    let cells: Vec<Cell> =
+        l.design.bodies.iter().zip(&surfaces).zip(&faces).map(|((b, surface), faces)| Cell { name: &b.name, surface, faces }).collect();
+    let m = match hex_mesh_cells(&cells, size, linear) {
+        Ok(m) => m,
+        Err(e) => {
+            return Ok(vec![
+                format!("no: {}", e.to_string().chars().take(REASON_CHARS).collect::<String>()),
+                String::new(),
+                String::new(),
+                String::new(),
+            ]);
+        }
+    };
+    let sj20 =
+        hex_mesh_cells(&cells, size, quadratic).map_or_else(|e| format!("error: {e}"), |q| format!("{:.3}", q.mesh.report.min_scaled_jacobian));
+    Ok(vec![
+        format!("blocks ({} cells, {} shared faces)", cells.len(), m.interfaces),
+        m.mesh.element_count().to_string(),
+        format!("{:.3}", m.mesh.report.min_scaled_jacobian),
+        sj20,
+    ])
+}
+
 /// The geometry's analytic face surface in the mesher's terms.
 pub fn analytic(s: FaceSurface) -> Surface {
     match s {
@@ -330,6 +370,17 @@ mod tests {
         }
     }
 
+    /// The partitioned bracket's 30 cells mesh into one conforming hex mesh.
+    #[test]
+    fn partitioned_bracket_meshes_as_blocks() {
+        let part = corpus::read_recipe(&corpus::corpus_dir().join("filleted-bracket-blocks")).unwrap();
+        let (loaded, _) = load(&part).unwrap();
+        assert_eq!(loaded.design.bodies.len(), 30);
+        let cells = StructuredHex.run(&loaded).unwrap();
+        assert!(cells[0].starts_with("blocks (30 cells"), "{cells:?}");
+        assert!(cells[2].parse::<f64>().unwrap() > 0.5, "{cells:?}");
+    }
+
     /// The constrained Delaunay tetrahedralization passes the gate on parts that exercise each
     /// of its stages: an assembly whose rounding needs a few-ulp repair, cavities refilled by
     /// gift-wrapping, and a holed plate.

crates/testkit/tests/bracket_blocks.rs

@@ -0,0 +1,96 @@
+//! The filleted bracket meshed all-hex by partitioning (`tests/corpus/filleted-bracket-blocks`):
+//! the upright's end clamped, the foot's end pushed down by 1 kN (steel). The conforming hex8 and
+//! hex20 meshes converge: the tip deflection, the von Mises stress at the middle of the fillet and
+//! the peak (at the upright's hole) agree between hex20 at 3 mm and 2 mm, and hex8 approaches them.
+//! For reference, a 48 016-element tet10 mesh of the unpartitioned part (curvature-refined, from
+//! the tet quality pipeline) gives a tip deflection of 0.5093 mm, 133.2 MPa at the fillet and a
+//! peak of 291.8 MPa; hex20 at 2 mm gives 0.5089 mm, 137.4 MPa and 297.0 MPa.
+#![allow(clippy::unwrap_used, clippy::expect_used, clippy::panic)]
+
+use fieldcraft_engine::Session;
+use fieldcraft_testkit::corpus;
+use serde_json::{Value, json};
+
+fn open(part: &str) -> (Session, std::path::PathBuf) {
+    let parts = corpus::parts().unwrap();
+    let p = parts.iter().find(|p| p.name == part).unwrap();
+    let dir = std::env::temp_dir().join(format!("fieldcraft-bracket-{}-{part}", std::process::id()));
+    std::fs::create_dir_all(&dir).unwrap();
+    let path = dir.join("part.solidcraft");
+    std::fs::write(&path, corpus::design_bytes(p).unwrap()).unwrap();
+    let mut s = Session::default();
+    s.execute("geometry.open", &json!({"path": path})).unwrap();
+    (s, dir)
+}
+
+/// Every (body, face) whose name is `face` (a face split into cells keeps its name in each).
+fn faces_named(s: &mut Session, face: &str) -> Vec<Value> {
+    let bodies = s.execute("geometry.bodies", &json!({})).unwrap();
+    let mut out = Vec::new();
+    for b in bodies.as_array().unwrap() {
+        for f in b["faces"].as_array().unwrap() {
+            if f["name"] == face {
+                out.push(json!({"body": b["name"], "face": face}));
+            }
+        }
+    }
+    out
+}
+
+/// Clamp the upright's top (y = 50), push the foot's end (x = 60) down by 1 kN; tip deflection,
+/// von Mises at the middle of the fillet and the largest von Mises.
+fn analyse(s: &mut Session) -> (f64, f64, f64) {
+    let top = faces_named(s, "F2:side:l5");
+    let tip = faces_named(s, "F2:side:l2");
+    assert!(!top.is_empty() && !tip.is_empty());
+    for (name, geometry) in [
+        (
+            "all",
+            json!(
+                s.execute("geometry.bodies", &json!({})).unwrap().as_array().unwrap().iter().map(|b| json!({"body": b["name"]})).collect::<Vec<_>>()
+            ),
+        ),
+        ("top", json!(top)),
+        ("tip", json!(tip)),
+    ] {
+        s.execute("model.set.create", &json!({"name": name, "geometry": geometry})).unwrap();
+    }
+    s.execute("model.material.create", &json!({"name": "Steel", "elastic": {"type": "isotropic", "young": 210000.0, "poisson": 0.3}})).unwrap();
+    s.execute("model.section.create", &json!({"name": "solid", "set": "all", "material": "Steel", "kind": {"type": "solid"}})).unwrap();
+    s.execute("model.boundary_condition.create", &json!({"name": "clamp", "type": "fixed", "target": {"set": "top"}})).unwrap();
+    s.execute("model.load.create", &json!({"name": "push", "type": "force", "target": {"set": "tip"}, "vector": [0.0, -1000.0, 0.0]})).unwrap();
+    s.execute("model.step.create", &json!({"name": "Load", "procedure": {"type": "static"}, "boundary_conditions": ["clamp"], "loads": ["push"]}))
+        .unwrap();
+    let r = s.execute("solve.static", &json!({})).unwrap();
+    let probe = s.execute("solve.probe", &json!({"step": "Load", "at": [[60.0, 4.0, 20.0], [9.757359312880714, 9.757359312880714, 20.0]]})).unwrap();
+    let nodes = probe["nodes"].as_array().unwrap();
+    let tip = nodes[0]["displacement"][1].as_f64().unwrap();
+    let fillet = nodes[1]["von_mises"].as_f64().unwrap();
+    (tip, fillet, r["max_von_mises"]["value"].as_f64().unwrap())
+}
+
+#[test]
+#[ignore = "heavy: cargo xtask test --heavy (CI heavy job)"]
+fn partitioned_bracket_converges_with_hex8_and_hex20() {
+    let mut results = Vec::new();
+    for (order, size) in [(1, 2.0), (2, 3.0), (2, 2.0)] {
+        let (mut s, dir) = open("filleted-bracket-blocks");
+        let m = s.execute("mesh.hex.blocks", &json!({"size": size, "order": order})).unwrap();
+        assert_eq!(m["cells"], 30, "{m}");
+        assert!(m["quality"]["min_scaled_jacobian"].as_f64().unwrap() > 0.6, "{m}");
+        let r = analyse(&mut s);
+        eprintln!("order {order}, size {size}: {} elements; tip {:.5} mm, fillet {:.2} MPa, peak {:.2} MPa", m["elements"], r.0, r.1, r.2);
+        results.push(r);
+        std::fs::remove_dir_all(dir).ok();
+    }
+    let close = |a: f64, b: f64, rel: f64| (a - b).abs() <= rel * b.abs();
+    let [hex8, coarse, fine] = [results[0], results[1], results[2]];
+    // The tip moves down; hex20 has converged on it, hex8 is within a few percent (stiffer).
+    assert!(fine.0 < 0.0);
+    assert!(close(coarse.0, fine.0, 0.005), "{coarse:?} {fine:?}");
+    assert!(close(hex8.0, fine.0, 0.04) && hex8.0.abs() < fine.0.abs(), "{hex8:?} {fine:?}");
+    // Stresses at the fillet and the peak.
+    assert!(close(coarse.1, fine.1, 0.06), "{coarse:?} {fine:?}");
+    assert!(close(coarse.2, fine.2, 0.02), "{coarse:?} {fine:?}");
+    assert!(close(hex8.1, fine.1, 0.15) && close(hex8.2, fine.2, 0.10), "{hex8:?} {fine:?}");
+}

crates/ui-egui/src/scene.rs

@@ -79,7 +79,11 @@ pub fn build(display: &DisplayMesh, style: Style, show: Show, t: &Tokens) -> Sce
             s.tris = display.tris.iter().map(|d| SceneTri { p: d.p, normal: d.normal, color: rgb(t.element(d.kind)) }).collect();
         }
         if !mesh_faces && show.geometry {
-            s.tris.extend(display.surface_tris.iter().map(|d| SceneTri { p: d.p, normal: d.normal, color: rgb(t.geometry) }));
+            s.tris.extend(display.surface_tris.iter().map(|d| SceneTri {
+                p: d.p,
+                normal: d.normal,
+                color: rgb(display.body_color(d, t.geometry, &t.bodies)),
+            }));
         }
     }
     if style != Style::Shaded && show.mesh {

crates/ui-egui/src/theme.rs

@@ -66,6 +66,8 @@ pub struct Tokens {
     pub geometry: Color32,
     /// Their face boundaries.
     pub geometry_edge: Color32,
+    /// The bodies of a design with several (the cells of a partitioned part), in turn.
+    pub bodies: [Color32; 8],
     // ---- results ----
     /// The contour colour scale, low to high (the same in both themes).
     pub contour: [Color32; 5],
@@ -137,6 +139,16 @@ impl Tokens {
             node: Color32::from_rgb(0xb4, 0x3c, 0x8c),
             geometry: Color32::from_rgb(176, 186, 198),
             geometry_edge: Color32::from_rgb(40, 44, 52),
+            bodies: [
+                Color32::from_rgb(0x8f, 0xb4, 0xd9),
+                Color32::from_rgb(0xe0, 0xa4, 0x6c),
+                Color32::from_rgb(0x8c, 0xc4, 0x8a),
+                Color32::from_rgb(0xd6, 0x8c, 0x96),
+                Color32::from_rgb(0xb0, 0x9c, 0xd8),
+                Color32::from_rgb(0xd8, 0xc8, 0x70),
+                Color32::from_rgb(0x7c, 0xc4, 0xc0),
+                Color32::from_rgb(0xc0, 0xa0, 0x84),
+            ],
             contour: CONTOUR,
             no_result: Color32::from_rgb(0x9a, 0xa0, 0xaa),
             glyph_constraint: Color32::from_rgb(0x00, 0x8a, 0xa8),
@@ -188,6 +200,16 @@ impl Tokens {
         t.node = Color32::from_rgb(0xf2, 0xe6, 0x6e);
         t.geometry = Color32::from_rgb(152, 156, 162);
         t.geometry_edge = Color32::from_rgb(16, 16, 18);
+        t.bodies = [
+            Color32::from_rgb(0x6c, 0x92, 0xb8),
+            Color32::from_rgb(0xbc, 0x84, 0x50),
+            Color32::from_rgb(0x6c, 0xa2, 0x6a),
+            Color32::from_rgb(0xb4, 0x6e, 0x78),
+            Color32::from_rgb(0x90, 0x7e, 0xb6),
+            Color32::from_rgb(0xb6, 0xa6, 0x54),
+            Color32::from_rgb(0x5c, 0xa2, 0x9e),
+            Color32::from_rgb(0xa0, 0x82, 0x66),
+        ];
         t.no_result = Color32::from_rgb(0x70, 0x76, 0x80);
         t.glyph_constraint = Color32::from_rgb(0x3c, 0xdc, 0xf0);
         t.glyph_load = Color32::from_rgb(0xff, 0x5a, 0x4c);

docs/features.md

@@ -6,7 +6,7 @@ Every entry is a toolbar button. **Ready** entries run their command; **coming s
 shown greyed with the roadmap milestone that brings them (ROADMAP.md). Scripts, the command line,
 the control channel and AI agents over MCP (docs/mcp.md) run the same commands.
 
-54 ready, 49 coming soon.
+55 ready, 49 coming soon.
 
 ## Geometry
 
@@ -35,6 +35,7 @@ the control channel and AI agents over MCP (docs/mcp.md) run the same commands.
 | Structured Hex | Quad Faces — All-quad meshes of faces on shared edge nodes | coming soon (M5) |  |
 | Structured Hex | Mapped Hex — Mapped hex mesh of a six-sided block (curved faces allowed) | ready | `mesh.hex.mapped`, `mesh.hex.generate` |
 | Structured Hex | Sweep Hex — Hex mesh swept along extrusions (several levels), around revolutions, or from a source face to a target | ready | `mesh.hex.sweep`, `mesh.hex.generate` |
+| Structured Hex | Block Hex — One conforming hex mesh of a part split into six-sided blocks: shared faces meshed once | ready | `mesh.hex.blocks` |
 | Structured Hex | Decompose — Find sweepable and mappable regions and propose the splits | coming soon (M5) |  |
 | Structured Hex | Transition — Conforming tet–hex interfaces with pyramids | coming soon (M5) |  |
 | Automatic Hex | All-Hex — Octree dual all-hex mesh of any part | coming soon (M6) |  |

docs/meshing/corpus.md

@@ -10,28 +10,29 @@ clock on the machine that ran it; every other column is deterministic.
 
 | part | categories | bodies | faces | triangles | tol. mm | build ms | enclosed mm³ | SolidCraft mm³ | rel. diff | exact mm³ | boundary | genus | feature edges | check ms | CDT tets | Steiner | cavities | tet gate | CDT ms | hex | hex elements | hex8 min SJ | hex20 min SJ | hex ms |
 |---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
-| bent-tube | sweepable, thin-wall, genus | 1 | 18 | 4296 | 0.121 | 2726 | 4213.1 | 4159.4 | 1.3e-2 | 4159.8 | ok, 2148 vertices | 1 | 268 | 19.8 | 15550 | 2338 | 0 | pass | 4070 | sweep | 192 | 0.859 | 0.997 | 535 |
-| bolted-assembly | assembly, blocky | 3 | 22 | 284 | 0.117 | 1141 | 33004.7 | 32999.0 | 1.7e-4 | 32998.9 | ok, 144 vertices | 0, 1, 1 | 104 | 1.6 | 690 | 181 | 0 | pass | 211 | extrude | 357 | 0.398 | 0.222 | 139 |
-| cantilever-i-beam | sweepable, thin-wall | 1 | 14 | 44 | 0.208 | 483 | 61600.0 | 61600.0 | 0.0e0 | 61600.0 | ok, 24 vertices | 0 | 36 | 0.2 | 90 | 24 | 0 | pass | 5 | extrude | 336 | 0.707 | 0.707 | 47 |
-| cross-cylinders | boolean, sharp | 1 | 10 | 1395 | 0.087 | 59348 | 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; … |  |  |  | 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; … |  |  |  |  | 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 | 3706 | 32162.3 | 32145.0 | 5.4e-4 | 32144.8 | ok, 268 vertices | 2 | 135 | 2.1 | 5582 | 1613 | 0 | pass | 1769 | no: Bracket: not supported yet: no structured hex method fits this body (extrude: cylindrical face 8 is not parallel to the direction; sweep: from face 1: linking face 6 has 2 boundary loops (periodic or holed linking faces are not swept); from face 6:  |  |  |  | 10 |
-| i-beam | sweepable, thin-wall | 1 | 14 | 44 | 0.208 | 280 | 61600.0 | 61600.0 | 0.0e0 | 61600.0 | ok, 24 vertices | 0 | 36 | 0.3 | 94 | 24 | 0 | pass | 7 | extrude | 336 | 0.707 | 0.707 | 52 |
-| loft-transition | organic | 1 | 6 | 292 | 0.075 | 397 | 56390.2 | 56450.6 | -1.1e-3 | – | ok, 148 vertices | 0 | 40 | 1.0 | 1329 | 256 | 1 | pass | 172 | sweep | 2873 | 0.189 | 0.189 | 2056 |
-| perforated-plate | genus, thin-wall, sweepable | 1 | 54 | 1644 | 0.144 | 2864 | 26822.2 | 26764.6 | 2.2e-3 | 26764.2 | ok, 776 vertices | 24 | 828 | 10.6 | 33611 | 9798 | 0 | pass | 6513 | extrude | 467 | 0.304 | 0.033 | 480 |
-| pinhole-block | small-feature, genus | 1 | 38 | 4044 | 0.073 | 2053 | 49919.4 | 49909.6 | 2.0e-4 | 49909.5 | ok, 1992 vertices | 16 | 1548 | 41.8 | 96411 | 30221 | 0 | pass | 28180 | extrude | 5236 | 0.408 | 0.392 | 3838 |
-| rounded-box | fillet | 1 | 26 | 2540 | 0.054 | 1069 | 23325.7 | 23340.5 | -6.3e-4 | 23340.8 | ok, 1272 vertices | 0 | 216 | 11.5 | 11653 | 1914 | 0 | pass | 2265 | no: Box: not supported yet: no structured hex method fits this body (extrude: cylindrical face 8 is not parallel to the direction; sweep: from face 0: 12 faces are not on the sweep path; from face 5: 12 faces are not on the sweep path; revolve: plan |  |  |  | 67 |
-| shallow-groove | near-tangent, boolean | 1 | 8 | 62 | 0.075 | 218 | 47795.2 | 47782.3 | 2.7e-4 | 47781.5 | ok, 33 vertices | 0 | 24 | 0.2 | 267 | 87 | 0 | pass | 24 | extrude | 1960 | 0.991 | 0.986 | 1078 |
-| sharp-wedge | sharp, thin-wall, sweepable | 1 | 5 | 8 | 0.083 | 72 | 3200.0 | 3200.0 | 0.0e0 | 3200.0 | ok, 6 vertices | 0 | 9 | 0.0 | 11 | 3 | 0 | pass | 1 | extrude | 294 | 0.050 | 0.050 | 87 |
-| sphere-on-plate | near-tangent, boolean | 1 | 12 | 2826 | 0.098 | 3775 | 69353.3 | 69509.0 | -2.2e-3 | 69509.7 | ok, 1415 vertices | 0 | 124 | 8.3 | 4898 | 309 | 9 | pass | 1265 | no: Plate: not supported yet: no structured hex method fits this body (extrude: face 6 is curved across the direction; sweep: from face 0: 6 faces are not on the sweep path; from face 5: linking face 11 has 2 boundary loops (periodic or holed linking  |  |  |  | 75 |
-| spline-vase | organic | 1 | 11 | 3890 | 0.106 | 1097 | 94891.2 | 94948.5 | -6.0e-4 | – | ok, 1947 vertices | 0 | 243 | 23.5 | 12461 | 1512 | 0 | pass | 2487 | sweep | 4554 | 0.649 | 0.649 | 3030 |
-| stepped-block | blocky, sharp | 1 | 14 | 96 | 0.078 | 438 | 47000.0 | 47000.0 | 0.0e0 | 47000.0 | ok, 48 vertices | 1 | 60 | 0.6 | 223 | 41 | 0 | pass | 11 | extrude | 1416 | 0.999 | 0.999 | 906 |
-| thin-ribs | thin-wall, blocky | 1 | 26 | 92 | 0.119 | 1085 | 32000.0 | 32000.0 | -4.5e-16 | 32000.0 | ok, 48 vertices | 0 | 72 | 0.7 | 171 | 33 | 0 | pass | 10 | extrude | 777 | 0.844 | 0.844 | 236 |
-| thin-wall-enclosure | thin-wall, blocky | 1 | 11 | 36 | 0.108 | 189 | 15564.0 | 15564.0 | 0.0e0 | 15564.0 | ok, 20 vertices | 0 | 28 | 0.2 | 191 | 47 | 2 | pass | 8 | extrude | 894 | 0.297 | 0.297 | 349 |
-| analytic-sphere | analytic | 1 | 6 | 2784 | 0.069 | 1340 | 33353.4 | 33507.9 | -4.6e-3 | 33510.3 | ok, 1394 vertices | 0 | 96 | 8.9 | 4211 | 58 | 0 | pass | 1430 | no: Sphere: not supported yet: no structured hex method fits this body (extrude: no direction to try; sweep: no planar face to sweep from; revolve: no axis to try; mapped: meshing failed: meshing failed: mesh validity gate: 8 inverted or degenerate ele |  |  |  | 159 |
-| analytic-torus | analytic, genus | 1 | 6 | 7230 | 0.109 | 1725 | 37644.1 | 37893.7 | -6.6e-3 | 37899.3 | ok, 3615 vertices | 1 | 180 | 31.7 | 19790 | 466 | 0 | pass | 4946 | revolve | 1134 | 0.706 | 0.706 | 1155 |
-| analytic-cylinder | analytic, sweepable | 1 | 4 | 1004 | 0.049 | 368 | 12533.1 | 12564.7 | -2.5e-3 | 12566.4 | ok, 504 vertices | 0 | 106 | 9.3 | 7437 | 1568 | 0 | pass | 1797 | extrude | 3200 | 0.707 | 0.707 | 2316 |
-| analytic-pipe | analytic, genus, sweepable | 1 | 6 | 1448 | 0.104 | 506 | 94150.9 | 94238.0 | -9.3e-4 | 94247.8 | ok, 724 vertices | 1 | 220 | 9.5 | 6238 | 1340 | 0 | pass | 1734 | extrude | 1288 | 0.998 | 0.998 | 1059 |
-| analytic-plate-hole | analytic, genus, sweepable | 1 | 8 | 112 | 0.224 | 524 | 98449.1 | 98429.5 | 2.0e-4 | 98429.2 | ok, 56 vertices | 1 | 62 | 0.6 | 1868 | 543 | 0 | pass | 78 | extrude | 321 | 0.320 | 0.142 | 248 |
+| bent-tube | sweepable, thin-wall, genus | 1 | 18 | 4296 | 0.121 | 428 | 4213.1 | 4159.4 | 1.3e-2 | 4159.8 | ok, 2148 vertices | 1 | 268 | 11.8 | 15550 | 2338 | 0 | pass | 598 | sweep | 192 | 0.859 | 0.997 | 62 |
+| bolted-assembly | assembly, blocky | 3 | 22 | 284 | 0.117 | 53 | 33004.7 | 32999.0 | 1.7e-4 | 32998.9 | ok, 144 vertices | 0, 1, 1 | 104 | 0.7 | 690 | 181 | 0 | pass | 33 | extrude | 357 | 0.398 | 0.222 | 40 |
+| cantilever-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 | 90 | 24 | 0 | pass | 3 | extrude | 336 | 0.707 | 0.707 | 33 |
+| cross-cylinders | boolean, sharp | 1 | 10 | 1395 | 0.087 | 10640 | 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; … |  |  |  | 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; … |  |  |  |  | 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 | 494 | 32162.3 | 32145.0 | 5.4e-4 | 32144.8 | ok, 268 vertices | 2 | 135 | 1.7 | 5582 | 1613 | 0 | pass | 233 | no: Bracket: not supported yet: no structured hex method fits this body (extrude: cylindrical face 8 is not parallel to the direction; sweep: from face 1: linking face 6 has 2 boundary loops (periodic or holed linking faces are not swept); from face 6:  |  |  |  | 5 |
+| filleted-bracket-blocks | fillet, genus, partitioned | 30 | 190 | 712 | 0.088 | 1057 | 32169.8 | 32145.1 | 7.7e-4 | 32144.8 | ok, 416 vertices | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 | 474 | 1.7 | 1423 | 349 | 0 | pass | 64 | blocks (30 cells, 45 shared faces) | 4986 | 0.705 | 0.705 | 1440 |
+| i-beam | sweepable, thin-wall | 1 | 14 | 44 | 0.208 | 234 | 61600.0 | 61600.0 | 0.0e0 | 61600.0 | ok, 24 vertices | 0 | 36 | 0.3 | 94 | 24 | 0 | pass | 8 | extrude | 336 | 0.707 | 0.707 | 74 |
+| loft-transition | organic | 1 | 6 | 292 | 0.075 | 493 | 56390.2 | 56450.6 | -1.1e-3 | – | ok, 148 vertices | 0 | 40 | 1.5 | 1329 | 256 | 1 | pass | 80 | sweep | 2873 | 0.189 | 0.189 | 892 |
+| perforated-plate | genus, thin-wall, sweepable | 1 | 54 | 1644 | 0.144 | 2861 | 26822.2 | 26764.6 | 2.2e-3 | 26764.2 | ok, 776 vertices | 24 | 828 | 10.6 | 33611 | 9798 | 0 | pass | 7636 | extrude | 467 | 0.304 | 0.033 | 555 |
+| pinhole-block | small-feature, genus | 1 | 38 | 4044 | 0.073 | 4821 | 49919.4 | 49909.6 | 2.0e-4 | 49909.5 | ok, 1992 vertices | 16 | 1548 | 29.3 | 96411 | 30221 | 0 | pass | 23460 | extrude | 5236 | 0.408 | 0.392 | 1853 |
+| rounded-box | fillet | 1 | 26 | 2540 | 0.054 | 1822 | 23325.7 | 23340.5 | -6.3e-4 | 23340.8 | ok, 1272 vertices | 0 | 216 | 15.1 | 11687 | 1921 | 0 | pass | 1081 | no: Box: not supported yet: no structured hex method fits this body (extrude: cylindrical face 6 is not parallel to the direction; sweep: from face 0: 12 faces are not on the sweep path; from face 5: 12 faces are not on the sweep path; revolve: plan |  |  |  | 37 |
+| shallow-groove | near-tangent, boolean | 1 | 8 | 62 | 0.075 | 266 | 47795.2 | 47782.3 | 2.7e-4 | 47781.5 | ok, 33 vertices | 0 | 24 | 0.4 | 267 | 87 | 0 | pass | 27 | extrude | 1960 | 0.991 | 0.986 | 500 |
+| sharp-wedge | sharp, thin-wall, sweepable | 1 | 5 | 8 | 0.083 | 143 | 3200.0 | 3200.0 | 0.0e0 | 3200.0 | ok, 6 vertices | 0 | 9 | 0.0 | 11 | 3 | 0 | pass | 1 | extrude | 294 | 0.050 | 0.050 | 67 |
+| sphere-on-plate | near-tangent, boolean | 1 | 12 | 2826 | 0.098 | 3307 | 69353.3 | 69509.0 | -2.2e-3 | 69509.7 | ok, 1415 vertices | 0 | 124 | 15.4 | 4898 | 309 | 9 | pass | 736 | no: Plate: not supported yet: no structured hex method fits this body (extrude: face 6 is curved across the direction; sweep: from face 0: 6 faces are not on the sweep path; from face 5: linking face 11 has 2 boundary loops (periodic or holed linking  |  |  |  | 55 |
+| spline-vase | organic | 1 | 11 | 3890 | 0.106 | 1175 | 94891.2 | 94948.5 | -6.0e-4 | – | ok, 1947 vertices | 0 | 243 | 20.5 | 12461 | 1512 | 0 | pass | 2093 | sweep | 4554 | 0.649 | 0.649 | 846 |
+| stepped-block | blocky, sharp | 1 | 14 | 96 | 0.078 | 32 | 47000.0 | 47000.0 | 0.0e0 | 47000.0 | ok, 48 vertices | 1 | 60 | 0.5 | 223 | 41 | 0 | pass | 8 | extrude | 1416 | 0.999 | 0.999 | 158 |
+| thin-ribs | thin-wall, blocky | 1 | 26 | 92 | 0.119 | 68 | 32000.0 | 32000.0 | -4.5e-16 | 32000.0 | ok, 48 vertices | 0 | 72 | 0.6 | 171 | 33 | 0 | pass | 8 | extrude | 777 | 0.844 | 0.844 | 86 |
+| 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 | 191 | 47 | 2 | pass | 6 | extrude | 894 | 0.297 | 0.297 | 100 |
+| analytic-sphere | analytic | 1 | 6 | 2784 | 0.069 | 284 | 33353.4 | 33507.9 | -4.6e-3 | 33510.3 | ok, 1394 vertices | 0 | 96 | 7.0 | 4211 | 58 | 0 | pass | 358 | no: Sphere: not supported yet: no structured hex method fits this body (extrude: no direction to try; sweep: no planar face to sweep from; revolve: no axis to try; mapped: meshing failed: meshing failed: mesh validity gate: 8 inverted or degenerate ele |  |  |  | 35 |
+| analytic-torus | analytic, genus | 1 | 6 | 7230 | 0.109 | 627 | 37644.1 | 37893.7 | -6.6e-3 | 37899.3 | ok, 3615 vertices | 1 | 180 | 30.0 | 19790 | 466 | 0 | pass | 1091 | revolve | 1134 | 0.706 | 0.706 | 210 |
+| analytic-cylinder | analytic, sweepable | 1 | 4 | 1004 | 0.049 | 32 | 12533.1 | 12564.7 | -2.5e-3 | 12566.4 | ok, 504 vertices | 0 | 106 | 3.4 | 7437 | 1568 | 0 | pass | 344 | extrude | 3200 | 0.707 | 0.707 | 409 |
+| analytic-pipe | analytic, genus, sweepable | 1 | 6 | 1448 | 0.104 | 51 | 94150.9 | 94238.0 | -9.3e-4 | 94247.8 | ok, 724 vertices | 1 | 220 | 4.5 | 6238 | 1340 | 0 | pass | 327 | extrude | 1288 | 0.998 | 0.998 | 176 |
+| analytic-plate-hole | analytic, genus, sweepable | 1 | 8 | 112 | 0.224 | 13 | 98449.1 | 98429.5 | 2.0e-4 | 98429.2 | ok, 56 vertices | 1 | 62 | 0.5 | 1868 | 543 | 0 | pass | 69 | extrude | 321 | 0.320 | 0.142 | 49 |
 
 ## Categories
 
@@ -39,8 +40,8 @@ clock on the machine that ran it; every other column is deterministic.
 |---|---|---|
 | thin-wall | walls thin against the part's size | 7 |
 | small-feature | features much smaller than the part | 1 |
-| fillet | blends between faces | 2 |
-| genus | through holes (genus > 0) | 7 |
+| fillet | blends between faces | 3 |
+| genus | through holes (genus > 0) | 8 |
 | sharp | sharp edges and acute dihedral angles | 3 |
 | near-tangent | faces meeting at a shallow angle | 2 |
 | boolean | intersection curves from booleans | 3 |
@@ -48,6 +49,7 @@ clock on the machine that ran it; every other column is deterministic.
 | blocky | axis-aligned blocks (polycube-friendly) | 4 |
 | organic | free-form, lofted or revolved spline surfaces | 2 |
 | assembly | several bodies placed by component transforms | 1 |
+| partitioned | split into cells that each mesh as a mapped block (conforming multi-block hex) | 1 |
 | analytic | closed-form geometry for solver verification | 5 |
 
 ## Parts
@@ -57,6 +59,7 @@ clock on the machine that ran it; every other column is deterministic.
 - `cantilever-i-beam`: The 200 mm I-beam of `i-beam` lying along X with its web upright (Z up): a sweep along a horizontal axis, and the cantilever of the README demo.
 - `cross-cylinders`: Two orthogonal R10 cylinders joined: equal radii, so the intersection curves cross at two singular points.
 - `filleted-bracket`: An 8 mm L-bracket with an R6 inside fillet and a hole through each leg: a blend beside holes.
+- `filleted-bracket-blocks`: The filleted bracket split into 30 cells with SolidCraft's Split Body (three slabs across its width, the corner halved along its diagonal, an O-grid of four cells round each hole) so that every cell is a mapped block: conforming multi-block hex meshing.
 - `i-beam`: A 200 mm I-beam with 3 mm flanges and a 2 mm web: a long sweepable prism with thin walls.
 - `loft-transition`: A 40 mm square lofted to an R15 circle 50 mm above: a ruled free-form surface between unlike sections.
 - `perforated-plate`: A 3 mm plate perforated by a 6 x 4 grid of 6 mm holes: genus 24 in a thin sheet.

docs/meshing/literature.md

@@ -880,7 +880,7 @@ input points), Steiner points kept half the local spacing from the boundary.
 
 ### B8. Mapped and swept hexahedral meshing (implemented in `fieldcraft-mesh::hex`)
 
-`mesh.hex.generate` (and `mesh.hex.sweep`, `mesh.hex.mapped`). Every method below is implemented from
+`mesh.hex.generate` (and `mesh.hex.sweep`, `mesh.hex.mapped`, `mesh.hex.blocks`). Every method below is implemented from
 its paper; no code was read. DOIs marked "checked" were confirmed on 2026-10-13 through a repository
 record or the publisher; the others are cited by venue only.
 
@@ -932,7 +932,17 @@ record or the publisher; the others are cited by venue only.
   Applied Mathematics and Computation 10–11:171–233, 1982.
 - **Use.** Interior nodes of a logical cube from its six face grids (`hex::mapped`).
 
-**B8.8 Support.** Generalized winding number for the inside test (A. Jacobson, L. Kavan, O. Sorkine-Hornung,
+**B8.8 Multi-block meshing of partitioned parts.** **Verdict: BUILD ON (conforming cells).**
+- **Citations.** J. F. Thompson, B. K. Soni, N. P. Weatherill (eds.). *Handbook of Grid Generation.* CRC Press,
+  1999, ch. 3 (transfinite and algebraic grids) and ch. 13 (multiblock grids). A. Sheffer, M. Bercovier,
+  T. Blacker, J. Clements. *Virtual topology operators for meshing.* IJCGA 10(3):309–331, 2000.
+- **Use.** `hex::cells` (`mesh.hex.blocks`): cells that meet face to face, each a logical cube after smoothly
+  meeting faces are joined into composite faces (virtual topology); vertices, edges and faces shared by
+  identity; interval assignment over all cells (B8.5); each edge divided once, each face's transfinite grid
+  made once and used by both cells; edges that leave a curved face graded from the spacing along it with a
+  growth of 1.2 per element; interiors by B8.7. The merged mesh is checked for free internal faces.
+
+**B8.9 Support.** Generalized winding number for the inside test (A. Jacobson, L. Kavan, O. Sorkine-Hornung,
 *Robust inside-outside segmentation using generalized winding numbers*, ACM TOG 32(4):33, 2013; solid angle
 of a triangle: Van Oosterom and Strackee 1983); closest point on a triangle (C. Ericson, *Real-Time
 Collision Detection*, 2005, §5.1.5).

tests/corpus/filleted-bracket-blocks/README.md

@@ -0,0 +1 @@
+The filleted bracket split into 30 cells with SolidCraft's Split Body (three slabs across its width, the corner halved along its diagonal, an O-grid of four cells round each hole) so that every cell is a mapped block: conforming multi-block hex meshing.

tests/corpus/filleted-bracket-blocks/recipe.json

@@ -0,0 +1,39 @@
+{"corpus": {"categories": ["fillet", "genus", "partitioned"], "volume_mm3": 32144.778925},
+ "commands": [
+  {"command": "sketch.create", "params": {"plane": "XY", "name": "Profile"}},
+  {"command": "sketch.line", "params": {"points": [[0, 0], [60, 0], [60, 8], [8, 8], [8, 50], [0, 50]], "closed": true}},
+  {"command": "sketch.finish", "params": {}},
+  {"command": "solid.extrude", "params": {"distance": 40, "body_name": "Bracket"}},
+  {"command": "solid.fillet", "params": {"edges": [[8, 8, 20]], "radius": 6}},
+  {"command": "solid.hole", "params": {"position": [40, 8, 20], "diameter": 8}},
+  {"command": "solid.hole", "params": {"position": [8, 32, 20], "diameter": 8}},
+  {"command": "solid.split_body", "params": {"body": "Bracket", "plane": {"origin": [0, 0, 30], "normal": [0, 0, 1]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1)", "plane": {"origin": [0, 0, 10], "normal": [0, 0, 1]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1) (1)", "plane": {"origin": [16, 0, 0], "normal": [1, 0, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1) (1) (1)", "plane": {"origin": [0, 16, 0], "normal": [0, 1, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1) (1) (1) (1)", "plane": {"origin": [0, 0, 0], "normal": [1, -1, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1) (1)", "plane": {"origin": [30, 0, 0], "normal": [1, 0, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1) (1)", "plane": {"origin": [50, 0, 0], "normal": [1, 0, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1) (1) (1)", "plane": {"origin": [0, 22, 0], "normal": [0, 1, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1) (1) (1)", "plane": {"origin": [0, 42, 0], "normal": [0, 1, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1)", "plane": {"origin": [16, 0, 0], "normal": [1, 0, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1) (1)'", "plane": {"origin": [0, 16, 0], "normal": [0, 1, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1) (1)' (1)", "plane": {"origin": [0, 0, 0], "normal": [1, -1, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1)", "plane": {"origin": [30, 0, 0], "normal": [1, 0, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1)", "plane": {"origin": [50, 0, 0], "normal": [1, 0, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1) (1)'", "plane": {"origin": [0, 22, 0], "normal": [0, 1, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1) (1)'", "plane": {"origin": [0, 42, 0], "normal": [0, 1, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket", "plane": {"origin": [16, 0, 0], "normal": [1, 0, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1)'", "plane": {"origin": [0, 16, 0], "normal": [0, 1, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1)' (1)", "plane": {"origin": [0, 0, 0], "normal": [1, -1, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket", "plane": {"origin": [30, 0, 0], "normal": [1, 0, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket", "plane": {"origin": [50, 0, 0], "normal": [1, 0, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1)'", "plane": {"origin": [0, 22, 0], "normal": [0, 1, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1)'", "plane": {"origin": [0, 42, 0], "normal": [0, 1, 0]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1) (1)'''", "plane": {"origin": [40, 0, 20], "normal": [1, 0, 1]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1) (1)'''", "plane": {"origin": [40, 0, 20], "normal": [1, 0, -1]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1) (1)''' (1)", "plane": {"origin": [40, 0, 20], "normal": [1, 0, -1]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1) (1)' (1)''", "plane": {"origin": [0, 32, 20], "normal": [0, 1, 1]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1) (1)' (1)''", "plane": {"origin": [0, 32, 20], "normal": [0, 1, -1]}}},
+  {"command": "solid.split_body", "params": {"body": "Bracket (1) (1)' (1)'' (1)", "plane": {"origin": [0, 32, 20], "normal": [0, 1, -1]}}}
+ ]}