implicit, torus: evaluates_as's bound derived term by term; the recognised torus's sweep frame from its form (#3185, #3199) (#13743)

Moved from tangency-labs/solidcore#3224.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

solidcore/crates/solidcore-convert/tests/flat_parity.rs

@@ -7,8 +7,8 @@
 //! kind, by what it became; a number within [`ROUNDING`] of its size (relative, no floor) is
 //! counted as rounding, with its count and the largest.
 //!
-//! On the whole corpus (1,332 bodies, at #3187): 1,256 come back bit for bit, 36 within
-//! rounding (the largest 4.9e-16 relative), and 40 differ, each in exactly the classes
+//! On the whole corpus (1,333 bodies, at #3199 with kernel case 11991 (#3231) in): 1,259 come
+//! back bit for bit, 36 within rounding (the largest 4.9e-16 relative), and 38 differ, each in exactly the classes
 //! [`pinned`] lists for it.
 #![allow(
     clippy::unwrap_used,
@@ -390,8 +390,8 @@ fn pinned() -> BTreeMap<String, BTreeSet<String>> {
             ],
         ),
         // a label its form is not proved on goes as the form alone (#2540): cylinders whose
-        // frames are about 1e-14 rad off their forms (two of them at the rounding of a derived
-        // bound, #3185's thread)
+        // frames are off their forms past the rounding `evaluates_as` derives (margins 1.23 to
+        // 2.03, tests/issue_3185_label_margins.rs)
         (
             label("cylinder -> rational_spline", &["radius"]),
             &[
@@ -401,18 +401,14 @@ fn pinned() -> BTreeMap<String, BTreeSet<String>> {
                 "fillets/2903-original-rim-round.json:b0",
             ],
         ),
-        // tori whose frame and fitted sweep are about 1e-14 rad off their forms, rightly
-        // refused (#3185)
+        // a torus label whose frame is off its form's own first sweep row past the rounding
+        // `evaluates_as` derives (margin 1.69, tests/issue_3185_label_margins.rs)
         (
             label(
                 "torus -> rational_spline",
                 &["major_radius", "minor_radius"],
             ),
-            &[
-                "sweeps/3897-sweep-ring-along-a-short-arc-joined-to-a-collar.json:b6",
-                "sweeps/4280-circular-pattern-of-sweep-cut-groove-on-hub.json:b46",
-                "sweeps/4280-circular-pattern-of-sweep-cut-groove-on-hub.json:b85",
-            ],
+            &["sweeps/4280-circular-pattern-of-sweep-cut-groove-on-hub.json:b46"],
         ),
         // an extruded circle the flat hop left a spline, proved a cylinder
         (
@@ -582,10 +578,12 @@ fn parity(stride: usize) {
             missing.is_empty(),
             "pinned bodies that no longer differ: {missing:#?}"
         );
+        // The corpus grows as cases are added, so the totals aren't pinned: every body that
+        // differs is pinned above, and the rest must be identical or within rounding.
         assert_eq!(
-            (t.bodies, t.identical, t.rounding_only),
-            (1332, 1256, 36),
-            "bodies, identical, within rounding"
+            (t.bodies - t.identical - t.rounding_only, t.rounding_only),
+            (pinned.len(), 36),
+            "bodies that differ, bodies within rounding"
         );
         // the largest rounding difference, relative (#3137: 4.885e-16)
         assert!(t.rounding_largest <= 4.9e-16, "{:e}", t.rounding_largest);

solidcore/crates/solidcore-convert/tests/issue_3185_label_margins.rs

@@ -0,0 +1,258 @@
+//! (#3185) An analytic label is kept on its NURBS form when the kernel's analytic surface,
+//! carrying the label's numbers, evaluates as the form within the rounding
+//! `AnalyticSupport::build` derives term by term (`surface::implicit`'s `evaluates_as`):
+//! [`AnalyticSupport::parametrisation_margin`] is the largest ratio of distance to bound over
+//! the form's Greville points, at most 1 where the label stands.
+//!
+//! - The corpus labels #3137's parity run refused under the old constant `gamma(16 (p + q +
+//!   2)) size`, each pinned to its margin as measured under the derived bound: the three tori
+//!   of b6 and b85 (0.02 to 0.50) stand; b46's torus (1.69) and the four cylinders (b94 1.23,
+//!   b43 1.40, the 2903 rims 1.52 and 2.03) do not: their frames are off their forms past
+//!   the rounding at their size. Since #3199 the kernel's recognised sweep frame is the form's
+//!   own (its end meridians), so a label is measured against the form, not against a fit
+//!   that shared the label's error.
+//! - A torus the kernel makes at the groove's sizes stands on its own form, and the same label
+//!   turned about its axis by four times what the bound allows is refused.
+#![allow(
+    clippy::unwrap_used,
+    clippy::expect_used,
+    clippy::panic,
+    clippy::indexing_slicing,
+    clippy::print_stderr,
+    clippy::disallowed_methods
+)] // test (the turn's cos and sin)
+
+use std::fmt::Write as _;
+use std::io::Read as _;
+use std::path::{Path, PathBuf};
+
+use serde_json::Value;
+use solidcore::blend::fillet_section::station::as_b_spline;
+use solidcore::geom::axis2_placement::Axis2Placement;
+use solidcore::geom::extent2d::Extent2d;
+use solidcore::geom::vector3d::{Point3, Vec3};
+use solidcore::kernel::{Kernel, KernelOptions, UpstreamFixes};
+use solidcore::surface::from_parameters::from_parameters;
+use solidcore::surface::implicit::AnalyticSupport;
+use solidcore::surface::surface::{Surface, SurfaceKind};
+use solidcore::surface::torus::{Torus, TorusOptions};
+use solidcore_record::brep::Brep as Record;
+
+fn corpus() -> PathBuf {
+    Path::new(env!("CARGO_MANIFEST_DIR")).join("../solidcore/tests/kernel_cases")
+}
+
+/// The corpus's kernel, its fixes on (as the parity run and SolidCraft run it).
+fn kernel() -> Kernel {
+    let v: Value =
+        serde_json::from_str(&std::fs::read_to_string(corpus().join("fixes.json")).unwrap())
+            .unwrap();
+    let mut fixes = UpstreamFixes::default();
+    for flag in v.as_array().unwrap() {
+        *fixes.flag_mut(flag.as_str().unwrap()).unwrap() = true;
+    }
+    let mut k = Kernel::new(&KernelOptions::default());
+    k.fixes = fixes;
+    k
+}
+
+fn body(case: &str, id: &str) -> Record {
+    let gz = corpus().join(format!("{case}.bodies.json.gz"));
+    let mut text = String::new();
+    flate2::read::GzDecoder::new(std::fs::File::open(&gz).unwrap())
+        .read_to_string(&mut text)
+        .unwrap();
+    let records: Value = serde_json::from_str(&text).unwrap();
+    serde_json::from_value(records[id].clone()).unwrap()
+}
+
+fn placement(f: &solidcore::eval::Frame) -> Axis2Placement {
+    let v = |a: [f64; 3]| Vec3::new(a[0], a[1], a[2]);
+    Axis2Placement::from_stored(
+        Point3::new(f.origin[0], f.origin[1], f.origin[2]),
+        v(f.x_axis),
+        v(f.y_axis),
+    )
+}
+
+/// A corpus face's label as the kernel's support, and its NURBS form.
+fn label_and_form(
+    k: &Kernel,
+    case: &str,
+    id: &str,
+    face: usize,
+) -> (
+    AnalyticSupport,
+    solidcore::nurbs::b_spline_surface::BSplineSurface,
+) {
+    let r = body(case, id);
+    let f = &r.faces[face];
+    let nurbs = solidcore_convert::surface(f.nurbs.as_ref().unwrap()).unwrap();
+    let (s, _) = from_parameters(k, &nurbs, f.domain).unwrap();
+    let SurfaceKind::BSpline(form) = s.kind() else {
+        panic!("{case} {id} {face}: the form is not a B-spline")
+    };
+    let support = match solidcore_convert::surface(&f.surface).unwrap() {
+        solidcore::eval::Surface::Torus {
+            frame,
+            major_radius,
+            minor_radius,
+            ..
+        } => AnalyticSupport::Torus {
+            placement: placement(&frame),
+            major_radius,
+            minor_radius,
+        },
+        solidcore::eval::Surface::Cylinder { frame, radius } => AnalyticSupport::Cylinder {
+            placement: placement(&frame),
+            radius,
+        },
+        other => panic!("{case} {id} {face}: {other:?}"),
+    };
+    (support, form.copy())
+}
+
+/// (case, body, face, the margin measured at #3185, whether the label stands)
+const CASES: &[(&str, &str, usize, f64, bool)] = &[
+    (
+        "sweeps/3897-sweep-ring-along-a-short-arc-joined-to-a-collar",
+        "b6",
+        1,
+        0.0201,
+        true,
+    ),
+    (
+        "sweeps/4280-circular-pattern-of-sweep-cut-groove-on-hub",
+        "b46",
+        1,
+        1.6943,
+        false,
+    ),
+    (
+        "sweeps/4280-circular-pattern-of-sweep-cut-groove-on-hub",
+        "b85",
+        15,
+        0.4995,
+        true,
+    ),
+    (
+        "sweeps/4280-circular-pattern-of-sweep-cut-groove-on-hub",
+        "b85",
+        21,
+        0.4245,
+        true,
+    ),
+    (
+        "booleans/3456-deep-cut-through-drilled-taper-leaves-base",
+        "b94",
+        3,
+        1.2339,
+        false,
+    ),
+    (
+        "crashes/2243-step-export-flanged-enclosure",
+        "b43",
+        1,
+        1.4008,
+        false,
+    ),
+    ("fillets/2903-original-rim-chamfer", "b0", 1, 1.5243, false),
+    ("fillets/2903-original-rim-chamfer", "b0", 3, 2.0275, false),
+    ("fillets/2903-original-rim-round", "b0", 1, 1.5243, false),
+    ("fillets/2903-original-rim-round", "b0", 3, 2.0275, false),
+];
+
+/// Each corpus label's margin, pinned to 1 % of what #3185 measured, on its side of 1, and
+/// `build` agreeing with it.
+#[test]
+fn the_corpus_labels_stand_or_fall_by_the_derived_bound() {
+    let k = kernel();
+    let mut report = String::new();
+    let mut wrong = Vec::new();
+    for &(case, id, face, pinned, stands) in CASES {
+        let (support, form) = label_and_form(&k, case, id, face);
+        let m = support.parametrisation_margin(&k, &form).unwrap().unwrap();
+        let _ = writeln!(report, "{case} {id} face {face}: margin {m:.4}");
+        let built = support.build(&k, form).is_ok();
+        if (m <= 1.0) != stands || built != stands || (m - pinned).abs() > 0.01 * pinned {
+            wrong.push(format!(
+                "{case} {id} face {face}: margin {m} (pinned {pinned}, stands {stands}), build {built}"
+            ));
+        }
+    }
+    eprintln!("{report}");
+    assert!(wrong.is_empty(), "{wrong:#?}");
+}
+
+/// A torus at the groove's sizes (#3185's b46: centre (19.375, -33.558, 37.5), R 33.657, r 2,
+/// a 43.6 degree sweep of the whole tube), made by the kernel: its own label stands on its
+/// form; turned about its axis by four times the turn the bound allows (the margin is linear
+/// in a small turn: a turn of `t` reads `m(t)`, so `4 t / m(t)` reads 4), it is refused.
+#[test]
+fn a_torus_label_turned_past_its_bound_is_refused() {
+    let k = kernel();
+    let o = Point3::new(19.374_999_999_999_975, -33.558_484_396_647_01, 37.5);
+    let x = Vec3::new(
+        -0.427_099_277_807_218_7,
+        0.739_757_649_038_078_7,
+        -0.519_946_946_895_745_4,
+    );
+    let y = Vec3::new(
+        0.259_973_473_447_872_36,
+        -0.450_287_264_631_874_2,
+        -0.854_198_555_614_438_5,
+    );
+    let domain = Extent2d::from_bounds(0.0, -180.0, 43.602_818_972_703_645, 180.0);
+    let made = Torus::new(
+        &k,
+        o,
+        x,
+        y,
+        domain,
+        33.657_280_044_590_65,
+        2.0,
+        TorusOptions::default(),
+    )
+    .unwrap();
+    let s = Surface::from(made);
+    let form = as_b_spline(&k, &s).unwrap().unwrap().into_owned();
+    let own = AnalyticSupport::of(&s).unwrap();
+    let m0 = own.parametrisation_margin(&k, &form).unwrap().unwrap();
+    assert!(m0 <= 1.0, "own label: {m0}");
+    assert!(own.build(&k, form.copy()).is_ok());
+    let turned = |t: f64| {
+        let AnalyticSupport::Torus {
+            placement,
+            major_radius,
+            minor_radius,
+        } = own
+        else {
+            panic!()
+        };
+        let (c, sn) = (t.cos(), t.sin());
+        let (px, py) = (placement.x_axis(), placement.y_axis());
+        AnalyticSupport::Torus {
+            placement: Axis2Placement::from_stored(
+                placement.origin(),
+                px * c + py * sn,
+                py * c - px * sn,
+            ),
+            major_radius,
+            minor_radius,
+        }
+    };
+    let probe = 1e-11;
+    let m1 = turned(probe)
+        .parametrisation_margin(&k, &form)
+        .unwrap()
+        .unwrap();
+    assert!(m1 > 1.0, "a turn of {probe} rad reads {m1}");
+    let four = 4.0 * probe / m1;
+    let m4 = turned(four)
+        .parametrisation_margin(&k, &form)
+        .unwrap()
+        .unwrap();
+    eprintln!("own {m0:.4}; turn {probe:e}: {m1:.3}; turn {four:e}: {m4:.3}");
+    assert!(m4 > 1.0, "a turn of {four:e} rad reads {m4}");
+    assert!(turned(four).build(&k, form).is_err());
+}

solidcore/crates/solidcore-ops/src/contact.rs

