boolean: a grazing edge-by-face vertex at its three faces' common point; certify the continued edge (#3062)

Rebased onto main (its net change from 9b559c4c07, three-way per file).

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

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

@@ -78,7 +78,10 @@ fn taken(r: f64) -> f64 {
 
 #[test]
 fn the_upright_round_runs_on_until_it_leaves_the_material() {
-    for r in [2.0, 1.0] {
+    // r 0.48 to 0.65 (#3062): where the round crosses the seam of the r 5 and r 3 rounds,
+    // near (60, 40 - r, 17), the r 5 round's curve with the upright meets the r 3 round at a
+    // grazing angle, so the vertex there is the three faces' common point
+    for r in [2.0, 1.0, 0.7, 0.65, 0.6, 0.55, 0.5, 0.48, 0.45, 0.42] {
         let b = rounded();
         let v0 = volume(&b);
         let e = edge_at(&b, [60.0, 40.0, 8.0]);
@@ -91,3 +94,62 @@ fn the_upright_round_runs_on_until_it_leaves_the_material() {
         assert!(rel(volume(&b), want) < 1e-8, "r {r}: {} {want}", volume(&b));
     }
 }
+
+/// #3062: at r 0.5 the edges at the vertex where the upright crosses the seam of the r 5 and
+/// r 3 rounds lie on both their faces' surfaces along their last stretch (the continued
+/// span, 1.2e-4 long, and on), within each edge's tolerance: measured at 200 points of the
+/// last 1e-3 of each edge's length, by dropping each point on each face's surface.
+#[test]
+fn the_edges_at_the_grazing_crossing_lie_on_their_faces() {
+    use solidcore::surface::surface::p4::NoObject;
+    use solidcore::surface::surface::p5::DropPointOptions;
+    use solidcore::topology::{edge, edgeuse, face, vertex};
+    let r = 0.5;
+    let b = rounded();
+    let e = edge_at(&b, [60.0, 40.0, 8.0]);
+    blend(&b, Blend::Round, &[e], &[r], &[], Corners::Filled).unwrap();
+    let (k, brep) = kernel_copy(&b);
+    let near = solidcore::math::Point3::new(60.0, 40.0 - r, 17.0);
+    let mut seen = 0;
+    for &e in brep.edge_list() {
+        let curve = edge::curve(&brep, e).unwrap().unwrap();
+        let ivl = edge::interval(&brep, e).unwrap();
+        let ends = [
+            (edge::start_vertex(&brep, e).unwrap(), ivl.min()),
+            (edge::end_vertex(&brep, e).unwrap(), ivl.max()),
+        ];
+        for (v, at) in ends {
+            let Some(v) = v else { continue };
+            let p = vertex::point(&brep, v).unwrap();
+            if p.distance(near) > 0.2 || (p.z - 17.0).abs() < 1e-6 {
+                continue;
+            }
+            seen += 1;
+            let tol = edge::tolerance(&brep, e).unwrap();
+            let span = 1e-3 * ivl.length();
+            for eu in edge::edgeuses(&brep, e, None, None).unwrap() {
+                let f = edgeuse::p1::face(&brep, eu).unwrap().unwrap();
+                let s = face::p2::surface(&brep, f).unwrap().unwrap();
+                let dom = s.natural_uv_domain_with_nurb(&k).unwrap();
+                let mut worst = 0.0_f64;
+                for i in 0..=200 {
+                    let t = if at == ivl.min() {
+                        at + span * f64::from(i) / 200.0
+                    } else {
+                        at - span * f64::from(i) / 200.0
+                    };
+                    let x = curve.evaluate_point(&k, t).unwrap();
+                    let d = s
+                        .drop_point::<NoObject>(&k, x, &dom, &DropPointOptions::default())
+                        .unwrap();
+                    worst = worst.max(d.gap);
+                }
+                assert!(
+                    worst <= tol,
+                    "{e:?} on {f:?}: {worst:e} off, tolerance {tol:e}"
+                );
+            }
+        }
+    }
+    assert!(seen >= 3, "the crossing's three edges, found {seen}");
+}

solidcore/crates/solidcore/src/error.rs

@@ -709,6 +709,17 @@ pub enum Reason {
     Offset(offset::OffsetFailure),
     /// A NURBS routine failed.
     Nurbs(nurbs::NurbsFailure),
+    /// An intersection edge whose end is continued to its vertex would leave its faces: the
+    /// continued span, `gap` long at its end, lies up to `bound` (certified) off them,
+    /// past the `tolerance` it must stay within (not upstream, #3062).
+    EdgeOffItsFaces {
+        /// How far the edge's end lay from its vertex.
+        gap: f64,
+        /// The certified bound of the continued span's distance from its faces.
+        bound: f64,
+        /// The tolerance it must stay within.
+        tolerance: f64,
+    },
 }
 
 impl Reason {
@@ -716,7 +727,7 @@ impl Reason {
     #[must_use]
     pub fn entity(&self) -> Option<EntityRef> {
         match self {
-            Self::Failed(_) | Self::Nurbs(_) => None,
+            Self::Failed(_) | Self::Nurbs(_) | Self::EdgeOffItsFaces { .. } => None,
             Self::Offset(o) => o.entity(),
         }
     }
@@ -725,7 +736,7 @@ impl Reason {
     /// comparison detail).
     pub(crate) fn comparison_status(&self) -> i32 {
         match self {
-            Self::Failed(_) => SM_ERR,
+            Self::Failed(_) | Self::EdgeOffItsFaces { .. } => SM_ERR,
             Self::Offset(o) => o.comparison_status(),
             Self::Nurbs(_) => nurbs::NurbsFailure::COMPARISON_STATUS,
         }
@@ -736,6 +747,7 @@ impl Reason {
     pub(crate) fn kind_given(&self, status: i32) -> ErrorKind {
         match self {
             Self::Failed(_) => ErrorKind::of_status(status),
+            Self::EdgeOffItsFaces { .. } => ErrorKind::NotWithinTolerance,
             Self::Offset(o) => o.kind(),
             Self::Nurbs(n) => n.kind(),
         }
@@ -746,6 +758,14 @@ impl core::fmt::Display for Reason {
     fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
         match self {
             Self::Failed(what) => f.write_str(what),
+            Self::EdgeOffItsFaces {
+                gap,
+                bound,
+                tolerance,
+            } => write!(
+                f,
+                "an intersection edge continued {gap:e} to its vertex would lie up to {bound:e} off its faces, past their tolerance {tolerance:e}"
+            ),
             Self::Offset(o) => o.fmt(f),
             Self::Nurbs(n) => n.fmt(f),
         }

solidcore/crates/solidcore/src/topology/proto_topology/p1/proto.rs

@@ -23,6 +23,7 @@ use crate::intersect::point_class::{
 use crate::kernel::Kernel;
 use crate::math::{max2, min2, Point2, Point3, Vec3, BIG_DOUBLE, EFF_ZERO};
 use crate::query::face::p2 as face_p2;
+use crate::surface::surface::Surface;
 use crate::topology::brep::p2 as brep_p2;
 use crate::topology::edgeuse::p1 as eu1;
 use crate::topology::topology as topo;
@@ -891,7 +892,121 @@ fn crossing_error(gap: f64, crossing_sine: f64) -> f64 {
     max2(gap, 0.0) / crossing_sine
 }
 
+/// The most Newton steps [`common_point`] takes: a bound on the loop; convergence is
+/// checked at each step, never assumed.
+const COMMON_POINT_STEPS: usize = 32;
+
+/// Not upstream (#3062): the value at `x` of the implicit equation of `a` and its gradient
+/// (a plane, a cylinder or cone, a sphere; `None` for a torus, whose equation is quartic).
+fn implicit_value(a: &crate::surface::implicit::Implicit, x: Point3) -> Option<(f64, Vec3)> {
+    use crate::surface::implicit::Implicit;
+    Some(match *a {
+        Implicit::Plane { origin, normal } => ((x - origin).dot(normal), normal),
+        Implicit::Cone {
+            origin,
+            axis,
+            radius,
+            tan_angle,
+        } => {
+            let u = x - origin;
+            let h = u.dot(axis);
+            let radial = u - axis * h;
+            let r = radius + tan_angle * h;
+            let f = radial.dot(radial) - r * r;
+            let g = (radial - axis * (r * tan_angle)) * 2.0;
+            (f, g)
+        }
+        Implicit::Sphere { centre, radius } => {
+            let u = x - centre;
+            (u.dot(u) - radius * radius, u * 2.0)
+        }
+        Implicit::Torus { .. } => return None,
+    })
+}
+
+/// Not upstream (#3062): the common point of three analytic surfaces nearest `start`, by
+/// Newton's method on their implicit equations: `None` where one is not a plane, cylinder,
+/// cone or sphere, the normals do not span space at a step, or it does not converge in
+/// [`COMMON_POINT_STEPS`] steps to within rounding of `start`'s size.
+fn common_point(surfaces: [&Surface; 3], start: Point3) -> Option<Point3> {
+    use crate::surface::implicit::AnalyticSupport;
+    let of = |s: &Surface| AnalyticSupport::of(s).map(|a| a.implicit());
+    let [s0, s1, s2] = surfaces;
+    let eq = [of(s0)?, of(s1)?, of(s2)?];
+    let size = 1.0 + start.max_dimension();
+    let mut x = start;
+    for _ in 0..COMMON_POINT_STEPS {
+        let mut rows = [(0.0, Vec3::ZERO); 3];
+        for (row, a) in rows.iter_mut().zip(&eq) {
+            *row = implicit_value(a, x)?;
+        }
+        let [(f0, g0), (f1, g1), (f2, g2)] = rows;
+        let det = g0.dot(g1.cross(g2));
+        let scale = g0.length() * g1.length() * g2.length();
+        if det.is_nan() || det.abs() <= EFF_ZERO * scale {
+            return None;
+        }
+        // Cramer's rule for g_i . dx = -f_i
+        let dx = (g1.cross(g2) * -f0 + g2.cross(g0) * -f1 + g0.cross(g1) * -f2) / det;
+        x += dx;
+        if dx.length() <= 4.0 * f64::EPSILON * size {
+            return Some(x);
+        }
+    }
+    None
+}
+
 impl<O: ClassObject> VertexDefinition<O> {
+    /// Not upstream (#3062): where this definition is a crossing of an edge of one B-rep
+    /// (between two faces) and a face of the other, all three on planes, cylinders, cones or
+    /// spheres, it is placed at their common point, solved on the surfaces from its
+    /// position: a curve crossing the edge at a grazing angle fixes the crossing along the
+    /// edge only to the tolerance over the angle's sine, and the curves meeting there each
+    /// come within tolerance of the third face at a different point. Moved only within the
+    /// definition's own tolerance. Returns whether it moved.
+    ///
+    /// # Errors
+    ///
+    /// Those of the topology walks.
+    pub fn place_on_three_surfaces(
+        &mut self,
+        k: &Kernel,
+        geometry: &impl ClassGeometry<O>,
+    ) -> Result<bool> {
+        let (edge_pc, face_pc) = match (self.pc1.point_class(), self.pc2.point_class()) {
+            (PointClassType::Edge, PointClassType::Face) => (self.pc1, self.pc2),
+            (PointClassType::Face, PointClassType::Edge) => (self.pc2, self.pc1),
+            _ => return Ok(false),
+        };
+        let (Some(eo), Some(fo)) = (edge_pc.object(), face_pc.object()) else {
+            return Ok(false);
+        };
+        let (Some((eb, EntityRef::Edge(e))), Some((fb, EntityRef::Face(f)))) =
+            (geometry.topology_entity(eo), geometry.topology_entity(fo))
+        else {
+            return Ok(false);
+        };
+        let [a, b] = edge::faces(eb, e, None)?[..] else {
+            return Ok(false);
+        };
+        let (Some(sa), Some(sb), Some(sf)) = (
+            crate::topology::face::p2::surface(eb, a)?,
+            crate::topology::face::p2::surface(eb, b)?,
+            crate::topology::face::p2::surface(fb, f)?,
+        ) else {
+            return Ok(false);
+        };
+        let Some(x) = common_point([sa, sb, sf], self.position) else {
+            return Ok(false);
+        };
+        let moved = x.distance(self.position);
+        if moved.is_nan() || moved > self.tolerance(k) {
+            return Ok(false);
+        }
+        self.position = x;
+        Ok(true)
+    }
+
     /// `SmVertexDefinition(pCtx, cpPC1, cpPC2, rPosition, dTol, dGap, pOwner)`: copies
     /// the classifications and computes their parameters at `position` when unset
     /// (statuses ignored); a non-positive `gap` becomes the distance between the two
@@ -2939,10 +3054,12 @@ impl<O: ClassObject> ProtoTopologyManager<O> {
         // Placed again by a crossing that fixes it more closely, it moves this
         // ProtoVertex, the average of its VertexDefinitions, with it.
         let old = self.vd(vd)?.position;
-        let moved = self
+        let def = self
             .vertex_def_mut(vd)
-            .ok_or(dangling("vertex definition is deleted"))?
-            .place_by_crossing(k, geometry, pc1, pc2, crossing_sine)?;
+            .ok_or(dangling("vertex definition is deleted"))?;
+        let crossed = def.place_by_crossing(k, geometry, pc1, pc2, crossing_sine)?;
+        // (#3062) and then on its three faces' common point, where it has one
+        let moved = def.place_on_three_surfaces(k, geometry)? || crossed;
         if moved {
             let (new, vd_tol) = {
                 let v = self.vd(vd)?;

solidcore/crates/solidcore/src/topology/proto_topology/p2/imprint_ends.rs

@@ -38,6 +38,7 @@
 use super::{curve_of, dangling, ignore_soft, EdgeDefId, ProtoBrep, ProtoTopologyManager};
 use super::{ProtoVertexId, Result};
 use crate::curve::curve::p2::{Curve, CurveKind};
+use crate::error::{Failure, Reason, SolidError};
 use crate::geom::extent1d::Extent1d;
 use crate::geom::global_solver::SolverOperation;
 use crate::geom::tol::are_same;
@@ -48,10 +49,281 @@ use crate::nurbs::b_spline_curve::p1::{BSplineCurveForm, CanonicalData, KnotType
 use crate::nurbs::b_spline_curve::BSplineCurve;
 use crate::nurbs::b_spline_curve_n_lib::TrimMode;
 
+const OP: &str = "topology::proto_topology::p2::imprint_ends";
+
 /// How far ahead of a curve's end, as a fraction of the interval, a proto vertex may lie
 /// for the curve to be continued to it.
 const MAX_EXTENSION_FRACTION: f64 = 0.05;
 
+/// Not upstream (#3062): a quadric's equation about a point `o`, `f(x) = y.M.y + 2 m.y + c`
+/// with `y = x - o`, and how a bound of `|f|` bounds the distance from the surface
+/// ([`QuadricKind`]).
+struct Quadric {
+    o: Vec3,
+    m: [[f64; 3]; 3],
+    lin: Vec3,
+    c: f64,
+    kind: QuadricKind,
+}
+
+/// A quadric's kind, for its distance bound: a plane's `|f|` (its normal unit), a sphere's
+/// `|f| / r`, a cylinder's or cone's `|f| / R`, `R` its least radius along the span.
+enum QuadricKind {
+    Plane,
+    Sphere {
+        radius: f64,
+    },
+    Cone {
+        axis: Vec3,
+        radius: f64,
+        tan_angle: f64,
+    },
+}
+
+/// The quadric of a plane, cylinder, cone or sphere; `None` for any other surface.
+fn quadric_of(surface: &crate::surface::surface::Surface) -> Option<Quadric> {
+    use crate::surface::implicit::{AnalyticSupport, Implicit};
+    Some(match AnalyticSupport::of(surface)?.implicit() {
+        Implicit::Plane { origin, normal } => Quadric {
+            o: origin,
+            m: [[0.0; 3]; 3],
+            lin: normal * 0.5,
+            c: 0.0,
+            kind: QuadricKind::Plane,
+        },
+        Implicit::Sphere { centre, radius } => Quadric {
+            o: centre,
+            m: [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
+            lin: Vec3::ZERO,
+            c: -radius * radius,
+            kind: QuadricKind::Sphere { radius },
+        },
+        // |y|^2 - (y.a)^2 - (r + t y.a)^2
+        Implicit::Cone {
+            origin,
+            axis,
+            radius,
+            tan_angle,
+        } => {
+            let a = [axis.x, axis.y, axis.z];
+            let s = 1.0 + tan_angle * tan_angle;
+            let mut m = [[0.0; 3]; 3];
+            for (row, ai) in m.iter_mut().zip(a) {
+                for (x, aj) in row.iter_mut().zip(a) {
+                    *x = -s * ai * aj;
+                }
+            }
+            for (i, row) in m.iter_mut().enumerate() {
+                if let Some(x) = row.get_mut(i) {
+                    *x += 1.0;
+                }
+            }
+            Quadric {
+                o: origin,
+                m,
+                lin: axis * (-radius * tan_angle),
+                c: -radius * radius,
+                kind: QuadricKind::Cone {
+                    axis,
+                    radius,
+                    tan_angle,
+                },
+            }
+        }
+        Implicit::Torus { .. } => return None,
+    })
+}
+
+/// `C(n, k)`.
+fn binomial(n: usize, k: usize) -> f64 {
+    let mut r = 1.0_f64;
+    for i in 0..k {
+        #[allow(clippy::cast_precision_loss)]
+        let f = (n - i) as f64 / (i + 1) as f64;
+        r *= f;
+    }
+    r
+}
+
+/// The product of two polynomials in Bernstein form over `[0, 1]`, of degrees
+/// `a.len() - 1` and `b.len() - 1`.
+fn bernstein_product(a: &[f64], b: &[f64]) -> Vec<f64> {
+    let (p, q) = (a.len().saturating_sub(1), b.len().saturating_sub(1));
+    let mut c = vec![0.0; p + q + 1];
+    for (i, x) in a.iter().enumerate() {
+        for (j, y) in b.iter().enumerate() {
+            if let Some(slot) = c.get_mut(i + j) {
+                *slot += binomial(p, i) * binomial(q, j) / binomial(p + q, i + j) * x * y;
+            }
+        }
+    }
+    c
+}
+
+/// De Casteljau's split of homogeneous control points at `u`: the part over `[u, 1]` when
+/// `upper`, else over `[0, u]`.
+fn de_casteljau(pts: &[[f64; 4]], u: f64, upper: bool) -> Vec<[f64; 4]> {
+    let mut row = pts.to_vec();
+    let mut lower = Vec::new();
+    let mut higher = Vec::new();
+    while let (Some(first), Some(last)) = (row.first().copied(), row.last().copied()) {
+        lower.push(first);
+        higher.push(last);
+        if row.len() == 1 {
+            break;
+        }
+        row = row
+            .windows(2)
+            .map(|w| {
+                let p = w.first().copied().unwrap_or_default();
+                let q = w.get(1).copied().unwrap_or_default();
+                [
+                    p[0] + u * (q[0] - p[0]),
+                    p[1] + u * (q[1] - p[1]),
+                    p[2] + u * (q[2] - p[2]),
+                    p[3] + u * (q[3] - p[3]),
+                ]
+            })
+            .collect();
+    }
+    if upper {
+        higher.reverse();
+        higher
+    } else {
+        lower
+    }
+}
+
+/// Not upstream (#3062): a certified bound of how far `curve` over `[from, to]` (inside one
+/// of its spans) lies from each of `surfaces` (planes, cylinders, cones, spheres), with no
+/// sampling. The span's Bézier form, its homogeneous control points `(X, W)` cut to the
+/// interval by de Casteljau, composes with each quadric to the polynomial
+/// `F = Y.M.Y + 2 m.Y W + c W^2`, `Y = X - o W`, in Bernstein form. `|f| = |F| / W^2` is at
+/// most the largest coefficient's size, with their rounding (`gamma`), over the least weight
+/// coefficient squared, and the distance follows ([`QuadricKind`]). `None` where a surface
+/// is of another kind, the curve is no B-spline, a weight coefficient is not positive, or a
+/// cone's radius along the span is not bounded away from zero.
+fn span_off_surfaces(
+    curve: &Curve,
+    from: f64,
+    to: f64,
+    surfaces: [&crate::surface::surface::Surface; 2],
+) -> Option<f64> {
+    let CurveKind::BSpline(bs) = curve.kind() else {
+        return None;
+    };
+    let (lo, hi) = (from.min(to), from.max(to));
+    let segments = bs.gw_nurb()?.decompose_bez().ok()?;
+    let segment = segments.iter().find(|seg| {
+        let kn = seg.knots();
+        kn.first()
+            .is_some_and(|a| *a <= lo + EFF_ZERO * (1.0 + lo.abs()))
+            && kn
+                .last()
+                .is_some_and(|b| *b >= hi - EFF_ZERO * (1.0 + hi.abs()))
+    })?;
+    let (a, b) = (*segment.knots().first()?, *segment.knots().last()?);
+    if !(b > a) {
+        return None;
+    }
+    let mut pts: Vec<[f64; 4]> = segment
+        .cpts()
+        .iter()
+        .map(|c| {
+            if c.w == crate::nurbs::math::NOW {
+                [c.x, c.y, c.z, 1.0]
+            } else {
+                [c.x, c.y, c.z, c.w]
+            }
+        })
+        .collect();
+    let (u0, u1) = ((lo - a) / (b - a), (hi - a) / (b - a));
+    pts = de_casteljau(&pts, u0.clamp(0.0, 1.0), true);
+    let u = if u0 < 1.0 {
+        ((u1 - u0) / (1.0 - u0)).clamp(0.0, 1.0)
+    } else {
+        1.0
+    };
+    pts = de_casteljau(&pts, u, false);
+    let w: Vec<f64> = pts.iter().map(|p| p[3]).collect();
+    let w_min = w.iter().copied().fold(f64::INFINITY, f64::min);
+    if !(w_min > 0.0) {
+        return None;
+    }
+    let n = pts.len().saturating_sub(1);
+    let mut worst = 0.0_f64;
+    for surface in surfaces {
+        let q = quadric_of(surface)?;
+        let y: Vec<[f64; 3]> = pts
+            .iter()
+            .map(|p| {
+                [
+                    p[0] - q.o.x * p[3],
+                    p[1] - q.o.y * p[3],
+                    p[2] - q.o.z * p[3],
+                ]
+            })
+            .collect();
+        let comp =
+            |i: usize| -> Vec<f64> { y.iter().map(|v| v.get(i).copied().unwrap_or(0.0)).collect() };
+        let mut f = vec![0.0; 2 * n + 1];
+        let mut f_abs = vec![0.0; 2 * n + 1];
+        let mut add = |terms: Vec<f64>, factor: f64| {
+            for ((slot, abs), t) in f.iter_mut().zip(f_abs.iter_mut()).zip(terms) {
+                *slot += factor * t;
+                *abs += (factor * t).abs();
+            }
+        };
+        for (i, row) in q.m.iter().enumerate() {
+            for (j, mij) in row.iter().enumerate() {
+                if *mij != 0.0 {
+                    add(bernstein_product(&comp(i), &comp(j)), *mij);
+                }
+            }
+        }
+        for (i, li) in [q.lin.x, q.lin.y, q.lin.z].iter().enumerate() {
+            if *li != 0.0 {
+                add(bernstein_product(&comp(i), &w), 2.0 * li);
+            }
+        }
+        if q.c != 0.0 {
+            add(bernstein_product(&w, &w), q.c);
+        }
+        let rounding =
+            crate::nurbs::math::gamma(8 * (n + 1) + 16) * f_abs.iter().copied().fold(0.0, f64::max);
+        let f_bound = f.iter().fold(0.0_f64, |m, x| m.max(x.abs())) + rounding;
+        let f_over = f_bound / (w_min * w_min);
+        let off = match q.kind {
+            QuadricKind::Plane => f_over,
+            QuadricKind::Sphere { radius } => f_over / radius,
+            QuadricKind::Cone {
+                axis,
+                radius,
+                tan_angle,
+            } => {
+                // the radius along the span, `r + t (x - o).a`, at least its least over the
+                // Euclidean control points (the span is their convex combination)
+                let along =
+                    |v: &[f64; 3], wi: f64| (v[0] * axis.x + v[1] * axis.y + v[2] * axis.z) / wi;
+                let mut h_lo = f64::INFINITY;
+                let mut h_hi = f64::NEG_INFINITY;
+                for (v, wi) in y.iter().zip(&w) {
+                    let h = along(v, *wi);
+                    h_lo = h_lo.min(h);
+                    h_hi = h_hi.max(h);
+                }
+                let r_min = (radius + tan_angle * h_lo).min(radius + tan_angle * h_hi);
+                if !(r_min > 0.0) {
+                    return None;
+                }
+                f_over / r_min
+            }
+        };
+        worst = worst.max(off);
+    }
+    Some(worst)
+}
+
 /// `ANALYTIC_TOL_SCALE` (`SmConfig.h`).
 const ANALYTIC_TOL_SCALE: f64 = 100.0;
 
@@ -366,6 +638,35 @@ impl<O: ClassObject> ProtoTopologyManager<O> {
             }
         }
 
+        // Not upstream (#3062): where both faces are planes, cylinders, cones or spheres the
+        // continued span is certified on both surfaces over its whole length, or the edge is
+        // refused, naming how far
+        let mut surfaces = Vec::new();
+        for which in [1, 2] {
+            let surface = d
+                .src_topo(which)
+                .filter(|_| d.src_topo_type(which) == PointClassType::Face)
+                .filter(|o| o.class_kind() == ClassObjectKind::Face)
+                .and_then(|f| brep.face_surface(f))
+                .flatten()
+                .and_then(|s| brep.surface(s));
+            surfaces.extend(surface);
+        }
+        if let [s1, s2] = surfaces[..] {
+            if let Some(bound) = span_off_surfaces(&new_3d, end_param, new_param, [s1, s2]) {
+                if bound.is_nan() || bound > tol {
+                    return Err(SolidError::Failed(Box::new(Failure::new(
+                        OP,
+                        Reason::EdgeOffItsFaces {
+                            gap,
+                            bound,
+                            tolerance: tol,
+                        },
+                    ))));
+                }
+            }
+        }
+
         let mut ivl = interval;
         let set = if at_end {
             ivl.set_max(new_param)
@@ -684,6 +985,33 @@ mod tests {
     #![allow(clippy::unwrap_used, clippy::indexing_slicing)] // test helpers
     use super::*;
 
+    /// #3062: the certified bound of a span off two planes is its distance from them, to
+    /// rounding, for a line; for a quarter circle at least that, and at most that over cos² 45°.
+    #[test]
+    fn a_span_is_bounded_off_its_surfaces() {
+        use crate::surface::plane::p1::Plane;
+        use crate::surface::surface::Surface;
+        let plane =
+            |z: f64, n: Vec3| Surface::from(Plane::from_normal(Point3::new(0.0, 0.0, z), n));
+        let (z, y) = (Vec3::new(0.0, 0.0, 1.0), Vec3::new(0.0, 1.0, 0.0));
+        let end = Point3::new(4.0, 0.0, 0.0);
+        let line = Curve::from(BSplineCurve::create_line_segment(3, Point3::ZERO, end).unwrap());
+        let ivl = line.natural_interval();
+        let (a, b) = (ivl.evaluate(0.25), ivl.evaluate(0.75));
+        let on = span_off_surfaces(&line, a, b, [&plane(0.0, z), &plane(0.0, y)]).unwrap();
+        assert!(on <= 1e-14, "{on:e}");
+        let off = span_off_surfaces(&line, a, b, [&plane(1e-3, z), &plane(0.0, y)]).unwrap();
+        assert!((off - 1e-3).abs() <= 1e-15, "{off:e}");
+        let arc = crate::primitive::api_curves::create_arc(Point3::ZERO, 2.0, 0.0, 90.0).unwrap();
+        let arc = Curve::from(arc);
+        let ivl = arc.natural_interval();
+        let (a, b) = (ivl.evaluate(0.1), ivl.evaluate(0.6));
+        let on = span_off_surfaces(&arc, a, b, [&plane(0.0, z), &plane(0.0, z)]).unwrap();
+        assert!(on <= 1e-14, "{on:e}");
+        let off = span_off_surfaces(&arc, a, b, [&plane(-2e-4, z), &plane(0.0, z)]).unwrap();
+        assert!((2e-4..=4e-4 * (1.0 + 1e-12)).contains(&off), "{off:e}");
+    }
+
     /// A clamped cubic with three spans (inner knots 1 and 2) over `[0, 3]`.
     fn cubic() -> Curve {
         let pts = [