feat(subdivide): four triangles where there was one

The last documented Surface operator. Every edge gets a midpoint, every
triangle becomes four, and the shape does not move.

It does NOT smooth, which the Houdini SOP of this name does. A subdivision that
also moved points would be two operations wearing one name, and the smoothing
one is already Remesh's relaxation — where it can be turned off, tuned, and
told to project back onto the surface it started from. Keeping them separate is
what makes it possible to say which one changed a mesh.

Attributes interpolate onto the midpoints through the remesher's own split, so
a field defined on a coarse mesh survives being refined, and the original
points keep their identities. Depth is capped at six: each level is four times
the triangles, so seven is sixteen thousand times the input and that is a hang
rather than a render.

Co-Authored-By: Claude Opus 5 <[email protected]>

 nodes/subdivide.json | 10 ++++++++++
 shapeshifter.md      |  7 +++++--
 src/geometry.rs      | 30 ++++++++++++++++++++++++++++
 src/main.rs          | 48 ++++++++++++++++++++++++++++++++++++++++++++
 src/remesh.rs        | 56 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 5 files changed, 149 insertions(+), 2 deletions(-)

diff --git a/nodes/subdivide.json b/nodes/subdivide.json
new file mode 100644
index 0000000..dddcb10
--- /dev/null
+++ b/nodes/subdivide.json
@@ -0,0 +1,10 @@
+{
+ "name": "Subdivide",
+ "type": "subdivide",
+ "inputs": 1,
+ "outputs": 1,
+ "params": [
+  { "name": "Input", "type": "text", "default": "" },
+  { "name": "Depth", "type": "spinbox", "default": "1", "min": 0.0, "max": 6.0, "step": 1.0 }
+ ]
+}
diff --git a/shapeshifter.md b/shapeshifter.md
index ee37425..a95490c 100644
--- a/shapeshifter.md
+++ b/shapeshifter.md
@@ -244,8 +244,11 @@ Touches: `geometry.rs`, `app.rs`, `render.rs`, `playbar.rs`.
 > parameter names carried a reading, these carry nothing. Say what they do and
 > they are a short job each.
 >
-> Remaining: `subdivide`, which is the one documented Surface operator still
-> missing and is the easy one.
+> `subdivide` is in: four triangles where there was one, attributes
+> interpolated onto the midpoints, and the shape left exactly where it was —
+> it refines, it does not smooth, which is what separates it from Remesh.
+>
+> Phase 3 is done but for the two operators nobody can describe.
 
 `Develop` is easy — displace along the normal by a development attribute.
 **Remesh is the hard one**, and it is load-bearing: without topology that keeps
diff --git a/src/geometry.rs b/src/geometry.rs
index ffe2195..410417b 100644
--- a/src/geometry.rs
+++ b/src/geometry.rs
@@ -740,6 +740,8 @@ pub fn generate_single_node_geometry_with_errors(
         resolve_neighbour_geometry_with_errors(root, target, visited, ocl_error, sim)
     } else if target.node_type.eq_ignore_ascii_case("time") {
         resolve_time_geometry_with_errors(root, target, visited, ocl_error, sim)
+    } else if target.node_type.eq_ignore_ascii_case("subdivide") {
+        resolve_subdivide_geometry_with_errors(root, target, visited, ocl_error, sim)
     } else if target.node_type.eq_ignore_ascii_case("detangle") {
         resolve_detangle_geometry_with_errors(root, target, visited, ocl_error, sim)
     } else if target.node_type.eq_ignore_ascii_case("suture") {
@@ -1296,6 +1298,24 @@ pub fn resolve_relax_geometry_with_errors(
     Some(geom)
 }
 
+/// The Subdivide node: four triangles where there was one.
+pub fn resolve_subdivide_geometry_with_errors(
+    root: &FsNode,
+    target: &FsNode,
+    visited: &mut Vec<String>,
+    ocl_error: &mut Option<String>,
+    sim: &mut EvalSim,
+) -> Option<Detail> {
+    let input_name = node_param_str(target, "Input", "");
+    if input_name.is_empty() {
+        return None;
+    }
+    let input_node = find_node_by_name(root, &input_name)?;
+    let geom = generate_single_node_geometry_with_errors(root, input_node, visited, ocl_error, sim)?;
+    let depth = node_param_f32(target, "Depth", 1.0).clamp(0.0, 6.0) as usize;
+    Some(crate::remesh::subdivide(&geom, depth))
+}
+
 /// The Detangle node: push a surface off itself.
 ///
 /// Growth folds a surface into its own neighbourhood long before it looks
@@ -3939,6 +3959,7 @@ pub fn is_geometry_node_type(node_type: &str) -> bool {
         || nt == "remesh"
         || nt == "suture"
         || nt == "detangle"
+        || nt == "subdivide"
         || nt == "attribute"
         || nt == "simnet"
 }
@@ -4086,6 +4107,15 @@ pub fn network_sphere_vertices_with_errors(
                     out.merge(&geom);
                 }
             }
+        } else if node.node_type.eq_ignore_ascii_case("subdivide") {
+            let _idx = *count;
+            *count += 1;
+            if is_visible {
+                let mut visited = Vec::new();
+                if let Some(geom) = resolve_subdivide_geometry_with_errors(root, node, &mut visited, ocl_error, sim) {
+                    out.merge(&geom);
+                }
+            }
         } else if node.node_type.eq_ignore_ascii_case("detangle") {
             let _idx = *count;
             *count += 1;
diff --git a/src/main.rs b/src/main.rs
index bd943fb..168cb3b 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -3623,6 +3623,54 @@ mod tests {
         }
     }
 
