diff --git a/crates/solidcore/src/blend/feature_executive.rs b/crates/solidcore/src/blend/feature_executive.rs index 7299c5bd5b..ff48a0c91e 100644 --- a/crates/solidcore/src/blend/feature_executive.rs +++ b/crates/solidcore/src/blend/feature_executive.rs @@ -96,7 +96,6 @@ use crate::geom::a_object::BooleanOperation; use crate::geom::context::{with_mark_lock, MarkType, ALL_MARKS}; use crate::geom::extent1d::Extent1d; use crate::geom::extent3d::Extent3d; -use crate::geom::geom_utility::p2::ContinuityType; use crate::geom::global_solver::SolverOperation; use crate::geom::tol::{self, SmallTopoType}; use crate::intersect::point_class::{PointClassType, PointClassification}; @@ -1624,7 +1623,13 @@ fn ext_surface_rebuild( .clone() .ok_or_else(|| SolidError::null_pointer(OP))?; let new_surf = surface - .create_extended_surface(k, ext, ContinuityType::G1) + .extended_sides( + k, + &crate::surface::side::BoxSide::ALL, + crate::surface::extend::Continuity::Tangent, + ext, + ) + .map(|e| e.surface) .map_err(|e| { if e.is_unported_or_cancelled() { e diff --git a/crates/solidcore/src/blend/fillet_corner/p2.rs b/crates/solidcore/src/blend/fillet_corner/p2.rs index 73f8c453a0..771f4da9d0 100644 --- a/crates/solidcore/src/blend/fillet_corner/p2.rs +++ b/crates/solidcore/src/blend/fillet_corner/p2.rs @@ -109,7 +109,6 @@ use crate::blend::fillet_standard_solver::p1::{FilletLaw, VariableRadiusFs}; use crate::curve::curve::p1::GlobalPointSolveOptions; use crate::curve::curve::{Curve, CurveKind}; use crate::error::{ignored, EntityRef, Result, SolidError, SM_ERR_INVALID_INPUT}; -use crate::geom::geom_utility::p2::ContinuityType; use crate::geom::global_solver::{SolutionRequested, SolverOperation}; use crate::geom::solution_array::SolutionObject; use crate::geom::vector2d::Point2; @@ -433,8 +432,14 @@ impl FilletCorner { .surface .as_ref() .ok_or_else(|| null_deref("side face without surface"))?; - let new_surface = - side_surface.create_extended_surface(k, 10.0 * gap_dist, ContinuityType::G1)?; + let new_surface = side_surface + .extended_sides( + k, + &crate::surface::side::BoxSide::ALL, + crate::surface::extend::Continuity::Tangent, + 10.0 * gap_dist, + )? + .surface; // Copy pNewSurface, when possible as an analytic surface. let new_analytic_surface = new_surface.copy_and_add_analytics(k)?; diff --git a/crates/solidcore/src/blend/fillet_corner/topology/mixed_3x2.rs b/crates/solidcore/src/blend/fillet_corner/topology/mixed_3x2.rs index 9cf8a81bad..be7270977b 100644 --- a/crates/solidcore/src/blend/fillet_corner/topology/mixed_3x2.rs +++ b/crates/solidcore/src/blend/fillet_corner/topology/mixed_3x2.rs @@ -37,7 +37,6 @@ use crate::blend::fillet_geom::p1::{ }; use crate::blend::fillet_solver::p1::FilletSolver; use crate::error::{EntityRef, Result, SolidError}; -use crate::geom::geom_utility::p2::ContinuityType; use crate::intersect::point_class::PointClassType; use crate::kernel::Kernel; use crate::nurbs::b_spline_surface_n_lib::p2::OrientType; @@ -361,7 +360,14 @@ impl FilletCorner { .surface .as_ref() .ok_or_else(|| SolidError::null_pointer(OP))?; - let new_surface = surface.create_extended_surface(k, ext_dist, ContinuityType::G1)?; + let new_surface = surface + .extended_sides( + k, + &crate::surface::side::BoxSide::ALL, + crate::surface::extend::Continuity::Tangent, + ext_dist, + )? + .surface; // Copy pNewSurface, when possible as an analytic surface let new_analytic_surface = new_surface.copy_and_add_analytics(k)?; diff --git a/crates/solidcore/src/blend/fillet_corner/topology/nx1_end.rs b/crates/solidcore/src/blend/fillet_corner/topology/nx1_end.rs index 971d73c07d..b2e0c67285 100644 --- a/crates/solidcore/src/blend/fillet_corner/topology/nx1_end.rs +++ b/crates/solidcore/src/blend/fillet_corner/topology/nx1_end.rs @@ -28,7 +28,6 @@ use crate::blend::fillet_geom::p1::{ FilletBrep, FilletCornerId, FilletEdgeType, FilletGeomId, FilletVertex, FilletVertexType, }; use crate::error::{EntityRef, Result, SolidError, SM_ERR, SM_ERR_INVALID_INPUT}; -use crate::geom::geom_utility::p2::ContinuityType; use crate::intersect::point_class::PointClassType; use crate::kernel::Kernel; use crate::math::{Point3, EFF_ZERO}; @@ -268,11 +267,13 @@ impl FilletCorner { // there depends on how its curvature differs from the rail's face across the // seam, which a linear (G1) extension would not keep (tangency-labs/smlib@856b54e). let extension = if mixed && (tangent_side[0] || tangent_side[1]) { - ContinuityType::CInfinity + crate::surface::extend::Continuity::Natural } else { - ContinuityType::G1 + crate::surface::extend::Continuity::Tangent }; - let new_surface = surface.create_extended_surface(k, ext_dist, extension)?; + let new_surface = surface + .extended_sides(k, &crate::surface::side::BoxSide::ALL, extension, ext_dist)? + .surface; // Copy pNewSurface, when possible as an analytic surface let new_extend_surface = new_surface.copy_and_add_analytics(k)?; // If AddExtendedOriginalSurface() deleted pNewExtendSurface, then diff --git a/crates/solidcore/src/blend/fillet_section/certificate.rs b/crates/solidcore/src/blend/fillet_section/certificate.rs index a95837c4ce..ae2bbeb229 100644 --- a/crates/solidcore/src/blend/fillet_section/certificate.rs +++ b/crates/solidcore/src/blend/fillet_section/certificate.rs @@ -1426,10 +1426,15 @@ pub struct GnVarying { pub r: [f64; 2], /// `pCells`: I's enclosures, computed once per cell. pub cells: GnAreaCells, + /// The thumbweight's proved enclosure, `[lo, hi]`, over every span certified with this + /// ([`gn_certify_span`]: the union of its pieces' [`GnRowTerms::w`] at each level whose + /// bound it returns); `[inf, -inf]` before any. Not upstream. + pub w: [f64; 2], } /// `sm_GnRowTerms`: the largest `|(sum_ik dg_ikj c'_ik)|` (`x`) and `|W_j|` (`y`), -/// the least `W_j` (`w_lo`) per row, and the largest `|S - O|` (`r`). +/// the least `W_j` (`w_lo`) per row, the largest `|S - O|` (`r`), and the equal area +/// thumbweight's proved enclosure over the pieces (`w`). #[derive(Clone, Debug, PartialEq)] pub struct GnRowTerms { /// `X`. @@ -1440,6 +1445,10 @@ pub struct GnRowTerms { pub w_lo: Vec, /// `R`. pub r: f64, + /// The equal area thumbweight over every piece folded in, `[lo, hi]`: the union of the + /// interval Newton roots `w0` (each holds every piece's root, for every s in it); `[inf, + /// -inf]` before any. Not upstream: the fork discards `w0` after the inversion. + pub w: [f64; 2], } /// `sm_Centred`: the coefficients `0..=SM_GN_TAY` of a series over a piece in their @@ -1548,6 +1557,8 @@ pub fn gn_row_terms_over( if !root { return false; } + let [w_lo, w_hi] = rt.w; + rt.w = [cmin(w_lo, w0.lo), cmax(w_hi, w0.hi)]; let over_w0 = vary.cells.over(k, w0); let sw = gn_invert(&tg, &over_w0, w0); let form = gn_form_normalized(n, sw); @@ -1906,6 +1917,7 @@ pub fn gn_certify_span( y: vec![0.0; ld], w_lo: vec![f64::INFINITY; ld], r: 0.0, + w: [f64::INFINITY, f64::NEG_INFINITY], }; for pair in ends.windows(2) { if !proved { @@ -2053,7 +2065,15 @@ pub fn gn_certify_span( continue; } // certified by each parameter's own term, never larger than the fork's - bound_out = cmin(bound_out, round_up(level_bound(&uv4) * (1.0 + 1.0e-12))); + let bound = round_up(level_bound(&uv4) * (1.0 + 1.0e-12)); + if bound < bound_out { + // the thumbweight over the span, proved at the level whose bound is returned + if let Some(v) = vary.as_deref_mut() { + let ([lo, hi], [rlo, rhi]) = (v.w, rows.w); + v.w = [cmin(lo, rlo), cmax(hi, rhi)]; + } + } + bound_out = cmin(bound_out, bound); if bound_out <= target { break; } diff --git a/crates/solidcore/src/blend/fillet_section/recertify.rs b/crates/solidcore/src/blend/fillet_section/recertify.rs index 4d41d556dc..0ea9f41519 100644 --- a/crates/solidcore/src/blend/fillet_section/recertify.rs +++ b/crates/solidcore/src/blend/fillet_section/recertify.rs @@ -453,16 +453,35 @@ fn pieces(pb: &ProceduralBlend, ld: usize) -> Result> { /// net's distance from its rows and the floor follow the round's size, not where it lies; /// the move's own rounding ([`frame_rounding`]) is added to the net's distance. pub fn recertify(k: &Kernel, pb: &ProceduralBlend) -> Result { + Ok(recertify_with_thumbweight(k, pb)?.0) +} + +/// [`recertify`], with the equal area thumbweight's proved enclosure `[lo, hi]` over every +/// span recertified ([`GnVarying::w`](super::certificate::GnVarying::w), from the certificate's interval Newton roots); `None` +/// for any other section (a constant thumbweight is its own). +/// +/// # Errors +/// +/// Those of [`recertify`]. +pub fn recertify_with_thumbweight( + k: &Kernel, + pb: &ProceduralBlend, +) -> Result<(f64, Option<[f64; 2]>)> { crate::profile::profile_scope!(k, "blend::certificate::recertify"); let centre = frame_centre(pb.net())?; let local = pb.translated(k, centre * -1.0)?; recertify_in_frame(k, &local, pb.bound(), frame_rounding(pb.net())?) } -/// [`recertify`] of a blend already in its own frame, against `bound`, its net's -/// distance from its rows taken as at least `shift` more. +/// [`recertify_with_thumbweight`] of a blend already in its own frame, against `bound`, +/// its net's distance from its rows taken as at least `shift` more. #[allow(clippy::too_many_lines)] -fn recertify_in_frame(k: &Kernel, pb: &ProceduralBlend, bound: f64, shift: f64) -> Result { +fn recertify_in_frame( + k: &Kernel, + pb: &ProceduralBlend, + bound: f64, + shift: f64, +) -> Result<(f64, Option<[f64; 2]>)> { let d = pb.definition(); let n = d.order as usize; let equal_area = n >= 2 && d.thumbweight == GN_THUMBWEIGHT_EQUAL_AREA; @@ -636,7 +655,8 @@ fn recertify_in_frame(k: &Kernel, pb: &ProceduralBlend, bound: f64, shift: f64) } achieved = achieved.max(total); } - Ok(achieved) + let [lo, hi] = vary.w; + Ok((achieved, (equal_area && lo <= hi).then_some([lo, hi]))) } /// The stations of a definition at `sigma`, its seeds' knots, where `needed`: each solved @@ -759,20 +779,18 @@ pub(crate) struct DefinitionSpan { pub(crate) stations: Vec, } -/// The span of a definition's net between its last two stations (`high`) or its first two, -/// solved and certified as [`net_of_definition`] does each. +/// The span of a definition's net from its station `first` to the next, solved and +/// certified as [`net_of_definition`] does each. /// /// # Errors /// /// Those of [`net_of_definition`]. -pub(crate) fn end_span_of_definition( +pub(crate) fn span_of_definition( k: &Kernel, pb: &ProceduralBlend, - high: bool, + first: usize, target: f64, ) -> Result { - let count = pb.definition().seeds[0].knots(None)?.0.len(); - let first = if high { count.saturating_sub(2) } else { 0 }; let (_, mut spans, certified, ld, stations) = spans_of_definition(k, pb, target, Some(first))?; let (homog, weights) = spans.pop().ok_or_else(|| SolidError::sm_err(OP))?; Ok(DefinitionSpan { diff --git a/crates/solidcore/src/blend/fillet_section/section.rs b/crates/solidcore/src/blend/fillet_section/section.rs index a309fe5952..15b63236e3 100644 --- a/crates/solidcore/src/blend/fillet_section/section.rs +++ b/crates/solidcore/src/blend/fillet_section/section.rs @@ -458,6 +458,7 @@ pub(crate) fn varying( n, r, cells: GnAreaCells::default(), + w: [f64::INFINITY, f64::NEG_INFINITY], }; if equal_area { let mut w_lo = f64::INFINITY; diff --git a/crates/solidcore/src/blend/fillet_solver/p2.rs b/crates/solidcore/src/blend/fillet_solver/p2.rs index 6468c00680..88d9e04881 100644 --- a/crates/solidcore/src/blend/fillet_solver/p2.rs +++ b/crates/solidcore/src/blend/fillet_solver/p2.rs @@ -117,7 +117,6 @@ use crate::error::{ignore_status, Result, SolidError}; use crate::geom::axis2_placement::Axis2Placement; use crate::geom::extent1d::{Extent1d, ExtentNd}; use crate::geom::extent2d::Extent2d; -use crate::geom::geom_utility::p2::ContinuityType; use crate::geom::global_solver::{SolutionRequested, SolverOperation}; use crate::geom::local_solve_nd::{BoundaryHandler, EvalNFunctions, LocalSolveNd}; use crate::geom::matrix::Matrix; @@ -2593,7 +2592,14 @@ impl FilletSolver { adj_surface } else { let adj = adj_surface.ok_or_else(nul)?; - let new_surface = adj.create_extended_surface(k, extension_dist, ContinuityType::G1)?; + let new_surface = adj + .extended_sides( + k, + &crate::surface::side::BoxSide::ALL, + crate::surface::extend::Continuity::Tangent, + extension_dist, + )? + .surface; // Copy pNewSurface, when possible as an analytic surface let copy = new_surface.copy_and_add_analytics(k)?; pseudo diff --git a/crates/solidcore/src/offset/construction/sheets.rs b/crates/solidcore/src/offset/construction/sheets.rs index 93c5ef2491..19983a60e4 100644 --- a/crates/solidcore/src/offset/construction/sheets.rs +++ b/crates/solidcore/src/offset/construction/sheets.rs @@ -53,7 +53,7 @@ use crate::error::{Result, SolidError}; use crate::kernel::Kernel; use crate::math::Vec3; use crate::nurbs::knots::{KnotVector, Side}; -use crate::nurbs::math::{CPoint, NOW, NOZ}; +use crate::nurbs::math::{gamma, CPoint, NOW, NOZ}; use crate::nurbs::nurbs_srf::decompose_srf; use crate::nurbs::srf_basic::Surface as Nurb; @@ -339,23 +339,38 @@ pub struct CertifiedOffset { /// The knots of `G` in one direction: the base's breakpoints `breaks` (an interior one /// of multiplicity `m` in the base's degree `deg` taken `3 - min(deg - m - 1, 2)` times, -/// so `G` keeps the normal's continuity) and the inserted `extra` knots (simple). -fn offset_knots(breaks: &[f64], mults: &[usize], deg: usize, extra: &[f64]) -> Result> { +/// so `G` keeps the normal's continuity) and the inserted `extra` knots (simple). An +/// interior breakpoint of multiplicity `deg` (the net only `C0` in its parameter) is a +/// crease unless `smooth` says its joint was proved tangent continuous ([`joint_meets`]): +/// then the normal is `C0` there and `G` takes it `3` times. +fn offset_knots( + breaks: &[f64], + mults: &[usize], + smooth: &[bool], + deg: usize, + extra: &[f64], +) -> Result> { let (Some(&lo), Some(&hi)) = (breaks.first(), breaks.last()) else { return Err(invalid("no breakpoints")); }; let mut interior: Vec<(f64, usize)> = Vec::new(); - for (&x, &m) in breaks + for ((&x, &m), at) in breaks .iter() .skip(1) .take(breaks.len().saturating_sub(2)) .zip(mults) + .zip(0..) { - // the normal is C^(deg - m - 1) there - let Some(cont) = deg.checked_sub(m + 1) else { - return Err(invalid( - "a crease inside the face: its offset is not one surface", - )); + // the normal is C^(deg - m - 1) there; at multiplicity deg C0 where the joint + // was proved tangent continuous + let cont = match deg.checked_sub(m + 1) { + Some(c) => c, + None if m == deg && smooth.get(at) == Some(&true) => 0, + None => { + return Err(invalid( + "a crease inside the face: its offset is not one surface", + )) + } }; interior.push((x, DEGREE - cont.min(DEGREE - 1))); } @@ -369,6 +384,146 @@ fn offset_knots(breaks: &[f64], mults: &[usize], deg: usize, extra: &[f64]) -> R Ok(k) } +/// The most pieces [`joint_meets`] splits one joint's edge into. +const MAX_JOINT_PIECES: usize = 256; + +/// Whether the joint of `ps` at the interior breakpoint `at` of the u (`along_u`) or v +/// breakpoints, a knot of multiplicity the base's degree (where its net is only `C0` in the +/// parameter), is tangent continuous closely enough that the offsets of its two sides meet +/// there within half the offset's tolerance `tol` at the distance `d`: then its offset is +/// one surface, `C0` across the joint, and [`certified_offset`] certifies it as any other +/// breakpoint. A knot's multiplicity is not the joint's geometry: an exact extension's joint +/// (Natural, Tangent, Ruled) is of multiplicity `p` and tangent continuous. +/// +/// The comparison is of the two one-sided normals along the joint's edge, `n_L` and `n_R` +/// (the patches' reduced normals [`reduced_normal`], of the one-sided derivatives across +/// the joint and the common one along it), as Bernstein polynomials of the edge's parameter: +/// `D = n_L . n_R`, `C = n_L x n_R` and `Y = |n_L|^2 |n_R|^2`. Where every coefficient of +/// `D` is positive the angle `theta` between the unit normals is below a right angle, and +/// where `|C| <= s sqrt(Y)` on the piece (the norm of `C`'s components' largest coefficients +/// against `Y`'s smallest, each with its rounding) `sin theta <= s`, so the offset points' +/// gap `|d| |n_L - n_R| = 2 |d| sin(theta / 2) <= sqrt 2 |d| sin theta`; `s = tol / (2 sqrt +/// 2 |d|)` holds the gap to `tol / 2`. An edge piece not proved so is halved (an explicit +/// stack, at most [`MAX_JOINT_PIECES`] pieces); a piece whose end shows the jump (`D` not +/// positive, or `|C| > s sqrt(Y)` there) ends the test: a crease. +fn joint_meets( + ps: &[Vec], + along_u: bool, + at: usize, + d: f64, + tol: f64, + b: &Binomials, +) -> Result { + let before = at + .checked_sub(1) + .ok_or_else(|| invalid("a joint at the domain's end"))?; + let s = tol / (2.0 * core::f64::consts::SQRT_2 * d.abs()); + if !(s.is_finite() && positive(s)) { + return Ok(false); + } + let pairs: Vec<(&Patch, &Patch)> = if along_u { + let (Some(left), Some(right)) = (ps.get(before), ps.get(at)) else { + return Err(invalid("no patch row at the joint")); + }; + left.iter().zip(right.iter()).collect() + } else { + ps.iter() + .map(|row| match (row.get(before), row.get(at)) { + (Some(l), Some(r)) => Ok((l, r)), + _ => Err(invalid("no patch column at the joint")), + }) + .collect::>()? + }; + for (left, right) in pairs { + let nl = reduced_normal(left, b)?; + let nr = reduced_normal(right, b)?; + // the edge: the left patch's high end, the right patch's low end, in the joint's + // direction (a polynomial constant in it, of the edge's parameter in the other) + let (el, er) = if along_u { + ( + restrict3(&nl, 1.0, 1.0, 0.0, 1.0), + restrict3(&nr, 0.0, 0.0, 0.0, 1.0), + ) + } else { + ( + restrict3(&nl, 0.0, 1.0, 1.0, 1.0), + restrict3(&nr, 0.0, 1.0, 0.0, 0.0), + ) + }; + let mut stack = vec![(el, er)]; + let mut pieces = 1_usize; + while let Some((a, c)) = stack.pop() { + let dd = dot(&a, &c, b)?; + let x = cross(&a, &c, b)?; + let y = dot(&a, &a, b)?.mul(&dot(&c, &c, b)?, b)?; + // every coefficient of these forms is a weighted mean of products of the normals' + // coefficients (the Bernstein product's weights sum to 1): at most 2 M for a + // component of the cross product, 3 M for D and 9 M^2 for Y, with M the product of + // the two normals' largest coefficients; each is formed in a few dozen operations + // and the splits are convex combinations: gamma_64 of that bounds their rounding + let m = a + .iter() + .map(Bernstein2::max_abs_coefficient) + .fold(0.0_f64, f64::max) + * c.iter() + .map(Bernstein2::max_abs_coefficient) + .fold(0.0_f64, f64::max); + let (rc, rd, ry) = ( + gamma(64) * 2.0 * m, + gamma(64) * 3.0 * m, + gamma(64) * 9.0 * m * m, + ); + // |n_L x n_R| over the piece: at most the norm of the components' largest + // coefficients (each component within the hull of its coefficients) + let cmax = x + .iter() + .map(|f| f.max_abs_coefficient().powi(2)) + .sum::() + .sqrt(); + let floor = |v: f64| (v - ry).max(0.0).sqrt(); + if dd.min_coefficient() > rd && cmax + rc <= s * floor(y.min_coefficient()) { + continue; + } + // a jump at an end of the piece (its end values are coefficients): a crease, no + // split proves it smooth + let first = |f: &Bernstein2| f.coefficients().first().copied().unwrap_or(0.0); + let last = |f: &Bernstein2| f.coefficients().last().copied().unwrap_or(0.0); + let norm = |pick: &dyn Fn(&Bernstein2) -> f64| { + x.iter().map(|f| pick(f).powi(2)).sum::().sqrt() + }; + let creased = |dv: f64, cv: f64, yv: f64| dv <= -rd || cv - rc > s * (yv + ry).sqrt(); + if creased(first(&dd), norm(&first), first(&y)) + || creased(last(&dd), norm(&last), last(&y)) + { + return Ok(false); + } + pieces = pieces.saturating_add(1); + if pieces > MAX_JOINT_PIECES { + return Ok(false); + } + let halves = |f: &Bernstein3| -> (Bernstein3, Bernstein3) { + let [x, yy, zz] = f; + let sp = |g: &Bernstein2| { + if along_u { + g.split_v(0.5) + } else { + g.split_u(0.5) + } + }; + let (x0, x1) = sp(x); + let (y0, y1) = sp(yy); + let (z0, z1) = sp(zz); + ([x0, y0, z0], [x1, y1, z1]) + }; + let (a0, a1) = halves(&a); + let (c0, c1) = halves(&c); + stack.push((a0, c0)); + stack.push((a1, c1)); + } + } + Ok(true) +} + /// The Greville abscissae of the knots `k` of degree [`DEGREE`]. fn greville(k: &[f64]) -> Vec { k.windows(DEGREE + 2) @@ -622,12 +777,28 @@ pub fn certified_offset(k: &Kernel, s: &Nurb, d: f64, tol: f64) -> Result Result> { + mults + .iter() + .zip(1..) + .map(|(&m, at)| { + if m == deg { + joint_meets(&ps, along_u, at, d, tol, &b) + } else { + Ok(false) + } + }) + .collect() + }; + let smooth_u = smooth(&sub_mults, p, true)?; + let smooth_v = smooth(&svb_mults, q, false)?; let mut extra_u: Vec = Vec::new(); let mut extra_v: Vec = Vec::new(); for _round in 0..MAX_ROUNDS { k.check_cancelled(OP)?; - let ku = offset_knots(&sub, &sub_mults, p, &extra_u)?; - let kv = offset_knots(&svb, &svb_mults, q, &extra_v)?; + let ku = offset_knots(&sub, &sub_mults, &smooth_u, p, &extra_u)?; + let kv = offset_knots(&svb, &svb_mults, &smooth_v, q, &extra_v)?; let gu = greville(&ku); let gv = greville(&kv); let mut grid: Vec> = Vec::new(); @@ -1043,6 +1214,100 @@ mod tests { worst } + /// `s` (a B-spline net) extended `how` past its u max by `e`: its net, the old end now an + /// interior knot of multiplicity `p` (the joint). + fn extended_net(s: &Nurb, how: crate::surface::extend::Continuity, e: f64) -> Nurb { + use crate::nurbs::b_spline_surface::BSplineSurface; + use crate::surface::side::BoxSide; + use crate::surface::surface::{Surface, SurfaceKind}; + let k = Kernel::default(); + let src = Surface::from(BSplineSurface::from_gw_nurb(s)); + let dom = src.natural_uv_domain_with_nurb(&k).unwrap(); + let out = src + .extended_by_distance(&k, &dom, BoxSide::UMax, how, e) + .unwrap() + .surface; + let SurfaceKind::BSpline(net) = out.kind() else { + panic!("not a B-spline: {out:?}"); + }; + net.gw_nurb().unwrap().clone() + } + + #[test] + fn a_tangent_continuous_extension_joint_at_multiplicity_p_is_no_crease() { + use crate::surface::extend::Continuity; + let s = parabola(false); + for how in [Continuity::Natural, Continuity::Tangent, Continuity::Ruled] { + let g = extended_net(&s, how, 6.0); + let (p, _) = g.degrees(); + // the joint: the old end, u = 1, now inside at multiplicity p + let at = g.u_knots().iter().filter(|&&x| x == 1.0).count(); + assert_eq!(at, p, "{how:?}: {:?}", g.u_knots()); + for d in [2.0, -2.0] { + let tol = 1e-7; + let made = certified_offset(&Kernel::default(), &g, d, tol).unwrap(); + assert!(made.bound <= tol, "{how:?} d {d}: bound {}", made.bound); + let err = sampled_error(&g, &made.surface, d); + assert!( + err <= made.bound * 1.000_001 + 1e-12, + "{how:?} d {d}: {err}" + ); + } + } + } + + #[test] + fn a_kink_at_multiplicity_p_is_a_crease() { + // two quadratic spans meeting at (0, 10) at multiplicity 2, their tangents there + // (10, 5) and (10, -5): a real crease + let pts = [ + (-10.0, 0.0), + (-5.0, 7.5), + (0.0, 10.0), + (5.0, 7.5), + (10.0, 0.0), + ]; + let net = pts + .iter() + .map(|&(x, y)| vec![CPoint::new(x, y, 0.0, NOW), CPoint::new(x, y, 20.0, NOW)]) + .collect(); + let s = Nurb::new( + net, + 2, + 1, + vec![0., 0., 0., 0.5, 0.5, 1., 1., 1.], + vec![0., 0., 1., 1.], + ) + .unwrap(); + let err = certified_offset(&Kernel::default(), &s, 2.0, 1e-7).unwrap_err(); + assert!( + matches!(err, SolidError::InvalidInput { reason, .. } if reason.starts_with("a crease")), + "{err:?}" + ); + // the same joint made smooth (its neighbours on the tangent through it) is none + let pts = [ + (-10.0, 0.0), + (-5.0, 10.0), + (0.0, 10.0), + (5.0, 10.0), + (10.0, 0.0), + ]; + let net = pts + .iter() + .map(|&(x, y)| vec![CPoint::new(x, y, 0.0, NOW), CPoint::new(x, y, 20.0, NOW)]) + .collect(); + let s = Nurb::new( + net, + 2, + 1, + vec![0., 0., 0., 0.5, 0.5, 1., 1., 1.], + vec![0., 0., 1., 1.], + ) + .unwrap(); + let made = certified_offset(&Kernel::default(), &s, 2.0, 1e-7).unwrap(); + assert!(made.bound <= 1e-7); + } + #[test] fn radius_classes() { assert_eq!(classify_radius(1.0, 1e-6), SheetClass::Regular); diff --git a/crates/solidcore/src/offset/offset_executive/p1.rs b/crates/solidcore/src/offset/offset_executive/p1.rs index ee15f5ac85..56a042d634 100644 --- a/crates/solidcore/src/offset/offset_executive/p1.rs +++ b/crates/solidcore/src/offset/offset_executive/p1.rs @@ -83,7 +83,6 @@ use crate::geom::context::{Context, MarkType}; use crate::geom::extent1d::Extent1d; use crate::geom::extent2d::Extent2d; use crate::geom::extent3d::Extent3d; -use crate::geom::geom_utility::p2::ContinuityType; use crate::geom::global_solver::{SolutionRequested, SolverOperation}; use crate::geom::object::{NotifyOperation, ObjectRef}; use crate::geom::solution_array::{Solution, SolutionArray}; @@ -93,7 +92,7 @@ use crate::heal::stitch::p2 as stitch; use crate::intersect::topology_solver::p1::TopologySolver; use crate::intersect::topology_solver::p2::TopologySolveOptions; use crate::kernel::Kernel; -use crate::math::{min2, Vec3, BIG_DOUBLE, EFF_ZERO, EFF_ZERO_SQRT, PI}; +use crate::math::{Vec3, BIG_DOUBLE, EFF_ZERO, EFF_ZERO_SQRT, PI}; use crate::nurbs::b_spline_surface_n_lib::p1::OffsetParamType; use crate::nurbs::b_spline_surface_n_lib::p2::OrientType; use crate::nurbs::knots::Side; @@ -101,6 +100,8 @@ use crate::offset::offset_geometry_creation::{OffsetGeometryCreation, SsiData}; use super::edgeuse::{CapNotifications, MarkSource}; use super::p2::OffsetNotifications; +use crate::surface::extend::Continuity; +use crate::surface::side::BoxSide; use crate::surface::surface::p4::NoObject; use crate::surface::surface::{Surface, SurfaceKind}; use crate::topology::brep::p1::{self as bp1, SetRegionVoidType}; @@ -184,222 +185,6 @@ pub struct VertexOffset { pub curve_orientations: Vec, } -/// Not upstream (WR-532; the fork's tangency-labs/smlib@b283ea2, `sm_ExtendConeAsItself`): -/// cone (or cylinder) `surface` extended by `dist` along its rulings and around its axis -/// as itself, each direction's end spans carried on (`SM_CT_CINFINITY`), which stays on -/// the circle where a tangent extension (`SM_CT_G1R`) leaves it. Around the axis it goes -/// no farther than half the turn the cone leaves open on each side, so it never laps -/// itself; a closed cone is extended along its rulings only. Around the axis each -/// side is extended on its own and each old end made a knot of multiplicity equal to -/// the degree again ([`full_knot`]; tangency-labs/smlib@3626d4b). -/// -/// # Errors -/// -/// Those of the extensions and knot insertions. -fn extend_cone_as_itself( - kernel: &Kernel, - surface: &Surface, - cone: &crate::surface::cone::Cone, - dist: f64, -) -> Result { - use crate::nurbs::b_spline_surface::SurfParamType; - let (around, along) = if cone.rev().swap_uv() { - (SurfParamType::V, SurfParamType::U) - } else { - (SurfParamType::U, SurfParamType::V) - }; - let along_srf = - surface.create_extended_surface_dir(kernel, along, dist, ContinuityType::CInfinity)?; - - // around: half the open turn on each side, in arc length on the smaller end circle - let sweep_rad = cone.rev().step_uv_domain().u_interval().length() * PI / 180.0; - let radius = min2(cone.bot_radius(), cone.top_radius()); - let open = 2.0 * PI - sweep_rad; - let around_dist = min2(dist, 0.5 * 0.5 * open * radius); - // (`SM_CAST_PTR(SmBSplineSurface, ..)`: the B-spline surface and its subclasses) - let is_bspline = !matches!( - along_srf.kind(), - SurfaceKind::OffsetSurface(_) - | SurfaceKind::CurveBounded(_) - | SurfaceKind::Step(_) - | SurfaceKind::SrfInVolume(_) - ); - if !is_bspline || open <= 0.0 || around_dist <= EFF_ZERO_SQRT * (1.0 + radius) { - return Ok(along_srf); - } - // (tangency-labs/smlib@3626d4b) each side on its own, each old end then a knot of - // multiplicity equal to the degree again: SMLib's own form for a circle - // (`SmCone::IsNurbSurfaceCone` takes only a triple knot at each end and double knots - // inside); both sides in one call extend the first side's extended span again and - // leave a single knot at the old end, and a face on that surface is no longer found - // to be a cone - let around_ivl = if around == SurfParamType::U { - along_srf.natural_uv_domain_with_nurb(kernel)?.u_interval() - } else { - along_srf.natural_uv_domain_with_nurb(kernel)?.v_interval() - }; - let (min_side, max_side) = if around == SurfParamType::U { - (SurfParamType::UMin, SurfParamType::UMax) - } else { - (SurfParamType::VMin, SurfParamType::VMax) - }; - let mut min_srf = along_srf.create_extended_surface_dir( - kernel, - min_side, - around_dist, - ContinuityType::CInfinity, - )?; - full_knot(&mut min_srf, around, around_ivl.min())?; - let mut both = min_srf.create_extended_surface_dir( - kernel, - max_side, - around_dist, - ContinuityType::CInfinity, - )?; - full_knot(&mut both, around, around_ivl.max())?; - Ok(both) -} - -/// Not upstream (WR-728, #2484): ball patch `patch` (the face's surface `surface`) -/// extended by `dist` as itself, a ball patch on the same sphere. The Ball's patch is the -/// sphere's triangle with one corner, its pole, on a side of the patch collapsed to a point -/// (a singularity, as the extension finds it); its iso-curves through the pole are great -/// circles, and the third side is a great circle's arc. The patch is carried on with -/// `SM_CT_CINFINITY`, the polynomial continuation of its homogeneous net, so the identity -/// `|X - C W|^2 = R^2 W^2` that puts it on the sphere holds on the extension too, where a -/// tangent extension (`SM_CT_G1R`) leaves the sphere at once. Around the pole (across -/// its two sides through the pole) each side goes no farther than a quarter of the turn -/// the third side's great circle leaves open, as a cone's turn is bounded -/// ([`extend_cone_as_itself`]); away from the pole (across the third side) no farther -/// than half the arc left from the third side's farthest point to the pole's antipode, -/// where the great circles through the pole meet again. The pole's side is not extended. -/// A patch with no singular side (not the Ball's triangle) is extended as its net is, a -/// plain B-spline. -/// -/// # Errors -/// -/// Those of the singularity tests, evaluations and extensions, and -/// [`SolidError::InvalidInput`] when the extended net has a weight that is not positive. -fn extend_ball_patch_as_itself( - kernel: &Kernel, - surface: &Surface, - patch: &crate::surface::ball_patch::BallPatch, - dist: f64, -) -> Result { - use crate::math::Point2; - use crate::nurbs::b_spline_surface::SurfParamType; - let net = patch.net(); - let dom = net.natural_uv_domain(); - let (lo, hi) = (dom.min(), dom.max()); - // the pole's side as the extension finds it (`CreateExtendedSurface`'s singularity - // test at the domain's corners): the pole, the third side's ends, the direction round - // the pole and the side across the third side - let at_lo = surface.is_singularity(kernel, lo, EFF_ZERO, false)?; - let at_hi = surface.is_singularity(kernel, hi, EFF_ZERO, false)?; - let (pole, ends, around, away) = if at_lo.singular && at_lo.direction == SurfParamType::U { - ( - lo, - [Point2::new(lo.x, hi.y), hi], - SurfParamType::U, - SurfParamType::VMax, - ) - } else if at_lo.singular { - ( - lo, - [Point2::new(hi.x, lo.y), hi], - SurfParamType::V, - SurfParamType::UMax, - ) - } else if at_hi.singular && at_hi.direction == SurfParamType::U { - ( - hi, - [lo, Point2::new(hi.x, lo.y)], - SurfParamType::U, - SurfParamType::VMin, - ) - } else if at_hi.singular { - ( - hi, - [lo, Point2::new(lo.x, hi.y)], - SurfParamType::V, - SurfParamType::UMin, - ) - } else { - return Ok(Surface::from(net.create_extended_surface( - kernel, - dist, - ContinuityType::G1R, - )?)); - }; - let centre = patch.centre(); - let radius = patch.radius(); - let [e0, e1] = ends; - let p = net.evaluate_point(pole)? - centre; - let a = net.evaluate_point(e0)? - centre; - let b = net.evaluate_point(e1)? - centre; - let angle = |x: Vec3, y: Vec3| x.cross(y).length().atan2(x.dot(y)); - // the third side's farthest point from the pole: the foot of the pole's antipode on - // the great circle's plane, where it is on the arc, else an end - let m = a.cross(b); - let foot = m * m.dot(p) - p * m.dot(m); - let mut farthest = angle(p, a).max(angle(p, b)); - if foot != Vec3::new(0.0, 0.0, 0.0) - && a.cross(foot).dot(m) >= 0.0 - && foot.cross(b).dot(m) >= 0.0 - { - farthest = farthest.max(angle(p, foot)); - } - let around_dist = min2(dist, 0.25 * (2.0 * PI - angle(a, b)) * radius); - let away_dist = min2(dist, 0.5 * (PI - farthest) * radius); - let round = - net.create_extended_surface_dir(kernel, around, around_dist, ContinuityType::CInfinity)?; - let both = - round.create_extended_surface_dir(kernel, away, away_dist, ContinuityType::CInfinity)?; - Ok(Surface::from( - crate::surface::ball_patch::BallPatch::create(&both, centre, radius)?, - )) -} - -/// Not upstream (WR-532; tangency-labs/smlib@3626d4b): makes the knot of B-spline -/// `surface` at `at` in direction `dir` of multiplicity equal to the degree -/// (`InsertOneKnot` of the missing count; a knot within `EFF_ZERO` of the interval's -/// length counts as the same). -/// -/// # Errors -/// -/// `SM_ERR_NULL_POINTER` for a surface that is not a B-spline surface (upstream's -/// `NER` of the cast); those of the knot query and insertion. -fn full_knot( - surface: &mut Surface, - dir: crate::nurbs::b_spline_surface::SurfParamType, - at: f64, -) -> Result<()> { - use crate::nurbs::b_spline_surface::SurfParamType; - const OP: &str = "offset::offset_executive::p1::extend_cone_as_itself"; - let SurfaceKind::BSpline(bs) = surface.kind_mut() else { - return Err(SolidError::null_pointer(OP)); - }; - let (knots, mults) = bs.knots(dir, None)?; - let domain = bs.natural_uv_domain(); - let ivl = if dir == SurfParamType::U { - domain.u_interval() - } else { - domain.v_interval() - }; - let same = EFF_ZERO * (1.0 + ivl.length()); - let mut have = 0; - for (&knot, &mult) in knots.iter().zip(&mults) { - if (knot - at).abs() <= same { - have = mult; - } - } - let degree = bs.degree(dir); - if have < degree { - bs.insert_one_knot(at, degree - have, dir)?; - } - Ok(()) -} - /// The 27-bit object namespace under which the lamina clean-up's `DeleteFaces` /// (merge p4's arena adapter) records its notifications; distinct from the Boolean /// kernels' namespaces, which share no attribute store with an offset executive. @@ -836,8 +621,17 @@ impl OffsetExecutive { Ok(out) } - /// Returns a cached extension, creating it with distance `offset_distance * 5.0` - /// and G1-ruled continuity. A hit never re-extends or replaces the cached surface. + /// Returns a cached extension, creating it with distance `offset_distance * 5.0`, every + /// side exactly ruled ([`Surface::extended_sides_keeping_parameters`], [`Continuity::Ruled`]): a plane, cone, + /// sphere, torus or ball patch as itself, a procedural blend by its definition, an + /// offset surface as the offset of its base extended. A hit never re-extends or replaces + /// the cached surface. + /// + /// Not upstream: upstream's `CreateExtendedSurface(.., SM_CT_G1R, ..)` fits a ruled + /// extension and removes knots within a tolerance (the WR-532 / WR-728 opt-ins + /// `offset_extends_cones_as_themselves` and `offset_extends_ball_patches_as_themselves` + /// extended cones and ball patches as themselves; the exact path always does, so the + /// flags no longer change this). /// /// `base_key` represents upstream pointer identity, including across Rust clones. /// Distinct source identities must use distinct keys. Failed extensions are not cached. @@ -853,19 +647,14 @@ impl OffsetExecutive { // The obsolete offset-specific `else` at cpp:2158 is unreachable: // NER(pBase) already established pOrigSurf == pBase != NULL. let dist = self.offset_distance * 5.0; - let extension = match base.kind() { - // Not upstream (WR-532), opt-in `offset_extends_cones_as_themselves` - SurfaceKind::Cone(cone) if kernel.fixes.offset_extends_cones_as_themselves => { - extend_cone_as_itself(kernel, base, cone, dist)? - } - // Not upstream (WR-728), opt-in `offset_extends_ball_patches_as_themselves` - SurfaceKind::BallPatch(patch) - if kernel.fixes.offset_extends_ball_patches_as_themselves => - { - extend_ball_patch_as_itself(kernel, base, patch, dist)? - } - _ => base.create_extended_surface(kernel, dist, ContinuityType::G1R)?, - }; + let extension = base + .extended_sides_keeping_parameters( + kernel, + &BoxSide::ALL, + Continuity::Ruled, + dist, + )? + .surface; Ok(entry.insert(extension)) } } @@ -3979,8 +3768,14 @@ mod tests { let mut executive = executive(); let base = plane(0.0); let expected = base - .create_extended_surface(&kernel, 0.125 * 5.0, ContinuityType::G1R) - .unwrap(); + .extended_sides_keeping_parameters( + &kernel, + &BoxSide::ALL, + Continuity::Ruled, + 0.125 * 5.0, + ) + .unwrap() + .surface; let first = executive .get_create_extended_surface(&kernel, 41, &base) .unwrap(); diff --git a/crates/solidcore/src/surface/extend.rs b/crates/solidcore/src/surface/extend.rs index 3029059997..e7f0fa9c7f 100644 --- a/crates/solidcore/src/surface/extend.rs +++ b/crates/solidcore/src/surface/extend.rs @@ -789,6 +789,12 @@ struct Local { /// A bound on how far each extension point's homogeneous coordinates are from the exact /// construction's (0 for the source's own points, which are kept bit for bit). err: f64, + /// The degree across the rows (the other direction's): `q + 1` consecutive rows make a + /// point of the surface. + q: usize, + /// The 3D bound of a construction made in a local frame ([`Self::natural`]), which + /// [`Self::bound`] returns in place of the one from [`Self::err`]. + frame: Option, } fn op_err() -> SolidError { @@ -806,6 +812,7 @@ impl Local { [c.x, c.y, c.z, c.w] } }; + let q = if side.is_u() { net.q() } else { net.p() }; let (rows, knots, p): (Vec>, &[f64], usize) = if side.is_u() { let cols = cpts.first().map_or(0, Vec::len); let rows = (0..cols) @@ -822,6 +829,8 @@ impl Local { p, rational, err: 0.0, + q, + frame: None, }; if !side.is_high() { l.reverse(); @@ -956,6 +965,8 @@ impl Local { p: self.p, rational: self.rational, err, + q: self.q, + frame: None, }) } @@ -977,6 +988,9 @@ impl Local { /// from the homogeneous bound [`Self::err`]: `S = sum N H / sum N w`, so a change of at /// most `err` in each coordinate moves it at most `err (sqrt 3 + |S|) / w_min`. fn bound(&self) -> f64 { + if let Some(b) = self.frame { + return b; + } if self.err == 0.0 { return 0.0; } @@ -1056,23 +1070,64 @@ impl Local { p: self.p, rational: self.rational, err: self.err, + q: self.q, + frame: None, }) } /// Natural: the last span's Bezier points (from a copy with its start knot at multiplicity /// `p`) continued past `b` by de Casteljau at `t = 1 + e / (b - c)`: the points /// `blossom(1^(p-i), t^i)` are the polynomial's Bezier points over `[b, b + e]`. + /// + /// In a local frame: each row is translated by its own origin `O_j` (its end point), `H' + /// = H - w O_j` (the weights unchanged), continued, and its new points translated back, + /// `H = H' + w O_j`; the source's rows are kept as they are. Translation commutes with the + /// continuation (both are linear in the homogeneous points and keep the weights), so the + /// construction is the same; only its rounding changes. The continuation amplifies its + /// input's and its own rounding by `a^p` (`a = |1 - t| + |t|`), which is now relative to + /// the rows' local size, not to their absolute coordinates; the translations round once + /// each, the way back not amplified. The 3D bound ([`Self::bound`]) is then that of the + /// local construction ([`Self::frame`]). fn natural(&self, e: f64) -> Result { let b = self.end()?; let c = self.before_end(0)?; let p = self.p; - let span = if c > self.start()? { - self.with_knot(c)? + // each row's origin: its end point (any point would do: it is exact by choice) + let origins: Vec<[f64; 3]> = self + .rows + .iter() + .map(|r| { + let x = r.last().ok_or_else(op_err)?; + let w = if self.rational { x[3] } else { 1.0 }; + if !(w.is_finite() && w > 0.0) { + return Err(SolidError::invalid_input( + OP, + "a weight that is not positive", + )); + } + Ok([x[0] / w, x[1] / w, x[2] / w]) + }) + .collect::>()?; + // the rows translated: one product and one difference a coordinate + let mut local = self.clone(); + let mut moved = 0.0_f64; + for (r, o) in local.rows.iter_mut().zip(&origins) { + for x in r.iter_mut() { + let w = x[3]; + for (v, &oc) in x.iter_mut().zip(o) { + moved = moved.max(v.abs() + (w * oc).abs()); + *v -= w * oc; + } + } + } + local.err = self.err + gamma(2) * moved; + let span = if c > local.start()? { + local.with_knot(c)? } else { - self.clone() + local.clone() }; let t = 1.0 + e / (b - c); - let ext: Vec> = span + let ext_local: Vec> = span .rows .iter() .map(|r| { @@ -1091,15 +1146,63 @@ impl Local { out }) .collect(); - let knots = vec![b + e; p + 1]; // de Casteljau past the span: each level `(1 - t) x + t y` scales the errors and the // points by `a = |1 - t| + |t|` and rounds by gamma_3 of that, over p levels (the // rounding of t itself only reparametrises the continuation, its points stay on it) let a = (1.0 - t).abs() + t.abs(); #[allow(clippy::cast_possible_truncation, clippy::cast_possible_wrap)] let amp = a.powi(p as i32); + let err_local = amp * (span.err + gamma(3 * p + 3) * span.size()); + // the local frame's extent: its points' largest l1 size over their weight, the + // smallest weight (the source's rows and the new ones), and the new rows translated + // back, each coordinate one product and one sum + let (mut wmin, mut pmax, mut back) = (f64::INFINITY, 0.0_f64, 0.0_f64); + for x in local + .rows + .iter() + .flatten() + .chain(ext_local.iter().flatten()) + { + let w = if self.rational { x[3] } else { 1.0 }; + wmin = wmin.min(w); + pmax = pmax.max((x[0].abs() + x[1].abs() + x[2].abs()) / w.abs()); + } + let ext: Vec> = ext_local + .into_iter() + .zip(&origins) + .map(|(r, o)| { + r.into_iter() + .map(|mut x| { + let w = x[3]; + for (v, &oc) in x.iter_mut().zip(o) { + back = back.max(v.abs() + (w * oc).abs()); + *v += w * oc; + } + x + }) + .collect() + }) + .collect(); + // rows of different origins mix at a point: the q + 1 rows of a span across, so the + // largest l1 distance between two origins within q + 1 consecutive rows + let mut diam = 0.0_f64; + for (i, oi) in origins.iter().enumerate() { + for oj in origins.iter().skip(i + 1).take(self.q) { + let dist = (oi[0] - oj[0]).abs() + (oi[1] - oj[1]).abs() + (oi[2] - oj[2]).abs(); + diam = diam.max(dist); + } + } + let knots = vec![b + e; p + 1]; let mut out = self.joined(ext, &knots)?; - out.err = amp * (span.err + gamma(3 * p + 3) * span.size()); + out.err = err_local; + // `S = sum N (H' + w O_j) / sum N w`: an error of at most `err_local` in each local + // coordinate and weight moves it at most `err_local (sqrt 3 + |S - O_j| ...)`, with + // `|S - O_j| <= pmax + diam` for every row j at the point; the way back's rounding + // `gamma_2 back` in each coordinate moves it at most `sqrt 3` of that over the weight + out.frame = Some( + (err_local * (2.0 + pmax + diam) * (1.0 + gamma(4)) + 2.0 * gamma(2) * back) / wmin + * (1.0 + gamma(8)), + ); Ok(out) } @@ -1331,6 +1434,8 @@ fn keeps(net: &Net, out: &Net, side: BoxSide) -> Result { Ok(points_kept && a.len() <= b.len()) } +mod analytic; +mod ball; mod blend; mod sides; pub use sides::{ExtendedSides, SideExtension}; diff --git a/crates/solidcore/src/surface/extend/analytic.rs b/crates/solidcore/src/surface/extend/analytic.rs new file mode 100644 index 0000000000..698a31b944 --- /dev/null +++ b/crates/solidcore/src/surface/extend/analytic.rs @@ -0,0 +1,596 @@ +//! A plane, cone, sphere or torus extended past the sides of its box by its own NURBS net +//! continued Natural ([`Surface::extended_sides`] on an analytic kind), still on its analytic +//! surface. +//! +//! Not upstream: upstream extends an analytic's NURB as any B-spline's. Here (owner, +//! 2026-11-04, the ball patch's rule) the analytic's OWN net is continued Natural past each +//! side (the homogeneous polynomial of its end span carried on), so the source's net, knots +//! and parameters stay bit for bit (a straight direction's end points to rounding, below): a face keeps its UV trims on the result (the shell's +//! extended offset replaces the face's surface with UV trim creation off). Widening the +//! analytic itself ([`Surface::continued`]) rebuilds its canonical form over the wider STEP +//! box, which moves every kept UV trim. The implicit equation composed with the net is a +//! polynomial identity, so it holds on its continuation: the result is certified on the +//! analytic on its net ([`lies_on`]), within its rounding. It is a B-spline, carrying the +//! source's kept exact round where it has one (#2487); an offset surface of an analytic is +//! the offset of its base so extended. +//! +//! A plane follows the same rule for a caller keeping its trims (owner, 2026-11-04; +//! [`Surface::extended_sides_keeping_parameters`]): its bilinear net continued, a B-spline +//! whose parameters over the source's box are the plane's own bit for bit, where widening the +//! plane ([`Surface::continued`]) handed back a plane kind the shell's offset takes other +//! branches for. Its form is straight both ways, so every continuity is its continuation. +//! +//! **Straight directions.** Across a side where the form is of degree 1 (a plane either way, +//! a cone along its rulings) the continuation is the end span's own line, so its joint knot +//! is removed exactly and the span is widened (`joint_removed`): the result is again one +//! straight span there, as the kind's canonical form is and its recognisers read it (a cone's +//! form must have its two end circles only), on the same parameters. The source's end control +//! points are then the widened span's points at the old ends, to their rounding (the plane's +//! corners, a ruling's ends); its other points, its knots and its parameters stay bit for +//! bit. +//! +//! A side is left where the analytic closes in its direction (a full turn), where the side +//! is not proved longer than a point (a pole or an apex), or a plane's side at +//! `+-INFINITE_PARAMETER`, as [`Surface::extended_sides`] leaves them; where no side is +//! extended the result is the source itself. Every other side asked is extended `distance` across it, capped (said so, +//! [`SideExtension::UpToLimit`]): +//! - **round the axis**, no farther than a quarter of the turn the source leaves open (the +//! shell's rule for a round, WR-532), so it never laps itself; +//! - **along the generator**, never past the cone's apex or the sphere's pole, and never past +//! the torus's full turn of its tube. +//! +//! **The distance** across a side is measured at the side's middle on the iso across it, in +//! closed form on the analytic: the iso is certified (on its homogeneous control points, within +//! their rounding) in a plane normal to the axis (round it: a parallel, a circle about the +//! axis), in a plane through the axis (a sphere's meridian, a great circle; a torus's tube, the +//! circle of the minor radius about its centre), or straight (a cone's ruling or a plane's +//! iso, of degree 1). +//! Its length is the circle's radius times the angle swept from the side through the +//! extension's middle point to its end, or the ruling's chord. The extension parameter is +//! solved for the distance (or the cap) by bisection on that closed form, met to `1e-12` of it. +//! +//! **Weights.** An extension whose weights are not all positive is refused +//! ([`ExtendRefusal::WeightVanishes`]) naming the longest distance found whose weights are. +//! +//! Deviations: none (no upstream counterpart). No recursion. + +use super::ball::{knot_box, reach}; +use super::sides::{round_the_axis, ExtendedSides, SideExtension}; +use super::{keeps, Continuity, ExtendRefusal, Local}; +use crate::error::{Result, SolidError}; +use crate::geom::vector2d::Point2; +use crate::geom::vector3d::Vec3; +use crate::kernel::Kernel; +use crate::math::INFINITE_PARAMETER; +use crate::nurbs::b_spline_surface::{BSplineSurface, SurfParamType}; +use crate::nurbs::math::{gamma, NOW}; +use crate::nurbs::srf_basic::{Direction, Surface as Net}; +use crate::nurbs::srf_tool::p1::extract_iso_crv; +use crate::surface::implicit::{lies_on, AnalyticSupport, Implicit}; +use crate::surface::offset_surface::{OffsetSurface, SS_UNKNOWN}; +use crate::surface::side::{Across, BoxSide, Bracket}; +use crate::surface::surface::{Surface, SurfaceKind}; +use std::f64::consts::PI; + +const OP: &str = "surface::extend::analytic::analytic_extended_sides"; + +/// A full turn, radians. +const TURN: f64 = 2.0 * PI; + +fn refused(refusal: ExtendRefusal, distance: Option) -> SolidError { + SolidError::ExtendRefused { + operation: OP, + refusal, + distance, + } +} + +/// The iso across a side, by the analytic's closed form. +#[derive(Debug, Clone, Copy)] +enum Iso { + /// A circle about the axis, in a plane normal to it. + Parallel, + /// A sphere's great circle through its poles, of the sphere's radius. + Meridian { radius: f64 }, + /// A torus's tube circle, of its minor radius about the point at its major radius. + Tube { major: f64, minor: f64 }, + /// A cone's straight ruling. + Ruling, +} + +/// The analytic's axis and the data its caps need. +#[derive(Debug, Clone, Copy)] +struct Frame { + /// A point of the axis (the sphere's or torus's centre). + origin: Vec3, + /// The unit axis. + axis: Vec3, + /// A cone's apex, if it has one (not a cylinder). + apex: Option, + /// The points' absolute size, for the rounding. + size: f64, + /// A quarter of the turn the source leaves open round the axis, radians. + quarter: f64, + /// The turn a torus's tube leaves open, radians (0 for the other kinds). + tube_left: f64, +} + +/// `s`, an analytic kind (cone, sphere, torus; a plane with `planes`) or an offset surface of one, extended past +/// the sides `asked` of its natural box by `distance` (module docs); `None` for any other +/// kind. +pub(super) fn analytic_extended_sides( + k: &Kernel, + s: &Surface, + asked: &[BoxSide], + distance: f64, + planes: bool, +) -> Result> { + let mut offsets = Vec::new(); + let mut carrier = s; + while let SurfaceKind::OffsetSurface(o) = carrier.kind() { + offsets.push(o.offset_distance()); + carrier = o.base_surface(); + } + let plane = matches!(carrier.kind(), SurfaceKind::Plane(_)); + if !(plane && planes + || matches!( + carrier.kind(), + SurfaceKind::Cone(_) | SurfaceKind::Sphere(_) | SurfaceKind::Torus(_) + )) + { + return Ok(None); + } + if !(distance.is_finite() && distance > 0.0) { + return Err(SolidError::invalid_input( + OP, + "the distance is not finite and positive", + )); + } + // the form proved on its label once; its continuation is certified at the end + let support = AnalyticSupport::check(k, carrier)?; + let implicit = support.implicit(); + let form = crate::blend::fillet_section::station::as_b_spline(k, carrier)? + .ok_or_else(|| SolidError::sm_err(OP))? + .into_owned(); + let step = carrier.step_uv_domain(); + let open = TURN - step.u_interval().length().to_radians(); + let (frame, along) = match implicit { + Implicit::Cone { + origin, + axis, + radius, + tan_angle, + } => ( + Frame { + origin, + axis, + apex: (tan_angle != 0.0).then(|| origin + axis * (-radius / tan_angle)), + size: 0.0, + quarter: 0.25 * open, + tube_left: 0.0, + }, + Iso::Ruling, + ), + Implicit::Sphere { centre, radius } => { + let AnalyticSupport::Sphere { placement, .. } = support else { + return Err(SolidError::sm_err(OP)); + }; + ( + Frame { + origin: centre, + axis: placement.x_axis().cross(placement.y_axis()), + apex: None, + size: radius, + quarter: 0.25 * open, + tube_left: 0.0, + }, + Iso::Meridian { radius }, + ) + } + Implicit::Torus { + centre, + axis, + major, + minor, + } => ( + Frame { + origin: centre, + axis, + apex: None, + size: major + minor, + quarter: 0.25 * open, + tube_left: TURN - step.v_interval().length().to_radians(), + }, + Iso::Tube { major, minor }, + ), + // straight both ways (its form is bilinear), never capped + Implicit::Plane { origin, normal } => ( + Frame { + origin, + axis: normal, + apex: None, + size: 0.0, + quarter: 0.0, + tube_left: 0.0, + }, + Iso::Ruling, + ), + }; + // the size of the net's points from the origin, for the rounding of the closed forms + let net0 = form.gw_nurb().ok_or_else(|| SolidError::sm_err(OP))?; + let mut size = frame.size; + for c in net0.cpts().iter().flatten() { + let w = if c.w == NOW { 1.0 } else { c.w }; + size = size.max(c.x.abs().max(c.y.abs()).max(c.z.abs()) / w.abs()); + } + let o = frame.origin; + let frame = Frame { + size: size + o.x.abs().max(o.y.abs()).max(o.z.abs()), + ..frame + }; + let dom_src = carrier.natural_uv_domain_with_nurb(k)?; + let mut cur = form; + let mut sides = [SideExtension::NotAsked; 4]; + let mut extension = 0.0_f64; + for (&side, done) in BoxSide::ALL.iter().zip(sides.iter_mut()) { + if !asked.contains(&side) { + continue; + } + let param = if side.is_u() { + SurfParamType::U + } else { + SurfParamType::V + }; + if carrier.is_closed(k, &dom_src, param, None, false)?.0 { + *done = SideExtension::Closed; + continue; + } + if carrier.side_diameter(k, &dom_src, side)?.lower <= 0.0 { + *done = SideExtension::Collapsed; + continue; + } + if matches!(implicit, Implicit::Plane { .. }) + && side.at(&dom_src).abs() >= INFINITE_PARAMETER + { + *done = SideExtension::Unbounded; + continue; + } + let iso = if round_the_axis(carrier, side) == Some(true) { + if open <= 0.0 { + *done = SideExtension::Closed; + continue; + } + Iso::Parallel + } else { + along + }; + let (grown, across, capped, bound) = extend_side(&cur, side, &frame, iso, distance)?; + cur = grown; + extension = extension.max(bound); + *done = if capped { + SideExtension::UpToLimit(across) + } else { + SideExtension::Extended(across) + }; + } + // nothing extended: the source itself, as [`Surface::extended_sides`] leaves it + if !sides + .iter() + .any(|x| matches!(x, SideExtension::Extended(_) | SideExtension::UpToLimit(_))) + { + return Ok(Some(ExtendedSides { + surface: s.copy(k)?, + sides, + extension, + })); + } + if !lies_on(&cur, &implicit)? { + return Err(refused(ExtendRefusal::NotExact, None)); + } + let mut surface = carrier.carry_round_on_continuation(Surface::from(cur)); + for &d in offsets.iter().rev() { + surface = Surface::from(OffsetSurface::new(k, d, surface, true, SS_UNKNOWN)?); + } + Ok(Some(ExtendedSides { + surface, + sides, + extension, + })) +} + +/// The unit vector along `v`, refused where it has no direction. +fn unit(v: Vec3) -> Result { + let l = v.length(); + if l.is_nan() || l <= 0.0 { + return Err(SolidError::sm_err(OP)); + } + Ok(v * (1.0 / l)) +} + +/// `x`'s part normal to the axis. +fn radial(f: &Frame, x: Vec3) -> Vec3 { + let d = x - f.origin; + d - f.axis * f.axis.dot(d) +} + +/// The angle between `a` and `b`, in radians. +fn angle(a: Vec3, b: Vec3) -> f64 { + a.cross(b).length().atan2(a.dot(b)) +} + +/// The circle an iso runs (its centre, its plane's normal and its radius) through `x0`; +/// `None` for a ruling. +fn circle(f: &Frame, iso: Iso, x0: Vec3) -> Result> { + let d = x0 - f.origin; + Ok(match iso { + Iso::Parallel => { + let centre = f.origin + f.axis * f.axis.dot(d); + Some((centre, f.axis, (x0 - centre).length())) + } + Iso::Meridian { radius } => Some((f.origin, unit(f.axis.cross(d))?, radius)), + Iso::Tube { major, minor } => { + let centre = f.origin + unit(radial(f, x0))? * major; + Some((centre, unit(f.axis.cross(d))?, minor)) + } + Iso::Ruling => None, + }) +} + +/// The length along the iso from `x0` through `xm` to `x1`, in closed form. +fn length(f: &Frame, iso: Iso, x0: Vec3, xm: Vec3, x1: Vec3) -> Result { + let Some((centre, n, radius)) = circle(f, iso, x0)? else { + return Ok((x1 - x0).length()); + }; + // the swept angle, the way round through the middle point + let n = if (xm - x0).cross(x1 - xm).dot(n) < 0.0 { + n * -1.0 + } else { + n + }; + let (u0, u1) = (x0 - centre, x1 - centre); + let mut phi = u0.cross(u1).dot(n).atan2(u0.dot(u1)); + if phi < 0.0 { + phi += TURN; + } + Ok(radius * phi) +} + +/// `cur` continued Natural past `side` by `distance` across it, capped (module docs): the +/// result, its across length, whether the cap stopped it, and its construction bound. +fn extend_side( + cur: &BSplineSurface, + side: BoxSide, + f: &Frame, + iso: Iso, + distance: f64, +) -> Result<(BSplineSurface, Across, bool, f64)> { + let net = cur.gw_nurb().ok_or_else(|| SolidError::sm_err(OP))?; + if matches!(iso, Iso::Ruling) && (if side.is_u() { net.p() } else { net.q() }) != 1 { + return Err(SolidError::Unsupported { + operation: OP, + reason: "a cone's form is not straight along its rulings", + }); + } + let source = Local::of(net, side)?; + let dom = knot_box(cur)?; + let straight = (if side.is_u() { net.p() } else { net.q() }) == 1; + let middle = side.point(&dom, 0.5); + let made = |e: f64| -> Result> { + let Some(g) = source.extended(Continuity::Natural, e)? else { + return Ok(None); + }; + if !g.is_finite() { + return Ok(None); + } + let bound = g.bound(); + let out = g.surface(net, side)?; + if !keeps(net, &out, side)? { + return Err(SolidError::port_invariant( + OP, + "the extension moved the source's net", + )); + } + let out = if straight { + joint_removed(&out, side, side.at(&dom))? + } else { + out + }; + Ok(Some((BSplineSurface::from_gw_nurb(&out), bound))) + }; + // the iso across the side at its middle, from the side (0) to the extension's end (1) + let from = side.at(&dom); + let at = |s: &BSplineSurface, t: f64| -> Point2 { + let end = knot_box(s).map_or(from, |d| side.at(&d)); + let x = from + t * (end - from); + if side.is_u() { + Point2::new(x, middle.y) + } else { + Point2::new(middle.x, x) + } + }; + let measure = |s: &BSplineSurface| -> Result { + let x0 = s.evaluate_point(at(s, 0.0))?; + let xm = s.evaluate_point(at(s, 0.5))?; + let x1 = s.evaluate_point(at(s, 1.0))?; + length(f, iso, x0, xm, x1) + }; + let value = |e: f64| -> Result> { made(e)?.map(|(s, _)| measure(&s)).transpose() }; + // the cap, from the side's middle point and the iso's other end inside the source + let x_side = cur.evaluate_point(middle)?; + let inside = { + let (lo, hi) = (dom.min(), dom.max()); + match side { + BoxSide::UMin => Point2::new(hi.x, middle.y), + BoxSide::UMax => Point2::new(lo.x, middle.y), + BoxSide::VMin => Point2::new(middle.x, hi.y), + BoxSide::VMax => Point2::new(middle.x, lo.y), + } + }; + let x_in = cur.evaluate_point(inside)?; + let cap = match iso { + Iso::Parallel => radial(f, x_side).length() * f.quarter, + Iso::Meridian { radius } => { + // toward the pole the latitude climbs to + let (h0, h_in) = (f.axis.dot(x_side - f.origin), f.axis.dot(x_in - f.origin)); + let pole = if h0 > h_in { f.axis } else { f.axis * -1.0 }; + radius * angle(x_side - f.origin, pole) + } + Iso::Tube { minor, .. } => minor * f.tube_left, + Iso::Ruling => match f.apex { + Some(apex) if (x_side - apex).length() < (x_in - apex).length() => { + (x_side - apex).length() + } + _ => f64::INFINITY, + }, + }; + let target = distance.min(cap); + let capped = target < distance; + let e = reach(&value, target, source.range()?)?; + let (grown, bound) = made(e)?.ok_or_else(|| SolidError::sm_err(OP))?; + let length = measure(&grown)?; + certify_iso(&grown, side, middle, f, iso)?; + // the closed form's rounding (the evaluations at the points' size, the products, the + // atan2), and the points' construction bound through the angle + let rounding = gamma(64) * (length + f.size) + 4.0 * bound; + let across = Across { + enclosure: Bracket { + lower: (length - rounding).max(0.0), + upper: length + rounding, + }, + truncation: 0.0, + rounding, + converged: true, + }; + Ok((grown, across, capped, bound)) +} + +/// `out`, a net of degree 1 across `side` continued Natural past it, with the joint knot +/// `joint` (the source's old end, now of multiplicity 1) removed: the span widened. Exact: +/// the continuation of a straight row is its own line, so the row at the joint lies on the +/// segment between its neighbours at the knots' ratio (on the homogeneous rows, so a rational +/// row too), and dropping it leaves the same surface on the same parameters. The analytic's +/// canonical form is one straight span (a plane's 2 x 2 net, a cone's two end circles), which +/// is what the kind's recognisers read; a joint left would make it two. +/// +/// # Errors +/// +/// `SM_ERR` where the joint is not a single interior knot of `out`; those of the net. +fn joint_removed(out: &Net, side: BoxSide, joint: f64) -> Result { + let (p, q) = out.degrees(); + let knots = if side.is_u() { + out.u_knots() + } else { + out.v_knots() + }; + let last = knots + .len() + .checked_sub(1) + .ok_or_else(|| SolidError::sm_err(OP))?; + let at: Vec = (0..knots.len()) + .filter(|&i| knots.get(i) == Some(&joint)) + .collect(); + let &[i] = at.as_slice() else { + return Err(SolidError::sm_err(OP)); + }; + // degree 1, clamped: an interior knot has two knots on each side, and the control + // point at knot index i is point i - 1 + if i < 2 || i.checked_add(2).is_none_or(|j| j > last) { + return Err(SolidError::sm_err(OP)); + } + let row = i.checked_sub(1).ok_or_else(|| SolidError::sm_err(OP))?; + let kept: Vec = knots + .iter() + .enumerate() + .filter(|&(j, _)| j != i) + .map(|(_, &x)| x) + .collect(); + let mut net = out.cpts(); + if side.is_u() { + if row >= net.len() { + return Err(SolidError::sm_err(OP)); + } + net = net + .into_iter() + .enumerate() + .filter(|&(j, _)| j != row) + .map(|(_, r)| r) + .collect(); + Net::new(net, p, q, kept, out.v_knots().to_vec()) + } else { + for r in &mut net { + if row >= r.len() { + return Err(SolidError::sm_err(OP)); + } + *r = std::mem::take(r) + .into_iter() + .enumerate() + .filter(|&(j, _)| j != row) + .map(|(_, x)| x) + .collect(); + } + Net::new(net, p, q, out.u_knots().to_vec(), kept) + } +} + +/// Whether the iso across `side` of `s` at the side's middle `middle` is the closed form's +/// curve: every homogeneous control point in a plane normal to the axis (a parallel) or +/// through it (a meridian, a tube), within its rounding; a ruling is straight by its +/// degree. Refused as not exact otherwise. +fn certify_iso( + s: &BSplineSurface, + side: BoxSide, + middle: Point2, + f: &Frame, + iso: Iso, +) -> Result<()> { + if matches!(iso, Iso::Ruling) { + return Ok(()); + } + let net = s.gw_nurb().ok_or_else(|| SolidError::sm_err(OP))?; + let (t, dir) = if side.is_u() { + (middle.y, Direction::U) + } else { + (middle.x, Direction::V) + }; + let curve = extract_iso_crv(net, t, dir)?; + let points: Vec<(Vec3, f64, f64)> = curve + .cpts() + .iter() + .map(|p| { + let w = if p.w == NOW { 1.0 } else { p.w }; + let h = Vec3::new(p.x, p.y, p.z); + // X W from the origin, and the size of its terms + let mag = p.x.abs().max(p.y.abs()).max(p.z.abs()) + f.size * w.abs(); + (h - f.origin * w, w, mag) + }) + .collect(); + let Some(&(x0, w0, _)) = points.first() else { + return Err(SolidError::sm_err(OP)); + }; + // the plane's normal and its offset from the origin, per unit weight + let (n, height) = if matches!(iso, Iso::Parallel) { + (f.axis, f.axis.dot(x0) / w0) + } else { + // through the axis and the iso: its farthest point from the axis sets it + let far = points + .iter() + .map(|&(x, w, _)| { + let d = x * (1.0 / w); + d - f.axis * f.axis.dot(d) + }) + .fold( + Vec3::ZERO, + |m, d| if d.length() > m.length() { d } else { m }, + ); + (unit(f.axis.cross(far))?, 0.0) + }; + for &(x, w, mag) in &points { + // the dot product's rounding and the plane's own (its normal and height from + // stored points), a few ulps of the terms' size each + let off = n.dot(x) - height * w; + if off.abs() > 8.0 * gamma(32) * (mag + height.abs() * w.abs()) { + return Err(refused(ExtendRefusal::NotExact, None)); + } + } + Ok(()) +} diff --git a/crates/solidcore/src/surface/extend/ball.rs b/crates/solidcore/src/surface/extend/ball.rs new file mode 100644 index 0000000000..8ba6bb20df --- /dev/null +++ b/crates/solidcore/src/surface/extend/ball.rs @@ -0,0 +1,495 @@ +//! A ball patch extended as itself ([`Surface::extended_sides`] on a [`BallPatch`]): its +//! own net continued Natural, still on its sphere. +//! +//! Not upstream: upstream extends a ball patch's net as any B-spline's (the base class), +//! which leaves the sphere. Here (owner, 2026-11-04) the patch's OWN net is continued +//! Natural past each side (the homogeneous polynomial of its end span carried on), so the +//! source's net and knots stay bit for bit (a face keeps its UV trims on it: the shell's +//! extended offset replaces the face's surface with UV trim creation off) and every point +//! stays on the sphere: the identity `|H - C w|^2 = R^2 w^2` that puts the net on the +//! sphere is a polynomial identity, so it holds on its continuation. The result is +//! certified on the sphere on its net ([`lies_on`]), within its rounding, and is the ball +//! patch on the same sphere. +//! +//! The patch is the sphere's triangle with one side collapsed to its pole `P` (proved: the +//! side's certified diameter's lower bound is 0), the two sides through the pole great +//! circles (meridians about the axis `C -> P`) and the third side a great circle's arc +//! across from the pole. The pole's side is collapsed (left). Each other side asked is +//! extended `distance` across it, capped (said so, [`SideExtension::UpToLimit`]): +//! - **round the pole** (a side through the pole), no farther than the sphere's radius +//! times a quarter of the turn the patch leaves open round its axis (the shell's rule, +//! WR-532 / WR-728: its gaps close on the rounds' ends there), and never past half that +//! turn in longitude, so the two sides never meet; +//! - **away from the pole** (the third side), no farther than half the arc left from the +//! third side's farthest point to the antipodal pole, where the meridians meet again +//! (WR-728). +//! +//! **The distance** across a side is measured at the side's middle on the iso across it, +//! in closed form: that iso lies on the sphere and (certified on its homogeneous control +//! points: every 4 x 4 determinant of three of them and a fourth zero within its rounding) +//! in a plane, so it is a circle, and its length from the side to the end is the circle's +//! radius times the angle swept (away from the pole the isos are meridians: the sphere's +//! radius times the angle along the great circle). The angle is the one from the side's +//! point through the extension's middle point to its end: the arc the iso runs. The Ball's +//! patch is the inverse stereographic image of a rational quadratic arc's sector (blend +//! corner `ball_triangle_patch`), so its isos run their circles monotonically, and Natural +//! carries on the same rational map. The extension parameter is solved for the distance +//! (or the cap) by bisection on that closed form, met to `1e-12` of it. +//! +//! **Weights.** An extension whose weights are not all positive is refused +//! ([`ExtendRefusal::WeightVanishes`]) naming the longest distance found whose weights are. +//! +//! Deviations: none (no upstream counterpart). No recursion. + +use super::sides::{ExtendedSides, SideExtension}; +use super::{keeps, Continuity, ExtendRefusal, Local}; +use crate::error::{Result, SolidError}; +use crate::geom::extent2d::Extent2d; +use crate::geom::vector2d::Point2; +use crate::geom::vector3d::Vec3; +use crate::kernel::Kernel; +use crate::nurbs::b_spline_surface::BSplineSurface; +use crate::nurbs::math::{gamma, NOW}; +use crate::nurbs::srf_basic::Direction; +use crate::nurbs::srf_tool::p1::extract_iso_crv; +use crate::surface::ball_patch::BallPatch; +use crate::surface::implicit::{lies_on, Implicit}; +use crate::surface::side::{Across, BoxSide, Bracket}; +use crate::surface::surface::Surface; +use std::f64::consts::PI; + +const OP: &str = "surface::extend::ball::ball_extended_sides"; + +/// Doublings of the extension parameter that bracket the distance. +const BRACKET_TRIES: usize = 60; +/// Halvings of the bracket: the length met to `MEET` of it, or the bracket a few ulps wide. +const BISECTIONS: usize = 200; +/// The relative accuracy the closed-form length is met to. +const MEET: f64 = 1e-12; + +/// The angle between `a` and `b`, in radians. +fn angle(a: Vec3, b: Vec3) -> f64 { + a.cross(b).length().atan2(a.dot(b)) +} + +fn refused(refusal: ExtendRefusal, distance: Option) -> SolidError { + SolidError::ExtendRefused { + operation: OP, + refusal, + distance, + } +} + +/// How far a side may go: round the pole by a length across it (the sphere's radius times +/// a quarter of the open turn) and no more than an angle of longitude (half the open +/// turn: the two sides never meet), or away from it by an arc along the meridians +/// (radians). +#[derive(Debug, Clone, Copy)] +enum Cap { + Round { length: f64, turn: f64 }, + Arc(f64), +} + +/// `patch` (the surface `s`) extended past the sides `asked` of its natural box by +/// `distance` as itself (module docs); `None` for a patch with no collapsed side (not the +/// sphere's triangle: extended as its net is). +pub(super) fn ball_extended_sides( + k: &Kernel, + s: &Surface, + patch: &BallPatch, + asked: &[BoxSide], + distance: f64, +) -> Result> { + if !(distance.is_finite() && distance > 0.0) { + return Err(SolidError::invalid_input( + OP, + "the distance is not finite and positive", + )); + } + // the net's knot box: the extension is past the net's ends, whatever box the patch keeps + let dom = knot_box(patch.net())?; + let mut pole_side = None; + for side in BoxSide::ALL { + if s.side_diameter(k, &dom, side)?.lower <= 0.0 { + pole_side = Some(side); + break; + } + } + let Some(pole_side) = pole_side else { + return Ok(None); + }; + let away = match pole_side { + BoxSide::UMin => BoxSide::UMax, + BoxSide::UMax => BoxSide::UMin, + BoxSide::VMin => BoxSide::VMax, + BoxSide::VMax => BoxSide::VMin, + }; + let net = patch.net(); + let (c, r) = (patch.centre(), patch.radius()); + let p = net.evaluate_point(pole_side.point(&dom, 0.5))? - c; + let a = net.evaluate_point(away.point(&dom, 0.0))? - c; + let b = net.evaluate_point(away.point(&dom, 1.0))? - c; + let z = p * (1.0 / p.length()); + let (ax, bx) = (a - z * a.dot(z), b - z * b.dot(z)); + if !(ax.length() > 0.0 && bx.length() > 0.0) { + return Err(SolidError::invalid_input( + OP, + "the third side's end is at the pole or its antipode", + )); + } + // the turn round the pole the patch spans, and the third side's farthest point from it + let lambda = angle(ax, bx); + let m = a.cross(b); + let foot = m * m.dot(p) - p * m.dot(m); + let mut theta = angle(p, a).max(angle(p, b)); + if foot.length() > 0.0 && a.cross(foot).dot(m) >= 0.0 && foot.cross(b).dot(m) >= 0.0 { + theta = theta.max(angle(p, foot)); + } + if !(lambda > 0.0 && lambda < PI && theta > 0.0 && theta < PI) { + return Err(SolidError::invalid_input( + OP, + "the patch is not a triangle under a half turn of its sphere", + )); + } + let mut cur = net.copy(); + let mut sides = [SideExtension::NotAsked; 4]; + let mut extension = 0.0_f64; + for (&side, done) in BoxSide::ALL.iter().zip(sides.iter_mut()) { + if !asked.contains(&side) { + continue; + } + if side == pole_side { + *done = SideExtension::Collapsed; + continue; + } + let cap = if side == away { + Cap::Arc(0.5 * (PI - theta)) + } else { + Cap::Round { + length: r * 0.25 * (2.0 * PI - lambda), + turn: 0.5 * (2.0 * PI - lambda), + } + }; + let (grown, across, capped, bound) = extend_side(&cur, side, c, r, z, cap, distance)?; + cur = grown; + extension = extension.max(bound); + *done = if capped { + SideExtension::UpToLimit(across) + } else { + SideExtension::Extended(across) + }; + } + // every point on the sphere, certified on the net + let on = Implicit::Sphere { + centre: c, + radius: r, + }; + if !lies_on(&cur, &on)? { + return Err(refused(ExtendRefusal::NotExact, None)); + } + Ok(Some(ExtendedSides { + surface: Surface::from(BallPatch::create(&cur, c, r)?), + sides, + extension, + })) +} + +/// The box of `s`'s knots (its B-spline form's ends). +pub(super) fn knot_box(s: &BSplineSurface) -> Result { + let net = s.gw_nurb().ok_or_else(|| SolidError::sm_err(OP))?; + let (u, v) = (net.u_knots(), net.v_knots()); + match (u.first(), u.last(), v.first(), v.last()) { + (Some(&u0), Some(&u1), Some(&v0), Some(&v1)) => Ok(Extent2d::from_bounds(u0, v0, u1, v1)), + _ => Err(SolidError::sm_err(OP)), + } +} + +/// What one ask's extension measures: its length across the side and its turn round the +/// pole's axis. +#[derive(Debug, Clone, Copy)] +struct Measure { + length: f64, + turn: f64, +} + +/// `cur` continued Natural past `side` by `distance` across it, capped by `cap`: the +/// result, its across length, whether the cap stopped it, and its construction bound. +fn extend_side( + cur: &BSplineSurface, + side: BoxSide, + c: Vec3, + r: f64, + z: Vec3, + cap: Cap, + distance: f64, +) -> Result<(BSplineSurface, Across, bool, f64)> { + let net = cur.gw_nurb().ok_or_else(|| SolidError::sm_err(OP))?; + let source = Local::of(net, side)?; + let dom = knot_box(cur)?; + let middle = side.point(&dom, 0.5); + let made = |e: f64| -> Result> { + let Some(g) = source.extended(Continuity::Natural, e)? else { + return Ok(None); + }; + if !g.is_finite() { + return Ok(None); + } + let bound = g.bound(); + let out = g.surface(net, side)?; + if !keeps(net, &out, side)? { + return Err(SolidError::port_invariant( + OP, + "the extension moved the source's net", + )); + } + Ok(Some((BSplineSurface::from_gw_nurb(&out), bound))) + }; + // the iso across the side at its middle, from the side (0) to the extension's end (1) + let from_end = side.at(&dom); + let at = |s: &BSplineSurface, f: f64| -> Point2 { + let end = knot_box(s).map_or(from_end, |d| side.at(&d)); + let from = side.at(&dom); + let x = from + f * (end - from); + if side.is_u() { + Point2::new(x, middle.y) + } else { + Point2::new(middle.x, x) + } + }; + let measure = |s: &BSplineSurface| -> Result { + let x0 = s.evaluate_point(at(s, 0.0))? - c; + let xm = s.evaluate_point(at(s, 0.5))? - c; + let x1 = s.evaluate_point(at(s, 1.0))? - c; + let n = (xm - x0).cross(x1 - xm); + let len = n.length(); + if len.is_nan() || len <= 0.0 { + return Err(SolidError::sm_err(OP)); + } + let n = n * (1.0 / len); + let cc = n * n.dot(x0); + let (u0, u1) = (x0 - cc, x1 - cc); + let mut phi = u0.cross(u1).dot(n).atan2(u0.dot(u1)); + if phi < 0.0 { + phi += 2.0 * PI; + } + let (t0, t1) = (x0 - z * x0.dot(z), x1 - z * x1.dot(z)); + Ok(Measure { + length: u0.length() * phi, + turn: angle(t0, t1), + }) + }; + let value = |e: f64, turn: bool| -> Result> { + Ok(match made(e)? { + Some((s, _)) => { + let m = measure(&s)?; + Some(if turn { m.turn } else { m.length }) + } + None => None, + }) + }; + let range = source.range()?; + let (e, capped) = match cap { + Cap::Arc(arc) => { + let target = distance.min(r * arc); + ( + reach(&|e| value(e, false), target, range)?, + target < distance, + ) + } + Cap::Round { length, turn } => { + // within the length cap and the turn's (none where the turn never gets there) + let at_turn = match solve(&|e| value(e, true), turn, range)? { + Reach::At(e_turn) => Some(( + e_turn, + value(e_turn, false)?.ok_or_else(|| SolidError::sm_err(OP))?, + )), + Reach::Short(_) => None, + }; + let target = distance.min(length); + match at_turn { + Some((e_turn, l_turn)) if l_turn <= target => (e_turn, true), + _ => ( + reach(&|e| value(e, false), target, range)?, + target < distance, + ), + } + } + }; + let (grown, bound) = made(e)?.ok_or_else(|| SolidError::sm_err(OP))?; + let length = measure(&grown)?.length; + // the closed form's rounding (the evaluations at the points' size, the cross products, + // the atan2), and the points' construction bound through the angle + let size = c.length() + r; + let rounding = gamma(64) * (length + size) + 4.0 * bound; + certify_circle(&grown, side, middle, c, size)?; + let across = Across { + enclosure: Bracket { + lower: (length - rounding).max(0.0), + upper: length + rounding, + }, + truncation: 0.0, + rounding, + converged: true, + }; + Ok((grown, across, capped, bound)) +} + +/// Where `f` reaches its target ([`solve`]). +pub(super) enum Reach { + /// At this extension parameter. + At(f64), + /// Never: `f` stops growing (its curve's continuation tends to a point) short of the + /// target, at this value. + Short(f64), +} + +/// The extension parameter where `f` (increasing from 0, `None` where a weight is not +/// positive) reaches `target`: bracketed by doubling from `e0`, then bisected. A target past +/// every positive ask is refused naming the longest value found. +pub(super) fn solve(f: &dyn Fn(f64) -> Result>, target: f64, e0: f64) -> Result { + let (mut lo, mut lo_v) = (0.0_f64, 0.0_f64); + let mut hi = None; + let mut bad = None; + let mut e = e0; + for _ in 0..BRACKET_TRIES { + match f(e)? { + Some(v) if v >= target => { + hi = Some(e); + break; + } + // no longer growing: the continuation tends to its limit point + Some(v) if lo > 0.0 && v - lo_v <= MEET * v => return Ok(Reach::Short(v)), + Some(v) => { + lo = e; + lo_v = v; + e *= 2.0; + } + None => { + bad = Some(e); + break; + } + } + } + let hi = match (hi, bad) { + (Some(h), _) => h, + (None, Some(mut b)) => { + // the longest positive ask below the first whose weights are not + for _ in 0..BISECTIONS { + let m = 0.5 * (lo + b); + if !(m > lo && m < b) { + break; + } + match f(m)? { + Some(v) if v >= target => return Ok(Reach::At(bisect(f, target, lo, m)?)), + Some(v) => { + lo = m; + lo_v = v; + } + None => b = m, + } + } + return Err(refused( + ExtendRefusal::WeightVanishes, + (lo_v > 0.0).then_some(lo_v), + )); + } + (None, None) => return Ok(Reach::Short(lo_v)), + }; + Ok(Reach::At(bisect(f, target, lo, hi)?)) +} + +/// The parameter where `f` reaches `target`, or refused as not exact naming the longest +/// value reached. +pub(super) fn reach(f: &dyn Fn(f64) -> Result>, target: f64, e0: f64) -> Result { + match solve(f, target, e0)? { + Reach::At(e) => Ok(e), + Reach::Short(v) => Err(refused(ExtendRefusal::NotExact, (v > 0.0).then_some(v))), + } +} + +/// Bisection of `[lo, hi]` (`f(lo) < target <= f(hi)`) until `f` is within `MEET` of the +/// target. +fn bisect(f: &dyn Fn(f64) -> Result>, target: f64, lo: f64, hi: f64) -> Result { + let (mut lo, mut hi) = (lo, hi); + for _ in 0..BISECTIONS { + let m = 0.5 * (lo + hi); + if !(m > lo && m < hi) { + break; + } + let v = f(m)?.ok_or_else(|| SolidError::sm_err(OP))?; + if (v - target).abs() <= MEET * target { + return Ok(m); + } + if v < target { + lo = m; + } else { + hi = m; + } + } + Ok(hi) +} + +/// Whether the iso across `side` of `s` at the side's middle `middle` lies in a plane: +/// every determinant of three spread homogeneous control points (from the centre `c`) and +/// a fourth is zero within its rounding (`size` the points' absolute size), and the three +/// span a plane. Refused as not exact otherwise. +fn certify_circle( + s: &BSplineSurface, + side: BoxSide, + middle: Point2, + c: Vec3, + size: f64, +) -> Result<()> { + let net = s.gw_nurb().ok_or_else(|| SolidError::sm_err(OP))?; + // the iso along the side's across direction: a u-curve across a u side + let (t, dir) = if side.is_u() { + (middle.y, Direction::U) + } else { + (middle.x, Direction::V) + }; + let iso = extract_iso_crv(net, t, dir)?; + let mut rows: Vec<[f64; 4]> = Vec::new(); + let mut mags: Vec = Vec::new(); + for p in iso.cpts() { + let w = if p.w == NOW { 1.0 } else { p.w }; + rows.push([p.x - c.x * w, p.y - c.y * w, p.z - c.z * w, w]); + mags.push(p.x.abs().max(p.y.abs()).max(p.z.abs()) + (size + 1.0) * w.abs()); + } + let n = rows.len(); + let (Some(&r0), Some(&rm), Some(&rl)) = (rows.first(), rows.get(n / 2), rows.last()) else { + return Err(SolidError::sm_err(OP)); + }; + let (Some(&m0), Some(&mm), Some(&ml)) = (mags.first(), mags.get(n / 2), mags.last()) else { + return Err(SolidError::sm_err(OP)); + }; + let g = |row: [f64; 4], x: usize| row.get(x).copied().unwrap_or(0.0); + // the 3 x 3 minor of the three spread rows on columns i, j, l + let minor = |i: usize, j: usize, l: usize| -> f64 { + g(r0, i) * (g(rm, j) * g(rl, l) - g(rm, l) * g(rl, j)) + - g(r0, j) * (g(rm, i) * g(rl, l) - g(rm, l) * g(rl, i)) + + g(r0, l) * (g(rm, i) * g(rl, j) - g(rm, j) * g(rl, i)) + }; + let three = m0 * mm * ml; + // the three span a plane: a minor well past its rounding + let spans = [(0, 1, 2), (0, 1, 3), (0, 2, 3), (1, 2, 3)] + .iter() + .any(|&(i, j, l)| minor(i, j, l).abs() > 1e3 * 6.0 * gamma(16) * three); + if !spans { + return Err(refused(ExtendRefusal::NotExact, None)); + } + let (m123, m023, m013, m012) = ( + minor(1, 2, 3), + minor(0, 2, 3), + minor(0, 1, 3), + minor(0, 1, 2), + ); + for (&row, &mag) in rows.iter().zip(&mags) { + // expansion along the fourth row; |fl(det) - det| <= gamma_K per(|A|) <= gamma_K 4! + // times the product of the rows' sizes + let det = -g(row, 0) * m123 + g(row, 1) * m023 - g(row, 2) * m013 + g(row, 3) * m012; + if det.abs() > gamma(32) * 24.0 * three * mag { + return Err(refused(ExtendRefusal::NotExact, None)); + } + } + Ok(()) +} diff --git a/crates/solidcore/src/surface/extend/blend.rs b/crates/solidcore/src/surface/extend/blend.rs index 8a6de24353..38e2b54c04 100644 --- a/crates/solidcore/src/surface/extend/blend.rs +++ b/crates/solidcore/src/surface/extend/blend.rs @@ -25,7 +25,13 @@ //! convex hull: every weight positive, as the Natural extension proves). The result's //! bound is `max(eps, eps G + the extension's rounding)`, the growth `G` rounded up. //! -//! A G^n section (order 2 and up) has no growth factor here and is refused across. +//! A G^n section (order `n >= 2`) is the same argument with `d_S = 3 n + 1` and its weight +//! `W = A B^n C^n` (`A = (1 - t)^(n+1) + t^(n+1)`, `B = 1 - w + w t`, `C = 1 - w t`), bounded +//! in closed form over a proved range of its thumbweight `w`: a constant one, or the equal +//! area one's enclosure from its certificate's interval Newton roots (the stations' own w's +//! padded would be no proof). `A` grows with `|t - 1/2|`; `B` and `C` are bilinear in `(t, +//! w)`, their extremes at the corners of the (t range) x (w range) box. A corner past which +//! `B` or `C` is not positive is a weight that vanishes: refused. //! //! **Along the spine (s).** Past the definition's end station (the side its last or first //! station, the net on its whole section `t in [0, 1]`, s along u): the definition is @@ -33,14 +39,19 @@ //! its seeds by a span to the new station, as the construction's seeds are the feet's //! Hermite curves on the stations), and the net gains that span of the construction's: //! the stations at the end and at `s_new` solved, the span's Hermite rows -//! ([`end_span_of_definition`], the construction's `NetForm::span`) certified against the definition ([`gn_certify_span`]'s bound), and -//! joined after the source's net at multiplicity 3 (the construction's own joints), scaled -//! to its weights. The source's net is kept bit for bit; the piece's first row is the -//! source's own, so its bound against the span's Hermite rows is the joint's row -//! difference ([`net_distance`]). The result's bound: `max(eps, certificate + joint + -//! three floors + rounding)`. The distance is solved as any extension's +//! ([`span_of_definition`], the construction's `NetForm::span`) certified against the +//! definition ([`gn_certify_span`]'s bound), and joined after the source's net at +//! multiplicity 3 (the construction's own joints), scaled to its weights. A side INSIDE the +//! spine (a piece of the blend cut across it) runs the definition's spans already there +//! first: each span the extension crosses solved and certified the same way, the part of +//! it past the side (or up to `s_new`) its Hermite rows' cubic Bezier cut there (blossoms), +//! each part joined at multiplicity 3 and scaled to the row before it; then past the end +//! as above. The source's net is kept bit for bit; each part's first row is the row before +//! it, so its bound against its span's rows is the joint's row difference +//! ([`net_distance`]). The result's bound: `max(eps, the largest part's certificate + +//! joint, + three floors + rounding)`. The distance is solved as any extension's //! ([`Surface::extended_by_distance`]), measured on the net; an ask whose station does not -//! solve, whose foot leaves its support's box, or whose span is not proved, works as one +//! solve, whose foot leaves a B-spline support's box (an analytic support continued is still itself), or whose span is not proved, works as one //! whose weights are not positive: the refusal names the longest that does. //! //! [`gn_certify_span`]: crate::blend::fillet_section::certificate::gn_certify_span @@ -48,10 +59,11 @@ //! Deviations: none (no upstream counterpart). No recursion. use super::{keeps, op_err, solve, Continuity, Extended, Local, Made, H, OP, REL}; +use crate::blend::fillet_section::recertify::recertify_with_thumbweight; use crate::blend::fillet_section::recertify::{ - end_span_of_definition, net_distance, precision_floor, DefinitionSpan, + net_distance, precision_floor, span_of_definition, DefinitionSpan, }; -use crate::blend::fillet_section::section::{param, param_in}; +use crate::blend::fillet_section::section::{param, param_in, GN_THUMBWEIGHT_EQUAL_AREA}; use crate::blend::fillet_section::station::GnStation; use crate::error::{Result, SolidError}; use crate::geom::extent2d::Extent2d; @@ -61,6 +73,7 @@ use crate::nurbs::b_spline_curve::BSplineCurve; use crate::nurbs::b_spline_surface::BSplineSurface; use crate::nurbs::math::{gamma, CPoint, NOW}; use crate::nurbs::srf_basic::Surface as Net; +use crate::surface::implicit::AnalyticSupport; use crate::surface::procedural_blend::ProceduralBlend; use crate::surface::side::BoxSide; use crate::surface::surface::{Surface, SurfaceKind}; @@ -79,15 +92,12 @@ pub(super) fn blend_extended( return along(k, pb, dom, side, distance); } let order = pb.order(); + let n = order as usize; + // the section's degree in t: the chord's 1, the circle's 2, a G^n section's 3 n + 1 let section_degree = match order { 0 => 1, 1 => 2, - _ => { - return Err(SolidError::Unsupported { - operation: OP, - reason: "a G^n section has no growth factor past its source's t range here", - }) - } + _ => 3 * n + 1, }; let source_net = pb.net.gw_nurb().ok_or_else(|| SolidError::sm_err(OP))?; let grown = Surface::from(pb.net.clone()).extended_by_distance( @@ -129,14 +139,19 @@ pub(super) fn blend_extended( let m = section_degree + net_degree; let x = (2.0 * te - t0 - t1).abs() / (t1 - t0); let cheb = chebyshev(m, x); - let (w_hi, w_lo) = if order == 0 { - (1.0, 1.0) - } else { - let off = |t: f64| (t * (t - 1.0)).max(0.0); - let hi = 1.0 + 2.0 * off(t0).max(off(t1)); - let (lo_t, hi_t) = (ts.min(te), ts.max(te)); - let lo = if lo_t < 1.0 && hi_t > 0.0 { 0.5 } else { 1.0 }; - (hi, lo) + let (w_hi, w_lo) = match order { + 0 => (1.0, 1.0), + 1 => { + let off = |t: f64| (t * (t - 1.0)).max(0.0); + let hi = 1.0 + 2.0 * off(t0).max(off(t1)); + let (lo_t, hi_t) = (ts.min(te), ts.max(te)); + let lo = if lo_t < 1.0 && hi_t > 0.0 { 0.5 } else { 1.0 }; + (hi, lo) + } + _ => { + let w = thumbweight(k, pb)?; + gn_weight_range(n, [t0, t1], [ts.min(te), ts.max(te)], w)? + } }; let (v_hi, _) = weights(source_net); let (_, v_lo) = weights(out_net); @@ -168,6 +183,87 @@ pub(super) fn blend_extended( }) } +/// The thumbweight of `pb`'s G^n section over its face, `[lo, hi]`: a constant one, or the +/// equal area one's proved enclosure from its certificate's interval Newton roots +/// ([`recertify_with_thumbweight`]). +fn thumbweight(k: &Kernel, pb: &ProceduralBlend) -> Result<[f64; 2]> { + let tw = pb.definition().thumbweight; + if tw != GN_THUMBWEIGHT_EQUAL_AREA { + return Ok([tw, tw]); + } + recertify_with_thumbweight(k, pb)? + .1 + .ok_or(SolidError::Unsupported { + operation: OP, + reason: "the equal area thumbweight is not enclosed over the blend", + }) +} + +/// A G^n section's weight `W(t) = A(t) B(t)^n C(t)^n` (`A = (1 - t)^(n+1) + t^(n+1)`, +/// `B = 1 - w + w t`, `C = 1 - w t`, `w` the thumbweight in `w`): its largest value over the +/// source's t range `src` and its least over the extension's `ext`, rounded outward. `A` +/// grows with `|t - 1/2|` (`A' = (n + 1) (t^n - (1 - t)^n)`, of the sign of `t - 1/2`), so +/// its extremes on an interval are at the point farthest from and nearest to 1/2; `B` and +/// `C` are bilinear in `(t, w)`, so their extremes on a box are at its corners. A corner of +/// the extension's box where `B` or `C` is not positive is a weight that vanishes there: +/// refused ([`super::ExtendRefusal::WeightVanishes`]). +fn gn_weight_range(n: usize, src: [f64; 2], ext: [f64; 2], w: [f64; 2]) -> Result<(f64, f64)> { + let exp = i32::try_from(n).map_err(|_| SolidError::sm_err(OP))?; + let a = |t: f64| (1.0 - t).powi(exp + 1) + t.powi(exp + 1); + let far = |[x, y]: [f64; 2]| { + if (x - 0.5).abs() >= (y - 0.5).abs() { + x + } else { + y + } + }; + let near = |[x, y]: [f64; 2]| { + if x <= 0.5 && y >= 0.5 { + 0.5 + } else if (x - 0.5).abs() <= (y - 0.5).abs() { + x + } else { + y + } + }; + let corners = |[t0, t1]: [f64; 2]| { + let mut b = [0.0_f64; 4]; + let mut c = [0.0_f64; 4]; + for (i, (&t, &wt)) in [t0, t0, t1, t1] + .iter() + .zip(&[w[0], w[1], w[0], w[1]]) + .enumerate() + { + if let (Some(bi), Some(ci)) = (b.get_mut(i), c.get_mut(i)) { + *bi = 1.0 - wt + wt * t; + *ci = 1.0 - wt * t; + } + } + (b, c) + }; + let (bs, cs) = corners(src); + let mag = |x: &[f64; 4]| x.iter().fold(0.0_f64, |m, v| m.max(v.abs())); + let hi = a(far(src)) * mag(&bs).powi(exp) * mag(&cs).powi(exp); + let (be, ce) = corners(ext); + let least = |x: &[f64; 4]| x.iter().fold(f64::INFINITY, |m, &v| m.min(v)); + let (bl, cl) = (least(&be), least(&ce)); + if !(bl > 0.0 && cl > 0.0) { + return Err(SolidError::ExtendRefused { + operation: OP, + refusal: super::ExtendRefusal::WeightVanishes, + distance: None, + }); + } + let lo = a(near(ext)) * bl.powi(exp) * cl.powi(exp); + // each a product of 3 n + 2 roundings at most: up and down by gamma_(4 n + 8) + let g = gamma(4 * n + 8); + let (hi, lo) = (hi * (1.0 + g), lo * (1.0 - g)); + if !(hi.is_finite() && lo > 0.0) { + return Err(SolidError::sm_err(OP)); + } + Ok((hi, lo)) +} + /// `T_m(x)` for `x >= 1` (`|x|` taken), by its recurrence (every term positive and /// increasing there), rounded up. fn chebyshev(m: usize, x: f64) -> f64 { @@ -222,11 +318,16 @@ fn along( let s_side = a[0] + b[0] * side.at(dom); // outward from the side, s grows (`up`) or falls let up = side.is_high() == (b[0] > 0.0); - let end = if up { sigma.last() } else { sigma.first() }; - if end.copied() != Some(s_side) { + // the side on the definition's spine: at its end station, or inside it (extended by + // the definition's spans already there, then past its end) + let (Some(&first), Some(&last)) = (sigma.first(), sigma.last()) else { + return Err(SolidError::sm_err(OP)); + }; + let end = if up { last } else { first }; + if !(s_side == end || (s_side > first && s_side < last)) { return Err(SolidError::Unsupported { operation: OP, - reason: "the side is not the definition's end station", + reason: "the side is not on the definition's spine", }); } let source_net = pb.net.gw_nurb().ok_or_else(|| SolidError::sm_err(OP))?; @@ -238,34 +339,66 @@ fn along( }); } let target = pb.bound(); - // the definition and its end span for the ask e (in the net's parameter) - let piece_of = |e: f64| -> Result> { - let ds = e * b[0].abs(); + let rate = b[0].abs(); + // the definition and its spans from the side to s_new for the ask e (in the net's + // parameter), outward, each with its part `[x, y]` in s + let pieces_of = |e: f64| -> Result)>> { + let ds = e * rate; let s_new = if up { s_side + ds } else { s_side - ds }; let next = continued_definition(pb, s_new, up)?; - // a station that does not solve, a foot off its support or a span not proved: an - // ask that does not work, as a weight not positive - let span = match end_span_of_definition(k, &next, up, target) { - Ok(span) if span.certified.is_finite() && feet_on(k, &next, &span)? => span, - Ok(_) | Err(SolidError::InvalidInput { .. } | SolidError::Upstream { .. }) => { - return Ok(None) + let sig = next + .definition() + .seeds + .first() + .ok_or_else(op_err)? + .knots(None)? + .0; + let (lo, hi) = if up { (s_side, s_new) } else { (s_new, s_side) }; + let mut pieces = Vec::new(); + for (idx, (sa, sb)) in sig.iter().zip(sig.iter().skip(1)).enumerate() { + let (x, y) = (sa.max(lo), sb.min(hi)); + if y <= x { + continue; } - Err(e) => return Err(e), - }; - Ok(Some((next, span))) - }; - let made = |e: f64| -> Result> { - let Some((_, span)) = piece_of(e)? else { - return Ok(None); - }; - let (ext, piece, scale) = piece_rows(&source, &span, up)?; - // the piece's own bound waits for the met ask: here only that it is finite - if !piece_distance(&piece, scale)?.is_finite() { + // a station that does not solve, a foot off its support or a span not proved: + // an ask that does not work, as a weight not positive + let span = match span_of_definition(k, &next, idx, target) { + Ok(span) if span.certified.is_finite() && feet_on(k, &next, &span)? => span, + Ok(_) | Err(SolidError::InvalidInput { .. } | SolidError::Upstream { .. }) => { + return Ok(None) + } + Err(e) => return Err(e), + }; + pieces.push(Piece { + span, + part: [(x - sa) / (sb - sa), (y - sa) / (sb - sa)], + length: (y - x) / rate, + }); + } + if !up { + pieces.reverse(); + } + if pieces.is_empty() { return Ok(None); } - let mut grown = source.joined(ext, &[source.end()? + e; 4])?; - // the scaled points' rounding, one product each - grown.err = gamma(1) * grown.size(); + Ok(Some((next, pieces))) + }; + let made = |e: f64| -> Result> { + let Some((_, pieces)) = pieces_of(e)? else { + return Ok(None); + }; + let (ext, rows) = pieces_rows(&source, &pieces, up)?; + // the pieces' own bounds wait for the met ask: here only that they are finite + for (piece, scale) in &rows { + if !piece_distance(piece, *scale)?.is_finite() { + return Ok(None); + } + } + let grown = source.joined(ext, &joints(source.end()?, e, &pieces))?; + let mut grown = grown; + // the scaled points' rounding, one product each, and each part's subdivision (three + // levels of convex combinations, the part's ends a few more) + grown.err = gamma(16) * grown.size(); let bound = grown.bound(); let out = grown.surface(source_net, side)?; Ok(Some(( @@ -294,14 +427,26 @@ fn along( "the extension moved the source's net", )); } - // the met ask's definition again (the same solve and certificate, deterministic), its - // seeds the feet's Hermite curves on its stations, and the piece's bound against it + // the met ask's definition again (the same solves and certificates, deterministic); + // past the end station its seeds the feet's Hermite curves on its new span's + // stations; each piece's bound against its span let grown_end = Local::of(out_net, side)?.end()?; let e = grown_end - source.end()?; - let (mut next, span) = piece_of(e)?.ok_or_else(op_err)?; - hermite_seeds(&mut next, &span.stations, up)?; - let (_, rows, scale) = piece_rows(&source, &span, up)?; - let piece = span.certified + piece_distance(&rows, scale)?; + let (mut next, pieces) = pieces_of(e)?.ok_or_else(op_err)?; + let past_end = if up { + s_side + e * rate > last + } else { + s_side - e * rate < first + }; + if past_end { + let outer = pieces.last().ok_or_else(op_err)?; + hermite_seeds(&mut next, &outer.span.stations, up)?; + } + let (_, rows) = pieces_rows(&source, &pieces, up)?; + let mut piece = 0.0_f64; + for (p, (r, scale)) in pieces.iter().zip(&rows) { + piece = piece.max(p.span.certified + piece_distance(r, *scale)?); + } let floor = 3.0 * precision_floor(k, next.supports(), next.order(), net)?; let mut out = ProceduralBlend { net: net.clone(), @@ -319,6 +464,123 @@ fn along( }) } +/// One span of a definition an extension along its spine runs: the span, the part of it +/// it runs (`[x, y]` of its own parameter in `[0, 1]`, in s's direction) and that part's +/// length in the net's parameter. +struct Piece { + span: DefinitionSpan, + part: [f64; 2], + length: f64, +} + +/// The knots after the source's joint for `pieces` from `end` by `e` in all: each part's +/// end at multiplicity 3 (the construction's own joints), the last clamped (4), at exactly +/// `end + e`. +fn joints(end: f64, e: f64, pieces: &[Piece]) -> Vec { + let mut knots = Vec::new(); + let mut at = end; + let n = pieces.len(); + for (i, p) in pieces.iter().enumerate() { + if i + 1 == n { + knots.extend([end + e; 4]); + } else { + at += p.length; + knots.extend([at; 3]); + } + } + knots +} + +/// The cubic Bezier `p` over `[t0, t1]` of its parameter: its blossoms `B(t0^(3-i), t1^i)`. +fn sub_bezier(p: [H; 4], t0: f64, t1: f64) -> [H; 4] { + let blossom = |u: [f64; 3]| -> H { + let mut q = p; + for (level, &t) in u.iter().enumerate() { + for i in 0..3 - level { + if let (Some(&a), Some(&b)) = (q.get(i), q.get(i + 1)) { + if let Some(x) = q.get_mut(i) { + *x = super::lin(1.0 - t, a, t, b); + } + } + } + } + q.first().copied().unwrap_or([f64::NAN; 4]) + }; + [ + blossom([t0, t0, t0]), + blossom([t0, t0, t1]), + blossom([t0, t1, t1]), + blossom([t1, t1, t1]), + ] +} + +/// The extension's rows past the source's joint for `pieces` (outward), each part's rows +/// scaled to the weight of the row before it (the source's joint row for the first), and +/// for each piece its rows against its span's ([`piece_distance`]) with its scale. +#[allow(clippy::type_complexity)] +fn pieces_rows( + source: &Local, + pieces: &[Piece], + up: bool, +) -> Result<(Vec>, Vec<(Vec<(Vec, Vec)>, f64)>)> { + let mut ext: Vec> = vec![Vec::new(); source.rows.len()]; + let mut out = Vec::new(); + // the row before each piece, across: the source's joint row first + let mut before: Vec = source + .rows + .iter() + .map(|r| r.last().copied().ok_or_else(op_err)) + .collect::>()?; + for p in pieces { + let ld = p.span.ld; + if source.rows.len() != ld { + return Err(SolidError::port_invariant( + OP, + "the net is not the section's width", + )); + } + // the span's part at across index j, outward from the joint + let part = |j: usize| -> Result<[H; 4]> { + let mut c = [[0.0; 4]; 4]; + for (i, ci) in c.iter_mut().enumerate() { + let at = i + .checked_mul(ld) + .and_then(|x| x.checked_add(j)) + .ok_or_else(op_err)?; + let (q, w) = ( + p.span.homog.get(at).ok_or_else(op_err)?, + p.span.weights.get(at).ok_or_else(op_err)?, + ); + *ci = [q.x, q.y, q.z, *w]; + } + let mut sub = sub_bezier(c, p.part[0], p.part[1]); + if !up { + sub.reverse(); + } + Ok(sub) + }; + let (Some(b0), Ok(h0)) = (before.first(), part(0)) else { + return Err(op_err()); + }; + let scale = b0[3] / h0[0][3]; + let mut rows = Vec::new(); + for (j, (e, bj)) in ext.iter_mut().zip(before.iter_mut()).enumerate() { + let h = part(j)?; + let (mut hs, mut bs) = (vec![h[0]], vec![*bj]); + for x in h.iter().skip(1) { + let y = x.map(|c| scale * c); + e.push(y); + hs.push(*x); + bs.push(y); + } + *bj = *bs.last().ok_or_else(op_err)?; + rows.push((bs, hs)); + } + out.push((rows, scale)); + } + Ok((ext, out)) +} + /// `|B - H|` over the piece: `piece[j]` the net's rows along s at across index j (the /// source's joint row, then the scaled span rows) and the span's own rows (scaled by /// `scale`, the same surface), as [`net_distance`] bounds two nets. @@ -381,6 +643,8 @@ fn continued_curve(c: &BSplineCurve, s_new: f64, up: bool) -> Result Result<(Vec>, Vec<(Vec, Vec)>, f64)> { - let ld = span.ld; - if source.rows.len() != ld { - return Err(SolidError::port_invariant( - OP, - "the net is not the section's width", - )); - } - // the span's row i outward from the joint at across index j - let row = |i: usize, j: usize| -> Option { - let i = if up { i } else { 3_usize.checked_sub(i)? }; - let at = i.checked_mul(ld)?.checked_add(j)?; - let (p, w) = (span.homog.get(at)?, span.weights.get(at)?); - Some([p.x, p.y, p.z, *w]) - }; - let (Some(joint), Some(src)) = (row(0, 0), source.rows.first().and_then(|r| r.last())) else { - return Err(op_err()); - }; - let scale = src[3] / joint[3]; - let mut ext = Vec::new(); - let mut piece = Vec::new(); - for (j, srow) in source.rows.iter().enumerate() { - let (Some(&sj), Some(hj)) = (srow.last(), row(0, j)) else { - return Err(op_err()); - }; - let (mut q, mut h, mut b) = (Vec::new(), vec![hj], vec![sj]); - for i in 1..4 { - let x = row(i, j).ok_or_else(op_err)?; - let y = x.map(|c| scale * c); - q.push(y); - h.push(x); - b.push(y); - } - ext.push(q); - piece.push((b, h)); - } - Ok((ext, piece, scale)) -} - -/// Whether every foot of the span's stations is on its support's natural box. +/// Whether every foot of the span's stations is on its support: on its natural box, or +/// anywhere on an analytic support (plane, cone, sphere, torus). The stations are solved on +/// the support's net continued past its knots ([`crate::blend::fillet_section::station`]); +/// an analytic's net continued is still on the analytic (its implicit equation composed +/// with the homogeneous net is a polynomial identity, so it holds on the continuation, +/// wherever the weight is not zero), so a foot off its box is still on the surface. A +/// B-spline's continuation is no longer its face's surface: a foot off its box is not on it. fn feet_on(k: &Kernel, pb: &ProceduralBlend, span: &DefinitionSpan) -> Result { for (f, support) in pb.supports().into_iter().enumerate() { + if AnalyticSupport::of(support).is_some() { + let finite = span.stations.iter().all(|st| { + param(st, 2 * f, 0, 0).is_finite() && param(st, 2 * f + 1, 0, 0).is_finite() + }); + if !finite { + return Ok(false); + } + continue; + } let b = support.natural_uv_domain_with_nurb(k)?; for st in &span.stations { let (u, v) = (param(st, 2 * f, 0, 0), param(st, 2 * f + 1, 0, 0)); diff --git a/crates/solidcore/src/surface/extend/sides.rs b/crates/solidcore/src/surface/extend/sides.rs index de34512675..c70cdd22ec 100644 --- a/crates/solidcore/src/surface/extend/sides.rs +++ b/crates/solidcore/src/surface/extend/sides.rs @@ -17,18 +17,26 @@ //! Two more rules (owner, 2026-11-04): //! - **Unbounded.** A plane's side at `+-INFINITE_PARAMETER` (the kernel's infinite plane) is //! left: the plane has no side there to extend past. -//! - **Analytic kinds as themselves.** A cone, sphere or torus is extended as its own exact -//! shape whatever the continuity asked ([`Surface::continued`], Natural): along the -//! generator up to its apex, pole or full turn where the distance asked runs past it (the -//! refusal's named distance, the longest that works), and said so. A surface of -//! revolution is extended along its generator as asked, on its B-spline form row by row: -//! its rows are one generator turned and weighted, every construction here commutes with -//! both, so the result is the revolution of the generator extended. Round the axis is a -//! closed direction for all of them, so a side round the axis is left. +//! - **Analytic kinds as themselves.** A cone, sphere or torus (or an offset surface of one; +//! a plane too for a caller keeping its trims, [`Surface::extended_sides_keeping_parameters`]; +//! otherwise a plane is widened as a plane) is extended by its own NURBS net continued Natural whatever the continuity asked, the +//! source's net and parameters bit for bit, still on its analytic surface (certified): +//! round its axis up to a quarter of the turn the source leaves open, along its generator +//! never past its apex, pole or full turn, and said so where that caps the distance +//! (`analytic`). A surface of revolution is extended along its generator as asked, on its +//! B-spline form row by row: its rows are one generator turned and weighted, every +//! construction here commutes with both, so the result is the revolution of the generator +//! extended; its side round its axis is left. +//! +//! - **Ball patches as themselves.** A ball patch (the sphere's triangle with a pole) is +//! extended by its own net continued Natural (the source bit for bit, every point still on +//! the sphere, certified), round the pole up to a quarter of the turn it leaves open and +//! away from it up to half the arc left to the antipode, its pole's side collapsed +//! (`ball`). //! //! Deviations: none (no upstream counterpart). No recursion. -use super::{Continuity, ExtendRefusal, Extended}; +use super::{Continuity, Extended}; use crate::error::{Result, SolidError}; use crate::kernel::Kernel; use crate::math::INFINITE_PARAMETER; @@ -59,7 +67,7 @@ pub enum SideExtension { /// For a cone, sphere, torus or surface of revolution, whether side `side` runs round its /// axis (its NURBS u is the sweep unless swapped); `None` for any other kind. -fn round_the_axis(s: &Surface, side: BoxSide) -> Option { +pub(super) fn round_the_axis(s: &Surface, side: BoxSide) -> Option { let swap = match s.kind() { SurfaceKind::Cone(c) => c.rev().swap_uv(), SurfaceKind::Sphere(sp) => sp.base.swap_uv, @@ -99,6 +107,48 @@ impl Surface { how: Continuity, distance: f64, ) -> Result { + self.extended_sides_with(k, asked, how, distance, false) + } + + /// [`Surface::extended_sides`] for a caller that keeps the face's UV trims on the result + /// (the shell's extended offset replaces a face's surface with UV trim creation off): a + /// plane too is extended by its own net, a B-spline on the plane's own parameters (owner, + /// 2026-11-04: the cone, sphere and torus rule), where [`Surface::extended_sides`] widens + /// it as a plane (its kind kept, which a caller intersecting with it reads). + /// + /// # Errors + /// + /// As [`Surface::extended_sides`]. + pub fn extended_sides_keeping_parameters( + &self, + k: &Kernel, + asked: &[BoxSide], + how: Continuity, + distance: f64, + ) -> Result { + self.extended_sides_with(k, asked, how, distance, true) + } + + fn extended_sides_with( + &self, + k: &Kernel, + asked: &[BoxSide], + how: Continuity, + distance: f64, + keeping: bool, + ) -> Result { + // a ball patch as its sphere (`ball`) + if let SurfaceKind::BallPatch(patch) = self.kind() { + if let Some(e) = super::ball::ball_extended_sides(k, self, patch, asked, distance)? { + return Ok(e); + } + } + // a cone, sphere or torus (a plane too when keeping) by its own net (`analytic`) + if let Some(e) = + super::analytic::analytic_extended_sides(k, self, asked, distance, keeping)? + { + return Ok(e); + } let mut current: Option = None; // the side of the current result each side of the source's box is on let mut on = BoxSide::ALL; @@ -135,44 +185,24 @@ impl Surface { *done = SideExtension::Unbounded; continue; } - let revolution = round_the_axis(carrier, side); - if revolution == Some(true) { + // a revolution's side round its axis is left; along it, its B-spline form row by + // row as asked, which is again a revolution (its rows are one generator turned + // and weighted, and every construction commutes with both) + if round_the_axis(carrier, side) == Some(true) { *done = SideExtension::AroundAxis; continue; } - // an analytic label is extended as itself (Natural, its continuation); a - // revolution's B-spline form row by row as asked, which is again a revolution - // (its rows are one generator turned and weighted, and every construction - // commutes with both) - let analytic = - revolution.is_some() && !matches!(carrier.kind(), SurfaceKind::SurfOfRevolution(_)); - let how = if analytic { Continuity::Natural } else { how }; - let (grown, capped) = match surface.extended_by_distance(k, &dom, side, how, distance) { - Ok(g) => (g, false), - // as itself up to its apex, pole or turn: the longest distance that works - Err(SolidError::ExtendRefused { - refusal: - ExtendRefusal::PastApex | ExtendRefusal::PastPole | ExtendRefusal::OverATurn, - distance: Some(limit), - .. - }) if analytic && limit > 0.0 && limit < distance => ( - surface.extended_by_distance(k, &dom, side, how, limit)?, - true, - ), - Err(e) => return Err(e), - }; + let grown = surface.extended_by_distance(k, &dom, side, how, distance)?; for s in &mut on { *s = grown.side_of(*s); } extension = extension.max(grown.extension); - *done = if capped { - SideExtension::UpToLimit(grown.across) - } else { - SideExtension::Extended(grown.across) - }; + *done = SideExtension::Extended(grown.across); current = Some(grown); } let surface = match current { + // a Natural continuation of the net is still on an exact round's S (#2487) + Some(e) if how == Continuity::Natural => self.carry_round_on_continuation(e.surface), Some(e) => e.surface, None => self.copy(k)?, }; diff --git a/crates/solidcore/src/surface/extend/tests.rs b/crates/solidcore/src/surface/extend/tests.rs index 413cac0f69..7d37fcf119 100644 --- a/crates/solidcore/src/surface/extend/tests.rs +++ b/crates/solidcore/src/surface/extend/tests.rs @@ -7,7 +7,7 @@ use crate::nurbs::srf_tool::p1::extract_iso_crv; use crate::surface::cone::Cone; use crate::surface::plane::Plane; use crate::surface::procedural_blend::ProceduralBlend; -use std::f64::consts::{FRAC_1_SQRT_2, FRAC_PI_2}; +use std::f64::consts::{FRAC_1_SQRT_2, FRAC_PI_2, PI}; fn k() -> Kernel { Kernel::default() @@ -817,7 +817,8 @@ fn every_side_extends_in_turn_and_closed_or_collapsed_sides_are_left() { } #[test] -fn an_infinite_planes_sides_are_unbounded_and_a_finite_planes_extend() { +fn an_infinite_planes_sides_are_unbounded_and_a_finite_planes_extend_its_own_net() { + use crate::surface::implicit::{lies_on, AnalyticSupport}; let kk = k(); let inf = Surface::from(Plane::from_normal(Vec3::new(0.0, 0.0, 3.0), Vec3::Z)); let e = inf @@ -837,9 +838,22 @@ fn an_infinite_planes_sides_are_unbounded_and_a_finite_planes_extend() { Vec2::new(1.0, 1.0), &Extent2d::from_bounds(0.0, 0.0, 10.0, 10.0), )); + // widened as a plane let e = p .extended_sides(&kk, &BoxSide::ALL, Continuity::Tangent, 4.33) .unwrap(); + assert!(matches!(e.surface.kind(), SurfaceKind::Plane(_))); + let kept = inf + .extended_sides_keeping_parameters(&kk, &BoxSide::ALL, Continuity::Ruled, 4.33) + .unwrap(); + assert!(kept + .sides + .iter() + .all(|x| matches!(x, SideExtension::Unbounded))); + // keeping its parameters + let e = p + .extended_sides_keeping_parameters(&kk, &BoxSide::ALL, Continuity::Tangent, 4.33) + .unwrap(); assert!( e.sides .iter() @@ -847,11 +861,45 @@ fn an_infinite_planes_sides_are_unbounded_and_a_finite_planes_extend() { "{:?}", e.sides ); - assert!(matches!(e.surface.kind(), SurfaceKind::Plane(_))); + // its own net continued: a B-spline on the plane's own parameters, on the plane + let SurfaceKind::BSpline(grown) = e.surface.kind() else { + panic!("not a B-spline: {:?}", e.surface.kind()); + }; + let on = AnalyticSupport::check(&kk, &p).unwrap().implicit(); + assert!(lies_on(grown, &on).unwrap()); + // one straight span each way (its joint knots removed): the plane's own 2 x 2 form + assert_eq!(grown.gw_nurb().unwrap().cpts().len(), 2); + assert_eq!(grown.gw_nurb().unwrap().u_knots().len(), 4); + for (u, v) in [(0.0, 0.0), (10.0, 0.0), (3.5, 7.25), (-4.0, 14.0)] { + let uv = Point2::new(u, v); + // the plane's own parameters: origin + u X + v Y at scale 1 + let at = Vec3::new(u, v, 0.0); + assert!( + (grown.evaluate_point(uv).unwrap() - at).length() <= 1e-13, + "{uv:?}" + ); + } + // the across length is the distance on every side, in the plane's own parameters + let d = grown.gw_nurb().unwrap(); + // (the extension parameter solved to 1e-12 of the distance) + let ends = |k: &[f64]| (k[0], *k.last().unwrap()); + for (lo, hi) in [ends(d.u_knots()), ends(d.v_knots())] { + assert!( + (lo + 4.33).abs() <= 1e-9 && (hi - 14.33).abs() <= 1e-9, + "{lo} {hi}" + ); + } + for a in e.sides { + let SideExtension::Extended(a) = a else { + unreachable!() + }; + assert!((a.enclosure.mid() - 4.33).abs() <= 1e-9 * 4.33, "{a:?}"); + } } #[test] -fn a_cone_extends_as_itself_along_up_to_its_apex_and_is_left_round_its_axis() { +fn a_cone_extends_its_own_net_along_up_to_its_apex_and_round_its_axis_a_quarter_of_its_open_turn() { + use crate::surface::implicit::{lies_on, AnalyticSupport}; let kk = k(); // radius 2 at the bottom, 1 at height 1: the apex at height 2, sqrt 2 slant past the top let pl = Axis2Placement::new(Vec3::ZERO, Vec3::Z, Vec3::X); @@ -859,31 +907,73 @@ fn a_cone_extends_as_itself_along_up_to_its_apex_and_is_left_round_its_axis() { c.adjust_step_uv_domain(&kk, &Extent2d::from_bounds(0.0, 0.0, 90.0, 1.0)) .unwrap(); let c = Surface::from(c); + let on = AnalyticSupport::check(&kk, &c).unwrap().implicit(); + let src = crate::blend::fillet_section::station::as_b_spline(&kk, &c) + .unwrap() + .unwrap() + .into_owned(); + let src_net = src.gw_nurb().unwrap().clone(); for how in [Continuity::Tangent, Continuity::Natural] { let e = c.extended_sides(&kk, &BoxSide::ALL, how, 5.0).unwrap(); - assert!(matches!(e.surface.kind(), SurfaceKind::Cone(_)), "{how:?}"); - let around = e - .sides - .iter() - .filter(|x| matches!(x, SideExtension::AroundAxis)) - .count(); - assert_eq!(around, 2, "{:?}", e.sides); - // one side along runs the whole 5, the other stops short of the apex - let mut lengths: Vec<(bool, f64)> = e - .sides - .iter() - .filter_map(|x| match x { - SideExtension::Extended(a) => Some((false, a.enclosure.mid())), - SideExtension::UpToLimit(a) => Some((true, a.enclosure.mid())), - _ => None, - }) - .collect(); + let SurfaceKind::BSpline(grown) = e.surface.kind() else { + panic!("{how:?}: not a B-spline"); + }; + // the source's circles' points and knots bit for bit (along the rulings the span is + // widened, one straight span: its end circles to rounding), the result on the cone + let out = grown.gw_nurb().unwrap(); + let (src_u, src_v) = (src_net.u_knots(), src_net.v_knots()); + assert!(src_u.iter().all(|x| out.u_knots().contains(x)), "{how:?}"); + assert_eq!(out.v_knots().len(), 4, "{how:?}: one straight span along"); + assert_eq!(out.cpts().len(), src_net.cpts().len() + 4, "{how:?}"); + for i in 0..=4 { + for j in 0..=4 { + let uv = Point2::new( + src_u[0] + f64::from(i) / 4.0 * (src_u.last().unwrap() - src_u[0]), + src_v[0] + f64::from(j) / 4.0 * (src_v.last().unwrap() - src_v[0]), + ); + let d = grown.evaluate_point(uv).unwrap() - src.evaluate_point(uv).unwrap(); + assert!(d.length() <= 1e-14, "{how:?} {uv:?}: {d:?}"); + } + } + assert!(lies_on(grown, &on).unwrap()); + // round the axis (u): each side a quarter of the 270 degrees left open, said so: the + // radius at the side's middle (1.5) times 3 pi / 8 + for i in [0, 1] { + let SideExtension::UpToLimit(a) = e.sides[i] else { + panic!("{:?}", e.sides); + }; + let want = 1.5 * 0.375 * PI; + assert!((a.enclosure.mid() - want).abs() <= 1e-9 * want, "{a:?}"); + } + // round the axis each end turned 67.5 degrees from the source's + let (d, sd) = (grown.gw_nurb().unwrap(), src_net.clone()); + let ends = + |n: &crate::nurbs::srf_basic::Surface| (n.u_knots()[0], *n.u_knots().last().unwrap()); + let vm = 0.5 * (sd.v_knots()[0] + *sd.v_knots().last().unwrap()); + let ((u0, u1), (s0, s1)) = (ends(d), ends(&sd)); + let crate::surface::implicit::Implicit::Cone { origin, axis, .. } = on else { + panic!("not a cone"); + }; + let flat = |s: &BSplineSurface, u: f64| { + let p = s.evaluate_point(Point2::new(u, vm)).unwrap() - origin; + p - axis * axis.dot(p) + }; + let deg = |a: Vec3, b: Vec3| a.cross(b).length().atan2(a.dot(b)).to_degrees(); + assert!((deg(flat(grown, u0), flat(&src, s0)) - 67.5).abs() <= 1e-6); + assert!((deg(flat(grown, u1), flat(&src, s1)) - 67.5).abs() <= 1e-6); + // along: one side runs the whole 5, the other stops at the apex + let along = |x: &SideExtension| match x { + SideExtension::Extended(a) => Some((false, a.enclosure.mid())), + SideExtension::UpToLimit(a) => Some((true, a.enclosure.mid())), + _ => None, + }; + let mut lengths: Vec<(bool, f64)> = e.sides[2..].iter().filter_map(along).collect(); lengths.sort_by(|a, b| a.1.total_cmp(&b.1)); let [(true, short), (false, long)] = lengths[..] else { panic!("{:?}", e.sides); }; - assert!((long - 5.0).abs() <= 1e-8, "{long}"); - assert!(short < 2f64.sqrt() && short > 2f64.sqrt() - 1e-6, "{short}"); + assert!((long - 5.0).abs() <= 1e-9 * 5.0, "{long}"); + assert!((short - 2f64.sqrt()).abs() <= 1e-9, "{short}"); } } @@ -1012,6 +1102,11 @@ fn ruled_with_end_weight_ratios_that_differ_is_c1_and_straight_up_to_its_reach() /// 1e-3); with `certified`, on its own net and bound ([`net_of_definition`]). The planes /// run to x = 20, so the feet have room along. fn round_blend(order: u32, certified: bool) -> (Kernel, ProceduralBlend) { + round_blend_with(order, certified, 0.5) +} + +/// [`round_blend`] with the G^n section's thumbweight rule `thumbweight`. +fn round_blend_with(order: u32, certified: bool, thumbweight: f64) -> (Kernel, ProceduralBlend) { use crate::blend::fillet_section::recertify::net_of_definition; use crate::blend::fillet_section::section::GnSectionFsg; use crate::geom::vector3d::Point3; @@ -1051,7 +1146,7 @@ fn round_blend(order: u32, certified: bool) -> (Kernel, ProceduralBlend) { offsets: [r, r], law: None, order, - thumbweight: if order >= 2 { 0.5 } else { 0.0 }, + thumbweight: if order >= 2 { thumbweight } else { 0.0 }, leader: &leader, seeds: [&seed, &seed], bound: 1e-3, @@ -1168,16 +1263,148 @@ fn a_round_continues_along_its_spine_by_a_span_of_its_definition() { } #[test] -fn a_gn_section_across_is_refused() { - let (kk, pb) = round_blend(2, false); - let s = Surface::from(pb); +fn a_side_inside_the_spine_runs_the_definitions_spans_then_past_its_end() { + use crate::nurbs::b_spline_surface::SurfParamType; + let (kk, pb) = round_blend(1, true); + // the half of the round up to x = 5 (u = 1/2, inside the spine's one span) + let (left, _) = pb.split_at(&kk, 0.5, SurfParamType::U).unwrap(); + let s = Surface::from(left.clone()); let d0 = s.natural_uv_domain_with_nurb(&kk).unwrap(); - match s.extended_by_distance(&kk, &d0, BoxSide::VMax, Continuity::Tangent, 0.1) { - Err(SolidError::Unsupported { .. }) => {} - other => panic!("{other:?}"), + assert!((d0.max().x - 0.5).abs() < 1e-15, "{d0:?}"); + for (d, ue, stations) in [(3.0, 0.8, 2), (8.0, 1.3, 3)] { + let e = s + .extended_by_distance(&kk, &d0, BoxSide::UMax, Continuity::Tangent, d) + .unwrap(); + meets(&e, d); + let SurfaceKind::ProceduralBlend(out) = e.surface.kind() else { + panic!("not a blend: {:?}", e.surface.kind()); + }; + let dd = e.surface.natural_uv_domain_with_nurb(&kk).unwrap(); + assert!((dd.max().x - ue).abs() < 1e-6, "{d}: {dd:?}"); + assert!( + out.bound() >= left.bound() && out.bound() < 1e-6, + "{}", + out.bound() + ); + for i in 0..=8 { + for j in 0..=8 { + let u = 0.5 + (dd.max().x - 0.5) * f64::from(i) / 8.0; + let v = f64::from(j) / 8.0; + let p = e.surface.evaluate_point(&kk, Vec2::new(u, v)).unwrap(); + let st = out.definition_parameters(Point2::new(u, v)); + let def = out.evaluate_definition(&kk, st.x, st.y, 0, 0).unwrap()[0][0]; + assert!(p.distance(def) <= out.bound(), "{u} {v}: {p:?} {def:?}"); + let rad = (p.y - 2.0).hypot(p.z - 2.0); + assert!((rad - 2.0).abs() <= out.bound() + 1e-12, "{p:?}"); + } + } + // inside the spine the definition is its own; past its end it gains a station + let sigma = out.definition().seeds[0].knots(None).unwrap().0; + assert_eq!(sigma.len(), stations, "{d}: {sigma:?}"); + let src = left.net.gw_nurb().unwrap(); + assert!(keeps(src, out.net.gw_nurb().unwrap(), BoxSide::UMax).unwrap()); } } +#[test] +fn a_gn_section_continues_past_a_foot_with_its_bound_grown_over_its_weight_in_closed_form() { + // the G^2 section at the constant thumbweight 1/2 + let (kk, pb) = round_blend(2, true); + let s = Surface::from(pb.clone()); + let d0 = s.natural_uv_domain_with_nurb(&kk).unwrap(); + let e = s + .extended_by_distance(&kk, &d0, BoxSide::VMax, Continuity::Tangent, 0.3) + .unwrap(); + meets(&e, 0.3); + let SurfaceKind::ProceduralBlend(out) = e.surface.kind() else { + panic!("not a blend: {:?}", e.surface.kind()); + }; + let d = e.surface.natural_uv_domain_with_nurb(&kk).unwrap(); + let te = d.max().y; + assert!(te > 1.0 && te < 2.0, "{d:?}"); + // the growth in closed form: T_m(2 te - 1), m = 3 n + 1 + the net's degree across, + // times the net's weight ratio, over W's least past the foot: A >= 1, B >= 1 and + // C = 1 - te / 2 at the end (W at most 1 on the source) + let src = pb.net.gw_nurb().unwrap(); + let m = 7 + src.q(); + let x = 2.0 * te - 1.0; + let (mut prev, mut cur) = (1.0_f64, x); + for _ in 1..m { + let next = 2.0 * x * cur - prev; + prev = cur; + cur = next; + } + let ws = |n: &crate::nurbs::srf_basic::Surface| { + let w: Vec = n + .cpts() + .iter() + .flatten() + .map(|c| if c.w == NOW { 1.0 } else { c.w }) + .collect(); + ( + w.iter().copied().fold(0.0_f64, f64::max), + w.iter().copied().fold(f64::INFINITY, f64::min), + ) + }; + let (vhi, _) = ws(src); + let (_, vlo) = ws(out.net.gw_nurb().unwrap()); + let c = 1.0 - 0.5 * te; + let want = pb.bound() * cur * vhi / vlo / (c * c); + assert!( + out.bound() >= want && out.bound() <= want * (1.0 + 1e-9) + 1e-15 + e.extension, + "{} vs {want}", + out.bound() + ); + // the continued part within the bound of the definition continued + for i in 0..=8 { + for j in 0..=8 { + let u = f64::from(i) / 8.0; + let v = 1.0 + (te - 1.0) * f64::from(j) / 8.0; + let p = e.surface.evaluate_point(&kk, Vec2::new(u, v)).unwrap(); + let st = out.definition_parameters(Point2::new(u, v)); + let def = out.evaluate_definition(&kk, st.x, st.y, 0, 0).unwrap()[0][0]; + assert!(p.distance(def) <= out.bound(), "{u} {v}: {p:?} {def:?}"); + } + } + // the source part is the source's net, bit for bit + assert!(keeps(src, out.net.gw_nurb().unwrap(), BoxSide::VMax).unwrap()); + // a long ask stays short of the section's weight's zero (C = 0 at t = 2): met there + // with the bound grown by its closed form, or refused + match s.extended_by_distance(&kk, &d0, BoxSide::VMax, Continuity::Tangent, 50.0) { + Ok(long) => { + let d = long.surface.natural_uv_domain_with_nurb(&kk).unwrap(); + assert!(d.max().y < 2.0, "{d:?}"); + } + Err(SolidError::ExtendRefused { .. }) => {} + Err(other) => panic!("{other:?}"), + } +} + +#[test] +fn an_equal_area_gn_section_continues_over_its_thumbweights_proved_enclosure() { + use crate::blend::fillet_section::recertify::recertify_with_thumbweight; + use crate::blend::fillet_section::section::GN_THUMBWEIGHT_EQUAL_AREA; + let (kk, pb) = round_blend_with(2, true, GN_THUMBWEIGHT_EQUAL_AREA); + // the certificate's interval Newton roots enclose the thumbweight: one w along this + // straight round, below 1 + let (_, w) = recertify_with_thumbweight(&kk, &pb).unwrap(); + let [lo, hi] = w.unwrap(); + assert!(lo <= hi && hi < 1.0 && hi - lo < 1e-6, "{lo} {hi}"); + let s = Surface::from(pb.clone()); + let d0 = s.natural_uv_domain_with_nurb(&kk).unwrap(); + let e = s + .extended_by_distance(&kk, &d0, BoxSide::VMax, Continuity::Tangent, 0.3) + .unwrap(); + meets(&e, 0.3); + let SurfaceKind::ProceduralBlend(out) = e.surface.kind() else { + panic!("not a blend"); + }; + assert!(out.bound().is_finite() && out.bound() >= pb.bound()); + // a constant thumbweight is its own + let (kk, pb) = round_blend(2, true); + assert!(recertify_with_thumbweight(&kk, &pb).unwrap().1.is_none()); +} + #[test] fn a_round_along_past_its_supports_is_refused_naming_a_distance_that_works() { // the planes run to x = 20 and the spine ends at x = 10: past x = 20 a foot leaves its @@ -1197,3 +1424,66 @@ fn a_round_along_past_its_supports_is_refused_naming_a_distance_that_works() { .unwrap(); meets(&e, named); } + +/// The octant ball patch of radius 2 about `(1, -2, 0.5)`, its pole on +x. +fn octant_ball() -> (Surface, Vec3, f64) { + let (c, r) = (Vec3::new(1.0, -2.0, 0.5), 2.0); + let dirs = [Vec3::X, Vec3::Y, Vec3::Z]; + let (net, _, _) = + crate::blend::blend_corner::ball::ball_triangle_patch(c, r, &dirs, 0, false).unwrap(); + let patch = crate::surface::ball_patch::BallPatch::create(&net, c, r).unwrap(); + (Surface::from(patch), c, r) +} + +#[test] +fn a_ball_patch_extends_its_own_net_on_its_sphere_up_to_its_caps() { + use crate::surface::implicit::{lies_on, Implicit}; + let kk = k(); + let (s, c, r) = octant_ball(); + let SurfaceKind::BallPatch(src) = s.kind() else { + panic!("not a ball patch"); + }; + let src_net = src.net().gw_nurb().unwrap().clone(); + // a short ask is met on every side but the pole's; a long one stops at the caps: round + // the pole r times a quarter of the 270 degrees left open, away half the 90 degrees + // left from the third side to the antipode + for (d, round, away) in [(0.1, 0.1, 0.1), (100.0, r * 0.375 * PI, r * 0.25 * PI)] { + let e = s + .extended_sides(&kk, &BoxSide::ALL, Continuity::Ruled, d) + .unwrap(); + let SurfaceKind::BallPatch(out) = e.surface.kind() else { + panic!("{d}: not a ball patch"); + }; + assert_eq!(e.sides[2], SideExtension::Collapsed, "{:?}", e.sides); + for (i, want) in [(0, round), (1, round), (3, away)] { + let a = match e.sides[i] { + SideExtension::Extended(a) if d < 1.0 => a, + SideExtension::UpToLimit(a) if d > 1.0 => a, + other => panic!("{d}: side {i}: {other:?}"), + }; + assert!( + (a.enclosure.mid() - want).abs() <= 1e-9 * want, + "{d}: {i}: {a:?}" + ); + } + // the source's net bit for bit, the whole result on the sphere + let out_net = out.net().gw_nurb().unwrap(); + let kept = src_net.cpts().iter().flatten().all(|p| { + out_net + .cpts() + .iter() + .flatten() + .any(|q| p.x == q.x && p.y == q.y && p.z == q.z && p.w == q.w) + }); + assert!(kept, "{d}: the source's net moved"); + assert!(lies_on( + out.net(), + &Implicit::Sphere { + centre: c, + radius: r + } + ) + .unwrap()); + assert!((out.radius() - r).abs() == 0.0 && out.centre() == c); + } +} diff --git a/crates/solidcore/src/surface/side.rs b/crates/solidcore/src/surface/side.rs index 9f49a7a404..eca0a182d9 100644 --- a/crates/solidcore/src/surface/side.rs +++ b/crates/solidcore/src/surface/side.rs @@ -18,7 +18,8 @@ use crate::geom::vector2d::Point2; use crate::geom::vector3d::Vec3; use crate::kernel::Kernel; use crate::nurbs::b_spline_surface::{BSplineSurface, IsoCurveOptions, SurfParamType}; -use crate::nurbs::math::{gamma, NOW}; +use crate::nurbs::math::{gamma, CPoint, NOW}; +use crate::nurbs::srf_basic::Surface as Net; use crate::surface::implicit::AnalyticSupport; use crate::surface::surface::{Surface, SurfaceKind}; @@ -137,6 +138,46 @@ impl Bracket { } } +/// `form`'s net translated by `P0` (the iso's first point, as computed), and the +/// translation's rounding: each homogeneous point `(X w - P0 w, w)`, its coordinates within +/// `gamma_2 (|x| + |P0| |w|)` of the exact translation's (a product and a difference; the +/// weights exact). A length is the same on the translated surface (`S - P0`, `P0` a fixed +/// vector); measured there, the points' and their rounding's magnitudes are the net's size +/// about `P0`, not its distance from the origin. +/// +/// # Errors +/// +/// `SM_ERR` for a form without a NURBS net; those of the net. +fn translated(form: &BSplineSurface, p0: Vec3, op: &'static str) -> Result<(BSplineSurface, f64)> { + let nurb = form.gw_nurb().ok_or_else(|| SolidError::sm_err(op))?; + let cpts = nurb.cpts(); + let size = p0.x.abs().max(p0.y.abs()).max(p0.z.abs()); + let g = gamma(2); + let mut shift = 0.0_f64; + let net: Vec> = cpts + .iter() + .map(|row| { + row.iter() + .map(|c| { + let w = if c.w == NOW { 1.0 } else { c.w }; + let mag = c.x.abs().max(c.y.abs()).max(c.z.abs()); + shift = shift.max(g * (mag + size * w.abs())); + CPoint::new(c.x - p0.x * w, c.y - p0.y * w, c.z - p0.z * w, c.w) + }) + .collect() + }) + .collect(); + let (p, q) = nurb.degrees(); + let out = Net::new(net, p, q, nurb.u_knots().to_vec(), nurb.v_knots().to_vec())?; + if !shift.is_finite() { + return Err(SolidError::invalid_input( + op, + "the net's translation is not bounded", + )); + } + Ok((BSplineSurface::from_gw_nurb(&out), shift)) +} + /// A rational Bezier piece in homogeneous coordinates `(w x, w y, w z, w)`. type Piece = Vec<[f64; 4]>; @@ -180,21 +221,156 @@ struct Rounding { homogeneous: f64, } +/// The largest knot magnitude and the least nonzero span of the knots of `kv` (the full +/// vector) from `margin` indices under the last knot below `a` to `margin` over the first +/// above `b`: what the knot insertions and the Bezier decomposition that make the pieces +/// over `[a, b]` read of `kv` (a segment's coefficients `(t - u_i) / (u_(i+p) - u_i)` take +/// `i` from its start group's last index, the one before it for its first points, to `p` +/// past it; the iso's extraction the knots within `q` of its span). +fn knot_window(kv: &[f64], a: f64, b: f64, margin: usize) -> (f64, f64) { + let lo = kv + .iter() + .position(|&x| x > a) + .unwrap_or(kv.len()) + .saturating_sub(margin); + let hi = kv + .iter() + .rposition(|&x| x < b) + .unwrap_or(0) + .saturating_add(margin); + let (mut kmax, mut smin) = (0.0_f64, f64::INFINITY); + let window: Vec = kv + .iter() + .enumerate() + .filter(|&(i, _)| i >= lo && i <= hi) + .map(|(_, &x)| x) + .collect(); + for (&x, &y) in window.iter().zip(window.iter().skip(1)) { + kmax = kmax.max(x.abs()).max(y.abs()); + if y > x { + smin = smin.min(y - x); + } + } + (kmax, smin) +} + +/// [`knot_window`] for the iso's extraction at `t` of a direction of degree `q` (full knot +/// vector `kv`): knot insertion at `t` in its span `[u_s, u_(s+1))` forms its coefficients +/// `(t - u_i) / (u_(i+q-j+1) - u_i)` from the knots `u_(s-q+1) .. u_(s+q)` only. +fn extraction_window(kv: &[f64], t: f64, q: usize) -> (f64, f64) { + let s = kv.iter().rposition(|&x| x <= t).unwrap_or(0); + let (lo, hi) = (s.saturating_sub(q), s.saturating_add(q)); + let window: Vec = kv + .iter() + .enumerate() + .filter(|&(i, _)| i >= lo && i <= hi) + .map(|(_, &x)| x) + .collect(); + let (mut kmax, mut smin) = (t.abs(), f64::INFINITY); + for (&x, &y) in window.iter().zip(window.iter().skip(1)) { + kmax = kmax.max(x.abs()).max(y.abs()); + if y > x { + smin = smin.min(y - x); + } + } + (kmax, smin) +} + +/// [`knot_window`] for the Bezier decomposition along the pieces over `[a, b]` of a curve of +/// degree `p` (full knot vector `kv`): the decomposition forms a coefficient only where it +/// inserts, at an interior knot of multiplicity below `p`, and each of its denominators +/// `u_(i+j) - u_i` there is at least the length of the segment that knot ends. So the least +/// span read is the least such segment's in the window; a window whose interior knots are all +/// of multiplicity `p` or more (Bezier joints already) forms none (`INFINITY`). +fn decomposed_window(kv: &[f64], a: f64, b: f64, p: usize) -> (f64, f64) { + let margin = p.saturating_add(2); + let (kmax, _) = knot_window(kv, a, b, margin); + let lo = kv + .iter() + .position(|&x| x > a) + .unwrap_or(kv.len()) + .saturating_sub(margin); + let hi = kv + .iter() + .rposition(|&x| x < b) + .unwrap_or(0) + .saturating_add(margin); + // the distinct knots in the window with their multiplicities in the whole vector + let mut groups: Vec<(f64, usize)> = Vec::new(); + for (i, &x) in kv.iter().enumerate() { + match groups.last_mut() { + Some((y, m)) if *y == x => *m = m.saturating_add(1), + _ if i >= lo && i <= hi => groups.push((x, 1)), + _ => {} + } + } + let (first, last) = (kv.first().copied(), kv.last().copied()); + let mut smin = f64::INFINITY; + for (&(x, _), &(y, my)) in groups.iter().zip(groups.iter().skip(1)) { + let interior = Some(y) != first && Some(y) != last; + if interior && my < p && y > x { + smin = smin.min(y - x); + } + } + (kmax, smin) +} + +/// The control indices of degree `deg` on the knots `k` whose basis function's support +/// `[k_i, k_(i + deg + 1)]` meets `[lo, hi]` (closed: a support ending at `lo` is kept). +fn support_meets(k: &[f64], deg: usize, lo: f64, hi: f64) -> std::ops::RangeInclusive { + let count = k.len().saturating_sub(deg + 1); + let meets = |i: usize| matches!((k.get(i), k.get(i + deg + 1)), (Some(&a), Some(&b)) if a <= hi && b >= lo); + let first = (0..count).find(|&i| meets(i)); + let last = (0..count).rev().find(|&i| meets(i)); + match (first, last) { + (Some(f), Some(l)) => f..=l, + // none: every point, as the whole net's bound + _ => 0..=count.saturating_sub(1), + } +} + impl Rounding { + /// [`Rounding::of`] for the pieces of the iso at `fixed` (constant v for `along_u`) + /// over `[from, to]` alone: the knots read are only those near them. Knot insertion and + /// the Bezier decomposition are local: a segment's points are combinations whose + /// coefficients take the knots from its start's segment's start to `p` past its end + /// (`(t - u_i) / (u_(i+p) - u_i)`, i from the segment's start index), and the iso's + /// extraction at `fixed` the knots within `q` of it. A span far from the pieces (a sliver + /// at the net's other end) does not touch them. + fn local( + form: &BSplineSurface, + along_u: bool, + fixed: f64, + from: f64, + to: f64, + op: &'static str, + ) -> Result { + let nurb = form.gw_nurb().ok_or_else(|| SolidError::sm_err(op))?; + let (p, q) = nurb.degrees(); + let (held, along, dp, dq) = if along_u { + (nurb.v_knots(), nurb.u_knots(), p, q) + } else { + (nurb.u_knots(), nurb.v_knots(), q, p) + }; + let (a, b) = (from.min(to), from.max(to)); + // each direction's coefficients read its own knots and parameters + let (k1, s1) = decomposed_window(along, a, b, dp); + let (k2, s2) = extraction_window(held, fixed, dq); + let ratio = (k1.max(from.abs()).max(to.abs()) / s1).max(k2 / s2); + // and only the control points the pieces are made of: those whose basis function's + // support meets `[from, to]` along and holds `fixed` across (a far row's size, the + // other end of a long spine, enters no piece) + let along_ix = support_meets(along, dp, a, b); + let held_ix = support_meets(held, dq, fixed, fixed); + let keep = |i: usize, j: usize| { + let (ia, ih) = if along_u { (i, j) } else { (j, i) }; + along_ix.contains(&ia) && held_ix.contains(&ih) + }; + Self::with_ratio_in(nurb, ratio, op, &keep) + } + fn of(form: &BSplineSurface, extra: &[f64], op: &'static str) -> Result { let nurb = form.gw_nurb().ok_or_else(|| SolidError::sm_err(op))?; - let (mut h, mut wmin, mut wmax) = (0.0_f64, f64::INFINITY, 0.0_f64); - for row in nurb.rows() { - for c in row { - let w = if c.w == NOW { 1.0 } else { c.w }; - if w.is_nan() || w <= 0.0 { - return Err(SolidError::invalid_input(op, "a weight is not positive")); - } - h = h.max(c.x.abs()).max(c.y.abs()).max(c.z.abs()); - wmin = wmin.min(w); - wmax = wmax.max(w); - } - } let (mut kmax, mut smin) = (0.0_f64, f64::INFINITY); for kv in [nurb.u_knots(), nurb.v_knots()] { for (&a, &b) in kv.iter().zip(kv.iter().skip(1)) { @@ -207,7 +383,45 @@ impl Rounding { for &x in extra { kmax = kmax.max(x.abs()); } - let level = gamma(3) * (1.0 + 6.0 * kmax / smin) * 2.0; + Self::with_ratio(nurb, kmax / smin, op) + } + + /// The rounding of a net whose coefficients' knots and parameters are at most `ratio` + /// times the least span they divide by (`K / S`). + fn with_ratio( + nurb: &crate::nurbs::srf_basic::Surface, + ratio: f64, + op: &'static str, + ) -> Result { + Self::with_ratio_in(nurb, ratio, op, &|_, _| true) + } + + /// [`Rounding::with_ratio`] over the control points `keep` names (`[i along u][j along + /// v]`): the pieces' points are combinations of those alone. + fn with_ratio_in( + nurb: &crate::nurbs::srf_basic::Surface, + ratio: f64, + op: &'static str, + keep: &dyn Fn(usize, usize) -> bool, + ) -> Result { + let (mut h, mut wmin, mut wmax) = (0.0_f64, f64::INFINITY, 0.0_f64); + for (i, row) in nurb.rows().enumerate() { + for c in row + .iter() + .enumerate() + .filter(|&(j, _)| keep(i, j)) + .map(|(_, c)| c) + { + let w = if c.w == NOW { 1.0 } else { c.w }; + if w.is_nan() || w <= 0.0 { + return Err(SolidError::invalid_input(op, "a weight is not positive")); + } + h = h.max(c.x.abs()).max(c.y.abs()).max(c.z.abs()); + wmin = wmin.min(w); + wmax = wmax.max(w); + } + } + let level = gamma(3) * (1.0 + 6.0 * ratio) * 2.0; let scale = h * wmax / (wmin * wmin); if level.is_finite() && scale.is_finite() { Ok(Self { @@ -565,10 +779,18 @@ impl Surface { converged: true, }); } - let rounding = Rounding::of(&form, &[fixed, from, to], OP_LENGTH)?; + let start = if side.is_u() { + Point2::new(from, fixed) + } else { + Point2::new(fixed, from) + }; + let (form, shift) = translated(&form, form.evaluate_point(start)?, OP_LENGTH)?; + let rounding = Rounding::local(&form, side.is_u(), fixed, from, to, OP_LENGTH)?; let (pieces, levels) = iso_pieces(k, &form, side.is_u(), fixed, from, to)?; - // the pieces' homogeneous points as computed, within this of the exact ones - let made = rounding.homogeneous(levels); + // the pieces' homogeneous points as computed, within this of the exact ones (of the + // net translated exactly), the translation's rounding carried through the levels' + // convex combinations unchanged + let made = rounding.homogeneous(levels) + shift; let (mut chords, mut polygons, mut leaves) = (0.0_f64, 0.0_f64, 0_usize); let (mut chord_cuts, mut polygon_cuts, mut moved) = (0.0_f64, 0.0_f64, 0.0_f64); let mut converged = true; diff --git a/crates/solidcore/src/surface/surface/p2.rs b/crates/solidcore/src/surface/surface/p2.rs index 22856bfe72..199e54a810 100644 --- a/crates/solidcore/src/surface/surface/p2.rs +++ b/crates/solidcore/src/surface/surface/p2.rs @@ -3826,6 +3826,13 @@ impl Surface { } } + /// An exact round's definition carried onto `continuation`, this surface's own NURBS + /// net continued Natural (`surface::extend`'s analytic kinds): still on S, as + /// [`Surface::create_extended_surface`]'s `CInfinity` carries it (#2487). + pub(crate) fn carry_round_on_continuation(&self, continuation: Self) -> Self { + self.carry_exact_round_on_extension(continuation, ContinuityType::CInfinity) + } + /// [`Surface::create_extended_surface`] of the kind alone. fn create_extended_surface_of_kind( &self, diff --git a/crates/solidcore/src/topology/brep/p2.rs b/crates/solidcore/src/topology/brep/p2.rs index f1c42597d7..98b33d34a7 100644 --- a/crates/solidcore/src/topology/brep/p2.rs +++ b/crates/solidcore/src/topology/brep/p2.rs @@ -2912,7 +2912,13 @@ pub fn mate_lamina_edge_to_surface( "negative surface extension distance", )); } - surf.create_extended_surface(k, ext_dist, ContinuityType::G1R)? + surf.extended_sides( + k, + &crate::surface::side::BoxSide::ALL, + crate::surface::extend::Continuity::Ruled, + ext_dist, + )? + .surface }; let ext_domain = ext_surface.natural_uv_domain_with_nurb(k)?; diff --git a/crates/solidcore/src/topology/brep/p3.rs b/crates/solidcore/src/topology/brep/p3.rs index 3e5e9cf507..faa80aa282 100644 --- a/crates/solidcore/src/topology/brep/p3.rs +++ b/crates/solidcore/src/topology/brep/p3.rs @@ -124,7 +124,6 @@ use crate::geom::context::{ }; use crate::geom::extent1d::Extent1d; use crate::geom::extent3d::Extent3d; -use crate::geom::geom_utility::p2::ContinuityType; use crate::geom::object::{self, NotifyOperation, ObjectRef}; use crate::geom::tol::{self, SmallTopoType}; use crate::kernel::Kernel; @@ -4264,21 +4263,16 @@ fn turn_to_nurbs_in( size = 1.0; } let extend_distance = 0.02 * size; - let SurfaceKind::BSpline(new_bss) = new_surface.kind() else { - return Err(ner()); - }; - let ext = new_bss.create_extended_surface(k, extend_distance, ContinuityType::G1)?; + let ext = new_surface + .extended_sides( + k, + &crate::surface::side::BoxSide::ALL, + crate::surface::extend::Continuity::Tangent, + extend_distance, + )? + .surface; let face_ref = (env.key_of)(EntityRef::Face(f)); - p2::replace_surface( - k, - b, - env, - f, - Surface::from(ext), - face_ref, - OldSurface::Delete, - true, - )?; + p2::replace_surface(k, b, env, f, ext, face_ref, OldSurface::Delete, true)?; continue; }; if surface.type_code() == BSplineSurface::TYPE { diff --git a/crates/solidcore/tests/issue_2484_shell_keeps_ball_patches.rs b/crates/solidcore/tests/issue_2484_shell_keeps_ball_patches.rs index a03417e29b..0ab0b2c0e5 100644 --- a/crates/solidcore/tests/issue_2484_shell_keeps_ball_patches.rs +++ b/crates/solidcore/tests/issue_2484_shell_keeps_ball_patches.rs @@ -4,8 +4,8 @@ //! eight corner faces is on the sphere of radius 1 concentric with its outer corner. //! The shell's extended offset (`SmOffsetExecutive::GetCreateExtendedSurface`) used to //! extend those patches tangentially (`SM_CT_G1R`) into plain B-splines, and 7 of the 8 -//! inner corners lost their sphere; with every fix on, each is a ball patch on its -//! sphere, the volume is the closed form's and the body checks clean. +//! inner corners lost their sphere; now each is a ball patch on its sphere (its own net +//! continued, flag or not), the volume is the closed form's and the body checks clean. #![allow( clippy::unwrap_used, clippy::expect_used, @@ -209,12 +209,13 @@ fn the_shell_has_the_closed_form_volume_and_checks_clean() { assert!(lines.is_empty(), "{lines:?}"); } -/// Without `offset_extends_ball_patches_as_themselves` the extended offset extends a ball -/// patch tangentially, a plain B-spline: only the inner corner whose trims stay inside -/// its patch's domain to `EFF_ZERO` (corner 1, at (27, 3, 3)) keeps its patch; the other -/// seven reach 1.6e-12 to 2.4e-11 past it in their trims' rounding and keep the B-spline. +/// `offset_extends_ball_patches_as_themselves` no longer changes the extended offset: it +/// always extends a ball patch as itself, its own net continued on its sphere +/// (`Surface::extended_sides`), so with the flag off every inner corner still keeps its +/// patch, on its sphere (before, a tangent extension into a plain B-spline kept only +/// corner 1's). #[test] -fn without_the_fix_seven_inner_corners_lose_their_sphere() { +fn without_the_fix_every_inner_corner_still_keeps_its_sphere() { let mut k = all_fixes_but("offset_extends_ball_patches_as_themselves"); let shell = shelled_box(&mut k); let mut kept = Vec::new(); @@ -222,9 +223,10 @@ fn without_the_fix_seven_inner_corners_lose_their_sphere() { let faces = corner_faces(&k, &shell, corner, R - WALL); assert_eq!(faces.len(), 1, "corner {corner}"); let s = face2::surface(&shell, faces[0]).unwrap().unwrap(); - if matches!(s.kind(), SurfaceKind::BallPatch(_)) { + if let SurfaceKind::BallPatch(p) = s.kind() { + assert!(off_sphere(p) <= 1e-12, "corner {corner}"); kept.push(corner); } } - assert_eq!(kept, [1]); + assert_eq!(kept, (0..8).collect::>()); } diff --git a/crates/solidcore/tests/issue_398_shell_box_rounds_chamfer_fillet.rs b/crates/solidcore/tests/issue_398_shell_box_rounds_chamfer_fillet.rs index a7ac3f2091..d51284ca95 100644 --- a/crates/solidcore/tests/issue_398_shell_box_rounds_chamfer_fillet.rs +++ b/crates/solidcore/tests/issue_398_shell_box_rounds_chamfer_fillet.rs @@ -256,8 +256,11 @@ fn issue_398_open_top_shell_with_edges_on_analytics_exact() { assert_eq!(got.face_list().len(), 51); let want = closed_form(Some(30.0)); assert!((v - want).abs() <= 1e-6 * want, "{v} vs {want}"); - // 2.3e-8 above the closed form; 13701.735330849975 (8.4e-8 above) with the flag off - let pinned = 13701.73449509504; + // 2.8e-9 above the closed form: the offset's face extensions are exact (the source bit + // for bit, no least squares fit or knot removal moving a face's own part); it was + // 13701.73449509504 (2.3e-8 above) on the fitted extensions, 13701.735330849975 (8.4e-8 + // above) with the flag off + let pinned = 13701.73421803366; assert!((v - pinned).abs() <= 1e-9 * pinned, "{v} vs {pinned}"); } diff --git a/crates/solidcore/tests/wr698_procedural_blend_definition.rs b/crates/solidcore/tests/wr698_procedural_blend_definition.rs index 23477e4c45..ad6456668f 100644 --- a/crates/solidcore/tests/wr698_procedural_blend_definition.rs +++ b/crates/solidcore/tests/wr698_procedural_blend_definition.rs @@ -422,12 +422,14 @@ fn partials_past_the_series_are_refused() { /// evaluators and box (fork 0136260's values on the same definitions, whose net is the /// unit square z = 0): the box is the net's grown by the bound, so the blend is not /// planar within 1e-6 although its net is; the degeneracy tests, the largest analytic -/// domain and the extension follow the fork's. +/// domain follow the fork's; the extension is the exact one (`extended_sides`). #[test] fn the_base_class_paths_take_the_blend_s_box_and_points() { + use solidcore::error::SolidError; use solidcore::geom::extent2d::Extent2d; - use solidcore::geom::geom_utility::p2::ContinuityType; use solidcore::math::EFF_ZERO_DEG; + use solidcore::surface::extend::{Continuity, SideExtension}; + use solidcore::surface::side::BoxSide; let k = kernel(); for (bumpy, n, w, law) in [ (false, 1, 0.0, false), @@ -456,13 +458,47 @@ fn the_base_class_paths_take_the_blend_s_box_and_points() { (m.min().x, m.min().y, m.max().x, m.max().y), (-100.0, -100.0, 101.0, 101.0) ); - let ext = s - .create_extended_surface(&k, 0.1, ContinuityType::G1) - .unwrap(); - let ed = ext.natural_uv_domain_with_nurb(&k).unwrap(); - assert!(close(ed.min().x, -0.1) && close(ed.max().y, 1.1)); - let p0 = ext.evaluate_point(&k, ed.min()).unwrap(); - assert!((p0 - Point3::new(-0.1, -0.1, 0.0)).length() < 1e-15); + // the exact extension (`extended_sides`): across (v, the section's t) the toy net + // continued with its bound grown in closed form; along (u) the toy net is not the + // construction's cubic Hermite net, so no span of the definition joins it + // (refused); the equal area G^n section's thumbweight is proved only on a net of + // its definition (refused) + let across = s.extended_sides( + &k, + &[BoxSide::VMin, BoxSide::VMax], + Continuity::Tangent, + 0.1, + ); + if n == 2 { + assert!( + matches!(across, Err(SolidError::InvalidInput { .. })), + "{across:?}" + ); + } else { + let ext = across.unwrap(); + assert!(matches!(ext.sides[2], SideExtension::Extended(_))); + assert!(matches!(ext.sides[3], SideExtension::Extended(_))); + let ed = ext.surface.natural_uv_domain_with_nurb(&k).unwrap(); + assert!(close(ed.min().y, -0.1) && close(ed.max().y, 1.1), "{ed:?}"); + assert!(ed.min().x == 0.0 && ed.max().x == 1.0, "{ed:?}"); + // the net continued: the unit square's corner moved by 0.1 across + let solidcore::surface::surface::SurfaceKind::ProceduralBlend(pb) = ext.surface.kind() + else { + panic!("not a blend"); + }; + let p0 = pb.net().evaluate_point(ed.min()).unwrap(); + assert!( + (p0 - Point3::new(0.0, -0.1, 0.0)).length() < 1e-12, + "{p0:?}" + ); + } + for side in [BoxSide::UMin, BoxSide::UMax] { + let along = s.extended_sides(&k, &[side], Continuity::Tangent, 0.1); + assert!( + matches!(along, Err(SolidError::Unsupported { .. })), + "{along:?}" + ); + } } } diff --git a/crates/solidcore/tests/wr_344_half_frustum_shell_keeps_its_rim.rs b/crates/solidcore/tests/wr_344_half_frustum_shell_keeps_its_rim.rs index 4d4d2b8583..f63e07fcd4 100644 --- a/crates/solidcore/tests/wr_344_half_frustum_shell_keeps_its_rim.rs +++ b/crates/solidcore/tests/wr_344_half_frustum_shell_keeps_its_rim.rs @@ -5,7 +5,8 @@ //! built as two half faces (`create_partial_cone` 0..180 and 180..360, united) is //! shelled inward, open at its top. With `offset_extended_edges_mated` the mated pairs //! are related and the shell has its seven faces, twelve edges and eight vertices; by -//! default these frustums come out with extra edges and vertices on the rim. +//! default these frustums came out with extra edges and vertices on the rim; on the exact +//! surface extensions the default keeps the rim too. //! //! The shell's volume is held to 1e-6 of the closed form: the frustum less the cavity //! bounded by the inner cone (the wall moved `t` along its normal, its radius smaller by @@ -81,13 +82,21 @@ fn wr_344_mated_rims_are_not_split_with_the_fix() { } } -/// Upstream intersects the mated pairs again and splits the rim. +/// Upstream intersects the mated pairs again and splits the rim. On the exact extensions +/// (a cone continued as its own net, its rulings one straight span, the offsets' extended +/// bases meeting along the mated edges) the default path no longer splits it: the same seven +/// faces, twelve edges and eight vertices as the fix, a manifold solid of the closed form's +/// volume to 1e-6. #[test] -fn wr_344_default_splits_the_rim() { +fn wr_344_default_keeps_the_rim_on_exact_extensions() { for case in CASES { - let (_, b) = shelled(false, case); - let (f, e, v) = counts(&b); - assert_eq!(f, 7, "{case:?}"); - assert!(e > 12 && v > 8, "{case:?}: {e} edges, {v} vertices"); + let (k, mut b) = shelled(false, case); + assert!(brep_is_manifold_solid(&b).unwrap(), "{case:?}"); + assert_eq!(counts(&b), (7, 12, 8), "{case:?}"); + let (v, exact) = (volume(&k, &mut b).unwrap(), closed_form(case)); + assert!( + ((v - exact) / exact).abs() < 1e-6, + "{case:?}: {v} against {exact}" + ); } }