diff --git a/crates/step/src/topology.rs b/crates/step/src/topology.rs index d6ac9dd71c..9328a55499 100644 --- a/crates/step/src/topology.rs +++ b/crates/step/src/topology.rs @@ -44,12 +44,13 @@ use solidcraft_geom::brep::{ Surface as BSurface, Vertex, }; -use crate::geom::{BSplineCurve, BSplineSurface, Curve, P3, Placement, Surface, dot, scale, sub}; +use crate::geom::{BSplineCurve, BSplineSurface, Curve, P3, Placement, dot, scale, sub}; use crate::p21::{Entity, Exchange, Param}; use crate::read::{R, Reader}; use crate::units; mod seams; +mod swept; /// Tolerance (mm) when the file states no uncertainty. pub const DEFAULT_TOLERANCE: f64 = 1e-3; @@ -553,8 +554,7 @@ impl<'b, 'a> Builder<'b, 'a> { if self.faces.len() >= MAX_FACES { return Err(format!("more than {MAX_FACES} faces")); } - let surface = self.reader.surface(rid(p.get(2), "face geometry")?).map_err(|m| format!("face #{id}: {m}"))?; - let surface = map_surface(&surface).map_err(|m| format!("face #{id}: {m}"))?; + let step_surface = self.reader.surface(rid(p.get(2), "face geometry")?).map_err(|m| format!("face #{id}: {m}"))?; let same_sense = flag(p.get(3), "face same_sense")?; let bounds = refs(p.get(1), "face bounds")?; if bounds.len() > MAX_BOUNDS { @@ -564,10 +564,7 @@ impl<'b, 'a> Builder<'b, 'a> { let mut outer = None; let mut edges = 0; for b in bounds { - let (mut l, is_outer) = self.bound(b, &mut edges).map_err(|m| format!("face #{id}: {m}"))?; - if !same_sense && let LoopB::Edges(uses) = &mut l { - reverse(uses); - } + let (l, is_outer) = self.bound(b, &mut edges).map_err(|m| format!("face #{id}: {m}"))?; if is_outer && outer.is_none() { outer = Some(loops.len()); } @@ -576,9 +573,53 @@ impl<'b, 'a> Builder<'b, 'a> { if !loops.iter().any(|l| matches!(l, LoopB::Edges(_))) { return Err(format!("face #{id}: {}", unsupported("a face without an edge loop"))); } + let extent = self.extent(&loops)?; + let (surface, reversed) = swept::map_surface(&step_surface, &extent).map_err(|m| format!("face #{id}: {m}"))?; + // Loops relative to the record's natural normal: the face's normal is STEP's surface + // normal when same_sense, and the record's is STEP's unless `reversed`. + let sense = same_sense != reversed; + if !sense { + for l in &mut loops { + if let LoopB::Edges(uses) = l { + reverse(uses); + } + } + } let i = index(self.faces.len())?; - self.faces.push(FaceB { step: id, surface, loops, outer, forward: same_sense }); - Ok((i, same_sense)) + self.faces.push(FaceB { step: id, surface, loops, outer, forward: sense }); + Ok((i, sense)) + } + + /// Points along a face's loops (each edge at nine points, and its vertex loops), and along + /// the edges its loops use twice (its seams). + fn extent(&self, loops: &[LoopB]) -> R { + let mut ext = swept::Extent::default(); + let mut count: HashMap = HashMap::new(); + for l in loops { + match l { + LoopB::Edges(es) => { + for (e, _) in es { + *count.entry(*e).or_default() += 1; + } + } + LoopB::Vertex(v) => ext.points.push(v3(self.point(*v)?)), + } + } + let mut edges: Vec = count.keys().copied().collect(); + edges.sort_unstable(); + for e in edges { + let edge = self.edges.get(e as usize).ok_or("edge reference")?; + let c = CurveEval::new(&edge.curve).map_err(|m| m.to_string())?; + let [lo, hi] = edge.interval; + for k in 0..=8 { + let p = c.point(lo + (hi - lo) * f64::from(k) / 8.0).ok_or("edge point")?; + ext.points.push(p); + if count.get(&e).is_some_and(|n| *n > 1) { + ext.seam.push(p); + } + } + } + Ok(ext) } fn shell(&mut self, id: u64) -> R> { @@ -1119,40 +1160,6 @@ fn map_bspline_surface(b: &BSplineSurface) -> R { }) } -/// A face's surface as a schema surface. -fn map_surface(s: &Surface) -> R { - map_surface_in(s, 0) -} - -fn map_surface_in(s: &Surface, depth: usize) -> R { - if depth > MAX_NESTING { - return Err("surfaces nested too deeply".into()); - } - Ok(match s { - Surface::Plane { position } => BSurface::Plane { origin: position.origin, u_axis: position.ref_direction, v_axis: position.y() }, - Surface::Cylinder { position, radius } => BSurface::Cylinder { frame: frame(position), radius: *radius }, - Surface::Cone { position, radius, semi_angle } => { - if !(*radius >= 0.0) || !(semi_angle.abs() > 1e-12 && semi_angle.abs() < std::f64::consts::FRAC_PI_2) { - return Err(unsupported(format!("a cone of radius {radius} and semi-angle {semi_angle}"))); - } - BSurface::Cone { frame: frame(position), radius_at_origin: *radius, half_angle: *semi_angle } - } - Surface::Sphere { position, radius } => BSurface::Sphere { frame: frame(position), radius: *radius }, - Surface::Torus { position, major_radius, minor_radius } => { - if !(major_radius > minor_radius) { - return Err(unsupported("a spindle or horn torus")); - } - BSurface::Torus { frame: frame(position), major_radius: *major_radius, minor_radius: *minor_radius } - } - Surface::DegenerateTorus { .. } => return Err(unsupported("DEGENERATE_TOROIDAL_SURFACE")), - Surface::BSpline(b) => map_bspline_surface(b)?, - Surface::RectangularTrimmed { basis, .. } => map_surface_in(basis, depth + 1)?, - Surface::Revolution { .. } => return Err(unsupported("SURFACE_OF_REVOLUTION")), - Surface::Extrusion { .. } => return Err(unsupported("SURFACE_OF_LINEAR_EXTRUSION")), - Surface::Offset { .. } => return Err(unsupported("OFFSET_SURFACE")), - }) -} - /// An edge on its schema curve. #[derive(Debug)] struct Mapped { @@ -1213,13 +1220,35 @@ fn map_edge(curve: &Curve, ps: P3, pe: P3, same_sense: bool, closed: bool, near: } open_or_closed(BCurve::Polyline { points: points.clone() }, ps, pe, same_sense, closed, tol) } - Curve::Hyperbola { .. } => Err(unsupported("HYPERBOLA")), - Curve::Parabola { .. } => Err(unsupported("PARABOLA")), + Curve::Hyperbola { position, semi_axis, semi_imag_axis } => { + open_conic(curve, ps, pe, same_sense, closed, |lo, hi| swept::hyperbola_spline(position, *semi_axis, *semi_imag_axis, lo, hi)) + } + Curve::Parabola { position, focal_dist } => { + open_conic(curve, ps, pe, same_sense, closed, |lo, hi| swept::parabola_spline(position, *focal_dist, lo, hi)) + } Curve::OffsetCurve3d { .. } => Err(unsupported("OFFSET_CURVE_3D")), Curve::Trimmed { .. } | Curve::SurfaceCurve { .. } => Err(unsupported("nested trimmed or surface curves")), } } +/// A hyperbola or parabola: the vertices' parameters in closed form, the edge running the way +/// they go, as an exact spline over just that interval. +fn open_conic(curve: &Curve, ps: P3, pe: P3, same_sense: bool, closed: bool, spline: impl Fn(f64, f64) -> BCurve) -> R { + if closed { + return Err("a closed edge on an open conic".into()); + } + let (ts, te) = swept::open_conic_params(curve, ps, pe)?; + if !((te - ts).abs() > 1e-12 * (1.0 + ts.abs().max(te.abs()))) { + return Err("a degenerate edge on an open conic".into()); + } + if (ts - te).abs() > 64.0 { + return Err(unsupported("a hyperbola or parabola edge that long")); + } + let with_step = ts < te; + let (lo, hi) = if with_step { (ts, te) } else { (te, ts) }; + Ok(Mapped { curve: spline(lo, hi), interval: [lo, hi], with_step, sense_mismatch: with_step != same_sense }) +} + /// A circle or ellipse: the vertices' angles, the arc between them in the edge's sense (a full /// turn when they coincide). fn conic(curve: BCurve, position: &Placement, radii: [f64; 2], ps: P3, pe: P3, same_sense: bool, near: bool) -> R { diff --git a/crates/step/src/topology/seams.rs b/crates/step/src/topology/seams.rs index 424fbb8bc4..cfcba40774 100644 --- a/crates/step/src/topology/seams.rs +++ b/crates/step/src/topology/seams.rs @@ -4,7 +4,7 @@ //! apex); schema 4 keeps every loop closed in the surface's (u, v), so such a face gets a seam: //! a ruling line on a cylinder or cone, a meridian arc on a sphere, from a vertex of one loop to //! the other loop (split there when it has no vertex at that angle) or to the pole the face -//! covers. The loops merge into one: the first, the seam out, the second (or nothing at a pole), +//! covers; on a surface of revolution, its profile turned to the vertex's angle. The loops merge into one: the first, the seam out, the second (or nothing at a pole), //! the seam back. A torus face winding without a seam is refused for now. use std::f64::consts::{PI, TAU}; @@ -37,10 +37,12 @@ struct Axis { poles: [Option; 2], } -#[derive(Clone, Copy, PartialEq)] +#[derive(Clone, PartialEq)] enum Kind { Ruled, Sphere(f64), + /// A placed revolution with u the angle: its profile. + Revolved(BCurve), } fn axis(s: &BSurface) -> Option { @@ -63,10 +65,53 @@ fn axis(s: &BSurface) -> Option { let (origin, z) = frame(f); Some(Axis { origin, z, kind: Kind::Sphere(*radius), poles: [Some(origin - z * *radius), Some(origin + z * *radius)] }) } + BSurface::PlacedRevolution { profile, origin, axis, matrix, forward: false } if *matrix == IDENTITY => { + let (origin, z) = (Vec3::from(*origin), Vec3::from(*axis)); + let c = CurveEval::new(profile).ok()?; + let solidcraft_geom::brep::Natural::Range(a, b) = profile.natural() else { return None }; + // A profile end on the axis is a pole. + let pole = |t: f64| { + c.point(t).filter(|p| { + let d = *p - origin; + (d - z * d.dot(z)).len() <= 1e-9 * (1.0 + d.len()) + }) + }; + Some(Axis { origin, z, kind: Kind::Revolved((**profile).clone()), poles: [pole(a), pole(b)] }) + } _ => None, } } +const IDENTITY: [[f64; 4]; 4] = [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]; + +/// `p` turned by `angle` about the axis through `o` along unit `k` (right-handed). +fn turned(p: Vec3, o: Vec3, k: Vec3, angle: f64) -> Vec3 { + let (s, c) = angle.sin_cos(); + let x = p - o; + o + x * c + k.cross(x) * s + k * (k.dot(x) * (1.0 - c)) +} + +/// A revolution's profile turned by `angle`: the meridian there, with the profile's parameter. +fn meridian(profile: &BCurve, o: Vec3, k: Vec3, angle: f64) -> R { + let turn = |h: &[f64; 4]| { + let w = h[3]; + let p = turned(Vec3::new(h[0] / w, h[1] / w, h[2] / w), o, k, angle); + [p.x * w, p.y * w, p.z * w, w] + }; + Ok(match profile { + BCurve::Spline { knots, points } => { + BCurve::Spline { knots: knots.clone(), points: points.iter().map(|p| arr(turned(Vec3::from(*p), o, k, angle))).collect() } + } + BCurve::RationalSpline { knots, homogeneous } => { + BCurve::RationalSpline { knots: knots.clone(), homogeneous: homogeneous.iter().map(turn).collect() } + } + BCurve::Line { start, end } => { + BCurve::Line { start: arr(turned(Vec3::from(*start), o, k, angle)), end: arr(turned(Vec3::from(*end), o, k, angle)) } + } + _ => return Err(unsupported("a revolution of this profile")), + }) +} + impl Builder<'_, '_> { /// Give every face that winds around its surface without a seam one. Returns how many /// seams were added. @@ -104,8 +149,8 @@ impl Builder<'_, '_> { }; match (up.as_slice(), down.as_slice()) { // The face lies to the left of its loops in (u, v): above a loop running +u. - ([l], []) => self.pole_seam(fi, *l, &ax, ax.poles[1])?, - ([], [l]) => self.pole_seam(fi, *l, &ax, ax.poles[0])?, + ([l], []) => self.pole_seam(fi, *l, &ax, ax.poles[1], &eval)?, + ([], [l]) => self.pole_seam(fi, *l, &ax, ax.poles[0], &eval)?, ([a], [b]) => self.band_seam(fi, *a, *b, &ax, &eval, &mut splits)?, _ => return Err(format!("face #{step}: {}", unsupported("a face with several loops around its axis"))), } @@ -156,12 +201,16 @@ impl Builder<'_, '_> { Ok(es.iter().skip(k).chain(es.iter().take(k)).copied().collect()) } - /// A seam edge from vertex `a` to vertex `b` on the surface: a ruling line, or a meridian - /// arc of a sphere. - fn seam_edge(&mut self, a: u32, b: u32, ax: &Axis, eval: Option<&SurfaceEval>) -> R { + /// A seam edge from vertex `a` to vertex `b` on the surface: a ruling line, a meridian arc + /// of a sphere, or a revolution's profile turned to their angle. Returns the edge and + /// whether it runs from `a` to `b`. + fn seam_edge(&mut self, a: u32, b: u32, ax: &Axis, eval: &SurfaceEval, check: bool) -> R<(u32, bool)> { let (pa, pb) = (self.vpoint(a)?, self.vpoint(b)?); - let (curve, interval) = match ax.kind { - Kind::Ruled => (BCurve::Line { start: arr(pa), end: arr(pb) }, [0.0, 1.0]), + if pa.dist(pb) <= self.tol { + return Err("a seam of no length".into()); + } + let (curve, interval, along) = match &ax.kind { + Kind::Ruled => (BCurve::Line { start: arr(pa), end: arr(pb) }, [0.0, 1.0], true), Kind::Sphere(r) => { let x = (pa - ax.origin).normalized().ok_or("a seam from the sphere's centre")?; let to = pb - ax.origin; @@ -170,37 +219,52 @@ impl Builder<'_, '_> { if !(theta > 0.0 && theta <= PI) { return Err("a degenerate seam".into()); } - (BCurve::Circle { frame: Frame { origin: arr(ax.origin), x_axis: arr(x), y_axis: arr(y) }, radius: r }, [0.0, theta]) + (BCurve::Circle { frame: Frame { origin: arr(ax.origin), x_axis: arr(x), y_axis: arr(y) }, radius: *r }, [0.0, theta], true) + } + Kind::Revolved(profile) => { + // a's angle; each vertex's profile parameter (a pole's from the profile's end). + let qa = eval.project(pa, None, f64::INFINITY).ok_or("a seam vertex off its surface")?; + let qb = eval.project(pb, Some(qa), f64::INFINITY).ok_or("a seam vertex off its surface")?; + let solidcraft_geom::brep::Natural::Range(lo, hi) = profile.natural() else { return Err("a periodic profile".into()) }; + let vb = match ax.poles { + [Some(p), _] if p.dist(pb) <= self.tol => lo, + [_, Some(p)] if p.dist(pb) <= self.tol => hi, + _ => qb[1], + }; + let curve = meridian(profile, ax.origin, ax.z, qa[0])?; + let along = qa[1] < vb; + (curve, if along { [qa[1], vb] } else { [vb, qa[1]] }, along) } }; - if pa.dist(pb) <= self.tol { - return Err("a seam of no length".into()); + if !(interval[0] < interval[1]) { + return Err("a degenerate seam".into()); } // On the surface: its middle within the tolerance (a ruling, a meridian). let c = CurveEval::new(&curve).map_err(|m| m.to_string())?; - if let Some(eval) = eval { + if check { let m = c.point(0.5 * (interval[0] + interval[1])).ok_or("seam point")?; let off = eval.project(m, None, f64::INFINITY).and_then(|q| eval.point(q)).map_or(f64::INFINITY, |s| s.dist(m)); if !(off <= self.tol) { return Err(unsupported("a seam off its surface: loops with no vertices at one angle")); } } - for (t, v) in interval.into_iter().zip([a, b]) { + let vertices = if along { [a, b] } else { [b, a] }; + for (t, v) in interval.into_iter().zip(vertices) { let gap = c.point(t).ok_or("seam point")?.dist(self.vpoint(v)?); if let Some(x) = self.vertices.get_mut(v as usize) { x.tolerance = x.tolerance.max(1.5 * gap); } } let i = index(self.edges.len())?; - self.edges.push(EdgeB { curve, interval, vertices: [a, b], tolerance: self.tol }); - Ok(i) + self.edges.push(EdgeB { curve, interval, vertices, tolerance: self.tol }); + Ok((i, along)) } /// Join loop `li` of face `fi` to the pole the face covers: seam out and back. - fn pole_seam(&mut self, fi: usize, li: usize, ax: &Axis, pole: Option) -> R<()> { + fn pole_seam(&mut self, fi: usize, li: usize, ax: &Axis, pole: Option, eval: &SurfaceEval) -> R<()> { let step = self.faces.get(fi).map_or(0, |f| f.step); let Some(p) = pole else { - return Err(format!("face #{step}: a face around an open cylinder or cone with one bound")); + return Err(format!("face #{step}: a face around an open surface of revolution with one bound")); }; let Some(LoopB::Edges(es)) = self.faces.get(fi).and_then(|f| f.loops.get(li)) else { return Err("loop reference".into()) }; let es = es.clone(); @@ -229,10 +293,10 @@ impl Builder<'_, '_> { v } }; - let seam = self.seam_edge(a, b, ax, None)?; + let (seam, along) = self.seam_edge(a, b, ax, eval, false)?; let mut merged = ring; - merged.push((seam, true)); - merged.push((seam, false)); + merged.push((seam, along)); + merged.push((seam, !along)); let face = self.faces.get_mut(fi).ok_or("face reference")?; if let Some(l) = face.loops.get_mut(li) { *l = LoopB::Edges(merged); @@ -306,11 +370,11 @@ impl Builder<'_, '_> { let b_ring = self.rotated(fi, lb, j)?; let va = self.start(*a_ring.first().ok_or("empty loop")?)?; let vb = self.start(*b_ring.first().ok_or("empty loop")?)?; - let seam = self.seam_edge(va, vb, ax, Some(eval))?; + let (seam, along) = self.seam_edge(va, vb, ax, eval, true)?; let mut merged = a_ring; - merged.push((seam, true)); + merged.push((seam, along)); merged.extend(b_ring); - merged.push((seam, false)); + merged.push((seam, !along)); let face = self.faces.get_mut(fi).ok_or("face reference")?; if let Some(l) = face.loops.get_mut(la) { *l = LoopB::Edges(merged); diff --git a/crates/step/src/topology/swept.rs b/crates/step/src/topology/swept.rs new file mode 100644 index 0000000000..aee6ee9e81 --- /dev/null +++ b/crates/step/src/topology/swept.rs @@ -0,0 +1,583 @@ +//! Swept, offset and open-conic geometry as exact schema 4 kinds. STEP's SURFACE_OF_REVOLUTION +//! becomes a `placed_revolution` (u the angle, v the profile's parameter, as STEP has them); +//! SURFACE_OF_LINEAR_EXTRUSION becomes a spline surface, the profile as an exact (rational) +//! B-spline times degree 1 along the extrusion vector (a circle extruded along its axis, and a +//! circle revolved in a plane through the axis, become the cylinder and ring torus they are, so +//! whole turns keep a periodic parameter); OFFSET_SURFACE of an analytic surface +//! becomes that surface moved along its normal; HYPERBOLA and PARABOLA become exact rational +//! quadratic and quadratic splines. +//! +//! Profiles are bounded where the face needs them: a line, a conic or an open conic over the +//! parameters of the face's boundary points (with a margin), a closed conic starting at the +//! face's seam or in the widest gap between them, so its seam lies outside the face. A +//! revolved profile stops where it meets the axis (beyond it the surface would fold over). The +//! surfaces' parameters may differ from STEP's (an extrusion's v, a conic's rational +//! parameter); pcurves and domains come from the records, never from STEP's parameters. + +use std::f64::consts::{FRAC_PI_2, TAU}; + +use solidcraft_geom::Vec3; +use solidcraft_geom::brep::{Curve as BCurve, Surface as BSurface}; + +use super::{MAX_NESTING, R, frame, map_bspline_curve, unsupported, v3}; +use crate::geom::{Axis1, Curve, P3, Placement, Surface}; + +/// Margin around a face's parameters on a profile, as a fraction of their span. +const MARGIN: f64 = 0.1; +/// Most sweep of one rational quadratic conic piece. +const CONIC_PIECE: f64 = FRAC_PI_2; +/// Most parameter step of one rational quadratic hyperbola piece (its weight `cosh(step / 2)` +/// stays near 1). +const HYPERBOLA_STEP: f64 = 1.0; +/// Most hyperbola pieces. +const MAX_PIECES: usize = 64; +/// Samples when looking for where a revolved profile meets its axis. +const AXIS_SAMPLES: usize = 128; + +/// Points of a face's boundary (its edges sampled) and of its seam edges (used twice). +#[derive(Default)] +pub(super) struct Extent { + pub points: Vec, + pub seam: Vec, +} + +fn p3(v: Vec3) -> P3 { + [v.x, v.y, v.z] +} + +/// The ellipse of `radii` in `p` from angle `t0` over `sweep` > 0 as rational quadratic pieces +/// whose knots are the pieces' end angles. +pub(super) fn conic_spline(p: &Placement, radii: [f64; 2], t0: f64, sweep: f64) -> BCurve { + let n = ((sweep / CONIC_PIECE).ceil() as usize).clamp(1, 8); + let dt = sweep / n as f64; + let w = (0.5 * dt).cos(); + let at = |c: f64, s: f64| p.at(radii[0] * c, radii[1] * s, 0.0); + let mut knots = vec![t0; 3]; + let first = at(t0.cos(), t0.sin()); + let mut homogeneous = vec![[first[0], first[1], first[2], 1.0]]; + for i in 0..n { + let ta = t0 + dt * i as f64; + let (tm, tb) = (ta + 0.5 * dt, ta + dt); + // The meeting point of the end tangents, weighted cos of the half sweep. + let m = at(tm.cos() / w, tm.sin() / w); + homogeneous.push([m[0] * w, m[1] * w, m[2] * w, w]); + let b = at(tb.cos(), tb.sin()); + homogeneous.push([b[0], b[1], b[2], 1.0]); + let k = if i + 1 == n { t0 + sweep } else { tb }; + knots.extend([k, k]); + } + knots.push(t0 + sweep); + BCurve::RationalSpline { knots, homogeneous } +} + +/// The hyperbola branch `a cosh t X + b sinh t Y` from `t0` to `t1` (t0 < t1) as rational +/// quadratic pieces whose knots are the pieces' end parameters. +pub(super) fn hyperbola_spline(p: &Placement, a: f64, b: f64, t0: f64, t1: f64) -> BCurve { + let n = (((t1 - t0) / HYPERBOLA_STEP).ceil() as usize).clamp(1, MAX_PIECES); + let dt = (t1 - t0) / n as f64; + let w = (0.5 * dt).cosh(); + let at = |c: f64, s: f64| p.at(a * c, b * s, 0.0); + let mut knots = vec![t0; 3]; + let first = at(t0.cosh(), t0.sinh()); + let mut homogeneous = vec![[first[0], first[1], first[2], 1.0]]; + for i in 0..n { + let ta = t0 + dt * i as f64; + let (tm, tb) = (ta + 0.5 * dt, ta + dt); + let m = at(tm.cosh() / w, tm.sinh() / w); + homogeneous.push([m[0] * w, m[1] * w, m[2] * w, w]); + let e = at(tb.cosh(), tb.sinh()); + homogeneous.push([e[0], e[1], e[2], 1.0]); + let k = if i + 1 == n { t1 } else { tb }; + knots.extend([k, k]); + } + knots.push(t1); + BCurve::RationalSpline { knots, homogeneous } +} + +/// The parabola `f t² X + 2 f t Y` from `t0` to `t1` (t0 < t1): one quadratic Bézier on knots +/// `t0, t1`, so its parameter is STEP's t exactly. +pub(super) fn parabola_spline(p: &Placement, f: f64, t0: f64, t1: f64) -> BCurve { + let at = |t: f64| p.at(f * t * t, 2.0 * f * t, 0.0); + // The end tangents meet at f·(t0·t1, t0 + t1). + let mid = p.at(f * t0 * t1, f * (t0 + t1), 0.0); + BCurve::Spline { knots: vec![t0, t0, t0, t1, t1, t1], points: vec![at(t0), mid, at(t1)] } +} + +/// The STEP line from `t0` to `t1` as a degree-1 spline on knots `t0, t1`: STEP's parameter. +fn line_spline(origin: P3, direction: P3, magnitude: f64, t0: f64, t1: f64) -> BCurve { + let at = |t: f64| [0, 1, 2].map(|i| origin[i] + direction[i] * magnitude * t); + BCurve::Spline { knots: vec![t0, t0, t1, t1], points: vec![at(t0), at(t1)] } +} + +/// A polyline as a degree-1 spline: segment i on [i, i + 1], as STEP has it. +fn polyline_spline(points: &[P3]) -> BCurve { + let n = points.len(); + let mut knots = vec![0.0]; + knots.extend((0..n).map(|i| i as f64)); + knots.push(n.saturating_sub(1) as f64); + BCurve::Spline { knots, points: points.to_vec() } +} + +/// A curve under its SURFACE_CURVE and TRIMMED_CURVE wrappers. +fn basis(c: &Curve) -> R<&Curve> { + let mut c = c; + for _ in 0..=MAX_NESTING { + match c { + Curve::SurfaceCurve { curve_3d, .. } => c = curve_3d, + Curve::Trimmed { basis, .. } => c = basis, + _ => return Ok(c), + } + } + Err("profile curves nested too deeply".into()) +} + +/// The conic angle of a point in the conic's plane. +fn conic_angle(p: &Placement, radii: [f64; 2], q: Vec3) -> f64 { + let d = q - v3(p.origin); + (d.dot(v3(p.y())) / radii[1]).atan2(d.dot(v3(p.ref_direction)) / radii[0]) +} + +/// The STEP parameter of `q` (a point on, or projected into the plane of, a planar curve): +/// angle, hyperbola or parabola parameter, or line parameter. +fn planar_param(c: &Curve, q: Vec3) -> Option { + let local = |p: &Placement| { + let d = q - v3(p.origin); + (d.dot(v3(p.ref_direction)), d.dot(v3(p.y()))) + }; + let t = match c { + Curve::Line { origin, direction, magnitude } => (q - v3(*origin)).dot(v3(*direction)) / magnitude, + Curve::Circle { position, radius } => conic_angle(position, [*radius; 2], q), + Curve::Ellipse { position, semi_axis_1, semi_axis_2 } => conic_angle(position, [*semi_axis_1, *semi_axis_2], q), + Curve::Hyperbola { position, semi_imag_axis, .. } => (local(position).1 / semi_imag_axis).asinh(), + Curve::Parabola { position, focal_dist } => local(position).1 / (2.0 * focal_dist), + _ => return None, + }; + t.is_finite().then_some(t) +} + +/// [min, max] of the parameters, or `None` without any. +fn range(ts: &[f64]) -> Option<(f64, f64)> { + let (lo, hi) = ts.iter().fold((f64::INFINITY, f64::NEG_INFINITY), |(a, b), t| (a.min(*t), b.max(*t))); + (lo.is_finite() && hi.is_finite()).then_some((lo, hi)) +} + +/// Middle-free description of the widest gap between angles: (where it starts, its width). +fn widest_gap(angles: &[f64]) -> Option<(f64, f64)> { + let mut a: Vec = angles.iter().map(|x| x.rem_euclid(TAU)).filter(|x| x.is_finite()).collect(); + a.sort_by(f64::total_cmp); + let (first, last) = (*a.first()?, *a.last()?); + let mut best = (last, first + TAU - last); + for w in a.windows(2) { + if let [x, y] = w + && y - x > best.1 + { + best = (*x, y - x); + } + } + Some(best) +} + +/// The angle all the given angles share (within 1e-6 rad), if there are any. +fn common_angle(angles: &[f64]) -> Option { + let first = *angles.first()?; + let off = |a: f64| ((a - first) + 0.5 * TAU).rem_euclid(TAU) - 0.5 * TAU; + angles.iter().all(|a| off(*a).abs() < 1e-6).then_some(first) +} + +/// A profile's parameter range: where the face is (`inner`) and the bounded range with its +/// margin (`outer`), or a whole closed conic from a start angle. +enum Span { + Open { inner: (f64, f64), outer: (f64, f64) }, + Turn { start: f64 }, +} + +/// The range a profile is built over, from the parameters of the face's points on it and of its +/// seam points; `closed` for a circle or ellipse. +fn span(ts: &[f64], seam: &[f64], closed: bool) -> R { + if closed { + if let Some(s) = common_angle(seam) { + return Ok(Span::Turn { start: s }); + } + let (gap_start, gap) = widest_gap(ts).ok_or("a face without boundary points")?; + let covered = TAU - gap; + if gap < 1e-3 { + return Ok(Span::Turn { start: gap_start + gap }); + } + let start = gap_start + gap; + let m = (MARGIN * covered + 1e-3).min(0.25 * gap); + return Ok(Span::Open { inner: (start, start + covered), outer: (start - m, start + covered + m) }); + } + let (lo, hi) = range(ts).ok_or("a face without boundary points")?; + let m = MARGIN * (hi - lo) + 1e-3 * (1.0 + lo.abs().max(hi.abs())); + Ok(Span::Open { inner: (lo, hi), outer: (lo - m, hi + m) }) +} + +/// The profile over its span as an exact spline (B-spline and polyline profiles whole). +fn bounded(c: &Curve, s: Span) -> R { + let (lo, hi) = match s { + Span::Open { outer, .. } => outer, + Span::Turn { start } => (start, start + TAU), + }; + if !(lo < hi && lo.is_finite() && hi.is_finite()) { + return Err("a profile without extent".into()); + } + Ok(match c { + Curve::Line { origin, direction, magnitude } => line_spline(*origin, *direction, *magnitude, lo, hi), + Curve::Circle { position, radius } => conic_spline(position, [*radius; 2], lo, hi - lo), + Curve::Ellipse { position, semi_axis_1, semi_axis_2 } => conic_spline(position, [*semi_axis_1, *semi_axis_2], lo, hi - lo), + Curve::Hyperbola { position, semi_axis, semi_imag_axis } => { + if hi - lo > HYPERBOLA_STEP * MAX_PIECES as f64 { + return Err(unsupported("a hyperbola profile that long")); + } + hyperbola_spline(position, *semi_axis, *semi_imag_axis, lo, hi) + } + Curve::Parabola { position, focal_dist } => parabola_spline(position, *focal_dist, lo, hi), + _ => return Err(unsupported("this profile curve")), + }) +} + +/// A profile that needs no bounds (a B-spline, a polyline), if it is one. +fn whole(c: &Curve) -> R> { + Ok(match c { + Curve::BSpline(b) => Some(map_bspline_curve(b)?), + Curve::Polyline(points) => { + if points.len() < 2 || points.len() > super::MAX_SPLINE_POINTS || points.windows(2).any(|w| w.first() == w.get(1)) { + return Err(unsupported(format!("a polyline profile of {} points", points.len()))); + } + Some(polyline_spline(points)) + } + Curve::Line { .. } | Curve::Circle { .. } | Curve::Ellipse { .. } | Curve::Hyperbola { .. } | Curve::Parabola { .. } => None, + Curve::OffsetCurve3d { .. } => return Err(unsupported("an OFFSET_CURVE_3D profile")), + Curve::Trimmed { .. } | Curve::SurfaceCurve { .. } => return Err(unsupported("nested profile curves")), + }) +} + +/// The middle parameter of a TRIMMED_CURVE trimmed by parameters (on a closed conic, the +/// middle of the arc its sense runs over). +fn trimmed_middle(c: &Curve) -> Option { + let Curve::Trimmed { basis, trim_1, trim_2, sense, .. } = c else { return None }; + let param = |t: &[crate::geom::Trim]| { + t.iter().find_map(|x| match x { + crate::geom::Trim::Parameter(p) => Some(*p), + crate::geom::Trim::Point(_) => None, + }) + }; + let (a, mut b) = (param(trim_1)?, param(trim_2)?); + if is_closed(basis) { + b = if *sense { a + (b - a).rem_euclid(TAU) } else { a - (a - b).rem_euclid(TAU) }; + } + let m = 0.5 * (a + b); + m.is_finite().then_some(m) +} + +fn is_closed(c: &Curve) -> bool { + matches!(c, Curve::Circle { .. } | Curve::Ellipse { .. }) +} + +/// SURFACE_OF_REVOLUTION: the profile turned about the axis, u the angle and v the profile's +/// parameter (`forward: false`), so the natural normal is ∂S/∂angle × ∂S/∂profile, STEP's. +fn revolution(profile: &Curve, axis: &Axis1, ext: &Extent) -> R { + let c = basis(profile)?; + let (o, k) = (v3(axis.origin), v3(axis.direction)); + let k = k.normalized().ok_or("a revolution about no direction")?; + let curve = match whole(c)? { + Some(curve) => curve, + None => { + // The meridian half-plane the profile lies in: the radial direction of its point + // farthest from the axis among a few. + let radial = |q: Vec3| { + let d = q - o; + d - k * d.dot(k) + }; + let probe = |t: f64| c.point(t).map(v3); + let mut e = None::; + let mut best = 0.0; + // A trimmed profile's own middle says which side of the axis it is on (a profile + // across the axis covers the surface twice, once each way round). + let probes = match trimmed_middle(profile) { + Some(t) => vec![t], + None => vec![-1.0, 0.0, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0], + }; + for t in probes { + if let Some(r) = probe(t).map(radial) + && r.len() > best + { + best = r.len(); + e = r.normalized(); + } + } + let e = e.ok_or_else(|| unsupported("a profile on its axis"))?; + // Each boundary point turned into the meridian plane, then its profile parameter: + // on whichever side of the axis the profile has a point where it lands. + let meridian = |q: Vec3| ((q - o).dot(k), radial(q).len()); + let param = |q: Vec3| { + let (h, r) = meridian(q); + [e, -e] + .into_iter() + .filter_map(|side| { + let at = o + k * h + side * r; + let t = planar_param(c, at)?; + Some((t, probe(t)?.dist(at))) + }) + .min_by(|a, b| a.1.total_cmp(&b.1)) + .map(|(t, _)| t) + }; + let ts: Vec = ext.points.iter().filter_map(|q| param(*q)).collect(); + let seam: Vec = ext.seam.iter().filter_map(|q| param(*q)).collect(); + // The profile's side of the axis where the face is. + let mid = ts.get(ts.len() / 2).copied().ok_or("a face without boundary points")?; + let e = probe(mid).map(radial).and_then(|r| r.normalized()).unwrap_or(e); + let s = span(&ts, &seam, is_closed(c))?; + // Stop where the profile crosses the axis outside the face. + let side = |t: f64| probe(t).map(|q| (q - o).dot(e)); + let s = match s { + Span::Open { inner, outer } => { + let lo = crossing(&side, inner.0, outer.0).unwrap_or(outer.0); + let hi = crossing(&side, inner.1, outer.1).unwrap_or(outer.1); + Span::Open { inner, outer: (lo, hi) } + } + turn => turn, + }; + bounded(c, s)? + } + }; + let identity = [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]; + Ok(BSurface::PlacedRevolution { profile: Box::new(curve), origin: p3(o), axis: p3(k), matrix: identity, forward: false }) +} + +/// Going from `from` (where `side` is positive or zero) towards `to`, the first parameter +/// where `side` reaches zero, if it does. +fn crossing(side: &dyn Fn(f64) -> Option, from: f64, to: f64) -> Option { + let s0 = side(from)?; + let mut prev = (from, s0); + for i in 1..=AXIS_SAMPLES { + let t = from + (to - from) * i as f64 / AXIS_SAMPLES as f64; + let s = side(t)?; + if s.signum() != prev.1.signum() || s == 0.0 { + let (mut a, mut b, fa) = (prev.0, t, prev.1); + for _ in 0..100 { + let m = 0.5 * (a + b); + let fm = side(m)?; + if fm.signum() == fa.signum() && fm != 0.0 { + a = m; + } else { + b = m; + } + } + return Some(b); + } + prev = (t, s); + } + None +} + +/// SURFACE_OF_LINEAR_EXTRUSION: the profile as an exact spline along u, degree 1 along the +/// vector in v, over the face's reach along it. +fn extrusion(profile: &Curve, direction: P3, magnitude: f64, ext: &Extent) -> R { + let c = basis(profile)?; + let d = v3(direction) * magnitude; + let dd = d.len2(); + if !(dd > 0.0 && dd.is_finite()) { + return Err("an extrusion of no length".into()); + } + let curve = match whole(c)? { + Some(curve) => curve, + None => { + // Each boundary point slid along the vector onto the profile. + let onto = |q: Vec3| -> Option { + match c { + Curve::Line { origin, direction: e, magnitude: m } => { + // q − origin = s·m·e + v·d, least squares in (s, v). + let (e, r) = (v3(*e) * *m, q - v3(*origin)); + let (a, b, cc) = (e.len2(), e.dot(d), dd); + let det = a * cc - b * b; + (det.abs() > 1e-12 * a * cc).then(|| (r.dot(e) * cc - r.dot(d) * b) / det) + } + Curve::Circle { position, .. } + | Curve::Ellipse { position, .. } + | Curve::Hyperbola { position, .. } + | Curve::Parabola { position, .. } => { + let n = v3(position.axis); + let dn = d.dot(n); + if !(dn.abs() > 1e-9 * d.len()) { + return None; + } + planar_param(c, q - d * ((q - v3(position.origin)).dot(n) / dn)) + } + _ => None, + } + }; + let ts: Vec = ext.points.iter().filter_map(|q| onto(*q)).collect(); + if ts.len() < ext.points.len() { + return Err(unsupported("an extrusion along its profile's plane")); + } + let seam: Vec = ext.seam.iter().filter_map(|q| onto(*q)).collect(); + bounded(c, span(&ts, &seam, is_closed(c))?)? + } + }; + let (knots, homogeneous, rational) = match curve { + BCurve::Spline { knots, points } => (knots, points.iter().map(|p| [p[0], p[1], p[2], 1.0]).collect::>(), false), + BCurve::RationalSpline { knots, homogeneous } => (knots, homogeneous, true), + _ => return Err("an extrusion profile".into()), + }; + // v: the face's points' reach along the vector, less the profile's (its control points + // bound it), with a margin. + let along = |p: Vec3| p.dot(d) / dd; + let (plo, phi) = range(&ext.points.iter().map(|p| along(*p)).collect::>()).ok_or("a face without boundary points")?; + let cs: Vec = homogeneous.iter().filter(|h| h[3] > 0.0).map(|h| along(Vec3::new(h[0] / h[3], h[1] / h[3], h[2] / h[3]))).collect(); + let (clo, chi) = range(&cs).ok_or("an extrusion profile")?; + let (lo, hi) = (plo - chi, phi - clo); + let m = MARGIN * (hi - lo) + 1e-3 / dd.sqrt(); + let (lo, hi) = (lo - m, hi + m); + let shift = |h: &[f64; 4], s: f64| [h[0] + d.x * s * h[3], h[1] + d.y * s * h[3], h[2] + d.z * s * h[3], h[3]]; + let rows: Vec> = homogeneous.iter().map(|h| vec![shift(h, lo), shift(h, hi)]).collect(); + let v_knots = vec![lo, lo, hi, hi]; + Ok(if rational { + BSurface::RationalSpline { u_knots: knots, v_knots, homogeneous: rows } + } else { + BSurface::Spline { u_knots: knots, v_knots, points: rows.into_iter().map(|r| r.into_iter().map(|h| [h[0], h[1], h[2]]).collect()).collect() } + }) +} + +/// An analytic surface moved `d` along its natural normal (every analytic kind's natural +/// normal points out of its axis, or along z × … for a plane): checked at each level. +fn offset(s: BSurface, d: f64) -> R { + let bad = |what: &str| Err(format!("an offset {what}")); + Ok(match s { + BSurface::Plane { origin, u_axis, v_axis } => { + let n = v3(u_axis).cross(v3(v_axis)).normalized().ok_or("a degenerate plane")?; + BSurface::Plane { origin: p3(v3(origin) + n * d), u_axis, v_axis } + } + BSurface::Cylinder { frame, radius } => { + if !(radius + d > 1e-9) { + return bad("cylinder that collapses or turns inside out"); + } + BSurface::Cylinder { frame, radius: radius + d } + } + BSurface::Sphere { frame, radius } => { + if !(radius + d > 1e-9) { + return bad("sphere that collapses or turns inside out"); + } + BSurface::Sphere { frame, radius: radius + d } + } + BSurface::Cone { frame, radius_at_origin, half_angle } => { + // The normal is (cos a, −sin a) in (radial, axial): the reference radius grows by + // d / cos a, the angle stays. + let r = radius_at_origin + d / half_angle.cos(); + if !(r >= 0.0 && r.is_finite()) { + return bad("cone whose reference radius turns negative"); + } + BSurface::Cone { frame, radius_at_origin: r, half_angle } + } + BSurface::Torus { frame, major_radius, minor_radius } => { + let r = minor_radius + d; + if !(r > 1e-9 && r < major_radius) { + return bad("torus that collapses or meets its axis"); + } + BSurface::Torus { frame, major_radius, minor_radius: r } + } + _ => return Err(unsupported("an OFFSET_SURFACE of a spline or swept surface")), + }) +} + +/// A face's surface as a schema surface, and whether the record's natural normal is opposite +/// STEP's surface normal (an extruded circle or revolved circle read as a cylinder or torus, +/// whose natural normals point away from their axes); `ext` bounds swept surfaces' profiles. +pub(super) fn map_surface(s: &Surface, ext: &Extent) -> R<(BSurface, bool)> { + let mut reversed = false; + let s = map_surface_in(s, ext, 0, &mut reversed)?; + Ok((s, reversed)) +} + +/// An extruded circle along its axis is a cylinder; a circle in a plane through the axis, +/// revolved, a ring torus. `None` for anything else. +fn canonical(s: &Surface, reversed: &mut bool) -> Option { + let c = match s { + Surface::Extrusion { profile, .. } | Surface::Revolution { profile, .. } => basis(profile).ok()?, + _ => return None, + }; + let Curve::Circle { position, radius } = c else { return None }; + let n = v3(position.axis); + match s { + Surface::Extrusion { direction, .. } => { + let d = v3(*direction); + if d.cross(n).len() > 1e-12 { + return None; + } + // STEP's normal is the profile's tangent × the vector: out of the axis when the + // vector runs along the circle's axis. + *reversed = d.dot(n) < 0.0; + Some(BSurface::Cylinder { frame: frame(position), radius: *radius }) + } + Surface::Revolution { axis, .. } => { + let (o, k) = (v3(axis.origin), v3(axis.direction).normalized()?); + let centre = v3(position.origin); + let scale = 1.0 + (centre - o).len() + radius; + if n.dot(k).abs() > 1e-12 || (centre - o).dot(n).abs() > 1e-12 * scale { + return None; + } + let foot = o + k * (centre - o).dot(k); + let big = (centre - foot).len(); + if !(big > *radius * (1.0 + 1e-9)) { + return None; + } + let x = (centre - foot).normalized()?; + // STEP's normal at the circle's angle 0: ∂angle × ∂profile; the torus's points out + // of the tube. + let (cx, cy) = (v3(position.ref_direction), v3(position.y())); + let at = centre + cx * *radius; + let step_normal = k.cross(at - o).cross(cy * *radius); + *reversed = step_normal.dot(cx) < 0.0; + Some(BSurface::Torus { + frame: solidcraft_geom::brep::Frame { origin: p3(foot), x_axis: p3(x), y_axis: p3(k.cross(x)) }, + major_radius: big, + minor_radius: *radius, + }) + } + _ => None, + } +} + +fn map_surface_in(s: &Surface, ext: &Extent, depth: usize, reversed: &mut bool) -> R { + if depth > MAX_NESTING { + return Err("surfaces nested too deeply".into()); + } + if let Some(c) = canonical(s, reversed) { + return Ok(c); + } + Ok(match s { + Surface::Plane { position } => BSurface::Plane { origin: position.origin, u_axis: position.ref_direction, v_axis: position.y() }, + Surface::Cylinder { position, radius } => BSurface::Cylinder { frame: frame(position), radius: *radius }, + Surface::Cone { position, radius, semi_angle } => { + if !(*radius >= 0.0) || !(semi_angle.abs() > 1e-12 && semi_angle.abs() < FRAC_PI_2) { + return Err(unsupported(format!("a cone of radius {radius} and semi-angle {semi_angle}"))); + } + BSurface::Cone { frame: frame(position), radius_at_origin: *radius, half_angle: *semi_angle } + } + Surface::Sphere { position, radius } => BSurface::Sphere { frame: frame(position), radius: *radius }, + Surface::Torus { position, major_radius, minor_radius } => { + if !(major_radius > minor_radius) { + return Err(unsupported("a spindle or horn torus")); + } + BSurface::Torus { frame: frame(position), major_radius: *major_radius, minor_radius: *minor_radius } + } + Surface::DegenerateTorus { .. } => return Err(unsupported("DEGENERATE_TOROIDAL_SURFACE")), + Surface::BSpline(b) => super::map_bspline_surface(b)?, + Surface::RectangularTrimmed { basis, .. } => map_surface_in(basis, ext, depth + 1, reversed)?, + Surface::Revolution { profile, axis } => revolution(profile, axis, ext)?, + Surface::Extrusion { profile, direction, magnitude } => extrusion(profile, *direction, *magnitude, ext)?, + Surface::Offset { basis, distance, .. } => { + // Along STEP's normal: the record's, or its opposite. + let b = map_surface_in(basis, ext, depth + 1, reversed)?; + offset(b, if *reversed { -*distance } else { *distance })? + } + }) +} + +/// The STEP parameters of `ps` and `pe` on a hyperbola or parabola. +pub(super) fn open_conic_params(c: &Curve, ps: P3, pe: P3) -> R<(f64, f64)> { + match (planar_param(c, v3(ps)), planar_param(c, v3(pe))) { + (Some(a), Some(b)) => Ok((a, b)), + _ => Err("no parameter on the conic".into()), + } +} diff --git a/crates/step/src/topology_tests.rs b/crates/step/src/topology_tests.rs index 702815d17d..838f584d20 100644 --- a/crates/step/src/topology_tests.rs +++ b/crates/step/src/topology_tests.rs @@ -687,9 +687,39 @@ fn vertex_parameters_on_lines_conics_and_splines() { let p = [4.0 * t.cos(), t.sin(), 0.0]; let m = map_edge(&ellipse, [4.0, 0.0, 0.0], p, true, false, false, tol).unwrap(); assert!(m.interval[0].abs() < 1e-15 && (m.interval[1] - t).abs() < 1e-14, "{:?}", m.interval); + // A hyperbola edge: rational quadratic pieces on the STEP parameters of its ends, through + // the hyperbola's points; a parabola: one Bézier whose parameter is STEP's. + let (a, b) = (3.0, 2.0); + let hyperbola = Curve::Hyperbola { position, semi_axis: a, semi_imag_axis: b }; + let at = |t: f64| [a * t.cosh(), b * t.sinh(), 0.0]; + let m = map_edge(&hyperbola, at(1.5), at(-2.5), true, false, false, tol).unwrap(); + assert_eq!(m.interval, [-2.5, 1.5]); + assert!(!m.with_step && m.sense_mismatch); + let c = CurveEval::new(&m.curve).unwrap(); + for k in 0..=40 { + let p = c.point(-2.5 + 4.0 * f64::from(k) / 40.0).unwrap(); + assert!(((p.x / a).powi(2) - (p.y / b).powi(2) - 1.0).abs() < 1e-12 && p.x > 0.0, "{p:?}"); + } + for t in [-2.5, -0.5, 0.5, 1.5] { + assert!(c.point(t).unwrap().dist(v3(at(t))) < 1e-12, "{t}"); + } + let parabola = Curve::Parabola { position, focal_dist: 0.5 }; + let at = |t: f64| [0.5 * t * t, t, 0.0]; + let m = map_edge(¶bola, at(-1.0), at(3.0), true, false, false, tol).unwrap(); + assert_eq!(m.interval, [-1.0, 3.0]); + let c = CurveEval::new(&m.curve).unwrap(); + for k in 0..=8 { + let t = -1.0 + 0.5 * f64::from(k); + assert!(c.point(t).unwrap().dist(v3(at(t))) < 1e-12, "{t}"); + } // Unsupported curves. - let hyperbola = Curve::Hyperbola { position, semi_axis: 1.0, semi_imag_axis: 1.0 }; - assert!(map_edge(&hyperbola, s, e, true, false, false, tol).unwrap_err().contains("not supported yet: HYPERBOLA")); + let offset = Curve::OffsetCurve3d { + basis: std::sync::Arc::new(line), + distance: 1.0, + self_intersect: crate::geom::Logical::False, + ref_direction: [0.0, 0.0, 1.0], + }; + assert!(map_edge(&offset, s, e, true, false, false, tol).unwrap_err().contains("not supported yet: OFFSET_CURVE_3D")); } use std::f64::consts::{FRAC_PI_2, PI}; @@ -865,3 +895,310 @@ fn corpus_solids_read_into_schema_4() { } } } + +/// One STEP surface (entity #10 of `data`, after the common placements) sampled over a box of +/// its own parameters: the mapped record passes through every sample, and its natural normal is +/// STEP's (or its opposite where the mapping says so). +#[track_caller] +fn same_surface(data: &str, us: (f64, f64), vs: (f64, f64)) -> BSurface { + const BASE: &str = "#1=CARTESIAN_POINT('',(0.,0.,0.));#2=DIRECTION('',(0.,0.,1.));#3=DIRECTION('',(1.,0.,0.));\ +#4=AXIS2_PLACEMENT_3D('',#1,#2,#3);#5=DIRECTION('',(0.,1.,0.));#6=AXIS1_PLACEMENT('',#1,#2);"; + let ex = parse(&format!("{BASE}{data}\n")); + let s = crate::read::Reader::new(&ex, crate::read::Units::default()).surface(10).unwrap(); + let n = 8; + let grid: Vec<(f64, f64)> = (0..=n) + .flat_map(|i| (0..=n).map(move |j| (us.0 + (us.1 - us.0) * i as f64 / n as f64, vs.0 + (vs.1 - vs.0) * j as f64 / n as f64))) + .collect(); + let ext = swept::Extent { points: grid.iter().map(|&(u, v)| v3(s.point(u, v).unwrap())).collect(), seam: Vec::new() }; + let (record, reversed) = swept::map_surface(&s, &ext).unwrap(); + record.natural(); + let eval = SurfaceEval::new(&record).unwrap(); + for &(u, v) in &grid { + let p = v3(s.point(u, v).unwrap()); + let q = eval.project(p, None, f64::INFINITY).unwrap_or_else(|| panic!("{data}: no parameters at ({u}, {v})")); + let (at, su, sv) = eval.derivs(q[0], q[1]).unwrap(); + assert!(at.dist(p) < 1e-8, "{data}: ({u}, {v}) {p:?} is {} off", at.dist(p)); + let mine = su.cross(sv).normalized().unwrap(); + let theirs = v3(s.normal(u, v).unwrap()); + let agree = mine.dot(theirs) * if reversed { -1.0 } else { 1.0 }; + assert!(agree > 0.999, "{data}: ({u}, {v}) normal {mine:?} vs {theirs:?} (reversed {reversed})"); + } + record +} + +#[test] +fn swept_and_offset_surfaces_are_the_files_own() { + // Revolutions: a slanted line (a cone's frustum), a ball's half circle, an ellipse off the + // axis, a B-spline profile; a circle in a plane through the axis is the ring torus it is. + let line = "#7=CARTESIAN_POINT('',(5.,0.,0.));#8=DIRECTION('',(-2.,0.,4.));#9=VECTOR('',#8,4.47213595499958);#11=LINE('',#7,#9);"; + let r = same_surface(&format!("{line}#10=SURFACE_OF_REVOLUTION('',#11,#6);"), (0.0, 6.0), (0.0, 1.0)); + assert!(matches!(r, BSurface::PlacedRevolution { forward: false, .. })); + // A circle across the axis covers the ball twice, once each way round: its trims say which + // half the surface is. + for (half, vs) in [ + ("(PARAMETER_VALUE(0.)),(PARAMETER_VALUE(3.14159265358979))", (0.2, 3.0)), + ("(PARAMETER_VALUE(3.14159265358979)),(PARAMETER_VALUE(0.))", (3.3, 6.0)), + ] { + let ball = format!( + "#7=AXIS2_PLACEMENT_3D('',#1,#5,#2);#12=CIRCLE('',#7,5.);#11=TRIMMED_CURVE('',#12,{half},.T.,.PARAMETER.);#10=SURFACE_OF_REVOLUTION('',#11,#6);" + ); + assert!(matches!(same_surface(&ball, (0.0, 6.0), vs), BSurface::PlacedRevolution { .. })); + } + let tilted = + "#7=CARTESIAN_POINT('',(10.,0.,1.));#8=AXIS2_PLACEMENT_3D('',#7,#5,#3);#11=ELLIPSE('',#8,3.,2.);#10=SURFACE_OF_REVOLUTION('',#11,#6);"; + same_surface(tilted, (-1.0, 2.0), (-1.0, 1.5)); + let spline = "#7=CARTESIAN_POINT('',(4.,0.,0.));#8=CARTESIAN_POINT('',(6.,0.,2.));#9=CARTESIAN_POINT('',(5.,0.,5.));\ +#11=B_SPLINE_CURVE_WITH_KNOTS('',2,(#7,#8,#9),.UNSPECIFIED.,.F.,.F.,(3,3),(0.,1.),.UNSPECIFIED.);#10=SURFACE_OF_REVOLUTION('',#11,#6);"; + same_surface(spline, (0.0, 6.0), (0.05, 0.95)); + for axis in ["#5", "#12"] { + let torus = format!( + "#12=DIRECTION('',(0.,-1.,0.));#7=CARTESIAN_POINT('',(10.,0.,0.));#8=AXIS2_PLACEMENT_3D('',#7,{axis},#3);#11=CIRCLE('',#8,2.);#10=SURFACE_OF_REVOLUTION('',#11,#6);" + ); + assert!(matches!(same_surface(&torus, (0.0, 6.0), (0.0, 6.0)), BSurface::Torus { .. })); + } + // Extrusions: an ellipse askew, a line (a plane strip), a B-spline, a parabola, a hyperbola; + // a circle along its axis either way is a cylinder. + let askew = "#7=ELLIPSE('',#4,9.,4.);#8=DIRECTION('',(0.5,0.25,1.));#9=VECTOR('',#8,2.);#10=SURFACE_OF_LINEAR_EXTRUSION('',#7,#9);"; + assert!(matches!(same_surface(askew, (-1.0, 1.0), (0.0, 5.0)), BSurface::RationalSpline { .. })); + let strip = "#7=VECTOR('',#3,2.);#8=LINE('',#1,#7);#9=VECTOR('',#5,1.);#10=SURFACE_OF_LINEAR_EXTRUSION('',#8,#9);"; + assert!(matches!(same_surface(strip, (-1.0, 3.0), (0.0, 5.0)), BSurface::Spline { .. })); + let spline = "#7=CARTESIAN_POINT('',(4.,0.,0.));#8=CARTESIAN_POINT('',(6.,1.,0.));#9=CARTESIAN_POINT('',(5.,3.,0.));\ +#11=B_SPLINE_CURVE_WITH_KNOTS('',2,(#7,#8,#9),.UNSPECIFIED.,.F.,.F.,(3,3),(0.,1.),.UNSPECIFIED.);#12=VECTOR('',#2,3.);#10=SURFACE_OF_LINEAR_EXTRUSION('',#11,#12);"; + same_surface(spline, (0.05, 0.95), (-1.0, 1.0)); + let parabola = "#7=PARABOLA('',#4,0.5);#8=VECTOR('',#2,1.);#10=SURFACE_OF_LINEAR_EXTRUSION('',#7,#8);"; + same_surface(parabola, (-2.0, 3.0), (0.0, 4.0)); + let hyperbola = "#7=HYPERBOLA('',#4,3.,2.);#8=VECTOR('',#2,1.);#10=SURFACE_OF_LINEAR_EXTRUSION('',#7,#8);"; + same_surface(hyperbola, (-2.0, 1.5), (0.0, 4.0)); + for (dir, reversed) in [("#2", false), ("#7", true)] { + let tube = format!("#7=DIRECTION('',(0.,0.,-1.));#8=CIRCLE('',#4,5.);#9=VECTOR('',{dir},3.);#10=SURFACE_OF_LINEAR_EXTRUSION('',#8,#9);"); + assert!(matches!(same_surface(&tube, (0.0, 6.0), (0.0, 2.0)), BSurface::Cylinder { radius: 5.0, .. })); + // Offset along STEP's normal: outward one way, inward the other. + let off = tube.replace("#10=SURFACE", "#11=SURFACE").replace("#10=", "#11=") + "#10=OFFSET_SURFACE('',#11,1.,.F.);"; + let want = if reversed { 4.0 } else { 6.0 }; + assert!(matches!(same_surface(&off, (0.0, 6.0), (0.0, 2.0)), BSurface::Cylinder { radius, .. } if (radius - want).abs() < 1e-12), "{dir}"); + } + // Offsets of each analytic kind, nested ones too. + for (basis, d) in [ + ("PLANE('',#4)", 2.0), + ("CYLINDRICAL_SURFACE('',#4,5.)", -2.0), + ("CONICAL_SURFACE('',#4,5.,0.3)", 1.5), + ("CONICAL_SURFACE('',#4,5.,-0.3)", -1.0), + ("SPHERICAL_SURFACE('',#4,5.)", 1.0), + ("TOROIDAL_SURFACE('',#4,10.,2.)", 1.0), + ] { + let one = format!("#11={basis};#10=OFFSET_SURFACE('',#11,{},.F.);", real(d)); + same_surface(&one, (0.1, 2.0), (0.1, 1.2)); + let two = format!("#12={basis};#11=OFFSET_SURFACE('',#12,{},.F.);#10=OFFSET_SURFACE('',#11,0.5,.U.);", real(d)); + same_surface(&two, (0.1, 2.0), (0.1, 1.2)); + } + // Refused: offsets that collapse or of other kinds, an extrusion in its profile's plane. + let ext = swept::Extent { points: vec![Vec3::new(1.0, 0.0, 0.0), Vec3::new(0.0, 1.0, 1.0)], seam: Vec::new() }; + for (data, why) in [ + ("#11=CYLINDRICAL_SURFACE('',#4,5.);#10=OFFSET_SURFACE('',#11,-5.,.F.);", "collapses"), + ("#11=TOROIDAL_SURFACE('',#4,10.,2.);#10=OFFSET_SURFACE('',#11,9.,.F.);", "meets its axis"), + ("#11=CYLINDRICAL_SURFACE('',#4,5.);#12=OFFSET_SURFACE('',#11,-6.,.F.);#10=OFFSET_SURFACE('',#12,6.,.F.);", "collapses"), + ("#7=ELLIPSE('',#4,9.,4.);#8=VECTOR('',#3,1.);#10=SURFACE_OF_LINEAR_EXTRUSION('',#7,#8);", "profile's plane"), + ( + "#7=ELLIPSE('',#4,9.,4.);#8=VECTOR('',#2,1.);#11=SURFACE_OF_LINEAR_EXTRUSION('',#7,#8);#10=OFFSET_SURFACE('',#11,1.,.F.);", + "not supported yet", + ), + ] { + let exch = parse(&format!( + "#1=CARTESIAN_POINT('',(0.,0.,0.));#2=DIRECTION('',(0.,0.,1.));#3=DIRECTION('',(1.,0.,0.));#4=AXIS2_PLACEMENT_3D('',#1,#2,#3);{data}\n" + )); + let s = crate::read::Reader::new(&exch, crate::read::Units::default()).surface(10).unwrap(); + let e = swept::map_surface(&s, &ext).unwrap_err(); + assert!(e.contains(why), "{data}: {e}"); + } +} + +/// A frustum (radius 5 at z = 0 to radius 3 at z = 4) whose side is a SURFACE_OF_REVOLUTION of +/// its slanted profile line, with the profile as a seam edge or without one, about +z or −z +/// (STEP's normal ∂angle × ∂profile then points in, and the face says same_sense .F.). +fn frustum(seam: bool, down: bool) -> p21::Exchange { + let mut w = W::default(); + let (pb, pt) = ([5.0, 0.0, 0.0], [3.0, 0.0, 4.0]); + let (vb, vt) = (w.vertex(pb), w.vertex(pt)); + let kb = w.circle([0.0; 3], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0], 5.0); + let kt = w.circle([0.0, 0.0, 4.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0], 3.0); + let cb = w.add(format!("EDGE_CURVE('',#{vb},#{vb},#{kb},.T.)")); + let ct = w.add(format!("EDGE_CURVE('',#{vt},#{vt},#{kt},.T.)")); + let (o, d) = (w.pt(pb), w.dir([-2.0, 0.0, 4.0])); + let vec = w.add(format!("VECTOR('',#{d},{})", real(20f64.sqrt()))); + let profile = w.add(format!("LINE('',#{o},#{vec})")); + let origin = w.pt([0.0; 3]); + let z = w.dir([0.0, 0.0, if down { -1.0 } else { 1.0 }]); + let axis = w.add(format!("AXIS1_PLACEMENT('',#{origin},#{z})")); + let surface = w.add(format!("SURFACE_OF_REVOLUTION('',#{profile},#{axis})")); + let bounds = if seam { + let s = w.add(format!("EDGE_CURVE('',#{vb},#{vt},#{profile},.T.)")); + let ring = [(cb, true), (s, true), (ct, false), (s, false)]; + let oes: Vec = ring.iter().map(|&(e, o)| w.oe(e, o)).collect(); + vec![w.bound(&oes, true, true)] + } else { + let (x, y) = (w.oe(cb, true), w.oe(ct, false)); + vec![w.bound(&[x], true, true), w.bound(&[y], false, true)] + }; + let side = w.face(&bounds, surface, !down); + let x = w.oe(cb, false); + let bb = w.bound(&[x], true, true); + let plane = w.plane([0.0; 3], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]); + let bottom = w.face(&[bb], plane, false); + let x = w.oe(ct, true); + let tb = w.bound(&[x], true, true); + let plane = w.plane([0.0, 0.0, 4.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]); + let top = w.face(&[tb], plane, true); + let shell = closed_shell(&mut w, &[side, bottom, top]); + let solid = w.add(format!("MANIFOLD_SOLID_BREP('frustum',#{shell})")); + w.file(&[solid]) +} + +#[test] +fn surfaces_of_revolution_read_as_placed_revolutions() { + for seam in [true, false] { + for down in [false, true] { + let s = one_solid(&frustum(seam, down)); + assert_eq!(counts(&s.brep), [2, 1, 3, 3, 2], "seam {seam} down {down}"); + assert!(s.brep.faces.iter().any(|f| matches!(f.surface, BSurface::PlacedRevolution { forward: false, .. }))); + assert_eq!(s.notes.iter().any(|n| n.contains("seam edges added")), !seam, "{:?}", s.notes); + faces_away(&s.brep, 0, |p| [0.0, 0.0, p[2].clamp(1.0, 3.0)]); + } + } + + // A ball: a half circle (trimmed from the north pole to the south through +x) revolved, + // bounded by its seam meridian alone. STEP's normal there points in: same_sense .F. + let mut w = W::default(); + let (vn, vs) = (w.vertex([0.0, 0.0, 5.0]), w.vertex([0.0, 0.0, -5.0])); + let k = w.circle([0.0; 3], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0], 5.0); + let half = w.add(format!("TRIMMED_CURVE('',#{k},(PARAMETER_VALUE(0.)),(PARAMETER_VALUE({})),.T.,.PARAMETER.)", real(PI))); + let e = w.add(format!("EDGE_CURVE('',#{vn},#{vs},#{k},.T.)")); + let (o1, o2) = (w.oe(e, true), w.oe(e, false)); + let bound = w.bound(&[o1, o2], true, true); + let (origin, z) = (w.pt([0.0; 3]), w.dir([0.0, 0.0, 1.0])); + let axis = w.add(format!("AXIS1_PLACEMENT('',#{origin},#{z})")); + let surface = w.add(format!("SURFACE_OF_REVOLUTION('',#{half},#{axis})")); + let face = w.face(&[bound], surface, false); + let shell = closed_shell(&mut w, &[face]); + let solid = w.add(format!("MANIFOLD_SOLID_BREP('ball',#{shell})")); + let s = one_solid(&w.file(&[solid])); + assert_eq!(counts(&s.brep), [2, 1, 1, 1, 2]); + faces_away(&s.brep, 0, |_| [0.0; 3]); +} + +#[test] +fn surfaces_of_extrusion_read_as_splines_or_cylinders() { + // A half elliptic prism: the ellipse (6 by 3) from −y round +x to +y, extruded 5 up, cut + // by the plane x = 0. + let mut w = W::default(); + let p = |y: f64, z: f64| [0.0, y, z]; + let v: Vec<(u64, V)> = [p(-3.0, 0.0), p(3.0, 0.0), p(-3.0, 5.0), p(3.0, 5.0)].iter().map(|q| (w.vertex(*q), *q)).collect(); + let ellipse = |w: &mut W, z: f64| { + let a = w.axis([0.0, 0.0, z], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]); + w.add(format!("ELLIPSE('',#{a},6.,3.)")) + }; + let (eb, et) = (ellipse(&mut w, 0.0), ellipse(&mut w, 5.0)); + let arc_b = w.add(format!("EDGE_CURVE('',#{},#{},#{eb},.T.)", v[0].0, v[1].0)); + let arc_t = w.add(format!("EDGE_CURVE('',#{},#{},#{et},.T.)", v[2].0, v[3].0)); + let up_m = w.line_edge(v[0], v[2], true, 1.0); + let up_p = w.line_edge(v[1], v[3], true, 1.0); + let low = w.line_edge(v[1], v[0], true, 1.0); + let high = w.line_edge(v[2], v[3], true, 1.0); + let z = w.dir([0.0, 0.0, 1.0]); + let vec = w.add(format!("VECTOR('',#{z},2.)")); + let surface = w.add(format!("SURFACE_OF_LINEAR_EXTRUSION('',#{eb},#{vec})")); + let face = |w: &mut W, ring: &[(u64, bool)], surface: u64, sense: bool| { + let oes: Vec = ring.iter().map(|&(e, o)| w.oe(e, o)).collect(); + let b = w.bound(&oes, true, true); + w.face(&[b], surface, sense) + }; + let curved = face(&mut w, &[(arc_b, true), (up_p, true), (arc_t, false), (up_m, false)], surface, true); + let flat_plane = w.plane([0.0; 3], [-1.0, 0.0, 0.0], [0.0, 1.0, 0.0]); + let flat = face(&mut w, &[(low, true), (up_m, true), (high, true), (up_p, false)], flat_plane, true); + let bottom_plane = w.plane([0.0; 3], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]); + let bottom = face(&mut w, &[(arc_b, false), (low, false)], bottom_plane, false); + let top_plane = w.plane([0.0, 0.0, 5.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]); + let top = face(&mut w, &[(arc_t, true), (high, false)], top_plane, true); + let shell = closed_shell(&mut w, &[curved, flat, bottom, top]); + let solid = w.add(format!("MANIFOLD_SOLID_BREP('half',#{shell})")); + let s = one_solid(&w.file(&[solid])); + assert_eq!(counts(&s.brep), [2, 1, 4, 6, 4]); + assert!(s.brep.faces.iter().any(|f| matches!(f.surface, BSurface::RationalSpline { .. }))); + faces_away(&s.brep, 0, |q| [q[0].clamp(0.5, 2.0), q[1].clamp(-1.0, 1.0), q[2].clamp(1.0, 4.0)]); + + // A cylinder whose side is its bottom circle extruded up (or down, STEP's normal then + // pointing in), or a narrower cylinder offset out by 1. + for variant in ["up", "down", "offset"] { + let mut w = W::default(); + let c = cylinder(&mut w, [1.0, 2.0], 5.0, [0.0, 10.0]); + let side_face = side(&mut w, &c, true); + let bottom_bound = { + let ob = w.oe(c.circles[0], false); + w.bound(&[ob], true, true) + }; + let bottom_plane = w.plane([1.0, 2.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]); + let bottom = w.face(&[bottom_bound], bottom_plane, false); + let top_bound = { + let ot = w.oe(c.circles[1], true); + w.bound(&[ot], true, true) + }; + let top_plane = w.plane([1.0, 2.0, 10.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]); + let top = w.face(&[top_bound], top_plane, true); + let k = w.circle([1.0, 2.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0], 5.0); + let d = w.dir([0.0, 0.0, if variant == "down" { -1.0 } else { 1.0 }]); + let vec = w.add(format!("VECTOR('',#{d},1.)")); + let a = w.axis([1.0, 2.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]); + let narrow = w.add(format!("CYLINDRICAL_SURFACE('',#{a},4.)")); + let shell = closed_shell(&mut w, &[bottom, side_face, top]); + let solid = w.add(format!("MANIFOLD_SOLID_BREP('',#{shell})")); + let old = format!("#{}=CYLINDRICAL_SURFACE", c.surface); + let text = w.s.lines().map(|l| { + if l.starts_with(&old) { + match variant { + "offset" => format!("#{}=OFFSET_SURFACE('',#{narrow},1.,.F.);", c.surface), + _ => format!("#{}=SURFACE_OF_LINEAR_EXTRUSION('',#{k},#{vec});", c.surface), + } + } else if variant == "down" && l.starts_with(&format!("#{side_face}=")) { + l.replace(",.T.);", ",.F.);") + } else { + l.to_string() + } + }); + w.s = text.collect::>().join("\n") + "\n"; + assert!(!w.s.contains(&old) && w.s.contains(if variant == "offset" { "OFFSET_SURFACE" } else { "LINEAR_EXTRUSION" })); + assert_eq!(w.s.contains(&format!("#{side_face}=ADVANCED_FACE('',(#{}),#{},.F.)", side_face - 1, c.surface)), variant == "down"); + let s = one_solid(&w.file(&[solid])); + assert_eq!(counts(&s.brep), [2, 1, 3, 3, 2], "{variant}"); + assert!(s.brep.faces.iter().any(|f| matches!(f.surface, BSurface::Cylinder { radius, .. } if (radius - 5.0).abs() < 1e-12)), "{variant}"); + faces_away(&s.brep, 0, |q| [1.0, 2.0, q[2].clamp(1.0, 9.0)]); + } +} + +#[test] +fn parabola_and_hyperbola_edges_bound_faces() { + for kind in ["PARABOLA('',#{a},0.5)", "HYPERBOLA('',#{a},1.,1.)"] { + let mut w = W::default(); + let at = |t: f64| if kind.starts_with("PARABOLA") { [0.5 * t * t, t, 0.0] } else { [t.cosh(), t.sinh(), 0.0] }; + let (p1, p2) = (at(-1.5), at(1.5)); + let (v1, v2) = (w.vertex(p1), w.vertex(p2)); + let a = w.axis([0.0; 3], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]); + let conic = w.add(kind.replace("{a}", &a.to_string())); + let e1 = w.add(format!("EDGE_CURVE('',#{v1},#{v2},#{conic},.T.)")); + let e2 = w.line_edge((v2, p2), (v1, p1), true, 1.0); + let (o1, o2) = (w.oe(e1, true), w.oe(e2, true)); + let bound = w.bound(&[o1, o2], true, true); + let plane = w.plane([0.0; 3], [0.0, 0.0, -1.0], [1.0, 0.0, 0.0]); + let face = w.face(&[bound], plane, true); + let shell = w.add(format!("OPEN_SHELL('',(#{face}))")); + let model = w.add(format!("SHELL_BASED_SURFACE_MODEL('',(#{shell}))")); + let s = one_solid(&w.file(&[model])); + assert_eq!(counts(&s.brep), [1, 1, 1, 2, 2], "{kind}"); + let e = s.brep.edges.iter().find(|e| !matches!(e.curve, BCurve::Line { .. })).unwrap(); + assert_eq!(e.interval, [-1.5, 1.5]); + let c = CurveEval::new(&e.curve).unwrap(); + for k in 0..=12 { + let p = c.point(-1.5 + 0.25 * f64::from(k)).unwrap(); + let on = if kind.starts_with("PARABOLA") { p.y * p.y - 2.0 * p.x } else { p.x * p.x - p.y * p.y - 1.0 }; + assert!(on.abs() < 1e-12, "{kind}: {p:?}"); + } + } +}