diff --git a/crates/difftest/src/units/curve_ellipse_p1.rs b/crates/difftest/src/units/curve_ellipse_p1.rs index 9f64cf6bac..7078bb2d77 100644 --- a/crates/difftest/src/units/curve_ellipse_p1.rs +++ b/crates/difftest/src/units/curve_ellipse_p1.rs @@ -33,11 +33,13 @@ use solidcore::geom::vector3d::Vec3; use solidcore::kernel::Kernel; use solidcore::nurbs::b_spline_curve::BSplineCurve; use solidcore::nurbs::b_spline_curve_n_lib::TrimMode; +use solidcore::nurbs::crv_basic::Curve as NurbCurve; +use solidcore::nurbs::math::{CPoint, NOW}; use solidcore::topology::database_io::{DatabaseReader, FileType}; use solidcore::SolidError; use solidcore_ffi::units::curve_ellipse_p1::{self as ffi}; -use crate::case::CaseContext; +use crate::case::{CaseContext, CaseSource}; use crate::runner::default_cases; use crate::{Case, HarnessError, Output, Rng, UnitDiff}; @@ -162,6 +164,22 @@ enum OtherJ { BsplineOfA, } +/// The curve tested by `SmCircle::IsNurbCurveCircle`. +#[derive(Debug, Clone, Serialize, Deserialize)] +#[serde(rename_all = "snake_case")] +enum SourceJ { + /// A itself. + A, + /// A plain B-spline with A's NURB. + BsplineOfA, + /// A raw 3D NURB: homogeneous `x y z w` points (`w` = `NOW` when non-rational). + Raw { + degree: usize, + cpts: Vec<[f64; 4]>, + knots: Vec, + }, +} + fn full() -> [f64; 2] { [0.0, 1.0] } @@ -239,6 +257,11 @@ enum OpJ { post: bool, }, IsArc, + IsNurbCircle { + source: SourceJ, + #[serde(default, skip_serializing_if = "Option::is_none")] + scale: Option, + }, } #[derive(Debug, Clone, Serialize, Deserialize)] @@ -422,6 +445,22 @@ fn ffi_op(op: &OpJ) -> ffi::Op { }, OpJ::Notify { post } => F::Notify { post: *post }, OpJ::IsArc => F::IsArc, + OpJ::IsNurbCircle { source, scale } => F::IsNurbCircle { + source: match source { + SourceJ::A => ffi::NurbSource::A, + SourceJ::BsplineOfA => ffi::NurbSource::BSplineOfA, + SourceJ::Raw { + degree, + cpts, + knots, + } => ffi::NurbSource::Raw { + degree: *degree, + cpts: cpts.clone(), + knots: knots.clone(), + }, + }, + scale: *scale, + }, } } @@ -975,6 +1014,36 @@ fn run_op(k: &Kernel, a: &mut Ellipse, op: &OpJ, o: &mut Vec) -> Result<(), o.extend([r, s0, s1]); } } + OpJ::IsNurbCircle { source, scale } => { + let test = match source { + SourceJ::A => Curve::from(a.clone()), + SourceJ::BsplineOfA => match a.base().gw_nurb() { + Some(n) => Curve::from(BSplineCurve::from_gw_nurb(3, n)), + None => return Err(harness("A without NURB").into()), + }, + SourceJ::Raw { + degree, + cpts, + knots, + } => { + let pw = cpts + .iter() + .map(|&[x, y, z, w]| CPoint::new(x, y, z, w)) + .collect(); + let cur = NurbCurve::new(pw, *degree, knots.clone()) + .map_err(|e| harness(format!("raw NURB: {e}")))?; + Curve::from(BSplineCurve::from_nurb(cur, 3, false)) + } + }; + if let Some(got) = st(o, Ellipse::is_nurb_curve_circle(k, &test, *scale))? { + o.push(b(got.is_some())); + if let Some(c) = got { + o.push(b(c.kind_of(CIRCLE_TYPE))); + o.push(c.radius()); + put_state(o, &c); + } + } + } } Ok(()) } @@ -1001,6 +1070,7 @@ fn op_code(op: &OpJ) -> f64 { OpJ::EditParam { .. } => 18.0, OpJ::Notify { .. } => 19.0, OpJ::IsArc => 20.0, + OpJ::IsNurbCircle { .. } => 21.0, } } @@ -1484,6 +1554,442 @@ fn gen_case(rng: &mut Rng) -> Option { Some(Inputs { build: bj, ops }) } +// ---- SmCircle::IsNurbCurveCircle cases ---------------------------------------------- + +/// Generated `IsNurbCurveCircle` cases (inline, fixed seed, independent of `--seeds` +/// so that the main generated set keeps its seeds). +const NURB_CIRCLE_GENERATED: usize = 400; +/// Seed of the generated `IsNurbCurveCircle` cases. +const NURB_CIRCLE_SEED: u64 = 0x1432_0c1c; + +fn inc(source: SourceJ, scale: Option) -> OpJ { + OpJ::IsNurbCircle { source, scale } +} + +fn unit_frame() -> [[f64; 3]; 3] { + [[0.0; 3], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]] +} + +/// The unit canonical circle: the host of raw-NURB ops (and the issue's repro). +fn host() -> BuildJ { + BuildJ::CircleCanonical { + frame: unit_frame(), + r: 1.0, + interval: None, + inside_out: None, + } +} + +fn circle_j(c: [f64; 3], x: [f64; 3], y: [f64; 3], domain: [f64; 2], r: f64) -> BuildJ { + BuildJ::Circle(CircleJ { + c, + x, + y, + domain, + r, + dim: 3, + inside_out: false, + nurb_from: None, + }) +} + +fn ellipse_j(c: [f64; 3], domain: [f64; 2], radii: [f64; 2]) -> BuildJ { + BuildJ::Ellipse(EllipseJ { + c, + x: [1.0, 0.0, 0.0], + y: [0.0, 1.0, 0.0], + domain, + radii, + dim: 3, + inside_out: false, + nurb_from: None, + }) +} + +fn add3(a: [f64; 3], b: [f64; 3]) -> [f64; 3] { + [a[0] + b[0], a[1] + b[1], a[2] + b[2]] +} + +fn scale3(a: [f64; 3], s: f64) -> [f64; 3] { + [a[0] * s, a[1] * s, a[2] * s] +} + +/// An orthonormal pair spanning a plane (from two directions, Gram-Schmidt). +fn plane(a: [f64; 3], b: [f64; 3]) -> ([f64; 3], [f64; 3]) { + let u = v3(a).unitize_or_z(); + let bv = v3(b); + let mut w = Vec3::new( + bv.x - u.x * u.dot(bv), + bv.y - u.y * u.dot(bv), + bv.z - u.z * u.dot(bv), + ); + if w.length() < 1e-3 { + w = u.cross(Vec3::new(0.3, -0.5, 0.8)); + } + let w = w.unitize_or_z(); + ([u.x, u.y, u.z], [w.x, w.y, w.z]) +} + +/// The rational quadratic NURB of the arc of radius `r` about `c` in the plane of +/// the orthonormal `u`, `v` from angle `a0` to `a1` (radians) in `spans` equal spans +/// (each below half a turn), weights times `wscale`, over `[t0, t0 + len]`. +#[allow(clippy::too_many_arguments)] +fn raw_arc( + c: [f64; 3], + (u, v): ([f64; 3], [f64; 3]), + r: f64, + [a0, a1]: [f64; 2], + spans: usize, + wscale: f64, + [t0, len]: [f64; 2], +) -> SourceJ { + let d = (a1 - a0) / spans as f64; + let pt = |a: f64, rr: f64| add3(c, add3(scale3(u, rr * a.cos()), scale3(v, rr * a.sin()))); + let hw = |p: [f64; 3], w: f64| { + let w = w * wscale; + [p[0] * w, p[1] * w, p[2] * w, w] + }; + let w1 = (0.5 * d).cos(); + let mut cpts = vec![hw(pt(a0, r), 1.0)]; + let mut knots = vec![t0; 3]; + for i in 0..spans { + let s0 = a0 + d * i as f64; + cpts.push(hw(pt(s0 + 0.5 * d, r / w1), w1)); + cpts.push(hw(pt(s0 + d, r), 1.0)); + let t = t0 + len * (i + 1) as f64 / spans as f64; + if i + 1 < spans { + knots.extend([t, t]); + } + } + knots.extend([t0 + len; 3]); + SourceJ::Raw { + degree: 2, + cpts, + knots, + } +} + +/// `raw_arc` in the XY plane about the origin over `[0, 1]`. +fn raw_xy(r: f64, angles: [f64; 2], spans: usize, wscale: f64) -> SourceJ { + raw_arc( + [0.0; 3], + ([1.0, 0.0, 0.0], [0.0, 1.0, 0.0]), + r, + angles, + spans, + wscale, + [0.0, 1.0], + ) +} + +/// `src` with control point `i` changed by `f` (homogeneous). +fn tweak(src: SourceJ, i: usize, f: impl Fn([f64; 4]) -> [f64; 4]) -> SourceJ { + match src { + SourceJ::Raw { + degree, + mut cpts, + knots, + } => { + if let Some(p) = cpts.get_mut(i) { + *p = f(*p); + } + SourceJ::Raw { + degree, + cpts, + knots, + } + } + other => other, + } +} + +/// The weight of control point `i` times `1 + e` (the Euclidean point kept). +fn tweak_weight(src: SourceJ, i: usize, e: f64) -> SourceJ { + tweak(src, i, |[x, y, z, w]| { + let s = 1.0 + e; + [x * s, y * s, z * s, w * s] + }) +} + +/// Control point `i` moved by `d` (Euclidean, the weight kept). +fn tweak_point(src: SourceJ, i: usize, d: [f64; 3]) -> SourceJ { + tweak(src, i, |[x, y, z, w]| { + [x + d[0] * w, y + d[1] * w, z + d[2] * w, w] + }) +} + +fn nurb_circle_fixed() -> Vec<(String, Inputs)> { + use std::f64::consts::{FRAC_PI_2, PI}; + let mut v: Vec<(String, Inputs)> = Vec::new(); + let mut push = |name: &str, build: BuildJ, ops: Vec| { + v.push((format!("inc_{name}"), Inputs { build, ops })); + }; + let both = |scale: Option| { + vec![ + inc(SourceJ::BsplineOfA, scale), + inc(SourceJ::A, scale), + ] + }; + // #1432: the NURB form of the full unit canonical circle as a plain B-spline + push("issue_unit_circle", host(), both(None)); + push( + "unit_circle_scales", + host(), + [0.0, 1e-3, 10.0, 1e4, -1.0] + .iter() + .map(|&s| inc(SourceJ::BsplineOfA, Some(s))) + .collect(), + ); + // arcs of every span count + let tilted = plane([0.3, -0.4, 0.8], [0.9, 0.1, -0.2]); + for (i, dom) in [ + [0.0, 90.0], + [30.0, 200.0], + [-90.0, 45.0], + [0.0, 359.9], + [0.0, 1e-3], + [0.0, 180.0], + [-180.0, 180.0], + [100.0, 370.0], + ] + .iter() + .enumerate() + { + push( + &format!("circle_arc_{i}"), + circle_j([1.0, 2.0, 3.0], tilted.0, tilted.1, *dom, 2.5), + both(None), + ); + } + // huge and tiny scales + for (i, (c, r)) in [ + ([1e5, -1e5, 5e4], 1e5), + ([4e5, 0.0, 0.0], 5e5), + ([0.0; 3], 1e-5), + ([1.0, 1.0, 1.0], 1e-7), + ([1e3, 0.0, 0.0], 1e-6), + ([0.0; 3], 1e-9), + ] + .iter() + .enumerate() + { + push( + &format!("circle_scale_{i}"), + circle_j(*c, [1.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 360.0], *r), + both(None), + ); + push( + &format!("circle_scale_arc_{i}"), + circle_j(*c, [0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [10.0, 75.0], *r), + both(Some(1e-2)), + ); + } + // SmEllipse (not SmCircle): equal radii is an arc, unequal is not + push( + "ellipse_equal_radii", + ellipse_j([1.0, 0.0, 0.0], [0.0, 360.0], [2.0, 2.0]), + both(None), + ); + push( + "ellipse_equal_radii_arc", + ellipse_j([1.0, 0.0, 0.0], [10.0, 100.0], [2.0, 2.0]), + both(None), + ); + push( + "ellipse_near_equal_radii", + ellipse_j([0.0; 3], [0.0, 360.0], [2.0, 2.0 * (1.0 + 1e-9)]), + both(None), + ); + push( + "ellipse_unequal_radii", + ellipse_j([0.0; 3], [0.0, 360.0], [2.0, 1.0]), + both(None), + ); + // raw NURBs + let raw = |name: &str, src: SourceJ, scale: Option| -> (String, Inputs) { + ( + format!("inc_raw_{name}"), + Inputs { + build: host(), + ops: vec![inc(src, scale)], + }, + ) + }; + let quarter = raw_xy(1.0, [0.0, FRAC_PI_2], 1, 1.0); + v.push(raw("quarter", quarter.clone(), None)); + v.push(raw("quarter_weights_x2", raw_xy(1.0, [0.0, FRAC_PI_2], 1, 2.0), None)); + v.push(raw("quarter_weights_x0.01", raw_xy(1.0, [0.0, FRAC_PI_2], 1, 0.01), None)); + v.push(raw("semicircle_2span", raw_xy(3.0, [0.0, PI], 2, 1.0), None)); + v.push(raw("full_4span", raw_xy(3.0, [0.0, 2.0 * PI], 4, 1.0), None)); + v.push(raw("full_3span", raw_xy(0.5, [-PI, PI], 3, 1.0), None)); + v.push(raw( + "tilted_120", + raw_arc([5.0, -2.0, 1.0], tilted, 0.75, [0.2, 0.2 + 2.0 * PI / 3.0], 1, 1.0, [2.0, 3.5]), + None, + )); + v.push(raw("wide_170", raw_xy(1.0, [0.0, 170f64.to_radians()], 1, 1.0), None)); + v.push(raw("tiny_angle", raw_xy(1.0, [0.0, 1e-6], 1, 1.0), None)); + for (i, e) in [1e-14, 1e-12, 1e-10, 1e-8, 1e-6, 1e-4].iter().enumerate() { + v.push(raw(&format!("weight_perturbed_{i}"), tweak(quarter.clone(), 1, |[x, y, z, w]| [x, y, z, w * (1.0 + e)]), None)); + v.push(raw(&format!("point_perturbed_{i}"), tweak_point(quarter.clone(), 1, [*e, -*e, 0.0]), None)); + v.push(raw(&format!("point_out_of_plane_{i}"), tweak_point(quarter.clone(), 2, [0.0, 0.0, *e]), None)); + v.push(raw(&format!("end_weight_{i}"), tweak_weight(quarter.clone(), 0, *e), None)); + } + // scaled tolerance at the limit: the same near-circle accepted or not by the scale + for (i, s) in [1e-3, 1.0, 1e3, 1e6].iter().enumerate() { + v.push(raw(&format!("tol_scale_{i}"), tweak_point(quarter.clone(), 1, [1e-9, 0.0, 0.0]), Some(*s))); + } + // not circles + v.push(raw("hyperbola_weight", tweak(quarter.clone(), 1, |[x, y, z, w]| [x / w * 2.0, y / w * 2.0, z / w * 2.0, 2.0]), None)); + v.push(raw("ellipse_weight", tweak(quarter.clone(), 1, |[x, y, z, w]| [x / w * 0.3, y / w * 0.3, z / w * 0.3, 0.3]), None)); + v.push(raw("unit_weights", tweak(quarter.clone(), 1, |[x, y, z, w]| [x / w, y / w, z / w, 1.0]), None)); + v.push(raw( + "non_rational", + SourceJ::Raw { + degree: 2, + cpts: vec![[1.0, 0.0, 0.0, NOW], [1.0, 1.0, 0.0, NOW], [0.0, 1.0, 0.0, NOW]], + knots: vec![0.0, 0.0, 0.0, 1.0, 1.0, 1.0], + }, + None, + )); + v.push(raw( + "rational_cubic", + SourceJ::Raw { + degree: 3, + cpts: vec![ + [1.0, 0.0, 0.0, 1.0], + [0.5, 0.25, 0.0, 0.5], + [0.25, 0.5, 0.0, 0.5], + [0.0, 1.0, 0.0, 1.0], + ], + knots: vec![0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0], + }, + None, + )); + v.push(raw( + "rational_linear", + SourceJ::Raw { + degree: 1, + cpts: vec![[1.0, 0.0, 0.0, 1.0], [0.0, 0.5, 0.0, 0.5]], + knots: vec![0.0, 0.0, 1.0, 1.0], + }, + None, + )); + v.push(raw( + "degenerate_point", + SourceJ::Raw { + degree: 2, + cpts: vec![[1.0, 2.0, 3.0, 1.0], [0.5, 1.0, 1.5, 0.5], [1.0, 2.0, 3.0, 1.0]], + knots: vec![0.0, 0.0, 0.0, 1.0, 1.0, 1.0], + }, + None, + )); + v.push(raw("degenerate_tiny_radius", raw_xy(1e-12, [0.0, FRAC_PI_2], 1, 1.0), None)); + v.push(raw("tiny_radius_1e-8", raw_xy(1e-8, [0.0, FRAC_PI_2], 1, 1.0), None)); + v.push(raw("huge_radius", raw_xy(1e5, [0.0, FRAC_PI_2], 1, 1.0), None)); + v.push(raw( + "huge_offset_small_radius", + raw_arc([5e5, 5e5, 0.0], ([1.0, 0.0, 0.0], [0.0, 1.0, 0.0]), 1e-2, [0.0, 1.0], 1, 1.0, [0.0, 1.0]), + None, + )); + v.push(raw( + "collinear", + SourceJ::Raw { + degree: 2, + cpts: vec![[0.0, 0.0, 0.0, 1.0], [0.5, 0.0, 0.0, 0.5], [2.0, 0.0, 0.0, 1.0]], + knots: vec![0.0, 0.0, 0.0, 1.0, 1.0, 1.0], + }, + None, + )); + v.push(raw( + "s_curve_2span", + tweak_point(raw_xy(1.0, [0.0, PI], 2, 1.0), 3, [0.0, 0.0, 0.5]), + None, + )); + v +} + +/// A generated `IsNurbCurveCircle` case. +fn gen_nurb_circle(rng: &mut Rng) -> Inputs { + use std::f64::consts::PI; + let scale = if rng.chance(0.25) { + Some(*rng.choose(&[0.0, 1e-3, 0.5, 2.0, 1e3, 1e6]).unwrap_or(&1.0)) + } else { + None + }; + let r = 10f64.powf(rng.uniform(-6.0, 5.0)); + let c = scale3(rng.point3(-1.0, 1.0), 10f64.powf(rng.uniform(-3.0, 5.0))); + match rng.below(4) { + 0 => { + let f = gen_frame(rng); + Inputs { + build: circle_j(c, f[1], f[2], gen_domain(rng), r), + ops: vec![inc(SourceJ::BsplineOfA, scale), inc(SourceJ::A, scale)], + } + } + 1 => { + let rr = 2.0 * (1.0 + if rng.chance(0.5) { 10f64.powf(rng.uniform(-14.0, -2.0)) } else { 0.0 }); + Inputs { + build: ellipse_j(c, gen_domain(rng), [2.0, rr]), + ops: vec![inc(SourceJ::A, scale), inc(SourceJ::BsplineOfA, scale)], + } + } + _ => { + let total = rng.uniform(1e-3, 2.0 * PI); + let min_spans = (total / (0.9 * PI)).ceil() as usize; + let spans = min_spans.max(1) + rng.below(2) as usize; + let a0 = rng.uniform(-PI, PI); + let pl = plane(gen_dir(rng), gen_dir(rng)); + let wscale = 10f64.powf(rng.uniform(-1.0, 1.0)); + let t0 = rng.uniform(-10.0, 10.0); + let len = 10f64.powf(rng.uniform(-2.0, 2.0)); + let mut src = raw_arc(c, pl, r, [a0, a0 + total], spans, wscale, [t0, len]); + let n = 2 * spans + 1; + let e = 10f64.powf(rng.uniform(-15.0, -2.0)); + let i = rng.below(n as u64) as usize; + src = match rng.below(4) { + 0 => tweak_weight(src, i, e), + 1 => { + let d = scale3(gen_dir(rng), e * r); + tweak_point(src, i, d) + } + _ => src, + }; + Inputs { + build: host(), + ops: vec![inc(src, scale)], + } + } + } +} + +/// The inline `IsNurbCurveCircle` cases: fixed, then generated (encodable only). +fn nurb_circle_cases(unit: &str) -> Vec { + let mut all = nurb_circle_fixed(); + let mut rng = Rng::new(NURB_CIRCLE_SEED); + let mut made = 0; + let mut tries = 0; + while made < NURB_CIRCLE_GENERATED && tries < 8 * NURB_CIRCLE_GENERATED { + tries += 1; + let i = gen_nurb_circle(&mut rng); + let (b, ops) = ffi_program(&i); + if ffi::encode(&b, &ops).is_ok() { + all.push((format!("incgen{made:04}"), i)); + made += 1; + } + } + all.into_iter() + .filter_map(|(name, i)| { + Some(Case { + unit: unit.to_string(), + name, + source: CaseSource::Inline, + inputs: serde_json::to_value(i).ok()?, + }) + }) + .collect() +} + impl UnitDiff for Runner { fn id(&self) -> &'static str { "curve.ellipse_p1" @@ -1507,7 +2013,9 @@ impl UnitDiff for Runner { seeds: ctx.seeds.max(MIN_GENERATED), ..ctx.clone() }; - default_cases(self, &ctx) + let mut v = default_cases(self, &ctx)?; + v.extend(nurb_circle_cases(self.id())); + Ok(v) } fn run_cpp(&self, case: &Case) -> Result { @@ -1518,3 +2026,56 @@ impl UnitDiff for Runner { rust(&case.parse()?) } } + +#[cfg(test)] +mod tests { + use super::*; + + /// The `IsNurbCurveCircle` flag of the first op of `case` (Rust side). + fn first_flag(case: &Case) -> f64 { + let i: Inputs = case.parse().unwrap(); + let out = rust(&i).unwrap(); + let rec = out.values.get("rec").unwrap(); + // construction status, op code 21, op status 0, flag + assert_eq!(&rec[..3], &[0.0, 21.0, 0.0], "{}", case.name); + rec[3] + } + + /// The inline `IsNurbCurveCircle` cases reach both outcomes (the circle is built + /// for the issue's repro and the clean raw arcs; not for non-circles). + #[test] + fn nurb_circle_cases_reach_both_outcomes() { + let cases = nurb_circle_cases("curve.ellipse_p1"); + let flag = |name: &str| { + first_flag(cases.iter().find(|c| c.name == name).unwrap()) + }; + for yes in [ + "inc_issue_unit_circle", + "inc_raw_quarter", + "inc_raw_quarter_weights_x2", + "inc_raw_full_4span", + "inc_raw_tilted_120", + "inc_circle_arc_1", + "inc_ellipse_equal_radii", + ] { + assert_eq!(flag(yes), 1.0, "{yes}"); + } + for no in [ + "inc_raw_non_rational", + "inc_raw_rational_cubic", + "inc_raw_hyperbola_weight", + "inc_raw_degenerate_point", + "inc_raw_collinear", + "inc_ellipse_unequal_radii", + ] { + assert_eq!(flag(no), 0.0, "{no}"); + } + let gen: Vec = cases + .iter() + .filter(|c| c.name.starts_with("incgen")) + .map(first_flag) + .collect(); + let yes = gen.iter().filter(|&&f| f == 1.0).count(); + assert!(yes > gen.len() / 5 && gen.len() - yes > gen.len() / 10, "{yes} of {}", gen.len()); + } +} diff --git a/crates/solidcore-ffi/src/units/curve_ellipse_p1.rs b/crates/solidcore-ffi/src/units/curve_ellipse_p1.rs index fa1f0c1bcf..dc512efdaa 100644 --- a/crates/solidcore-ffi/src/units/curve_ellipse_p1.rs +++ b/crates/solidcore-ffi/src/units/curve_ellipse_p1.rs @@ -34,6 +34,7 @@ //! | `EditParam` | `EditParameterization(SmExtent1d(d0, d1), notify)` | | //! | `Notify` | `Notify(SM_NO_PRE_EDIT / SM_NO_POST_EDIT, this, NULL, NULL)` | | //! | `IsArc` | `SmCircle::IsArc(5, 1e-6, ..)` (circles only) | `flag O X Y r start end` | +//! | `IsNurbCircle` | `SmCircle::IsNurbCurveCircle(ctx, curve, circle, scale?)` on [`NurbSource`] | `flag` (`TRUE`: `IsKindOf(SmCircle_TYPE)`, `GetRadius`, state of the circle) | //! //! Intervals are normalised parameters `[f0, f1]` of A's natural interval //! (extrapolated outside `[0, 1]`); a NURB parameter is raw or normalised. Solutions are @@ -45,7 +46,7 @@ use std::fmt::Write as _; -use crate::validate::bad; +use crate::validate::{bad, check_knots, NOW}; use crate::{call_with_buffer, FfiError}; /// Largest magnitude of coordinates, lengths, parameters and tolerances. @@ -56,6 +57,12 @@ pub const MAX_ANGLE: f64 = 720.0; pub const MAX_OPS: usize = 64; /// Most types of [`Op::State`]. pub const MAX_TYPES: usize = 64; +/// Most control points of a raw NURB of [`NurbSource::Raw`]. +pub const MAX_RAW_POINTS: usize = 65; +/// Largest degree of a raw NURB of [`NurbSource::Raw`]. +pub const MAX_RAW_DEGREE: usize = 8; +/// Weight range of a rational raw NURB of [`NurbSource::Raw`]. +pub const RAW_WEIGHTS: [f64; 2] = [1e-3, 1e3]; /// Largest last control point index and polar table length of a database record. pub const MAX_DB_INDEX: i64 = 64; @@ -220,6 +227,25 @@ pub enum JoinOther { BSplineOfA, } +/// The curve tested by [`Op::IsNurbCircle`]. +#[derive(Debug, Clone, PartialEq)] +pub enum NurbSource { + /// A itself (an `SmCircle` is copied; an `SmEllipse` dispatches to its overrides). + A, + /// `SmBSplineCurve(3, A's NURB)`. + BSplineOfA, + /// `SmBSplineCurve(gw_CURVE, 3, FALSE)` of a clamped 3D NURB. + Raw { + /// Degree (`1..=MAX_RAW_DEGREE`). + degree: usize, + /// Homogeneous control points `x y z w` (`w` positive in [`RAW_WEIGHTS`] for + /// every point, or `NOW` for every point: non-rational). + cpts: Vec<[f64; 4]>, + /// Clamped knots (`cpts.len() + degree + 1`). + knots: Vec, + }, +} + /// An operation on curve A (see the module docs). #[derive(Debug, Clone, PartialEq)] pub enum Op { @@ -328,6 +354,13 @@ pub enum Op { }, /// `SmCircle::IsArc` (circles only). IsArc, + /// `SmCircle::IsNurbCurveCircle` (static) of a curve. + IsNurbCircle { + /// The tested curve. + source: NurbSource, + /// `dToleranceScale` (`None`: the default `1.0`). + scale: Option, + }, } fn finite(x: f64) -> Result { @@ -544,6 +577,55 @@ impl Enc { } self.num(x) } + /// `n p m cpts knots` of a clamped 3D NURB. + fn raw_nurb( + &mut self, + degree: usize, + cpts: &[[f64; 4]], + knots: &[f64], + ) -> Result<(), FfiError> { + if cpts.len() < 2 || cpts.len() > MAX_RAW_POINTS { + return Err(bad("raw NURB: control point count")); + } + let n = cpts.len() - 1; + if degree == 0 || degree > MAX_RAW_DEGREE || degree > n { + return Err(bad("raw NURB: need 1 <= p <= min(n, MAX_RAW_DEGREE)")); + } + if knots.len() != n + degree + 2 { + return Err(bad("raw NURB: need m = n + p + 1")); + } + if knots.iter().any(|k| !k.is_finite() || k.abs() > MAX_COORD) { + return Err(bad("raw NURB: non-finite or huge knot")); + } + check_knots(knots, degree, "raw NURB")?; + let rational = cpts.first().is_some_and(|p| p[3] != NOW); + for &[x, y, z, w] in cpts { + let w_e = if rational { + if !(RAW_WEIGHTS[0]..=RAW_WEIGHTS[1]).contains(&w) { + return Err(bad("raw NURB: weight outside RAW_WEIGHTS or mixed rationality")); + } + w + } else { + if w != NOW { + return Err(bad("raw NURB: mixed rationality")); + } + 1.0 + }; + for c in [x, y, z] { + if !c.is_finite() || (c / w_e).abs() > MAX_COORD { + return Err(bad("raw NURB: non-finite or huge point")); + } + } + } + self.int(n); + self.int(degree); + self.int(n + degree + 1); + for p in cpts { + self.0.extend_from_slice(p); + } + self.0.extend_from_slice(knots); + Ok(()) + } fn int(&mut self, x: usize) { self.0.push(x as f64); } @@ -873,6 +955,23 @@ impl Enc { self.flag(*post); } Op::IsArc => self.int(20), + Op::IsNurbCircle { source, scale } => { + self.int(21); + match source { + NurbSource::A => self.int(0), + NurbSource::BSplineOfA => self.int(1), + NurbSource::Raw { + degree, + cpts, + knots, + } => { + self.int(2); + self.raw_nurb(*degree, cpts, knots)?; + } + } + self.flag(scale.is_some()); + self.num(scale.unwrap_or(0.0))?; + } } Ok(()) } diff --git a/crates/solidcore/src/curve/ellipse/p1.rs b/crates/solidcore/src/curve/ellipse/p1.rs index 24dee72777..25db7c4dad 100644 --- a/crates/solidcore/src/curve/ellipse/p1.rs +++ b/crates/solidcore/src/curve/ellipse/p1.rs @@ -22,6 +22,7 @@ //! | [`Ellipse::is_circle_class`], [`Ellipse::kind_of`] | `SmCircle_TYPE` objects, `SmCircle::IsKindOf` | //! | [`Ellipse::type_code`], [`Ellipse::type_string`], [`Ellipse::class_string`] (in `p2`), [`CIRCLE_CLASS_TYPE`] | `SmCircle::GetType`, `GetTypeString`, `GetClassString`, `GetClassType` | //! | [`Ellipse::is_arc`] | `SmCircle::IsArc` | +//! | [`Ellipse::is_nurb_curve_circle`] | `SmCircle::IsNurbCurveCircle` | //! | [`Ellipse::create_simple_offset`] | `SmCircle::CreateSimpleOffset` | //! | [`Ellipse::assert_valid`] | `SmEllipse::AssertValid`, `SmCircle::AssertValid` | //! | [`Ellipse::read_from_db`] | `SmEllipse::ReadFromDB`, `SmCircle::ReadFromDB` | @@ -88,6 +89,11 @@ //! - **Back-edges.** `SmCurve::IsArc` of a non-circular ellipse (inside //! `SmBSplineCurve::CreateSimpleOffset`) is the base class of unit `curve.curve_p5` //! ([`crate::curve::curve::p5::is_arc_base`]) on the ellipse's NURB. +//! - **`IsNurbCurveCircle`.** The `SmCircle*&` out-param is the returned `Option` +//! (`None` for `FALSE`), the defaulted `dToleranceScale` an `Option` (`None`: `1.0`); +//! the `SmBSplineCurve*` is a [`Curve`] whose kind must derive from `SmBSplineCurve` +//! (B-spline, ellipse, circle, line, parabola, hyperbola), others are +//! [`SolidError::InvalidInput`]. The ignored `EvaluatePoint` status propagates. use core::ops::Div; @@ -136,6 +142,9 @@ const OP_CANONICAL: &str = "curve::ellipse::p1::Ellipse::create_canonical"; const OP_CIRCLE_CANONICAL: &str = "curve::ellipse::p1::Ellipse::create_canonical_circle"; const OP_JOIN: &str = "curve::ellipse::p1::Ellipse::join_with"; const OP_OFFSET: &str = "curve::ellipse::p1::Ellipse::create_simple_offset"; +const OP_IS_NURB_CURVE_CIRCLE: &str = "curve::ellipse::p1::Ellipse::is_nurb_curve_circle"; +/// `ANALYTIC_TOL_SCALE` (`SmConfig.h`). +const ANALYTIC_TOL_SCALE: f64 = 100.0; /// `SmCircle_TYPE` (`CURV_BASE_TYPE + 21`), the class type of a circle. pub const CIRCLE_CLASS_TYPE: i64 = CIRCLE_TYPE; @@ -1014,6 +1023,96 @@ impl Ellipse { ) } + /// The circle a B-spline curve is, when its geometry and representation are those + /// of a circle (`SmCircle::IsNurbCurveCircle`; `tolerance_scale` defaults to `1.0`): + /// a copy of a circle; `None` for a curve that is not of degree 2 and rational, is + /// degenerate (`IsDegenerate` to `scale * ANALYTIC_TOL_SCALE * SM_EFF_ZERO * (1 + + /// |P|max)`, `P` the point at 0.4567 of the natural interval) or is not an arc + /// (`IsArc(5, ..)` to the same tolerance); else the circle of the arc's frame, + /// radius and angular domain, dimension 3, with a copy of the curve's NURB + /// (`pOptNurb`). + /// + /// # Errors + /// + /// [`SolidError::InvalidInput`] for a curve that is not an `SmBSplineCurve` (a + /// Hermite, composite, conic or procedural curve); those of the curve's + /// evaluator, [`Curve::is_degenerate`], [`Curve::is_arc`] and + /// [`Ellipse::new_circle`]. + pub fn is_nurb_curve_circle( + k: &Kernel, + test_curve: &Curve, + tolerance_scale: Option, + ) -> Result> { + let tolerance_scale = tolerance_scale.unwrap_or(1.0); + let line_base; + let test_base: &BSplineCurve = match test_curve.kind() { + CurveKind::BSpline(c) => c, + CurveKind::Ellipse(e) => { + // low work - already a circle + if e.is_circle_class() { + // copy the circle + return Ok(Some(e.copy(k)?)); + } + e.base() + } + CurveKind::Parabola(c) => c.base(), + CurveKind::Hyperbola(c) => c.base(), + CurveKind::Line(l) => { + line_base = l.base_with_nurb()?; + &line_base + } + _ => { + return Err(SolidError::invalid_input( + OP_IS_NURB_CURVE_CIRCLE, + "the curve is not an SmBSplineCurve", + )) + } + }; + + // no work - not a degree 2 rational curve + let deg = test_curve.degree()?; + if deg != 2 || !test_curve.is_rational()? { + return Ok(None); + } + + // scale zero value + let mid_point = + test_curve.evaluate_point(k, test_curve.natural_interval().evaluate(0.4567))?; + let scaled_zero = + tolerance_scale * ANALYTIC_TOL_SCALE * EFF_ZERO * (1.0 + mid_point.max_dimension()); + + // skip degenerate curves + if test_curve.is_degenerate(k, scaled_zero, None)? { + return Ok(None); + } + + // Get circle reference frame, radius and angle extent + let Some((ref_frame, radius, start_ang_deg, end_ang_deg)) = + test_curve.is_arc(k, 5, scaled_zero)? + else { + // when isoParameterCurve is not a circular Arc + return Ok(None); + }; + let anal_domain = Extent1d::from_min_max(start_ang_deg, end_ang_deg); + + // make new circle + let circle = Self::new_circle( + k, + ref_frame.origin(), + ref_frame.x_axis(), + ref_frame.y_axis(), + anal_domain, + radius, + &EllipseOptions { + dim: 3, + nurb: Some(test_base), + inside_out: false, + }, + )?; + // all done + Ok(Some(circle)) + } + /// The offset of the curve by `offset_distance` in the plane with normal /// `offset_plane_normal` within `approx_tol_3d`, and its maximum gap /// (`CreateSimpleOffset`). For a circle (`SmCircle::CreateSimpleOffset`) whose @@ -2213,6 +2312,39 @@ mod tests { assert!(c.polar_converter().is_polar_conversion_possible()); } + /// #1432: a full unit circle given as its rational quadratic NURBS (9 points, + /// weights `1, s, 1, s, ..`, a plain B-spline) is recognised as that circle. + #[test] + fn is_nurb_curve_circle_full_unit_circle() { + let k = kernel(); + let circle = + Ellipse::create_canonical_circle(&k, &Axis2Placement::default(), 1.0, None, None) + .unwrap(); + let nurbs = Curve::from(circle.base().copy()); + assert!(matches!(nurbs.kind(), CurveKind::BSpline(_))); + let got = Ellipse::is_nurb_curve_circle(&k, &nurbs, None) + .unwrap() + .expect("a circle"); + assert!(got.is_circle_class()); + assert!((got.radius() - 1.0).abs() < 1e-12); + let dom = got.step_interval(); + assert!((dom.max() - dom.min() - 360.0).abs() < 1e-9); + assert!(got.base().gw_nurb().is_some()); + + // already a circle: a copy + let same = Ellipse::is_nurb_curve_circle(&k, &Curve::from(circle.copy(&k).unwrap()), None) + .unwrap() + .expect("a copy"); + assert!(same.is_circle_class()); + assert_eq!(same.radius(), 1.0); + + // a non-circular ellipse is rational degree 2 but not an arc + let ell = Ellipse::create_canonical(&k, &Axis2Placement::default(), 2.0, 1.0, None, None) + .unwrap(); + let none = Ellipse::is_nurb_curve_circle(&k, &Curve::from(ell.base().copy()), None); + assert!(none.unwrap().is_none()); + } + #[test] fn simple_offset_and_join() { let k = kernel(); diff --git a/harness/cpp/shim/curve_ellipse_p1.cpp b/harness/cpp/shim/curve_ellipse_p1.cpp index b9752b5610..84aab747fa 100644 --- a/harness/cpp/shim/curve_ellipse_p1.cpp +++ b/harness/cpp/shim/curve_ellipse_p1.cpp @@ -22,6 +22,7 @@ #include #include #include +#include #include #include #include @@ -62,7 +63,8 @@ enum : int32_t OP_EDIT_PARAM = 18, OP_NOTIFY = 19, OP_IS_ARC = 20, - OP_LAST = 20, + OP_IS_NURB_CIRCLE = 21, + OP_LAST = 21, }; // Constructions. @@ -556,6 +558,56 @@ void run_op(int32_t op, Rd& rd, const SmContext& ctx, SmEllipse* a, std::vector< o.push_back(s1); break; } + case OP_IS_NURB_CIRCLE: + { + // the tested curve: A itself (virtual dispatch: a circle is copied), a plain + // SmBSplineCurve with A's NURB, or a raw 3D NURB `n p m cpts(x y z w) knots` + // (bounds and validity checked by the FFI) + int64_t src = rd.idx(0, 2); + const SmBSplineCurve* test = a; + SmBSplineCurve* owned = NULL; + if (src == 1) + owned = new (ctx) SmBSplineCurve(3, a->GetGwNurbPointer()); + else if (src == 2) + { + gw_INDEX n = (gw_INDEX)rd.idx(1, 64); + gw_DEGREE p = (gw_DEGREE)rd.idx(1, 8); + gw_INDEX m = (gw_INDEX)rd.idx(2, 73); + if (m != n + p + 1) + throw Bad(); + gw_CURVE* cur = sm_AllocateNurbCurve(n, p, m); + for (gw_INDEX i = 0; i <= n; i++) + { + cur->pol->Pw[i].x = rd.num(); + cur->pol->Pw[i].y = rd.num(); + cur->pol->Pw[i].z = rd.num(); + cur->pol->Pw[i].w = rd.num(); + } + for (gw_INDEX i = 0; i <= m; i++) + cur->knt->U[i] = rd.num(); + owned = new (ctx) SmBSplineCurve(cur, 3, FALSE); + } + if (owned != NULL) + test = owned; + bool has_scale = rd.flag(); + double scale = rd.num(); + SmCircle* circle = NULL; + SmBoolean ok = has_scale ? SmCircle::IsNurbCurveCircle(ctx, test, circle, scale) + : SmCircle::IsNurbCurveCircle(ctx, test, circle); + o.push_back(0.0); + o.push_back(flg(ok)); + if (ok) + { + if (circle == NULL) + throw Bad(); + o.push_back(flg(circle->IsKindOf(SmCircle_TYPE))); + o.push_back(circle->GetRadius()); + put_state(o, circle); + } + delete circle; + delete owned; + break; + } default: throw Bad(); } } diff --git a/plan/units/curve.ellipse_p1.toml b/plan/units/curve.ellipse_p1.toml index 604c4933b1..eaa69a7c0f 100644 --- a/plan/units/curve.ellipse_p1.toml +++ b/plan/units/curve.ellipse_p1.toml @@ -9,4 +9,4 @@ tests = ["prog_test:my_test_extent3d", "prog_test:my_test_vector2d", "prog_test: status = "green" claimed_by = "swarm-curve.ellipse_p1-a1" claimed_at = "2026-10-13T04:34:25Z" -notes = "main CI e33985411 (2026-10-23): PASS curve.ellipse_p1 cases=2031 passed=2031 failed=0 deferred=0. GREEN (promoted by curve.curve_p5: 2028 cases, 0 deferred after its stubs were wired); PARTIAL: difftest 2028 cases (28 fixture: my_test_analytic_curves circle/ellipse (CreateCanonical, AdjustSTEPInterval, GlobalPointSolveSTEP, ReverseParameterization), smSurfaceTest cylinder circle inside-out, smTopoTest CreateCanonical + AdjustSTEPInterval / EditParameterization circles, smCurveTest circle copy constructor, edge cases: non-positive radius, collinear / 2D / closed 3-point circles, joins (overlap, flipped normal, non-ellipse), circle offsets, skewed ConvertTo2D, 2D construction, Trim modes; + 2000 generated programs over constructors, CreateCanonical and ReadFromDB records), 1896 pass bit-identical (max_rel_err 0), 132 deferred (curve.curve_p5 132: LocalPointSolve / LocalPropertyAnalysis / SmCurve::IsArc); --perturb SENSITIVE; SmEllipse and SmCircle join Curve as CurveKind::Ellipse (Curve methods whose ellipse override rebuilds the polar converter take &Kernel); wired UNPORTED(curve.ellipse_p1) in curve.curve_p2 (ReadFromDB), curve.ellipse_p2 (GlobalPointSolve, GlobalCoincidenceChecker), nurbs.b_spline_curve_p2 (CreateAnalyticCurve SmCircle), nurbs.b_spline_surface_p2 (sweep circle). 39 fns; largest: SmEllipse::MakeNurb, SmPolarConversion::ConvertToPolarParameter, SmEllipse::JoinWith, SmEllipse::Trim; cycle back-edges (call into later units; stub or co-port): surface.surf_of_extrusion_p1(1), surface.surf_of_revolution_p1(1); unit-SCC #22 (220 units); 319 upstream test entry points reach this unit" +notes = "SmCircle::IsNurbCurveCircle ported (#1432) as Ellipse::is_nurb_curve_circle (shim op IsNurbCircle), +476 inline cases (76 fixed: NURBS circles/arcs, raw rational arcs, weight/point perturbations, tolerance scales, non-circles, degenerate, huge/tiny scales; 400 generated, fixed seed), all bit-identical. main CI e33985411 (2026-10-23): PASS curve.ellipse_p1 cases=2031 passed=2031 failed=0 deferred=0. GREEN (promoted by curve.curve_p5: 2028 cases, 0 deferred after its stubs were wired); PARTIAL: difftest 2028 cases (28 fixture: my_test_analytic_curves circle/ellipse (CreateCanonical, AdjustSTEPInterval, GlobalPointSolveSTEP, ReverseParameterization), smSurfaceTest cylinder circle inside-out, smTopoTest CreateCanonical + AdjustSTEPInterval / EditParameterization circles, smCurveTest circle copy constructor, edge cases: non-positive radius, collinear / 2D / closed 3-point circles, joins (overlap, flipped normal, non-ellipse), circle offsets, skewed ConvertTo2D, 2D construction, Trim modes; + 2000 generated programs over constructors, CreateCanonical and ReadFromDB records), 1896 pass bit-identical (max_rel_err 0), 132 deferred (curve.curve_p5 132: LocalPointSolve / LocalPropertyAnalysis / SmCurve::IsArc); --perturb SENSITIVE; SmEllipse and SmCircle join Curve as CurveKind::Ellipse (Curve methods whose ellipse override rebuilds the polar converter take &Kernel); wired UNPORTED(curve.ellipse_p1) in curve.curve_p2 (ReadFromDB), curve.ellipse_p2 (GlobalPointSolve, GlobalCoincidenceChecker), nurbs.b_spline_curve_p2 (CreateAnalyticCurve SmCircle), nurbs.b_spline_surface_p2 (sweep circle). 39 fns; largest: SmEllipse::MakeNurb, SmPolarConversion::ConvertToPolarParameter, SmEllipse::JoinWith, SmEllipse::Trim; cycle back-edges (call into later units; stub or co-port): surface.surf_of_extrusion_p1(1), surface.surf_of_revolution_p1(1); unit-SCC #22 (220 units); 319 upstream test entry points reach this unit"