@@ -0,0 +1,140 @@
+//! Contact between bodies: whether two bodies meet (share material or touch within a
+//! tolerance) and the distance between their boundaries, each decided by the kernel's
+//! certified bracket ([`solidcore::query::contact`]). Ops passes the tolerance on as the
+//! caller gives it and reads the kernel's verdict; it computes no distance of its own.
+
+use solidcore::geom::certified::Verdict;
+use solidcore::kernel::Kernel;
+use solidcore::query::contact::{self as c, Contact, Prepared};
+use solidcore::topology::brep::p3::brep_copy;
+use solidcore::topology::Brep;
+use solidcore::Point3;
+
+use crate::body::Body;
+use crate::session::{run_on, Answer};
+use crate::{Error, Result};
+
+impl Answer for Contact {}
+impl Answer for Vec<Contact> {}
+
+/// The bracket of the distance between two bodies' boundaries, at most the width asked.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct Gap {
+    /// No point of either boundary is nearer the other than this.
+    pub lower: f64,
+    /// The distance of the witness points, rounded up.
+    pub upper: f64,
+    /// A point on each body `upper` apart, when one is known (the first on this body).
+    pub witness: Option<[[f64; 3]; 2]>,
+}
+
+fn array(p: Point3) -> [f64; 3] {
+    [p.x, p.y, p.z]
+}
+
+/// A copy of `body`'s kernel body, the body released at once: the contact test takes its
+/// operands mutably (for its classifications' caches) and never edits the caller's, and a
+/// body tested against itself is locked only once at a time.
+fn copy_of(k: &mut Kernel, body: &Body) -> solidcore::Result<Brep> {
+    let b = body.lock();
+    brep_copy(k, &b)
+}
+
+/// The verdict as an answer: `Yes` meets, `No` does not, `Undecided` refused with the
+/// bracket, which names the tolerance that would decide.
+fn decided(contact: &Contact, other: usize) -> Result<bool> {
+    match contact.verdict {
+        Verdict::Yes => Ok(true),
+        Verdict::No => Ok(false),
+        Verdict::Undecided => Err(Error::ContactUndecided {
+            other,
+            lower: contact.lower,
+            upper: contact.upper,
+        }),
+    }
+}
+
+fn poisoned(bodies: &[&Body]) -> Result<()> {
+    if bodies.iter().any(|b| b.is_poisoned()) {
+        return Err(Error::Internal(
+            "contact: a body was left unfinished by an earlier panic".to_owned(),
+        ));
+    }
+    Ok(())
+}
+
+impl Body {
+    /// Whether this body and `other` meet at `tol` (the caller's: ops invents none): share
+    /// material, or have boundaries within `tol` of each other
+    /// ([`solidcore::query::contact::meets`]). Neither body changes; `other` may be this body.
+    ///
+    /// # Errors
+    /// [`Error::ContactUndecided`] when the kernel's bracket does not decide it (`other` 0):
+    /// the bracket, never read as either answer; [`Error::Status`] `SM_ERR_INVALID_INPUT` for a
+    /// tolerance not finite and positive or a body without faces; those of the session.
+    pub fn meets(&self, other: &Body, tol: f64) -> Result<bool> {
+        poisoned(&[self, other])?;
+        let contact = run_on("contact::meets", &[("a", self), ("b", other)], |k| {
+            let mut a = copy_of(k, self)?;
+            let mut b = copy_of(k, other)?;
+            c::meets(k, &mut a, &mut b, tol)
+        })?;
+        decided(&contact, 0)
+    }
+
+    /// [`Body::meets`] against each of `others`, each with its own tolerance, this body
+    /// prepared once (its faces' pieces about the centre of its box,
+    /// [`Prepared::about`]).
+    ///
+    /// # Errors
+    /// [`Error::ContactUndecided`] for the first of `others` the kernel's bracket does not
+    /// decide, by its index; those of [`Body::meets`].
+    pub fn meets_each(&self, others: &[(&Body, f64)]) -> Result<Vec<bool>> {
+        let mut bodies: Vec<&Body> = vec![self];
+        bodies.extend(others.iter().map(|(b, _)| *b));
+        poisoned(&bodies)?;
+        let inputs: Vec<(&'static str, &Body)> = bodies.iter().map(|b| ("brep", *b)).collect();
+        let contacts = run_on("contact::meets_each", &inputs, |k| {
+            let mut a = copy_of(k, self)?;
+            let prepared = Prepared::about(k, &a)?;
+            let mut contacts = Vec::new();
+            for (other, tol) in others {
+                let mut b = copy_of(k, other)?;
+                contacts.push(c::meets_prepared(k, &mut a, &prepared, &mut b, *tol)?);
+            }
+            Ok(contacts)
+        })?;
+        contacts
+            .iter()
+            .enumerate()
+            .map(|(i, contact)| decided(contact, i))
+            .collect()
+    }
+
+    /// The distance between this body's boundary and `other`'s, bracketed to at most `width`
+    /// (relative to its upper end, absolute below 1;
+    /// [`solidcore::query::contact::boundary_distance`]).
+    ///
+    /// # Errors
+    /// [`Error::ContactUndecided`] (`other` 0) with the narrowest bracket reached, when the
+    /// kernel's work runs out before `width`; [`Error::Status`] `SM_ERR_INVALID_INPUT` for a
+    /// width not finite and positive or a body without faces; those of the session.
+    pub fn boundary_distance(&self, other: &Body, width: f64) -> Result<Gap> {
+        poisoned(&[self, other])?;
+        let contact = run_on(
+            "contact::boundary_distance",
+            &[("a", self), ("b", other)],
+            |k| {
+                let mut a = copy_of(k, self)?;
+                let mut b = copy_of(k, other)?;
+                c::boundary_distance(k, &mut a, &mut b, width)
+            },
+        )?;
+        decided(&contact, 0)?;
+        Ok(Gap {
+            lower: contact.lower,
+            upper: contact.upper,
+            witness: contact.witness.map(|[p, q]| [array(p), array(q)]),
+        })
+    }
+}

solidcore/crates/solidcore-ops/src/error.rs

@@ -134,6 +134,19 @@ pub enum Error {
     /// The kernel found its own state inconsistent part way through a call, or was called
     /// before the preparation it needs.
     Exception,
+    /// A contact test the kernel's certified bracket of the boundary distance does not decide
+    /// (`solidcore::query::contact`): the tolerance lies inside `[lower, upper]` (or the
+    /// boundaries are apart but one body's enclosure of the other is not decided), or a
+    /// distance's bracket is still wider than asked. Never read as either answer. `other` is
+    /// the index of the body tested among the call's others.
+    ContactUndecided {
+        /// The other body's index among the call's (0 for one).
+        other: usize,
+        /// No point of either boundary is nearer the other than this.
+        lower: f64,
+        /// A pair of points, one on each boundary, is this far apart.
+        upper: f64,
+    },
     /// The application cancelled the call.
     Cancelled,
     /// Ops was called from inside a kernel call or from the failure callback.
@@ -325,6 +338,14 @@ impl std::fmt::Display for Error {
             ),
             Error::TooFar { largest } => write!(f, "too far to extend (largest {largest:?})"),
             Error::Exception => write!(f, "solidcore raised an internal exception"),
+            Error::ContactUndecided {
+                other,
+                lower,
+                upper,
+            } => write!(
+                f,
+                "contact with body {other} is undecided: the boundary distance is in [{lower}, {upper}]"
+            ),
             Error::Cancelled => write!(f, "cancelled"),
             Error::Reentered => write!(
                 f,

solidcore/crates/solidcore-ops/src/lib.rs

@@ -19,8 +19,8 @@
 //! - [`provenance`]: the hook an application's construction history rides on.
 //!
 //! The families: [`measure`] (mass properties, area, box, point class, gap samples, the picks'
-//! distances and normals), [`boolean`], [`stitch`], [`prism`] (planar region booleans and
-//! extrusion, prism sections).
+//! distances and normals), [`contact`] (whether bodies meet, their boundary distance),
+//! [`boolean`], [`stitch`], [`prism`] (planar region booleans and extrusion, prism sections).
 //!
 //! Operations change the body they are given, as the kernel does. An application that keeps
 //! its input copies it first ([`Body::copy`]).
@@ -32,6 +32,7 @@ pub mod blend;
 pub mod body;
 pub mod boolean;
 pub mod check;
+pub mod contact;
 pub mod error;
 pub mod extend;
 pub mod file;
@@ -56,6 +57,7 @@ pub use blend::{
 pub use body::{make_box, make_cone, make_cylinder, make_sphere, make_torus, Body};
 pub use boolean::BooleanOp;
 pub use check::{CheckLevel, CheckReport};
+pub use contact::Gap;
 pub use error::{Error, FitKind, Result, RoomReason, Unblendable};
 pub use measure::{
     point_distance, Distance, EdgePoint, FaceCentroid, GapSample, MassPropertiesEstimate, Normal,
@@ -65,6 +67,8 @@ pub use offset::{draft, offset_faces, offset_sheet, shell, thicken};
 pub use order::Order;
 pub use step_assembly::{placed_bodies, PlacedBody};
 
+#[cfg(test)]
+mod noisy_round_tests;
 #[cfg(test)]
 mod test_corpus;
 #[cfg(test)]

solidcore/crates/solidcore-ops/src/noisy_round_tests.rs

@@ -0,0 +1,199 @@
+//! #3066: a round on supports read from a record with rounding in their frames keeps its exact
+//! definition. The class is decided within the body's tolerance; the round lies on a torus of
+//! its own, its net certified there at bound 0, and that torus's distance from the class's
+//! torus from the supports is held to the body's tolerance, so the round answers
+//! `procedural_blend()`.
+#![allow(
+    clippy::unwrap_used,
+    clippy::expect_used,
+    clippy::panic,
+    clippy::indexing_slicing,
+    clippy::disallowed_methods
+)] // test
+
+use std::f64::consts::PI;
+
+use solidcore::surface::surface::SurfaceKind;
+use solidcore::topology::face::p2 as face;
+use solidcore_record::brep::Brep as Record;
+
+use crate::blend::{blend, Blend, Corners};
+use crate::boolean::boolean;
+use crate::{make_box, make_cylinder, session, Body, BooleanOp, CheckLevel};
+
+fn fixes() {
+    let fixes: Vec<String> = serde_json::from_str(include_str!(
+        "../../solidcore/tests/kernel_cases/fixes.json"
+    ))
+    .unwrap();
+    session::set_fixes(fixes.iter().map(String::as_str)).unwrap();
+}
+
+/// The faces of `b` that are no plane, as (with a definition, without one).
+fn definitions(b: &Body) -> (usize, usize) {
+    let brep = b.lock();
+    let (mut with, mut without) = (0, 0);
+    for &f in brep.face_list() {
+        let s = face::surface(&brep, f).unwrap().unwrap();
+        if matches!(s.kind(), SurfaceKind::Plane(_)) {
+            continue;
+        }
+        if s.procedural_blend().is_some() {
+            with += 1;
+        } else {
+            without += 1;
+        }
+    }
+    (with, without)
+}
+
+/// The edge of `b` through `p`.
+fn edge_through(b: &Body, p: [f64; 3]) -> usize {
+    let k = solidcore::Kernel::default();
+    let brep = b.lock();
+    let q = solidcore::Vec3::new(p[0], p[1], p[2]);
+    brep.edge_list()
+        .iter()
+        .position(|&e| {
+            solidcore::query::api_queries::edge_closest_point(&k, &brep, e, q)
+                .unwrap()
+                .closest
+                .distance
+                < 1e-6
+        })
+        .unwrap()
+}
+
+/// A boss R 8 standing 10 out of a 40 x 40 x 10 block, its axis turned 1e-10 off square to the
+/// block's top, its root rounded 2: a circle-class round on supports only that close to the
+/// class (the error the 1e-10 turn makes over the part is far inside the body's tolerance, and
+/// far outside the rounding the net is certified to). The round is concave: the corner between
+/// the wall and the top, less the quarter disc, turned about the axis. By Pappus, it adds
+/// 2 pi (R + x) r^2 (1 - pi / 4), x = r (10 - 3 pi) / (12 - 3 pi) the corner's centroid from the
+/// wall.
+#[test]
+fn a_round_on_a_boss_tilted_off_the_class_keeps_its_definition() {
+    fixes();
+    let block = make_box([0.0; 3], [40.0, 40.0, 10.0]).unwrap();
+    // the boss from z = 5 to 20, its axis turned 1e-10 about x off square to the block's top
+    let boss = make_cylinder([0.0; 3], 8.0, 15.0).unwrap();
+    let (c, s) = (1e-10_f64.cos(), 1e-10_f64.sin());
+    boss.transform(&[
+        [1.0, 0.0, 0.0, 0.0],
+        [0.0, c, s, 0.0],
+        [0.0, -s, c, 0.0],
+        [20.0, 20.0, 5.0, 1.0],
+    ])
+    .unwrap();
+    let b = boolean(block, boss, BooleanOp::Union).unwrap().unwrap();
+    let e = edge_through(&b, [28.0, 20.0, 10.0]);
+    let (r, big_r) = (2.0, 8.0);
+    blend(&b, Blend::Round, &[e], &[r], &[], Corners::Filled).unwrap();
+    let report = b.check(CheckLevel::All).unwrap();
+    assert!(report.is_valid(), "{}", report.text);
+    // the boss's wall and the round: the round keeps its definition
+    assert_eq!(definitions(&b), (1, 1));
+    let x = r * (10.0 - 3.0 * PI) / (12.0 - 3.0 * PI);
+    let want =
+        16000.0 + PI * big_r * big_r * 10.0 + 2.0 * PI * (big_r + x) * r * r * (1.0 - PI / 4.0);
+    let got = b.volume(1e-10).unwrap();
+    assert!((got - want).abs() <= 1e-9 * want, "{got} {want}");
+}
+
+/// SolidCraft's raked wedge (`tests/data/raked_wedge_noisy_cylinders.json`): its end outline's
+/// rounds keep their definitions, the two circle-class rounds on the read-back cylinders among
+/// them; the operand's own faces, read from the record, have none.
+#[test]
+fn the_raked_wedges_rounds_keep_their_definitions() {
+    fixes();
+    let record: Record = serde_json::from_str(include_str!(
+        "../tests/data/raked_wedge_noisy_cylinders.json"
+    ))
+    .unwrap();
+    let (b, _) = Body::from_record(&record, true).unwrap();
+    let (_, before) = definitions(&b);
+    let top = 60.0 - 60.0 * 10f64.to_radians().tan();
+    let mut seeds = Vec::new();
+    for x in [-69.0, 69.0] {
+        for p in [
+            [x, -41.0, 22.0],
+            [x, 8.0, 41.0],
+            [x, (38.0 + top) / 2.0, 60.0],
+        ] {
+            seeds.push(edge_through(&b, p));
+        }
+    }
+    let chain = crate::blend::tangent_chain(&b, &seeds, 0.5_f64.to_radians()).unwrap();
+    blend(&b, Blend::Round, &chain, &[0.854], &[], Corners::Filled).unwrap();
+    let report = b.check(CheckLevel::All).unwrap();
+    assert!(report.is_valid(), "{}", report.text);
+    let (with, without) = definitions(&b);
+    // as many of the new rounds keep their definitions as on the same wedge with exact
+    // cylinders, and the rest are as bare there
+    let (exact_with, exact_bare) = exact_wedge_rounds();
+    assert_eq!((with, without - before), (exact_with, exact_bare));
+    assert_eq!(with + without - before, chain.len());
+}
+
+/// The same wedge built here, its cylinders exact, its end outline rounded as SolidCraft's: of the
+/// rounds that last call makes, how many keep a definition and how many don't.
+fn exact_wedge_rounds() -> (usize, usize) {
+    use crate::sweep::{extrude, Curve, DraftCorner};
+    let top = 60.0 - 60.0 * 10f64.to_radians().tan();
+    let profile = [
+        [-60.0, 0.0],
+        [60.0, 0.0],
+        [top, 60.0],
+        [38.0, 60.0],
+        [-22.0, 22.0],
+        [-60.0, 22.0],
+    ];
+    let at = |p: [f64; 2]| [-69.0, p[0], p[1]];
+    let curves: Vec<Curve> = (0..6)
+        .map(|i| Curve::Line {
+            start: at(profile[i]),
+            end: at(profile[(i + 1) % 6]),
+        })
+        .collect();
+    let b = extrude(
+        &curves,
+        [1.0, 0.0, 0.0],
+        138.0,
+        0.0,
+        DraftCorner::Curvature,
+        true,
+    )
+    .unwrap();
+    let round = |pts: &[[f64; 3]], r: f64, chain: bool| {
+        let seeds: Vec<usize> = pts.iter().map(|&p| edge_through(&b, p)).collect();
+        let edges = if chain {
+            crate::blend::tangent_chain(&b, &seeds, 0.5_f64.to_radians()).unwrap()
+        } else {
+            seeds
+        };
+        blend(&b, Blend::Round, &edges, &[r], &[], Corners::Filled).unwrap();
+    };
+    round(&[[0.0, -60.0, 22.0]], 1.2, true);
+    round(&[[0.0, -22.0, 22.0]], 1.5, true);
+    round(&[[0.0, 38.0, 60.0]], 1.5, true);
+    let back = (60.0 + top) / 2.0;
+    let four = [
+        [-69.0, -60.0, 12.0],
+        [69.0, -60.0, 12.0],
+        [-69.0, back, 30.0],
+        [69.0, back, 30.0],
+    ];
+    round(&four, 4.0, false);
+    let (with_before, bare_before) = definitions(&b);
+    let mut ends = Vec::new();
+    for x in [-69.0, 69.0] {
+        ends.extend([
+            [x, -41.0, 22.0],
+            [x, 8.0, 41.0],
+            [x, (38.0 + top) / 2.0, 60.0],
+        ]);
+    }
+    round(&ends, 0.854, true);
+    let (with, bare) = definitions(&b);
+    (with - with_before, bare - bare_before)
+}

solidcore/crates/solidcore-ops/tests/contact.rs

@@ -0,0 +1,95 @@
+//! Contact between bodies on closed forms: meets, meets against several, the boundary
+//! distance, and the refusals.
+#![allow(clippy::unwrap_used)]
+
+use solidcore_ops::session;
+use solidcore_ops::{make_box, make_sphere, Error};
+
+fn all_fixes() {
+    session::set_thread_fixes(solidcore::kernel::UpstreamFixes::NAMES.iter().copied()).unwrap();
+}
+
+#[test]
+fn boxes_sharing_a_face_meet() {
+    all_fixes();
+    let a = make_box([0.0, 0.0, 0.0], [1.0, 1.0, 1.0]).unwrap();
+    let b = make_box([1.0, 0.0, 0.0], [1.0, 1.0, 1.0]).unwrap();
+    assert!(a.meets(&b, 1e-6).unwrap());
+    assert!(b.meets(&a, 1e-6).unwrap());
+}
+
+#[test]
+fn boxes_one_apart_do_not_meet_and_their_distance_is_one() {
+    all_fixes();
+    let a = make_box([0.0, 0.0, 0.0], [1.0, 1.0, 1.0]).unwrap();
+    let b = make_box([2.0, 0.0, 0.0], [1.0, 1.0, 1.0]).unwrap();
+    assert!(!a.meets(&b, 1e-3).unwrap());
+    let g = a.boundary_distance(&b, 1e-9).unwrap();
+    assert!(g.lower <= 1.0 && 1.0 <= g.upper, "{g:?}");
+    assert!(g.upper - g.lower <= 1e-9 * g.upper.max(1.0), "{g:?}");
+    let [p, q] = g.witness.unwrap();
+    assert!(
+        (p[0] - 1.0).abs() <= 1e-12 && (q[0] - 2.0).abs() <= 1e-12,
+        "{g:?}"
+    );
+}
+
+#[test]
+fn a_body_meets_itself_and_one_it_encloses() {
+    all_fixes();
+    let a = make_box([0.0, 0.0, 0.0], [3.0, 3.0, 3.0]).unwrap();
+    let inner = make_box([1.0, 1.0, 1.0], [1.0, 1.0, 1.0]).unwrap();
+    assert!(a.meets(&a, 1e-6).unwrap());
+    assert!(a.meets(&inner, 1e-6).unwrap());
+    assert!(inner.meets(&a, 1e-6).unwrap());
+    // neither body changes
+    assert_eq!(a.counts().unwrap(), [6, 12, 8]);
+    assert_eq!(inner.counts().unwrap(), [6, 12, 8]);
+}
+
+#[test]
+fn meets_each_answers_each_with_its_own_tolerance() {
+    all_fixes();
+    let a = make_sphere([0.0, 0.0, 0.0], 1.0).unwrap();
+    let touching = make_box([1.0, -1.0, -1.0], [1.0, 2.0, 2.0]).unwrap();
+    let near = make_sphere([2.5, 0.0, 0.0], 1.0).unwrap(); // 0.5 apart
+    let far = make_box([5.0, 0.0, 0.0], [1.0, 1.0, 1.0]).unwrap();
+    let got = a
+        .meets_each(&[(&touching, 1e-6), (&near, 0.1), (&near, 0.6), (&far, 1.0)])
+        .unwrap();
+    assert_eq!(got, vec![true, false, true, false]);
+}
+
+/// The refusal is the kernel's own bracket: a tolerance equal to a gap the kernel brackets
+/// to its rounding is inside the bracket, and neither answer is given.
+#[test]
+fn a_tolerance_inside_the_bracket_is_refused_with_it() {
+    all_fixes();
+    let a = make_box([0.0, 0.0, 0.0], [1.0, 1.0, 1.0]).unwrap();
+    let b = make_box([2.0, 0.0, 0.0], [1.0, 1.0, 1.0]).unwrap();
+    let far = make_box([5.0, 0.0, 0.0], [1.0, 1.0, 1.0]).unwrap();
+    match a.meets(&b, 1.0) {
+        Err(Error::ContactUndecided {
+            other: 0,
+            lower,
+            upper,
+        }) => assert!(lower <= 1.0 && 1.0 <= upper, "[{lower}, {upper}]"),
+        r => panic!("{r:?}"),
+    }
+    match a.meets_each(&[(&far, 1e-3), (&b, 1.0)]) {
+        Err(Error::ContactUndecided { other: 1, .. }) => {}
+        r => panic!("{r:?}"),
+    }
+}
+
+#[test]
+fn a_tolerance_or_width_not_finite_and_positive_is_refused() {
+    all_fixes();
+    let a = make_box([0.0, 0.0, 0.0], [1.0, 1.0, 1.0]).unwrap();
+    let b = make_box([2.0, 0.0, 0.0], [1.0, 1.0, 1.0]).unwrap();
+    let invalid = Error::Status(solidcore::error::SM_ERR_INVALID_INPUT);
+    for bad in [0.0, -1.0, f64::NAN, f64::INFINITY] {
+        assert_eq!(a.meets(&b, bad).unwrap_err(), invalid, "{bad}");
+        assert_eq!(a.boundary_distance(&b, bad).unwrap_err(), invalid, "{bad}");
+    }
+}

solidcore/crates/solidcore/src/blend/fillet_exact.rs

@@ -2521,46 +2521,11 @@ fn keep_definition(
     let Some(surface) = pseudo.geom_mut(g)?.fillet_surface_mut(None)? else {
         return Err(SolidError::null_pointer(OP));
     };
-    keep_round_if_proved(k, surface, round, inputs.tol)?;
+    round.recertify(k, surface)?;
+    surface.set_exact_round(Some(round));
     Ok(())
 }
 
-/// `round` kept on `net` where it is proved ([`ExactRound::recertify`]), returning whether it
-/// was. The class is read within `tol`, the body's tolerance, and the certificate to rounding.
-/// Supports only that close to the class (a cylinder's axis 1e-13 off square to a plane, as
-/// read from a record) make a net that is the round but not provably on the class's surface:
-/// no bound of 0, so it keeps no definition. A net further off the
-/// class than `tol` ([`ExactRound::miss`]) is built wrongly, and stays an error.
-///
-/// # Errors
-///
-/// The certificate's refusal of a net more than `tol` off the class, or of supports that
-/// give the class no surface; those of the certificate and of [`ExactRound::miss`].
-pub fn keep_round_if_proved(
-    k: &Kernel,
-    net: &mut Surface,
-    round: ExactRound,
-    tol: f64,
-) -> Result<bool> {
-    match round.recertify(k, net) {
-        Ok(_) => {
-            net.set_exact_round(Some(round));
-            Ok(true)
-        }
-        Err(e) => {
-            let off = matches!(
-                &e,
-                SolidError::InvalidInput { operation, .. }
-                    if *operation == crate::surface::implicit::OP_EXACT
-            );
-            match off.then(|| round.miss(k, net)).transpose()?.flatten() {
-                Some(miss) if miss <= tol => Ok(false),
-                _ => Err(e),
-            }
-        }
-    }
-}
-
 /// [`calc_exact_fillet_geom`]'s build by the class (`CalcExactFilletGeom`), without the
 /// definition kept.
 #[allow(clippy::too_many_arguments)]
@@ -2719,7 +2684,8 @@ fn exact_fillet_geom_of_class(
 mod tests {
     use super::*;
 
-    /// A box's edge rounded 2: the round's net, without its definition, and the definition.
+    /// A cylinder R 10, 20 high, its top rim rounded 2 (a round of the circle class): the
+    /// round's net, without its definition, and the definition.
     fn kept_round() -> (Kernel, Surface, ExactRound) {
         use crate::blend::fillet_executive::session::FilletSession;
         use crate::kernel::UpstreamFixes;
@@ -2729,24 +2695,51 @@ mod tests {
             fixes: UpstreamFixes::all(),
             ..Kernel::default()
         };
-        let mut b = crate::primitive::api_primitives::create_box(
+        let mut b = crate::primitive::api_primitives::create_cylinder(
             &mut k,
             Point3::new(0.0, 0.0, 0.0),
-            30.0,
-            20.0,
             10.0,
+            20.0,
         )
         .unwrap();
-        let edge = b.edge_list()[0];
+        // the top rim: the edge through (10, 0, 20)
+        let rim = *b
+            .edge_list()
+            .iter()
+            .min_by(|&&x, &&y| {
+                let d = |e| {
+                    crate::query::api_queries::edge_closest_point(
+                        &k,
+                        &b,
+                        e,
+                        Point3::new(10.0, 0.0, 20.0),
+                    )
+                    .unwrap()
+                    .closest
+                    .distance
+                        + crate::query::api_queries::edge_closest_point(
+                            &k,
+                            &b,
+                            e,
+                            Point3::new(-10.0, 0.0, 20.0),
+                        )
+                        .unwrap()
+                        .closest
+                        .distance
+                };
+                d(x).total_cmp(&d(y))
+            })
+            .unwrap();
         {
             let mut s = FilletSession::new(&mut k, &mut b).unwrap();
-            s.prepare_constant_radius_edges(&mut k, &[edge], 2.0, 1, 1, 1.0)
+            s.prepare_constant_radius_edges(&mut k, &[rim], 2.0, 1, 1, 1.0)
                 .unwrap();
             s.do_fillet(&mut k).unwrap();
         }
         for &f in b.face_list() {
             let s = face::surface(&b, f).unwrap().unwrap();
             if let Some(BlendRef::Exact { round, .. }) = s.procedural_blend() {
+                assert_eq!(round.class, ExactClass::Circle);
                 let round = round.clone();
                 let mut net = s.clone();
                 net.set_exact_round(None);
@@ -2769,23 +2762,21 @@ mod tests {
         net
     }
 
-    /// The class is read within 2e-5 here. A net on it keeps its definition; one off it by
-    /// less keeps none (the round of supports that close to the class); one further off is
-    /// built wrongly, and stays an error.
+    /// A circle-class round's definition is the torus its net is on, a stated distance from
+    /// the supports' class torus within the tolerance the class was read at (2e-5 here): the
+    /// net as built is at 0, the net moved by (1, 1, 1) 1e-7 is at sqrt 3 1e-7, and the net
+    /// moved by (1, 1, 1) 1e-3 is refused.
     #[test]
-    fn a_net_off_its_class_keeps_no_definition_within_the_tolerance_only() {
-        let tol = 2e-5;
+    fn a_circle_rounds_net_is_certified_on_its_own_torus_near_the_class() {
         let (k, net, round) = kept_round();
-        let mut on = net.clone();
-        assert!(keep_round_if_proved(&k, &mut on, round.clone(), tol).unwrap());
-        assert!(on.procedural_blend().is_some());
-
-        let mut near = moved(&k, &net, 1e-7);
-        assert!(!keep_round_if_proved(&k, &mut near, round.clone(), tol).unwrap());
-        assert!(near.procedural_blend().is_none());
-
-        let mut far = moved(&k, &net, 1e-3);
-        let refused = keep_round_if_proved(&k, &mut far, round, tol);
+        assert_eq!(round.recertify(&k, &net).unwrap(), 0.0);
+        assert!(round.supports_distance(&k, &net).unwrap() < 1e-12);
+        let near = moved(&k, &net, 1e-7);
+        assert_eq!(round.recertify(&k, &near).unwrap(), 0.0);
+        let d = round.supports_distance(&k, &near).unwrap();
+        assert!((d - 3.0_f64.sqrt() * 1e-7).abs() < 1e-12, "{d}");
+        let far = moved(&k, &net, 1e-3);
+        let refused = round.recertify(&k, &far);
         assert!(
             matches!(refused, Err(SolidError::InvalidInput { .. })),
             "{refused:?}"

solidcore/crates/solidcore/src/blend/fillet_exact/rational.rs

@@ -1304,7 +1304,7 @@ fn rational_spine_parts(
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::blend::fillet_exact::{classify_supports, keep_round_if_proved, ClassInputs};
+    use crate::blend::fillet_exact::{classify_supports, ClassInputs};
     use crate::curve::curve::CurveKind;
     use crate::geom::axis2_placement::Axis2Placement;
     use crate::surface::cone::Cone;
@@ -1779,10 +1779,10 @@ mod tests {
                             inputs,
                         )
                         .unwrap();
-                        assert!(
-                            keep_round_if_proved(&k, &mut net, round, inputs.tol).unwrap(),
-                            "{offsets:?}, {contact:?}"
-                        );
+                        if let Err(e) = round.recertify(&k, &net) {
+                            panic!("{offsets:?}, {contact:?}: {e}");
+                        }
+                        net.set_exact_round(Some(round));
                         assert_eq!(
                             net.procedural_blend().unwrap().exact().unwrap().bound(),
                             0.0

solidcore/crates/solidcore/src/query.rs

@@ -6,6 +6,7 @@
 
 pub mod api_queries;
 pub mod brep;
+pub mod contact;
 pub mod distance_bound;
 pub mod edge;
 pub mod edge_turn;

solidcore/crates/solidcore/src/query/contact.rs

@@ -0,0 +1,732 @@
+//! Not upstream: whether two bodies meet (share material or touch within a tolerance),
+//! decided only where it is proved: [`meets`] and its distance, [`boundary_distance`].
+//!
+//! | Rust | Upstream |
+//! |------|----------|
+//! | [`meets`], [`boundary_distance`], [`Prepared`] | none (`SmApiBrepRelationship`, ported as [`crate::query::api_queries::brep_relationship`], is a Newton solve: its gap is a point pair, never a lower bound, and its containment counts samples) |
+//!
+//! **The bracket.** The distance between the two bodies' boundaries is enclosed in `[lower,
+//! upper]`:
+//!
+//! - `lower`: the least, over pairs of pieces of the two bodies' faces, of the distance
+//!   between their hulls' boxes. The pieces are those of
+//!   [`crate::query::distance_bound`] (each face's NURBS form over its UV box, its Bézier
+//!   patches about one common point, every coefficient's error carried and each box widened
+//!   by it), so no point of either boundary is nearer the other than `lower`. Pairs wait on a
+//!   heap, nearest first; the nearest is split (its larger piece, in four), so `lower`
+//!   rises. A piece whose UV rectangle is proved off its face (no trim curve's box meets it,
+//!   and its centre classifies outside the face) is dropped: a face's hole is no material.
+//! - `upper` and its witness: a pair of points, one on each body, each on its face's surface
+//!   at parameters the face classifies decisively inside (`face_classify_uv`, the face
+//!   itself: neither on its boundary nor undecided), their distance rounded up. The
+//!   candidates are the centres of the pairs taken from the heap, each moved to the other
+//!   surface's closest point and back ([`crate::query::api_queries::surface_closest_point`],
+//!   a few times), so two surfaces that cross reach a distance of their rounding.
+//!
+//! **Meets.** At a tolerance `tol`: `Yes` when `upper <= tol` (with the witness); when
+//! `lower > tol`, the boundaries do not touch, and the bodies meet only if one encloses the
+//! other: a point of each boundary is classified in the other body, `Yes` when one is
+//! inside its material, `No` when both are outside. A classification on the other's boundary
+//! or not decided is no answer, and neither is a bracket that still straddles `tol` when
+//! the pairs reach [`MAX_PAIRS`]: [`Verdict::Undecided`], with the bracket, which names
+//! the size that would decide.
+//!
+//! Deviations: none (no upstream counterpart). No recursion: the pairs wait on a heap.
+
+use std::cmp::Ordering;
+use std::collections::BinaryHeap;
+
+use crate::geom::certified::{gamma, Verdict};
+use crate::geom::extent1d::Extent1d;
+use crate::geom::tol::{brep_zone_tol_3d, SmallTopoType};
+use crate::geom::vector2d::Point2;
+use crate::intersect::point_class::PointClassType;
+use crate::kernel::Kernel;
+use crate::mass::integrator::properties::region_material;
+use crate::math::Point3;
+use crate::query::api_queries::{
+    edge_closest_point, face_bounding_box, face_classify_uv, surface_closest_point,
+};
+use crate::query::brep::BrepClassifier;
+use crate::query::distance_bound::{pieces_about, Piece};
+use crate::surface::surface::Surface;
+use crate::topology::{edgeuse::p1 as eu1, loopuse, Brep, FaceId};
+use crate::{EntityRef, Result, SolidError};
+
+const OP: &str = "query::contact::meets";
+
+/// The most piece pairs made for one answer.
+pub const MAX_PAIRS: usize = 1 << 16;
+
+/// The times a witness is moved to the other surface's closest point and back.
+const PROJECTIONS: usize = 6;
+
+/// How two bodies stand: the verdict at the tolerance asked, and the boundary distance's
+/// bracket.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct Contact {
+    /// `Yes` (they share material or their boundaries come within the tolerance), `No`, or
+    /// `Undecided` (the bracket straddles the tolerance at the work's limit, or a
+    /// classification was not decided).
+    pub verdict: Verdict,
+    /// A certified lower bound of the distance between the boundaries.
+    pub lower: f64,
+    /// An upper bound (the witness's distance, rounded up); infinite without a witness.
+    pub upper: f64,
+    /// Two points, one on each boundary, `upper` apart.
+    pub witness: Option<[Point3; 2]>,
+}
+
+/// A face's pieces and what decides whether a piece is on it.
+struct FacePieces {
+    face: FaceId,
+    surface: Surface,
+    /// Its form's certified distance from its definition, lost from every bound.
+    slack: f64,
+    pieces: Vec<Piece>,
+    /// The boxes of its own UV trim curves `[u0, u1, v0, v1]`; `None` when a trim curve or
+    /// its box is missing (then no piece is dropped).
+    trims: Option<Vec<[f64; 4]>>,
+}
+
+/// A body made ready for contact tests: its faces' pieces about a common point.
+pub struct Prepared {
+    origin: Point3,
+    faces: Vec<FacePieces>,
+}
+
+/// A pair of pieces, one of each body, with its lower bound.
+struct Pair {
+    lower: f64,
+    a: (usize, Piece),
+    b: (usize, Piece),
+}
+
+impl PartialEq for Pair {
+    fn eq(&self, o: &Self) -> bool {
+        self.lower.total_cmp(&o.lower) == Ordering::Equal
+    }
+}
+impl Eq for Pair {}
+impl PartialOrd for Pair {
+    fn partial_cmp(&self, o: &Self) -> Option<Ordering> {
+        Some(self.cmp(o))
+    }
+}
+impl Ord for Pair {
+    // the heap pops the smallest bound first
+    fn cmp(&self, o: &Self) -> Ordering {
+        o.lower.total_cmp(&self.lower)
+    }
+}
+
+/// The box of body `b`'s faces, `(lo, hi)`.
+fn body_box(k: &Kernel, b: &Brep) -> Result<(Point3, Point3)> {
+    let mut lo = Point3::new(f64::INFINITY, f64::INFINITY, f64::INFINITY);
+    let mut hi = Point3::new(f64::NEG_INFINITY, f64::NEG_INFINITY, f64::NEG_INFINITY);
+    for &f in b.face_list() {
+        let e = face_bounding_box(k, b, f, false)?;
+        let (l, h) = (e.min(), e.max());
+        lo = Point3::new(lo.x.min(l.x), lo.y.min(l.y), lo.z.min(l.z));
+        hi = Point3::new(hi.x.max(h.x), hi.y.max(h.y), hi.z.max(h.z));
+    }
+    if !(lo.x <= hi.x && lo.y <= hi.y && lo.z <= hi.z) {
+        return Err(SolidError::invalid_input(OP, "a body without faces"));
+    }
+    Ok((lo, hi))
+}
+
+/// The common point two bodies' pieces are made about: the centre of their boxes.
+fn common_origin(k: &Kernel, a: &Brep, b: &Brep) -> Result<Point3> {
+    let ((al, ah), (bl, bh)) = (body_box(k, a)?, body_box(k, b)?);
+    let mid = |x: f64, y: f64, z: f64, w: f64| 0.5 * (x.min(z) + y.max(w));
+    Ok(Point3::new(
+        mid(al.x, ah.x, bl.x, bh.x),
+        mid(al.y, ah.y, bl.y, bh.y),
+        mid(al.z, ah.z, bl.z, bh.z),
+    ))
+}
+
+/// The parts each Bézier segment of a trim curve is cut into for its boxes: a rectangle
+/// that no part's box meets is proved clear of the curve.
+const TRIM_PARTS: usize = 16;
+
+/// The boxes, `[u0, u1, v0, v1]`, of a UV trim curve's NURBS form: each Bézier segment cut in
+/// [`TRIM_PARTS`] by de Casteljau on its homogeneous points, each part in the box of its
+/// Euclidean points (its hull: positive weights), widened by the subdivision's rounding;
+/// `None` for a form whose weights are not positive.
+fn nurb_boxes(c: &crate::nurbs::crv_basic::Curve) -> Result<Option<Vec<[f64; 4]>>> {
+    use crate::nurbs::math::NOW;
+    let mut out = Vec::new();
+    for seg in c.decompose_bez()? {
+        let pts: Vec<[f64; 3]> = seg
+            .cpts()
+            .iter()
+            .map(|q| {
+                let w = if q.w == NOW { 1.0 } else { q.w };
+                [q.x, q.y, w]
+            })
+            .collect();
+        if pts
+            .iter()
+            .any(|q| !(q[2] > 0.0 && q.iter().all(|v| v.is_finite())))
+        {
+            return Ok(None);
+        }
+        let size = pts.iter().flatten().fold(0.0_f64, |m, v| m.max(v.abs()));
+        let wmin = pts.iter().fold(f64::INFINITY, |m, q| m.min(q[2]));
+        // the right part from t: de Casteljau at t on the homogeneous points
+        let cut = |p: &[[f64; 3]], t: f64| -> (Vec<[f64; 3]>, Vec<[f64; 3]>) {
+            let mut level = p.to_vec();
+            let (mut left, mut right) = (Vec::new(), Vec::new());
+            while let (Some(&f), Some(&l)) = (level.first(), level.last()) {
+                left.push(f);
+                right.push(l);
+                level = level
+                    .windows(2)
+                    .filter_map(|w| match w {
+                        [a, b] => Some([
+                            (1.0 - t) * a[0] + t * b[0],
+                            (1.0 - t) * a[1] + t * b[1],
+                            (1.0 - t) * a[2] + t * b[2],
+                        ]),
+                        _ => None,
+                    })
+                    .collect();
+            }
+            right.reverse();
+            (left, right)
+        };
+        let degree = pts.len().saturating_sub(1);
+        // TRIM_PARTS levels of convex combinations at most: their rounding, of the size
+        let err = gamma(3 * degree * TRIM_PARTS + 3) * size;
+        let mut rest = pts;
+        for i in 0..TRIM_PARTS {
+            let t = 1.0 / (TRIM_PARTS - i) as f64;
+            let (part, next) = if i + 1 == TRIM_PARTS {
+                (rest.clone(), Vec::new())
+            } else {
+                cut(&rest, t)
+            };
+            rest = next;
+            let (mut lo, mut hi) = ([f64::INFINITY; 2], [f64::NEG_INFINITY; 2]);
+            for q in &part {
+                for (k, v) in [q[0] / q[2], q[1] / q[2]].into_iter().enumerate() {
+                    if let (Some(l), Some(h)) = (lo.get_mut(k), hi.get_mut(k)) {
+                        *l = l.min(v);
+                        *h = h.max(v);
+                    }
+                }
+            }
+            // a quotient of coefficients within err: within err (1 + |x|) / (w - err)
+            let widen = err * (1.0 + size / wmin) / (wmin - err).max(f64::MIN_POSITIVE)
+                + gamma(2) * size / wmin;
+            out.push([lo[0] - widen, hi[0] + widen, lo[1] - widen, hi[1] + widen]);
+        }
+    }
+    Ok(Some(out))
+}
+
+/// The boxes of face `f`'s own UV trim curves (each its NURBS form's parts,
+/// [`nurb_boxes`], or the curve's whole box where it has none); `None` where a trim curve or
+/// its box is missing.
+fn trim_boxes(k: &Kernel, b: &Brep, f: FaceId) -> Result<Option<Vec<[f64; 4]>>> {
+    use crate::curve::curve::CurveKind;
+    let Some(fu) = b.face(f)?.faceuse else {
+        return Ok(None);
+    };
+    let mut out = Vec::new();
+    for lu in eu1::faceuse_loopuses(b, fu, false, None)? {
+        for eu in loopuse::edgeuses(b, lu, None, None)? {
+            // its own, or its mate's on the face's other side (the same face's parameters)
+            let Some(c) = eu1::uv_trim_curve(b, eu)?.0 else {
+                return Ok(None);
+            };
+            let nurb = match c.kind() {
+                CurveKind::BSpline(x) => x.gw_nurb(),
+                CurveKind::Ellipse(x) => x.base().gw_nurb(),
+                CurveKind::Hyperbola(x) => x.base().gw_nurb(),
+                CurveKind::Parabola(x) => x.base().gw_nurb(),
+                CurveKind::Line(x) => x.base().gw_nurb(),
+                _ => None,
+            };
+            if let Some(n) = nurb {
+                if let Some(boxes) = nurb_boxes(n)? {
+                    out.extend(boxes);
+                    continue;
+                }
+            }
+            let interval: Extent1d = c.natural_interval();
+            let boxes = c.calculate_bounding_box(
+                k,
+                &interval,
+                crate::nurbs::b_spline_curve::BoundingBoxRequest {
+                    normal: true,
+                    ..Default::default()
+                },
+                true,
+            )?;
+            let Some(e) = boxes.normal else {
+                return Ok(None);
+            };
+            let (l, h) = (e.min(), e.max());
+            // the box's own rounding, a few ulps of its size
+            let r = 16.0 * f64::EPSILON * (l.x.abs() + h.x.abs() + l.y.abs() + h.y.abs());
+            out.push([l.x - r, h.x + r, l.y - r, h.y + r]);
+        }
+    }
+    Ok(Some(out))
+}
+
+impl Prepared {
+    /// Body `b`'s faces' pieces about `origin`.
+    ///
+    /// # Errors
+    ///
+    /// `InvalidInput` for a body without faces or a face without a surface; those of the
+    /// pieces ([`crate::query::distance_bound`]) and the trim curves' boxes.
+    pub fn new(k: &Kernel, b: &Brep, origin: Point3) -> Result<Self> {
+        let mut faces = Vec::new();
+        for &f in b.face_list() {
+            let face = b.face(f)?;
+            let surface = face
+                .surface
+                .as_ref()
+                .ok_or_else(|| SolidError::invalid_input(OP, "a face without a surface"))?
+                .clone();
+            let (pieces, slack) = pieces_about(k, &surface, Some(&face.uv_domain), origin)?;
+            let trims = trim_boxes(k, b, f)?;
+            faces.push(FacePieces {
+                face: f,
+                surface,
+                slack,
+                pieces,
+                trims,
+            });
+        }
+        if faces.is_empty() {
+            return Err(SolidError::invalid_input(OP, "a body without faces"));
+        }
+        Ok(Self { origin, faces })
+    }
+
+    /// Body `b`'s faces' pieces about the centre of its own box: the preparation of one
+    /// body for tests against several ([`meets_prepared`]).
+    ///
+    /// # Errors
+    ///
+    /// Those of [`Prepared::new`] and of the faces' boxes.
+    pub fn about(k: &Kernel, b: &Brep) -> Result<Self> {
+        Self::new(k, b, common_origin(k, b, b)?)
+    }
+}
+
+/// The certified lower bound of the distance between two pieces about the same point: the
+/// larger of the distance of their widened boxes and, along a unit direction `axis` (the
+/// witness's, when there is one), the gap between their hulls' extents
+/// ([`Piece::extent_along`]), each rounded down. The axis steers which bound is tight, never
+/// whether either holds: any direction separates the hulls by at most their distance.
+fn pair_lower(a: &Piece, b: &Piece, axis: Option<[f64; 3]>) -> Result<f64> {
+    let by_box = box_lower(a, b)?;
+    let Some(u) = axis else {
+        return Ok(by_box);
+    };
+    let len = (u[0] * u[0] + u[1] * u[1] + u[2] * u[2]).sqrt() * (1.0 + gamma(4));
+    if !(len > 0.0 && len.is_finite()) {
+        return Ok(by_box);
+    }
+    let ((a0, a1), (b0, b1)) = (a.extent_along(u)?, b.extent_along(u)?);
+    let gap = (b0 - a1).max(a0 - b1).max(0.0);
+    // (the difference rounded down by gamma_1 of its terms, then over the length)
+    let gap = (gap - gamma(1) * (a0.abs() + a1.abs() + b0.abs() + b1.abs())).max(0.0);
+    let by_axis = gap / len * (1.0 - gamma(2));
+    Ok(by_box.max(by_axis))
+}
+
+/// The distance of two pieces' widened boxes, rounded down.
+fn box_lower(a: &Piece, b: &Piece) -> Result<f64> {
+    let ((al, ah), (bl, bh)) = (a.widened_box()?, b.widened_box()?);
+    let mut sum = 0.0_f64;
+    for i in 0..3 {
+        let (Some(&x0), Some(&x1), Some(&y0), Some(&y1)) =
+            (al.get(i), ah.get(i), bl.get(i), bh.get(i))
+        else {
+            return Err(SolidError::sm_err(OP));
+        };
+        let gap = (y0 - x1).max(x0 - y1).max(0.0);
+        // (the difference rounded up by at most gamma_1 of its terms)
+        let gap = (gap - gamma(1) * (x0.abs() + x1.abs() + y0.abs() + y1.abs())).max(0.0);
+        sum += gap * gap;
+    }
+    Ok(sum.sqrt() * (1.0 - gamma(5)))
+}
+
+/// Where a piece stands against its face: `Some(true)` proved on it and decisively inside
+/// at its centre, `Some(false)` proved off it (no trim curve's box meets its rectangle and
+/// its centre is outside), `None` not decided.
+fn piece_on_face(k: &Kernel, b: &mut Brep, fp: &FacePieces, p: &Piece) -> Result<Option<bool>> {
+    let [u0, u1, v0, v1] = p.uv;
+    let crosses = match &fp.trims {
+        None => true,
+        Some(boxes) => boxes
+            .iter()
+            .any(|t| !(t[1] < u0 || t[0] > u1 || t[3] < v0 || t[2] > v1)),
+    };
+    let centre = Point2::new(0.5 * (u0 + u1), 0.5 * (v0 + v1));
+    let class = match face_classify_uv(k, b, fp.face, centre) {
+        Ok(c) => c,
+        // a classification that fails decides nothing
+        Err(e) if e.is_cancelled() => return Err(e),
+        Err(_) => return Ok(None),
+    };
+    Ok(match class {
+        PointClassType::Unknown if !crosses => Some(false),
+        PointClassType::Face => Some(true),
+        _ => None,
+    })
+}
+
+/// A point of body `b` for face `fp` at `uv`: the surface's point when the face classifies
+/// `uv` decisively inside it; when it classifies it on the face's boundary (within its
+/// tolerance, either side), the nearest point of the face's edges to the surface's point,
+/// which is on the body whichever side `uv` was; `None` outside or not decided.
+fn on_face(k: &Kernel, b: &mut Brep, fp: &FacePieces, uv: Point2) -> Result<Option<Point3>> {
+    match face_classify_uv(k, b, fp.face, uv) {
+        Ok(PointClassType::Face) => Ok(Some(fp.surface.evaluate_point(k, uv)?)),
+        Ok(PointClassType::Edge | PointClassType::Vertex) => {
+            let x = fp.surface.evaluate_point(k, uv)?;
+            nearest_edge_point(k, b, fp.face, x)
+        }
+        Err(e) if e.is_cancelled() => Err(e),
+        _ => Ok(None),
+    }
+}
+
+/// The nearest point to `x` of face `f`'s edges (each solved on its edge, so on the body);
+/// `None` when no edge's solve succeeds.
+fn nearest_edge_point(k: &Kernel, b: &Brep, f: FaceId, x: Point3) -> Result<Option<Point3>> {
+    let Some(fu) = b.face(f)?.faceuse else {
+        return Ok(None);
+    };
+    let mut best: Option<(f64, Point3)> = None;
+    for lu in eu1::faceuse_loopuses(b, fu, false, None)? {
+        for eu in loopuse::edgeuses(b, lu, None, None)? {
+            let Some(e) = b.edgeuse(eu)?.edge else {
+                continue;
+            };
+            let Ok(c) = edge_closest_point(k, b, e, x) else {
+                continue;
+            };
+            let d = c.closest.point.distance(x);
+            if best.is_none_or(|(m, _)| d < m) {
+                best = Some((d, c.closest.point));
+            }
+        }
+    }
+    Ok(best.map(|(_, p)| p))
+}
+
+/// The distance of two points rounded up: `|p - q|` within `gamma_5` of the terms' size.
+fn distance_up(p: Point3, q: Point3) -> f64 {
+    let size = p.x.abs() + p.y.abs() + p.z.abs() + q.x.abs() + q.y.abs() + q.z.abs();
+    p.distance(q) * (1.0 + gamma(5)) + gamma(5) * size
+}
+
+/// The closest point of `fp`'s surface to `x`, kept when its face classifies it decisively
+/// inside; `None` for a solve that fails or a point off the face.
+fn closest_on(
+    k: &Kernel,
+    b: &mut Brep,
+    fp: &FacePieces,
+    x: Point3,
+) -> Result<Option<(Point3, Point2)>> {
+    let Ok(c) = surface_closest_point(k, &fp.surface, x) else {
+        return Ok(None);
+    };
+    let Some(&uv) = c.uvs.first() else {
+        return Ok(None);
+    };
+    Ok(on_face(k, b, fp, uv)?.map(|p| (p, uv)))
+}
+
+/// The bracket between `a` and `b` (as [`Prepared`]), refined until `stop` says done or the
+/// pairs reach [`MAX_PAIRS`]: `stop(lower, upper)`.
+fn bracket(
+    k: &Kernel,
+    a: &mut Brep,
+    pa: &Prepared,
+    b: &mut Brep,
+    pb: &Prepared,
+    seed: f64,
+    stop: &dyn Fn(f64, f64) -> bool,
+) -> Result<(f64, f64, Option<[Point3; 2]>, bool)> {
+    let mut heap = BinaryHeap::new();
+    let mut made = 0_usize;
+    for (ia, fa) in pa.faces.iter().enumerate() {
+        for (ib, fb) in pb.faces.iter().enumerate() {
+            for x in &fa.pieces {
+                for y in &fb.pieces {
+                    let lower = (pair_lower(x, y, None)? - fa.slack - fb.slack).max(0.0);
+                    heap.push(Pair {
+                        lower,
+                        a: (ia, x.clone()),
+                        b: (ib, y.clone()),
+                    });
+                    made += 1;
+                }
+            }
+        }
+    }
+    let (mut upper, mut witness) = (f64::INFINITY, None);
+    // witnesses first from every face pair whose bound is within `seed`, nearest first (the
+    // pairs tie closely where two bodies lie along each other; each is visited once here)
+    let mut first: Vec<&Pair> = heap.iter().filter(|p| p.lower <= seed).collect();
+    first.sort_by(|x, y| x.lower.total_cmp(&y.lower));
+    for pair in first {
+        let (Some(fa), Some(fb)) = (pa.faces.get(pair.a.0), pb.faces.get(pair.b.0)) else {
+            return Err(SolidError::sm_err(OP));
+        };
+        if pair.lower < upper {
+            witness_from(
+                k,
+                (&mut *a, fa, &pair.a.1),
+                (&mut *b, fb, &pair.b.1),
+                &mut upper,
+                &mut witness,
+            )?;
+        }
+        if stop(0.0, upper) {
+            break;
+        }
+    }
+    loop {
+        k.check_cancelled(OP)?;
+        let Some(pair) = heap.pop() else {
+            // every piece proved off its face: no boundary is there to bound
+            return Ok((f64::INFINITY, upper, witness, true));
+        };
+        let lower = pair.lower.min(upper);
+        if stop(lower, upper) {
+            return Ok((lower, upper, witness, true));
+        }
+        if made >= MAX_PAIRS {
+            return Ok((lower, upper, witness, false));
+        }
+        let (ia, x) = &pair.a;
+        let (ib, y) = &pair.b;
+        let (Some(fa), Some(fb)) = (pa.faces.get(*ia), pb.faces.get(*ib)) else {
+            return Err(SolidError::sm_err(OP));
+        };
+        // a witness from the pair's centres (none can be nearer than the pair's bound, so a
+        // pair past the witness so far is skipped)
+        if pair.lower < upper {
+            witness_from(
+                k,
+                (&mut *a, fa, x),
+                (&mut *b, fb, y),
+                &mut upper,
+                &mut witness,
+            )?;
+        }
+        // split the larger piece
+        let size = |p: &Piece| -> Result<f64> {
+            let (l, h) = p.widened_box()?;
+            Ok((0..3)
+                .map(|i| h.get(i).copied().unwrap_or(0.0) - l.get(i).copied().unwrap_or(0.0))
+                .fold(0.0_f64, f64::max))
+        };
+        let split_a = size(x)? >= size(y)?;
+        let (from, keep) = if split_a { (x, y) } else { (y, x) };
+        // the witness's direction, for the hulls' extents along it
+        let axis = witness.map(|[p, q]| [q.x - p.x, q.y - p.y, q.z - p.z]);
+        for q in from.quarters()? {
+            let (fp, brep) = if split_a {
+                (fa, &mut *a)
+            } else {
+                (fb, &mut *b)
+            };
+            if piece_on_face(k, brep, fp, &q)? == Some(false) {
+                continue;
+            }
+            let (qa, qb) = if split_a { (&q, keep) } else { (keep, &q) };
+            let lower = (pair_lower(qa, qb, axis)? - fa.slack - fb.slack).max(0.0);
+            heap.push(Pair {
+                lower,
+                a: (*ia, qa.clone()),
+                b: (*ib, qb.clone()),
+            });
+            made += 1;
+        }
+    }
+}
+
+/// A witness from a pair's centres, each moved to the other's surface's closest point and
+/// back ([`PROJECTIONS`] times), kept where both points are decisively on their faces and
+/// nearer than `upper`.
+fn witness_from(
+    k: &Kernel,
+    (a, fa, x): (&mut Brep, &FacePieces, &Piece),
+    (b, fb, y): (&mut Brep, &FacePieces, &Piece),
+    upper: &mut f64,
+    witness: &mut Option<[Point3; 2]>,
+) -> Result<()> {
+    let (_, ua) = x.centre();
+    let (_, ub) = y.centre();
+    let mut pa_pt = on_face(k, a, fa, Point2::new(ua[0], ua[1]))?;
+    let mut pb_pt = on_face(k, b, fb, Point2::new(ub[0], ub[1]))?;
+    for _ in 0..PROJECTIONS {
+        if let Some(p) = pa_pt {
+            if let Some((q, _)) = closest_on(k, b, fb, p)? {
+                pb_pt = Some(q);
+            }
+        }
+        if let Some(q) = pb_pt {
+            if let Some((p, _)) = closest_on(k, a, fa, q)? {
+                pa_pt = Some(p);
+            }
+        }
+        if let (Some(p), Some(q)) = (pa_pt, pb_pt) {
+            let d = distance_up(p, q) + fa.slack + fb.slack;
+            if d < *upper {
+                *upper = d;
+                *witness = Some([p, q]);
+            }
+        }
+    }
+    Ok(())
+}
+
+/// A point of `inner`'s boundary classified in `outer`: `Some(true)` inside its material,
+/// `Some(false)` outside it, `None` when none of `inner`'s points tried (its pieces'
+/// centres, where their faces hold them: [`on_face`]) is proved more than twice the
+/// classification's tolerance from `outer`'s boundary
+/// ([`crate::query::distance_bound::face_distance_lower_bound`] over `outer`'s faces), or
+/// its classification is not decided.
+fn point_in(
+    k: &Kernel,
+    inner: &mut Brep,
+    p_inner: &Prepared,
+    outer: &Brep,
+) -> Result<Option<bool>> {
+    use crate::query::distance_bound::face_distance_lower_bound;
+    let tol = brep_zone_tol_3d(k, outer, SmallTopoType::NormalSize);
+    let clear = 2.0 * tol;
+    for fp in &p_inner.faces {
+        for piece in &fp.pieces {
+            let (_, uv) = piece.centre();
+            let Some(p) = on_face(k, inner, fp, Point2::new(uv[0], uv[1]))? else {
+                continue;
+            };
+            let mut far = true;
+            for &g in outer.face_list() {
+                if face_distance_lower_bound(k, outer, g, p, clear)? <= clear {
+                    far = false;
+                    break;
+                }
+            }
+            if !far {
+                continue;
+            }
+            let mut classifier = BrepClassifier::of_brep(outer);
+            let pc = match classifier.point_3d_classify(k, p, tol, true) {
+                Ok(pc) => pc,
+                Err(e) if e.is_cancelled() => return Err(e),
+                Err(_) => return Ok(None),
+            };
+            return Ok(match (pc.point_class(), pc.object()) {
+                (PointClassType::Region, Some(EntityRef::Region(r))) => {
+                    Some(region_material(classifier.brep(), r)?)
+                }
+                _ => None,
+            });
+        }
+    }
+    Ok(None)
+}
+
+/// Whether `a` and `b` meet at `tol` (module docs), `a` prepared once for several `b`s.
+///
+/// # Errors
+///
+/// `InvalidInput` for a tolerance not finite and positive, or a body without faces; those
+/// of the pieces, the classifications' cancellation and the closest-point solves' setup.
+pub fn meets_prepared(
+    k: &Kernel,
+    a: &mut Brep,
+    pa: &Prepared,
+    b: &mut Brep,
+    tol: f64,
+) -> Result<Contact> {
+    if !(tol.is_finite() && tol > 0.0) {
+        return Err(SolidError::invalid_input(
+            OP,
+            "the tolerance is finite and positive",
+        ));
+    }
+    let pb = Prepared::new(k, b, pa.origin)?;
+    let stop = |lower: f64, upper: f64| upper <= tol || lower > tol;
+    let (lower, upper, witness, done) = bracket(k, a, pa, b, &pb, tol, &stop)?;
+    let verdict = if upper <= tol {
+        Verdict::Yes
+    } else if !done {
+        Verdict::Undecided
+    } else {
+        // the boundaries are more than tol apart: they meet only if one encloses the other
+        match (point_in(k, b, &pb, a)?, point_in(k, a, pa, b)?) {
+            (Some(true), _) | (_, Some(true)) => Verdict::Yes,
+            (Some(false), Some(false)) => Verdict::No,
+            _ => Verdict::Undecided,
+        }
+    };
+    Ok(Contact {
+        verdict,
+        lower,
+        upper,
+        witness,
+    })
+}
+
+/// [`meets_prepared`] for one pair: `a` prepared about the centre of both bodies' boxes.
+///
+/// # Errors
+///
+/// Those of [`meets_prepared`].
+pub fn meets(k: &Kernel, a: &mut Brep, b: &mut Brep, tol: f64) -> Result<Contact> {
+    let origin = common_origin(k, a, b)?;
+    let pa = Prepared::new(k, a, origin)?;
+    meets_prepared(k, a, &pa, b, tol)
+}
+
+/// The bracket of the distance between `a`'s and `b`'s boundaries, refined until it is at
+/// most `width` wide (relative to the upper end, or absolute below 1) or the pairs reach
+/// [`MAX_PAIRS`]; the verdict is `Yes` when it is that narrow, `Undecided` when the work
+/// ran out first. The bracket always holds the distance; it can stay wider than `width`
+/// where a nearest point lies on a face's trim, since a piece straddling the trim carries
+/// surface off the face and is cut finer only along the trim.
+///
+/// # Errors
+///
+/// `InvalidInput` for a width not finite and positive; those of [`meets`].
+pub fn boundary_distance(k: &Kernel, a: &mut Brep, b: &mut Brep, width: f64) -> Result<Contact> {
+    if !(width.is_finite() && width > 0.0) {
+        return Err(SolidError::invalid_input(
+            OP,
+            "the width is finite and positive",
+        ));
+    }
+    let origin = common_origin(k, a, b)?;
+    let pa = Prepared::new(k, a, origin)?;
+    let pb = Prepared::new(k, b, origin)?;
+    let stop = |lower: f64, upper: f64| upper - lower <= width * upper.max(1.0);
+    let (lower, upper, witness, done) = bracket(k, a, &pa, b, &pb, f64::INFINITY, &stop)?;
+    Ok(Contact {
+        verdict: if done {
+            Verdict::Yes
+        } else {
+            Verdict::Undecided
+        },
+        lower,
+        upper,
+        witness,
+    })
+}

solidcore/crates/solidcore/src/query/distance_bound.rs

@@ -44,15 +44,18 @@ pub const MAX_PATCHES: usize = 4096;
 /// A piece of a face's surface: its homogeneous coordinates about the point, as Bernstein
 /// polynomials on the unit square, and the absolute error each coefficient may carry.
 #[derive(Debug, Clone)]
-struct Piece {
+pub(crate) struct Piece {
     /// `X W - p W`, `Y W - p W`, `Z W - p W`, `W`.
-    h: [Bernstein2; 4],
+    pub(crate) h: [Bernstein2; 4],
     /// The bound on every coefficient's error.
-    err: f64,
+    pub(crate) err: f64,
     /// The certified lower bound of the piece's distance from the point.
-    lower: f64,
+    pub(crate) lower: f64,
     /// Its subdivision depth.
-    depth: usize,
+    pub(crate) depth: usize,
+    /// Its rectangle in the surface's NURB parameters, `[u0, u1, v0, v1]` (to the rounding
+    /// of the splits that made it).
+    pub(crate) uv: [f64; 4],
 }
 
 impl PartialEq for Piece {
@@ -74,21 +77,62 @@ impl Ord for Piece {
 }
 
 impl Piece {
-    /// The piece with its bound, the coefficients `h` within `err` each.
-    fn new(h: [Bernstein2; 4], err: f64, depth: usize) -> Result<Self> {
+    /// The piece with its bound, the coefficients `h` within `err` each, over `uv`.
+    fn new(h: [Bernstein2; 4], err: f64, depth: usize, uv: [f64; 4]) -> Result<Self> {
         let lower = hull_bound(&h, err)?;
         Ok(Self {
             h,
             err,
             lower,
             depth,
+            uv,
         })
     }
 
+    /// The box of its Euclidean control points about the point, each coordinate widened
+    /// by the error a quotient of coefficients within `err` can carry: every point of the
+    /// piece lies in it (the hull, with positive weights).
+    pub(crate) fn widened_box(&self) -> Result<([f64; 3], [f64; 3])> {
+        let (lo, hi, widen) = euclid_box(&self.h, self.err)?;
+        Ok((lo.map(|l| l - widen), hi.map(|u| u + widen)))
+    }
+
+    /// Its centre's Euclidean point about the point, and the centre's parameters.
+    pub(crate) fn centre(&self) -> ([f64; 3], [f64; 2]) {
+        let [x, y, z, w] = &self.h;
+        let wc = w.eval(0.5, 0.5);
+        let [u0, u1, v0, v1] = self.uv;
+        (
+            [
+                x.eval(0.5, 0.5) / wc,
+                y.eval(0.5, 0.5) / wc,
+                z.eval(0.5, 0.5) / wc,
+            ],
+            [0.5 * (u0 + u1), 0.5 * (v0 + v1)],
+        )
+    }
+
+    /// The certified range of `x . u` over the piece, for a direction `u`: the least and
+    /// largest of its Euclidean control points' (its hull), each moved out by the widening
+    /// along `u` and the products' rounding.
+    pub(crate) fn extent_along(&self, u: [f64; 3]) -> Result<(f64, f64)> {
+        let (lo, hi, widen) = euclid_box(&self.h, self.err)?;
+        let size = lo.iter().chain(&hi).fold(0.0_f64, |m, v| m.max(v.abs())) + widen;
+        let norm1 = u[0].abs() + u[1].abs() + u[2].abs();
+        let (mut least, mut most) = (f64::INFINITY, f64::NEG_INFINITY);
+        for e in euclid_points(&self.h) {
+            let d = e[0] * u[0] + e[1] * u[1] + e[2] * u[2];
+            least = least.min(d);
+            most = most.max(d);
+        }
+        let slack = widen * norm1 + gamma(3) * size * norm1;
+        Ok((least - slack, most + slack))
+    }
+
     /// Its four quarters: one de Casteljau level each way, which adds the rounding of a
     /// convex combination per level (`gamma` of twice the degree plus two) of the largest
     /// coefficient.
-    fn quarters(&self) -> Result<[Self; 4]> {
+    pub(crate) fn quarters(&self) -> Result<[Self; 4]> {
         let largest = self
             .h
             .iter()
@@ -108,14 +152,23 @@ impl Piece {
             split(&self.h[2]),
             split(&self.h[3]),
         ];
+        let [u0, u1, v0, v1] = self.uv;
+        let (um, vm) = (0.5 * (u0 + u1), 0.5 * (v0 + v1));
+        let rects = [
+            [u0, um, v0, vm],
+            [u0, um, vm, v1],
+            [um, u1, v0, vm],
+            [um, u1, vm, v1],
+        ];
         let mut out = Vec::new();
-        for i in 0..4 {
+        for (i, uv) in rects.into_iter().enumerate() {
             let pick =
                 |a: &[Bernstein2; 4]| a.get(i).cloned().ok_or_else(|| SolidError::sm_err(OP));
             out.push(Self::new(
                 [pick(&x)?, pick(&y)?, pick(&z)?, pick(&w)?],
                 err,
                 self.depth + 1,
+                uv,
             )?);
         }
         out.try_into().map_err(|_| SolidError::sm_err(OP))
@@ -130,42 +183,15 @@ impl Piece {
 /// least the least P . u, and |X| >= X . u / |u|), which, unlike the box, tightens with the
 /// square of the piece's size on a piece tilted to the axes.
 fn hull_bound(h: &[Bernstein2; 4], err: f64) -> Result<f64> {
-    let [x, y, z, w] = h;
-    let (xs, ys, zs, ws) = (
-        x.coefficients(),
-        y.coefficients(),
-        z.coefficients(),
-        w.coefficients(),
-    );
-    let w_min = ws.iter().copied().fold(f64::INFINITY, f64::min);
-    if !above(w_min - err, 0.0) {
-        return Err(SolidError::Unsupported {
-            operation: OP,
-            reason: "a patch's weights are not proved positive",
-        });
-    }
-    let mut points = Vec::new();
-    let mut lo = [f64::INFINITY; 3];
-    let mut hi = [f64::NEG_INFINITY; 3];
-    let mut size = 0.0_f64;
+    let (lo, hi, widen) = euclid_box(h, err)?;
+    let points = euclid_points(h);
+    let size = lo.iter().chain(&hi).fold(0.0_f64, |m, v| m.max(v.abs()));
     let mut sum_points = [0.0_f64; 3];
-    for (((&a, &b), &c), &wi) in xs.iter().zip(ys).zip(zs).zip(ws) {
-        let e = [a / wi, b / wi, c / wi];
-        for (k, v) in e.into_iter().enumerate() {
-            size = size.max(v.abs());
-            if let (Some(l), Some(u), Some(t)) =
-                (lo.get_mut(k), hi.get_mut(k), sum_points.get_mut(k))
-            {
-                *l = l.min(v);
-                *u = u.max(v);
-                *t += v;
-            }
+    for e in &points {
+        for (t, v) in sum_points.iter_mut().zip(e) {
+            *t += v;
         }
-        points.push(e);
     }
-    // |(v + dv) / (w + dw) - v / w| <= err (1 + |v / w|) / (w - err), and the quotient's
-    // own rounding
-    let widen = err * (1.0 + size) / (w_min - err) + gamma(2) * size;
     let mut sum = 0.0_f64;
     for (l, u) in lo.into_iter().zip(hi) {
         let out = (l - widen).max(-(u + widen)).max(0.0);
@@ -202,6 +228,61 @@ fn hull_bound(h: &[Bernstein2; 4], err: f64) -> Result<f64> {
     Ok(bound.max(0.0))
 }
 
+/// The Euclidean control points of the homogeneous coefficients `h`.
+fn euclid_points(h: &[Bernstein2; 4]) -> Vec<[f64; 3]> {
+    let [x, y, z, w] = h;
+    x.coefficients()
+        .iter()
+        .zip(y.coefficients())
+        .zip(z.coefficients())
+        .zip(w.coefficients())
+        .map(|(((&a, &b), &c), &wi)| [a / wi, b / wi, c / wi])
+        .collect()
+}
+
+/// The box of the Euclidean control points of `h`, and how far each coordinate is to be
+/// widened for coefficients within `err` (and the quotients' own rounding): `(lo, hi,
+/// widen)`.
+fn euclid_box(h: &[Bernstein2; 4], err: f64) -> Result<([f64; 3], [f64; 3], f64)> {
+    let [x, y, z, w] = h;
+    let (xs, ys, zs, ws) = (
+        x.coefficients(),
+        y.coefficients(),
+        z.coefficients(),
+        w.coefficients(),
+    );
+    let w_min = ws.iter().copied().fold(f64::INFINITY, f64::min);
+    if !above(w_min - err, 0.0) {
+        return Err(SolidError::Unsupported {
+            operation: OP,
+            reason: "a patch's weights are not proved positive",
+        });
+    }
+    let mut lo = [f64::INFINITY; 3];
+    let mut hi = [f64::NEG_INFINITY; 3];
+    let mut size = 0.0_f64;
+    for (((&a, &b), &c), &wi) in xs.iter().zip(ys).zip(zs).zip(ws) {
+        let e = [a / wi, b / wi, c / wi];
+        for (k, v) in e.into_iter().enumerate() {
+            size = size.max(v.abs());
+            if let (Some(l), Some(u)) = (lo.get_mut(k), hi.get_mut(k)) {
+                *l = l.min(v);
+                *u = u.max(v);
+            }
+        }
+    }
+    // |(v + dv) / (w + dw) - v / w| <= err (1 + |v / w|) / (w - err), and the quotient's
+    // own rounding
+    let widen = err * (1.0 + size) / (w_min - err) + gamma(2) * size;
+    if !(widen.is_finite() && lo.iter().chain(&hi).all(|v| v.is_finite())) {
+        return Err(SolidError::Unsupported {
+            operation: OP,
+            reason: "a patch's box is not finite",
+        });
+    }
+    Ok((lo, hi, widen))
+}
+
 /// A certified lower bound of the distance from `p` to face `f` of `b`: no point of the face
 /// is nearer to `p`. The face is bounded by its surface over the face's UV box (its trims
 /// lie inside it): [`surface_distance_lower_bound`] on them.
@@ -320,6 +401,51 @@ pub fn surface_distance_lower_bound(
             "the bound asked for is finite and not negative",
         ));
     }
+    let (pieces, slack) = pieces_about(k, surf, uv_box, p)?;
+    let want_net = want + slack;
+    let mut made = pieces.len();
+    let mut heap: BinaryHeap<Piece> = pieces.into_iter().collect();
+    let off = |lower: f64| (lower - slack).max(0.0);
+    loop {
+        // (the box is inside the patches' range, so some patch overlaps it)
+        let nearest = heap.pop().ok_or_else(|| SolidError::sm_err(OP))?;
+        if nearest.lower >= want_net || made >= MAX_PATCHES {
+            return Ok(off(nearest.lower));
+        }
+        // the piece's centre is a point of the net over the box: nearer than the target, no
+        // bound reaches it
+        let (centre, _) = nearest.centre();
+        if (centre[0] * centre[0] + centre[1] * centre[1] + centre[2] * centre[2]).sqrt() < want_net
+        {
+            return Ok(off(nearest.lower));
+        }
+        for q in nearest.quarters()? {
+            heap.push(q);
+        }
+        made += 4;
+    }
+}
+
+/// The pieces of `surf` over `uv_box` (as [`surface_distance_lower_bound`] takes it) about
+/// `p`: its NURBS form's Bézier patches restricted to the box, translated by `-p`, each
+/// with its coefficients' error (the form's against its definition, the extraction's and
+/// the translation's, module docs); and the slack every distance from them is to lose (a
+/// net's certified bound from its definition, 0 for a surface that evaluates on its net).
+/// The one construction of the face pieces that [`surface_distance_lower_bound`] and the
+/// contact test between bodies ([`crate::query::contact`]) bound distances with.
+///
+/// # Errors
+///
+/// As [`surface_distance_lower_bound`]'s, but for `want`.
+pub(crate) fn pieces_about(
+    k: &Kernel,
+    surf: &Surface,
+    uv_box: Option<&Extent2d>,
+    p: Point3,
+) -> Result<(Vec<Piece>, f64)> {
+    if !(p.x.is_finite() && p.y.is_finite() && p.z.is_finite()) {
+        return Err(SolidError::invalid_input(OP, "point not finite"));
+    }
     if !measurable(k, surf)? {
         return Err(SolidError::Unsupported {
             operation: OP,
@@ -424,9 +550,7 @@ pub fn surface_distance_lower_bound(
             ));
         }
     }
-    let want_net = want + slack;
-    let mut heap = BinaryHeap::new();
-    let mut made = 0_usize;
+    let mut pieces = Vec::new();
     for patch in patches.iter().flatten() {
         let (pu, pv) = patch.degrees();
         let (uk, vk) = (patch.u_knots(), patch.v_knots());
@@ -496,33 +620,15 @@ pub fn surface_distance_lower_bound(
         } else {
             h
         };
-        heap.push(Piece::new(h, err, 0)?);
-        made += 1;
-    }
-    let off = |lower: f64| (lower - slack).max(0.0);
-    loop {
-        // (the box is inside the patches' range, so some patch overlaps it)
-        let nearest = heap.pop().ok_or_else(|| SolidError::sm_err(OP))?;
-        if nearest.lower >= want_net || made >= MAX_PATCHES {
-            return Ok(off(nearest.lower));
-        }
-        // the piece's centre is a point of the net over the box: nearer than the target, no
-        // bound reaches it
-        let [x, y, z, w] = &nearest.h;
-        let wc = w.eval(0.5, 0.5);
-        let centre = (
-            x.eval(0.5, 0.5) / wc,
-            y.eval(0.5, 0.5) / wc,
-            z.eval(0.5, 0.5) / wc,
-        );
-        if (centre.0 * centre.0 + centre.1 * centre.1 + centre.2 * centre.2).sqrt() < want_net {
-            return Ok(off(nearest.lower));
-        }
-        for q in nearest.quarters()? {
-            heap.push(q);
-        }
-        made += 4;
+        let uv = [
+            u0 + a0 * (u1 - u0),
+            u0 + a1 * (u1 - u0),
+            v0 + b0 * (v1 - v0),
+            v0 + b1 * (v1 - v0),
+        ];
+        pieces.push(Piece::new(h, err, 0, uv)?);
     }
+    Ok((pieces, slack))
 }
 
 #[cfg(test)]
@@ -890,7 +996,10 @@ mod tests {
                     .unwrap()
                 });
                 let exact = h.each_ref().map(exact_quarters);
-                let quarters = Piece::new(h, 0.0, 0).unwrap().quarters().unwrap();
+                let quarters = Piece::new(h, 0.0, 0, [0.0, 1.0, 0.0, 1.0])
+                    .unwrap()
+                    .quarters()
+                    .unwrap();
                 for (q, quarter) in quarters.iter().enumerate() {
                     for (axis, component) in quarter.h.iter().enumerate() {
                         for (rounded, reference) in

solidcore/crates/solidcore/src/surface/implicit.rs

@@ -11,7 +11,6 @@
 use crate::error::{Result, SolidError};
 use crate::geom::axis2_placement::Axis2Placement;
 use crate::geom::extent2d::Extent2d;
-use crate::geom::vector2d::Point2;
 use crate::kernel::Kernel;
 use crate::math::{Vec2, Vec3, PI};
 use crate::nurbs::b_spline_curve::BSplineCurve;
@@ -69,42 +68,6 @@ pub enum Implicit {
     },
 }
 
-impl Implicit {
-    /// How far `x` is from this surface.
-    #[must_use]
-    pub fn distance(&self, x: Vec3) -> f64 {
-        // the axial and radial parts of `x` about an axis through `o`
-        let about = |o: Vec3, axis: Vec3| {
-            let v = x - o;
-            let h = v.dot(axis);
-            (h, (v - axis * h).length())
-        };
-        match *self {
-            Self::Plane { origin, normal } => (x - origin).dot(normal).abs(),
-            Self::Cone {
-                origin,
-                axis,
-                radius,
-                tan_angle,
-            } => {
-                let (h, rho) = about(origin, axis);
-                (rho - (radius + tan_angle * h)).abs() / (1.0 + tan_angle * tan_angle).sqrt()
-            }
-            Self::Sphere { centre, radius } => ((x - centre).length() - radius).abs(),
-            Self::Torus {
-                centre,
-                axis,
-                major,
-                minor,
-            } => {
-                let (h, rho) = about(centre, axis);
-                let (a, b) = (rho - major, h);
-                ((a * a + b * b).sqrt() - minor).abs()
-            }
-        }
-    }
-}
-
 /// A polynomial in Bernstein form over a patch, degrees `(p, q)`, `c[i (q + 1) + j]`.
 #[derive(Debug, Clone, PartialEq)]
 struct BPoly {
@@ -343,7 +306,7 @@ fn bezier_on<'a>(
     f.c.iter().all(within)
 }
 
-pub(crate) const OP_EXACT: &str = "surface::implicit::check_exact_net";
+const OP_EXACT: &str = "surface::implicit::check_exact_net";
 
 /// A quadric `x.A.x + 2 b.x + c = 0` as `(A, b, c)`.
 pub type QuadricEquation = ([[f64; 3]; 3], [f64; 3], f64);
@@ -439,46 +402,6 @@ pub fn check_exact_net(k: &Kernel, net: &Surface, check: &ExactNetCheck) -> Resu
     Ok(())
 }
 
-/// The number of intervals [`net_miss`] samples a net at, each way.
-const MISS_SAMPLES: u32 = 16;
-
-/// How far `net` strays from the class's `surface`, sampled over its domain. A guard, not a
-/// certificate: [`check_exact_net`] is the proof, and a net that passes only this keeps no
-/// definition. It tells a net that the class's rounding leaves off the surface from one
-/// built wrongly.
-///
-/// # Errors
-///
-/// Those of the evaluations.
-pub fn net_miss(k: &Kernel, net: &Surface, surface: &ClassSurface) -> Result<f64> {
-    let dom = net.natural_uv_domain_with_nurb(k)?;
-    let (lo, hi) = (dom.min(), dom.max());
-    let mut worst = 0.0_f64;
-    for i in 0..=MISS_SAMPLES {
-        for j in 0..=MISS_SAMPLES {
-            let (fu, fv) = (
-                f64::from(i) / f64::from(MISS_SAMPLES),
-                f64::from(j) / f64::from(MISS_SAMPLES),
-            );
-            let uv = Point2::new(lo.x + (hi.x - lo.x) * fu, lo.y + (hi.y - lo.y) * fv);
-            let x = net.evaluate_point(k, uv)?;
-            let miss = match surface {
-                ClassSurface::Implicit(a) => a.distance(x),
-                ClassSurface::Canal {
-                    spine, map, radius, ..
-                } => {
-                    let w = if map.along_v { uv.y } else { uv.x };
-                    let ivl = spine.natural_interval();
-                    let s = (map.a + map.b * w).max(ivl.min()).min(ivl.max());
-                    ((x - spine.evaluate_point(s)?).length() - *radius).abs()
-                }
-            };
-            worst = worst.max(miss);
-        }
-    }
-    Ok(worst)
-}
-
 /// A Bezier segment of a rational curve, Bernstein polynomials of degree `(p, 0)`: the
 /// homogeneous `(Y - O V)` components and the weight `V`.
 type CurveSegment = ([BPoly; 3], BPoly);
@@ -1200,12 +1123,37 @@ impl AnalyticSupport {
         // a surface of revolution: the one the kernel recognises in the form, which must be
         // this one, take the form as its NURB and evaluate as the form at the form's own
         // parameters
-        let not_form = || {
-            SolidError::invalid_input(
-                OP_BUILD,
-                "a support's analytic surface is not its NURBS form's parametrisation",
-            )
-        };
+        let surface = self.on_form(k, &form)?;
+        if !evaluates_as(k, &surface, &form)? {
+            return Err(not_form());
+        }
+        Ok(surface)
+    }
+
+    /// The largest ratio, over the form's Greville points, of this support's analytic
+    /// surface's distance from `form` to the rounding [`evaluates_as`] derives there (#3185):
+    /// at most 1 where [`Self::build`]'s parametrisation test passes. The surface is the one
+    /// `build` makes (the kernel's recognition of `form`, carrying this support's numbers);
+    /// `lies_on` and the placement's axes are not checked here. `None` for a plane, or a
+    /// surface the bound does not cover.
+    ///
+    /// # Errors
+    ///
+    /// `InvalidInput` where the kernel recognises no surface of this kind in `form`; those of
+    /// the recognisers, the conversions and the evaluations.
+    pub fn parametrisation_margin(&self, k: &Kernel, form: &BSplineSurface) -> Result<Option<f64>> {
+        if matches!(self, Self::Plane { .. }) {
+            return Ok(None);
+        }
+        let surface = self.on_form(k, form)?;
+        evaluation_margin(k, &surface, form)
+    }
+
+    /// The surface of revolution of this support's kind the kernel recognises in `form`,
+    /// carrying this support's own numbers and `form` as its NURB ([`Self::build`] after its
+    /// checks of the point set, before the parametrisation's).
+    fn on_form(&self, k: &Kernel, form: &BSplineSurface) -> Result<Surface> {
+        let form = form.copy();
         let mut surface = recognise(k, self.kind(), &form)?.ok_or_else(not_form)?;
         // the recognised surface must be of this support's kind; the numbers recognition
         // recomputed are not compared with the support's, because it recomputes them from
@@ -1255,13 +1203,18 @@ impl AnalyticSupport {
             }
             _ => return Err(not_form()),
         }
-        if !evaluates_as(k, &surface, &form)? {
-            return Err(not_form());
-        }
         Ok(surface)
     }
 }
 
+/// [`AnalyticSupport::build`]'s refusal of a surface that is not its form's parametrisation.
+fn not_form() -> SolidError {
+    SolidError::invalid_input(
+        OP_BUILD,
+        "a support's analytic surface is not its NURBS form's parametrisation",
+    )
+}
+
 /// A cylinder's or cone's recognised parametrisation carried by the support's own
 /// `placement`: its heights are measured along the axis from the placement's origin, so the
 /// recognised heights move by how far along the axis the recognised origin lies from it.
@@ -1336,41 +1289,314 @@ fn recognise(k: &Kernel, kind: SupportKind, form: &BSplineSurface) -> Result<Opt
 }
 
 /// Whether `surface` evaluates as `form` at the form's Greville parameters (each control
-/// point's), within the rounding of the two evaluations: `gamma_K` times the form's size,
-/// `K = 16 (p + q + 2)` for the rational evaluation of degrees `(p, q)` and the analytic
-/// one's conversion of the parameters and its trigonometry, a few ulps per operation.
+/// point's): at each, the analytic surface's point and the form's agree within the bound
+/// [`Rounding::at`] derives for that point, term by term below. The analytic side is the
+/// **STEP** evaluation of the converted parameters ([`Surface::convert_uv_from_nurbs_to_step`]
+/// and [`Surface::evaluate_step_point`]), which a cone, sphere, torus and surface of
+/// revolution compute from their placement, radii, semi-angle and generator.
+/// [`Surface::evaluate_point`] would not do: a sphere, a torus and a surface of revolution
+/// take it from their NURBS base, and that base is the `form` the caller has just installed,
+/// so it would compare the form with itself and read nothing of what the caller imposed
+/// (#2588).
 ///
-/// The analytic side is the **STEP** evaluation of the converted parameters
-/// ([`Surface::convert_uv_from_nurbs_to_step`] and [`Surface::evaluate_step_point`]),
-/// which a cone, sphere, torus and surface of revolution compute from their placement,
-/// radii, semi-angle and generator. [`Surface::evaluate_point`] would not do: a sphere, a
-/// torus and a surface of revolution take it from their NURBS base, and that base is the
-/// `form` the caller has just installed, so it would compare the form with itself and read
-/// nothing of what the caller imposed (#2588; the cone already converted and evaluated
-/// analytically, which is why only the other kinds went unchecked).
+/// What the bound forgives and what it tests (#3185): it forgives the rounding of both
+/// evaluations and of the two parameter conversions, nothing else. It does not forgive the
+/// label's frame turned about the axis against the converter's, its axis, origin or radii,
+/// a form whose parameter is not the analytic surface's (a reparametrised arc), nor the
+/// form's distance from the analytic surface: those are what the comparison proves.
+///
+/// `u` is the unit roundoff, `gamma_n = n u / (1 - n u)`. Stated assumption: the libm
+/// `sin`, `cos`, `tan` and `atan2` are within one ulp (relative `2 u`).
+///
+/// 1. **The form** (`Surface::eval`, Cox-de Boor then two blends): each basis value is a
+///    nonnegative combination over `p` levels of at most 6 operations (two knot differences,
+///    a sum, a division, a product and a sum), so relative `gamma_6p`; the blends sum
+///    `q + 1` then `p + 1` products. With positive weights the homogeneous point's error is
+///    at most `gamma_(7 (p + q) + 4) W size` and the weight's `gamma_(7 (p + q) + 4) W`, so
+///    the Euclidean point's, after the division, `E_form = (2 gamma_(7 (p + q) + 4) + u)
+///    size` (`gamma_(7 (p + q) + 4) size` when polynomial), `size` the net's largest point.
+/// 2. **The STEP point**, `O + m (cos u X + sin u Y) + h Z`, `m` and `h` the generator's
+///    radial and axial parts (the torus `R + r cos v`, `r sin v`; the sphere `r cos v`,
+///    `r sin v`; the cone `r0 + v tan a`, `v`). An angle `x` in degrees goes to radians in
+///    two operations of the rounded `PI`: `gamma_3 |x|`; with one ulp of the cosine and sine,
+///    each is off by `dc(x) = 2 u + gamma_3 |x|`. Per coordinate (unit axes) the error is at
+///    most `gamma_3 (|O| + M + H) + M (gamma_2 + 2 dc(u)) + dm + dh`, `M`, `H` the sizes of
+///    `m`, `h` and `dm`, `dh` their errors: the torus `gamma_3 (R + r) + 2 r dc(v)`, the
+///    sphere `gamma_1 r + 2 r dc(v)`, the cone `gamma_2 (|r0| + |v t|) + |v| dt` (`dt =
+///    (1 + t^2) gamma_3 |a| + 2 u |t|` for the tangent of the half angle `a`); a general
+///    surface of revolution rotates its generator's point `G` (a B-spline: item 1's bound of
+///    its own net, `E_gen`) about the axis, `gamma_8 (|O| + 3 D) + 2 D dc(u) + E_gen`, `D`
+///    the net's reach from `O`. `E_step` is `sqrt 3` times the coordinate bound.
+/// 3. **The generator's parameter**: a circle (torus, sphere) inverts its NURBS point by
+///    `atan2` about its centre `c`: at its radius `r` the error is that point's (`E_form`)
+///    plus `gamma_4 (size + |c|)` and `gamma_4 r |v|` of the angle's own rounding; a line
+///    (cone) maps it affinely in 6 operations, `gamma_6 2 (|v0| + |v1|)` along the axis and
+///    its slope, times `1 + |t|`; a B-spline generator's parameter is its NURBS's (none).
+/// 4. **The sweep's parameter**: the polar converter evaluates its circle's NURBS at the
+///    parameter (item 1's bound of that NURBS's own net, `E_polar`) and inverts the point by
+///    `atan2` about the circle's centre `c` in its frame: at the circle's radius `rc` the
+///    angle is off by `(E_polar + gamma_4 (rc + 2 |c|)) / rc`, plus `gamma_4 |s|` of the
+///    angle's own rounding and the degrees' period (`gamma_1 2 pi`), times the STEP point's
+///    distance from that axis (the lever). The converter's ends are its domain's, which
+///    recognition takes from the form's own end meridians (#3199), so they are the inverted
+///    ends to the same rounding.
+///
+/// A surface of another kind, or a general surface of revolution whose generator is not a
+/// B-spline, has no bound here and is not proved.
 ///
 /// # Errors
 ///
 /// Those of the parameter conversion and the evaluations and of the form's knots.
 fn evaluates_as(k: &Kernel, surface: &Surface, form: &BSplineSurface) -> Result<bool> {
+    Ok(evaluation_margin(k, surface, form)?.is_some_and(|m| m <= 1.0))
+}
+
+/// [`evaluates_as`]'s largest ratio of distance to bound over the Greville points (infinite
+/// for a distance that is not a number); `None` for a surface the bound does not cover.
+fn evaluation_margin(k: &Kernel, surface: &Surface, form: &BSplineSurface) -> Result<Option<f64>> {
     let (net, p, q) = euclidean_net(form)?;
     let size = net.iter().flatten().fold(0.0_f64, |a, x| a.max(x.length()));
-    let tol = gamma(16 * (p + q + 2)) * size;
+    let Some(r) = Rounding::of(surface, form.is_rational(), p, q, size, &net)? else {
+        return Ok(None);
+    };
+    let mut worst = 0.0_f64;
     for u in greville(form, true)? {
         for v in greville(form, false)? {
             let uv = Vec2::new(u, v);
             let step_uv = surface.convert_uv_from_nurbs_to_step(k, uv)?;
-            let off =
-                (surface.evaluate_step_point(k, step_uv)? - form.evaluate_point(uv)?).length();
-            if !matches!(
-                off.partial_cmp(&tol),
-                Some(core::cmp::Ordering::Less | core::cmp::Ordering::Equal)
-            ) {
-                return Ok(false);
-            }
+            let s = surface.evaluate_step_point(k, step_uv)?;
+            let f = form.evaluate_point(uv)?;
+            let ratio = (s - f).length() / r.at(step_uv, s, f);
+            worst = if ratio.is_nan() {
+                f64::INFINITY
+            } else {
+                worst.max(ratio)
+            };
         }
     }
-    Ok(true)
+    Ok(Some(worst))
+}
+
+/// A surface of revolution's generator, for [`Rounding`]'s item 2 and 3.
+#[derive(Debug, Clone, Copy)]
+enum Generator {
+    /// `R + r cos v`, `r sin v` about the circle of centre `c`.
+    Torus { major: f64, minor: f64, c: f64 },
+    /// `r cos v`, `r sin v` about the circle of centre `c`.
+    Sphere { radius: f64, c: f64 },
+    /// `r0 + v t`, `v` (the half angle `a` in radians), its STEP `v` over `[v0, v1]`.
+    Cone {
+        r0: f64,
+        t: f64,
+        a: f64,
+        v0: f64,
+        v1: f64,
+    },
+    /// A B-spline generator, its own evaluation's bound.
+    Spline { e_gen: f64 },
+}
+
+/// The pieces of [`evaluates_as`]'s bound that do not depend on the point.
+#[derive(Debug, Clone, Copy)]
+struct Rounding {
+    e_form: f64,
+    origin: Vec3,
+    /// The net's reach from the origin.
+    reach: f64,
+    generator: Generator,
+    /// The sweep converter's circle: centre, axes.
+    polar: (Vec3, Vec3, Vec3, f64, f64),
+}
+
+/// `|v|`'s largest coordinate.
+fn max_coord(v: Vec3) -> f64 {
+    v.x.abs().max(v.y.abs()).max(v.z.abs())
+}
+
+/// A rational (or polynomial) evaluation's bound of degrees summing to `n` over points of
+/// size `size` ([`evaluates_as`] item 1).
+fn eval_bound(rational: bool, n: usize, size: f64) -> f64 {
+    let g = gamma(7 * n + 4);
+    if rational {
+        (2.0 * g + crate::nurbs::math::UNIT) * size
+    } else {
+        g * size
+    }
+}
+
+impl Rounding {
+    /// The bound's data for `surface` (a cone, sphere, torus or surface of revolution with a
+    /// B-spline generator) on its form of degrees `(p, q)`, net `net`; `None` for any other.
+    fn of(
+        surface: &Surface,
+        rational: bool,
+        p: usize,
+        q: usize,
+        size: f64,
+        net: &[Vec<Vec3>],
+    ) -> Result<Option<Self>> {
+        use crate::curve::curve::CurveKind;
+        let e_form = eval_bound(rational, p + q, size);
+        let circle_centre = |rev: &crate::surface::surf_of_revolution::SurfOfRevolution| match rev
+            .gen_curve
+            .as_ref()
+            .map(crate::curve::curve::Curve::kind)
+        {
+            Some(CurveKind::Ellipse(e)) => Some(e.position.origin().length()),
+            _ => None,
+        };
+        let (rev, generator) = match surface.kind() {
+            SurfaceKind::Torus(t) => {
+                let Some(c) = circle_centre(&t.base) else {
+                    return Ok(None);
+                };
+                (
+                    &t.base,
+                    Generator::Torus {
+                        major: t.major_radius,
+                        minor: t.minor_radius,
+                        c,
+                    },
+                )
+            }
+            SurfaceKind::Sphere(sp) => {
+                let Some(c) = circle_centre(&sp.base) else {
+                    return Ok(None);
+                };
+                (
+                    &sp.base,
+                    Generator::Sphere {
+                        radius: sp.radius,
+                        c,
+                    },
+                )
+            }
+            SurfaceKind::Cone(c) => {
+                let a = c.semi_angle_deg * PI / 180.0;
+                let d = c.rev.anal_uv_domain;
+                (
+                    &c.rev,
+                    Generator::Cone {
+                        r0: c.base_radius,
+                        t: c.tan_semi_angle(),
+                        a,
+                        v0: d.v_min(),
+                        v1: d.v_max(),
+                    },
+                )
+            }
+            SurfaceKind::SurfOfRevolution(rev) => {
+                let Some(CurveKind::BSpline(g)) =
+                    rev.gen_curve.as_ref().map(crate::curve::curve::Curve::kind)
+                else {
+                    return Ok(None);
+                };
+                let pts = g.control_points()?;
+                let gsize = pts.iter().fold(0.0_f64, |m, c| {
+                    let w = if c.w == NOW { 1.0 } else { c.w };
+                    m.max((Vec3::new(c.x, c.y, c.z) * (1.0 / w)).length())
+                });
+                let rational = pts.iter().any(|c| c.w != NOW);
+                (
+                    &**rev,
+                    Generator::Spline {
+                        e_gen: eval_bound(rational, g.degree()?, gsize),
+                    },
+                )
+            }
+            _ => return Ok(None),
+        };
+        let Some(pc) = rev.polar_curve() else {
+            return Ok(None);
+        };
+        let e_polar = {
+            let pts = pc.base().control_points()?;
+            let psize = pts.iter().fold(0.0_f64, |m, c| {
+                let w = if c.w == NOW { 1.0 } else { c.w };
+                m.max((Vec3::new(c.x, c.y, c.z) * (1.0 / w)).length())
+            });
+            eval_bound(pts.iter().any(|c| c.w != NOW), pc.base().degree()?, psize)
+        };
+        let origin = rev.position.origin();
+        let reach = net
+            .iter()
+            .flatten()
+            .fold(0.0_f64, |m, x| m.max((*x - origin).length()));
+        Ok(Some(Self {
+            e_form,
+            origin,
+            reach,
+            generator,
+            polar: (
+                pc.position.origin(),
+                pc.position.x_axis(),
+                pc.position.y_axis(),
+                pc.radius_x,
+                e_polar,
+            ),
+        }))
+    }
+
+    /// The bound at the STEP parameters `step` (degrees, and the generator's), the STEP
+    /// point `s` and the form's point `f` ([`evaluates_as`] items 1 to 4).
+    fn at(&self, step: Vec2, s: Vec3, f: Vec3) -> f64 {
+        let u = crate::nurbs::math::UNIT;
+        let rad = |deg: f64| deg * PI / 180.0;
+        let dc = |deg: f64| 2.0 * u + gamma(3) * rad(deg).abs();
+        let o = max_coord(self.origin);
+        let size_f = f.length();
+        // item 2 (per coordinate) and item 3 (as a length)
+        let (coord, generator) = match self.generator {
+            Generator::Torus { major, minor, c } => {
+                let (m, h) = (major + minor, minor);
+                let dm = gamma(3) * (major + minor) + 2.0 * minor * dc(step.y);
+                let dh = gamma(1) * minor + 2.0 * minor * dc(step.y);
+                (
+                    gamma(3) * (o + m + h) + m * (gamma(2) + 2.0 * dc(step.x)) + dm + dh,
+                    self.e_form + gamma(4) * (size_f + c) + gamma(4) * minor * rad(step.y).abs(),
+                )
+            }
+            Generator::Sphere { radius, c } => {
+                let dm = gamma(1) * radius + 2.0 * radius * dc(step.y);
+                (
+                    gamma(3) * (o + 2.0 * radius)
+                        + radius * (gamma(2) + 2.0 * dc(step.x))
+                        + 2.0 * dm,
+                    self.e_form + gamma(4) * (size_f + c) + gamma(4) * radius * rad(step.y).abs(),
+                )
+            }
+            Generator::Cone { r0, t, a, v0, v1 } => {
+                let v = step.y;
+                let m = r0.abs() + (v * t).abs();
+                let dt = (1.0 + t * t) * gamma(3) * a.abs() + 2.0 * u * t.abs();
+                let dm = gamma(2) * m + v.abs() * dt;
+                (
+                    gamma(3) * (o + m + v.abs()) + m * (gamma(2) + 2.0 * dc(step.x)) + dm,
+                    gamma(6) * 2.0 * (v0.abs() + v1.abs()) * (1.0 + t.abs()),
+                )
+            }
+            Generator::Spline { e_gen } => (
+                gamma(8) * (o + 3.0 * self.reach) + 2.0 * self.reach * dc(step.x) + e_gen,
+                0.0,
+            ),
+        };
+        let e_step = 3.0_f64.sqrt() * coord;
+        // item 4: the polar converter's inversion of its own circle's point
+        let (c, x, y, rc, e_polar) = self.polar;
+        let z = x.cross(y);
+        let along = |v: Vec3| v - z * v.dot(z);
+        let lever = along(s - c).length();
+        let sweep = if rc > 0.0 {
+            lever
+                * ((e_polar + gamma(4) * (rc + 2.0 * c.length())) / rc
+                    + gamma(4) * rad(step.x).abs()
+                    + gamma(1) * 2.0 * PI)
+        } else {
+            // a circle of no radius has no angle: any, at the STEP point's lever
+            lever * PI
+        };
+        self.e_form + e_step + generator + sweep
+    }
 }
 
 const OP_DERIVE: &str = "surface::implicit::AnalyticSupport::derive";

solidcore/crates/solidcore/src/surface/procedural_blend/exact.rs

@@ -39,11 +39,62 @@ use crate::geom::vector3d::{Point3, Vec3};
 use crate::kernel::Kernel;
 use crate::math::EFF_ZERO;
 use crate::nurbs::b_spline_curve::BSplineCurve;
-use crate::surface::implicit::{check_exact_net, net_miss, ClassSurface, ExactNetCheck, LeaderMap};
+use crate::surface::implicit::{
+    check_exact_net, lies_on, AnalyticSupport, ClassSurface, ExactNetCheck, Implicit, LeaderMap,
+};
 use crate::surface::procedural_blend::chamfer::ChamferDefinition;
 use crate::surface::surface::Surface;
 
 const OP_TRANSFORM: &str = "surface::procedural_blend::exact::ExactRound::transform";
+const OP_CIRCLE: &str = "surface::procedural_blend::exact::ExactRound::supports_distance";
+
+/// The torus a round of the circle class is: a surface of revolution's, about its
+/// placement's axis, of the circle its generating arc is on (the centre of a rational quadratic
+/// arc from its first span's control points: the midpoint of the chord pushed away from the
+/// middle point by `a^2 / h`, `a` the half chord, `h` the middle point's distance from the
+/// chord's midpoint); or a torus label's own. `None` for any other surface.
+///
+/// # Errors
+///
+/// Those of reading the arc's control polygon.
+fn own_torus(net: &Surface) -> Result<Option<Implicit>> {
+    use crate::curve::curve::{Curve, CurveKind};
+    use crate::surface::surface::SurfaceKind;
+    match net.kind() {
+        SurfaceKind::Torus(_) => Ok(AnalyticSupport::of(net)
+            .as_ref()
+            .map(AnalyticSupport::implicit)),
+        SurfaceKind::SurfOfRevolution(rev) => {
+            let Some(CurveKind::BSpline(arc)) = rev.gen_curve().map(Curve::kind) else {
+                return Ok(None);
+            };
+            let (points, _) = arc.control_polygon()?;
+            let [p0, p1, p2, ..] = points.as_slice() else {
+                return Ok(None);
+            };
+            let mid = (*p0 + *p2) * 0.5;
+            let (a, h) = ((*p0 - mid).length(), (mid - *p1).length());
+            if h.partial_cmp(&0.0) != Some(core::cmp::Ordering::Greater) {
+                return Ok(None);
+            }
+            let centre = mid + (mid - *p1) * (a * a / (h * h));
+            let position = rev.position();
+            let (o, axis) = (
+                position.origin(),
+                position.x_axis().cross(position.y_axis()),
+            );
+            let along = (centre - o).dot(axis);
+            let across = (centre - o) - axis * along;
+            Ok(Some(Implicit::Torus {
+                centre: o + axis * along,
+                axis,
+                major: across.length(),
+                minor: (*p0 - centre).length(),
+            }))
+        }
+        _ => Ok(None),
+    }
+}
 
 /// The exact class of a round whose net is its definition.
 #[derive(Debug, Clone, Copy, PartialEq, Eq)]
@@ -263,7 +314,8 @@ impl ExactRound {
 
     /// Re-certifies `net` against this definition: the supports classified again are of
     /// this class, and the net lies on the class's surface derived from them
-    /// ([`check_exact_net`]). Returns the distance proved, 0.
+    /// ([`check_exact_net`]). A round of the circle class is certified on its own torus
+    /// instead ([`ExactRound::supports_distance`]). Returns the distance proved, 0.
     ///
     /// # Errors
     ///
@@ -271,6 +323,10 @@ impl ExactRound {
     /// surface from them, a net off it; those of reading the supports' quadrics and the
     /// net's NURBS form.
     pub fn recertify(&self, k: &Kernel, net: &Surface) -> Result<f64> {
+        if self.class == ExactClass::Circle {
+            self.supports_distance(k, net)?;
+            return Ok(0.0);
+        }
         let [a, b] = &self.supports;
         let blend = classify_supports(k, [a, b], &self.inputs)?;
         // the chord (#2486): a line spine between planes, of sizes that give its offsets
@@ -302,19 +358,72 @@ impl ExactRound {
         Ok(0.0)
     }
 
-    /// How far `net` strays from the surface of the class its supports give
-    /// ([`net_miss`]); `None` where they give it none.
+    /// A round of the circle class: the torus `net` is ([`own_torus`]: its revolution's
+    /// axis and generating circle, or its torus label), proved to hold its form to rounding
+    /// ([`lies_on`]),
+    /// and how far that torus is from the class's torus derived from the supports
+    /// ([`class_surface`]): returned, refused past the tolerance the class was read within.
+    ///
+    /// The class is read within the body's tolerance. Supports only that close to it (a plane
+    /// 1e-13 off square to a cylinder's axis, read from a record) have no torus touching both
+    /// exactly, so the round the solver builds, tangent to the supports as they are, is on a
+    /// torus of its own. Its definition is that torus, at bound 0, a stated distance from the
+    /// supports' class.
     ///
     /// # Errors
     ///
-    /// Those of reading the supports' quadrics and of the evaluations.
-    pub fn miss(&self, k: &Kernel, net: &Surface) -> Result<Option<f64>> {
+    /// `InvalidInput` for supports of another class or with no torus, a net of no torus, a
+    /// net off its torus, or a torus further from the class's than the tolerance; those of
+    /// reading the supports' quadrics, the net's NURBS form and [`lies_on`].
+    pub fn supports_distance(&self, k: &Kernel, net: &Surface) -> Result<f64> {
+        let refused = |why| SolidError::invalid_input(OP_CIRCLE, why);
         let [a, b] = &self.supports;
         let blend = classify_supports(k, [a, b], &self.inputs)?;
-        match class_surface(&blend, self.leader.as_ref(), &self.map) {
-            Some(surface) => net_miss(k, net, &surface).map(Some),
-            None => Ok(None),
+        if ExactClass::of_spine(blend.spine) != Some(ExactClass::Circle) {
+            return Err(refused(
+                "the supports classified again are not of the class the definition names",
+            ));
         }
+        let Some(ClassSurface::Implicit(Implicit::Torus {
+            centre,
+            axis,
+            major,
+            minor,
+        })) = class_surface(&blend, self.leader.as_ref(), &self.map)
+        else {
+            return Err(refused("the supports give the class no torus"));
+        };
+        let form = crate::blend::fillet_section::station::as_b_spline(k, net)?
+            .ok_or_else(|| refused("the net has no NURBS form to check"))?;
+        let own = own_torus(net)?.ok_or_else(|| refused("the net is no torus"))?;
+        if !lies_on(&form, &own)? {
+            return Err(refused("the net is off its own torus"));
+        }
+        let Implicit::Torus {
+            centre: c,
+            axis: d,
+            major: big_r,
+            minor: r,
+        } = own
+        else {
+            return Err(refused("the net's own surface is no torus"));
+        };
+        // every point of the net's torus moves to the class's torus by at most: the turn of
+        // its axis (an angle, over the torus's reach from its centre), the move of its
+        // centre, and the change of its radii
+        let turn = d.cross(axis).length().atan2(d.dot(axis).abs());
+        let distance =
+            turn * (big_r + r) + (c - centre).length() + (big_r - major).abs() + (r - minor).abs();
+        let within = matches!(
+            distance.partial_cmp(&self.inputs.tol),
+            Some(core::cmp::Ordering::Less | core::cmp::Ordering::Equal)
+        );
+        if !within {
+            return Err(refused(
+                "the net's torus is further from the class's than the tolerance",
+            ));
+        }
+        Ok(distance)
     }
 
     /// Moves the definition with its net (`Transform`): the supports, the leader and the

solidcore/crates/solidcore/src/surface/torus.rs

@@ -1061,7 +1061,6 @@ impl Torus {
         let ref_frame = ref_iso.reference_frame.ok_or_else(missing)?;
         let other_frame = other_iso.reference_frame.ok_or_else(missing)?;
         let min_rad = other_iso.radius.ok_or_else(missing)?;
-        let sweep_ivl = ref_iso.angles.ok_or_else(missing)?;
         let common_point = other_iso.common_line_point.ok_or_else(missing)?;
         let common_vec = other_iso.common_line_vec.ok_or_else(missing)?;
 
@@ -1069,6 +1068,17 @@ impl Torus {
         let torus_origin =
             common_point + common_vec.dot(other_frame.origin() - common_point) * common_vec;
 
+        // not upstream (#3199): the x axis and the sweep interval from the form's own first
+        // sweep row's ends about the axis (`sweep_frame`), not from `is_arc`'s fit of it
+        let Some((torus_x_axis, torus_y_axis, sweep_ivl)) = sweep_frame(
+            torus_origin,
+            ref_frame.z_axis(),
+            fitted_sides(&ref_frame, ref_iso.angles.ok_or_else(missing)?),
+            meridians(test, swap_uv)?,
+        ) else {
+            return Ok(None);
+        };
+
         // get torus orientation and size
         // not upstream (WR-501, #1470, opt-in; tangency-labs/smlib@6136a04): the major
         // radius runs from the torus origin, not from the first sweep circle's centre,
@@ -1079,9 +1089,7 @@ impl Torus {
             other_frame.origin() - ref_frame.origin()
         };
         let maj_rad = major_vec.length();
-        let torus_x_axis = ref_frame.x_axis(); // the frame where all circles are parallel
-        let torus_y_axis = ref_frame.y_axis();
-        let torus_z_axis = ref_frame.z_axis();
+        let torus_z_axis = ref_frame.z_axis(); // the frame where all circles are parallel
 
         // torus must have positive radii
         if min_rad <= EFF_ZERO || maj_rad <= EFF_ZERO {
@@ -1294,7 +1302,6 @@ impl Torus {
         let mut ref_frame = ref_iso.reference_frame.ok_or_else(missing)?;
         let other_frame = other_iso.reference_frame.ok_or_else(missing)?;
         let min_rad = other_iso.radius.ok_or_else(missing)?;
-        let sweep_ivl = ref_iso.angles.ok_or_else(missing)?;
         let common_point = other_iso.common_line_point.ok_or_else(missing)?;
         let common_vec = other_iso.common_line_vec.ok_or_else(missing)?;
 
@@ -1302,8 +1309,18 @@ impl Torus {
         let torus_origin =
             common_point + common_vec.dot(other_frame.origin() - common_point) * common_vec;
 
+        // not upstream (#3199): the x axis and the sweep interval from the form's own first
+        // sweep row's ends about the axis (`sweep_frame`), not from `is_arc`'s fit of it
+        let (x_axis, y_axis, sweep_ivl) = sweep_frame(
+            torus_origin,
+            ref_frame.z_axis(),
+            fitted_sides(&ref_frame, ref_iso.angles.ok_or_else(missing)?),
+            meridians(srf, swap_uv)?,
+        )
+        .ok_or_else(|| SolidError::sm_err(OP_REBUILD))?;
+
         // Install the origin into pRefFrame.  // [B328]
-        let _ = ref_frame.set_canonical(torus_origin, ref_frame.x_axis(), ref_frame.y_axis());
+        let _ = ref_frame.set_canonical(torus_origin, x_axis, y_axis);
 
         // get torus orientation and size
         let major_vec = other_frame.origin() - ref_frame.origin();
@@ -1587,6 +1604,82 @@ impl Revolved for Torus {
 
 /// `SmTorus::CreateIsoParametricCurve` of the surface of revolution `rev` of a torus
 /// of `shape` (also called from the base class's rebuilds, where it is virtual).
+/// The control points, Euclidean, of the form's first and last meridians: the columns of its
+/// net at the sweep's two ends (the sweep is along u, or along v with `swap_uv`). A clamped
+/// end's column is that boundary curve's own control polygon, and a meridian lies in a plane
+/// through the axis, so every one of its control points has the meridian's angle about the
+/// axis.
+fn meridians(srf: &BSplineSurface, swap_uv: bool) -> Result<[Vec<Point3>; 2]> {
+    let nurb = srf
+        .gw_nurb()
+        .ok_or_else(|| SolidError::null_pointer(OP_NURB_TORUS))?;
+    let (n, m) = (nurb.n(), nurb.m());
+    let column = |end: bool| -> Result<Vec<Point3>> {
+        if swap_uv {
+            let j = if end { m } else { 0 };
+            (0..=n).map(|i| Ok(nurb.cpt(i, j)?.to_euclid())).collect()
+        } else {
+            let i = if end { n } else { 0 };
+            (0..=m).map(|j| Ok(nurb.cpt(i, j)?.to_euclid())).collect()
+        }
+    };
+    Ok([column(false)?, column(true)?])
+}
+
+/// The directions, about the axis, of the fitted first sweep row's ends (`is_arc`'s frame and
+/// angles): [`sweep_frame`] reads only their sides from them.
+fn fitted_sides(frame: &Axis2Placement, angles: Extent1d) -> [Vec3; 2] {
+    let at = |deg: f64| {
+        let a = deg * PI / 180.0;
+        a.cos() * frame.x_axis() + a.sin() * frame.y_axis()
+    };
+    [at(angles.min()), at(angles.max())]
+}
+
+/// Not upstream (#3199): a recognised torus's x axis, y axis and sweep interval (degrees) from
+/// its form's first and last meridians (`meridians`: their control points, exact data), about
+/// the axis through `origin` along the unit `z`. A meridian lies in a plane through the axis,
+/// so its control points share its angle on the side the tube's centre is on, and the
+/// opposite angle past the axis (a spindle torus's tube crosses it); `sides` (the fitted
+/// row's end directions, [`fitted_sides`]) says which side that is, and nothing else. Each
+/// meridian's angle is read at its control point farthest from the axis on that side: x is
+/// the first meridian's direction square to the axis (so the interval starts at 0), y `z x
+/// x`, and the interval's end the last meridian's angle from x about z (`atan2` in that
+/// frame), a closed sweep (within `SM_EFF_ZERO_SQRT` of 0) at 360, as `is_arc` maps its end.
+/// `is_arc` fitted the first sweep row's circle through three samples, so on a short arc its
+/// centre, and with it the x axis and the end angle, was off the axis by its conditioning
+/// (1e-14 rad on #3185's grooves); the converter centres the sweep circle on the axis, so its
+/// ends and its inside now agree to rounding. `None` when a meridian has no control point off
+/// the axis on its side.
+fn sweep_frame(
+    origin: Point3,
+    z: Vec3,
+    sides: [Vec3; 2],
+    [first, last]: [Vec<Point3>; 2],
+) -> Option<(Vec3, Vec3, Extent1d)> {
+    let radial = |p: &Point3| {
+        let s = *p - origin;
+        s - s.dot(z) * z
+    };
+    let farthest = |ps: &[Point3], side: Vec3| {
+        ps.iter()
+            .map(radial)
+            .filter(|r| r.dot(side) > 0.0)
+            .max_by(|a, b| a.length().total_cmp(&b.length()))
+    };
+    let [first_side, last_side] = sides;
+    let x = farthest(&first, first_side)?.unitize().ok()?;
+    let y = z.cross(x);
+    let e = farthest(&last, last_side)?;
+    let mut end_deg = e.dot(y).atan2(e.dot(x)) * 180.0 / PI;
+    if end_deg.abs() < EFF_ZERO_SQRT {
+        end_deg = 360.0;
+    } else if end_deg < 0.0 {
+        end_deg += 360.0;
+    }
+    Some((x, y, Extent1d::from_min_max(0.0, end_deg)))
+}
+
 fn iso_curve(
     rev: &SurfOfRevolution,
     shape: Shape,

solidcore/crates/solidcore/tests/contact_between_bodies.rs

@@ -0,0 +1,200 @@
+//! #3144: whether two bodies meet, decided by a certified bracket of their boundary distance
+//! (`query::contact`), against closed forms: boxes touching on a face, apart by 1e-9,
+//! overlapping in a face, an edge resting on a face, interpenetrating, one inside the other,
+//! far apart; spheres apart and crossing; a cylinder crossing a box's face and one against a
+//! side along a line; a peg through a hole that it clears.
+
+#![allow(clippy::unwrap_used, clippy::panic, clippy::float_cmp)] // test
+
+use solidcore::boolean::api_brep::boolean_difference;
+use solidcore::geom::certified::Verdict;
+use solidcore::kernel::Kernel;
+use solidcore::primitive::api_primitives::{create_box, create_cylinder, create_sphere};
+use solidcore::query::contact::{boundary_distance, meets};
+use solidcore::topology::brep::mirror::transform_matrix;
+use solidcore::{Brep, Point3};
+
+fn k() -> Kernel {
+    Kernel {
+        fixes: solidcore::kernel::UpstreamFixes::all(),
+        ..Kernel::default()
+    }
+}
+
+fn cube(k: &mut Kernel, at: [f64; 3], size: f64) -> Brep {
+    create_box(k, Point3::new(at[0], at[1], at[2]), size, size, size).unwrap()
+}
+
+/// `b` turned by `angle` about the x axis, then moved by `by`.
+fn turned_about_x(k: &mut Kernel, b: &mut Brep, angle: f64, by: [f64; 3]) {
+    let (c, s) = (angle.cos(), angle.sin());
+    // row vectors (`p' = p A + t`): A's rows are the images of x, y and z, t the fourth row
+    let m = [
+        [1.0, 0.0, 0.0, 0.0],
+        [0.0, c, s, 0.0],
+        [0.0, -s, c, 0.0],
+        [by[0], by[1], by[2], 1.0],
+    ];
+    transform_matrix(k, b, &m, true).unwrap();
+}
+
+/// Asserts the verdict at `tol`, the bracket shown on failure.
+fn verdict_is(k: &Kernel, a: &mut Brep, b: &mut Brep, tol: f64, want: Verdict) {
+    let c = meets(k, a, b, tol).unwrap();
+    assert!(c.lower <= c.upper, "{c:?}");
+    assert_eq!(c.verdict, want, "at {tol}: {c:?}");
+}
+
+/// The bracket of the boundary distance, asserted to hold `exact`.
+fn holds(k: &Kernel, a: &mut Brep, b: &mut Brep, exact: f64) {
+    let c = boundary_distance(k, a, b, 1e-6).unwrap();
+    assert_eq!(c.verdict, Verdict::Yes, "{c:?}");
+    assert!(c.lower <= exact && exact <= c.upper, "{exact} not in {c:?}");
+    assert!(c.upper - c.lower <= 1e-6 * c.upper.max(1.0), "{c:?}");
+}
+
+#[test]
+fn boxes_touching_on_a_face_meet_at_any_tolerance() {
+    let mut k = k();
+    let mut a = cube(&mut k, [0.0; 3], 1.0);
+    let mut b = cube(&mut k, [1.0, 0.0, 0.0], 1.0);
+    let c = meets(&k, &mut a, &mut b, 1e-12).unwrap();
+    assert_eq!(c.verdict, Verdict::Yes, "{c:?}");
+    let [p, q] = c.witness.unwrap();
+    assert!(
+        p.distance(q) <= 1e-12 && (p.x - 1.0).abs() <= 1e-12,
+        "{c:?}"
+    );
+}
+
+#[test]
+fn boxes_apart_by_a_nanometre_meet_only_at_a_tolerance_past_it() {
+    let mut k = k();
+    // the gap as the doubles have it
+    let gap = (1.0 + 1e-9) - 1.0;
+    let mut a = cube(&mut k, [0.0; 3], 1.0);
+    let mut b = cube(&mut k, [1.0 + 1e-9, 0.0, 0.0], 1.0);
+    verdict_is(&k, &mut a, &mut b, 1e-8, Verdict::Yes);
+    verdict_is(&k, &mut a, &mut b, 1e-10, Verdict::No);
+    holds(&k, &mut a, &mut b, gap);
+}
+
+#[test]
+fn boxes_overlapping_in_a_face_meet() {
+    let mut k = k();
+    let mut a = cube(&mut k, [0.0; 3], 1.0);
+    let mut b = cube(&mut k, [1.0, 0.5, 0.5], 1.0);
+    verdict_is(&k, &mut a, &mut b, 1e-12, Verdict::Yes);
+}
+
+#[test]
+fn an_edge_resting_on_a_face_meets_it() {
+    let mut k = k();
+    let mut a = cube(&mut k, [0.0; 3], 1.0);
+    // a cube centred on the x axis, turned 45 degrees about it (either way: it is square), is
+    // a diamond whose lowest edge is at z = -sqrt(2)/2; moved onto a's top face
+    let mut b = cube(&mut k, [0.0, -0.5, -0.5], 1.0);
+    let half = std::f64::consts::FRAC_1_SQRT_2;
+    turned_about_x(
+        &mut k,
+        &mut b,
+        std::f64::consts::FRAC_PI_4,
+        [0.0, 0.5, 1.0 + half],
+    );
+    let c = meets(&k, &mut a, &mut b, 1e-12).unwrap();
+    assert_eq!(c.verdict, Verdict::Yes, "{c:?}");
+    let [p, _] = c.witness.unwrap();
+    assert!((p.z - 1.0).abs() <= 1e-9, "{c:?}");
+}
+
+#[test]
+fn interpenetrating_boxes_meet() {
+    let mut k = k();
+    let mut a = cube(&mut k, [0.0; 3], 1.0);
+    let mut b = cube(&mut k, [0.5; 3], 1.0);
+    verdict_is(&k, &mut a, &mut b, 1e-12, Verdict::Yes);
+}
+
+#[test]
+fn a_box_inside_another_meets_it_and_their_walls_are_one_apart() {
+    let mut k = k();
+    let mut a = cube(&mut k, [0.0; 3], 3.0);
+    let mut b = cube(&mut k, [1.0; 3], 1.0);
+    verdict_is(&k, &mut a, &mut b, 1e-3, Verdict::Yes);
+    verdict_is(&k, &mut b, &mut a, 1e-3, Verdict::Yes);
+    holds(&k, &mut a, &mut b, 1.0);
+}
+
+#[test]
+fn boxes_far_apart_do_not_meet_and_their_gap_is_bracketed() {
+    let mut k = k();
+    let mut a = cube(&mut k, [0.0; 3], 1.0);
+    let mut b = cube(&mut k, [5.0, 0.0, 0.0], 1.0);
+    verdict_is(&k, &mut a, &mut b, 0.1, Verdict::No);
+    holds(&k, &mut a, &mut b, 4.0);
+}
+
+#[test]
+fn spheres_apart_are_their_centres_gap_less_their_radii() {
+    let mut k = k();
+    let mut a = create_sphere(&mut k, Point3::new(0.0, 0.0, 0.0), 1.0).unwrap();
+    let mut b = create_sphere(&mut k, Point3::new(3.0, 0.0, 0.0), 1.0).unwrap();
+    verdict_is(&k, &mut a, &mut b, 0.5, Verdict::No);
+    verdict_is(&k, &mut a, &mut b, 2.0, Verdict::Yes);
+    holds(&k, &mut a, &mut b, 1.0);
+}
+
+#[test]
+fn crossing_spheres_meet_where_their_surfaces_cross() {
+    let mut k = k();
+    let mut a = create_sphere(&mut k, Point3::new(0.0, 0.0, 0.0), 1.0).unwrap();
+    let mut b = create_sphere(&mut k, Point3::new(1.5, 0.0, 0.0), 1.0).unwrap();
+    let c = meets(&k, &mut a, &mut b, 1e-9).unwrap();
+    assert_eq!(c.verdict, Verdict::Yes, "{c:?}");
+    // the witness is on the circle where they cross: x = 0.75, |p| = 1
+    let [p, _] = c.witness.unwrap();
+    assert!((p.x - 0.75).abs() <= 1e-6, "{c:?}");
+}
+
+#[test]
+fn a_cylinder_crossing_a_box_face_meets_the_box() {
+    let mut k = k();
+    let mut a = cube(&mut k, [0.0; 3], 1.0);
+    let mut b = create_cylinder(&mut k, Point3::new(0.5, 0.5, 0.5), 0.3, 1.0).unwrap();
+    verdict_is(&k, &mut a, &mut b, 1e-9, Verdict::Yes);
+}
+
+/// A cylinder touching a box side along a line meets it; moved 0.25 off, its bracket holds
+/// 0.25 and the verdicts either side of it are exact.
+///
+/// The bracket is not 1e-6 wide here: trimmed rim pieces stall `lower` at 0.25 - 0.0013 (the
+/// end caps' pieces at the rim straddle the circular trim, and each hull carries a sliver of
+/// the cap's plane outside the disk toward the box; their subdivision runs out of pairs
+/// first). The bracket contains the truth; verdicts outside [0.2487, 0.25] are exact. The
+/// width asserted is what the pair budget certifies, not a target.
+#[test]
+fn a_cylinder_against_a_box_side_meets_it_along_a_line() {
+    let mut k = k();
+    let mut a = cube(&mut k, [0.0; 3], 1.0);
+    let mut b = create_cylinder(&mut k, Point3::new(1.5, 0.5, 0.0), 0.5, 1.0).unwrap();
+    verdict_is(&k, &mut a, &mut b, 1e-9, Verdict::Yes);
+    let mut far = create_cylinder(&mut k, Point3::new(1.75, 0.5, 0.0), 0.5, 1.0).unwrap();
+    verdict_is(&k, &mut a, &mut far, 0.24, Verdict::No);
+    verdict_is(&k, &mut a, &mut far, 0.26, Verdict::Yes);
+    let c = boundary_distance(&k, &mut a, &mut far, 1e-6).unwrap();
+    assert!(c.lower <= 0.25 && 0.25 <= c.upper, "{c:?}");
+    assert!(c.lower >= 0.2487 && c.upper - 0.25 <= 1e-12, "{c:?}");
+}
+
+#[test]
+fn a_peg_through_a_hole_it_clears_does_not_meet_the_plate() {
+    let mut k = k();
+    // a 4 x 4 x 1 plate with a hole of radius 1 through its middle, and a peg of radius 0.75
+    // down the hole's axis: the peg clears the hole's wall by 0.25
+    let plate = create_box(&mut k, Point3::new(0.0, 0.0, 0.0), 4.0, 4.0, 1.0).unwrap();
+    let hole = create_cylinder(&mut k, Point3::new(2.0, 2.0, -1.0), 1.0, 3.0).unwrap();
+    let mut plate = boolean_difference(&mut k, plate, hole).unwrap();
+    let mut peg = create_cylinder(&mut k, Point3::new(2.0, 2.0, -2.0), 0.75, 5.0).unwrap();
+    verdict_is(&k, &mut plate, &mut peg, 0.1, Verdict::No);
+    holds(&k, &mut plate, &mut peg, 0.25);
+}

solidcore/docs/definitions-and-certificates.md

@@ -126,6 +126,17 @@ give the same class, and the net lies on that class's surface derived from the s
 parametrisation's: the net reversed, swapped, split or trimmed is still on S, so those edits
 carry it.
 
+A round of the `Circle` class is certified on its own torus. The class is read within the body's
+tolerance, so supports only that close to it (a plane 1e-13 off square to a cylinder's axis, as a
+record can carry them) have no torus that touches both exactly, and the round the solver builds,
+tangent to the supports as they are, lies on a torus of its own. Its definition is that torus: the
+torus of the net's own form (`own_torus`: a revolution's axis and generating arc, or a torus's
+label), which the form is proved to lie on (`lies_on`), at bound 0; and the distance from it to the
+class's torus derived from the supports, stated (`ExactRound::supports_distance`: the axes' angle
+times the reach, the centres' distance and the radii's differences) and refused past the tolerance
+the class was read within. Supports exactly of the class give a distance at rounding. The supports
+are never moved to fit.
+
 ### 3.4 The blend definition
 
 `ProceduralBlend` holds the net (`net`, a `BSplineSurface`) and the definition