Merge pull request #1456 from brett19/fix/face-props-split-at-knots

fix: precise face properties split at knot spans (fork 08d9a14) as opt-in flag (fixes #1435)

crates/difftest/src/solidcraft.rs

@@ -1428,6 +1428,9 @@ fn enable_fix(k: &mut Kernel, name: &str) -> Result<(), String> {
         "cone_parallel_line_midpoint" => &mut k.fixes.cone_parallel_line_midpoint,
         "analytic_copy_keeps_parameterization" => &mut k.fixes.analytic_copy_keeps_parameterization,
         "precise_properties_body_length_scale" => &mut k.fixes.precise_properties_body_length_scale,
+        "precise_properties_split_at_surface_knots" => {
+            &mut k.fixes.precise_properties_split_at_surface_knots
+        }
         "revolve_nested_shells_as_cavities" => &mut k.fixes.revolve_nested_shells_as_cavities,
         "sphere_small_circle_not_touching" => &mut k.fixes.sphere_small_circle_not_touching,
         _ => return Err(format!("unknown `fixes` entry `{name}`")),

crates/difftest/src/units/mass_integrator.rs

@@ -293,7 +293,8 @@ fn geometry_cpp(g: &Geometry) -> Result<Output, HarnessError> {
         .map(output)
         .map_err(|e| HarnessError::new(e.to_string()))
 }
-/// The shim's face surface: kind 3 is a quadratic B-spline patch, every other kind an
+/// The shim's face surface: kind 3 is a quadratic B-spline patch, kind 5 a multi-span
+/// rational quadratic patch (u knot 0.4, double v knot 0.5; #1435), every other kind an
 /// `SmPlane` over `domain` (kind 4 with rotated axes).
 fn face_surface(
     kind: i32,
@@ -330,6 +331,34 @@ fn face_surface(
                 &[0., 0., 0., 1., 1., 1.],
             )?,
         )
+    } else if kind == 5 {
+        let cp: Vec<Vec<_>> = (0..4)
+            .map(|i: i32| {
+                (0..5)
+                    .map(|j: i32| {
+                        CPoint::new(
+                            p[0] + f64::from(i) * p[3] / 3.0,
+                            p[1] + f64::from(j) * p[4] / 4.0,
+                            p[2] + f64::from(i * (3 - i) * j * (4 - j)) * 0.125,
+                            if j == 2 { 0.8 } else { 1.0 },
+                        )
+                    })
+                    .collect()
+            })
+            .collect();
+        solidcore::surface::surface::Surface::from(
+            solidcore::nurbs::b_spline_surface::BSplineSurface::from_data(
+                3,
+                4,
+                2,
+                2,
+                6,
+                7,
+                &cp,
+                &[0., 0., 0., 0.4, 1., 1., 1.],
+                &[0., 0., 0., 0.5, 0.5, 1., 1., 1.],
+            )?,
+        )
     } else {
         solidcore::surface::surface::Surface::from(
             solidcore::surface::plane::Plane::new(

crates/solidcore-ffi/src/units/mass_integrator.rs

@@ -148,7 +148,8 @@ pub fn edge(
     })
 }
 /// Face kinds: rectangular plane (0), disk (1), annulus (2), quadratic B-spline
-/// patch (3), rotated plane (4). Disk radii are 0.5 and 0.2.
+/// patch (3), rotated plane (4), multi-span rational quadratic B-spline patch (5: u knot
+/// 0.4, double v knot 0.5). Disk radii are 0.5 and 0.2.
 /// Precise face properties (0), actual oriented trim export (1), API area
 /// (2), centroid (3), B-rep area (4), volume (5), mass properties (6), rectangular surface area (7).
 /// Input: plane origin xyz, UV scales, domain u0 v0 u1 v1, computation origin xyz,
@@ -156,7 +157,7 @@ pub fn edge(
 pub fn face(op: i32, kind: i32, p: &[f64; 16]) -> Result<Vec<f64>, FfiError> {
     let [_, _, _, sx, sy, u0, v0, u1, v1, _, _, _, _, _, orientation, flags] = *p;
     if !(0..=7).contains(&op)
-        || !(0..=4).contains(&kind)
+        || !(0..=5).contains(&kind)
         || !p.iter().all(|v| v.is_finite())
         || ((kind == 1 || kind == 2)
             && (sx != 1.0 || sy != 1.0 || u0 > -0.5 || v0 > -0.5 || u1 < 0.5 || v1 < 0.5))

crates/solidcore/src/kernel.rs

@@ -283,6 +283,8 @@ pub struct UpstreamFixes {
     pub fillet_through_extended_rail_end: bool,
     /// `SmFace::ComputePreciseProperties` scales the volume and moment error bounds by the body size (ten times the root of the estimated area, or the origin's distance when larger), not 1.0.1's tenth of it, so the integration converges at tight accuracies instead of returning `SM_ERR` (WR-472; the SolidCraft fork's tangency-labs/smlib@073d4740).
     pub precise_properties_body_length_scale: bool,
+    /// `SmFace::ComputePreciseProperties` integrates each UV trim curve piece by piece between its knots and its crossings of the surface's u knot lines, and each iso-u integral of `SmODEAreaIFEO` piece by piece between the surface's v knots, so no Runge-Kutta step crosses a knot of a B-spline surface (a NURB cylinder's, sphere's or torus's rational-circle double knots) where the integrand's derivatives jump and the step's error estimate does not bound its error; the requested accuracy then holds at 1e-8 (WR-494, #1435; the SolidCraft fork's tangency-labs/smlib@08d9a14).
+    pub precise_properties_split_at_surface_knots: bool,
     /// An offset surface evaluated beyond its base's domain extends the base itself on every side that can grow at once, each by the domain size times the least power of two that reaches the point and never less than the extension it replaces, so two copies of one offset give the same points whatever was evaluated first (WR-458, #419; the fork's tangency-labs/smlib@2d7c7ed).
     pub offset_extension_from_base: bool,
     /// `N_FillNSidedHole` (through `SmBSplineSurface::CreateCornerBlend`) shortens the cross-boundary derivative across an inner curve where its cubic blend from the opposite boundary would leave that boundary's tangent plane, so a corner patch (a fillet setback corner) stays inside planar faces it meets, with G1 kept (WR-461, #1170). Curved neighbours are not covered.

crates/solidcore/src/mass/integrator/properties.rs

@@ -18,7 +18,9 @@
 //! B-rep enumeration, bounding boxes, and region classification remain calls into
 //! their unported units. `boundary_properties` exposes the integration after trim
 //! construction for callers already holding oriented trim curves. Nested integration
-//! has a fixed two-level call structure, not input-driven recursion.
+//! has a fixed two-level call structure, not input-driven recursion. With
+//! `UpstreamFixes::precise_properties_split_at_surface_knots` (WR-494) both levels split
+//! at the surface's knots, as the SolidCraft fork's tangency-labs/smlib@08d9a14 does.
 use super::{OdeIntegrator, OdeOptions};
 use crate::curve::curve::Curve;
 use crate::geom::attribute::{ids, AttributeStore, AttributeValue, ObjectKey};
@@ -31,7 +33,7 @@ use crate::nurbs::b_spline_surface_n_lib::p2::OrientType;
 use crate::query::brep::{self as qbrep, BoundingBoxOptions};
 use crate::surface::surface::Surface;
 use crate::topology::face::p2::{
-    OdeArea, PrecisePropertyLayout, PPF_ALL, PPF_AREA, PPF_AREA_FIRST_MOMENTS,
+    trim_curve_breaks, OdeArea, PrecisePropertyLayout, PPF_ALL, PPF_AREA, PPF_AREA_FIRST_MOMENTS,
 };
 use crate::topology::{edge, vertex, Brep, EdgeId, FaceId, FaceuseId, RegionId, VertexId};
 use crate::{EntityRef, Result, SolidError};
@@ -333,7 +335,16 @@ pub fn boundary_properties(
         // Integrate each knot span of the natural interval separately: the RK error
         // estimate assumes a smooth integrand, so a step across a knot can hide a large
         // error.
-        let (knots, _) = trim.knots(k, Some(&ivl))?;
+        //
+        // With `fixes.precise_properties_split_at_surface_knots` (WR-494; the SolidCraft
+        // fork's tangency-labs/smlib@08d9a14) the pieces are those between the trim
+        // curve's knots and its crossings of the surface's knot lines u = k, where the
+        // integrand's derivatives jump too.
+        let knots = if k.fixes.precise_properties_split_at_surface_knots {
+            trim_curve_breaks(k, trim, surface)?
+        } else {
+            trim.knots(k, Some(&ivl))?.0
+        };
 
         // A coarse first Runge-Kutta step can trip the integrator's abort guard for one
         // traversal direction, so a reversed cap loop fails where its twin succeeds. Try
@@ -345,6 +356,12 @@ pub fn boundary_properties(
             let &[a, b] = pair else {
                 continue;
             };
+            // (`!(dB > dA)`: also an unordered pair)
+            if k.fixes.precise_properties_split_at_surface_knots
+                && b.partial_cmp(&a) != Some(std::cmp::Ordering::Greater)
+            {
+                continue; // a crossing at a knot
+            }
             let mut status = Err(SolidError::sm_err(OP));
             for divisor in [2.0, 8.0, 64.0, 512.0, 4096.0] {
                 status = OdeIntegrator::default()

crates/solidcore/src/topology/face/p2.rs

@@ -20,6 +20,8 @@
 //! | [`PrecisePropertyLayout`]               | `SmPrecisePropertyLayout` (constructor, `Add`)    |
 //! | [`OdeIsoCurveLength`]                   | `SmODEIsoCurveLengthIFEO` (constructor, `Evaluate`) |
 //! | [`OdeArea`]                             | `SmODEAreaIFEO::Evaluate`                         |
+//! | [`surface_knots_inside`]                | `sm_SurfaceKnotsInside` (static; the SolidCraft fork's tangency-labs/smlib@08d9a14) |
+//! | [`trim_curve_breaks`]                   | `sm_TrimCurveBreaks` (static; the SolidCraft fork's tangency-labs/smlib@08d9a14) |
 //! | [`create_curves_by_trimming`], [`TrimmedCurves`] | `SmFace::CreateCurvesByTrimming`         |
 //! | [`create_curves_from_face`], [`FaceCurves`] | `SmFace::CreateCurvesFromFace`                |
 //! | [`fix_uv_trim_curves`], [`fix_uv_trim_curves_notified`] | `sm_FixUVTrimCurves`              |
@@ -2128,18 +2130,49 @@ impl OdeArea<'_> {
             let num_components = self.y_scale.len();
             let mut dy = Vec::new();
             dy.resize(num_components, 0.0);
-            let step = p0.y - self.min_v;
+            if k.fixes.precise_properties_split_at_surface_knots {
+                // Along the iso-u line the integrand is smooth only within each knot
+                // span of the surface: integrated span by span, as the trim curve is
+                // (WR-494; the SolidCraft fork's tangency-labs/smlib@08d9a14).
+                let mut breaks =
+                    surface_knots_inside(k, self.surface, SurfParamType::V, self.min_v, p0.y)?;
+                breaks.insert(0, self.min_v);
+                breaks.push(p0.y);
+                let mut piece_y = Vec::new();
+                piece_y.resize(num_components, 0.0);
+                for pair in breaks.windows(2) {
+                    let &[a, b] = pair else {
+                        continue;
+                    };
+                    piece_y.fill(0.0);
+                    // A failed inner integration is not a zero-area contribution.
+                    runge_kutta_ode_integrate(
+                        &eval,
+                        &mut piece_y,
+                        a,
+                        b,
+                        self.iso_curve_accuracy,
+                        (b - a) / 2.0,
+                        self.y_scale,
+                    )?;
+                    for (total, part) in dy.iter_mut().zip(piece_y.iter()) {
+                        *total += *part;
+                    }
+                }
+            } else {
+                let step = p0.y - self.min_v;
 
-            // A failed inner integration is not a zero-area contribution.
-            runge_kutta_ode_integrate(
-                &eval,
-                &mut dy,
-                self.min_v,
-                p0.y,
-                self.iso_curve_accuracy,
-                step / 2.0,
-                self.y_scale,
-            )?;
+                // A failed inner integration is not a zero-area contribution.
+                runge_kutta_ode_integrate(
+                    &eval,
+                    &mut dy,
+                    self.min_v,
+                    p0.y,
+                    self.iso_curve_accuracy,
+                    step / 2.0,
+                    self.y_scale,
+                )?;
+            }
             dydx.clear();
             dydx.push(0.0);
             dydx.extend_from_slice(&dy);
@@ -2153,6 +2186,131 @@ impl OdeArea<'_> {
     }
 }
 
+/// `sm_SurfaceKnotsInside` (the SolidCraft fork's tangency-labs/smlib@08d9a14): the
+/// unique knots of a B-spline-kind surface in `param`'s direction strictly inside
+/// `(lo, hi)`, ascending; none for a surface not of the B-spline kind (smooth
+/// everywhere). The knots are read as upstream's virtual `GetKnots` does, with the NURB
+/// of a lazy analytic surface made ([`Surface::knots_with_nurb`]).
+///
+/// # Errors
+///
+/// Those of [`Surface::is_kind_of`] and [`Surface::knots_with_nurb`].
+pub fn surface_knots_inside(
+    k: &Kernel,
+    surface: &Surface,
+    param: SurfParamType,
+    lo: f64,
+    hi: f64,
+) -> Result<Vec<f64>> {
+    let mut knots = Vec::new();
+    if !surface.is_kind_of(crate::surface::surface::p2::B_SPLINE_SURFACE_TYPE)? {
+        return Ok(knots);
+    }
+    let (all, _) = surface.knots_with_nurb(k, param, None)?;
+    for &knot in &all {
+        if knot > lo && knot < hi {
+            knots.push(knot);
+        }
+    }
+    Ok(knots)
+}
+
+/// `sm_TrimCurveBreaks` (the SolidCraft fork's tangency-labs/smlib@08d9a14): where the
+/// precise-property integrand is not smooth along the UV trim curve `uv_curve` on
+/// `surface`: the curve's own knots over its natural interval, and the parameters where
+/// it crosses a knot line `u = k` of a B-spline surface (found per knot span by a sign
+/// change of `u(t) - k` between `4 (degree + 1)` evenly spaced points, then at most 64
+/// bisections). Ascending, the curve's ends included; breaks the integrator cannot step
+/// between are merged into one, a curve knot kept over a crossing.
+///
+/// # Errors
+///
+/// `SM_ERR` for fewer than two curve knots; those of the curve's knots, bounding box,
+/// degree and evaluation and of [`surface_knots_inside`].
+pub fn trim_curve_breaks(k: &Kernel, uv_curve: &Curve, surface: &Surface) -> Result<Vec<f64>> {
+    let mut breaks = Vec::new();
+    let ivl = uv_curve.natural_interval();
+    let (knots, _) = uv_curve.knots(k, Some(&ivl))?;
+    if knots.len() < 2 {
+        return Err(SolidError::sm_err("topology::face::trim_curve_breaks"));
+    }
+    let req = crate::nurbs::b_spline_curve::p2::BoundingBoxRequest {
+        normal: true,
+        pseudo: false,
+        polar: false,
+    };
+    let bbox = uv_curve
+        .calculate_bounding_box(k, &uv_curve.natural_interval(), req, true)?
+        .normal
+        .ok_or_else(|| null("trim_curve_breaks: bounding box not computed"))?;
+    let u_lines = surface_knots_inside(k, surface, SurfParamType::U, bbox.min().x, bbox.max().x)?;
+
+    let samples = 4 * (uv_curve.degree()? + 1);
+    let u = |t: f64| -> Result<f64> { Ok(uv_curve.evaluate_point(k, t)?.x) };
+    for pair in knots.windows(2) {
+        let &[a, b] = pair else {
+            continue;
+        };
+        breaks.push(a);
+        for &line in &u_lines {
+            let mut t0 = a;
+            let mut u0 = u(t0)?;
+            for i in 1..=samples {
+                // (`(dB - dA) * iS / lSamples`: an unsigned count converted to double)
+                let t1 = a + (b - a) * (i as f64) / (samples as f64);
+                let u1 = u(t1)?;
+                if (u0 - line) * (u1 - line) < 0.0 {
+                    let mut lo = t0;
+                    let mut hi = t1;
+                    let mut u_lo = u0;
+                    let mut bisections = 0;
+                    while bisections < 64 && hi - lo > EFF_ZERO * (1.0 + hi.abs()) {
+                        let mid = 0.5 * (lo + hi);
+                        let u_mid = u(mid)?;
+                        if (u_lo - line) * (u_mid - line) <= 0.0 {
+                            hi = mid;
+                        } else {
+                            lo = mid;
+                            u_lo = u_mid;
+                        }
+                        bisections += 1;
+                    }
+                    breaks.push(0.5 * (lo + hi));
+                } else if u1 == line && i < samples {
+                    breaks.push(t1);
+                }
+                t0 = t1;
+                u0 = u1;
+            }
+        }
+    }
+    if let Some(&last) = knots.last() {
+        breaks.push(last);
+    }
+    breaks.sort_by(|x, y| x.partial_cmp(y).unwrap_or(std::cmp::Ordering::Equal));
+
+    // A crossing found next to a knot (a corner of the trim curve on a knot line) is the
+    // knot: breaks the integrator cannot step between are one, the curve's own knots and
+    // ends kept. The integrator takes no step when its ends are `SM_ARE_SAME` or its
+    // first step is `SM_ARE_SAME` as `SM_EFF_ZERO`, and the first step is half the span.
+    let unsteppable = |a: f64, b: f64| {
+        are_same(a, b) || are_same(0.5 * (b - a), EFF_ZERO) || 0.5 * (b - a) < EFF_ZERO
+    };
+    let mut distinct: Vec<f64> = Vec::new();
+    for &brk in &breaks {
+        if let Some(prev) = distinct.last_mut() {
+            if unsteppable(*prev, brk) {
+                if knots.contains(&brk) {
+                    *prev = brk;
+                }
+                continue;
+            }
+        }
+        distinct.push(brk);
+    }
+    Ok(distinct)
+}
+
 // ---- face curves ----------------------------------------------------------------------------
 
 /// `new SmBSplineCurve(*SM_CAST_PTR(SmBSplineCurve, pCurve))`: the copy constructor of

crates/solidcore/tests/issue_1435_precise_properties_split_at_knots.rs

@@ -0,0 +1,147 @@
+//! #1435 (WR-494): with `UpstreamFixes::precise_properties_split_at_surface_knots`
+//! (the SolidCraft fork's tangency-labs/smlib@08d9a14) `SmFace::ComputePreciseProperties`
+//! splits both integration levels at the surface's knots, so curved primitives' areas and
+//! inertia hold 1e-9 of their closed forms at a requested accuracy of 1e-8.
+#![allow(clippy::unwrap_used, clippy::indexing_slicing)]
+use solidcore::geom::attribute::{AttributeStore, ObjectKey};
+use solidcore::mass::integrator::properties::{
+    brep_area, brep_mass_properties, face_area, PropertyAttributes,
+};
+use solidcore::math::Vec3;
+use solidcore::primitive::api_primitives::{
+    create_cone, create_cylinder, create_sphere, create_torus,
+};
+use solidcore::topology::brep::p1::faces;
+use solidcore::Kernel;
+use std::f64::consts::PI;
+
+const ACCURACY: f64 = 1e-8;
+const REL: f64 = 1e-9;
+
+fn fixed() -> Kernel {
+    let mut k = Kernel::default();
+    k.fixes.precise_properties_split_at_surface_knots = true;
+    k
+}
+
+fn assert_rel(name: &str, got: f64, exact: f64) {
+    let err = (got - exact).abs() / exact.abs();
+    assert!(err <= REL, "{name}: {got} vs {exact} ({err:e} relative)");
+}
+
+#[test]
+fn flag_is_off_by_default() {
+    assert!(
+        !Kernel::default()
+            .fixes
+            .precise_properties_split_at_surface_knots
+    );
+}
+
+/// The issue's repro: the cylindrical face of `create_cylinder`.
+#[test]
+fn cylinder_side_face_area_holds_its_accuracy() {
+    let mut k = fixed();
+    let mut b = create_cylinder(&mut k, Vec3::new(1.0, 2.0, -3.0), 5.0, 12.0).unwrap();
+    let side = faces(&b, None).unwrap()[0];
+    let a = face_area(&k, &mut b, side, ACCURACY).unwrap();
+    assert_rel("cylinder side", a, 2.0 * PI * 5.0 * 12.0);
+}
+
+/// The issue's table: the whole brep's area of each curved primitive.
+#[test]
+fn curved_primitive_areas_hold_their_accuracy() {
+    let mut k = fixed();
+    let origin = Vec3::new(1.0, 2.0, -3.0);
+    let s_cone = (36.0_f64 + 81.0).sqrt();
+    let s_frustum = (16.0_f64 + 81.0).sqrt();
+    let cases: [(&str, solidcore::topology::Brep, f64); 5] = [
+        (
+            "cylinder r5 h12",
+            create_cylinder(&mut k, origin, 5.0, 12.0).unwrap(),
+            2.0 * PI * 25.0 + 2.0 * PI * 5.0 * 12.0,
+        ),
+        (
+            "sphere r7",
+            create_sphere(&mut k, origin, 7.0).unwrap(),
+            4.0 * PI * 49.0,
+        ),
+        (
+            "cone r6 h9",
+            create_cone(&mut k, origin, 6.0, 0.0, 9.0).unwrap(),
+            PI * (36.0 + 6.0 * s_cone),
+        ),
+        (
+            "frustum r6..2 h9",
+            create_cone(&mut k, origin, 6.0, 2.0, 9.0).unwrap(),
+            PI * (36.0 + 4.0 + 8.0 * s_frustum),
+        ),
+        (
+            "torus R10 r3",
+            create_torus(&mut k, origin, 10.0, 3.0).unwrap(),
+            4.0 * PI * PI * 30.0,
+        ),
+    ];
+    for (name, mut b, exact) in cases {
+        let a = brep_area(&k, &mut b, ACCURACY).unwrap();
+        assert_rel(name, a, exact);
+    }
+}
+
+/// The issue's `brep_properties` case: a torus R10 r3 at (1, -1, 0.5); its axial
+/// inertia about the centre against V (R^2/2 + 5 r^2/8).
+#[test]
+fn torus_mass_properties_hold_their_accuracy() {
+    let mut k = fixed();
+    let centre = Vec3::new(1.0, -1.0, 0.5);
+    let (rr, r) = (10.0, 3.0);
+    let mut b = create_torus(&mut k, centre, rr, r).unwrap();
+    let store = AttributeStore::default();
+    let key_of = |_| ObjectKey(0);
+    let attributes = PropertyAttributes {
+        store: &store,
+        key_of: &key_of,
+    };
+    let m = brep_mass_properties(&k, &mut b, ACCURACY, 1.0, centre, &attributes).unwrap();
+    let volume = 2.0 * PI * PI * rr * r * r;
+    assert_rel("torus area", m.properties.area, 4.0 * PI * PI * rr * r);
+    assert_rel("torus volume", m.properties.volume, volume);
+    let ixx = volume * (rr * rr / 2.0 + 5.0 * r * r / 8.0);
+    let izz = volume * (rr * rr + 3.0 * r * r / 4.0);
+    assert_rel("torus Ixx", m.inertia.x, ixx);
+    assert_rel("torus Iyy", m.inertia.y, ixx);
+    assert_rel("torus Izz", m.inertia.z, izz);
+}
+
+/// Flag off (upstream v1.0.1, the default): the issue's repro misses the accuracy it
+/// asks for, as the C++ reference does: the cylinder's side face 5.2e-8 off 2πrh.
+#[test]
+fn upstream_cylinder_side_face_area_misses_its_accuracy() {
+    let mut k = Kernel::default();
+    let mut b = create_cylinder(&mut k, Vec3::new(1.0, 2.0, -3.0), 5.0, 12.0).unwrap();
+    let side = faces(&b, None).unwrap()[0];
+    let a = face_area(&k, &mut b, side, ACCURACY).unwrap();
+    let exact = 2.0 * PI * 5.0 * 12.0;
+    let err = (a - exact).abs() / exact;
+    assert!(err > ACCURACY, "{a} vs {exact} ({err:e} relative)");
+}
+
+/// Flag off: the issue's torus inertia, 1.7e-9 off V (R^2/2 + 5 r^2/8).
+#[test]
+fn upstream_torus_inertia_misses_1e_9() {
+    let mut k = Kernel::default();
+    let centre = Vec3::new(1.0, -1.0, 0.5);
+    let (rr, r) = (10.0, 3.0);
+    let mut b = create_torus(&mut k, centre, rr, r).unwrap();
+    let store = AttributeStore::default();
+    let key_of = |_| ObjectKey(0);
+    let attributes = PropertyAttributes {
+        store: &store,
+        key_of: &key_of,
+    };
+    let m = brep_mass_properties(&k, &mut b, ACCURACY, 1.0, centre, &attributes).unwrap();
+    let volume = 2.0 * PI * PI * rr * r * r;
+    let ixx = volume * (rr * rr / 2.0 + 5.0 * r * r / 8.0);
+    let err = (m.inertia.x - ixx).abs() / ixx;
+    assert!(err > REL, "{} vs {ixx} ({err:e} relative)", m.inertia.x);
+}

harness/cpp/shim/mass_integrator.cpp

@@ -141,11 +141,12 @@ SC_EXPORT int32_t sc_mass_integrator_edge(int32_t op,int32_t kind,const double*
 
 // p: plane origin[3], UV scales[2], UV domain[4], computation origin[3], accuracy,
 // thickness, orientation (0 SAME, 1 OPPOSITE), property flags.
-// kind 0 rectangle, 1 disk, 2 annulus, 3 quadratic B-spline patch, 4 rotated plane.
+// kind 0 rectangle, 1 disk, 2 annulus, 3 quadratic B-spline patch, 4 rotated plane,
+// 5 multi-span rational quadratic B-spline patch (u knot 0.4, double v knot 0.5).
 // op 0 properties, 1 export input UV domain/bbox scale/oriented trims (NURB, then 1 + SmLine members or 0), 2 API area, 3 API centroid,
 // 4 API Brep area, 5 API Brep volume, 6 API Brep mass properties (sheet rejection), 7 Surface::Area.
 SC_EXPORT int32_t sc_mass_integrator_face(int32_t op,int32_t kind,const double* p,double* out,int64_t cap,int64_t* len){
-    if(op<0||op>7||kind<0||kind>4)return SC_ERR_ARG;
+    if(op<0||op>7||kind<0||kind>5)return SC_ERR_ARG;
     SC_TRY
     SmContext ctx;
     std::unique_ptr<SmBrep> brep(new(ctx)SmBrep());
@@ -163,6 +164,20 @@ SC_EXPORT int32_t sc_mass_integrator_face(int32_t op,int32_t kind,const double*
         }
         for(int i=0;i<6;++i){s->knu->U[i]=i<3?0:1;s->knv->U[i]=i<3?0:1;}
         surface=new(ctx)SmBSplineSurface(s,FALSE,NULL);
+    }else if(kind==5){
+        // 4 x 5 control points: knot spans [0,.4],[.4,1] in u and [0,.5],[.5,1] in v,
+        // the v knot double (C0, as a rational circle's), weights .8 on the middle row.
+        gw_SURFACE* s=sm_AllocateNurbSurface(3,4,2,2,6,7);
+        for(int i=0;i<4;++i)for(int j=0;j<5;++j){
+            s->net->Pw[i][j].x=p[0]+i*p[3]/3.0;
+            s->net->Pw[i][j].y=p[1]+j*p[4]/4.0;
+            s->net->Pw[i][j].z=p[2]+(double)(i*(3-i)*j*(4-j))*0.125;
+            s->net->Pw[i][j].w=j==2?0.8:1.0;
+        }
+        const double ku[7]={0,0,0,.4,1,1,1},kv[8]={0,0,0,.5,.5,1,1,1};
+        for(int i=0;i<7;++i)s->knu->U[i]=ku[i];
+        for(int i=0;i<8;++i)s->knv->U[i]=kv[i];
+        surface=new(ctx)SmBSplineSurface(s,FALSE,NULL);
     }else surface=new(ctx)SmPlane(SmPoint3d(p[0],p[1],p[2]),xaxis,yaxis,SmVector2d(p[3],p[4]),dom,&ctx);
     SmFace* face=nullptr;
     SmStatus st;

plan/upstream-issues/WR-494--precise-properties-step-across-surface-knots.md

@@ -0,0 +1,9 @@
+#### WR-494. `SmFace::ComputePreciseProperties` steps across the surface's knots, so curved faces miss the requested accuracy (1.0.1)
+- **Severity:** Wrong result
+- **Area:** Mass properties
+- **Location:** `source/SMLib/src/SmFace.cpp:6233` at `3fd5b95` (`SER(pBSC->GetKnots(sKnots, NULL, &sIvl));`, the outer integration's pieces), `SmFace::ComputePreciseProperties` (def. `:6144`); `:6078` (`SER(sIntegrator.RungeKuttaODEIntegrate(dY,m_dMinV, sPV[0].y, ...))`, the inner iso-u integral), `SmODEAreaIFEO::Evaluate` (def. `:6042`).
+- **What happens:** the face's properties are a double integral: along each UV trim curve (the outer adaptive Runge-Kutta run) of the integral along the iso-u line from the face's lowest v (the inner run, `SmODEAreaIFEO`). 1.0.1 splits the outer run at the trim curve's own knots only. The integrand's derivatives also jump where the trim curve crosses a knot line `u = k` of a B-spline surface, and the inner integrand's where the iso-u line crosses a v knot. A NURB cylinder, cone, sphere or torus has a rational circle's double knots in both directions. A Runge-Kutta step across one is outside the method's error model: its error estimate does not bound its error, so the accuracy asked for does not hold.
+- **Trigger:** `SmApiFaceComputeArea` of the cylindrical face of `SmApiCreateCylinder` (origin (1, 2, -3), radius 5, height 12) at accuracy 1e-8: 376.99109865933315 against 2πrh = 376.99111843077515 (5.2e-8 relative). `SmApiBrepComputeArea` at 1e-8: cylinder r5 h12 3.7e-8, sphere r7 1.2e-8, cone r6 h9 4.2e-8, frustum r6..2 h9 3.5e-8, torus R10 r3 6.0e-9 off their closed forms. `SmApiBrepComputeMassProperties` of a torus R10 r3 at (1, -1, 0.5): Ixx about the centre 98819.4138975709 against V (R²/2 + 5r²/8) = 98819.4140659072 (1.7e-9).
+- **How found:** wiring solidcore into SolidCraft as a second kernel (tangency-labs/solidcore#1435): the SolidCraft fork of SMLib gives 1.5e-10 (cylinder area) and 1.8e-11 (torus inertia) for the same call sequence.
+- **Confirmed:** the Rust port's default (difftest-gated parity with the C++ reference) reproduces the numbers above exactly (`crates/solidcore/tests/issue_1435_precise_properties_split_at_knots.rs` with the flag off); the SolidCraft fork measured the same defect in its C++ build ("at 1e-8 a cylinder's area 5.6e-8", tangency-labs/smlib@08d9a14) and fixed it there.
+- **Fix (opt-in, not the default):** `Kernel::fixes.precise_properties_split_at_surface_knots` ports the fork's tangency-labs/smlib@08d9a14. The outer run integrates each trim curve piece by piece between its knots and its crossings of the surface's u knot lines (`sm_TrimCurveBreaks` -> `topology::face::p2::trim_curve_breaks`: a sign change of `u(t) - k` between `4 (degree + 1)` points of each knot span, then bisection; breaks the integrator cannot step between merged), skipping an empty piece. Each inner iso-u integral is split at the surface's v knots inside `(min v, v(t))` (`sm_SurfaceKnotsInside` -> `topology::face::p2::surface_knots_inside`), each piece from a half-length first step, the pieces summed. A surface not of the B-spline kind has no knots. With it on, the cases above hold 1e-9 of their closed forms at accuracy 1e-8 (the test above with the flag on).