+    #[test]
+    fn test_subdivide_multiplies_triangles_without_moving_the_shape() {
+        use crate::remesh::subdivide;
+        let mut sphere = sphere_detail(Vec3::ZERO, 1.0, 6, 8);
+        sphere.points_mut().create("mass", AttribValue::Float(0.0));
+        for p in 0..sphere.num_points() {
+            let y = sphere.pos(p).y;
+            sphere.points_mut().set_value("mass", p, AttribValue::Float(y)).unwrap();
+        }
+        let (before_prims, before_bounds) = (sphere.num_prims(), sphere.bounds().unwrap());
+
+        let once = subdivide(&sphere, 1);
+        // Every triangle becomes four. The sphere's quad bands fan to two
+        // triangles each on the way in, so the count is against THAT.
+        let tri_count = |d: &Detail| d.triangulate(|_, _| ()).len() / 3;
+        assert_eq!(tri_count(&once), tri_count(&sphere) * 4);
+
+        // The shape does not move: this refines, it does not smooth. A
+        // subdivision that also moved points would be two operations wearing
+        // one name.
+        let after_bounds = once.bounds().unwrap();
+        assert!((after_bounds.0 - before_bounds.0).length() < 1e-5, "{:?}", after_bounds.0);
+        assert!((after_bounds.1 - before_bounds.1).length() < 1e-5, "{:?}", after_bounds.1);
+        // Original points keep their exact positions AND identities.
+        for p in 0..sphere.num_points() {
+            assert!(once.ids().contains(&sphere.id(p).unwrap()), "point {p} lost its identity");
+        }
+
+        // Attributes interpolate onto the midpoints, so a field defined on a
+        // coarse mesh survives being refined.
+        for p in 0..once.num_points() {
+            let v = once.points().value("mass", p).unwrap().as_f32();
+            assert!((v - once.pos(p).y).abs() < 0.12, "point {p}: {v} vs y {}", once.pos(p).y);
+        }
+
+        // Depth compounds, and zero is a pass-through.
+        assert_eq!(tri_count(&subdivide(&sphere, 2)), tri_count(&sphere) * 16);
+        assert_eq!(subdivide(&sphere, 0).num_prims(), before_prims);
+
+        // The result is still a closed, usable mesh.
+        for prim in 0..once.num_prims() {
+            assert_eq!(once.prim_points(prim).len(), 3);
+        }
+        for p in 0..once.num_points() {
+            assert!(!once.point_prims(p).is_empty(), "point {p} belongs to nothing");
+        }
+    }
+
     #[test]
     fn test_the_projection_pass_stops_a_remeshed_surface_creeping() {
         use crate::spatial::TriGrid;
diff --git a/src/remesh.rs b/src/remesh.rs
index a622ecc..e4bc318 100644
--- a/src/remesh.rs
+++ b/src/remesh.rs
@@ -589,6 +589,62 @@ fn project_pass(m: &mut Mesh, rest: &crate::spatial::TriGrid) {
     }
 }
 
+/// Subdivide every triangle into four, `depth` times.
+///
+/// Distinct from remeshing, and deliberately so: this makes a predictable,
+/// uniform refinement of the mesh it is given — every edge gets a midpoint,
+/// every triangle becomes four, and the shape does not move. Remesh steers
+/// toward a length and rearranges topology to get there; Subdivide multiplies
+/// what is already there.
+///
+/// It does NOT smooth, which the Houdini SOP of this name does. A subdivision
+/// that also moved points would be two operations wearing one name, and the
+/// smoothing one is already available as Remesh's relaxation.
+///
+/// Attributes interpolate onto the midpoints, the same way a remesh split
+/// does, so a field defined on a coarse mesh survives being refined.
+pub fn subdivide(input: &Detail, depth: usize) -> Detail {
+    if input.num_prims() == 0 || depth == 0 {
+        return input.clone();
+    }
+    let mut m = Mesh::from_detail(input);
+    // Capped because this is exponential: each level is four times the
+    // triangles, so six levels is four thousand times the input and anything
+    // past that is a hang rather than a render.
+    for _ in 0..depth.min(6) {
+        let mut mids: HashMap<[u32; 2], u32> = HashMap::new();
+        for (e, _) in m.edges() {
+            let mid = m.split_point(e[0], e[1]);
+            mids.insert(e, mid);
+        }
+        let key = |a: u32, b: u32| [a.min(b), a.max(b)];
+        // Snapshotted, because the loop adds triangles as it goes and the new
+        // ones are already subdivided.
+        let count = m.tris.len();
+        for t in 0..count {
+            if m.dead_tri[t] {
+                continue;
+            }
+            let tri = m.tris[t];
+            let (Some(&ab), Some(&bc), Some(&ca)) = (
+                mids.get(&key(tri[0], tri[1])),
+                mids.get(&key(tri[1], tri[2])),
+                mids.get(&key(tri[2], tri[0])),
+            ) else {
+                continue;
+            };
+            m.dead_tri[t] = true;
+            // Three corner triangles and the middle one, each keeping the
+            // original winding.
+            m.add_tri([tri[0], ab, ca]);
+            m.add_tri([ab, tri[1], bc]);
+            m.add_tri([ca, bc, tri[2]]);
+            m.add_tri([ab, bc, ca]);
+        }
+    }
+    m.into_detail()
+}
+
 /// Remesh toward `settings.target` edge length.
 pub fn remesh(input: &Detail, settings: Settings) -> Detail {
     if input.num_prims() == 0 || settings.target <= 0.0 {