query + ops: whether two bodies meet, on a certified bracket of their boundary distance (#3144) (#13554)

query::contact decides whether two bodies meet from a certified bracket [lower, upper] of the distance between their boundaries; solidcore-ops exposes it as Body::meets, meets_each and boundary_distance. Undecided carries the bracket and is never read as either answer.

Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

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,

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/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/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);
+}