From ed611ed7ba148f512ea10e33f6255d1039a05b51 Mon Sep 17 00:00:00 2001 From: bherbruck Date: Wed, 7 Oct 2026 20:26:35 +0000 Subject: [PATCH] STEP: reader Own ISO 10303-21 parser and AP203/AP214/AP242 geometry reader in the kernel (step_import): planes, cylinders, cones, spheres, tori, B-spline and rational surfaces, surfaces of revolution/extrusion; lines, conics, B-spline curves; breps with voids, faceted breps, shell models; units (SI prefixes, conversion based), product names, colours, assemblies (NAUO + item-defined transformations, mapped items). Hostile input is capped and guarded; faces that cannot be read are left out and reported. Faces that wrap around a closed non-periodic direction (pipe tori, revolved B-spline bands) are split for truck's mesher. Bodies are meshed per shell so open imports never panic. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01AS5Tn3xieewUSg3uxrMVdL --- crates/kernel/src/body.rs | 27 +- crates/kernel/src/lib.rs | 4 +- crates/kernel/src/step_in/brep.rs | 366 +++++++++++++ crates/kernel/src/step_in/geom.rs | 835 +++++++++++++++++++++++++++++ crates/kernel/src/step_in/mod.rs | 638 ++++++++++++++++++++++ crates/kernel/src/step_in/p21.rs | 551 +++++++++++++++++++ crates/kernel/src/step_in/split.rs | 409 ++++++++++++++ crates/kernel/src/step_in/tests.rs | 285 ++++++++++ 8 files changed, 3109 insertions(+), 6 deletions(-) create mode 100644 crates/kernel/src/step_in/brep.rs create mode 100644 crates/kernel/src/step_in/geom.rs create mode 100644 crates/kernel/src/step_in/mod.rs create mode 100644 crates/kernel/src/step_in/p21.rs create mode 100644 crates/kernel/src/step_in/split.rs create mode 100644 crates/kernel/src/step_in/tests.rs diff --git a/crates/kernel/src/body.rs b/crates/kernel/src/body.rs index d3ea3094..9cd6e5ea 100644 --- a/crates/kernel/src/body.rs +++ b/crates/kernel/src/body.rs @@ -34,6 +34,23 @@ pub struct FaceInfo { pub plane_normal: Option, } +type MeshedShell = truck_topology::Shell, Option>; + +/// Triangulate each shell (not the solid: truck's solid constructor panics on the open or +/// non-manifold shells an import can produce). Faces whose boundary is not exactly on the +/// surface (imported geometry within the file's tolerance) fail the exact parameter search; +/// such shells are meshed again projecting onto the surface. +fn mesh_shells(solid: &Solid, tol: f64) -> Vec { + solid + .boundaries() + .iter() + .map(|sh| { + let m = sh.triangulation(tol); + if m.face_iter().any(|f| f.surface().is_none()) { sh.robust_triangulation(tol) } else { m } + }) + .collect() +} + pub(crate) fn p3(v: Vec3) -> mt::Point3 { mt::Point3::new(v.x, v.y, v.z) } @@ -65,7 +82,7 @@ impl Body { return Err(KernelError::Failed("empty result".into())); } // Orient outward: a closed solid must have positive volume. - let v = guard("orient", || Ok(solid.triangulation(1.0).to_polygon().volume()))?; + let v = guard("orient", || Ok(mesh_shells(&solid, 1.0).iter().map(|s| s.to_polygon().volume()).sum::()))?; if v < 0.0 { solid.not(); } @@ -98,9 +115,9 @@ impl Body { let tol = if tol.is_finite() && tol > 0.0 { tol } else { 0.05 }; let tol = tol.max(self.size() * 1e-6); guard("tessellate", || { - let meshed = self.solid.triangulation(tol); + let meshed = mesh_shells(&self.solid, tol); let mut out = Mesh::default(); - for (fi, face) in meshed.face_iter().enumerate() { + for (fi, face) in meshed.iter().flat_map(|s| s.face_iter()).enumerate() { let Some(pm) = face.surface() else { continue }; let flip = !face.orientation(); let base = u32::try_from(out.positions.len()).unwrap_or(u32::MAX); @@ -142,7 +159,7 @@ impl Body { } } let mut seen = std::collections::HashMap::new(); - for e in meshed.edge_iter() { + for e in meshed.iter().flat_map(|s| s.edge_iter()) { if seen.contains_key(&e.id()) { continue; } @@ -152,7 +169,7 @@ impl Body { } // Which faces each edge bounds. out.edge_faces = vec![Vec::new(); out.edges.len()]; - for (fi, face) in meshed.face_iter().enumerate() { + for (fi, face) in meshed.iter().flat_map(|s| s.face_iter()).enumerate() { for e in face.edge_iter() { if let Some(slot) = seen.get(&e.id()).and_then(|i| out.edge_faces.get_mut(*i)) { let f = u32::try_from(fi).unwrap_or(u32::MAX); diff --git a/crates/kernel/src/lib.rs b/crates/kernel/src/lib.rs index ba76df67..0ccf300a 100644 --- a/crates/kernel/src/lib.rs +++ b/crates/kernel/src/lib.rs @@ -2,7 +2,7 @@ //! //! Everything SolveCraft does with solids goes through this crate: building bodies from planar //! regions (extrude, revolve) and primitives, booleans, fillets and chamfers, tessellation, -//! measurement and STEP export. The implementation uses the `truck` crates (Apache-2.0; see +//! measurement and STEP export and import. The implementation uses the `truck` crates (Apache-2.0; see //! `plan/adr/0001-geometry-kernel.md`) plus our own local operations (edge fillets and //! chamfers in `blend.rs`), but no truck type appears in the public API, so the kernel can be //! replaced without touching the rest of the workspace. @@ -20,6 +20,7 @@ mod ops; mod polybool; mod polyhedron; mod step; +mod step_in; mod topo; pub use blend::{chamfer, fillet}; @@ -30,6 +31,7 @@ pub use ops::{BoolOp, boolean, split_by_plane, transform, transform_matrix}; pub use polybool::planar_boolean; pub use polyhedron::{HalfSpace, convex_polyhedron, draft, shell}; pub use step::step_export; +pub use step_in::{ImportedBody, StepImport, StepNode, step_import}; pub use topo::{TopoCounts, merged_topology, seam_flags}; #[derive(Debug, Clone, PartialEq, thiserror::Error)] diff --git a/crates/kernel/src/step_in/brep.rs b/crates/kernel/src/step_in/brep.rs new file mode 100644 index 00000000..dae0cd79 --- /dev/null +++ b/crates/kernel/src/step_in/brep.rs @@ -0,0 +1,366 @@ +//! STEP topology (ISO 10303-42 topology entities) to kernel solids: breps, shells, faces, loops, +//! edges and vertices. A face that cannot be built is left out and reported; the rest of the +//! body is still imported (as an open body when a face is missing). + +use std::collections::HashMap; + +use mt::{BoundedCurve, EuclideanSpace, InnerSpace, MetricSpace, ParametricCurve}; +use truck_modeling as mt; + +use super::Ctx; +use super::geom::{self, R, entity}; +use super::p21::Param; +use crate::guard; + +/// Most faces in one body. +const MAX_FACES: usize = 500_000; +/// Most edges in one loop. +const MAX_LOOP: usize = 100_000; + +/// A body read from a brep entity. +pub(crate) struct BrepOut { + pub solid: mt::Solid, + /// Faces left out (with the reason). + pub failed: Vec, + /// The shells are closed (a solid) rather than open (missing faces, or a surface model). + pub closed: bool, + /// Faces read from the file (before band faces are split for the kernel). + pub file_faces: usize, +} + +#[derive(Default)] +struct Builder { + vertices: HashMap, + edges: HashMap>, + poly_edges: HashMap<(u64, u64), mt::Edge>, + faces: usize, +} + +fn refs(p: Option<&Param>) -> Vec { + p.and_then(Param::as_list).map(|v| v.iter().filter_map(Param::as_ref_id).collect()).unwrap_or_default() +} + +fn rid(p: Option<&Param>, what: &str) -> R { + p.and_then(Param::as_ref_id).ok_or_else(|| format!("{what}: expected a reference")) +} + +impl Builder { + fn vertex(&mut self, cx: &Ctx, id: u64) -> R { + if let Some(v) = self.vertices.get(&id) { + return Ok(v.clone()); + } + let e = entity(cx, id)?; + if e.name() != "VERTEX_POINT" { + return Err(format!("#{id}: expected VERTEX_POINT, found {}", e.name())); + } + let p = geom::point(cx, rid(e.params().get(1), "VERTEX_POINT")?)?; + let v = mt::Vertex::new(p); + self.vertices.insert(id, v.clone()); + Ok(v) + } + + /// Kernel edges (start to end) of an EDGE_CURVE. + fn edge_curve(&mut self, cx: &Ctx, id: u64) -> R> { + if let Some(e) = self.edges.get(&id) { + return Ok(e.clone()); + } + let e = entity(cx, id)?; + if e.name() != "EDGE_CURVE" { + return Err(format!("#{id}: unsupported edge {}", e.name())); + } + let p = e.params(); + let (vs_id, ve_id) = (rid(p.get(1), "edge start")?, rid(p.get(2), "edge end")?); + let vs = self.vertex(cx, vs_id)?; + let ve = self.vertex(cx, ve_id)?; + let sense = p.get(4).and_then(Param::as_bool).unwrap_or(true); + let g = geom::curve(cx, rid(p.get(3), "edge geometry")?, 0)?; + let closed = vs_id == ve_id; + let pieces = guard("step edge", || geom::edge_pieces(&g, vs.point(), ve.point(), sense, closed, cx.tol).map_err(crate::KernelError::Failed)) + .map_err(|e| e.to_string())?; + let n = pieces.len(); + let mut out = Vec::with_capacity(n); + let mut cur = vs.clone(); + for (i, c) in pieces.into_iter().enumerate() { + let next = if i + 1 == n { ve.clone() } else { mt::Vertex::new(c.back()) }; + out.push(mt::Edge::new_unchecked(&cur, &next, c)); + cur = next; + } + self.edges.insert(id, out.clone()); + Ok(out) + } + + /// Edges of an ORIENTED_EDGE (or a bare EDGE_CURVE) in loop order. + fn oriented_edge(&mut self, cx: &Ctx, id: u64, depth: usize) -> R> { + if depth > 4 { + return Err("oriented edges nested too deeply".into()); + } + let e = entity(cx, id)?; + match e.name() { + "ORIENTED_EDGE" => { + let p = e.params(); + let inner = self.oriented_edge(cx, rid(p.get(3), "ORIENTED_EDGE")?, depth + 1)?; + if p.get(4).and_then(Param::as_bool).unwrap_or(true) { Ok(inner) } else { Ok(inner.iter().rev().map(mt::Edge::inverse).collect()) } + } + _ => self.edge_curve(cx, id), + } + } + + fn poly_loop(&mut self, cx: &Ctx, pts: &[u64]) -> R> { + let n = pts.len(); + if !(3..=MAX_LOOP).contains(&n) { + return Err(format!("polygon loop with {n} points")); + } + let mut verts = Vec::with_capacity(n); + for &id in pts { + if let Some(v) = self.vertices.get(&id) { + verts.push(v.clone()); + } else { + let v = mt::Vertex::new(geom::point(cx, id)?); + self.vertices.insert(id, v.clone()); + verts.push(v); + } + } + let mut out = Vec::with_capacity(n); + for i in 0..n { + let (Some(&a), Some(&b)) = (pts.get(i), pts.get((i + 1) % n)) else { continue }; + if let Some(e) = self.poly_edges.get(&(b, a)) { + out.push(e.inverse()); + continue; + } + let (Some(va), Some(vb)) = (verts.get(i), verts.get((i + 1) % n)) else { continue }; + if va.point().distance(vb.point()) < 1e-12 { + return Err("polygon loop with repeated points".into()); + } + let e = mt::Edge::try_new(va, vb, mt::Curve::Line(mt::Line(va.point(), vb.point()))).map_err(|e| format!("polygon edge: {e}"))?; + self.poly_edges.insert((a, b), e.clone()); + out.push(e); + } + Ok(out) + } + + /// A face bound as a wire (oriented by the bound's orientation), or `None` for a vertex loop. + fn bound(&mut self, cx: &Ctx, id: u64) -> R> { + let b = entity(cx, id)?; + if !matches!(b.name(), "FACE_BOUND" | "FACE_OUTER_BOUND") { + return Err(format!("#{id}: unsupported face bound {}", b.name())); + } + let lp_id = rid(b.params().get(1), "face bound loop")?; + let orient = b.params().get(2).and_then(Param::as_bool).unwrap_or(true); + let lp = entity(cx, lp_id)?; + let edges = match lp.name() { + "EDGE_LOOP" => { + let ids = refs(lp.params().get(1)); + if ids.is_empty() || ids.len() > MAX_LOOP { + return Err(format!("edge loop with {} edges", ids.len())); + } + let mut v = Vec::new(); + for oe in ids { + v.extend(self.oriented_edge(cx, oe, 0)?); + } + v + } + "POLY_LOOP" => self.poly_loop(cx, &refs(lp.params().get(1)))?, + "VERTEX_LOOP" => return Ok(None), + n => return Err(format!("unsupported loop {n}")), + }; + let mut w: mt::Wire = edges.into(); + if !w.is_closed() { + return Err("face boundary is not a closed loop".into()); + } + if !orient { + w.invert(); + } + Ok(Some(w)) + } + + /// A face (ADVANCED_FACE, FACE_SURFACE, ORIENTED_FACE or a polygonal FACE). + fn face(&mut self, cx: &Ctx, id: u64, depth: usize) -> R { + if depth > 4 { + return Err("oriented faces nested too deeply".into()); + } + self.faces += 1; + if self.faces > MAX_FACES { + return Err(format!("more than {MAX_FACES} faces")); + } + let e = entity(cx, id)?; + let p = e.params(); + if e.name() == "ORIENTED_FACE" { + let mut f = self.face(cx, rid(p.get(2), "ORIENTED_FACE")?, depth + 1)?; + if !p.get(3).and_then(Param::as_bool).unwrap_or(true) { + f.invert(); + } + return Ok(f); + } + if !matches!(e.name(), "ADVANCED_FACE" | "FACE_SURFACE" | "FACE") { + return Err(format!("unsupported face {}", e.name())); + } + let mut wires = Vec::new(); + for b in refs(p.get(1)) { + if let Some(w) = self.bound(cx, b)? { + wires.push(w); + } + } + if wires.is_empty() { + return Err("face without boundary edges".into()); + } + let extent: Vec = wires + .iter() + .flat_map(|w| w.edge_iter()) + .flat_map(|e| { + let c = e.curve(); + let (t0, t1) = c.range_tuple(); + (0..8).map(move |i| c.subs(t0 + (t1 - t0) * i as f64 / 8.0)).collect::>() + }) + .collect(); + // Seam edges: used twice by the face's loops. + let mut uses: HashMap = HashMap::new(); + for e in wires.iter().flat_map(|w| w.edge_iter()) { + *uses.entry(e.id()).or_default() += 1; + } + let mut seen = std::collections::HashSet::new(); + let seam: Vec = wires + .iter() + .flat_map(|w| w.edge_iter()) + .filter(|e| uses.get(&e.id()).copied().unwrap_or(0) > 1 && seen.insert(e.id())) + .flat_map(|e| { + let c = e.curve(); + let (t0, t1) = c.range_tuple(); + (0..=8).map(move |i| c.subs(t0 + (t1 - t0) * i as f64 / 8.0)).collect::>() + }) + .collect(); + let surface = if e.name() == "FACE" { + plane_through(&extent).ok_or("polygonal face without a plane")? + } else { + let sid = rid(p.get(2), "face geometry")?; + let sense = p.get(3).and_then(Param::as_bool).unwrap_or(true); + guard("step surface", || geom::surface(cx, sid, &extent, &seam, sense, 0).map_err(crate::KernelError::Failed)) + .map_err(|e| e.to_string())? + }; + match mt::Face::try_new(wires.clone(), surface.clone()) { + Ok(f) => Ok(f), + // Faces on periodic surfaces run along a seam edge twice: valid in STEP. + Err(truck_topology::errors::Error::NotSimpleWire) => Ok(mt::Face::new_unchecked(wires, surface)), + Err(e) => Err(format!("face: {e}")), + } + } + + fn shell(&mut self, cx: &Ctx, id: u64, failed: &mut Vec) -> R<(mt::Shell, bool, usize)> { + let e = entity(cx, id)?; + let p = e.params(); + let (faces, flip, closed) = match e.name() { + "CLOSED_SHELL" => (refs(p.get(1)), false, true), + "OPEN_SHELL" => (refs(p.get(1)), false, false), + "ORIENTED_CLOSED_SHELL" | "ORIENTED_OPEN_SHELL" => { + let inner = entity(cx, rid(p.get(2), "oriented shell")?)?; + (refs(inner.params().get(1)), !p.get(3).and_then(Param::as_bool).unwrap_or(true), e.name() == "ORIENTED_CLOSED_SHELL") + } + n => return Err(format!("unsupported shell {n}")), + }; + if faces.len() > MAX_FACES { + return Err(format!("shell with {} faces", faces.len())); + } + let mut out = Vec::with_capacity(faces.len()); + for f in faces { + match self.face(cx, f, 0) { + Ok(mut face) => { + if flip { + face.invert(); + } + out.push(face); + } + Err(err) => { + if self.faces > MAX_FACES { + return Err(err); + } + failed.push(format!("face #{f}: {err}")); + } + } + } + if out.is_empty() { + return Err("no face could be imported".into()); + } + let read = out.len(); + let out = super::split::split_bands(out, failed); + Ok((out.into(), closed, read)) + } +} + +/// Plane through a polygon (Newell's normal). +fn plane_through(pts: &[mt::Point3]) -> Option { + let n = pts.len(); + let mut nrm = mt::Vector3::new(0.0, 0.0, 0.0); + let mut c = mt::Vector3::new(0.0, 0.0, 0.0); + for i in 0..n { + let (a, b) = (pts.get(i)?, pts.get((i + 1) % n)?); + nrm.x += (a.y - b.y) * (a.z + b.z); + nrm.y += (a.z - b.z) * (a.x + b.x); + nrm.z += (a.x - b.x) * (a.y + b.y); + c += a.to_homogeneous().truncate(); + } + let len = nrm.magnitude(); + if !(len > 1e-12) || n == 0 { + return None; + } + let z = nrm / len; + let o = mt::Point3::from_vec(c / n as f64); + let a = if z.x.abs() < 0.9 { mt::Vector3::unit_x() } else { mt::Vector3::unit_y() }; + let x = (a - z * a.dot(z)).normalize(); + let y = z.cross(x); + Some(mt::Surface::Plane(mt::Plane::new(o, o + x, o + y))) +} + +/// Read a solid-like representation item: MANIFOLD_SOLID_BREP, BREP_WITH_VOIDS, FACETED_BREP, +/// or a SHELL_BASED_SURFACE_MODEL (closed shells give solids, open ones open bodies). +pub(crate) fn read_item(cx: &Ctx, id: u64) -> R { + let e = entity(cx, id)?; + let p = e.params(); + let mut b = Builder::default(); + let mut failed = Vec::new(); + let mut shells = Vec::new(); + let mut closed = true; + let mut file_faces = 0; + match e.name() { + "MANIFOLD_SOLID_BREP" | "FACETED_BREP" | "BREP_WITH_VOIDS" => { + let (s, c, n) = b.shell(cx, rid(p.get(1), "brep outer shell")?, &mut failed)?; + closed &= c; + file_faces += n; + shells.push(s); + if e.name() == "BREP_WITH_VOIDS" { + for v in refs(p.get(2)) { + match b.shell(cx, v, &mut failed) { + Ok((s, c, n)) => { + closed &= c; + file_faces += n; + shells.push(s); + } + Err(err) => failed.push(format!("void shell #{v}: {err}")), + } + } + } + } + "SHELL_BASED_SURFACE_MODEL" => { + for s in refs(p.get(1)) { + match b.shell(cx, s, &mut failed) { + Ok((s, c, n)) => { + closed &= c; + file_faces += n; + shells.push(s); + } + Err(err) => failed.push(format!("shell #{s}: {err}")), + } + } + } + n => return Err(format!("unsupported item {n}")), + } + if shells.is_empty() { + return Err(failed.first().cloned().unwrap_or_else(|| "no shells".into())); + } + let solid = match guard("step solid", || mt::Solid::try_new(shells.clone()).map_err(|e| crate::KernelError::Failed(e.to_string()))) { + Ok(s) => s, + Err(_) => { + closed = false; + mt::Solid::new_unchecked(shells) + } + }; + Ok(BrepOut { solid, failed, closed, file_faces }) +} diff --git a/crates/kernel/src/step_in/geom.rs b/crates/kernel/src/step_in/geom.rs new file mode 100644 index 00000000..c128baf8 --- /dev/null +++ b/crates/kernel/src/step_in/geom.rs @@ -0,0 +1,835 @@ +//! STEP geometry (ISO 10303-42 entities) to kernel curves and surfaces. +//! +//! Points, placements, lines, conics, B-splines (plain and rational), elementary surfaces and +//! swept surfaces. Edge curves are trimmed to their vertices here, so every kernel edge runs +//! exactly from its start vertex to its end vertex. Lengths are scaled to millimetres while +//! reading. + +use std::f64::consts::{FRAC_PI_2, PI, TAU}; + +use mt::{BoundedCurve, Cut, InnerSpace, Invertible, MetricSpace, ParametricCurve, ParametricSurface3D, SearchNearestParameter}; +use truck_modeling as mt; + +use super::Ctx; +use super::p21::{Entity, Param}; + +/// Geometry nesting we follow (trimmed curve of a surface curve of …). +const MAX_GEOM_DEPTH: usize = 8; +/// Most control points in one B-spline (per direction product for surfaces). +const MAX_CTRL: usize = 1_000_000; +const MAX_COORD: f64 = 1e9; + +pub(crate) type R = std::result::Result; + +/// An orthonormal frame (AXIS2_PLACEMENT_3D). +#[derive(Clone, Copy, Debug)] +pub(crate) struct Frame { + pub o: mt::Point3, + pub x: mt::Vector3, + pub y: mt::Vector3, + pub z: mt::Vector3, +} + +impl Frame { + pub fn matrix(&self) -> mt::Matrix4 { + mt::Matrix4::from_cols(self.x.extend(0.0), self.y.extend(0.0), self.z.extend(0.0), self.o.to_homogeneous()) + } + fn at(&self, cx: f64, cy: f64, cz: f64) -> mt::Point3 { + self.o + self.x * cx + self.y * cy + self.z * cz + } +} + +/// A curve before trimming. +#[derive(Clone, Debug)] +pub(crate) enum CurveGeo { + Line { + p: mt::Point3, + d: mt::Vector3, + }, + /// Ellipse (`a == b` for circles) in the frame's XY plane. + Conic { + f: Frame, + a: f64, + b: f64, + }, + BSpline(mt::BSplineCurve), + Nurbs(mt::NurbsCurve), +} + +pub(crate) fn entity<'a>(cx: &Ctx<'a>, id: u64) -> R<&'a Entity> { + cx.ex.get(id).ok_or_else(|| format!("missing entity #{id}")) +} + +fn ent_of<'a>(cx: &Ctx<'a>, p: Option<&Param>, what: &str) -> R<(u64, &'a Entity)> { + let id = p.and_then(Param::as_ref_id).ok_or_else(|| format!("{what}: expected a reference"))?; + Ok((id, entity(cx, id)?)) +} + +fn num(p: Option<&Param>, what: &str) -> R { + let x = p.and_then(Param::as_f64).ok_or_else(|| format!("{what}: expected a number"))?; + if x.is_finite() && x.abs() < MAX_COORD { Ok(x) } else { Err(format!("{what}: {x} out of range")) } +} + +fn length(cx: &Ctx, p: Option<&Param>, what: &str) -> R { + Ok(num(p, what)? * cx.len) +} + +pub(crate) fn point(cx: &Ctx, id: u64) -> R { + let e = entity(cx, id)?; + if e.name() != "CARTESIAN_POINT" { + return Err(format!("#{id}: expected CARTESIAN_POINT, found {}", e.name())); + } + let c = e.params().get(1).and_then(Param::as_list).ok_or("CARTESIAN_POINT without coordinates")?; + let g = |i: usize| -> R { if i < c.len() { num(c.get(i), "coordinate") } else { Ok(0.0) } }; + if c.is_empty() || c.len() > 3 { + return Err(format!("#{id}: point with {} coordinates", c.len())); + } + Ok(mt::Point3::new(g(0)?, g(1)?, g(2)?) * cx.len) +} + +fn direction(cx: &Ctx, id: u64) -> R { + let e = entity(cx, id)?; + if e.name() != "DIRECTION" { + return Err(format!("#{id}: expected DIRECTION, found {}", e.name())); + } + let c = e.params().get(1).and_then(Param::as_list).ok_or("DIRECTION without ratios")?; + let g = |i: usize| -> R { if i < c.len() { num(c.get(i), "direction") } else { Ok(0.0) } }; + let v = mt::Vector3::new(g(0)?, g(1)?, g(2)?); + let m = v.magnitude(); + if m > 1e-12 && m.is_finite() { Ok(v / m) } else { Err(format!("#{id}: zero direction")) } +} + +fn any_perpendicular(z: mt::Vector3) -> mt::Vector3 { + let a = if z.x.abs() < 0.9 { mt::Vector3::unit_x() } else { mt::Vector3::unit_y() }; + let x = a - z * a.dot(z); + x / x.magnitude() +} + +/// AXIS2_PLACEMENT_3D (or 2D, as a frame in the XY plane). +pub(crate) fn frame(cx: &Ctx, id: u64) -> R { + let e = entity(cx, id)?; + let p = e.params(); + let o = point(cx, p.get(1).and_then(Param::as_ref_id).ok_or("placement without location")?)?; + let (z, rf) = match e.name() { + "AXIS2_PLACEMENT_3D" => { + let z = match p.get(2).and_then(Param::as_ref_id) { + Some(d) => direction(cx, d)?, + None => mt::Vector3::unit_z(), + }; + (z, p.get(3).and_then(Param::as_ref_id)) + } + "AXIS2_PLACEMENT_2D" => (mt::Vector3::unit_z(), p.get(2).and_then(Param::as_ref_id)), + n => return Err(format!("#{id}: expected AXIS2_PLACEMENT_3D, found {n}")), + }; + let x = match rf { + Some(d) => { + let r = direction(cx, d)?; + let x = r - z * r.dot(z); + if x.magnitude() > 1e-9 { x / x.magnitude() } else { any_perpendicular(z) } + } + None => any_perpendicular(z), + }; + Ok(Frame { o, x, y: z.cross(x), z }) +} + +/// AXIS1_PLACEMENT: origin and unit axis. +fn axis1(cx: &Ctx, id: u64) -> R<(mt::Point3, mt::Vector3)> { + let e = entity(cx, id)?; + let p = e.params(); + let o = point(cx, p.get(1).and_then(Param::as_ref_id).ok_or("axis without location")?)?; + let d = match p.get(2).and_then(Param::as_ref_id) { + Some(d) => direction(cx, d)?, + None => mt::Vector3::unit_z(), + }; + Ok((o, d)) +} + +fn knot_vec(mults: Option<&Param>, knots: Option<&Param>) -> R { + let m: Vec = mults + .and_then(Param::as_list) + .ok_or("B-spline without knot multiplicities")? + .iter() + .map(|x| x.as_i64().filter(|v| (1..=64).contains(v)).map(|v| v as usize).ok_or_else(|| "bad knot multiplicity".to_string())) + .collect::>()?; + let k: Vec = knots.and_then(Param::as_list).ok_or("B-spline without knots")?.iter().map(|x| num(Some(x), "knot")).collect::>()?; + if m.len() != k.len() || k.is_empty() || k.windows(2).any(|w| w.first() > w.get(1)) { + return Err("inconsistent knot vector".into()); + } + mt::KnotVec::from_single_multi(k, m).map_err(|e| format!("knot vector: {e}")) +} + +fn generated_knots(kind: &str, n_ctrl: usize, degree: usize) -> R { + if n_ctrl <= degree { + return Err("too few control points".into()); + } + let n = n_ctrl + degree + 1; + let v: Vec = match kind { + "BEZIER" => { + if n_ctrl != degree + 1 { + // Piecewise Bézier: segments of `degree` with full interior multiplicity. + let segs = (n_ctrl - 1) / degree.max(1); + if segs * degree + 1 != n_ctrl { + return Err("Bezier control point count".into()); + } + let mut v = vec![0.0; degree + 1]; + for s in 1..segs { + v.extend(std::iter::repeat_n(s as f64, degree)); + } + v.extend(std::iter::repeat_n(segs as f64, degree + 1)); + v + } else { + let mut v = vec![0.0; degree + 1]; + v.extend(std::iter::repeat_n(1.0, degree + 1)); + v + } + } + "UNIFORM" => (0..n).map(|i| i as f64 - degree as f64).collect(), + _ => { + // Quasi-uniform: clamped ends, uniform interior. + let inner = n_ctrl - degree; + let mut v = vec![0.0; degree + 1]; + v.extend((1..inner).map(|i| i as f64)); + v.extend(std::iter::repeat_n(inner as f64, degree + 1)); + v + } + }; + Ok(mt::KnotVec::from(v)) +} + +fn weights_list(p: Option<&Param>) -> R> { + p.and_then(Param::as_list) + .ok_or("rational B-spline without weights")? + .iter() + .map(|w| num(Some(w), "weight").and_then(|w| if w > 0.0 { Ok(w) } else { Err("non-positive weight".into()) })) + .collect() +} + +fn bspline_curve(cx: &Ctx, e: &Entity) -> R { + // Common part (B_SPLINE_CURVE): degree, control points, form, closed, self-intersect. + let (base, knots, kind): (Vec, Option<(Option<&Param>, Option<&Param>)>, &str) = if e.is_complex() { + let b = e.record("B_SPLINE_CURVE").ok_or("complex curve without B_SPLINE_CURVE")?; + let mut base = vec![Param::Str(String::new())]; + base.extend(b.iter().cloned()); + let k = e.record("B_SPLINE_CURVE_WITH_KNOTS").map(|k| (k.first(), k.get(1))); + let kind = if e.has("BEZIER_CURVE") { + "BEZIER" + } else if e.has("UNIFORM_CURVE") { + "UNIFORM" + } else { + "QUASI" + }; + (base, k, kind) + } else { + let p = e.params(); + let k = (e.name() == "B_SPLINE_CURVE_WITH_KNOTS").then(|| (p.get(6), p.get(7))); + let kind = match e.name() { + "BEZIER_CURVE" => "BEZIER", + "UNIFORM_CURVE" => "UNIFORM", + _ => "QUASI", + }; + (p.to_vec(), k, kind) + }; + let degree = base.get(1).and_then(Param::as_i64).filter(|d| (1..=32).contains(d)).ok_or("bad B-spline degree")? as usize; + let ctrl_ids = base.get(2).and_then(Param::as_list).ok_or("B-spline without control points")?; + if ctrl_ids.len() > MAX_CTRL { + return Err("B-spline with too many control points".into()); + } + let ctrl: Vec = ctrl_ids.iter().map(|c| point(cx, c.as_ref_id().ok_or("bad control point")?)).collect::>()?; + let kv = match knots { + Some((m, k)) => knot_vec(m, k)?, + None => generated_knots(kind, ctrl.len(), degree)?, + }; + if kv.len() != ctrl.len() + degree + 1 { + return Err(format!("B-spline knot count {} does not match {} control points of degree {degree}", kv.len(), ctrl.len())); + } + let bsp = mt::BSplineCurve::try_new(kv, ctrl).map_err(|e| format!("B-spline: {e}"))?; + match e.record("RATIONAL_B_SPLINE_CURVE") { + Some(r) => { + let w = weights_list(r.first())?; + let n = mt::NurbsCurve::::try_from_bspline_and_weights(bsp, w).map_err(|e| format!("rational B-spline: {e}"))?; + Ok(CurveGeo::Nurbs(n)) + } + None => Ok(CurveGeo::BSpline(bsp)), + } +} + +/// A STEP curve entity (unbounded or bounded). +pub(crate) fn curve(cx: &Ctx, id: u64, depth: usize) -> R { + if depth > MAX_GEOM_DEPTH { + return Err("curve definitions nested too deeply".into()); + } + let e = entity(cx, id)?; + let p = e.params(); + if e.is_complex() { + return bspline_curve(cx, e); + } + match e.name() { + "LINE" => { + let pt = point(cx, p.get(1).and_then(Param::as_ref_id).ok_or("LINE without point")?)?; + let (_, v) = ent_of(cx, p.get(2), "LINE vector")?; + let d = direction(cx, v.params().get(1).and_then(Param::as_ref_id).ok_or("VECTOR without direction")?)?; + Ok(CurveGeo::Line { p: pt, d }) + } + "CIRCLE" => { + let f = frame(cx, p.get(1).and_then(Param::as_ref_id).ok_or("CIRCLE without placement")?)?; + let r = length(cx, p.get(2), "radius")?; + if r <= 0.0 { + return Err("circle radius must be positive".into()); + } + Ok(CurveGeo::Conic { f, a: r, b: r }) + } + "ELLIPSE" => { + let f = frame(cx, p.get(1).and_then(Param::as_ref_id).ok_or("ELLIPSE without placement")?)?; + let (a, b) = (length(cx, p.get(2), "semi axis")?, length(cx, p.get(3), "semi axis")?); + if a <= 0.0 || b <= 0.0 { + return Err("ellipse axes must be positive".into()); + } + Ok(CurveGeo::Conic { f, a, b }) + } + "B_SPLINE_CURVE_WITH_KNOTS" | "BEZIER_CURVE" | "UNIFORM_CURVE" | "QUASI_UNIFORM_CURVE" => bspline_curve(cx, e), + "POLYLINE" => { + let pts: Vec = p + .get(1) + .and_then(Param::as_list) + .ok_or("POLYLINE without points")? + .iter() + .map(|c| point(cx, c.as_ref_id().ok_or("bad polyline point")?)) + .collect::>()?; + if pts.len() < 2 || pts.len() > MAX_CTRL { + return Err("polyline point count".into()); + } + let mut k = vec![0.0]; + k.extend((0..pts.len()).map(|i| i as f64)); + k.push((pts.len() - 1) as f64); + Ok(CurveGeo::BSpline(mt::BSplineCurve::try_new(mt::KnotVec::from(k), pts).map_err(|e| format!("polyline: {e}"))?)) + } + "TRIMMED_CURVE" => curve(cx, p.get(1).and_then(Param::as_ref_id).ok_or("TRIMMED_CURVE without basis")?, depth + 1), + "SURFACE_CURVE" | "SEAM_CURVE" | "INTERSECTION_CURVE" | "BOUNDED_SURFACE_CURVE" => { + curve(cx, p.get(1).and_then(Param::as_ref_id).ok_or("surface curve without 3D curve")?, depth + 1) + } + n => Err(format!("unsupported curve {n}")), + } +} + +/// Rational quadratic NURBS arc of an ellipse from angle `t0` over `sweep > 0`. +pub(crate) fn conic_arc(f: &Frame, a: f64, b: f64, t0: f64, sweep: f64) -> mt::NurbsCurve { + let n = ((sweep / FRAC_PI_2).ceil() as usize).clamp(1, 8); + let dt = sweep / n as f64; + let w = (dt * 0.5).cos(); + let at = |c: f64, s: f64| f.at(a * c, b * s, 0.0); + let mut ctrl = Vec::with_capacity(2 * n + 1); + let p0 = at(t0.cos(), t0.sin()); + ctrl.push(p0.to_homogeneous()); + let mut knots = vec![0.0, 0.0, 0.0]; + for i in 0..n { + let ta = t0 + dt * i as f64; + let tm = ta + dt * 0.5; + let tb = ta + dt; + let m = at(tm.cos() / w, tm.sin() / w); + ctrl.push(mt::Vector4::new(m.x * w, m.y * w, m.z * w, w)); + ctrl.push(at(tb.cos(), tb.sin()).to_homogeneous()); + let k = (i + 1) as f64 / n as f64; + knots.push(k); + knots.push(k); + } + // Last knot pair becomes the clamped end (multiplicity 3). + knots.push(1.0); + mt::NurbsCurve::new(mt::BSplineCurve::new(mt::KnotVec::from(knots), ctrl)) +} + +fn conic_angle(f: &Frame, a: f64, b: f64, p: mt::Point3) -> f64 { + let d = p - f.o; + (d.dot(f.y) / b).atan2(d.dot(f.x) / a) +} + +/// Parameter of the curve point nearest to `p`, and its distance. +fn nearest_param(c: &C, p: mt::Point3) -> Option<(f64, f64)> +where + C: ParametricCurve + BoundedCurve + SearchNearestParameter, +{ + let (t0, t1) = c.range_tuple(); + let n = 256; + let mut best = (t0, f64::INFINITY); + for i in 0..=n { + let t = t0 + (t1 - t0) * i as f64 / n as f64; + let d = c.subs(t).distance(p); + if d < best.1 { + best = (t, d); + } + } + if let Some(t) = c.search_nearest_parameter(p, best.0, 64) { + let t = t.clamp(t0, t1); + let d = c.subs(t).distance(p); + if d < best.1 { + best = (t, d); + } + } + best.1.is_finite().then_some(best) +} + +/// Sub-curve on `[a, b]` (inside the curve's range). +fn sub_curve(c: &C, a: f64, b: f64) -> C { + let (t0, t1) = c.range_tuple(); + let eps = (t1 - t0).abs() * 1e-9; + let mut out = c.clone(); + if a > t0 + eps { + out = out.cut(a); + } + if b < t1 - eps { + let _ = out.cut(b); + } + out +} + +/// Pieces of a spline edge from `ps` to `pe` (two pieces for a closed edge). +fn spline_pieces(c: &C, ps: mt::Point3, pe: mt::Point3, sense: bool, closed: bool, tol: f64) -> R> +where + C: ParametricCurve + BoundedCurve + SearchNearestParameter + Cut + Clone + Invertible, +{ + let (t0, t1) = c.range_tuple(); + let (ts, ds) = nearest_param(c, ps).ok_or("spline evaluation failed")?; + let (te, de) = nearest_param(c, pe).ok_or("spline evaluation failed")?; + if ds > tol || de > tol { + return Err(format!("edge vertices are {:.3e} mm off their curve", ds.max(de))); + } + let span = t1 - t0; + let near = |a: f64, b: f64| (a - b).abs() <= span * 1e-7; + let mut pieces = if closed { + // Whole closed curve starting at the vertex. + if near(ts, t0) || near(ts, t1) { + let m = t0 + span * 0.5; + vec![sub_curve(c, t0, m), sub_curve(c, m, t1)] + } else { + vec![sub_curve(c, ts, t1), sub_curve(c, t0, ts)] + } + } else { + let (mut a, mut b) = if sense { (ts, te) } else { (te, ts) }; + // A closed curve traversed across its seam. + if a >= b && c.subs(t0).distance(c.subs(t1)) <= tol { + if near(b, t0) { + b = t1; + } else if near(a, t1) { + a = t0; + } else { + let mut v = vec![sub_curve(c, a, t1), sub_curve(c, t0, b)]; + if !sense { + v.reverse(); + v.iter_mut().for_each(|x| x.invert()); + } + return Ok(v); + } + } + if b - a <= span * 1e-9 { + return Err("degenerate spline edge".into()); + } + vec![sub_curve(c, a, b)] + }; + if !sense { + pieces.reverse(); + pieces.iter_mut().for_each(|x| x.invert()); + } + Ok(pieces) +} + +/// Kernel curves for an edge from `ps` to `pe`: one piece, or two for a closed edge (so no +/// kernel edge starts and ends at the same vertex). +pub(crate) fn edge_pieces(g: &CurveGeo, ps: mt::Point3, pe: mt::Point3, sense: bool, closed: bool, tol: f64) -> R> { + match g { + CurveGeo::Line { .. } => { + if closed { + return Err("closed straight edge".into()); + } + Ok(vec![mt::Curve::Line(mt::Line(ps, pe))]) + } + CurveGeo::Conic { f, a, b } => { + let (a, b) = (*a, *b); + let ts = conic_angle(f, a, b, ps); + let te = conic_angle(f, a, b, pe); + let mut sweep = if sense { te - ts } else { ts - te }; + sweep = sweep.rem_euclid(TAU); + if closed || sweep < 1e-12 { + sweep = TAU; + } + let pieces: Vec<(f64, f64)> = if sweep > PI * 1.5 { vec![(0.0, sweep * 0.5), (sweep * 0.5, sweep)] } else { vec![(0.0, sweep)] }; + Ok(pieces + .into_iter() + .map(|(s0, s1)| { + if sense { + mt::Curve::NurbsCurve(conic_arc(f, a, b, ts + s0, s1 - s0)) + } else { + let mut c = conic_arc(f, a, b, ts - s1, s1 - s0); + c.invert(); + mt::Curve::NurbsCurve(c) + } + }) + .collect()) + } + CurveGeo::BSpline(c) => Ok(spline_pieces(c, ps, pe, sense, closed, tol)?.into_iter().map(mt::Curve::BSplineCurve).collect()), + CurveGeo::Nurbs(c) => Ok(spline_pieces(c, ps, pe, sense, closed, tol)?.into_iter().map(mt::Curve::NurbsCurve).collect()), + } +} + +/// A bounded kernel curve for a swept surface's profile (lines are bounded by `extent`). +fn bounded_curve(g: &CurveGeo, extent: &[mt::Point3]) -> R { + Ok(match g { + CurveGeo::Line { p, d } => { + let (lo, hi) = span_along(extent, *p, *d); + mt::Curve::Line(mt::Line(*p + *d * lo, *p + *d * hi)) + } + CurveGeo::Conic { f, a, b } => mt::Curve::NurbsCurve(conic_arc(f, *a, *b, 0.0, TAU)), + CurveGeo::BSpline(c) => mt::Curve::BSplineCurve(c.clone()), + CurveGeo::Nurbs(c) => mt::Curve::NurbsCurve(c.clone()), + }) +} + +/// Range of `(q - o)·d` over the points, widened by a margin. +fn span_along(pts: &[mt::Point3], o: mt::Point3, d: mt::Vector3) -> (f64, f64) { + let (mut lo, mut hi) = (f64::INFINITY, f64::NEG_INFINITY); + for q in pts { + let t = (q - o).dot(d); + lo = lo.min(t); + hi = hi.max(t); + } + if !(lo.is_finite() && hi.is_finite()) { + return (-1.0, 1.0); + } + let m = (hi - lo) * 0.25 + 1e-3; + (lo - m, hi + m) +} + +/// The angle shared by all the given angles (within 1e-4 rad), if there are any. +fn common_angle(angles: impl Iterator) -> Option { + let (mut sx, mut sy, mut n) = (0.0, 0.0, 0usize); + let mut all = Vec::new(); + for a in angles { + sx += a.cos(); + sy += a.sin(); + n += 1; + all.push(a); + } + if n == 0 { + return None; + } + let m = sy.atan2(sx); + let close = all.iter().all(|a| { + let d = (a - m).rem_euclid(TAU); + d.min(TAU - d) < 1e-4 + }); + close.then_some(m) +} + +/// Middle and width of the largest gap between angles (radians, any range). +fn largest_gap(angles: &[f64]) -> Option<(f64, f64)> { + let mut a: Vec = angles.iter().map(|x| x.rem_euclid(TAU)).filter(|x| x.is_finite()).collect(); + if a.is_empty() { + return None; + } + a.sort_by(f64::total_cmp); + let (mut best, mut mid) = (0.0, 0.0); + for i in 0..a.len() { + let (Some(&x), Some(&y)) = (a.get(i), a.get((i + 1) % a.len())) else { continue }; + let gap = if i + 1 == a.len() { y + TAU - x } else { y - x }; + if gap > best { + best = gap; + mid = x + gap * 0.5; + } + } + Some((mid, best)) +} + +/// The frame turned about its Z axis so that its X axis (the parameter seam of a revolved +/// surface) lies in the widest angular gap of the face boundary. +fn seam_frame(f: &Frame, extent: &[mt::Point3], seam: &[mt::Point3]) -> Frame { + // A seam edge of the face along the rotation: the parameter seam must be there. + if let Some(m) = common_angle(seam.iter().filter_map(|p| { + let d = p - f.o; + let (x, y) = (d.dot(f.x), d.dot(f.y)); + (x.hypot(y) > 1e-6 * (1.0 + d.magnitude())).then(|| y.atan2(x)) + })) { + let x = f.x * m.cos() + f.y * m.sin(); + return Frame { o: f.o, x, y: f.z.cross(x), z: f.z }; + } + let az: Vec = extent + .iter() + .filter_map(|p| { + let d = p - f.o; + let (x, y) = (d.dot(f.x), d.dot(f.y)); + (x.hypot(y) > 1e-9).then(|| y.atan2(x)) + }) + .collect(); + match largest_gap(&az) { + Some((m, gap)) if gap > 1e-3 => { + let x = f.x * m.cos() + f.y * m.sin(); + Frame { o: f.o, x, y: f.z.cross(x), z: f.z } + } + _ => *f, + } +} + +/// Polar axis for a sphere: far from every boundary point, and (for a face concentrated on +/// one side) with both poles outside the face, so every loop closes in parameter space. +fn sphere_axis(f: &Frame, extent: &[mt::Point3]) -> mt::Vector3 { + let dirs: Vec = extent + .iter() + .filter_map(|p| { + let d = p - f.o; + let m = d.magnitude(); + (m > 1e-12).then(|| d / m) + }) + .collect(); + if dirs.is_empty() { + return f.z; + } + let mean = dirs.iter().fold(mt::Vector3::new(0.0, 0.0, 0.0), |a, d| a + d) / dirs.len() as f64; + let concentrated = mean.magnitude() > 0.3; + let mhat = if concentrated { mean / mean.magnitude() } else { mean }; + let n = 400; + let golden = PI * (3.0 - 5f64.sqrt()); + let mut best = (f64::NEG_INFINITY, f.z); + // Half of a Fibonacci sphere (an axis and its opposite are the same poles). + for i in 0..n { + let z = 1.0 - (i as f64 + 0.5) / n as f64; + let r = (1.0 - z * z).max(0.0).sqrt(); + let t = golden * i as f64; + let a = mt::Vector3::new(r * t.cos(), r * t.sin(), z); + if concentrated && a.dot(mhat).abs() > 0.35 { + continue; + } + let score = dirs.iter().map(|d| 1.0 - d.dot(a).abs()).fold(f64::INFINITY, f64::min); + if score > best.0 { + best = (score, a); + } + } + best.1 +} + +/// Tube angle of a point on a torus. +fn tube_angle(f: &Frame, big: f64, p: mt::Point3) -> f64 { + let d = p - f.o; + let z = d.dot(f.z); + let rho = (d - f.z * z).magnitude(); + z.atan2(rho - big) +} + +fn revolved(curve: mt::Curve, o: mt::Point3, axis: mt::Vector3) -> mt::Surface { + mt::Surface::RevolutedCurve(mt::Processor::new(mt::RevolutedCurve::by_revolution(curve, o, axis))) +} + +/// A surface of revolution whose normal at the profile point `(u, 0)` points along `expected`. +/// The rotation is reversed rather than the surface inverted: truck's inverted `Processor` +/// swaps the parameters but not the search hints, which breaks meshing. +fn revolved_oriented(curve: mt::Curve, o: mt::Point3, axis: mt::Vector3, u: f64, expected: mt::Vector3) -> mt::Surface { + let s = revolved(curve.clone(), o, axis); + if s.normal(u, 0.0).dot(expected) >= 0.0 { + return s; + } + revolved(curve, o, -axis) +} + +fn bspline_surface(cx: &Ctx, e: &Entity) -> R { + let (base, knots, kind): (Vec, Option<[Option<&Param>; 4]>, &str) = if e.is_complex() { + let b = e.record("B_SPLINE_SURFACE").ok_or("complex surface without B_SPLINE_SURFACE")?; + let mut base = vec![Param::Str(String::new())]; + base.extend(b.iter().cloned()); + let k = e.record("B_SPLINE_SURFACE_WITH_KNOTS").map(|k| [k.first(), k.get(1), k.get(2), k.get(3)]); + let kind = if e.has("BEZIER_SURFACE") { + "BEZIER" + } else if e.has("UNIFORM_SURFACE") { + "UNIFORM" + } else { + "QUASI" + }; + (base, k, kind) + } else { + let p = e.params(); + let k = (e.name() == "B_SPLINE_SURFACE_WITH_KNOTS").then(|| [p.get(8), p.get(9), p.get(10), p.get(11)]); + let kind = match e.name() { + "BEZIER_SURFACE" => "BEZIER", + "UNIFORM_SURFACE" => "UNIFORM", + _ => "QUASI", + }; + (p.to_vec(), k, kind) + }; + let deg = |i: usize| base.get(i).and_then(Param::as_i64).filter(|d| (1..=32).contains(d)).map(|d| d as usize).ok_or("bad B-spline degree"); + let (du, dv) = (deg(1)?, deg(2)?); + let rows = base.get(3).and_then(Param::as_list).ok_or("B-spline surface without control points")?; + let mut ctrl: Vec> = Vec::with_capacity(rows.len()); + let mut total = 0usize; + for r in rows { + let r = r.as_list().ok_or("bad control point row")?; + total += r.len(); + if total > MAX_CTRL { + return Err("B-spline surface with too many control points".into()); + } + ctrl.push(r.iter().map(|c| point(cx, c.as_ref_id().ok_or("bad control point")?)).collect::>()?); + } + let nu = ctrl.len(); + let nv = ctrl.first().map(Vec::len).unwrap_or(0); + if nu == 0 || ctrl.iter().any(|r| r.len() != nv) { + return Err("ragged B-spline control net".into()); + } + let (ku, kv) = match knots { + Some([mu, mv, ku, kv]) => (knot_vec(mu, ku)?, knot_vec(mv, kv)?), + None => (generated_knots(kind, nu, du)?, generated_knots(kind, nv, dv)?), + }; + if ku.len() != nu + du + 1 || kv.len() != nv + dv + 1 { + return Err("B-spline surface knot counts do not match the control net".into()); + } + let bsp = mt::BSplineSurface::try_new((ku, kv), ctrl).map_err(|e| format!("B-spline surface: {e}"))?; + match e.record("RATIONAL_B_SPLINE_SURFACE") { + Some(r) => { + let w: Vec> = r + .first() + .and_then(Param::as_list) + .ok_or("rational surface without weights")? + .iter() + .map(|row| weights_list(Some(row))) + .collect::>()?; + if w.len() != nu || w.iter().any(|r| r.len() != nv) { + return Err("weights do not match the control net".into()); + } + let n = mt::NurbsSurface::::try_from_bspline_and_weights(bsp, w).map_err(|e| format!("rational surface: {e}"))?; + Ok(mt::Surface::NurbsSurface(n)) + } + None => Ok(mt::Surface::BSplineSurface(bsp)), + } +} + +/// A STEP surface entity as a kernel surface whose normal is the STEP surface normal. +/// `extent` holds points of the face boundary (bounds the analytic surfaces' parameter ranges). +/// `sense` false gives the opposite normal (a face with `same_sense = .F.`). +/// `seam` holds points of the face's seam edges (edges its loops use twice). +pub(crate) fn surface(cx: &Ctx, id: u64, extent: &[mt::Point3], seam: &[mt::Point3], sense: bool, depth: usize) -> R { + let sign = if sense { 1.0 } else { -1.0 }; + let flip = |mut s: mt::Surface| { + if !sense { + s.invert(); + } + s + }; + if depth > MAX_GEOM_DEPTH { + return Err("surface definitions nested too deeply".into()); + } + let e = entity(cx, id)?; + if e.is_complex() { + return bspline_surface(cx, e).map(flip); + } + let p = e.params(); + let placement = || -> R { frame(cx, p.get(1).and_then(Param::as_ref_id).ok_or("surface without placement")?) }; + match e.name() { + "PLANE" => { + let f = placement()?; + Ok(flip(mt::Surface::Plane(mt::Plane::new(f.o, f.o + f.x, f.o + f.y)))) + } + "CYLINDRICAL_SURFACE" => { + let f = seam_frame(&placement()?, extent, seam); + let r = length(cx, p.get(2), "radius")?; + if r <= 0.0 { + return Err("cylinder radius must be positive".into()); + } + let (lo, hi) = span_along(extent, f.o, f.z); + let line = mt::Curve::Line(mt::Line(f.at(r, 0.0, lo), f.at(r, 0.0, hi))); + Ok(revolved_oriented(line, f.o, f.z, 0.5, f.x * sign)) + } + "CONICAL_SURFACE" => { + let f = seam_frame(&placement()?, extent, seam); + let r = length(cx, p.get(2), "radius")?; + let a = num(p.get(3), "semi angle")? * cx.ang; + if r < 0.0 || !(a.abs() > 1e-9 && a.abs() < FRAC_PI_2 - 1e-9) { + return Err("bad cone".into()); + } + let t = a.tan(); + let (mut lo, mut hi) = span_along(extent, f.o, f.z); + // Stop at the apex. + let apex = -r / t; + if t > 0.0 { + lo = lo.max(apex); + } else { + hi = hi.min(apex); + } + if hi - lo <= 1e-9 { + return Err("cone face has no extent".into()); + } + let line = mt::Curve::Line(mt::Line(f.at(r + lo * t, 0.0, lo), f.at(r + hi * t, 0.0, hi))); + Ok(revolved_oriented(line, f.o, f.z, 0.5, (f.x * a.cos() - f.z * a.sin()) * sign)) + } + "SPHERICAL_SURFACE" => { + let f = placement()?; + let r = length(cx, p.get(2), "radius")?; + if r <= 0.0 { + return Err("sphere radius must be positive".into()); + } + // With seam edges the file's own parameterisation fits the loops; else pick poles + // away from the boundary. + let f = if seam.is_empty() { + let z = sphere_axis(&f, extent); + let x = any_perpendicular(z); + seam_frame(&Frame { o: f.o, x, y: z.cross(x), z }, extent, seam) + } else { + seam_frame(&f, extent, seam) + }; + let half = Frame { o: f.o, x: f.x, y: f.z, z: -f.y }; + let arc = conic_arc(&half, r, r, -FRAC_PI_2, PI); + let um = (arc.range_tuple().0 + arc.range_tuple().1) * 0.5; + Ok(revolved_oriented(mt::Curve::NurbsCurve(arc), f.o, f.z, um, f.x * sign)) + } + "TOROIDAL_SURFACE" => { + let f = placement()?; + let (big, small) = (length(cx, p.get(2), "major radius")?, length(cx, p.get(3), "minor radius")?); + if small <= 0.0 || big < 0.0 { + return Err("bad torus radii".into()); + } + // Seam edges either run around the axis (constant tube angle) or around the tube + // (constant azimuth); each fixes one parameter seam. + let tube_seam: Vec = + seam.iter().copied().filter(|p| common_angle(std::iter::once(tube_angle(&f, big, *p))).is_some()).collect(); + let phis_seam = common_angle(tube_seam.iter().map(|p| tube_angle(&f, big, *p))); + let az_seam: Vec = if phis_seam.is_some() { Vec::new() } else { seam.to_vec() }; + let f = seam_frame(&f, extent, &az_seam); + let phis: Vec = extent.iter().map(|p| tube_angle(&f, big, *p)).collect(); + let phi0 = match (phis_seam, largest_gap(&phis)) { + (Some(m), _) => m, + (None, Some((m, _))) => m, + (None, None) => 0.0, + }; + let tube = Frame { o: f.at(big, 0.0, 0.0), x: f.x, y: f.z, z: -f.y }; + let circle = conic_arc(&tube, small, small, phi0, TAU); + let u0 = circle.range_tuple().0; + let outward = f.x * phi0.cos() + f.z * phi0.sin(); + Ok(revolved_oriented(mt::Curve::NurbsCurve(circle), f.o, f.z, u0, outward * sign)) + } + "B_SPLINE_SURFACE_WITH_KNOTS" | "BEZIER_SURFACE" | "UNIFORM_SURFACE" | "QUASI_UNIFORM_SURFACE" => bspline_surface(cx, e).map(flip), + "SURFACE_OF_REVOLUTION" => { + let g = curve(cx, p.get(1).and_then(Param::as_ref_id).ok_or("revolution without profile")?, depth + 1)?; + let (o, axis) = axis1(cx, p.get(2).and_then(Param::as_ref_id).ok_or("revolution without axis")?)?; + // Reversed sense: turn the other way round the axis. + Ok(revolved(bounded_curve(&g, extent)?, o, axis * sign)) + } + "SURFACE_OF_LINEAR_EXTRUSION" => { + let g = curve(cx, p.get(1).and_then(Param::as_ref_id).ok_or("extrusion without profile")?, depth + 1)?; + let (_, v) = ent_of(cx, p.get(2), "extrusion vector")?; + let d = direction(cx, v.params().get(1).and_then(Param::as_ref_id).ok_or("VECTOR without direction")?)?; + let nurbs: mt::NurbsCurve = match bounded_curve(&g, extent)? { + mt::Curve::Line(l) => mt::BSplineCurve::new(mt::KnotVec::bezier_knot(1), vec![l.0, l.1]).into(), + mt::Curve::BSplineCurve(c) => c.into(), + mt::Curve::NurbsCurve(c) => c, + mt::Curve::IntersectionCurve(_) => return Err("unsupported extrusion profile".into()), + }; + let (t0, t1) = nurbs.range_tuple(); + let profile: Vec = (0..=32).map(|i| nurbs.subs(t0 + (t1 - t0) * i as f64 / 32.0)).collect(); + let (plo, phi) = span_along(&profile, mt::Point3::new(0.0, 0.0, 0.0), d); + let (elo, ehi) = span_along(extent, mt::Point3::new(0.0, 0.0, 0.0), d); + let (lo, hi) = (elo - phi, ehi - plo); + let bsp = nurbs.into_non_rationalized(); + let shift = |w: &mt::Vector4, s: f64| mt::Vector4::new(w.x + d.x * s * w.w, w.y + d.y * s * w.w, w.z + d.z * s * w.w, w.w); + let ctrl: Vec> = bsp.control_points().iter().map(|w| vec![shift(w, lo), shift(w, hi)]).collect(); + let s = + mt::BSplineSurface::try_new((bsp.knot_vec().clone(), mt::KnotVec::bezier_knot(1)), ctrl).map_err(|e| format!("extrusion: {e}"))?; + Ok(flip(mt::Surface::NurbsSurface(mt::NurbsSurface::new(s)))) + } + "RECTANGULAR_TRIMMED_SURFACE" | "CURVE_BOUNDED_SURFACE" => { + surface(cx, p.get(1).and_then(Param::as_ref_id).ok_or("trimmed surface without basis")?, extent, seam, sense, depth + 1) + } + n => Err(format!("unsupported surface {n}")), + } +} diff --git a/crates/kernel/src/step_in/mod.rs b/crates/kernel/src/step_in/mod.rs new file mode 100644 index 00000000..06d9314b --- /dev/null +++ b/crates/kernel/src/step_in/mod.rs @@ -0,0 +1,638 @@ +//! STEP import (ISO 10303-21 files, AP203/AP214/AP242 geometry). +//! +//! Our own reader, written from the public standards: [`p21`] parses the exchange structure, +//! [`geom`] and [`brep`] turn geometry and topology into kernel solids, and this module walks +//! the product structure (products, shape representations, assemblies with their placements, +//! units, colours). STEP files are untrusted: every count and nesting depth is capped, every +//! truck call runs under [`guard`], and what cannot be read becomes a warning rather than a +//! failure wherever the rest of the file is still usable. + +mod brep; +mod geom; +pub(crate) mod p21; +mod split; + +use std::collections::{BTreeMap, HashMap, HashSet}; + +use serde::{Deserialize, Serialize}; +use truck_modeling as mt; + +use crate::body::Body; +use crate::{KernelError, Result, guard}; +use geom::R; +use p21::{Entity, Exchange, Param}; + +/// Deepest assembly nesting followed. +const MAX_ASSEMBLY_DEPTH: usize = 32; +/// Most body instances in one import. +const MAX_BODIES: usize = 10_000; + +/// A body read from a STEP file. +#[derive(Clone, Debug)] +pub struct ImportedBody { + /// Solid name (the brep's label), else the product name. + pub name: String, + pub body: Body, + /// Surface colour (linear-ish RGB 0..1) from the file's presentation styles. + pub color: Option<[f32; 3]>, + /// Product names from the root to the part that holds the body. + pub path: Vec, + /// The body is a closed solid (false: faces were left out, or the file had a surface model). + pub closed: bool, + /// Faces read from the file. The body can have more: a face that wraps all the way around + /// a closed surface is split in two for the kernel. + pub file_faces: usize, +} + +/// One product occurrence in the file's assembly tree. +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] +pub struct StepNode { + pub name: String, + /// Placement in the parent (column-major 4×4, millimetres). + pub transform: [[f64; 4]; 4], + /// Indices into [`StepImport::bodies`] of this product's own bodies. + #[serde(default, skip_serializing_if = "Vec::is_empty")] + pub bodies: Vec, + #[serde(default, skip_serializing_if = "Vec::is_empty")] + pub children: Vec, +} + +/// Result of [`step_import`]. +#[derive(Clone, Debug)] +pub struct StepImport { + pub bodies: Vec, + /// Root products (an assembly tree; a plain part is a single root). + pub tree: Vec, + /// What could not be read (unsupported entities, faces left out…). + pub warnings: Vec, + /// FILE_SCHEMA (e.g. `AUTOMOTIVE_DESIGN`, `CONFIG_CONTROL_DESIGN`, `AP242_…`). + pub schema: String, + /// FILE_NAME originating system. + pub system: String, +} + +/// Reading context: the file and the units of the representation being read. +pub(crate) struct Ctx<'a> { + pub ex: &'a Exchange, + /// Millimetres per file length unit. + pub len: f64, + /// Radians per file angle unit. + pub ang: f64, + /// Vertex-on-curve tolerance (mm). + pub tol: f64, +} + +#[derive(Clone, Copy)] +struct Units { + len: f64, + ang: f64, + tol: f64, +} + +const DEFAULT_UNITS: Units = Units { len: 1.0, ang: 1.0, tol: 0.01 }; + +fn si_prefix(p: Option<&Param>) -> f64 { + match p.and_then(Param::as_enum) { + Some("EXA") => 1e18, + Some("PETA") => 1e15, + Some("TERA") => 1e12, + Some("GIGA") => 1e9, + Some("MEGA") => 1e6, + Some("KILO") => 1e3, + Some("HECTO") => 1e2, + Some("DECA") => 1e1, + Some("DECI") => 1e-1, + Some("CENTI") => 1e-2, + Some("MILLI") => 1e-3, + Some("MICRO") => 1e-6, + Some("NANO") => 1e-9, + Some("PICO") => 1e-12, + Some("FEMTO") => 1e-15, + Some("ATTO") => 1e-18, + _ => 1.0, + } +} + +/// A unit's size: millimetres for lengths, radians for plane angles. +fn unit_factor(ex: &Exchange, id: u64, depth: usize) -> Option<(&'static str, f64)> { + if depth > 4 { + return None; + } + let e = ex.get(id)?; + let kind = if e.has("LENGTH_UNIT") { + "length" + } else if e.has("PLANE_ANGLE_UNIT") { + "angle" + } else { + "other" + }; + if let Some(si) = e.record("SI_UNIT") { + let base = match si.get(1).and_then(Param::as_enum) { + Some("METRE") => 1000.0, + Some("RADIAN") => 1.0, + _ => return Some((kind, 1.0)), + }; + return Some((kind, si_prefix(si.first()) * base)); + } + if let Some(cb) = e.record("CONVERSION_BASED_UNIT") { + let m = ex.get(cb.get(1)?.as_ref_id()?)?; + let mp = m.params(); + let v = mp.first()?.as_f64()?; + let (_, f) = unit_factor(ex, mp.get(1)?.as_ref_id()?, depth + 1)?; + let x = v * f; + return (x.is_finite() && x > 0.0).then_some((kind, x)); + } + Some((kind, 1.0)) +} + +/// Units of a representation context (GLOBAL_UNIT_ASSIGNED_CONTEXT, uncertainty). +fn context_units(ex: &Exchange, ctx: Option, cache: &mut HashMap) -> Units { + let Some(id) = ctx else { return DEFAULT_UNITS }; + if let Some(u) = cache.get(&id) { + return *u; + } + let mut u = DEFAULT_UNITS; + let mut tol_file = None; + if let Some(e) = ex.get(id) { + if let Some(units) = e.record("GLOBAL_UNIT_ASSIGNED_CONTEXT").and_then(|r| r.first()).and_then(Param::as_list) { + for uid in units.iter().filter_map(Param::as_ref_id) { + match unit_factor(ex, uid, 0) { + Some(("length", f)) if f.is_finite() && f > 0.0 => u.len = f, + Some(("angle", f)) if f.is_finite() && f > 0.0 => u.ang = f, + _ => {} + } + } + } + if let Some(unc) = e.record("GLOBAL_UNCERTAINTY_ASSIGNED_CONTEXT").and_then(|r| r.first()).and_then(Param::as_list) { + for uid in unc.iter().filter_map(Param::as_ref_id) { + if let Some(m) = ex.get(uid) + && let Some(v) = m.params().first().and_then(Param::as_f64) + { + let f = m.params().get(1).and_then(Param::as_ref_id).and_then(|x| unit_factor(ex, x, 0)).map(|x| x.1); + tol_file = Some(v * f.unwrap_or(u.len)); + } + } + } + } + u.tol = tol_file.filter(|t| t.is_finite() && *t > 0.0).map(|t| (t * 10.0).max(1e-3)).unwrap_or(u.len * 1e-2).clamp(1e-4, 1.0); + cache.insert(id, u); + u +} + +/// Items and context of a representation entity. +fn representation(e: &Entity) -> Option<(Vec, Option)> { + let p = if e.is_complex() { + e.record("REPRESENTATION") + .or_else(|| e.records.iter().find(|r| r.name.ends_with("REPRESENTATION") && r.params.len() >= 3).map(|r| r.params.as_slice()))? + } else if e.name().ends_with("REPRESENTATION") { + e.params() + } else { + return None; + }; + let items = p.get(1)?.as_list()?.iter().filter_map(Param::as_ref_id).collect(); + Some((items, p.get(2).and_then(Param::as_ref_id))) +} + +fn mat_mul(a: &mt::Matrix4, b: &mt::Matrix4) -> mt::Matrix4 { + a * b +} + +fn to_array(m: &mt::Matrix4) -> [[f64; 4]; 4] { + [[m.x.x, m.x.y, m.x.z, m.x.w], [m.y.x, m.y.y, m.y.z, m.y.w], [m.z.x, m.z.y, m.z.z, m.z.w], [m.w.x, m.w.y, m.w.z, m.w.w]] +} + +/// Inverse of a rigid placement matrix (orthonormal rotation + translation). +fn rigid_inverse(m: &mt::Matrix4) -> mt::Matrix4 { + use mt::SquareMatrix; + m.invert().unwrap_or_else(mt::Matrix4::identity) +} + +struct Reader<'a> { + ex: &'a Exchange, + units: HashMap, + warnings: BTreeMap, + /// Representation → its items and context. + reps: HashMap, Option)>, + /// Plain representation relationships (same product, no placement). + rep_links: HashMap>, + /// Product definition → its shape representations. + pd_reps: HashMap>, + /// Product definition → child occurrences (NAUO ids). + children: HashMap>, + /// NAUO → placement relationship entity. + placements: HashMap, + /// Styled item → colour. + colors: HashMap, + /// Converted items (brep id → solid, or the error). + solids: HashMap>, + used_items: HashSet, + /// Items already instanced for the product being read. + node_items: HashSet, + bodies: Vec, +} + +impl<'a> Reader<'a> { + fn warn(&mut self, msg: impl Into) { + *self.warnings.entry(msg.into()).or_default() += 1; + } + + fn index(&mut self) { + let ex = self.ex; + let mut pds: HashMap = HashMap::new(); // PDS → definition + let mut sdrs: Vec<(u64, u64)> = Vec::new(); // (PDS, rep) + let mut cdsrs: Vec<(u64, u64)> = Vec::new(); // (relationship, PDS) + let mut nauos: Vec<(u64, u64, u64)> = Vec::new(); // (id, parent, child) + let mut ids: Vec<&u64> = ex.entities.keys().collect(); + ids.sort(); + for &id in ids { + let Some(e) = ex.get(id) else { continue }; + if let Some(r) = representation(e) { + self.reps.insert(id, r); + } + let p = e.params(); + let rf = |i: usize| p.get(i).and_then(Param::as_ref_id); + match e.name() { + "PRODUCT_DEFINITION_SHAPE" | "PROPERTY_DEFINITION" => { + if let Some(d) = rf(2) { + pds.insert(id, d); + } + } + "SHAPE_DEFINITION_REPRESENTATION" => { + if let (Some(d), Some(r)) = (rf(0), rf(1)) { + sdrs.push((d, r)); + } + } + "CONTEXT_DEPENDENT_SHAPE_REPRESENTATION" => { + if let (Some(r), Some(d)) = (rf(0), rf(1)) { + cdsrs.push((r, d)); + } + } + "NEXT_ASSEMBLY_USAGE_OCCURRENCE" | "ASSEMBLY_COMPONENT_USAGE" | "SPECIFIED_HIGHER_USAGE_OCCURRENCE" => { + if let (Some(a), Some(b)) = (rf(3), rf(4)) { + nauos.push((id, a, b)); + } + } + "SHAPE_REPRESENTATION_RELATIONSHIP" | "REPRESENTATION_RELATIONSHIP" if !e.is_complex() => { + if let (Some(a), Some(b)) = (rf(2), rf(3)) { + self.rep_links.entry(a).or_default().push(b); + self.rep_links.entry(b).or_default().push(a); + } + } + "STYLED_ITEM" | "OVER_RIDING_STYLED_ITEM" => { + if let Some(item) = rf(2) + && let Some(c) = self.find_colour(p.get(1)) + { + self.colors.entry(item).or_insert(c); + } + } + _ => {} + } + } + for (d, r) in sdrs { + if let Some(&pd) = pds.get(&d) { + self.pd_reps.entry(pd).or_default().push(r); + } + } + for (r, d) in cdsrs { + if let Some(&nauo) = pds.get(&d) { + self.placements.insert(nauo, r); + } + } + for (id, parent, child) in nauos { + if parent != child { + self.children.entry(parent).or_default().push(id); + } + } + } + + /// First colour reachable from a style list (bounded search). + fn find_colour(&self, styles: Option<&Param>) -> Option<[f32; 3]> { + let mut queue: Vec = styles?.as_list()?.iter().filter_map(Param::as_ref_id).collect(); + let mut seen = HashSet::new(); + while let Some(id) = queue.pop() { + if seen.len() > 64 || !seen.insert(id) { + continue; + } + let e = self.ex.get(id)?; + if e.name() == "COLOUR_RGB" { + let p = e.params(); + let c = |i: usize| p.get(i).and_then(Param::as_f64).map(|x| x.clamp(0.0, 1.0) as f32); + return Some([c(1)?, c(2)?, c(3)?]); + } + if e.name() == "DRAUGHTING_PRE_DEFINED_COLOUR" { + return Some(match e.params().first().and_then(Param::as_str).unwrap_or("") { + "red" => [1.0, 0.0, 0.0], + "green" => [0.0, 1.0, 0.0], + "blue" => [0.0, 0.0, 1.0], + "yellow" => [1.0, 1.0, 0.0], + "magenta" => [1.0, 0.0, 1.0], + "cyan" => [0.0, 1.0, 1.0], + "black" => [0.0, 0.0, 0.0], + _ => [1.0, 1.0, 1.0], + }); + } + fn walk(p: &Param, out: &mut Vec, depth: usize) { + match p { + Param::Ref(r) => out.push(*r), + Param::List(v) | Param::Typed(_, v) if depth < 4 => v.iter().for_each(|x| walk(x, out, depth + 1)), + _ => {} + } + } + for r in &e.records { + for p in &r.params { + walk(p, &mut queue, 0); + } + } + } + None + } + + fn units_of(&mut self, rep: u64) -> Units { + let ctx = self.reps.get(&rep).and_then(|r| r.1); + context_units(self.ex, ctx, &mut self.units) + } + + fn placement_matrix(&mut self, item: u64, u: Units) -> R { + let cx = Ctx { ex: self.ex, len: u.len, ang: u.ang, tol: u.tol }; + let e = geom::entity(&cx, item)?; + if e.name() == "CARTESIAN_TRANSFORMATION_OPERATOR_3D" { + let p = e.params(); + let o = geom::point(&cx, p.get(4).and_then(Param::as_ref_id).ok_or("transformation without origin")?)?; + return Ok(mt::Matrix4::from_translation(o.to_homogeneous().truncate())); + } + Ok(geom::frame(&cx, item)?.matrix()) + } + + /// Representations of a product (its own shape representations and those linked to them). + fn rep_closure(&self, roots: &[u64]) -> Vec { + let mut seen: Vec = Vec::new(); + let mut stack: Vec = roots.to_vec(); + while let Some(r) = stack.pop() { + if seen.contains(&r) || seen.len() > 1024 { + continue; + } + seen.push(r); + if let Some(l) = self.rep_links.get(&r) { + stack.extend(l.iter().copied()); + } + } + seen + } + + fn item_solid(&mut self, item: u64, u: Units) -> Option<(mt::Solid, bool, usize)> { + if let Some(r) = self.solids.get(&item) { + return r.as_ref().ok().cloned(); + } + let cx = Ctx { ex: self.ex, len: u.len, ang: u.ang, tol: u.tol }; + let r = match guard("step brep", || brep::read_item(&cx, item).map_err(KernelError::Failed)) { + Ok(out) => { + let label = self.label(item); + for f in &out.failed { + self.warn(format!("{label}: {f}")); + } + if !out.closed { + self.warn(format!("{label}: imported as an open body (the shell is not closed)")); + } + Ok((out.solid, out.closed, out.file_faces)) + } + Err(e) => { + let label = self.label(item); + let msg = e.to_string().trim_start_matches("the operation failed: ").to_string(); + self.warn(format!("{label}: not imported: {msg}")); + Err(msg) + } + }; + self.solids.insert(item, r.clone()); + r.ok() + } + + fn label(&self, item: u64) -> String { + let name = self.ex.get(item).and_then(|e| e.params().first()).and_then(Param::as_str).unwrap_or("").trim().to_string(); + if name.is_empty() { format!("solid #{item}") } else { format!("solid '{name}'") } + } + + /// Bodies of a representation's items (following mapped items), placed by `m`. + fn rep_bodies(&mut self, rep: u64, m: &mt::Matrix4, path: &[String], product: &str, out: &mut Vec, depth: usize) { + if depth > MAX_ASSEMBLY_DEPTH { + self.warn("mapped items nested too deeply"); + return; + } + let u = self.units_of(rep); + let items = self.reps.get(&rep).map(|r| r.0.clone()).unwrap_or_default(); + for item in items { + let Some(e) = self.ex.get(item) else { continue }; + match e.name() { + "MANIFOLD_SOLID_BREP" | "BREP_WITH_VOIDS" | "FACETED_BREP" | "SHELL_BASED_SURFACE_MODEL" => { + self.used_items.insert(item); + if !self.node_items.insert(item) { + continue; + } + if self.bodies.len() >= MAX_BODIES { + self.warn(format!("more than {MAX_BODIES} bodies: the rest are left out")); + return; + } + let Some((solid, closed, file_faces)) = self.item_solid(item, u) else { continue }; + let identity = *m == mt::Matrix4::from_scale(1.0); + let placed = if identity { + Ok(solid) + } else { + let mm = *m; + guard("step placement", || Ok(mt::builder::transformed(&solid, mm))) + }; + let body = placed.and_then(Body::new); + match body { + Ok(body) => { + let own = e.params().first().and_then(Param::as_str).unwrap_or("").trim().to_string(); + let name = if own.is_empty() { product.to_string() } else { own }; + out.push(self.bodies.len()); + self.bodies.push(ImportedBody { + name, + body, + color: self.colors.get(&item).copied(), + path: path.to_vec(), + closed, + file_faces, + }); + } + Err(err) => { + let label = self.label(item); + self.warn(format!("{label}: not imported: {err}")); + } + } + } + "MAPPED_ITEM" => { + let p = e.params(); + let (Some(map), Some(target)) = (p.get(1).and_then(Param::as_ref_id), p.get(2).and_then(Param::as_ref_id)) else { continue }; + let Some(me) = self.ex.get(map) else { continue }; + let (Some(origin), Some(mapped)) = + (me.params().first().and_then(Param::as_ref_id), me.params().get(1).and_then(Param::as_ref_id)) + else { + continue; + }; + let mu = self.units_of(mapped); + match (self.placement_matrix(target, u), self.placement_matrix(origin, mu)) { + (Ok(t), Ok(o)) => { + let mm = mat_mul(m, &mat_mul(&t, &rigid_inverse(&o))); + for r in self.rep_closure(&[mapped]) { + self.rep_bodies(r, &mm, path, product, out, depth + 1); + } + } + (Err(e), _) | (_, Err(e)) => self.warn(format!("mapped item #{item}: {e}")), + } + } + _ => {} + } + } + } + + fn product_name(&self, pd: u64) -> String { + let name = (|| { + let pdf = self.ex.get(self.ex.get(pd)?.params().get(2)?.as_ref_id()?)?; + let prod = self.ex.get(pdf.params().get(2)?.as_ref_id()?)?; + let p = prod.params(); + let n = p.get(1).and_then(Param::as_str).filter(|s| !s.trim().is_empty()).or_else(|| p.first().and_then(Param::as_str))?; + Some(n.trim().to_string()) + })(); + name.filter(|n| !n.is_empty()).unwrap_or_else(|| "Import".into()) + } + + /// Placement of a child occurrence in its parent. + fn occurrence_matrix(&mut self, nauo: u64, child_reps: &[u64]) -> mt::Matrix4 { + let Some(&rel) = self.placements.get(&nauo) else { return mt::Matrix4::from_scale(1.0) }; + let r = (|| -> R { + let e = self.ex.get(rel).ok_or("missing placement")?; + let rr = e.record("REPRESENTATION_RELATIONSHIP").ok_or("placement without representations")?; + let (r1, r2) = + (rr.get(2).and_then(Param::as_ref_id).ok_or("bad placement")?, rr.get(3).and_then(Param::as_ref_id).ok_or("bad placement")?); + let op = e + .record("REPRESENTATION_RELATIONSHIP_WITH_TRANSFORMATION") + .and_then(|x| x.first()) + .and_then(Param::as_ref_id) + .ok_or("placement without transformation")?; + let ope = self.ex.get(op).ok_or("missing transformation")?; + let (u1, u2) = (self.units_of(r1), self.units_of(r2)); + if ope.name() == "ITEM_DEFINED_TRANSFORMATION" { + let p = ope.params(); + let (i1, i2) = ( + p.get(2).and_then(Param::as_ref_id).ok_or("bad transformation")?, + p.get(3).and_then(Param::as_ref_id).ok_or("bad transformation")?, + ); + let m1 = self.placement_matrix(i1, u1)?; + let m2 = self.placement_matrix(i2, u2)?; + // rep_1 is normally the component: map its item onto the assembly's item. + let child_is_1 = child_reps.contains(&r1) || !child_reps.contains(&r2); + Ok(if child_is_1 { mat_mul(&m2, &rigid_inverse(&m1)) } else { mat_mul(&m1, &rigid_inverse(&m2)) }) + } else { + self.placement_matrix(op, u2) + } + })(); + match r { + Ok(m) => m, + Err(e) => { + self.warn(format!("component placement #{rel}: {e} (placed at the origin)")); + mt::Matrix4::from_scale(1.0) + } + } + } + + fn node(&mut self, pd: u64, local: mt::Matrix4, world: mt::Matrix4, path: &[String], depth: usize, stack: &mut Vec) -> StepNode { + let name = self.product_name(pd); + let mut node = StepNode { name: name.clone(), transform: to_array(&local), bodies: Vec::new(), children: Vec::new() }; + if depth > MAX_ASSEMBLY_DEPTH || stack.contains(&pd) { + self.warn("assembly nested too deeply or recursive: some components were left out"); + return node; + } + stack.push(pd); + let mut path = path.to_vec(); + path.push(name.clone()); + let own = self.pd_reps.get(&pd).cloned().unwrap_or_default(); + let reps = self.rep_closure(&own); + let mut bodies = Vec::new(); + self.node_items.clear(); + for r in &reps { + self.rep_bodies(*r, &world, &path, &name, &mut bodies, 0); + } + node.bodies = bodies; + for nauo in self.children.get(&pd).cloned().unwrap_or_default() { + let Some(child) = self.ex.get(nauo).and_then(|e| e.params().get(4)).and_then(Param::as_ref_id) else { continue }; + let child_own = self.pd_reps.get(&child).cloned().unwrap_or_default(); + let child_reps = self.rep_closure(&child_own); + let m = self.occurrence_matrix(nauo, &child_reps); + let w = mat_mul(&world, &m); + let c = self.node(child, m, w, &path, depth + 1, stack); + node.children.push(c); + } + stack.pop(); + node + } +} + +/// Read a STEP file: its solids (as bodies, placed by the assembly structure) and the tree. +pub fn step_import(text: &str) -> Result { + let ex = p21::parse(text).map_err(KernelError::Invalid)?; + let mut rd = Reader { + ex: &ex, + units: HashMap::new(), + warnings: BTreeMap::new(), + reps: HashMap::new(), + rep_links: HashMap::new(), + pd_reps: HashMap::new(), + children: HashMap::new(), + placements: HashMap::new(), + colors: HashMap::new(), + solids: HashMap::new(), + used_items: HashSet::new(), + node_items: HashSet::new(), + bodies: Vec::new(), + }; + rd.index(); + // Roots: product definitions that are nobody's component. + let child_pds: HashSet = + rd.children.values().flatten().filter_map(|n| ex.get(*n).and_then(|e| e.params().get(4)).and_then(Param::as_ref_id)).collect(); + let mut roots: Vec = + ex.entities.iter().filter(|(id, e)| e.name() == "PRODUCT_DEFINITION" && !child_pds.contains(id)).map(|(id, _)| *id).collect(); + roots.sort(); + let identity = mt::Matrix4::from_scale(1.0); + let mut tree = Vec::new(); + for pd in roots { + let n = rd.node(pd, identity, identity, &[], 0, &mut Vec::new()); + if !n.bodies.is_empty() || !n.children.is_empty() { + tree.push(n); + } + } + // Solids outside any product structure (bare geometry files). + let mut loose: Vec = rd + .reps + .iter() + .flat_map(|(r, (items, _))| items.iter().map(move |i| (*r, *i))) + .filter(|(_, i)| !rd.used_items.contains(i)) + .filter(|(_, i)| { + ex.get(*i).is_some_and(|e| matches!(e.name(), "MANIFOLD_SOLID_BREP" | "BREP_WITH_VOIDS" | "FACETED_BREP" | "SHELL_BASED_SURFACE_MODEL")) + }) + .map(|(r, _)| r) + .collect(); + loose.sort(); + loose.dedup(); + if !loose.is_empty() { + let mut bodies = Vec::new(); + rd.node_items = rd.used_items.clone(); + for r in loose { + rd.rep_bodies(r, &identity, &[], "Import", &mut bodies, 0); + } + if !bodies.is_empty() { + tree.push(StepNode { name: "Import".into(), transform: to_array(&identity), bodies, children: Vec::new() }); + } + } + if let Some(first) = ex.duplicates.first() { + rd.warn(format!("{} entity ids are defined twice (first #{first}); the first definitions were used", ex.duplicates.len())); + } + let warnings: Vec = rd.warnings.iter().map(|(m, n)| if *n > 1 { format!("{m} (×{n})") } else { m.clone() }).collect(); + if rd.bodies.is_empty() { + let why = warnings.first().cloned().unwrap_or_else(|| "the file has no solids".into()); + return Err(KernelError::Failed(format!("STEP import: no bodies ({why})"))); + } + Ok(StepImport { bodies: rd.bodies, tree, warnings, schema: ex.schemas.join(", "), system: ex.originating_system.clone() }) +} + +#[cfg(test)] +mod tests; diff --git a/crates/kernel/src/step_in/p21.rs b/crates/kernel/src/step_in/p21.rs new file mode 100644 index 00000000..996583cb --- /dev/null +++ b/crates/kernel/src/step_in/p21.rs @@ -0,0 +1,551 @@ +//! ISO 10303-21 ("Part 21") exchange structure parser: the HEADER and the DATA section's entity +//! instances, simple (`#1=NAME(...)`) and complex (`#2=(A(...) B(...))`). Written from the +//! public standard. Input is hostile: sizes, entity counts and nesting depth are capped and every +//! error is a message, never a panic. + +use std::collections::HashMap; + +/// Largest STEP text we parse. +pub const MAX_BYTES: usize = 512 << 20; +/// Most entity instances in one file. +pub const MAX_ENTITIES: usize = 20_000_000; +/// Deepest parameter list nesting. +const MAX_DEPTH: usize = 64; + +/// A parameter value. +#[derive(Clone, Debug, PartialEq)] +pub enum Param { + Int(i64), + Real(f64), + Str(String), + /// `.NAME.` (`.T.`/`.F.` included). + Enum(String), + Ref(u64), + List(Vec), + /// `NAME(params)` (a typed value such as `LENGTH_MEASURE(1.)`). + Typed(String, Vec), + Binary(String), + /// `$` + Null, + /// `*` + Derived, +} + +impl Param { + pub fn as_f64(&self) -> Option { + match self { + Param::Real(x) => Some(*x).filter(|x| x.is_finite()), + Param::Int(i) => Some(*i as f64), + Param::Typed(_, v) if v.len() == 1 => v.first().and_then(Param::as_f64), + _ => None, + } + } + pub fn as_i64(&self) -> Option { + match self { + Param::Int(i) => Some(*i), + Param::Real(x) if x.fract() == 0.0 && x.abs() < 9e15 => Some(*x as i64), + _ => None, + } + } + pub fn as_ref_id(&self) -> Option { + match self { + Param::Ref(r) => Some(*r), + _ => None, + } + } + pub fn as_str(&self) -> Option<&str> { + match self { + Param::Str(s) => Some(s), + _ => None, + } + } + pub fn as_enum(&self) -> Option<&str> { + match self { + Param::Enum(s) => Some(s), + _ => None, + } + } + /// `.T.` → true, `.F.` → false (`.U.` and anything else → None). + pub fn as_bool(&self) -> Option { + match self.as_enum()? { + "T" => Some(true), + "F" => Some(false), + _ => None, + } + } + pub fn as_list(&self) -> Option<&[Param]> { + match self { + Param::List(v) => Some(v), + _ => None, + } + } +} + +/// One record of an entity instance: a simple instance has one, a complex one several. +#[derive(Clone, Debug, PartialEq)] +pub struct Record { + pub name: String, + pub params: Vec, +} + +/// An entity instance. +#[derive(Clone, Debug, PartialEq)] +pub struct Entity { + pub records: Vec, +} + +impl Entity { + pub fn is_complex(&self) -> bool { + self.records.len() > 1 + } + /// Name of a simple instance (the first record of a complex one). + pub fn name(&self) -> &str { + self.records.first().map(|r| r.name.as_str()).unwrap_or("") + } + /// Parameters of the record called `name`. + pub fn record(&self, name: &str) -> Option<&[Param]> { + self.records.iter().find(|r| r.name == name).map(|r| r.params.as_slice()) + } + pub fn has(&self, name: &str) -> bool { + self.records.iter().any(|r| r.name == name) + } + /// Parameters of a simple instance. + pub fn params(&self) -> &[Param] { + self.records.first().map(|r| r.params.as_slice()).unwrap_or(&[]) + } +} + +/// A parsed exchange file. +#[derive(Clone, Debug, Default)] +pub struct Exchange { + /// FILE_SCHEMA names. + pub schemas: Vec, + /// FILE_NAME originating system. + pub originating_system: String, + pub entities: HashMap, + /// Ids defined more than once (the first definition is kept). + pub duplicates: Vec, +} + +impl Exchange { + pub fn get(&self, id: u64) -> Option<&Entity> { + self.entities.get(&id) + } +} + +struct Lexer<'a> { + s: &'a [u8], + i: usize, +} + +fn err(msg: impl Into) -> Result { + Err(msg.into()) +} + +impl<'a> Lexer<'a> { + fn line(&self) -> usize { + self.s.get(..self.i).map(|b| b.iter().filter(|c| **c == b'\n').count() + 1).unwrap_or(0) + } + fn fail(&self, msg: &str) -> Result { + err(format!("STEP syntax error at line {}: {msg}", self.line())) + } + fn peek(&self) -> Option { + self.s.get(self.i).copied() + } + /// Skip whitespace and comments. + fn ws(&mut self) -> Result<(), String> { + loop { + match self.peek() { + Some(c) if c.is_ascii_whitespace() => self.i += 1, + Some(b'/') if self.s.get(self.i + 1) == Some(&b'*') => { + let rest = self.s.get(self.i + 2..).unwrap_or(&[]); + match rest.windows(2).position(|w| w == b"*/") { + Some(p) => self.i += p + 4, + None => return self.fail("unterminated comment"), + } + } + _ => return Ok(()), + } + } + } + fn eat(&mut self, c: u8) -> Result { + self.ws()?; + if self.peek() == Some(c) { + self.i += 1; + Ok(true) + } else { + Ok(false) + } + } + fn expect(&mut self, c: u8) -> Result<(), String> { + if self.eat(c)? { Ok(()) } else { self.fail(&format!("expected `{}`", c as char)) } + } + fn keyword(&mut self) -> Result { + self.ws()?; + let st = self.i; + if self.peek() == Some(b'!') { + self.i += 1; + } + while let Some(c) = self.peek() { + if c.is_ascii_alphanumeric() || c == b'_' || c == b'-' { + self.i += 1; + } else { + break; + } + } + if self.i == st { + return self.fail("expected a keyword"); + } + Ok(String::from_utf8_lossy(self.s.get(st..self.i).unwrap_or(&[])).to_ascii_uppercase()) + } + fn number(&mut self) -> Result { + let st = self.i; + while let Some(c) = self.peek() { + if c.is_ascii_digit() || matches!(c, b'+' | b'-' | b'.' | b'E' | b'e') { + self.i += 1; + } else { + break; + } + } + let t = std::str::from_utf8(self.s.get(st..self.i).unwrap_or(&[])).unwrap_or(""); + if !t.contains(['.', 'E', 'e']) + && let Ok(i) = t.parse::() + { + return Ok(Param::Int(i)); + } + // STEP allows `1.` and `1.E-3`. + let fixed = t.replace(".E", ".0E").replace(".e", ".0e"); + let fixed = if fixed.ends_with('.') { format!("{fixed}0") } else { fixed }; + match fixed.parse::() { + Ok(x) => Ok(Param::Real(x)), + Err(_) => self.fail(&format!("bad number `{t}`")), + } + } + fn string(&mut self) -> Result { + // Opening quote consumed by the caller. + let mut raw = Vec::new(); + loop { + match self.peek() { + None => return self.fail("unterminated string"), + Some(b'\'') => { + if self.s.get(self.i + 1) == Some(&b'\'') { + raw.push(b'\''); + self.i += 2; + } else { + self.i += 1; + break; + } + } + Some(b'\n') | Some(b'\r') => self.i += 1, + Some(c) => { + raw.push(c); + self.i += 1; + } + } + } + Ok(decode_string(&raw)) + } + fn param(&mut self, depth: usize) -> Result { + if depth > MAX_DEPTH { + return self.fail("parameters nested too deeply"); + } + self.ws()?; + let Some(c) = self.peek() else { return self.fail("unexpected end of file") }; + match c { + b'$' => { + self.i += 1; + Ok(Param::Null) + } + b'*' => { + self.i += 1; + Ok(Param::Derived) + } + b'\'' => { + self.i += 1; + Ok(Param::Str(self.string()?)) + } + b'"' => { + self.i += 1; + let st = self.i; + while self.peek().is_some_and(|c| c != b'"') { + self.i += 1; + } + let v = String::from_utf8_lossy(self.s.get(st..self.i).unwrap_or(&[])).to_string(); + self.expect(b'"')?; + Ok(Param::Binary(v)) + } + b'#' => { + self.i += 1; + let st = self.i; + while self.peek().is_some_and(|c| c.is_ascii_digit()) { + self.i += 1; + } + let t = std::str::from_utf8(self.s.get(st..self.i).unwrap_or(&[])).unwrap_or(""); + t.parse::().map(Param::Ref).or_else(|_| self.fail("bad entity reference")) + } + b'.' if self.s.get(self.i + 1).is_some_and(|c| c.is_ascii_alphabetic() || *c == b'_') => { + self.i += 1; + let st = self.i; + while self.peek().is_some_and(|c| c != b'.') { + self.i += 1; + } + let v = String::from_utf8_lossy(self.s.get(st..self.i).unwrap_or(&[])).trim().to_ascii_uppercase(); + self.expect(b'.')?; + Ok(Param::Enum(v)) + } + b'(' => { + self.i += 1; + Ok(Param::List(self.list_rest(depth + 1)?)) + } + c if c.is_ascii_digit() || matches!(c, b'+' | b'-' | b'.') => self.number(), + c if c.is_ascii_alphabetic() || c == b'!' => { + let name = self.keyword()?; + self.expect(b'(')?; + Ok(Param::Typed(name, self.list_rest(depth + 1)?)) + } + _ => self.fail(&format!("unexpected `{}`", c as char)), + } + } + /// Parameters after an opening `(` up to and including the `)`. + fn list_rest(&mut self, depth: usize) -> Result, String> { + let mut v = Vec::new(); + if self.eat(b')')? { + return Ok(v); + } + loop { + v.push(self.param(depth)?); + if self.eat(b',')? { + continue; + } + self.expect(b')')?; + return Ok(v); + } + } + fn record(&mut self) -> Result { + let name = self.keyword()?; + self.expect(b'(')?; + Ok(Record { name, params: self.list_rest(1)? }) + } +} + +/// Decode the Part 21 string escapes (`\X\hh`, `\X2\…\X0\`, `\X4\…\X0\`, `\S\c`, `\\`). +fn decode_string(raw: &[u8]) -> String { + let mut out = String::new(); + let mut i = 0; + let hex = |b: &[u8]| std::str::from_utf8(b).ok().and_then(|t| u32::from_str_radix(t, 16).ok()); + while let Some(&c) = raw.get(i) { + if c == b'\\' { + let rest = raw.get(i..).unwrap_or(&[]); + if rest.starts_with(b"\\\\") { + out.push('\\'); + i += 2; + continue; + } + if rest.starts_with(b"\\X\\") + && let Some(v) = rest.get(3..5).and_then(hex) + { + out.push(char::from_u32(v).unwrap_or('?')); + i += 5; + continue; + } + for (tag, width) in [(&b"\\X2\\"[..], 4usize), (&b"\\X4\\"[..], 8)] { + if rest.starts_with(tag) + && let Some(end) = rest.windows(4).position(|w| w == b"\\X0\\") + { + let body = rest.get(4..end).unwrap_or(&[]); + let units: Vec = body.chunks(width).filter_map(hex).collect(); + if width == 4 { + let u16s: Vec = units.iter().map(|u| *u as u16).collect(); + out.push_str(&String::from_utf16_lossy(&u16s)); + } else { + out.extend(units.iter().map(|u| char::from_u32(*u).unwrap_or('?'))); + } + i += end + 4; + break; + } + } + if i < raw.len() && raw.get(i) == Some(&b'\\') { + if rest.starts_with(b"\\S\\") + && let Some(&ch) = rest.get(3) + { + out.push(char::from_u32(ch as u32 + 128).unwrap_or('?')); + i += 4; + continue; + } + out.push('\\'); + i += 1; + } + continue; + } + // Plain bytes: keep UTF-8 runs intact. + let st = i; + while raw.get(i).is_some_and(|c| *c != b'\\') { + i += 1; + } + out.push_str(&String::from_utf8_lossy(raw.get(st..i).unwrap_or(&[]))); + } + out +} + +/// Parse a Part 21 file. +pub fn parse(text: &str) -> Result { + if text.len() > MAX_BYTES { + return err(format!("STEP file too large ({} MB, limit {} MB)", text.len() >> 20, MAX_BYTES >> 20)); + } + let mut lx = Lexer { s: text.as_bytes(), i: 0 }; + let magic = lx.keyword()?; + if magic != "ISO-10303-21" { + return err("not a STEP file (no ISO-10303-21 header)"); + } + lx.expect(b';')?; + let mut ex = Exchange::default(); + loop { + lx.ws()?; + if lx.peek().is_none() { + return err("STEP file ends without END-ISO-10303-21"); + } + let kw = lx.keyword()?; + match kw.as_str() { + "HEADER" => { + lx.expect(b';')?; + loop { + let k = lx.keyword()?; + if k == "ENDSEC" { + lx.expect(b';')?; + break; + } + lx.expect(b'(')?; + let params = lx.list_rest(1)?; + lx.expect(b';')?; + match k.as_str() { + "FILE_SCHEMA" => { + if let Some(Param::List(v)) = params.first() { + ex.schemas = v.iter().filter_map(|p| p.as_str().map(str::to_string)).collect(); + } + } + "FILE_NAME" => { + ex.originating_system = params.get(5).and_then(Param::as_str).unwrap_or("").to_string(); + } + _ => {} + } + } + } + "DATA" => { + // Optional `DATA('name', ('schema'));` + if lx.eat(b'(')? { + lx.list_rest(1)?; + } + lx.expect(b';')?; + loop { + lx.ws()?; + if lx.peek() != Some(b'#') { + let k = lx.keyword()?; + if k == "ENDSEC" { + lx.expect(b';')?; + break; + } + return lx.fail(&format!("unexpected `{k}` in DATA")); + } + lx.i += 1; + let st = lx.i; + while lx.peek().is_some_and(|c| c.is_ascii_digit()) { + lx.i += 1; + } + let id: u64 = match std::str::from_utf8(lx.s.get(st..lx.i).unwrap_or(&[])).ok().and_then(|t| t.parse().ok()) { + Some(id) => id, + None => return lx.fail("bad entity id"), + }; + lx.expect(b'=')?; + let records = if lx.eat(b'(')? { + let mut rs = Vec::new(); + while !lx.eat(b')')? { + if rs.len() > 64 { + return lx.fail("complex entity with too many parts"); + } + rs.push(lx.record()?); + } + rs + } else { + vec![lx.record()?] + }; + lx.expect(b';')?; + if ex.entities.len() >= MAX_ENTITIES { + return err(format!("STEP file has more than {MAX_ENTITIES} entities")); + } + // A duplicate id is invalid; some writers emit them for geometry nobody + // reads. Keep the first definition and report it. + match ex.entities.entry(id) { + std::collections::hash_map::Entry::Occupied(_) => ex.duplicates.push(id), + std::collections::hash_map::Entry::Vacant(v) => { + v.insert(Entity { records }); + } + } + } + } + "END-ISO-10303-21" => { + lx.expect(b';')?; + return Ok(ex); + } + // ANCHOR / REFERENCE / SIGNATURE sections (Part 21 edition 3): skip to ENDSEC. + "ANCHOR" | "REFERENCE" | "SIGNATURE" => { + let rest = lx.s.get(lx.i..).unwrap_or(&[]); + match rest.windows(7).position(|w| w == b"ENDSEC;") { + Some(p) => lx.i += p + 7, + None => return lx.fail("unterminated section"), + } + } + _ => return lx.fail(&format!("unexpected `{kw}`")), + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn parses_simple_and_complex() { + let t = "ISO-10303-21;\nHEADER;FILE_DESCRIPTION((''),'2;1');\nFILE_NAME('a','',(''),(''),'','sys','');\nFILE_SCHEMA(('AUTOMOTIVE_DESIGN'));ENDSEC;\nDATA;\n#1=CARTESIAN_POINT('p',(0.,1.5E1,-2));\n/* c */ #2=(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT(.MILLI.,.METRE.));\n#3=X('it''s','\\X2\\00E9\\X0\\',$,#1,LENGTH_MEASURE(1.E-3),\"0F\");\nENDSEC;\nEND-ISO-10303-21;\n"; + let ex = parse(t).unwrap(); + assert_eq!(ex.schemas, vec!["AUTOMOTIVE_DESIGN"]); + assert_eq!(ex.originating_system, "sys"); + let p = ex.get(1).unwrap(); + assert_eq!(p.name(), "CARTESIAN_POINT"); + let xyz: Vec = p.params()[1].as_list().unwrap().iter().filter_map(Param::as_f64).collect(); + assert_eq!(xyz, vec![0.0, 15.0, -2.0]); + let u = ex.get(2).unwrap(); + assert!(u.is_complex() && u.has("SI_UNIT")); + assert_eq!(u.record("SI_UNIT").unwrap()[0].as_enum(), Some("MILLI")); + let x = ex.get(3).unwrap().params(); + assert_eq!(x[0].as_str(), Some("it's")); + assert_eq!(x[1].as_str(), Some("\u{e9}")); + assert_eq!(x[2], Param::Null); + assert_eq!(x[3].as_ref_id(), Some(1)); + assert_eq!(x[4].as_f64(), Some(1e-3)); + } + + #[test] + fn hostile_text_is_an_error() { + for t in [ + "", + "garbage", + "ISO-10303-21;", + "ISO-10303-21;DATA;#1=A(", + "ISO-10303-21;DATA;#1=A('x", + "ISO-10303-21;DATA;#1=A(/*", + "ISO-10303-21;DATA;#=A();ENDSEC;END-ISO-10303-21;", + "ISO-10303-21;DATA;#99999999999999999999999=A();ENDSEC;END-ISO-10303-21;", + "ISO-10303-21;DATA;#1=A(1e999999.);ENDSEC;END-ISO-10303-21;", + ] { + let _ = parse(t); + } + let deep = format!("ISO-10303-21;DATA;#1=A({}{});ENDSEC;END-ISO-10303-21;", "(".repeat(10_000), ")".repeat(10_000)); + assert!(parse(&deep).is_err()); + let dup = parse("ISO-10303-21;DATA;#1=A(1);#1=B();ENDSEC;END-ISO-10303-21;").unwrap(); + assert_eq!(dup.duplicates, vec![1]); + assert_eq!(dup.get(1).unwrap().name(), "A"); + assert_eq!(decode_string(b"\\X\\E9"), "\u{e9}"); + assert_eq!(decode_string(b"a\\X2\\00E900E8\\X0\\b"), "a\u{e9}\u{e8}b"); + assert_eq!(decode_string(b"\\"), "\\"); + } +} diff --git a/crates/kernel/src/step_in/split.rs b/crates/kernel/src/step_in/split.rs new file mode 100644 index 00000000..985746fa --- /dev/null +++ b/crates/kernel/src/step_in/split.rs @@ -0,0 +1,409 @@ +//! Splitting band faces that wrap all the way around a closed, non-periodic parameter direction. +//! +//! A STEP face may run once around a closed surface and close its loop with a seam edge used +//! twice (a pipe segment of a torus, a revolved B-spline band). Truck meshes faces in parameter +//! space and only treats the rotation of a revolved surface as periodic, so such a face is cut +//! in two along the parameter line opposite its seam; the ring edges it crosses are cut in every +//! face that uses them. + +use std::collections::HashMap; + +use mt::{BoundedCurve, BoundedSurface, InnerSpace, Invertible, MetricSpace, ParametricCurve, ParametricSurface, SearchNearestParameter}; +use truck_modeling as mt; + +use super::geom::{Frame, R, conic_arc}; + +/// Samples per edge when looking for the split crossing. +const SAMPLES: usize = 48; + +/// Which surface parameter wraps around. +#[derive(Clone, Copy, Debug, PartialEq)] +enum Dir { + U, + V, +} + +fn uv_of(s: &mt::Surface, p: mt::Point3) -> Option<(f64, f64)> { + s.search_nearest_parameter(p, None, 100).filter(|(u, v)| u.is_finite() && v.is_finite()) +} + +fn coord(d: Dir, uv: (f64, f64)) -> f64 { + match d { + Dir::U => uv.0, + Dir::V => uv.1, + } +} + +/// The closed, non-periodic direction of a surface with its range. +fn closed_dirs(s: &mt::Surface) -> Vec<(Dir, f64, f64)> { + match s { + mt::Surface::RevolutedCurve(p) => { + let c = p.entity().entity_curve(); + let (t0, t1) = c.range_tuple(); + if c.front().distance(c.back()) < 1e-7 * (1.0 + c.front().to_homogeneous().truncate().magnitude()) { + vec![(Dir::U, t0, t1)] + } else { + vec![] + } + } + mt::Surface::BSplineSurface(_) | mt::Surface::NurbsSurface(_) => { + let ((u0, u1), (v0, v1)) = match s { + mt::Surface::BSplineSurface(b) => b.range_tuple(), + mt::Surface::NurbsSurface(n) => n.range_tuple(), + _ => return vec![], + }; + let size = (0..=4) + .flat_map(|i| (0..=4).map(move |j| (i, j))) + .map(|(i, j)| s.subs(u0 + (u1 - u0) * i as f64 / 4.0, v0 + (v1 - v0) * j as f64 / 4.0).to_homogeneous().truncate().magnitude()) + .fold(1.0, f64::max); + let tol = size * 1e-7; + let mut out = Vec::new(); + if (0..=8).all(|k| { + let v = v0 + (v1 - v0) * k as f64 / 8.0; + s.subs(u0, v).distance(s.subs(u1, v)) < tol + }) { + out.push((Dir::U, u0, u1)); + } + if (0..=8).all(|k| { + let u = u0 + (u1 - u0) * k as f64 / 8.0; + s.subs(u, v0).distance(s.subs(u, v1)) < tol + }) { + out.push((Dir::V, v0, v1)); + } + out + } + mt::Surface::Plane(_) => vec![], + } +} + +fn samples(c: &mt::Curve, n: usize) -> Vec<(f64, mt::Point3)> { + let (t0, t1) = c.range_tuple(); + (0..=n).map(|i| t0 + (t1 - t0) * i as f64 / n as f64).map(|t| (t, c.subs(t))).collect() +} + +/// Where a ring chain crosses `w = mid`: (edge index in the chain, parameter on the edge's +/// absolute curve, point). +fn crossing(s: &mt::Surface, d: Dir, chain: &[mt::Edge], mid: f64, period: f64) -> R<(usize, f64, mt::Point3)> { + for (k, e) in chain.iter().enumerate() { + let c = e.curve(); + let pts = samples(&c, SAMPLES); + let ws: Vec> = pts.iter().map(|(_, p)| uv_of(s, *p).map(|uv| coord(d, uv))).collect(); + for i in 0..SAMPLES { + let (Some(Some(wa)), Some(Some(wb)), Some((ta, _)), Some((tb, _))) = (ws.get(i), ws.get(i + 1), pts.get(i), pts.get(i + 1)) else { + continue; + }; + if (wb - wa).abs() > period * 0.5 || (wa - mid) * (wb - mid) > 0.0 { + continue; + } + // Bisection on the curve parameter. + let (mut lo, mut hi, mut flo) = (*ta, *tb, wa - mid); + for _ in 0..60 { + let m = 0.5 * (lo + hi); + let Some(w) = uv_of(s, c.subs(m)).map(|uv| coord(d, uv)) else { break }; + let f = w - mid; + if (f - flo).abs() > period * 0.5 { + break; + } + if f * flo <= 0.0 { + hi = m; + } else { + lo = m; + flo = f; + } + } + let t = 0.5 * (lo + hi); + return Ok((k, t, c.subs(t))); + } + } + Err("no crossing with the split line".into()) +} + +/// Iso-parameter edge curve `w = mid` from `q1` to `q2`. +fn iso_curve(s: &mt::Surface, d: Dir, mid: f64, q1: mt::Point3, q2: mt::Point3, seam: &[mt::Point3]) -> R { + match s { + mt::Surface::RevolutedCurve(p) => { + // A circle around the axis, turning opposite to the seam (which runs from the q2 + // side to the q1 side). + let (o, a) = (p.entity().origin(), p.entity().axis()); + let az = |x: mt::Point3, f: &Frame| { + let v = x - f.o; + v.dot(f.y).atan2(v.dot(f.x)) + }; + let c = o + a * (q1 - o).dot(a); + let r = (q1 - c).magnitude(); + if !(r > 1e-9) { + return Err("split line on the axis".into()); + } + let x = (q1 - c) / r; + let f = Frame { o: c, x, y: a.cross(x), z: a }; + let mut seam_turn = 0.0; + for w in seam.windows(2) { + if let (Some(p0), Some(p1)) = (w.first(), w.get(1)) { + let d = az(*p1, &f) - az(*p0, &f); + seam_turn += (d + std::f64::consts::PI).rem_euclid(std::f64::consts::TAU) - std::f64::consts::PI; + } + } + let sweep = -seam_turn; + if sweep.abs() < 1e-9 { + return Err("degenerate split line".into()); + } + let f = if sweep > 0.0 { f } else { Frame { o: c, x, y: -a.cross(x), z: -a } }; + let arc = conic_arc(&f, r, r, 0.0, sweep.abs()); + if arc.back().distance(q2) > 1e-6 * (1.0 + r) { + return Err("split line does not reach the other ring".into()); + } + Ok(mt::Curve::NurbsCurve(arc)) + } + mt::Surface::BSplineSurface(_) | mt::Surface::NurbsSurface(_) => { + let hom: mt::BSplineSurface = match s { + mt::Surface::BSplineSurface(b) => mt::BSplineSurface::new( + b.knot_vecs().clone(), + b.control_points().iter().map(|row| row.iter().map(|p| p.to_homogeneous()).collect()).collect(), + ), + mt::Surface::NurbsSurface(n) => n.non_rationalized().clone(), + _ => return Err("unsupported surface".into()), + }; + let iso = mt::Curve::NurbsCurve(iso_bsp(&hom, mid, d == Dir::V)?); + let near = |p: mt::Point3| { + let (t0, t1) = iso.range_tuple(); + let mut best = (t0, f64::INFINITY); + for i in 0..=256 { + let t = t0 + (t1 - t0) * i as f64 / 256.0; + let dd = iso.subs(t).distance(p); + if dd < best.1 { + best = (t, dd); + } + } + if let Some(t) = iso.search_nearest_parameter(p, best.0, 64) { + let t = t.clamp(t0, t1); + if iso.subs(t).distance(p) <= best.1 { + best = (t, iso.subs(t).distance(p)); + } + } + best + }; + let ((ta, da), (tb, db)) = (near(q1), near(q2)); + let scale = 1.0 + q1.to_homogeneous().truncate().magnitude(); + if da > 1e-5 * scale || db > 1e-5 * scale || (ta - tb).abs() < 1e-12 { + return Err("split line does not reach the rings".into()); + } + let (lo, hi) = if ta < tb { (ta, tb) } else { (tb, ta) }; + let (t0, t1) = iso.range_tuple(); + let mut c = iso; + let eps = (t1 - t0) * 1e-9; + if lo > t0 + eps { + c = mt::Cut::cut(&mut c, lo); + } + if hi < t1 - eps { + let _ = mt::Cut::cut(&mut c, hi); + } + if ta > tb { + c.invert(); + } + Ok(c) + } + mt::Surface::Plane(_) => Err("planes do not wrap".into()), + } +} + +/// The iso-curve of a homogeneous B-spline surface at `v = w` (`along_u`) or `u = w`. +fn iso_bsp(b: &mt::BSplineSurface, w: f64, along_u: bool) -> R> { + let (ku, kv) = b.knot_vecs(); + let ctrl = b.control_points(); + let nu = ctrl.len(); + let nv = ctrl.first().map(Vec::len).unwrap_or(0); + let zero = mt::Vector4::new(0.0, 0.0, 0.0, 0.0); + let pts: Vec = if along_u { + let deg = kv.len().checked_sub(nv + 1).ok_or("bad knot vector")?; + let basis = kv.try_bspline_basis_functions(deg, w).map_err(|e| e.to_string())?; + ctrl.iter().map(|row| row.iter().zip(&basis).fold(zero, |acc, (p, f)| acc + *p * *f)).collect() + } else { + let deg = ku.len().checked_sub(nu + 1).ok_or("bad knot vector")?; + let basis = ku.try_bspline_basis_functions(deg, w).map_err(|e| e.to_string())?; + (0..nv).map(|j| ctrl.iter().zip(&basis).fold(zero, |acc, (row, f)| acc + row.get(j).copied().unwrap_or(zero) * *f)).collect() + }; + let knots = if along_u { ku.clone() } else { kv.clone() }; + mt::BSplineCurve::try_new(knots, pts).map(mt::NurbsCurve::new).map_err(|e| e.to_string()) +} + +type Cut = (mt::Edge, (mt::Edge, mt::Edge)); + +/// Split a chain at parameter `t` of its edge `k` (absolute curve): the edges before and after +/// the split point, the edge cut (unless the point is an existing vertex) and the split vertex. +fn split_chain(chain: &[mt::Edge], k: usize, t: f64, q: mt::Point3) -> R<(Vec, Vec, Option, mt::Vertex)> { + let e = chain.get(k).ok_or("bad crossing")?; + let (t0, t1) = e.curve().range_tuple(); + let eps = (t1 - t0).abs() * 1e-7; + let at_front = (t - t0).abs() <= eps; + let at_back = (t - t1).abs() <= eps; + if at_front || at_back { + // The split point is a vertex already: split the chain between edges. + let v = if at_front { e.absolute_front().clone() } else { e.absolute_back().clone() }; + let first_end = at_front != e.orientation(); // the point is this edge's end in loop order + let cut_at = if first_end { k + 1 } else { k }; + let before = chain.iter().take(cut_at).cloned().collect(); + let after = chain.iter().skip(cut_at).cloned().collect(); + return Ok((before, after, None, v)); + } + let v = mt::Vertex::new(q); + let (a, b) = e.cut_with_parameter(&v, t).ok_or("cannot cut a ring edge at the split line")?; + let mut before: Vec = chain.iter().take(k).cloned().collect(); + before.push(a.clone()); + let mut after = vec![b.clone()]; + after.extend(chain.iter().skip(k + 1).cloned()); + Ok((before, after, Some((e.clone(), (a, b))), v)) +} + +struct Plan { + /// Edge cuts: (edge as used in the band's loop, its two pieces in that loop's order). + cuts: Vec, + faces: (mt::Face, mt::Face), +} + +fn plan(face: &mt::Face) -> R> { + let s = face.surface(); + let wires = face.absolute_boundaries().clone(); + let mut uses: HashMap = HashMap::new(); + for e in wires.iter().flat_map(|w| w.edge_iter()) { + *uses.entry(e.id()).or_default() += 1; + } + let is_seam = |e: &mt::Edge| uses.get(&e.id()).copied().unwrap_or(0) > 1; + if !wires.iter().flat_map(|w| w.edge_iter()).any(is_seam) { + return Ok(None); + } + // The closed direction whose boundary carries the seam. + let seam_pts: Vec = + wires.iter().flat_map(|w| w.edge_iter()).filter(|e| is_seam(e)).flat_map(|e| samples(&e.curve(), 4).into_iter().map(|x| x.1)).collect(); + let Some((d, w0, w1)) = closed_dirs(&s).into_iter().find(|(d, w0, w1)| { + let p = w1 - w0; + seam_pts.iter().all(|q| uv_of(&s, *q).map(|uv| coord(*d, uv)).is_some_and(|w| (w - w0).abs().min((w - w1).abs()) < p * 1e-4)) + }) else { + return Ok(None); + }; + let (period, mid) = (w1 - w0, 0.5 * (w0 + w1)); + let outer_idx = wires.iter().position(|w| w.edge_iter().any(is_seam)).ok_or("no seam loop")?; + if wires.iter().enumerate().any(|(i, w)| i != outer_idx && w.edge_iter().any(is_seam)) { + return Err("seam edges in several loops".into()); + } + let outer: Vec = wires.get(outer_idx).ok_or("no seam loop")?.edge_iter().cloned().collect(); + let n = outer.len(); + let start = (0..n) + .find(|&i| outer.get(i).is_some_and(is_seam) && outer.get((i + n - 1) % n).is_some_and(|e| !is_seam(e))) + .ok_or("loop made only of seam edges")?; + let rot: Vec = (0..n).filter_map(|i| outer.get((start + i) % n).cloned()).collect(); + let mut chains: Vec<(bool, Vec)> = Vec::new(); + for e in rot { + let seam = is_seam(&e); + match chains.last_mut() { + Some((s0, c)) if *s0 == seam => c.push(e), + _ => chains.push((seam, vec![e])), + } + } + let [(true, s1), (false, r1), (true, s2), (false, r2)] = chains.as_slice() else { + return Err(format!("unsupported seam loop layout ({} runs)", chains.len())); + }; + let s2_ids: Vec = s2.iter().rev().map(mt::Edge::id).collect(); + if s1.iter().map(mt::Edge::id).collect::>() != s2_ids { + return Err("seam edges are not used back and forth".into()); + } + let (k1, t1, q1) = crossing(&s, d, r1, mid, period)?; + let (k2, t2, q2) = crossing(&s, d, r2, mid, period)?; + let seam_line: Vec = s1.iter().flat_map(|e| samples(&e.oriented_curve(), 16).into_iter().map(|x| x.1)).collect(); + let iso = iso_curve(&s, d, mid, q1, q2, &seam_line)?; + let (r1a, r1b, c1, v1) = split_chain(r1, k1, t1, q1)?; + let (r2a, r2b, c2, v2) = split_chain(r2, k2, t2, q2)?; + let e_split = mt::Edge::new_unchecked(&v1, &v2, iso); + // Face A: split, rest of ring 2, seam, start of ring 1. Face B: the other half. + let mut wa: Vec = vec![e_split.clone()]; + wa.extend(r2b); + wa.extend(s1.iter().cloned()); + wa.extend(r1a.iter().cloned()); + let mut wb: Vec = vec![e_split.inverse()]; + wb.extend(r1b); + wb.extend(s2.iter().cloned()); + wb.extend(r2a); + let wire_a: mt::Wire = wa.into(); + let wire_b: mt::Wire = wb.into(); + if !wire_a.is_closed() || !wire_b.is_closed() { + return Err("split loops do not close".into()); + } + // Holes go to the half that holds them. + let side_a = r1a + .iter() + .flat_map(|e| samples(&e.oriented_curve(), 4)) + .filter_map(|(_, p)| uv_of(&s, p).map(|uv| coord(d, uv))) + .find(|w| (w - mid).abs() > period * 1e-3) + .map(|w| w < mid) + .ok_or("cannot place the split halves")?; + let (mut ha, mut hb) = (vec![wire_a], vec![wire_b]); + for (i, w) in wires.iter().enumerate() { + if i == outer_idx { + continue; + } + let ws: Vec = w.edge_iter().flat_map(|e| samples(&e.curve(), 4)).filter_map(|(_, p)| uv_of(&s, p).map(|uv| coord(d, uv))).collect(); + if ws.iter().all(|x| (*x < mid) == side_a) { + ha.push(w.clone()); + } else if ws.iter().all(|x| (*x < mid) != side_a) { + hb.push(w.clone()); + } else { + return Err("a hole crosses the split line".into()); + } + } + let mut fa = mt::Face::new_unchecked(ha, s.clone()); + let mut fb = mt::Face::new_unchecked(hb, s); + if !face.orientation() { + fa.invert(); + fb.invert(); + } + Ok(Some(Plan { cuts: c1.into_iter().chain(c2).collect(), faces: (fa, fb) })) +} + +/// Replace cut edges in a face's loops. +fn recut(face: &mt::Face, cuts: &[(mt::Edge, (mt::Edge, mt::Edge))]) -> mt::Face { + let touches = face.absolute_boundaries().iter().flat_map(|w| w.edge_iter()).any(|e| cuts.iter().any(|(c, _)| c.id() == e.id())); + if !touches { + return face.clone(); + } + let wires: Vec = face + .absolute_boundaries() + .iter() + .map(|w| { + let mut out: Vec = Vec::new(); + for e in w.edge_iter() { + match cuts.iter().find(|(c, _)| c.id() == e.id()) { + Some((c, (a, b))) if c.orientation() == e.orientation() => out.extend([a.clone(), b.clone()]), + Some((_, (a, b))) => out.extend([b.inverse(), a.inverse()]), + None => out.push(e.clone()), + } + } + out.into() + }) + .collect(); + let mut f = mt::Face::new_unchecked(wires, face.surface()); + if !face.orientation() { + f.invert(); + } + f +} + +/// Split every band face of a shell (see the module docs). Faces that cannot be split stay as +/// they are and are reported. +pub(crate) fn split_bands(mut faces: Vec, warnings: &mut Vec) -> Vec { + let n = faces.len(); + for i in 0..n { + let Some(face) = faces.get(i).cloned() else { continue }; + match crate::guard("step band split", || plan(&face).map_err(crate::KernelError::Failed)) { + Ok(Some(p)) => { + faces = faces.iter().map(|f| recut(f, &p.cuts)).collect(); + if let Some(slot) = faces.get_mut(i) { + *slot = p.faces.0; + } + faces.push(p.faces.1); + } + Ok(None) => {} + Err(e) => warnings + .push(format!("face wrapping around a closed surface kept whole: {}", e.to_string().trim_start_matches("the operation failed: "))), + } + } + faces +} diff --git a/crates/kernel/src/step_in/tests.rs b/crates/kernel/src/step_in/tests.rs new file mode 100644 index 00000000..9387104d --- /dev/null +++ b/crates/kernel/src/step_in/tests.rs @@ -0,0 +1,285 @@ +use solvecraft_geom::{Loop2, Plane, Region2, Vec2, Vec3}; + +use super::*; +use crate::{BoolOp, boolean, box_solid, cylinder, fillet, measure, revolve, sphere, step_export, torus}; + +fn rel(a: f64, b: f64) -> f64 { + (a - b).abs() / b.abs().max(1e-9) +} + +/// Export → import → same volume, area and face count. +fn round_trip(b: &Body) -> StepImport { + let text = step_export(&[b], "test").unwrap(); + let imp = step_import(&text).unwrap(); + assert_eq!(imp.bodies.len(), 1, "{:?}", imp.warnings); + // truck's writer repeats entity ids for the surfaces of intersection curves (unused here). + assert!(imp.warnings.iter().all(|w| w.contains("defined twice")), "{:?}", imp.warnings); + let (m0, m1) = (measure(b).unwrap(), measure(&imp.bodies[0].body).unwrap()); + assert!(rel(m1.volume, m0.volume) < 1e-3, "volume {} vs {}", m1.volume, m0.volume); + assert!(rel(m1.area, m0.area) < 1e-3, "area {} vs {}", m1.area, m0.area); + assert_eq!(imp.bodies[0].file_faces, b.face_count()); + assert!(imp.bodies[0].closed); + imp +} + +#[test] +fn round_trip_primitives() { + round_trip(&box_solid(Vec3::ZERO, Vec3::new(10.0, 20.0, 30.0)).unwrap()); + round_trip(&cylinder(Vec3::new(1.0, 2.0, 3.0), Vec3::new(0.0, 1.0, 1.0), 5.0, 12.0).unwrap()); + round_trip(&sphere(Vec3::new(1.0, 2.0, 3.0), 10.0).unwrap()); + round_trip(&torus(Vec3::ZERO, 20.0, 4.0).unwrap()); +} + +#[test] +fn round_trip_features() { + let plate = box_solid(Vec3::ZERO, Vec3::new(40.0, 30.0, 10.0)).unwrap(); + let pin = cylinder(Vec3::new(20.0, 15.0, -1.0), Vec3::Z, 5.0, 12.0).unwrap(); + let holed = boolean(&plate, &pin, BoolOp::Cut).unwrap().unwrap(); + round_trip(&holed); + let edge = Vec3::new(20.0, 0.0, 10.0); + round_trip(&fillet(&plate, &[edge], 3.0).unwrap()); + let ring = Region2 { outer: Loop2::circle(Vec2::new(10.0, 0.0), 3.0), holes: vec![] }; + let donut = revolve(&Plane::XZ, &[ring], Vec2::ZERO, Vec2::Y, 1.5).unwrap(); + round_trip(&donut[0]); +} + +#[test] +fn several_bodies_and_names() { + let a = box_solid(Vec3::ZERO, Vec3::new(10.0, 10.0, 10.0)).unwrap(); + let b = box_solid(Vec3::new(20.0, 0.0, 0.0), Vec3::new(30.0, 10.0, 10.0)).unwrap(); + let text = step_export(&[&a, &b], "test").unwrap(); + let imp = step_import(&text).unwrap(); + assert_eq!(imp.bodies.len(), 2); + assert!(imp.bodies.iter().all(|b| !b.name.is_empty())); + let vol: f64 = imp.bodies.iter().map(|b| measure(&b.body).unwrap().volume).sum(); + assert!(rel(vol, 2000.0) < 1e-6); +} + +#[test] +fn imported_bodies_take_features() { + let plate = box_solid(Vec3::ZERO, Vec3::new(40.0, 30.0, 10.0)).unwrap(); + let imp = step_import(&step_export(&[&plate], "t").unwrap()).unwrap(); + let body = &imp.bodies[0].body; + // Fillet an edge, cut a hole through it. + let rounded = fillet(body, &[Vec3::new(20.0, 0.0, 10.0)], 2.0).unwrap(); + let v_fillet = 40.0 * (4.0 - std::f64::consts::PI); + assert!(rel(measure(&rounded).unwrap().volume, 12000.0 - v_fillet) < 1e-3); + let pin = cylinder(Vec3::new(20.0, 15.0, -1.0), Vec3::Z, 5.0, 12.0).unwrap(); + let holed = boolean(&rounded, &pin, BoolOp::Cut).unwrap().unwrap(); + let want = 12000.0 - v_fillet - std::f64::consts::PI * 25.0 * 10.0; + assert!(rel(measure(&holed).unwrap().volume, want) < 1e-3); +} + +const HEAD: &str = "ISO-10303-21;\nHEADER;\nFILE_DESCRIPTION((''),'2;1');\nFILE_NAME('t','',(''),(''),'','hand','');\nFILE_SCHEMA(('AUTOMOTIVE_DESIGN { 1 0 10303 214 1 1 1 1 }'));\nENDSEC;\nDATA;\n"; +const TAIL: &str = "ENDSEC;\nEND-ISO-10303-21;\n"; + +/// A hand-written unit cube brep (`#100`…) in a representation with the given unit and items. +fn cube_entities() -> String { + let mut s = String::new(); + let pts = [(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0), (0, 0, 1), (1, 0, 1), (1, 1, 1), (0, 1, 1)]; + for (i, (x, y, z)) in pts.iter().enumerate() { + s += &format!("#{}=CARTESIAN_POINT('',({x}.,{y}.,{z}.));\n#{}=VERTEX_POINT('',#{});\n", 200 + i, 210 + i, 200 + i); + } + // 12 edges as lines (direction vector filler; trimmed by vertices). + let edges = [(0, 1), (1, 2), (2, 3), (3, 0), (4, 5), (5, 6), (6, 7), (7, 4), (0, 4), (1, 5), (2, 6), (3, 7)]; + s += "#290=DIRECTION('',(1.,0.,0.));\n#291=VECTOR('',#290,1.);\n"; + for (k, (a, b)) in edges.iter().enumerate() { + s += &format!("#{}=LINE('',#{},#291);\n#{}=EDGE_CURVE('',#{},#{},#{},.T.);\n", 300 + k, 200 + a, 320 + k, 210 + a, 210 + b, 300 + k); + } + // Faces: (edge, forward?) loops counter-clockwise seen from outside, and plane frames. + let faces: [(&[(usize, bool)], (i32, i32, i32), (i32, i32, i32)); 6] = [ + (&[(3, false), (2, false), (1, false), (0, false)], (0, 0, 0), (0, 0, -1)), + (&[(4, true), (5, true), (6, true), (7, true)], (0, 0, 1), (0, 0, 1)), + (&[(0, true), (9, true), (4, false), (8, false)], (0, 0, 0), (0, -1, 0)), + (&[(1, true), (10, true), (5, false), (9, false)], (1, 0, 0), (1, 0, 0)), + (&[(2, true), (11, true), (6, false), (10, false)], (1, 1, 0), (0, 1, 0)), + (&[(3, true), (8, true), (7, false), (11, false)], (0, 1, 0), (-1, 0, 0)), + ]; + let mut face_ids = Vec::new(); + for (f, (lp, o, n)) in faces.iter().enumerate() { + let base = 400 + f * 20; + let oes: Vec = lp + .iter() + .enumerate() + .map(|(j, (e, fwd))| { + s += &format!("#{}=ORIENTED_EDGE('',*,*,#{},{});\n", base + j, 320 + e, if *fwd { ".T." } else { ".F." }); + format!("#{}", base + j) + }) + .collect(); + s += &format!("#{}=EDGE_LOOP('',({}));\n#{}=FACE_OUTER_BOUND('',#{},.T.);\n", base + 5, oes.join(","), base + 6, base + 5); + s += &format!("#{}=CARTESIAN_POINT('',({}.,{}.,{}.));\n#{}=DIRECTION('',({}.,{}.,{}.));\n", base + 7, o.0, o.1, o.2, base + 8, n.0, n.1, n.2); + s += &format!("#{}=AXIS2_PLACEMENT_3D('',#{},#{},$);\n#{}=PLANE('',#{});\n", base + 9, base + 7, base + 8, base + 10, base + 9); + s += &format!("#{}=ADVANCED_FACE('',(#{}),#{},.T.);\n", base + 11, base + 6, base + 10); + face_ids.push(format!("#{}", base + 11)); + } + s += &format!("#100=CLOSED_SHELL('',({}));\n#101=MANIFOLD_SOLID_BREP('Cube',#100);\n", face_ids.join(",")); + s +} + +fn context(id: u64, unit: &str) -> String { + format!( + "#{id}=(GEOMETRIC_REPRESENTATION_CONTEXT(3)GLOBAL_UNCERTAINTY_ASSIGNED_CONTEXT((#{}))GLOBAL_UNIT_ASSIGNED_CONTEXT((#{},#{}))REPRESENTATION_CONTEXT('',''));\n#{}=UNCERTAINTY_MEASURE_WITH_UNIT(LENGTH_MEASURE(1.E-6),#{},'','');\n#{}={unit};\n#{}=(NAMED_UNIT(*)PLANE_ANGLE_UNIT()SI_UNIT($,.RADIAN.));\n", + id + 1, + id + 2, + id + 3, + id + 1, + id + 2, + id + 2, + id + 3 + ) +} + +fn product(pd: u64, name: &str) -> String { + format!( + "#{}=PRODUCT('{name}','{name}','',(#{}));\n#{}=PRODUCT_CONTEXT('',#9,'mechanical');\n#{}=PRODUCT_DEFINITION_FORMATION('','',#{});\n#{pd}=PRODUCT_DEFINITION('design','',#{},#8);\n#{}=PRODUCT_DEFINITION_SHAPE('','',#{pd});\n", + pd + 1, + pd + 2, + pd + 2, + pd + 3, + pd + 1, + pd + 3, + pd + 4 + ) +} + +#[test] +fn hand_written_part_in_inches_with_colour() { + let mut t = String::from(HEAD); + t += "#8=PRODUCT_DEFINITION_CONTEXT('part definition',#9,'design');\n#9=APPLICATION_CONTEXT('');\n"; + t += &cube_entities(); + t += &context( + 20, + "(CONVERSION_BASED_UNIT('INCH',#30)LENGTH_UNIT()NAMED_UNIT(*));\n#30=LENGTH_MEASURE_WITH_UNIT(LENGTH_MEASURE(25.4),#31);\n#31=(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT(.MILLI.,.METRE.))", + ); + t += &product(50, "Widget"); + t += "#60=ADVANCED_BREP_SHAPE_REPRESENTATION('',(#101),#20);\n#61=SHAPE_DEFINITION_REPRESENTATION(#54,#62);\n#62=SHAPE_REPRESENTATION('',(),#20);\n#63=SHAPE_REPRESENTATION_RELATIONSHIP('','',#62,#60);\n"; + t += "#70=STYLED_ITEM('',(#71),#101);\n#71=PRESENTATION_STYLE_ASSIGNMENT((#72));\n#72=SURFACE_STYLE_USAGE(.BOTH.,#73);\n#73=SURFACE_SIDE_STYLE('',(#74));\n#74=SURFACE_STYLE_FILL_AREA(#75);\n#75=FILL_AREA_STYLE('',(#76));\n#76=FILL_AREA_STYLE_COLOUR('',#77);\n#77=COLOUR_RGB('',1.,0.5,0.);\n"; + t += TAIL; + let imp = step_import(&t).unwrap(); + assert!(imp.warnings.is_empty(), "{:?}", imp.warnings); + assert_eq!(imp.bodies.len(), 1); + let b = &imp.bodies[0]; + assert_eq!(b.name, "Cube"); + assert_eq!(b.path, vec!["Widget".to_string()]); + assert_eq!(b.color, Some([1.0, 0.5, 0.0])); + let m = measure(&b.body).unwrap(); + assert!(rel(m.volume, 25.4f64.powi(3)) < 1e-9, "{}", m.volume); + assert_eq!(m.faces, 6); + assert_eq!(imp.tree.len(), 1); + assert_eq!(imp.tree[0].name, "Widget"); +} + +#[test] +fn hand_written_assembly() { + let mut t = String::from(HEAD); + t += "#8=PRODUCT_DEFINITION_CONTEXT('part definition',#9,'design');\n#9=APPLICATION_CONTEXT('');\n"; + t += &cube_entities(); + t += &context(20, "(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT(.CENTI.,.METRE.))"); + t += &context(30, "(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT(.MILLI.,.METRE.))"); + // Part (cube, in cm) and assembly (mm) with two placed instances of it. + t += &product(50, "Block"); + t += &product(70, "Assembly"); + t += "#60=ADVANCED_BREP_SHAPE_REPRESENTATION('',(#101,#65),#20);\n#61=SHAPE_DEFINITION_REPRESENTATION(#54,#60);\n"; + t += "#65=AXIS2_PLACEMENT_3D('',#66,$,$);\n#66=CARTESIAN_POINT('',(0.,0.,0.));\n"; + t += "#80=SHAPE_REPRESENTATION('',(#81,#84),#30);\n#81=AXIS2_PLACEMENT_3D('',#82,$,$);\n#82=CARTESIAN_POINT('',(0.,0.,0.));\n"; + t += "#84=AXIS2_PLACEMENT_3D('',#85,#86,#87);\n#85=CARTESIAN_POINT('',(100.,0.,0.));\n#86=DIRECTION('',(0.,0.,1.));\n#87=DIRECTION('',(0.,1.,0.));\n"; + t += "#88=SHAPE_DEFINITION_REPRESENTATION(#74,#80);\n"; + for (k, item) in [(0u64, 81u64), (1, 84)] { + let n = 600 + k * 10; + t += &format!("#{n}=NEXT_ASSEMBLY_USAGE_OCCURRENCE('{k}','Block:{k}','',#70,#50,$);\n"); + t += &format!("#{}=PRODUCT_DEFINITION_SHAPE('','',#{n});\n", n + 1); + t += &format!("#{}=ITEM_DEFINED_TRANSFORMATION('','',#65,#{item});\n", n + 2); + t += &format!( + "#{}=(REPRESENTATION_RELATIONSHIP('','',#60,#80)REPRESENTATION_RELATIONSHIP_WITH_TRANSFORMATION(#{})SHAPE_REPRESENTATION_RELATIONSHIP());\n", + n + 3, + n + 2 + ); + t += &format!("#{}=CONTEXT_DEPENDENT_SHAPE_REPRESENTATION(#{},#{});\n", n + 4, n + 3, n + 1); + } + t += TAIL; + let imp = step_import(&t).unwrap(); + assert!(imp.warnings.is_empty(), "{:?}", imp.warnings); + assert_eq!(imp.bodies.len(), 2); + let boxes: Vec<_> = imp.bodies.iter().map(|b| measure(&b.body).unwrap()).collect(); + for m in &boxes { + assert!(rel(m.volume, 1000.0) < 1e-9, "{}", m.volume); + } + // First at the origin, second moved 100 mm along X and turned 90° about Z. + assert!(boxes[0].bbox.min.dist(Vec3::ZERO) < 1e-6 && boxes[0].bbox.max.dist(Vec3::new(10.0, 10.0, 10.0)) < 1e-6); + assert!(boxes[1].bbox.min.dist(Vec3::new(90.0, 0.0, 0.0)) < 1e-6, "{:?}", boxes[1].bbox); + assert!(boxes[1].bbox.max.dist(Vec3::new(100.0, 10.0, 10.0)) < 1e-6, "{:?}", boxes[1].bbox); + assert_eq!(imp.tree.len(), 1); + let root = &imp.tree[0]; + assert_eq!(root.name, "Assembly"); + assert_eq!(root.children.len(), 2); + assert!(root.children.iter().all(|c| c.name == "Block" && c.bodies.len() == 1)); + assert_eq!(imp.bodies[1].path, vec!["Assembly".to_string(), "Block".to_string()]); + assert!((root.children[1].transform[3][0] - 100.0).abs() < 1e-9); +} + +#[test] +fn unsupported_face_is_left_out_and_reported() { + let mut t = String::from(HEAD); + let cube = cube_entities().replace("#410=PLANE('',#409);", "#410=OFFSET_SURFACE('',#409,1.,.F.);"); + t += &cube; + t += "#60=ADVANCED_BREP_SHAPE_REPRESENTATION('',(#101),#20);\n"; + t += &context(20, "(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT(.MILLI.,.METRE.))"); + t += TAIL; + let imp = step_import(&t).unwrap(); + assert_eq!(imp.bodies.len(), 1); + assert!(!imp.bodies[0].closed); + assert_eq!(imp.bodies[0].body.face_count(), 5); + assert!(imp.warnings.iter().any(|w| w.contains("OFFSET_SURFACE")), "{:?}", imp.warnings); + assert!(imp.warnings.iter().any(|w| w.contains("open body")), "{:?}", imp.warnings); +} + +#[test] +fn hostile_files_never_panic() { + let good = step_export(&[&box_solid(Vec3::ZERO, Vec3::new(5.0, 5.0, 5.0)).unwrap()], "t").unwrap(); + let mut cases: Vec = vec![ + String::new(), + "ISO-10303-21;".into(), + format!("{HEAD}{TAIL}"), + format!("{HEAD}#1=MANIFOLD_SOLID_BREP('',#1);\n#2=ADVANCED_BREP_SHAPE_REPRESENTATION('',(#1),#3);\n{TAIL}"), + format!( + "{HEAD}#1=MANIFOLD_SOLID_BREP('',#2);\n#2=CLOSED_SHELL('',(#3));\n#3=ADVANCED_FACE('',(#4),#5,.T.);\n#4=FACE_BOUND('',#6,.T.);\n#6=EDGE_LOOP('',(#7));\n#7=ORIENTED_EDGE('',*,*,#7,.T.);\n#5=PLANE('',#5);\n#9=SHAPE_REPRESENTATION('',(#1),#10);\n{TAIL}" + ), + format!( + "{HEAD}#1=CARTESIAN_POINT('',(1.E308,1.E308,1.E308,1.,2.));\n#2=MANIFOLD_SOLID_BREP('',#1);\n#3=SHAPE_REPRESENTATION('',(#2),$);\n{TAIL}" + ), + format!( + "{HEAD}#1=NEXT_ASSEMBLY_USAGE_OCCURRENCE('','','',#2,#2,$);\n#2=PRODUCT_DEFINITION('','',#3,#3);\n#3=PRODUCT_DEFINITION_FORMATION('','',#3);\n{TAIL}" + ), + ]; + // Truncations and byte flips of a real file. + for cut in (0..good.len()).step_by(97) { + cases.push(good.get(..cut).unwrap_or("").to_string()); + } + for (i, from, to) in [(0, "1.0", "-1.0"), (1, "#1", "#9999"), (2, ".T.", ".F."), (3, "(0.0", "(nan"), (4, "B_SPLINE", "X_SPLINE")] { + let _ = i; + cases.push(good.replacen(from, to, 3)); + } + cases.push(good.replace("CLOSED_SHELL", "OPEN_SHELL")); + for c in &cases { + let r = std::panic::catch_unwind(|| step_import(c)); + assert!(r.is_ok(), "panicked on {:?}", c.get(..c.len().min(200))); + } +} + +#[test] +fn assembly_cycles_are_bounded() { + let mut t = String::from(HEAD); + t += "#8=PRODUCT_DEFINITION_CONTEXT('part definition',#9,'design');\n#9=APPLICATION_CONTEXT('');\n"; + t += &cube_entities(); + t += &context(20, "(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT(.MILLI.,.METRE.))"); + t += &product(50, "A"); + t += &product(70, "B"); + t += "#60=ADVANCED_BREP_SHAPE_REPRESENTATION('',(#101),#20);\n#61=SHAPE_DEFINITION_REPRESENTATION(#54,#60);\n"; + // A contains B contains A. + t += "#90=NEXT_ASSEMBLY_USAGE_OCCURRENCE('1','','',#50,#70,$);\n#91=NEXT_ASSEMBLY_USAGE_OCCURRENCE('2','','',#70,#50,$);\n"; + t += TAIL; + let r = step_import(&t); + // Every product is somebody's component: the solid is still imported (outside the tree). + let imp = r.unwrap(); + assert_eq!(imp.bodies.len(), 1); +}