From 87710a83ab4747f6ceac595f4513122fecca6c33 Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Thu, 15 Oct 2026 22:10:24 -0400 Subject: [PATCH] io: write element faces as binary and ASCII STL Model::skin_faces gives the mesh's outer skin: the exterior faces of the solids and the top side of every shell. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo --- crates/io/src/lib.rs | 1 + crates/io/src/stl.rs | 225 ++++++++++++++++++++++++++++++++++++ crates/model/src/contact.rs | 42 +++++++ 3 files changed, 268 insertions(+) create mode 100644 crates/io/src/stl.rs diff --git a/crates/io/src/lib.rs b/crates/io/src/lib.rs index 753b755..63d8c88 100644 --- a/crates/io/src/lib.rs +++ b/crates/io/src/lib.rs @@ -27,6 +27,7 @@ pub mod memory; pub mod mesh; pub mod recovery; pub mod roots; +pub mod stl; use std::collections::{BTreeMap, BTreeSet}; use std::io::{Cursor, Read, Write}; diff --git a/crates/io/src/stl.rs b/crates/io/src/stl.rs new file mode 100644 index 0000000..ab1eb36 --- /dev/null +++ b/crates/io/src/stl.rs @@ -0,0 +1,225 @@ +//! STL: element faces written as triangles, for mesh viewers, slicers and tools that read only +//! triangles. +//! +//! The layout follows the StereoLithography Interface Specification (3D Systems, 1988/1989). +//! Binary: an 80-byte header that does not start with `solid` (readers take that as ASCII), the +//! triangle count as a little-endian `u32`, then 50 bytes per triangle: the normal and the three +//! vertices as little-endian `f32` triples, and a `u16` attribute of 0. ASCII: `solid name`, +//! one `facet normal … outer loop … endloop endfacet` per triangle, `endsolid name`. +//! +//! Faces are triangulated through their corner nodes (a quadrilateral as 0-1-2 and 0-2-3), so the +//! vertices run counter-clockwise seen from outside when the faces do ([`ElementKind::faces`]). +//! Each facet's normal comes from its own vertices. Coordinates are rounded to single precision. +//! +//! [`ElementKind::faces`]: fieldcraft_model::ElementKind::faces + +use fieldcraft_model::{ElementFace, Model}; + +use crate::{Error, Result}; + +/// The most triangles one file holds: two per face of a model at its element cap is far below. +pub const MAX_TRIANGLES: usize = 200_000_000; + +const HEADER: &[u8] = b"FieldCraft STL"; + +type Triangle = [[f64; 3]; 3]; + +/// The triangles of `faces`, in their order (each quadrilateral as two). +pub fn triangles(model: &Model, faces: &[ElementFace]) -> Result> { + let mut out = Vec::new(); + for &face in faces { + let nodes = model.face_nodes(face).ok_or_else(|| Error::Invalid(format!("element {} has no face {}", face.element, face.face)))?; + let corners = if nodes.len() >= 6 { nodes.len() / 2 } else { nodes.len() }; + let mut p = Vec::with_capacity(corners); + for label in nodes.iter().take(corners) { + let n = model.node(*label).ok_or_else(|| Error::Invalid(format!("element {} uses missing node {label}", face.element)))?; + p.push(n.xyz); + } + match p.as_slice() { + [a, b, c] => out.push([*a, *b, *c]), + [a, b, c, d] => out.extend([[*a, *b, *c], [*a, *c, *d]]), + _ => return Err(Error::Invalid(format!("face {} of element {} is not a triangle or a quadrilateral", face.face, face.element))), + } + if out.len() > MAX_TRIANGLES { + return Err(Error::Limit(format!("an STL file here holds at most {MAX_TRIANGLES} triangles"))); + } + } + Ok(out) +} + +/// The unit normal of a triangle from its vertices (right-hand rule); zero when degenerate. +fn normal(t: &Triangle) -> [f64; 3] { + let [a, b, c] = t; + let u = [b[0] - a[0], b[1] - a[1], b[2] - a[2]]; + let v = [c[0] - a[0], c[1] - a[1], c[2] - a[2]]; + let n = [u[1] * v[2] - u[2] * v[1], u[2] * v[0] - u[0] * v[2], u[0] * v[1] - u[1] * v[0]]; + let len = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt(); + if len > 0.0 && len.is_finite() { [n[0] / len, n[1] / len, n[2] / len] } else { [0.0; 3] } +} + +/// `faces` as binary STL. +pub fn write_binary(model: &Model, faces: &[ElementFace]) -> Result> { + let tris = triangles(model, faces)?; + let count = u32::try_from(tris.len()).map_err(|_| Error::Limit(format!("binary STL holds at most {} triangles", u32::MAX)))?; + let mut out = Vec::with_capacity(84 + 50 * tris.len()); + out.extend_from_slice(HEADER); + out.resize(80, 0); + out.extend_from_slice(&count.to_le_bytes()); + for t in &tris { + for v in std::iter::once(normal(t)).chain(t.iter().copied()) { + for x in v { + out.extend_from_slice(&(x as f32).to_le_bytes()); + } + } + out.extend_from_slice(&0u16.to_le_bytes()); + } + Ok(out) +} + +/// `faces` as ASCII STL, the solid called `name` (spaces and line breaks become `_`). +pub fn write_ascii(model: &Model, faces: &[ElementFace], name: &str) -> Result { + use std::fmt::Write; + let tris = triangles(model, faces)?; + let name: String = name.chars().take(80).map(|c| if c.is_whitespace() || c.is_control() { '_' } else { c }).collect(); + let mut out = String::new(); + let _ = writeln!(out, "solid {name}"); + let xyz = |v: [f64; 3]| format!("{} {} {}", v[0] as f32, v[1] as f32, v[2] as f32); + for t in &tris { + let _ = writeln!(out, " facet normal {}", xyz(normal(t))); + out.push_str(" outer loop\n"); + for v in t { + let _ = writeln!(out, " vertex {}", xyz(*v)); + } + out.push_str(" endloop\n endfacet\n"); + } + let _ = writeln!(out, "endsolid {name}"); + Ok(out) +} + +#[cfg(test)] +mod tests { + use super::*; + use fieldcraft_model::ElementKind; + + /// An `n`×`n`×`n` block of hex8 of side `size`, node labels from 1. + fn block(n: u32, size: f64) -> Model { + let mut m = Model::default(); + let id = |i: u32, j: u32, k: u32| 1 + i + (n + 1) * (j + (n + 1) * k); + let h = size / f64::from(n); + for k in 0..=n { + for j in 0..=n { + for i in 0..=n { + m.set_node(id(i, j, k), [f64::from(i) * h, f64::from(j) * h, f64::from(k) * h]).unwrap(); + } + } + } + let mut label = 1; + for k in 0..n { + for j in 0..n { + for i in 0..n { + let b = [id(i, j, k), id(i + 1, j, k), id(i + 1, j + 1, k), id(i, j + 1, k)]; + let t = [id(i, j, k + 1), id(i + 1, j, k + 1), id(i + 1, j + 1, k + 1), id(i, j + 1, k + 1)]; + m.set_element(label, ElementKind::Hex8, b.into_iter().chain(t).collect()).unwrap(); + label += 1; + } + } + } + m + } + + /// The triangles read back from binary STL: (normal, vertices). + fn read(bytes: &[u8]) -> Vec<([f32; 3], [[f32; 3]; 3])> { + let count = u32::from_le_bytes(bytes[80..84].try_into().unwrap()) as usize; + assert_eq!(bytes.len(), 84 + 50 * count); + let f = |o: usize| f32::from_le_bytes(bytes[o..o + 4].try_into().unwrap()); + let v = |o: usize| [f(o), f(o + 4), f(o + 8)]; + (0..count) + .map(|t| { + let o = 84 + 50 * t; + assert_eq!(&bytes[o + 48..o + 50], &[0, 0]); + (v(o), [v(o + 12), v(o + 24), v(o + 36)]) + }) + .collect() + } + + /// The volume a closed triangulated surface encloses (divergence theorem); positive when + /// every normal points out. + fn volume(tris: &[([f32; 3], [[f32; 3]; 3])]) -> f64 { + tris.iter() + .map(|(_, [a, b, c])| { + let [a, b, c] = [a, b, c].map(|p| p.map(f64::from)); + (a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + a[2] * (b[0] * c[1] - b[1] * c[0])) / 6.0 + }) + .sum() + } + + fn check_closed(m: &Model, expected_tris: usize, expected_volume: f64, centre: [f64; 3]) { + let bytes = write_binary(m, &m.skin_faces()).unwrap(); + assert!(!bytes.starts_with(b"solid")); + let tris = read(&bytes); + assert_eq!(tris.len(), expected_tris); + assert!((volume(&tris) - expected_volume).abs() < 1e-5, "{}", volume(&tris)); + for (n, t) in &tris { + let fc: Vec = (0..3).map(|i| t.iter().map(|p| f64::from(p[i])).sum::() / 3.0 - centre[i]).collect(); + let dot: f64 = (0..3).map(|i| f64::from(n[i]) * fc[i]).sum(); + assert!(dot > 0.0, "normal {n:?} of {t:?} points in"); + } + assert_eq!(bytes, write_binary(m, &m.skin_faces()).unwrap(), "deterministic"); + } + + #[test] + fn a_hex_block_gives_two_triangles_per_boundary_quad_and_its_volume() { + // 2×2×2 hex8 of side 2: 24 boundary quads. + check_closed(&block(2, 2.0), 48, 8.0, [1.0; 3]); + } + + #[test] + fn tet10_and_wedge6_close_with_outward_normals() { + let mut t = Model::default(); + let p = [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.], [0.5, 0., 0.], [0.5, 0.5, 0.], [0., 0.5, 0.], [0., 0., 0.5], [0.5, 0., 0.5], [0., 0.5, 0.5]]; + for (i, xyz) in p.iter().enumerate() { + t.set_node(i as u32 + 1, *xyz).unwrap(); + } + t.set_element(1, ElementKind::Tet10, (1..=10).collect()).unwrap(); + check_closed(&t, 4, 1.0 / 6.0, [0.25; 3]); + + let mut w = Model::default(); + for (i, xyz) in [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [0., 1., 1.]].iter().enumerate() { + w.set_node(i as u32 + 1, *xyz).unwrap(); + } + w.set_element(1, ElementKind::Wedge6, (1..=6).collect()).unwrap(); + // Two triangles and three quads: 8 facets, volume 1/2. + check_closed(&w, 8, 0.5, [1.0 / 3.0, 1.0 / 3.0, 0.5]); + } + + #[test] + fn a_shell_plate_writes_one_pair_per_quad_along_its_normal() { + let mut m = Model::default(); + for (i, xy) in [[0., 0.], [1., 0.], [2., 0.], [0., 1.], [1., 1.], [2., 1.]].iter().enumerate() { + m.set_node(i as u32 + 1, [xy[0], xy[1], 0.0]).unwrap(); + } + m.set_element(1, ElementKind::Quad4, vec![1, 2, 5, 4]).unwrap(); + m.set_element(2, ElementKind::Quad4, vec![2, 3, 6, 5]).unwrap(); + let tris = read(&write_binary(&m, &m.skin_faces()).unwrap()); + assert_eq!(tris.len(), 4); + assert!(tris.iter().all(|(n, _)| *n == [0.0, 0.0, 1.0])); + } + + #[test] + fn ascii_holds_the_same_facets() { + let m = block(1, 1.0); + let text = write_ascii(&m, &m.skin_faces(), "my part").unwrap(); + assert!(text.starts_with("solid my_part\n") && text.ends_with("endsolid my_part\n"), "{text}"); + assert_eq!(text.matches("facet normal").count(), 12); + assert_eq!(text.matches("vertex").count(), 36); + assert!(text.contains("facet normal 0 0 -1"), "{text}"); + } + + #[test] + fn faces_that_are_not_there_are_refused() { + let m = block(1, 1.0); + assert!(write_binary(&m, &[ElementFace { element: 9, face: 1 }]).is_err()); + assert!(write_binary(&m, &[ElementFace { element: 1, face: 7 }]).is_err()); + assert_eq!(write_binary(&m, &[]).unwrap().len(), 84); + } +} diff --git a/crates/model/src/contact.rs b/crates/model/src/contact.rs index f7a6a3d..eab4586 100644 --- a/crates/model/src/contact.rs +++ b/crates/model/src/contact.rs @@ -232,6 +232,16 @@ impl Model { self.unshared(3, |kind| kind.faces(), 4, |element, face| crate::ElementFace { element, face }) } + /// The mesh's outer skin, in element order: the [`exterior_faces`](Model::exterior_faces) + /// of the solids and face 1 (the side its normal points to) of every 2D element. Each face's + /// corners run counter-clockwise seen from outside ([`crate::ElementKind::faces`]). + pub fn skin_faces(&self) -> Vec { + let shells = self.elements.iter().filter(|(_, e)| e.kind.dimension() == 2).map(|(&element, _)| crate::ElementFace { element, face: 1 }); + let mut out: Vec = self.exterior_faces().into_iter().chain(shells).collect(); + out.sort_unstable(); + out + } + /// The sides (`sides` of each element kind) of the elements of `dimension` that no other /// such element shares, each known by its sorted corners (at most `corners`; triangles of /// quadratic kinds by their first three), made by `make` from the element and the 1-based @@ -376,6 +386,38 @@ mod tests { assert!(Model::default().exterior_faces().is_empty()); } + #[test] + fn skin_faces_count_solid_and_shell_sides() { + let cube = |m: &mut Model, base: u32, z: f64| { + for (i, xy) in [[0., 0.], [1., 0.], [1., 1.], [0., 1.]].iter().enumerate() { + m.set_node(base + i as u32, [xy[0], xy[1], z]).unwrap(); + } + }; + let mut m = Model::default(); + cube(&mut m, 1, 0.0); + cube(&mut m, 5, 1.0); + m.set_element(1, crate::ElementKind::Hex8, (1..=8).collect()).unwrap(); + assert_eq!(m.skin_faces().len(), 6, "one hex8"); + cube(&mut m, 9, 2.0); + m.set_element(2, crate::ElementKind::Hex8, (5..=12).collect()).unwrap(); + assert_eq!(m.skin_faces().len(), 10, "two hex8 sharing a face"); + + let mut t = Model::default(); + for (i, p) in [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.], [0.5, 0., 0.], [0.5, 0.5, 0.], [0., 0.5, 0.], [0., 0., 0.5], [0.5, 0., 0.5], [0., 0.5, 0.5]] + .iter() + .enumerate() + { + t.set_node(i as u32 + 1, *p).unwrap(); + } + t.set_element(1, crate::ElementKind::Tet10, (1..=10).collect()).unwrap(); + assert_eq!(t.skin_faces().len(), 4, "a tet10"); + + let mut s = Model::default(); + cube(&mut s, 1, 0.0); + s.set_element(7, crate::ElementKind::Quad4, vec![1, 2, 3, 4]).unwrap(); + assert_eq!(s.skin_faces(), [crate::ElementFace { element: 7, face: 1 }], "a quad4 shell, its top side"); + } + #[test] fn exterior_edges_are_those_of_one_2d_element() { // Two quads sharing an edge and a triangle on the second's top edge.