diff --git a/crates/difftest/src/units/surface_surface_point_solve_p1.rs b/crates/difftest/src/units/surface_surface_point_solve_p1.rs index 4fdd57f9f2..87aa466c0b 100644 --- a/crates/difftest/src/units/surface_surface_point_solve_p1.rs +++ b/crates/difftest/src/units/surface_surface_point_solve_p1.rs @@ -240,7 +240,7 @@ fn cpp(i: &Inputs) -> Result { // ---- Rust side ---------------------------------------------------------------------- -fn wrapped_of(wrap: bool, s: &SurfIn) -> Result { +pub(crate) fn wrapped_of(wrap: bool, s: &SurfIn) -> Result { let s = surface_in(s)?; Ok(if wrap { Surface::from(StepSurface::new(s, true)) @@ -249,15 +249,15 @@ fn wrapped_of(wrap: bool, s: &SurfIn) -> Result { }) } -fn p2(uv: [f64; 2]) -> Point2 { +pub(crate) fn p2(uv: [f64; 2]) -> Point2 { Point2::new(uv[0], uv[1]) } -fn push2(o: &mut Vec, v: Point2) { +pub(crate) fn push2(o: &mut Vec, v: Point2) { o.extend([v.x, v.y]); } -fn solver_op(code: u8) -> Result { +pub(crate) fn solver_op(code: u8) -> Result { SolverOperation::ALL .get(usize::from(code)) .copied() @@ -485,10 +485,10 @@ fn rust(i: &Inputs) -> Result { // ---- generated cases ---------------------------------------------------------------- /// The supported solver operations. -const SOLVERS: [u8; 5] = [0, 3, 4, 5, 10]; +pub(crate) const SOLVERS: [u8; 5] = [0, 3, 4, 5, 10]; /// A point of `d`: a corner, a point on a side, or inside. -fn gen_uv(rng: &mut Rng, d: [f64; 4]) -> [f64; 2] { +pub(crate) fn gen_uv(rng: &mut Rng, d: [f64; 4]) -> [f64; 2] { let pick = |rng: &mut Rng, a: f64, b: f64| -> f64 { match rng.below(6) { 0 => a, @@ -500,7 +500,7 @@ fn gen_uv(rng: &mut Rng, d: [f64; 4]) -> [f64; 2] { } /// A point in `d` or, half the time, somewhat outside it. -fn gen_uv_any(rng: &mut Rng, d: [f64; 4]) -> [f64; 2] { +pub(crate) fn gen_uv_any(rng: &mut Rng, d: [f64; 4]) -> [f64; 2] { let [u, v] = gen_uv(rng, d); if rng.chance(0.5) { return [u, v]; @@ -514,7 +514,7 @@ fn gen_uv_any(rng: &mut Rng, d: [f64; 4]) -> [f64; 2] { /// A test point for the surface: on it, slightly or well off it along its normal, or /// anywhere around it. -fn gen_test(rng: &mut Rng, s: &Surface, d: [f64; 4], size: f64) -> [f64; 3] { +pub(crate) fn gen_test(rng: &mut Rng, s: &Surface, d: [f64; 4], size: f64) -> [f64; 3] { let uv = gen_uv(rng, d); let Ok((p, su, sv)) = s.evaluate_1st_derivatives(p2(uv), false, false) else { return gvec(rng, size * 2.0); @@ -539,7 +539,7 @@ fn gen_test(rng: &mut Rng, s: &Surface, d: [f64; 4], size: f64) -> [f64; 3] { /// A test point on the surface normal at `uv`: on the surface, or off it on either /// side by a small to large distance. -fn test_on_normal(rng: &mut Rng, s: &Surface, uv: [f64; 2], size: f64) -> [f64; 3] { +pub(crate) fn test_on_normal(rng: &mut Rng, s: &Surface, uv: [f64; 2], size: f64) -> [f64; 3] { let Ok((p, su, sv)) = s.evaluate_1st_derivatives(p2(uv), false, false) else { return gvec(rng, size); }; @@ -555,7 +555,7 @@ fn test_on_normal(rng: &mut Rng, s: &Surface, uv: [f64; 2], size: f64) -> [f64; } /// A closed B-spline: a full rational circle (or a closed polygon) in u swept in v. -fn gen_closed(rng: &mut Rng) -> InlineSrf { +pub(crate) fn gen_closed(rng: &mut Rng) -> InlineSrf { let r = rng.uniform(0.5, 5.0); let h = rng.uniform(0.5, 5.0); let v_rows = 2 + rng.below(2) as usize; @@ -623,7 +623,7 @@ fn gen_closed(rng: &mut Rng) -> InlineSrf { } /// A bilinear (2 x 2) surface, sometimes twisted. -fn gen_bilinear(rng: &mut Rng) -> InlineSrf { +pub(crate) fn gen_bilinear(rng: &mut Rng) -> InlineSrf { let tw = if rng.chance(0.5) { 0.0 } else { rng.uniform(-1.0, 1.0) }; let net = vec![ vec![[0.0, 0.0, 0.0, NOW], [0.0, 1.0, 0.0, NOW]], @@ -675,7 +675,7 @@ fn gen_zigzag(rng: &mut Rng) -> InlineSrf { } } -fn gen_surface(rng: &mut Rng) -> SurfIn { +pub(crate) fn gen_surface(rng: &mut Rng) -> SurfIn { let bss = |s: InlineSrf| BssIn { oob: false, form: default_form(), @@ -690,7 +690,7 @@ fn gen_surface(rng: &mut Rng) -> SurfIn { } /// The control-net counts of a surface input (a plane's net is 2 x 2). -fn net_counts(s: &SurfIn) -> Result<(usize, usize), HarnessError> { +pub(crate) fn net_counts(s: &SurfIn) -> Result<(usize, usize), HarnessError> { match s { SurfIn::Plane(_) => Ok((2, 2)), SurfIn::Bspline(b) | SurfIn::CurveBounded(b) => { @@ -845,7 +845,7 @@ fn gen_case(rng: &mut Rng) -> Result { /// The test points of `my_test_surface_closest_point` on `Surface.sms` /// (`tests/prog_test/src/smSurfaceTest.cpp:4117`). -const SURFACE_SMS_POINTS: [[f64; 3]; 15] = [ +pub(crate) const SURFACE_SMS_POINTS: [[f64; 3]; 15] = [ [-1.0616429349692744, -20.167126109553802, 20.020823739437560], [-1.1583434957685239, -19.953438722680168, 20.005921768970701], [-4.5941574701649142, -9.6528635542995485, 19.290949275641431], diff --git a/crates/difftest/src/units/surface_surface_point_solve_p2.rs b/crates/difftest/src/units/surface_surface_point_solve_p2.rs new file mode 100644 index 0000000000..a4a68d9c3f --- /dev/null +++ b/crates/difftest/src/units/surface_surface_point_solve_p2.rs @@ -0,0 +1,948 @@ +//! Runner for unit `surface.surface_point_solve_p2`: the global solver's guesses, +//! refinements, seam solutions and `Solve`, and the local stepper's step adjustment, +//! norms, convergence and boundary tests, steps and line searches +//! (`SmSurfacePointSolve.cpp`). +//! +//! C++: `solidcore_ffi::units::surface_surface_point_solve_p2` (the shim compiles the +//! upstream file, whose classes are file-local). Rust: +//! `solidcore::surface::surface_point_solve::p2`. +//! +//! A case is a problem as in the `surface.surface_point_solve_p1` runner (surface `surf`, +//! wrapped in a STEP surface when `wrap`; solver domain `dom`, operation `solver` as an +//! `SmSolverOperationType` code, test point `test`, tolerance `tol`) and an `op` with its +//! parameters. Both sides produce the flat record documented in the ffi module; the +//! runner reports its first entry as the `status` (`0` = `SM_SUCCESS`, else the +//! `SmStatus`; `-1` for a Rust error upstream cannot produce) and the rest as the `out` +//! values. Solutions are `n` then `u v value objects variables object0_is_surface` each. +//! +//! Cases: inline cases on the upstream `Surface.sms` with the test points of +//! `my_test_surface_closest_point` (`tests/prog_test/src/smSurfaceTest.cpp:4117`), and +//! generated cases (planes, B-splines incl. closed, rational, bilinear, zigzag and +//! huge / tiny ones, curve-bounded and STEP-wrapped surfaces; sub-domains; test points on, +//! near and far from the surface; stepper points inside, on and near the boundary). +//! +//! Deferred cases: `SmSurface::LocalPointSolve` (`surface.surface_p5`; every `Solve` that +//! gets a rough guess), `SmSurface::IsOnSeam` (`surface.surface_p4`; `RefineGuess` on any +//! surface but a plane), `SmSurfaceEval::Normal` (`surface.surface_p4`; convergence at a +//! singularity), `SmSurface::InvertPointNearSingularity` (`surface.surface_p8`; a step +//! from a singularity), and the bounding box of a curve-bounded or STEP surface for the +//! directed operation (`surface.surface_p6`, as in p1). + +use serde::{Deserialize, Serialize}; +use serde_json::Value; + +use solidcore::geom::extent2d::Extent2d; +use solidcore::geom::global_solver::SolutionRequested; +use solidcore::geom::solution_array::{ObjectKind, Solution, SolutionArray, SolutionObject}; +use solidcore::kernel::Kernel; +use solidcore::surface::surface::Surface; +use solidcore::surface::surface_point_solve::p1::ErrorNormType; +use solidcore::surface::surface_point_solve::{SurfPtGlobalSolver, SurfPtSolveStepper}; +use solidcore_ffi::units::surface_surface_point_solve_p2::{self as ffi, Op, MAX_GRID, MAX_HITS}; + +use super::surface_surface_cache_p1::{dom_of, flag, gen_sub_dom, harness, srf, status, v3, BssIn, SrfIn}; +use super::surface_surface_cache_p3::{push_surf, surf_domain, surf_size, surface_in, SurfIn}; +use super::surface_surface_point_solve_p1::{ + gen_surface, gen_test, gen_uv, gen_uv_any, net_counts, p2, push2, solver_op, test_on_normal, + wrapped_of, SOLVERS, SURFACE_SMS_POINTS, +}; +use crate::case::{CaseContext, CaseSource}; +use crate::runner::default_cases; +use crate::{Case, HarnessError, Output, Rng, UnitDiff}; + +/// Generated cases per run, whatever `--seeds` says (it can only raise this). +pub const MIN_GENERATED: usize = 900; + +/// `surface.surface_point_solve_p2` runner. +pub struct Runner; + +// ---- inputs ----------------------------------------------------------------------- + +#[derive(Debug, Clone, Serialize, Deserialize)] +#[serde(tag = "op", rename_all = "snake_case")] +enum OpIn { + Rough, + Poly { + iu: usize, + iv: usize, + }, + Refine { + uv: [f64; 2], + }, + Patch { + sub: [f64; 4], + nu: i32, + nv: i32, + uv: [f64; 2], + best: f64, + }, + Grid { + uv: [f64; 2], + iu: usize, + iv: usize, + }, + Seams { + uv: [f64; 2], + value: f64, + }, + Solve, + Adjust { + uv0: [f64; 2], + tight: f64, + uv1: [f64; 2], + hu: i32, + hv: i32, + }, + Norms { + uv0: [f64; 2], + tight: f64, + }, + Conv { + uv0: [f64; 2], + tight: f64, + ctol: f64, + }, + Push { + uv0: [f64; 2], + tight: f64, + ptol: f64, + }, + Snap { + uv0: [f64; 2], + tight: f64, + stol: f64, + }, + Step { + uv0: [f64; 2], + tight: f64, + #[serde(default)] + plane: bool, + }, + Line { + uv0: [f64; 2], + tight: f64, + uv1: [f64; 2], + }, + LineB { + uv0: [f64; 2], + tight: f64, + uv1: [f64; 2], + norm: u8, + }, +} + +impl OpIn { + fn ffi(&self) -> Op { + match self { + OpIn::Rough => Op::Rough, + OpIn::Poly { .. } => Op::Poly, + OpIn::Refine { .. } => Op::Refine, + OpIn::Patch { .. } => Op::Patch, + OpIn::Grid { .. } => Op::Grid, + OpIn::Seams { .. } => Op::Seams, + OpIn::Solve => Op::Solve, + OpIn::Adjust { .. } => Op::Adjust, + OpIn::Norms { .. } => Op::Norms, + OpIn::Conv { .. } => Op::Conv, + OpIn::Push { .. } => Op::Push, + OpIn::Snap { .. } => Op::Snap, + OpIn::Step { .. } => Op::Step, + OpIn::Line { .. } => Op::Line, + OpIn::LineB { .. } => Op::LineB, + } + } +} + +#[derive(Debug, Clone, Serialize, Deserialize)] +struct Inputs { + /// Wrap the surface in an `SmSTEPSurface`. + #[serde(default)] + wrap: bool, + surf: SurfIn, + /// The solver domain. + dom: [f64; 4], + /// The `SmSolverOperationType` code. + solver: u8, + test: [f64; 3], + tol: f64, + #[serde(flatten)] + op: OpIn, +} + +// ---- C++ side ----------------------------------------------------------------------- + +fn cpp_input(i: &Inputs) -> Result<(Op, Vec), HarnessError> { + let mut o = vec![flag(i.wrap)]; + push_surf(&mut o, &i.surf)?; + o.extend(i.dom); + o.push(f64::from(i.solver)); + o.extend(i.test); + o.push(i.tol); + match &i.op { + OpIn::Rough | OpIn::Solve => {} + OpIn::Poly { iu, iv } => o.extend([*iu as f64, *iv as f64]), + OpIn::Refine { uv } => o.extend(uv), + OpIn::Patch { + sub, + nu, + nv, + uv, + best, + } => { + o.extend(sub); + o.extend([f64::from(*nu), f64::from(*nv)]); + o.extend(uv); + o.push(*best); + } + OpIn::Grid { uv, iu, iv } => { + o.extend(uv); + o.extend([*iu as f64, *iv as f64]); + } + OpIn::Seams { uv, value } => { + o.extend(uv); + o.push(*value); + } + OpIn::Adjust { + uv0, + tight, + uv1, + hu, + hv, + } => { + o.extend(uv0); + o.push(*tight); + o.extend(uv1); + o.extend([f64::from(*hu), f64::from(*hv)]); + } + OpIn::Norms { uv0, tight } => { + o.extend(uv0); + o.push(*tight); + } + OpIn::Conv { + uv0, + tight, + ctol: t, + } + | OpIn::Push { + uv0, + tight, + ptol: t, + } + | OpIn::Snap { + uv0, + tight, + stol: t, + } => { + o.extend(uv0); + o.extend([*tight, *t]); + } + OpIn::Step { uv0, tight, plane } => { + o.extend(uv0); + o.extend([*tight, flag(*plane)]); + } + OpIn::Line { uv0, tight, uv1 } => { + o.extend(uv0); + o.push(*tight); + o.extend(uv1); + } + OpIn::LineB { + uv0, + tight, + uv1, + norm, + } => { + o.extend(uv0); + o.push(*tight); + o.extend(uv1); + o.push(f64::from(*norm)); + } + } + Ok((i.op.ffi(), o)) +} + +fn to_output(flat: &[f64]) -> Output { + let (status, rest) = flat + .split_first() + .map_or((f64::NAN, &[][..]), |(s, r)| (*s, r)); + Output::new() + .values("status", vec![status]) + .values("out", rest.to_vec()) +} + +fn cpp(i: &Inputs) -> Result { + let (op, input) = cpp_input(i)?; + let flat = ffi::eval(op, &input).map_err(|e| harness(e.to_string()))?; + Ok(to_output(&flat)) +} + +// ---- Rust side ---------------------------------------------------------------------- + +/// The solutions' handle for the case's surface: it reports the surface's natural +/// domain (`GetNaturalUVDomain` of the solution object, for `IsEqual`). +#[derive(Debug, Clone, Copy, PartialEq)] +struct SrfObj(Extent2d); + +impl SolutionObject for SrfObj { + fn kind(&self) -> ObjectKind { + ObjectKind::Surface + } + + fn natural_uv_domain(&self) -> Option { + Some(self.0) + } +} + +fn push_solutions(o: &mut Vec, sols: &SolutionArray, obj: SrfObj) { + o.push(sols.len() as f64); + for s in sols.as_slice() { + o.extend([ + s.start.parameters[0], + s.start.parameters[1], + s.start.solution_value, + s.num_objects as f64, + s.num_variables as f64, + flag(s.objects[0] == Some(obj)), + ]); + } +} + +fn norm_of(code: u8) -> ErrorNormType { + match code { + 0 => ErrorNormType::SumAbs, + 1 => ErrorNormType::Sum, + 2 => ErrorNormType::F, + 3 => ErrorNormType::G, + 4 => ErrorNormType::FAbs, + 5 => ErrorNormType::GAbs, + 6 => ErrorNormType::Max, + 7 => ErrorNormType::MaxAbs, + _ => ErrorNormType::Gap, + } +} + +const ALL_NORMS: [ErrorNormType; 9] = [ + ErrorNormType::SumAbs, + ErrorNormType::Sum, + ErrorNormType::F, + ErrorNormType::G, + ErrorNormType::FAbs, + ErrorNormType::GAbs, + ErrorNormType::Max, + ErrorNormType::MaxAbs, + ErrorNormType::Gap, +]; + +/// `UV BoundHitsU BoundHitsV NRStepSuspect IterCount` of a stepper. +fn push_state(o: &mut Vec, s: &SurfPtSolveStepper<'_>) { + push2(o, s.uv()); + o.extend([ + f64::from(s.bound_hits_u()), + f64::from(s.bound_hits_v()), + flag(s.nr_step_suspect()), + f64::from(s.iter_count()), + ]); +} + +/// `HitBound(which, sign)` `|n|` times. +fn hits(s: &mut SurfPtSolveStepper<'_>, which: solidcore::nurbs::b_spline_surface::SurfParamType, n: i32) { + let sign = if n < 0 { -1 } else { 1 }; + for _ in 0..n.unsigned_abs() { + s.hit_bound(which, sign); + } +} + +/// An error status record (`-1` for an error upstream cannot produce), or the deferral. +fn err_record(e: &solidcore::SolidError) -> Result, solidcore::SolidError> { + if e.is_unported() { + return Err(e.clone()); + } + Ok(vec![status(e)]) +} + +#[allow(clippy::too_many_lines)] +fn rust_op(i: &Inputs) -> Result>, HarnessError> { + use solidcore::nurbs::b_spline_surface::SurfParamType; + let k = Kernel::default(); + let surface: Surface = wrapped_of(i.wrap, &i.surf)?; + let dom = dom_of(i.dom); + let op = solver_op(i.solver)?; + let test = v3(i.test); + let tol = i.tol; + let obj = SrfObj(surface.natural_uv_domain()); + let global = || { + SurfPtGlobalSolver::new( + &k, + &surface, + &dom, + op, + test, + tol, + None, + SolutionRequested::Single, + ) + }; + let stepper = |uv: [f64; 2], tight: f64| { + SurfPtSolveStepper::new(&k, &surface, &dom, test, p2(uv), tight, tol, op) + }; + let run = || -> solidcore::Result> { + let mut o = Vec::new(); + match &i.op { + OpIn::Rough => { + let mut g = global()?; + match g.find_rough_guess() { + Ok(r) => { + o.push(0.0); + push2(&mut o, r.uv); + o.extend([r.best_idx_u as f64, r.best_idx_v as f64]); + } + Err(e) => return err_record(&e), + } + } + OpIn::Poly { iu, iv } => { + let g = global()?; + match g.find_polygon_guess((*iu, *iv)) { + Ok(r) => { + o.push(0.0); + push2(&mut o, r.uv); + o.extend([r.best_idx_u as f64, r.best_idx_v as f64, flag(r.suspicious)]); + } + Err(e) => return err_record(&e), + } + } + OpIn::Refine { uv } => { + let mut g = global()?; + let mut uv = p2(*uv); + match g.refine_guess(&mut uv) { + Ok(()) => o.push(0.0), + Err(e) => o = err_record(&e)?, + } + push2(&mut o, uv); + } + OpIn::Patch { + sub, + nu, + nv, + uv, + best, + } => { + let g = global()?; + let mut uv = p2(*uv); + let mut best = *best; + g.refine_patch_by_grid(&dom_of(*sub), *nu, *nv, &mut uv, &mut best)?; + o.push(0.0); + push2(&mut o, uv); + o.push(best); + } + OpIn::Grid { uv, iu, iv } => { + let g = global()?; + let r = g.refine_guess_by_grid(p2(*uv), *iu, *iv)?; + o.push(0.0); + push2(&mut o, r); + } + OpIn::Seams { uv, value } => { + let mut g = global()?; + let mut sol: Solution = Solution::new(); + sol.num_variables = 2; + sol.num_objects = 1; + sol.objects[0] = Some(obj); + sol.start.parameters[0] = uv[0]; + sol.start.parameters[1] = uv[1]; + sol.start.solution_value = *value; + let tight = g.base.tight_tol(); + let sols = g.check_seams(&k, tight, &sol, Some(obj))?; + o.push(0.0); + push_solutions(&mut o, &sols, obj); + } + OpIn::Solve => { + let mut g = global()?; + match g.solve(&k, Some(obj)) { + Ok(sols) => { + o.push(0.0); + push_solutions(&mut o, &sols, obj); + } + Err(e) => return err_record(&e), + } + } + OpIn::Adjust { + uv0, + tight, + uv1, + hu, + hv, + } => { + let mut prev = stepper(*uv0, *tight)?; + hits(&mut prev, SurfParamType::U, *hu); + hits(&mut prev, SurfParamType::V, *hv); + let mut next = SurfPtSolveStepper::copy_of(&prev); + next.set_uv(p2(*uv1)); + match next.adjust_step(&prev) { + Ok((adjusted, scale)) => o.extend([0.0, flag(adjusted), scale]), + Err(e) => { + o = err_record(&e)?; + // upstream's outputs on SM_ERR + o.extend([flag(true), 0.0]); + } + } + push_state(&mut o, &next); + } + OpIn::Norms { uv0, tight } => { + let mut s = stepper(*uv0, *tight)?; + o.push(0.0); + for n in ALL_NORMS { + o.push(s.error_norm(n)?); + } + let (b, which) = s.is_singular()?; + o.extend([flag(b), f64::from(which as u8)]); + } + OpIn::Conv { uv0, tight, ctol } => { + let mut s = stepper(*uv0, *tight)?; + o.push(0.0); + o.push(flag(s.converged(*ctol)?)); + o.push(flag(s.converged_with(true)?)); + o.push(flag(s.converged_with(false)?)); + } + OpIn::Push { uv0, tight, ptol } => { + let mut s = stepper(*uv0, *tight)?; + let mut c = SurfPtSolveStepper::copy_of(&s); + o.push(0.0); + let (b, which) = s.pushing_boundary(*ptol)?; + o.extend([flag(b), f64::from(which as u8)]); + let (b, which) = c.converged_on_boundary(*ptol)?; + o.extend([flag(b), f64::from(which as u8)]); + } + OpIn::Snap { uv0, tight, stol } => { + let mut s = stepper(*uv0, *tight)?; + let b = s.snap_to_boundaries(*stol)?; + o.extend([0.0, flag(b)]); + push_state(&mut o, &s); + } + OpIn::Step { uv0, tight, plane } => { + let mut prev = stepper(*uv0, *tight)?; + let mut next = SurfPtSolveStepper::copy_of(&prev); + let r = next.step(&mut prev, *plane)?; + o.extend([0.0, f64::from(r.code())]); + push_state(&mut o, &next); + } + OpIn::Line { uv0, tight, uv1 } => { + let mut prev = stepper(*uv0, *tight)?; + let mut next = SurfPtSolveStepper::copy_of(&prev); + next.set_uv(p2(*uv1)); + match next.line_search(&mut prev) { + Ok(()) => o.push(0.0), + Err(e) => o = err_record(&e)?, + } + push_state(&mut o, &next); + } + OpIn::LineB { + uv0, + tight, + uv1, + norm, + } => { + let mut prev = stepper(*uv0, *tight)?; + let mut next = SurfPtSolveStepper::copy_of(&prev); + next.set_uv(p2(*uv1)); + match next.line_search_bracketted(&mut prev, norm_of(*norm)) { + Ok(()) => o.push(0.0), + Err(e) => o = err_record(&e)?, + } + push_state(&mut o, &next); + } + } + Ok(o) + }; + Ok(run()) +} + +fn rust(i: &Inputs) -> Result { + Ok(match rust_op(i)? { + Ok(o) => to_output(&o), + Err(e) => Output::unported(&e).unwrap_or_else(|| to_output(&[status(&e)])), + }) +} + +// ---- generated cases ---------------------------------------------------------------- + +/// A point on the boundary of `d` (a side or a corner), or just inside or outside it. +fn gen_uv_edge(rng: &mut Rng, d: [f64; 4]) -> [f64; 2] { + let [mut u, mut v] = gen_uv(rng, d); + let (w, h) = (d[2] - d[0], d[3] - d[1]); + match rng.below(4) { + 0 => u = if rng.chance(0.5) { d[0] } else { d[2] }, + 1 => v = if rng.chance(0.5) { d[1] } else { d[3] }, + 2 => { + // near a side, by a tiny to small amount + let e = 10f64.powf(rng.uniform(-12.0, -2.0)); + u = if rng.chance(0.5) { d[0] + e * w } else { d[2] - e * w }; + } + _ => { + let e = 10f64.powf(rng.uniform(-12.0, -2.0)); + v = if rng.chance(0.5) { d[1] + e * h } else { d[3] - e * h }; + } + } + [u, v] +} + +/// A point of `d` clipped into the surface domain `sd`. +fn clip(uv: [f64; 2], sd: [f64; 4]) -> [f64; 2] { + [uv[0].clamp(sd[0], sd[2]), uv[1].clamp(sd[1], sd[3])] +} + +#[allow(clippy::too_many_lines)] +fn gen_op(rng: &mut Rng, s: &SurfIn, d: [f64; 4], sd: [f64; 4], size: f64) -> Result { + let (nu, nv) = net_counts(s)?; + let idx = |rng: &mut Rng, n: usize| -> usize { + if rng.chance(0.3) { + if rng.chance(0.5) { + 0 + } else { + n.saturating_sub(1) + } + } else { + rng.below(n.max(1) as u64) as usize + } + }; + let small_tol = |rng: &mut Rng| -> f64 { + match rng.below(5) { + 0 => 1e-10, + 1 => size * 1e-6, + 2 => -size * 1e-6, + 3 => 0.0, + _ => size * rng.uniform(1e-4, 0.1), + } + }; + // stepper points: inside, on or near the boundary of the solver domain + let st_uv = |rng: &mut Rng| -> [f64; 2] { + if rng.chance(0.5) { + gen_uv(rng, sd) + } else { + clip(gen_uv_edge(rng, sd), d) + } + }; + let tight = |rng: &mut Rng| -> f64 { + match rng.below(3) { + 0 => 1e-12, + 1 => size * 1e-10, + _ => size * 1e-7, + } + }; + Ok(match rng.below(40) { + 0..=2 => OpIn::Rough, + 3..=5 => OpIn::Poly { + iu: idx(rng, nu), + iv: idx(rng, nv), + }, + 6..=7 => OpIn::Refine { + uv: clip(gen_uv_any(rng, sd), d), + }, + 8..=10 => { + let sub = if rng.chance(0.3) { + sd + } else { + gen_sub_dom(rng, d) + }; + let g = |rng: &mut Rng| -> i32 { + if rng.chance(0.2) { + rng.range_i64(-2, 2) as i32 + } else { + rng.range_i64(2, MAX_GRID.min(9)) as i32 + } + }; + OpIn::Patch { + sub, + nu: g(rng), + nv: g(rng), + uv: gen_uv(rng, sd), + best: match rng.below(3) { + 0 => 1e300, + 1 => size * rng.uniform(0.0, 2.0), + _ => -size, + }, + } + } + 11..=13 => OpIn::Grid { + uv: gen_uv(rng, sd), + iu: idx(rng, nu), + iv: idx(rng, nv), + }, + 14..=17 => OpIn::Seams { + uv: if rng.chance(0.7) { + clip(gen_uv_edge(rng, sd), d) + } else { + gen_uv(rng, sd) + }, + value: size * rng.uniform(0.0, 1.0), + }, + 18 => OpIn::Solve, + 19..=21 => { + let uv0 = st_uv(rng); + let uv1 = match rng.below(4) { + 0 => uv0, + 1 => [uv0[0] + 1e-13, uv0[1]], + _ => gen_uv_any(rng, sd), + }; + OpIn::Adjust { + uv0, + tight: tight(rng), + uv1, + hu: rng.range_i64(-MAX_HITS, MAX_HITS) as i32, + hv: rng.range_i64(-MAX_HITS, MAX_HITS) as i32, + } + } + 22..=24 => OpIn::Norms { + uv0: st_uv(rng), + tight: tight(rng), + }, + 25..=26 => OpIn::Conv { + uv0: st_uv(rng), + tight: tight(rng), + ctol: small_tol(rng).abs().max(1e-12), + }, + 27..=29 => OpIn::Push { + uv0: st_uv(rng), + tight: tight(rng), + ptol: small_tol(rng), + }, + 30..=31 => OpIn::Snap { + uv0: st_uv(rng), + tight: tight(rng), + stol: { + let w = (sd[2] - sd[0]).min(sd[3] - sd[1]); + w * 10f64.powf(rng.uniform(-10.0, -1.0)) + }, + }, + 32..=35 => OpIn::Step { + uv0: st_uv(rng), + tight: tight(rng), + plane: rng.chance(0.5), + }, + 36..=37 => OpIn::Line { + uv0: st_uv(rng), + tight: tight(rng), + uv1: gen_uv(rng, d), + }, + _ => OpIn::LineB { + uv0: st_uv(rng), + tight: tight(rng), + uv1: gen_uv(rng, d), + norm: rng.below(9) as u8, + }, + }) +} + +fn gen_case(rng: &mut Rng) -> Result { + let surf = gen_surface(rng); + let wrap = rng.chance(0.15); + let d = surf_domain(&surf)?; + let size = surf_size(&surf)?; + let sd = if rng.chance(0.5) { d } else { gen_sub_dom(rng, d) }; + let solver = *rng.choose(&SOLVERS).unwrap_or(&3); + let s = surface_in(&surf)?; + let test = if solver == 10 { + match rng.below(4) { + 0 => [1.0, 0.0, 0.0], + 1 => [0.0, -1.0, 0.0], + _ => super::surface_surface_cache_p1::gvec(rng, 1.0), + } + } else { + gen_test(rng, &s, sd, size) + }; + let tol = match rng.below(4) { + 0 => 1e-6, + 1 => size * 1e-8, + 2 => size * rng.uniform(1e-5, 1e-2), + _ => 1e-10, + }; + let op = gen_op(rng, &surf, d, sd, size)?; + let mut test = test; + // converged / boundary tests: test points on the normal at the stepper point + if solver != 10 && rng.chance(0.5) { + if let OpIn::Conv { uv0, .. } | OpIn::Push { uv0, .. } | OpIn::Seams { uv: uv0, .. } = &op { + test = test_on_normal(rng, &s, *uv0, size); + } + } + Ok(Inputs { + wrap, + surf, + dom: sd, + solver, + test, + tol, + op, + }) +} + +// ---- inline cases from the upstream tests ---------------------------------------------------- + +fn upstream_cases(ctx: &CaseContext, id: &str) -> Result, HarnessError> { + let path = ctx.upstream_test_file("Surface.sms")?; + let sms = SrfIn::File { + sms: path.to_string_lossy().into_owned(), + }; + let inline = srf(&sms)?; + let surf = SurfIn::Bspline(BssIn { + oob: false, + form: super::surface_surface_cache_p1::default_form(), + srf: sms, + }); + let d = [ + inline.u.first().copied().unwrap_or(0.0), + inline.v.first().copied().unwrap_or(0.0), + inline.u.last().copied().unwrap_or(1.0), + inline.v.last().copied().unwrap_or(1.0), + ]; + let (nu, nv) = (inline.net.len(), inline.net.first().map_or(0, Vec::len)); + let mid = [0.5 * (d[0] + d[2]), 0.5 * (d[1] + d[3])]; + let at = |fu: f64, fv: f64| [d[0] + (d[2] - d[0]) * fu, d[1] + (d[3] - d[1]) * fv]; + let mut out = Vec::new(); + for (pi, test) in SURFACE_SMS_POINTS.iter().enumerate() { + let f = pi as f64 / 14.0; + let ops = [ + OpIn::Rough, + OpIn::Poly { iu: 0, iv: 0 }, + OpIn::Grid { + uv: at(f, 1.0 - f), + iu: pi % nu, + iv: (pi * 3) % nv, + }, + OpIn::Patch { + sub: d, + nu: 7, + nv: 7, + uv: mid, + best: 1e300, + }, + OpIn::Seams { + uv: at(0.0, f), + value: 0.0, + }, + OpIn::Adjust { + uv0: mid, + tight: 1e-12, + uv1: at(1.2, f), + hu: 0, + hv: 0, + }, + OpIn::Norms { + uv0: at(f, 0.5), + tight: 1e-12, + }, + OpIn::Conv { + uv0: at(f, f), + tight: 1e-12, + ctol: 1e-6, + }, + OpIn::Push { + uv0: at(1.0, f), + tight: 1e-12, + ptol: 1e-9, + }, + OpIn::Step { + uv0: mid, + tight: 1e-12, + plane: false, + }, + OpIn::Step { + uv0: at(f, 0.25), + tight: 1e-12, + plane: true, + }, + OpIn::Line { + uv0: at(0.1, 0.1), + tight: 1e-12, + uv1: at(0.9, f), + }, + ]; + for (oi, op) in ops.into_iter().enumerate() { + let inputs = Inputs { + wrap: false, + surf: surf.clone(), + dom: d, + solver: 3, + test: *test, + tol: 1e-6, + op, + }; + out.push(Case { + unit: id.to_string(), + name: format!("surface_sms_pt{:02}_op{oi}", pi + 1), + source: CaseSource::Inline, + inputs: serde_json::to_value(inputs).map_err(|e| harness(e.to_string()))?, + }); + } + } + Ok(out) +} + +impl UnitDiff for Runner { + fn id(&self) -> &'static str { + "surface.surface_point_solve_p2" + } + + fn generate(&self, rng: &mut Rng) -> Option { + for _ in 0..64 { + let Ok(i) = gen_case(rng) else { + continue; + }; + let Ok((op, input)) = cpp_input(&i) else { + continue; + }; + if ffi::validate(op, &input).is_ok() { + return serde_json::to_value(i).ok(); + } + } + None + } + + fn cases(&self, ctx: &CaseContext) -> Result, HarnessError> { + let mut v = upstream_cases(ctx, self.id())?; + let ctx = CaseContext { + seeds: ctx.seeds.max(MIN_GENERATED), + ..ctx.clone() + }; + v.extend(default_cases(self, &ctx)?); + Ok(v) + } + + fn run_cpp(&self, case: &Case) -> Result { + cpp(&case.parse()?) + } + + fn run_rust(&self, case: &Case) -> Result { + rust(&case.parse()?) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use super::super::surface_surface_cache_p1::{default_form, InlineSrf}; + use super::super::surface_surface_point_solve_p1::gen_bilinear; + + #[test] + fn op_codes_round_trip() { + let bil: InlineSrf = gen_bilinear(&mut Rng::new(1)); + let i = Inputs { + wrap: false, + surf: SurfIn::Bspline(BssIn { + oob: false, + form: default_form(), + srf: SrfIn::Inline(bil), + }), + dom: [0.0, 0.0, 1.0, 1.0], + solver: 3, + test: [0.5, 0.5, 1.0], + tol: 1e-6, + op: OpIn::Patch { + sub: [0.0, 0.0, 1.0, 1.0], + nu: MAX_GRID as i32, + nv: 1, + uv: [0.5, 0.5], + best: 1e300, + }, + }; + let (op, input) = cpp_input(&i).unwrap(); + assert_eq!(op, Op::Patch); + assert!(ffi::validate(op, &input).is_ok()); + } +} diff --git a/crates/solidcore-ffi/src/units/surface_surface_point_solve_p1.rs b/crates/solidcore-ffi/src/units/surface_surface_point_solve_p1.rs index da7e355fcd..406bde585d 100644 --- a/crates/solidcore-ffi/src/units/surface_surface_point_solve_p1.rs +++ b/crates/solidcore-ffi/src/units/surface_surface_point_solve_p1.rs @@ -108,7 +108,7 @@ impl Op { pub const MAX_HIT_STEPS: usize = 64; /// The solver operations the file's classes accept. -const OPS: [usize; 5] = [0, 3, 4, 5, 10]; +pub(crate) const OPS: [usize; 5] = [0, 3, 4, 5, 10]; mod sys { extern "C" { @@ -124,13 +124,13 @@ mod sys { } /// The problem `P`: the surface domain and the control-net counts (`2 x 2` for a plane). -struct Problem { - srf: Srf, - num_u: usize, - num_v: usize, +pub(crate) struct Problem { + pub(crate) srf: Srf, + pub(crate) num_u: usize, + pub(crate) num_v: usize, } -fn problem(rd: &mut Reader<'_>) -> Result { +pub(crate) fn problem(rd: &mut Reader<'_>) -> Result { rd.flag()?; let at = rd.at; let (kind, srf, _) = sized_surface(rd)?; @@ -151,7 +151,7 @@ fn problem(rd: &mut Reader<'_>) -> Result { } /// A point `uv` inside the surface domain. -fn uv_in(rd: &mut Reader<'_>, srf: &Srf) -> Result<(), FfiError> { +pub(crate) fn uv_in(rd: &mut Reader<'_>, srf: &Srf) -> Result<(), FfiError> { let u = rd.finite()?; let v = rd.finite()?; let d = srf.dom; @@ -162,7 +162,7 @@ fn uv_in(rd: &mut Reader<'_>, srf: &Srf) -> Result<(), FfiError> { } /// A point `uv` within a domain size of the surface domain. -fn uv_near(rd: &mut Reader<'_>, srf: &Srf) -> Result<(), FfiError> { +pub(crate) fn uv_near(rd: &mut Reader<'_>, srf: &Srf) -> Result<(), FfiError> { let u = rd.finite()?; let v = rd.finite()?; let d = srf.dom; @@ -174,7 +174,7 @@ fn uv_near(rd: &mut Reader<'_>, srf: &Srf) -> Result<(), FfiError> { } /// Control-net indices `i j` with `i + du < num_u`, `j + dv < num_v`. -fn net_index(rd: &mut Reader<'_>, p: &Problem, d: usize) -> Result<(), FfiError> { +pub(crate) fn net_index(rd: &mut Reader<'_>, p: &Problem, d: usize) -> Result<(), FfiError> { let i = rd.index(1 << 20)?; let j = rd.index(1 << 20)?; if i + d >= p.num_u || j + d >= p.num_v { diff --git a/crates/solidcore-ffi/src/units/surface_surface_point_solve_p2.rs b/crates/solidcore-ffi/src/units/surface_surface_point_solve_p2.rs new file mode 100644 index 0000000000..038f4c2bee --- /dev/null +++ b/crates/solidcore-ffi/src/units/surface_surface_point_solve_p2.rs @@ -0,0 +1,219 @@ +//! FFI wrappers for unit `surface.surface_point_solve_p2`. +//! +//! Shim: `harness/cpp/shim/surface_surface_point_solve_p2.cpp` +//! (`sc_surface_surface_point_solve_p2_eval`, which compiles the upstream file with its +//! file-local classes renamed). Upstream: `source/SMLib/src/SmSurfacePointSolve.cpp`. +//! +//! Notation as in [`super::surface_surface_point_solve_p1`]: every op starts with the +//! problem `P = W S D op test(3) tol`; `uv` is a point `u v`, `D` a domain +//! `umin vmin umax vmax`. The stepper ops continue with `uv0 tight`: the stepper +//! `SmSurfPtSolveStepper(S, D, test, uv0, tight, tol, op)`. +//! +//! | [`Op`] | upstream | parameters | +//! |------------|----------------------------------------------------------|------------------------| +//! | `Rough` | `SmSurfPtGlobalSolver::FindRoughGuess` | `P` | +//! | `Poly` | `FindPolygonGuess` (from best indices `iu iv`) | `P iu iv` | +//! | `Refine` | `RefineGuess` | `P uv` | +//! | `Patch` | `RefinePatchByGrid` | `P sub nu nv uv best` | +//! | `Grid` | `RefineGuessByGrid` | `P uv iu iv` | +//! | `Seams` | `CheckSeams` of the solution `(u, v)` with `value` | `P uv value` | +//! | `Solve` | `Solve` | `P` | +//! | `Adjust` | `AdjustStep` of a copy moved to `uv1` (after `hu`/`hv` boundary hits) | `P uv0 tight uv1 hu hv` | +//! | `Norms` | the nine `ErrorNorm`s, `IsSingular` | `P uv0 tight` | +//! | `Conv` | `Converged(ctol)`, `Converged(TRUE)`, `Converged(FALSE)` | `P uv0 tight ctol` | +//! | `Push` | `PushingBoundary(ptol)`, `ConvergedOnBoundary(ptol)` | `P uv0 tight ptol` | +//! | `Snap` | `SnapToBoundaries(stol)` | `P uv0 tight stol` | +//! | `Step` | `Step` of a copy, Newton or plane | `P uv0 tight plane` | +//! | `Line` | `LineSearch` of a copy moved to `uv1` | `P uv0 tight uv1` | +//! | `LineB` | `LineSearchBracketted` of a copy moved to `uv1` | `P uv0 tight uv1 norm` | +//! +//! The output is the op status (`0` = `SM_SUCCESS`, else the `SmStatus`) and the op's +//! outputs (see the shim). [`eval`] accepts only what upstream handles without undefined +//! behaviour: the problems of [`super::surface_surface_point_solve_p1`], points inside the +//! surface domain wherever a surface evaluator reads them (the step target of `Adjust`, +//! which is only clamped, within a domain size of it), a grid sub-domain inside the +//! surface domain with grid sizes up to [`MAX_GRID`], control-net indices inside the net, +//! boundary-hit counts up to 3, an `SmSrfPtErrorNormType`, finite values. + +use super::surface_surface_cache_p1::{bad, Reader}; +use super::surface_surface_point_solve_p1::{net_index, problem, uv_in, uv_near}; +use crate::{call_with_buffer, FfiError}; + +/// The upstream routine run by the shim. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum Op { + /// `FindRoughGuess`. + Rough, + /// `FindPolygonGuess`. + Poly, + /// `RefineGuess`. + Refine, + /// `RefinePatchByGrid`. + Patch, + /// `RefineGuessByGrid`. + Grid, + /// `CheckSeams`. + Seams, + /// `Solve`. + Solve, + /// `AdjustStep`. + Adjust, + /// `ErrorNorm`, `IsSingular`. + Norms, + /// `Converged`. + Conv, + /// `PushingBoundary`, `ConvergedOnBoundary`. + Push, + /// `SnapToBoundaries`. + Snap, + /// `Step`. + Step, + /// `LineSearch`. + Line, + /// `LineSearchBracketted`. + LineB, +} + +impl Op { + /// The shim's op code. + #[must_use] + pub fn code(self) -> i32 { + match self { + Op::Rough => 0, + Op::Poly => 1, + Op::Refine => 2, + Op::Patch => 3, + Op::Grid => 4, + Op::Seams => 5, + Op::Solve => 6, + Op::Adjust => 7, + Op::Norms => 8, + Op::Conv => 9, + Op::Push => 10, + Op::Snap => 11, + Op::Step => 12, + Op::Line => 13, + Op::LineB => 14, + } + } +} + +/// Largest grid size of `Patch` (sizes below 2 count as 2 upstream). +pub const MAX_GRID: i64 = 24; +/// Largest boundary-hit count of `Adjust` (either sign). +pub const MAX_HITS: i64 = 3; +/// Largest `SmSrfPtErrorNormType` (`SM_SPN_GAP`). +pub const MAX_NORM: usize = 8; + +mod sys { + extern "C" { + pub fn sc_surface_surface_point_solve_p2_eval( + op: i32, + input: *const f64, + n_in: i64, + out: *mut f64, + cap: i64, + out_len: *mut i64, + ) -> i32; + } +} + +/// An integer in `-max..=max`. +fn small_int(rd: &mut Reader<'_>, max: i64, what: &str) -> Result<(), FfiError> { + let v = rd.num()?; + if v.fract() != 0.0 || v.abs() > max as f64 { + return Err(bad(format!("{what}: bad integer {v}"))); + } + Ok(()) +} + +fn validate_op(op: Op, rd: &mut Reader<'_>) -> Result<(), FfiError> { + let p = problem(rd)?; + match op { + Op::Rough | Op::Solve => {} + Op::Poly => net_index(rd, &p, 0)?, + Op::Refine => uv_in(rd, &p.srf)?, + Op::Patch => { + rd.domain_in(&p.srf, false)?; + small_int(rd, MAX_GRID, "grid size")?; + small_int(rd, MAX_GRID, "grid size")?; + rd.finite()?; + rd.finite()?; + if rd.num()?.is_nan() { + return Err(bad("best value is NaN")); + } + } + Op::Grid => { + uv_in(rd, &p.srf)?; + net_index(rd, &p, 0)?; + } + Op::Seams => { + uv_in(rd, &p.srf)?; + rd.finite()?; + } + Op::Adjust | Op::Norms | Op::Conv | Op::Push | Op::Snap | Op::Step | Op::Line | Op::LineB => { + uv_in(rd, &p.srf)?; + rd.finite()?; + match op { + Op::Adjust => { + uv_near(rd, &p.srf)?; + small_int(rd, MAX_HITS, "boundary hits")?; + small_int(rd, MAX_HITS, "boundary hits")?; + } + Op::Conv | Op::Push | Op::Snap => { + rd.finite()?; + } + Op::Step => { + rd.flag()?; + } + Op::Line => uv_in(rd, &p.srf)?, + Op::LineB => { + uv_in(rd, &p.srf)?; + rd.index(MAX_NORM)?; + } + _ => {} + } + } + } + Ok(()) +} + +/// Checks `input` for `op` (layouts in the module docs). +/// +/// # Errors +/// +/// [`FfiError::InvalidArgument`] for input the shim must not see. +pub fn validate(op: Op, input: &[f64]) -> Result<(), FfiError> { + let mut rd = Reader { input, at: 0 }; + validate_op(op, &mut rd)?; + rd.done() +} + +/// Runs `op` on `input` (layouts in the module docs); returns the shim's output record. +/// +/// # Errors +/// +/// [`FfiError::InvalidArgument`] for input the shim must not see; [`FfiError::Status`] +/// for a shim failure. +pub fn eval(op: Op, input: &[f64]) -> Result, FfiError> { + validate(op, input)?; + call_with_buffer(4096, |buf, len| { + // SAFETY: input is valid for input.len() reads and holds exactly what the shim + // reads for `op` (validated above: well-formed clamped surfaces, planes with + // independent axes, a solver domain inside the surface domain, supported + // operations, evaluated points and grid sub-domains inside the surface domain, + // control-net indices inside the net, bounded grid sizes and hit counts, a valid + // norm type); buf is valid for buf.len() writes and the shim honours the + // capacity; len is valid. + unsafe { + sys::sc_surface_surface_point_solve_p2_eval( + op.code(), + input.as_ptr(), + input.len() as i64, + buf.as_mut_ptr(), + buf.len() as i64, + &raw mut *len, + ) + } + }) +} diff --git a/crates/solidcore/src/surface/surface_point_solve.rs b/crates/solidcore/src/surface/surface_point_solve.rs index 9df1464395..5b2f89fdad 100644 --- a/crates/solidcore/src/surface/surface_point_solve.rs +++ b/crates/solidcore/src/surface/surface_point_solve.rs @@ -8,5 +8,6 @@ //! `_p3`; each lives in its own submodule and adds methods to the classes of `p1`. pub mod p1; +pub mod p2; pub use p1::{SurfPtGlobalSolver, SurfPtSolveStepper, SurfacePointSolver}; diff --git a/crates/solidcore/src/surface/surface_point_solve/p1.rs b/crates/solidcore/src/surface/surface_point_solve/p1.rs index bef7ede086..ba678fa069 100644 --- a/crates/solidcore/src/surface/surface_point_solve/p1.rs +++ b/crates/solidcore/src/surface/surface_point_solve/p1.rs @@ -150,14 +150,7 @@ impl<'s> SurfacePointSolver<'s> { _solution_requested: SolutionRequested, ) -> Result { // We need an SmBSplineSurface for some operations. - // (An SmOffsetSurface would work from its base surface: not ported.) - let b_spline = match surface.kind() { - SurfaceKind::BSpline(b) => Some(b), - SurfaceKind::Plane(p) => Some(p.base()), - SurfaceKind::SurfOfRevolution(s) => Some(s.base()), - SurfaceKind::CurveBounded(_) | SurfaceKind::Step(_) => None, - } - .filter(|b| b.gw_nurb().is_some()); + let b_spline = b_spline_part(surface); let mut solver = Self { surface, @@ -255,6 +248,21 @@ impl<'s> SurfacePointSolver<'s> { self.b_spline } + /// The same problem over `surface` (the surface data re-read from it; every other + /// field copied), for a solve on a copy of the surface. + pub(super) fn with_surface<'t>(&self, surface: &'t Surface) -> SurfacePointSolver<'t> { + SurfacePointSolver { + surface, + domain: self.domain, + test_pt: self.test_pt, + dir_vec: self.dir_vec, + solver_op: self.solver_op, + b_spline: b_spline_part(surface), + tight_tol: self.tight_tol, + loose_tol: self.loose_tol, + } + } + /// `crGap` reduced to its component along the direction for the directed /// operation (`dDot * m_sDirVec`), else unchanged. fn directed(&self, gap: Vec3) -> Vec3 { @@ -267,11 +275,23 @@ impl<'s> SurfacePointSolver<'s> { } } +/// `m_pBSplSrf`: the surface as an `SmBSplineSurface` with a NURB, if it is one (an +/// `SmOffsetSurface` would work from its base surface: not ported). +fn b_spline_part(surface: &Surface) -> Option<&BSplineSurface> { + match surface.kind() { + SurfaceKind::BSpline(b) => Some(b), + SurfaceKind::Plane(p) => Some(p.base()), + SurfaceKind::SurfOfRevolution(s) => Some(s.base()), + SurfaceKind::CurveBounded(_) | SurfaceKind::Step(_) => None, + } + .filter(|b| b.gw_nurb().is_some()) +} + // ---- control-net helpers ------------------------------------------------------------------ /// The control point `(i, j)` of a net stored row by row with `num_v` points per row, /// divided by its weight when the net is rational. -fn control_point( +pub(super) fn control_point( pts: &[CPoint], num_v: usize, i: usize, @@ -498,6 +518,22 @@ impl<'s> SurfPtGlobalSolver<'s> { }) } + /// The same solver over `surface` (see [`SurfacePointSolver::with_surface`]), + /// keeping the closure computed so far. + pub(super) fn with_surface<'t>(&self, surface: &'t Surface) -> SurfPtGlobalSolver<'t> { + SurfPtGlobalSolver { + base: self.base.with_surface(surface), + closed_in_u: self.closed_in_u, + closed_in_v: self.closed_in_v, + } + } + + /// Takes over the closure `other` has computed (the same surface geometry). + pub(super) fn copy_closure_from(&mut self, other: &SurfPtGlobalSolver<'_>) { + self.closed_in_u = other.closed_in_u; + self.closed_in_v = other.closed_in_v; + } + /// Whether the surface is closed over the domain in `which_dir` (`SM_SP_U` or /// `SM_SP_V`; `false` for any other value), computed once (`ClosedSurface`). /// diff --git a/crates/solidcore/src/surface/surface_point_solve/p2.rs b/crates/solidcore/src/surface/surface_point_solve/p2.rs new file mode 100644 index 0000000000..3a1376a352 --- /dev/null +++ b/crates/solidcore/src/surface/surface_point_solve/p2.rs @@ -0,0 +1,2221 @@ +//! Surface-point solving, part 2: the global solver's rough guess (control net and +//! grid), its refinements (tangent-plane steps, grids), the seam solutions and the +//! main `Solve`; the local stepper's step adjustment, error norms, convergence and +//! singularity tests, Newton / plane step, line searches and boundary tests. +//! +//! Upstream: `source/SMLib/src/SmSurfacePointSolve.cpp`. +//! +//! | Rust | Upstream | +//! |---------------------------------------------------|--------------------------------------------------| +//! | [`SurfPtGlobalSolver::find_rough_guess`], [`RoughGuess`] | `SmSurfPtGlobalSolver::FindRoughGuess` | +//! | [`SurfPtGlobalSolver::find_polygon_guess`], [`PolygonGuess`] | `SmSurfPtGlobalSolver::FindPolygonGuess` | +//! | [`SurfPtGlobalSolver::refine_guess`] | `SmSurfPtGlobalSolver::RefineGuess` | +//! | [`SurfPtGlobalSolver::refine_patch_by_grid`] | `SmSurfPtGlobalSolver::RefinePatchByGrid` | +//! | [`SurfPtGlobalSolver::refine_guess_by_grid`] | `SmSurfPtGlobalSolver::RefineGuessByGrid` | +//! | [`SurfPtGlobalSolver::check_seams`] | `SmSurfPtGlobalSolver::CheckSeams` | +//! | [`SurfPtGlobalSolver::solve`] | `SmSurfPtGlobalSolver::Solve` | +//! | [`SurfPtSolveStepper::copy_of`] (in `p1`) | `SmSurfPtSolveStepper` copy constructor | +//! | [`SurfPtSolveStepper::adjust_step`] | `SmSurfPtSolveStepper::AdjustStep` | +//! | [`SurfPtSolveStepper::error_norm`] | `SmSurfPtSolveStepper::ErrorNorm` | +//! | [`SurfPtSolveStepper::converged`], [`SurfPtSolveStepper::converged_with`] | `Converged(double)`, `Converged(SmBoolean)` | +//! | [`SurfPtSolveStepper::is_singular`] | `SmSurfPtSolveStepper::IsSingular` | +//! | [`SurfPtSolveStepper::step`] | `SmSurfPtSolveStepper::Step` | +//! | [`SurfPtSolveStepper::line_search`] | `SmSurfPtSolveStepper::LineSearch` | +//! | [`SurfPtSolveStepper::line_search_bracketted`] | `SmSurfPtSolveStepper::LineSearchBracketted` | +//! | [`SurfPtSolveStepper::converged_on_boundary`] | `SmSurfPtSolveStepper::ConvergedOnBoundary` | +//! | [`SurfPtSolveStepper::pushing_boundary`] | `SmSurfPtSolveStepper::PushingBoundary` | +//! | [`SurfPtSolveStepper::snap_to_boundaries`] | `SmSurfPtSolveStepper::SnapToBoundaries` | +//! +//! Deviations: +//! - **Out-params** become return values ([`RoughGuess`], [`PolygonGuess`], +//! `(adjusted, scale)` of the step adjustment, `(bool, SurfParamType)` of the +//! boundary and singularity tests, the solution array). The in/out guesses of +//! `RefineGuess` and `RefinePatchByGrid` (and the latter's best value) stay `&mut` +//! arguments: upstream changes them before some of its error returns. The unused +//! index outputs of `RefineGuess` are dropped. +//! - **Status returns.** `AdjustStep`'s `SM_ERR` (no step and outside the domain) +//! is [`SolidError`] `SM_ERR`; upstream's outputs are then `adjusted = TRUE`, +//! `scale = 0`, and the stepper is clamped to the domain. A boundary test that is +//! not on a boundary reports [`SurfParamType::Neither`] where upstream leaves the +//! caller's variable unset. +//! - **Ignored statuses** keep upstream's values: an evaluation that fails leaves +//! the point undefined ([`Point3::UNDEF`], upstream uninitialised) or at its previous +//! value; `smgu_VecLinCombTwoVectors` writes `(0, 0)` before failing. +//! [`SolidError::Unported`] always propagates. +//! - **Out-of-bounds evaluation.** `Solve` enables out-of-bounds evaluation of an +//! `SmBSplineSurface` (a B-spline, plane or surface of revolution) for the duration +//! of the solve (`SmTemporaryChangeValue`). The solver borrows its surface, so the +//! port solves on a copy with the flag set; solutions carry the caller's object +//! handle, so the copy is not visible. +//! - **Solutions' objects** are the caller's handle for the surface (`O`, upstream +//! the `SmSurface*`). +//! - **Unported callees.** `SmSurface::LocalPointSolve` (`surface.surface_p5`) for +//! every surface, `SmSurface::InvertPointNearSingularity` (`surface.surface_p8`) +//! and `SmSurfaceEval::Normal` (`surface.surface_p4`, through `p1`) return +//! [`SolidError::Unported`]. +//! - The stepper's `IsLocalMin` / `IsLocalMax` are not used upstream and are not +//! ported. Debug-only drawing and `SM_ASSERT`s are omitted. + +use crate::error::{Result, SolidError}; +use crate::geom::extent2d::Extent2d; +use crate::geom::geom_utility::p1::{line_closest_point, vec_lin_comb_two_vectors}; +use crate::geom::global_solver::SolverOperation; +use crate::geom::local_solve_nd::TerminationReason; +use crate::geom::solution_array::{Solution, SolutionArray, SolutionObject, SortKey}; +use crate::geom::tol::is_zero; +use crate::geom::vector2d::{Point2, Vec2}; +use crate::geom::vector3d::Point3; +use crate::kernel::Kernel; +use crate::math::{min2, BIG_DOUBLE, EFF_ZERO, EFF_ZERO_PARAM, EFF_ZERO_SQ, EFF_ZERO_SQRT}; +use crate::nurbs::b_spline_surface::SurfParamType; +use crate::surface::surface::p2::SurfaceEval; +use crate::surface::surface::{Surface, SurfaceKind}; + +use super::p1::{control_point, ErrorNormType, SurfPtGlobalSolver, SurfPtSolveStepper}; + +const OP_ROUGH: &str = "surface::surface_point_solve::p2::SurfPtGlobalSolver::find_rough_guess"; +const OP_POLY: &str = "surface::surface_point_solve::p2::SurfPtGlobalSolver::find_polygon_guess"; +const OP_ADJUST: &str = "surface::surface_point_solve::p2::SurfPtSolveStepper::adjust_step"; +const OP_LINE_B: &str = + "surface::surface_point_solve::p2::SurfPtSolveStepper::line_search_bracketted"; + +/// The result of [`SurfPtGlobalSolver::find_rough_guess`]. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct RoughGuess { + /// The UV guess (`rUVGuess`). + pub uv: Point2, + /// The best control point's u index (`rlBestIdxU`; `0` without a control net). + pub best_idx_u: usize, + /// The best control point's v index (`rlBestIdxV`). + pub best_idx_v: usize, +} + +/// The result of [`SurfPtGlobalSolver::find_polygon_guess`]. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct PolygonGuess { + /// The UV guess (`rUVGuess`; undefined when the control net could not be read + /// and the guess is not suspicious, as upstream). + pub uv: Point2, + /// The best control point's u index (`rlBestIdxU`). + pub best_idx_u: usize, + /// The best control point's v index (`rlBestIdxV`). + pub best_idx_v: usize, + /// The control net may be far from the surface (`rbSuspicious`). + pub suspicious: bool, +} + +impl SurfPtGlobalSolver<'_> { + /// A first rough guess (`FindRoughGuess`): the control-net guess, and when that + /// fails or is suspicious, the best of a 7 x 7, 5 x 5 and 3 x 3 grid search over + /// ever smaller sub-domains around the best grid point (fewer points for a net + /// with fewer control points); with both, the better of the two. + /// + /// # Errors + /// + /// `SM_ERR` when neither guess has a defined objective value; those of the + /// objective function and the grid search. + pub fn find_rough_guess(&mut self) -> Result { + // init output + let mut best_idx_u = 0; + let mut best_idx_v = 0; + let mut uv = Point2::UNDEF; + + // locals + let mut suspicious = false; + + // Find a guess from a BSplineSurface's control polygon + let poly_ok = match self.find_polygon_guess((best_idx_u, best_idx_v)) { + Ok(g) => { + uv = g.uv; + best_idx_u = g.best_idx_u; + best_idx_v = g.best_idx_v; + suspicious = g.suspicious; + true + } + Err(e) if e.is_unported() => return Err(e), + Err(_) => false, + }; + + // when FindPolygonGuess() worked without setting the bSuspicious flag + if poly_ok && !suspicious { + return Ok(RoughGuess { + uv, + best_idx_u, + best_idx_v, + }); + } + + // Next: Find the best GuessUV out of a Grid of UV samples. + // Grid 7x7, then 5x5, then 3x3. [080718; B46; B51; B394] + let mut grid_size_u: i32 = 7; // 5 is too small [B394] + let mut grid_size_v: i32 = 7; + let num_iters = 3; + + // If it's a single linear patch, don't subdivide. + // Similar for just a few control points. + let surface = self.base.surface(); + let num_cpts = + i32::try_from(surface.number_control_points(SurfParamType::U)).unwrap_or(i32::MAX); + if num_cpts < grid_size_u && num_cpts > 0 { + grid_size_u = num_cpts; + } + let num_cpts = + i32::try_from(surface.number_control_points(SurfParamType::V)).unwrap_or(i32::MAX); + if num_cpts < grid_size_v && num_cpts > 0 { + grid_size_v = num_cpts; + } + + let domain = *self.base.domain(); + let mut domain_size = domain.size(); + + let mut sub_domain = domain; + let mut grid_uv_guess = Point2::UNDEF; + let mut grid_ok = false; + + let mut best_func_val = BIG_DOUBLE; + + // for a sequence of ever smaller sub-domains centered on the last iterations + // best guess, search a grid of UV samples for the best UV guess point + for _ in 0..num_iters { + self.refine_patch_by_grid( + &sub_domain, + grid_size_u, + grid_size_v, + &mut grid_uv_guess, + &mut best_func_val, + )?; + grid_ok = true; + if grid_size_u < 2 || grid_size_v < 2 { + break; // Just to be safe. + } + + let fraction_u = 1.0 / f64::from(grid_size_u - 1) + 0.01; + let fraction_v = 1.0 / f64::from(grid_size_v - 1) + 0.01; + domain_size.x *= fraction_u; + domain_size.y *= fraction_v; + let uv_min = domain.clamp_point(grid_uv_guess - domain_size); + let uv_max = domain.clamp_point(grid_uv_guess + domain_size); + // (status ignored, as upstream) + let _ = sub_domain.set_min_max(uv_min, uv_max); + + // Decrease grid size for the next iteration: 7-5-3. + if grid_size_u > 4 { + grid_size_u -= 1; + } + if grid_size_v > 4 { + grid_size_v -= 1; + } + if grid_size_u > 2 { + grid_size_u -= 1; + } + if grid_size_v > 2 { + grid_size_v -= 1; + } + } + + // If guesses from both FindPolygonGuess() and RefinePatchByGrid(), pick the + // better one. + if poly_ok && grid_ok { + let poly_dist = self.calc_func_val(uv)?; + let grid_dist = self.calc_func_val(grid_uv_guess)?; + + // when either guess is uninitialized - quit + if poly_dist >= BIG_DOUBLE && grid_dist >= BIG_DOUBLE { + return Err(SolidError::sm_err(OP_ROUGH)); + } + + // when the GridGuess is better than the PolyGuess, use it + if grid_dist < poly_dist { + uv = grid_uv_guess; + } + } else if grid_ok { + uv = grid_uv_guess; + } else { + // both failed + return Err(SolidError::sm_err(OP_ROUGH)); + } + + Ok(RoughGuess { + uv, + best_idx_u, + best_idx_v, + }) + } + + /// A UV guess from the control net (`FindPolygonGuess`): the control-net + /// triangle guess, unless it is suspicious (a suspicious guess is accepted when + /// both the net and the surface come within 1% of the average control-point + /// spacing); else the closest control point (the extreme one for maximise and + /// the directed operation), refined along its polygon legs, trying every point of + /// a row or column of coincident points. `best_idx` is the caller's best control + /// point, kept when no net point qualifies (upstream's in/out indices). + /// + /// # Errors + /// + /// `SM_ERR` without a B-spline surface with a NURB; those of the seam check, the + /// leg refinement and the Greville abscissae. + pub fn find_polygon_guess(&self, best_idx: (usize, usize)) -> Result { + let b = &self.base; + // We need a SmBSplineSurface to work with control points. + let Some(bspl) = b.b_spline() else { + return Err(SolidError::sm_err(OP_POLY)); + }; + let op = b.solver_op(); + let test_pt = b.test_point(); + + let (mut best_idx_u, mut best_idx_v) = best_idx; + let mut uv = Point2::UNDEF; + + // find closest ControlPoint triangle with point closest to target + let suspicious = match self.find_polygon_triangle_guess((best_idx_u, best_idx_v)) { + Ok(t) => { + uv = t.uv; + best_idx_u = t.best_idx_u; + best_idx_v = t.best_idx_v; + let mut suspicious = t.suspicious; + + // When worried ControlPoint Mesh may be far from Surface + if suspicious && op != SolverOperation::Signed3dDirectedMaximize { + // when the BestDist is small compared to the Average spacing + if t.best_dist < t.ave_cpt_dist / 100.0 { + // Also check distance to the surface. [091002] + let srf_pt = eval_point(b.surface(), uv, Point3::UNDEF)?; + let surf_dist = srf_pt.distance(test_pt); + + if surf_dist < t.ave_cpt_dist / 100.0 { + suspicious = false; + } + } + } + suspicious + } + Err(e) if e.is_unported() => return Err(e), + // (the control net could not be read: only the rational warning is set) + Err(_) => bspl.is_rational(), + }; + + // If it looks like the UVGuess is good - return + if !suspicious { + return Ok(PolygonGuess { + uv, + best_idx_u, + best_idx_v, + suspicious, + }); + } + + // Next try: Iter every ControlPoint and find the one closest to the target + let rational = bspl.is_rational(); + let (num_u, num_v, cpts) = bspl + .control_points() + .map_err(|_| SolidError::sm_err(OP_POLY))?; + + let mut tied: Option<(usize, usize)> = None; + let mut num_ties: usize = 0; + let mut best_dist = BIG_DOUBLE; + if op == SolverOperation::Maximize { + best_dist = -best_dist; + } + + for ii in 0..num_u { + for jj in 0..num_v { + let this_pt = control_point(&cpts, num_v, ii, jj, rational)?; + + // the best dThisDist has the smallest signed value + let this_dist = if op == SolverOperation::Maximize { + -test_pt.distance(this_pt) + } else if op == SolverOperation::Signed3dDirectedMaximize { + this_pt.dot(b.dir_vec()) + } else { + test_pt.distance(this_pt) + }; + + // save best found case + if this_dist < best_dist - EFF_ZERO { + best_dist = this_dist; + best_idx_u = ii; + best_idx_v = jj; + // No tie; reset. + tied = None; + num_ties = 0; + } + // else if two equivalent choices are found + else if this_dist < best_dist + EFF_ZERO { + // 'Equal': record a tie. + num_ties += 1; + tied = Some((best_idx_u, best_idx_v)); + best_idx_u = ii; + best_idx_v = jj; + } + } + } + + // Have to check here for seams (status ignored). + let (iu, iv) = + self.check_guess_for_seams_cpts(best_idx_u, best_idx_v, num_u, num_v, &cpts, rational)?; + best_idx_u = iu; + best_idx_v = iv; + + // Refine along the legs of the control polygon. + let mut uv_guess_poly = Point2::UNDEF; + let mut leg_ok = match self.refine_guess_polygon_leg(best_idx_u, best_idx_v) { + Ok(g) => { + uv_guess_poly = g; + true + } + Err(e) if e.is_unported() => return Err(e), + Err(_) => false, + }; + + // Check for ties -- coincident control points. + let mut best_dist = BIG_DOUBLE; + let mut srf_pt = Point3::UNDEF; + let tied_u = tied.map(|t| t.0); + let tied_v = tied.map(|t| t.1); + + if num_ties >= 2 && num_ties == num_u.wrapping_sub(1) && tied_v == Some(best_idx_v) { + // Adjacent u-points with the same value (presumably coincident). + for ii in 0..num_u { + if let Some(guess) = self.leg_or_none(ii, best_idx_v)? { + srf_pt = eval_point(b.surface(), guess, srf_pt)?; + let this_dist = srf_pt.distance(test_pt); + if this_dist < best_dist { + best_dist = this_dist; + uv_guess_poly = guess; + best_idx_u = ii; + + // Quit if we get a direct hit. + if this_dist < b.loose_tol() { + break; + } + } + } + } + leg_ok = best_dist < BIG_DOUBLE / 2.0; + } else if num_ties >= 2 && num_ties == num_v.wrapping_sub(1) && tied_u == Some(best_idx_u) { + // Adjacent v-points. + for ii in 0..num_v { + if let Some(guess) = self.leg_or_none(best_idx_u, ii)? { + srf_pt = eval_point(b.surface(), guess, srf_pt)?; + let this_dist = srf_pt.distance(test_pt); + if this_dist < best_dist { + best_dist = this_dist; + uv_guess_poly = guess; + best_idx_v = ii; + + if this_dist < b.loose_tol() { + break; + } + } + } + } + leg_ok = best_dist < BIG_DOUBLE / 2.0; + } + + if leg_ok { + uv = uv_guess_poly; + } else { + // Use the Greville abscissae for the control point indices. + uv = Point2::new( + bspl.greville_abscissa(SurfParamType::U, best_idx_u)?, + bspl.greville_abscissa(SurfParamType::V, best_idx_v)?, + ); + } + + Ok(PolygonGuess { + uv: b.domain().clamp_point(uv), + best_idx_u, + best_idx_v, + suspicious, + }) + } + + /// `RefineGuessPolygonLeg` with its status: the guess, `None` on failure. + fn leg_or_none(&self, idx_u: usize, idx_v: usize) -> Result> { + match self.refine_guess_polygon_leg(idx_u, idx_v) { + Ok(g) => Ok(Some(g)), + Err(e) if e.is_unported() => Err(e), + Err(_) => Ok(None), + } + } + + /// Refines `uv_guess` to where Newton iteration should work (`RefineGuess`): up to + /// three tangent-plane steps (each at most 10% of the domain; for the directed + /// operation steps along the descent direction, from 5% of the domain halving), + /// with seam checks, stopping on a negligible objective value or step; the result + /// is the better of the last two points. + /// + /// On an evaluation error `uv_guess` holds the seam-checked input, as upstream. + /// + /// # Errors + /// + /// Those of the surface derivatives (upstream's `SER`), the objective function + /// and the seam checks. + pub fn refine_guess(&mut self, uv_guess: &mut Point2) -> Result<()> { + // Note: don't use a Newton step here, we're presumably not close enough. + // Try a couple of plane steps. + *uv_guess = self.seams_status_ignored(*uv_guess)?; // Since we do this on subsequent guesses... + + let mut this_uv = *uv_guess; + let mut next_uv = *uv_guess; + let mut this_f_val; + let mut next_f_val = BIG_DOUBLE; + + let domain = *self.base.domain(); + let u_extent = domain.u_interval().length(); + let v_extent = domain.v_interval().length(); + + // Don't ever take very big steps here. + let max_u_shift = u_extent * 0.10; + let max_v_shift = v_extent * 0.10; + + // For directed operations, we have to guess a step size. + // Start at 1/20 of the domain, and shrink each iteration. + let domain_size = domain.size().length(); + let mut step_size_guess = 0.05 * domain_size; + + let op = self.base.solver_op(); + + for _ in 0..3 { + this_f_val = next_f_val; + + next_f_val = self.calc_func_val(next_uv)?; + + if next_f_val.abs() < EFF_ZERO_SQRT { + *uv_guess = next_uv; + return Ok(()); // very close, definitely good enough guess. + } + + if next_f_val >= this_f_val { + // if no movement, then we're close enough. + if (next_uv.x - this_uv.x).abs() < u_extent / 1000.0 + && (next_uv.y - this_uv.y).abs() < v_extent / 1000.0 + { + break; + } + + if op == SolverOperation::Signed3dDirectedMaximize { + next_uv = (this_uv + next_uv) / 2.0; + continue; + } + break; // have to try something different. + } + + // Take a plane step. + let (pt, su, sv) = self + .base + .surface() + .evaluate_1st_derivatives(next_uv, false, false)?; + + let (u_shift, v_shift); + if op == SolverOperation::Signed3dDirectedMaximize { + // Plane step doesn't work here. All we can get is a direction. + let mut uv_step = vec_lin_comb_two_vectors(su, sv, -self.base.dir_vec()) + .unwrap_or(Vec2::new(0.0, 0.0)); + let mag = uv_step.length(); + if mag < EFF_ZERO_SQRT { + break; + } + + uv_step *= step_size_guess / mag; + step_size_guess /= 2.0; // shrink on each iteration. + + u_shift = uv_step.x; + v_shift = uv_step.y; + } else { + // 'Point' cases. + let gap = pt - self.base.test_point(); + + let f = gap.dot(su); // These are the right-hand side (b-vector) + let g = gap.dot(sv); + + // Solve 2x2 to get step + let fu = su.dot(su); + let fv = su.dot(sv); + let gu = fv; + let gv = sv.dot(sv); + + let det = fu * gv - fv * gu; + let mut us = -(f * gv - fv * g); // numerators + let mut vs = -(fu * g - f * gu); // to be divided by dDet. + + // For our purposes, if singular, just do what we can. + if is_zero(det) || is_zero(us) { + us = 0.0; + } else if det.abs() > us.abs() * EFF_ZERO_SQRT { + us /= det; + } else { + us = 0.0; + } + + if is_zero(det) || is_zero(vs) { + vs = 0.0; + } else if det.abs() > vs.abs() * EFF_ZERO_SQRT { + vs /= det; + } else { + vs = 0.0; + } + + // if no movement, then we're close enough. + if us.abs() < u_extent / 1000.0 && vs.abs() < v_extent / 1000.0 { + *uv_guess = next_uv; + return Ok(()); + } + + // Don't take big steps. + if us.abs() > max_u_shift { + us = if us > 0.0 { max_u_shift } else { -max_u_shift }; + } + if vs.abs() > max_v_shift { + vs = if vs > 0.0 { max_v_shift } else { -max_v_shift }; + } + u_shift = us; + v_shift = vs; + } + + // Set up for next iter. + this_uv = next_uv; + + next_uv = Point2::new(this_uv.x + u_shift, this_uv.y + v_shift); + + next_uv = domain.clamp_point(next_uv); + + next_uv = self.seams_status_ignored(next_uv)?; + } + + // We didn't return with success during the loop. + this_f_val = next_f_val; + next_f_val = self.calc_func_val(next_uv)?; + + *uv_guess = if next_f_val < this_f_val { + next_uv + } else { + this_uv + }; + Ok(()) + } + + /// `CheckGuessForSeams(rUVGuess)` with its status ignored: the input on failure. + fn seams_status_ignored(&mut self, uv: Point2) -> Result { + match self.check_guess_for_seams(uv) { + Ok(u) => Ok(u), + Err(e) if e.is_unported() => Err(e), + Err(_) => Ok(uv), + } + } + + /// Grids `sub_domain` with `num_u` x `num_v` points (at least 2 x 2) and keeps the + /// best point to the test point (`RefinePatchByGrid`), trying a plane step from + /// each grid point too ([`SurfPtGlobalSolver::check_bad_guess_stride`]). + /// `best_func_val` carries the best objective value over sequential calls: + /// `uv_guess` changes only when a point beats it. + /// + /// # Errors + /// + /// Those of the objective function and the stride check. + pub fn refine_patch_by_grid( + &self, + sub_domain: &Extent2d, + num_u: i32, + num_v: i32, + uv_guess: &mut Point2, + best_func_val: &mut f64, + ) -> Result<()> { + let surface = self.base.surface(); + let mut this_grid_pt = SurfaceEval::new(Some(surface)); + + let num_u = num_u.max(2); + let num_v = num_v.max(2); + + for ii in 0..num_u { + let frac_u = f64::from(ii) / f64::from(num_u - 1); + + for jj in 0..num_v { + let frac_v = f64::from(jj) / f64::from(num_v - 1); + + let this_uv = sub_domain.evaluate(frac_u, frac_v); + + this_grid_pt.set_uv(this_uv); + + let mut this_val = self.calc_func_val(this_uv)?; + + let mut best_uv = this_uv; // Best of: this uv, or one of the refined ones. + + // Check for too-big strides. [B46, B51] + if let Some((temp_srf_pt, temp_val)) = + self.check_bad_guess_stride(&mut this_grid_pt)? + { + if temp_val < this_val { + this_val = temp_val; + best_uv = temp_srf_pt.uv(); + } + } + + if this_val < *best_func_val { + *uv_guess = best_uv; + *best_func_val = this_val; + } + } + } + + Ok(()) + } + + /// A third refinement of `uv_guess` by brute-force gridding (`RefineGuessByGrid`): + /// 6 x 6 grids over the control-polygon neighbourhood of the best control point + /// `(best_idx_u, best_idx_v)` (a quarter of the domain each way around the guess + /// without a control net), then twice over a grid cell around the best point; + /// keeps the input when no grid point is better. + /// + /// # Errors + /// + /// Those of the objective function, the grid search and the Greville abscissae. + pub fn refine_guess_by_grid( + &self, + uv_guess: Point2, + best_idx_u: usize, + best_idx_v: usize, + ) -> Result { + let mut uv_guess = uv_guess; + let domain = *self.base.domain(); + + // Sometimes the original guess is better than any gridded value. + let mut best_uv_guess = uv_guess; + let mut best_val = self.calc_func_val(best_uv_guess)?; + + // Find a subdomain around rUVGuess. + let (mut min_uv, mut max_uv); + + if let Some(bspl) = self.base.b_spline() { + let num_u = bspl.number_control_points(SurfParamType::U); + let num_v = bspl.number_control_points(SurfParamType::V); + + let lo_u = if best_idx_u > 0 { best_idx_u - 1 } else { 0 }; + let hi_u = if best_idx_u < num_u.wrapping_sub(1) { + best_idx_u + 1 + } else { + num_u.wrapping_sub(1) + }; + let lo_v = if best_idx_v > 0 { best_idx_v - 1 } else { 0 }; + let hi_v = if best_idx_v < num_v.wrapping_sub(1) { + best_idx_v + 1 + } else { + num_v.wrapping_sub(1) + }; + + min_uv = Point2::new( + bspl.greville_abscissa(SurfParamType::U, lo_u)?, + bspl.greville_abscissa(SurfParamType::V, lo_v)?, + ); + max_uv = Point2::new( + bspl.greville_abscissa(SurfParamType::U, hi_u)?, + bspl.greville_abscissa(SurfParamType::V, hi_v)?, + ); + } else { + // No B-Spline surface. Use 1/4 of the domain. + let domain_size = domain.size(); + min_uv = uv_guess - domain_size / 4.0; + max_uv = uv_guess + domain_size / 4.0; + } + + min_uv = domain.clamp_point(min_uv); + max_uv = domain.clamp_point(max_uv); + + let mut sub_domain = Extent2d::from_points(min_uv, max_uv); + + let grid_size: i32 = 6; + + let mut this_val = BIG_DOUBLE; + + self.refine_patch_by_grid( + &sub_domain, + grid_size, + grid_size, + &mut uv_guess, + &mut this_val, + )?; + + // Check whether it improved. + if this_val < best_val { + best_val = this_val; + best_uv_guess = uv_guess; + } + + // Smaller grids, twice. + for last in [false, true] { + let delta_u = max_uv.x - min_uv.x; + let delta_v = max_uv.y - min_uv.y; + + min_uv = Point2::new( + uv_guess.x - delta_u / f64::from(grid_size - 1), + uv_guess.y - delta_v / f64::from(grid_size - 1), + ); + max_uv = Point2::new( + uv_guess.x + delta_u / f64::from(grid_size - 1), + uv_guess.y + delta_v / f64::from(grid_size - 1), + ); + + min_uv = domain.clamp_point(min_uv); + max_uv = domain.clamp_point(max_uv); + + // (status ignored, as upstream) + let _ = sub_domain.set_min_max(min_uv, max_uv); + + self.refine_patch_by_grid( + &sub_domain, + grid_size, + grid_size, + &mut uv_guess, + &mut this_val, + )?; + + if last { + // Check whether it improved. + if this_val > best_val { + uv_guess = best_uv_guess; + } + } else if this_val < best_val { + best_val = this_val; + best_uv_guess = uv_guess; + } + } + + Ok(uv_guess) + } + + /// The solutions for `sol` across the seams it is on (`CheckSeams`): `sol` (moved + /// onto a seam when that is no worse) and, on a closed surface, the point on the + /// other side of each seam (and of both, for a torus) when it is a solution too, + /// sorted by parameters. `surface_obj` is the caller's handle for the surface, + /// stored in the added solutions. `_tol` is not used (as upstream). + /// + /// # Errors + /// + /// Those of the boundary and closure tests, the evaluations and the sorted + /// insertion. + pub fn check_seams( + &mut self, + k: &Kernel, + _tol: f64, + sol: &Solution, + surface_obj: Option, + ) -> Result> { + let mut solutions = SolutionArray::new(); + solutions.add(sol); + + let mut uv = Point2::new(sol.start.param(0)?, sol.start.param(1)?); + + let b = &self.base; + let domain = *b.domain(); + let surface = b.surface(); + + // Check with a looser tol, because we will check whether the cross-seam + // point is a solution, and reject it if not. + let domain_tol = 0.001 * domain.size().length(); + let (on_boundary, which_bound, _, _) = + surface.is_on_boundary(k, &domain, uv, Some(domain_tol))?; + if !on_boundary { + return Ok(solutions); // can't be on a seam. + } + + let mut other_sol = *sol; + + // Make two steppers, to evaluate and compare. + let mut this_eval = SurfPtSolveStepper::new( + k, + surface, + &domain, + b.test_point(), + uv, + b.tight_tol(), + b.loose_tol(), + b.solver_op(), + )?; + let mut other_eval = SurfPtSolveStepper::copy_of(&this_eval); + let tight_tol = b.tight_tol(); + + // For comparison: + let norm_type = if b.solver_op() == SolverOperation::Normalize { + ErrorNormType::SumAbs + } else { + ErrorNormType::Gap + }; + + // Do U direction. + if self.closed_surface(SurfParamType::U)? + && (which_bound == SurfParamType::U || which_bound == SurfParamType::Both) + { + let dom = domain.u_interval(); + let (near_seam_u, far_seam_u) = if uv.x < dom.mid() { + (dom.min(), dom.max()) + } else { + (dom.max(), dom.min()) + }; + // If not quite on the seam, see whether right on is better. + if (uv.x - near_seam_u).abs() > tight_tol { + other_eval.set_uv(Point2::new(near_seam_u, uv.y)); + if other_eval.error_norm(norm_type)? <= this_eval.error_norm(norm_type)? { + uv.x = near_seam_u; + let gap_len = other_eval.gap_length()?; + if let Some(first) = solutions.get_mut(0) { + *first.start.param_mut(0)? = near_seam_u; + first.start.solution_value = gap_len; + } + } + } + + // Now check across the seam: is it a solution too? + other_eval.set_uv(Point2::new(far_seam_u, uv.y)); + if other_eval.error_norm(norm_type)? < this_eval.error_norm(norm_type)? + EFF_ZERO_SQRT + { + *other_sol.start.param_mut(0)? = far_seam_u; + *other_sol.start.param_mut(1)? = uv.y; + other_sol.start.solution_value = other_eval.gap_length()?; + other_sol.num_objects = 1; + set_first_object(&mut other_sol, surface_obj); + solutions.add_sorted_solution(&other_sol, SortKey::ByFirstTwoParameters)?; + } + } + + // Now do V direction. + if self.closed_surface(SurfParamType::V)? + && (which_bound == SurfParamType::V || which_bound == SurfParamType::Both) + { + let dom = domain.v_interval(); + let (near_seam_v, far_seam_v) = if uv.y < dom.mid() { + (dom.min(), dom.max()) + } else { + (dom.max(), dom.min()) + }; + // If not quite on the seam, see whether right on is better. + if (uv.y - near_seam_v).abs() > tight_tol { + other_eval.set_uv(Point2::new(uv.x, near_seam_v)); + if other_eval.error_norm(norm_type)? <= this_eval.error_norm(norm_type)? { + uv.y = near_seam_v; + let gap_len = other_eval.gap_length()?; + if let Some(first) = solutions.get_mut(0) { + *first.start.param_mut(1)? = near_seam_v; + first.start.solution_value = gap_len; + } + } + } + + // Now check across the seam: is it a solution too? + other_eval.set_uv(Point2::new(uv.x, far_seam_v)); + if other_eval.error_norm(norm_type)? < this_eval.error_norm(norm_type)? + EFF_ZERO_SQRT + { + *other_sol.start.param_mut(0)? = uv.x; + *other_sol.start.param_mut(1)? = far_seam_v; + other_sol.start.solution_value = other_eval.gap_length()?; + other_sol.num_objects = 1; + set_first_object(&mut other_sol, surface_obj); + solutions.add_sorted_solution(&other_sol, SortKey::ByFirstTwoParameters)?; + } + } + + // Finally, check both directions (torus). + if self.closed_surface(SurfParamType::U)? + && self.closed_surface(SurfParamType::V)? + && which_bound == SurfParamType::Both + { + let dom = domain.u_interval(); + let far_seam_u = if uv.x < dom.mid() { + dom.max() + } else { + dom.min() + }; + let dom = domain.v_interval(); + let far_seam_v = if uv.y < dom.mid() { + dom.max() + } else { + dom.min() + }; + + // Now check across both seams: is it a solution too? + other_eval.set_uv(Point2::new(far_seam_u, far_seam_v)); + if other_eval.error_norm(norm_type)? < this_eval.error_norm(norm_type)? + EFF_ZERO_SQRT + { + *other_sol.start.param_mut(0)? = far_seam_u; + *other_sol.start.param_mut(1)? = far_seam_v; + other_sol.start.solution_value = other_eval.gap_length()?; + solutions.add_sorted_solution(&other_sol, SortKey::ByFirstTwoParameters)?; + } + } + + Ok(solutions) + } + + /// The main global solve (`Solve`): a rough guess, refined and handed to the + /// local solver over each crease sub-domain (the whole domain: crease splitting is + /// off upstream); a result that is the wrong kind of extremum (or not on the + /// surface for intersect, or pushing a boundary for normalise) is retried from a + /// grid-refined guess; the best solution is returned with its across-seam + /// partners. Out-of-bounds evaluation is enabled for the solve. `surface_obj` is + /// the caller's handle for the surface. + /// + /// # Errors + /// + /// `SM_ERR` without a rough guess; [`SolidError::Unported`]: + /// `SmSurface::LocalPointSolve` (`surface.surface_p5`); those of the guesses and + /// the tests. + pub fn solve( + &mut self, + k: &Kernel, + surface_obj: Option, + ) -> Result> { + // We have to enable out-of-bounds evaluations. Otherwise, the evaluators + // will silently clamp, leading to wrong answers. + if let Some(oob_surface) = with_out_of_bounds_enabled(self.base.surface()) { + let mut solver = self.with_surface(&oob_surface); + let solutions = solver.solve_body(k, surface_obj); + // (the closure is a property of the geometry: keep what was computed) + self.copy_closure_from(&solver); + return solutions; + } + self.solve_body(k, surface_obj) + } + + /// The body of [`Self::solve`], once out-of-bounds evaluation is enabled. + fn solve_body( + &mut self, + k: &Kernel, + surface_obj: Option, + ) -> Result> { + let mut solutions = SolutionArray::new(); + + // First find a rough guess. + let mut sol: Solution = Solution::new(); + sol.num_variables = 2; + let mut found_answer = false; + + let rough = self.find_rough_guess()?; + let (best_idx_u, best_idx_v) = (rough.best_idx_u, rough.best_idx_v); + + // Check for creases. Note: this is currently not used. + let domain_list = self.crease_subdomains()?; + + let mut sol_array: SolutionArray = SolutionArray::new(); + let mut stat_ok = true; + let op = self.base.solver_op(); + let loose_tol = self.base.loose_tol(); + let surface = self.base.surface(); + + for this_domain in &domain_list { + // Use only those guesses that are in or near this subdomain. + if !this_domain.contains_point(rough.uv, EFF_ZERO_SQRT) { + continue; + } + + // Refine the guess (status ignored, as upstream). + let mut uv_guess = rough.uv; + if let Err(e) = self.refine_guess(&mut uv_guess) { + if e.is_unported() { + return Err(e); + } + } + + // We have to check here for seams: start on the proper side. + uv_guess = self.seams_status_ignored(uv_guess)?; + + // Do the local solve. + let pt_for_local = if op == SolverOperation::Signed3dDirectedMaximize { + -self.base.dir_vec() + } else { + self.base.test_point() + }; + match local_point_solve(surface, this_domain, op, pt_for_local, uv_guess, &sol) { + Ok((found, s)) => { + found_answer = found; + sol = s; + stat_ok = true; + } + Err(e) if e.is_unported() => return Err(e), + Err(_) => stat_ok = false, + } + + // Double-check this solution; see whether we should try it again with a + // more refined guess. + let mut try_again = false; + let mut surf_eval = + SurfaceEval::with_uv(Some(surface), sol.start.param(0)?, sol.start.param(1)?); + + match op { + SolverOperation::Minimize | SolverOperation::Signed3dDirectedMaximize => { + if found_answer && self.is_local_max(&mut surf_eval, loose_tol)? { + found_answer = false; + } + if !found_answer { + try_again = true; + } + } + SolverOperation::Maximize => { + if found_answer && self.is_local_min(&mut surf_eval, loose_tol)? { + found_answer = false; + } + if !found_answer { + try_again = true; + } + } + SolverOperation::Intersect => { + if found_answer && sol.start.solution_value > loose_tol { + found_answer = false; + } + if !found_answer + && !self.is_local_min(&mut surf_eval, loose_tol)? + && !self.pushing_boundary(&mut surf_eval, loose_tol)?.0 + { + try_again = true; + } + } + // Note: Even if it's pushing a boundary, try again anyway. + SolverOperation::Normalize if !found_answer => try_again = true, + _ => {} + } + + if try_again { + let save_guess = uv_guess; + uv_guess = self.refine_guess_by_grid(uv_guess, best_idx_u, best_idx_v)?; + + if !uv_guess.closer_than(EFF_ZERO_SQRT, save_guess) { + match local_point_solve(surface, this_domain, op, pt_for_local, uv_guess, &sol) + { + Ok((found, s)) => { + found_answer = found; + sol = s; + stat_ok = true; + if found_answer { + surf_eval + .set_uv(Point2::new(sol.start.param(0)?, sol.start.param(1)?)); + } + } + Err(e) if e.is_unported() => return Err(e), + Err(_) => stat_ok = false, + } + } + } + + // check the final result to see whether it's really a solution. + if found_answer { + match op { + SolverOperation::Minimize | SolverOperation::Signed3dDirectedMaximize => { + if self.is_local_max(&mut surf_eval, loose_tol)? { + found_answer = false; + } + } + SolverOperation::Maximize => { + if self.is_local_min(&mut surf_eval, loose_tol)? { + found_answer = false; + } + } + SolverOperation::Intersect => { + if sol.start.solution_value > loose_tol { + found_answer = false; + } + } + SolverOperation::Normalize + if self.pushing_boundary(&mut surf_eval, loose_tol)?.0 => + { + found_answer = false; + } + _ => {} + } + } + + if stat_ok && found_answer { + if sol.num_variables == 0 { + sol.num_variables = 2; + } + sol_array.add(&sol); + } + } + + // Sort out multiple solutions. + let mut best_sol_value = BIG_DOUBLE; + for s in sol_array.as_slice() { + if s.start.solution_value < best_sol_value { + sol = *s; + best_sol_value = sol.start.solution_value; + found_answer = true; + } + } + + // If on a seam, add across-seam solutions. + if stat_ok && found_answer { + // TightTol: [080520] + solutions = self.check_seams(k, self.base.tight_tol(), &sol, surface_obj)?; + } + + Ok(solutions) + } +} + +/// `sOtherSol.m_apObjects[0] = m_pSrf`. +fn set_first_object(sol: &mut Solution, surface_obj: Option) { + if let Some(slot) = sol.objects.first_mut() { + *slot = surface_obj; + } +} + +/// A copy of `surface` with out-of-bounds evaluation enabled, when it is an +/// `SmBSplineSurface` (B-spline, plane, surface of revolution) without it; `None` +/// when nothing changes (`SmTemporaryChangeValue` on `GetOutOfBoundsEnabled()`). +fn with_out_of_bounds_enabled(surface: &Surface) -> Option { + let enabled = match surface.kind() { + SurfaceKind::BSpline(b) => b.is_out_of_bounds_enabled(), + SurfaceKind::Plane(p) => p.base().is_out_of_bounds_enabled(), + SurfaceKind::SurfOfRevolution(s) => s.base().is_out_of_bounds_enabled(), + SurfaceKind::CurveBounded(_) | SurfaceKind::Step(_) => return None, + }; + if enabled { + return None; + } + let mut copy = surface.clone(); + match copy.kind_mut() { + SurfaceKind::BSpline(b) => *b.out_of_bounds_enabled_mut() = true, + SurfaceKind::Plane(p) => *p.out_of_bounds_enabled_mut() = true, + SurfaceKind::SurfOfRevolution(s) => *s.base_mut().out_of_bounds_enabled_mut() = true, + SurfaceKind::CurveBounded(_) | SurfaceKind::Step(_) => {} + } + Some(copy) +} + +/// `pSurf->EvaluatePoint(uv, rPt)` with its status ignored: `prev` on failure. +fn eval_point(surface: &Surface, uv: Point2, prev: Point3) -> Result { + match surface.evaluate_point(uv) { + Ok(p) => Ok(p), + Err(e) if e.is_unported() => Err(e), + Err(_) => Ok(prev), + } +} + +/// `m_pSrf->LocalPointSolve(domain, op, point, guess, bFoundAnswer, sSol)`: whether +/// an answer was found and the solution. +fn local_point_solve( + _surface: &Surface, + _domain: &Extent2d, + _op: SolverOperation, + _point: Point3, + _guess: Point2, + _sol: &Solution, +) -> Result<(bool, Solution)> { + // UNPORTED(surface.surface_p5): SmSurface::LocalPointSolve + Err(SolidError::unported( + "SmSurface::LocalPointSolve", + "surface.surface_p5", + )) +} + +/// `m_pSrf->InvertPointNearSingularity(testPt, domain, uv, tightTol, eReason)`: the +/// inverted UV. +fn invert_point_near_singularity( + _surface: &Surface, + _test_pt: Point3, + _domain: &Extent2d, + _uv: Point2, + _tol: f64, +) -> Result { + // UNPORTED(surface.surface_p8): SmSurface::InvertPointNearSingularity + Err(SolidError::unported( + "SmSurface::InvertPointNearSingularity", + "surface.surface_p8", + )) +} + +// ---- SmSurfPtSolveStepper ----------------------------------------------------------------- + +impl SurfPtSolveStepper<'_> { + /// Adjusts the step from `prev` to this stepper's UV (`AdjustStep`): limits each + /// parameter's step to [`SurfPtSolveStepper::max_step`] of the domain, clips at + /// the domain boundary (the better of the clamped and the proportionally scaled + /// step), records boundary hits, and moves there. Returns whether the step was + /// adjusted and the final step length as a fraction of the original. + /// + /// # Errors + /// + /// `SM_ERR` when there is no step and the UV is outside the domain (it is clamped + /// to the domain: upstream's outputs are `adjusted`, scale `0`); those of the + /// evaluations. + pub fn adjust_step(&mut self, prev: &Self) -> Result<(bool, f64)> { + let mut adjusted = false; + let mut scale = 1.0; + let domain = *self.base.domain(); + + let mut step = self.uv() - prev.uv(); + let orig_step_size = step.length(); + + if orig_step_size < EFF_ZERO_SQ { + if domain.contains_point(self.uv(), 0.0) { + return Ok((adjusted, scale)); + } + self.set_uv(domain.clamp_point(self.uv())); + return Err(SolidError::sm_err(OP_ADJUST)); + } + + // Never take too big a step. + // Scale the two parameters individually [ Fillet 2:234 ] + let mut size_scale = 1.0; + + let max_step = self.max_step(); + + let max_u_step = max_step * domain.u_interval().length(); + if step.x.abs() > max_u_step { + adjusted = true; + size_scale = max_u_step / step.x.abs(); + step.x *= size_scale; + } + + let max_v_step = max_step * domain.v_interval().length(); + if step.y.abs() > max_v_step { + adjusted = true; + let scale_y = min2(max_v_step / step.y.abs(), size_scale); + step.y *= scale_y; + size_scale = min2(size_scale, scale_y); + } + + // Apply the step and check out of bounds. + let curr_uv = prev.uv(); + let next_uv = curr_uv + step; + + let bc = Self::check_bounds(curr_uv, next_uv, &domain); + + self.hit_bound(SurfParamType::U, bc.hits_u); + self.hit_bound(SurfParamType::V, bc.hits_v); + + if bc.out_of_bounds { + let clamped_uv = bc.clamped_uv.unwrap_or(next_uv); + // Clamped. Use the better of scaled and unscaled. + if bc.step_scale > 1.0 - EFF_ZERO || bc.step_scale < EFF_ZERO { + // both are the same, or scaled would be zero step. + self.set_uv(clamped_uv); + } else { + // Evaluate both and use the better one. + self.set_uv(clamped_uv); + + let mut scaled_eval = Self::copy_of(self); + scaled_eval.set_uv(bc.clamped_uv_scaled.unwrap_or(next_uv)); + + if scaled_eval.error_norm(ErrorNormType::SumAbs)? + < self.error_norm(ErrorNormType::SumAbs)? + { + self.assign(&scaled_eval); + // (upstream then resets its local hit flags, which are not read + // again) + } + } + } else { + self.set_uv(next_uv); + } + + if size_scale < 1.0 || bc.out_of_bounds { + let final_step_size = (self.uv() - prev.uv()).length(); + scale = final_step_size / orig_step_size; + adjusted = true; + } + + Ok((adjusted, scale)) + } + + /// The measure of goodness `norm_type` of the evaluation (`ErrorNorm`): from the + /// unitised-derivative check values (`F / |Su|`, `G / |Sv|`), or the gap length. + /// + /// # Errors + /// + /// [`SolidError::Unported`] from the evaluation. + pub fn error_norm(&mut self, norm_type: ErrorNormType) -> Result { + let v = self.eval()?; // Make sure we're initialized. + + Ok(match norm_type { + ErrorNormType::SumAbs => v.sum_check, + ErrorNormType::Sum => v.f_check + v.g_check, + ErrorNormType::F => v.f_check, + ErrorNormType::G => v.g_check, + ErrorNormType::FAbs => v.f_check.abs(), + ErrorNormType::GAbs => v.g_check.abs(), + ErrorNormType::Max => { + if v.f_check.abs() >= v.g_check.abs() { + v.f_check + } else { + v.g_check + } + } + ErrorNormType::MaxAbs => { + if v.f_check.abs() >= v.g_check.abs() { + v.f_check.abs() + } else { + v.g_check.abs() + } + } + ErrorNormType::Gap => v.gap_len, + }) + } + + /// Whether the evaluation has converged to the tight (`use_tight_tol`) or loose + /// tolerance (`Converged(SmBoolean)`). + /// + /// # Errors + /// + /// Those of [`Self::converged`]. + pub fn converged_with(&mut self, use_tight_tol: bool) -> Result { + let tol = if use_tight_tol { + self.base.tight_tol() + } else { + self.base.loose_tol() + }; + self.converged(tol) + } + + /// Whether the evaluation has converged within `tol` (`Converged(double)`): on + /// the test point, or on the surface normal (by the unitised check values away + /// from a singularity; at one, the normal lined up with the direction, or passing + /// within `tol` (scaled by the distance) of the test point). + /// + /// # Errors + /// + /// [`SolidError::Unported`] from the evaluation and, at a singularity, + /// `SmSurfaceEval::Normal` (`surface.surface_p4`). + pub fn converged(&mut self, tol: f64) -> Result { + self.eval()?; + + // If surface point equals test point, we're done. + if self.error_norm(ErrorNormType::Gap)? < tol { + return Ok(true); + } + + // Check on the normal. Not reliable for singularities however. + if !self.is_singular()?.0 { + return Ok(self.error_norm(ErrorNormType::SumAbs)? < tol); + } + + // We're at a singularity. + let norm = self.normal()?; + + if norm.length_squared() < EFF_ZERO_SQ { + return Ok(false); + } + + // For this case, converged is when the normals line up. + if self.base.solver_op() == SolverOperation::Signed3dDirectedMaximize { + let dir = self.base.dir_vec(); + let diff = if norm.dot(dir) > 0.0 { + norm - dir + } else { + norm + dir + }; + return Ok(diff.length_squared() <= tol * tol); + } + + // All other cases: check distance from test pt to surface normal. + let test_pt = self.base.test_point(); + let pos = self.pos()?; + let line_t = line_closest_point(pos, norm, test_pt); + + let line_pnt = pos + line_t * norm; + let diff = line_pnt - test_pt; + + // Scale by how far the point is from the surface. + let mut dist_to_surf = line_pnt.distance(self.pos()?); + if dist_to_surf < 1.0 { + dist_to_surf = 1.0; + } + + Ok(diff.length() < tol * dist_to_surf) + } + + /// Whether the evaluation is near a singularity (`IsSingular`): one first + /// derivative much smaller (by a factor 10000) than the other and below + /// `SM_EFF_ZERO_PARAM` in squared length; which one ([`SurfParamType::U`], + /// [`SurfParamType::V`] or [`SurfParamType::Neither`]). + /// + /// # Errors + /// + /// [`SolidError::Unported`] from the evaluation. + pub fn is_singular(&mut self) -> Result<(bool, SurfParamType)> { + self.eval()?; + + let ratio_limit_squared = 10000.0 * 10000.0; + let abs_mag_limit = EFF_ZERO_PARAM; + + let mut local_which = SurfParamType::Neither; + + let su = self.su()?; + let sv = self.sv()?; + let mag_su_sq = su.dot(su); + let mag_sv_sq = sv.dot(sv); + + if mag_su_sq < mag_sv_sq / ratio_limit_squared && mag_su_sq < abs_mag_limit { + local_which = SurfParamType::U; + } + + if mag_sv_sq < mag_su_sq / ratio_limit_squared && mag_sv_sq < abs_mag_limit { + local_which = SurfParamType::V; + } + + Ok((local_which != SurfParamType::Neither, local_which)) + } + + /// Takes a Newton (or, with `do_plane_step`, tangent-plane) step from `prev` + /// and moves this stepper there (`Step`), using unitised first derivatives and + /// limiting the second-derivative terms on the concave side (marking the Newton + /// step suspect). For the directed operation a plane step only gets a direction + /// (a step of 5% of the domain). Near a singularity the point is inverted there. + /// Returns [`TerminationReason::StillIterating`] (a step was taken), + /// [`TerminationReason::FoundAnswerConverged`] (inverted at a singularity) or + /// [`TerminationReason::UnableToConverge`] (no step). + /// + /// # Errors + /// + /// [`SolidError::Unported`] from the evaluations and, near a singularity, + /// `SmSurface::InvertPointNearSingularity` (`surface.surface_p8`). + pub fn step(&mut self, prev: &mut Self, do_plane_step: bool) -> Result { + const SECOND_DERIV_TERM_LIMIT: f64 = 0.90; + + self.eval()?; + prev.eval()?; + + let op = self.base.solver_op(); + + // Special case: Plane step with directed operations: + // Can get direction but not step size. + if do_plane_step && op == SolverOperation::Signed3dDirectedMaximize { + let mut step = vec_lin_comb_two_vectors(prev.su()?, prev.sv()?, -self.base.dir_vec()) + .unwrap_or(Vec2::new(0.0, 0.0)); + let mag = step.length(); + if mag > EFF_ZERO_SQRT / 100.0 { + const STEP_SIZE_GUESS: f64 = 0.05; + let domain_size = self.base.domain().size().length(); + step *= STEP_SIZE_GUESS * domain_size / mag; + } + self.set_uv(prev.uv() + step); + + return Ok(TerminationReason::StillIterating); + } + + let mut got_step = false; + + // Check for near singularity: if one deriv is much smaller than the other. + // Unfortunately, this isn't implemented for signed-directed case. + if prev.is_singular()?.0 && op != SolverOperation::Signed3dDirectedMaximize { + // Near singularity, use the surface's own domain. [B50] + let surface = self.base.surface(); + let srf_domain = surface.natural_uv_domain(); + + match invert_point_near_singularity( + surface, + self.base.test_point(), + &srf_domain, + prev.uv(), + self.base.tight_tol(), + ) { + Ok(sing_uv) => { + self.set_uv(sing_uv); + return Ok(TerminationReason::FoundAnswerConverged); // end iteration. + } + Err(e) if e.is_unported() => return Err(e), + Err(_) => {} + } + } + + // We use unitized first derivatives. + let su = prev.su()?; + let sv = prev.sv()?; + let mu = su.length(); + let mv = sv.length(); + // We already checked for singularity, so this is unlikely. + if mu < EFF_ZERO || mv < EFF_ZERO { + return Ok(TerminationReason::UnableToConverge); + } + + let su_n = su / mu; + let sv_n = sv / mv; + + // Now the function values. + let pv = prev.eval()?; + let f = pv.gap.dot(su_n); + let g = pv.gap.dot(sv_n); + + // The four derivatives of the unitized first derivatives. + let suu = prev.suu()?; + let suv = prev.suv()?; + let svv = prev.svv()?; + let su_nu = su.unitized_derivative(suu); + let su_nv = su.unitized_derivative(suv); + let sv_nu = sv.unitized_derivative(suv); + let sv_nv = sv.unitized_derivative(svv); + + // Keep the two terms of Fu and Gv separate. + let fu1 = pv.gap_u.dot(su_n); + let mut fu2 = pv.gap.dot(su_nu); + + // Limit the magnitude of dFu2, if it's negative (see upstream's notes). + + if fu2 < -fu1 * SECOND_DERIV_TERM_LIMIT + && (op == SolverOperation::Minimize + || op == SolverOperation::Intersect + || op == SolverOperation::Signed3dDirectedMaximize) + { + self.set_nr_suspect(true); + fu2 = -fu1 * SECOND_DERIV_TERM_LIMIT; + } + + let gv1 = pv.gap_v.dot(sv_n); + let mut gv2 = pv.gap.dot(sv_nv); + + if gv2 < -gv1 * SECOND_DERIV_TERM_LIMIT + && (op == SolverOperation::Minimize || op == SolverOperation::Intersect) + { + self.set_nr_suspect(true); + gv2 = -gv1 * SECOND_DERIV_TERM_LIMIT; + } + + let (fu, gv, fv, gu) = if do_plane_step { + // planar step: ignore 2nd deriv terms. + (fu1, gv1, pv.gap_v.dot(su_n), pv.gap_u.dot(sv_n)) + } else { + // Newton step. + ( + fu1 + fu2, + gv1 + gv2, + pv.gap_v.dot(su_n) + pv.gap.dot(su_nv), + pv.gap_u.dot(sv_n) + pv.gap.dot(sv_nu), + ) + }; + + // 2x2 for the steps: Cramer's rule + let det = fu * gv - gu * fv; + let numer_u = -(f * gv - g * fv); + let numer_v = -(fu * g - gu * f); + + if det.abs() <= EFF_ZERO * numer_u.abs() || det.abs() <= EFF_ZERO * numer_v.abs() { + self.set_nr_suspect(true); + } else { + // Calculate the step. + got_step = true; + let step = Vec2::new(numer_u / det, numer_v / det); + self.set_uv(prev.uv() + step); + } + + Ok(if got_step { + TerminationReason::StillIterating + } else { + TerminationReason::UnableToConverge + }) + } + + /// Looks for a better point between `prev` and this stepper (`LineSearch`): a + /// bracketed secant search on `F` or `G` when one dominates and changes sign; + /// else samples at 10%, 20%, 40% and 70% of the way (returning the first that + /// halves the error), and beyond either end when the error looks monotone, and + /// moves to the best (clamped to the domain). + /// + /// # Errors + /// + /// Those of [`Self::line_search_bracketted`]; [`SolidError::Unported`] from the + /// evaluations. + pub fn line_search(&mut self, prev: &mut Self) -> Result<()> { + const VALS: [f64; 4] = [0.1, 0.2, 0.4, 0.7]; // was: { 0.5, 0.2, 0.8, 0.4, 0.6, 0.1, 0.9 } + + // If one function value (F or G) is much bigger than the other, and it + // changes sign, iterate on that one only. + if prev.error_norm(ErrorNormType::F)? * self.error_norm(ErrorNormType::F)? < 0.0 + && prev.error_norm(ErrorNormType::FAbs)? + > 10.0 * prev.error_norm(ErrorNormType::GAbs)? + && self.error_norm(ErrorNormType::FAbs)? + > 10.0 * self.error_norm(ErrorNormType::GAbs)? + { + return self.line_search_bracketted(prev, ErrorNormType::F); + } + if prev.error_norm(ErrorNormType::G)? * self.error_norm(ErrorNormType::G)? < 0.0 + && prev.error_norm(ErrorNormType::GAbs)? + > 10.0 * prev.error_norm(ErrorNormType::FAbs)? + && self.error_norm(ErrorNormType::GAbs)? + > 10.0 * self.error_norm(ErrorNormType::FAbs)? + { + return self.line_search_bracketted(prev, ErrorNormType::G); + } + + // Since we're likely in a problem area, use 3d gap for the error, + // unless we're Normalizing or doing a directed operation. + let op = self.base.solver_op(); + let norm_type = if op == SolverOperation::Normalize + || op == SolverOperation::Signed3dDirectedMaximize + { + ErrorNormType::SumAbs + } else { + ErrorNormType::Gap + }; + + let mut third_eval = Self::copy_of(self); // init to this: bnds counts, etc. + + // Set dBestErr and dBestVal. + let err0 = prev.error_norm(norm_type)?; + let err1 = self.error_norm(norm_type)?; + + let mut best_err = err0; + let mut best_val = 0.0; + if err1 < best_err { + best_err = err1; + best_val = 1.0; + } + let best_init_err = best_err; + + for fraction in VALS { + third_eval.set_uv((1.0 - fraction) * prev.uv() + fraction * self.uv()); + + let this_err = third_eval.error_norm(norm_type)?; + if this_err < best_err { + // If good enough, return it, else record it and keep looking. + if this_err < best_init_err * 0.50 { + self.assign(&third_eval); + return Ok(()); + } + + best_err = this_err; + best_val = fraction; + } + } + + // If no minimum between the two inputs, try outside the range. + if best_val < 0.01 && err0 < err1 { + // Appears to be monotonically increasing. + let fraction = -1.0; + third_eval.set_uv((1.0 - fraction) * prev.uv() + fraction * self.uv()); + // Need to clamp this. + third_eval.clamp(); + + let this_err = third_eval.error_norm(norm_type)?; + if this_err < best_err { + best_err = this_err; + best_val = fraction; + } + } + + if best_val > 0.99 && err1 < err0 { + // Appears to be monotonically decreasing. + let fraction = 2.0; + third_eval.set_uv((1.0 - fraction) * prev.uv() + fraction * self.uv()); + third_eval.clamp(); + + // (the best error is not read again) + if third_eval.error_norm(norm_type)? < best_err { + best_val = fraction; + } + } + + self.set_uv((1.0 - best_val) * prev.uv() + best_val * self.uv()); + + self.clamp(); // In case not between the two inputs. + + Ok(()) + } + + /// A secant search for the zero of the signed norm `norm_type` between `prev` + /// and this stepper, which bracket it (`LineSearchBracketted`): up to 5 secant + /// steps (clamped to the domain), extended by 5 at a time up to 30 with a + /// bisection when the value stops improving; moves to the last point. + /// + /// # Errors + /// + /// `SM_ERR` for a norm that cannot be negative, or when the two values do not + /// bracket zero; [`SolidError::Unported`] from the evaluations. + pub fn line_search_bracketted( + &mut self, + prev: &mut Self, + norm_type: ErrorNormType, + ) -> Result<()> { + if matches!( + norm_type, + ErrorNormType::SumAbs + | ErrorNormType::FAbs + | ErrorNormType::GAbs + | ErrorNormType::MaxAbs + | ErrorNormType::Gap + ) { + return Err(SolidError::sm_err(OP_LINE_B)); + } + + let mut y0 = prev.error_norm(norm_type)?; + let mut y1 = self.error_norm(norm_type)?; + + if y0 * y1 > 0.0 { + return Err(SolidError::sm_err(OP_LINE_B)); + } + + let best_input_y = min2(y0.abs(), y1.abs()); + + let mut max_iter = 5; // Don't spend too much time on a line search. + + // Call the initial points the unit interval: + let mut x0 = 0.0; + let mut x1 = 1.0; + + let mut mid_eval = Self::copy_of(prev); + let domain = *self.base.domain(); + + let mut ii = 0; + while ii < max_iter { + // The secant step: find the zero on the line between the evals. + let mut x_mid = x0 - y0 * (x1 - x0) / (y1 - y0); + + mid_eval.set_uv((1.0 - x_mid) * prev.uv() + x_mid * self.uv()); + + // Due to numerical noise the interpolation can leave the domain. So: + mid_eval.set_uv(domain.clamp_point(mid_eval.uv())); + + let mut y_mid = mid_eval.error_norm(norm_type)?; + + if y_mid.abs() < EFF_ZERO_SQRT { + self.assign(&mid_eval); + return Ok(()); + } + // Another check: the fcn value can get worse with this method. + // If it's not improving, try something different: a bisection. + if ii >= max_iter - 1 && y_mid.abs() >= best_input_y { + max_iter += 5; + if max_iter > 30 { + break; + } + x_mid = (x0 + x1) / 2.0; + + mid_eval.set_uv((1.0 - x_mid) * prev.uv() + x_mid * self.uv()); + y_mid = mid_eval.error_norm(norm_type)?; + } + + if y_mid * y0 < 0.0 { + x1 = x_mid; + y1 = y_mid; + } else { + x0 = x_mid; + y0 = y_mid; + } + ii += 1; + } + + self.assign(&mid_eval); + + Ok(()) + } + + /// Whether the evaluation is on a domain boundary and converged there + /// (`ConvergedOnBoundary`): pushing the boundary (right-on counts), and + /// perpendicular to the boundary curve within `tol` or pushing off its end at a + /// corner; which boundary ([`SurfParamType::Neither`] when not pushing one). + /// + /// # Errors + /// + /// [`SolidError::Unported`] from the evaluation. + pub fn converged_on_boundary(&mut self, tol: f64) -> Result<(bool, SurfParamType)> { + // Neg tol: if right on, do the check. + let (pushing, ret_val) = self.pushing_boundary(-EFF_ZERO)?; + if !pushing { + return Ok((false, SurfParamType::Neither)); + } + + self.eval()?; + let domain = *self.base.domain(); + + if matches!( + ret_val, + SurfParamType::UMin | SurfParamType::UMax | SurfParamType::Both + ) { + if self.error_norm(ErrorNormType::GAbs)? < tol { + return Ok((true, ret_val)); // perpendicular to boundary curve. + } + + let param = self.uv().y; + let crv_domain = domain.v_interval(); + if param < crv_domain.min() + tol { + return Ok((self.error_norm(ErrorNormType::G)? > -tol, ret_val)); + } else if param > crv_domain.max() - tol { + return Ok((self.error_norm(ErrorNormType::G)? < tol, ret_val)); + } + } + if matches!( + ret_val, + SurfParamType::VMin | SurfParamType::VMax | SurfParamType::Both + ) { + if self.error_norm(ErrorNormType::FAbs)? < tol { + return Ok((true, ret_val)); // perpendicular to boundary curve. + } + + let param = self.uv().x; + let crv_domain = domain.u_interval(); + if param < crv_domain.min() + tol { + return Ok((self.error_norm(ErrorNormType::F)? > -tol, ret_val)); + } else if param > crv_domain.max() - tol { + return Ok((self.error_norm(ErrorNormType::F)? < tol, ret_val)); + } + } + + Ok((false, ret_val)) + } + + /// Whether the UV is on a domain boundary and the plane step towards the test + /// point leaves the domain there by more than `tol` (`PushingBoundary`; a + /// negative `tol` makes a right-on point count), and which boundary: + /// [`SurfParamType::UMin`], [`SurfParamType::UMax`], [`SurfParamType::VMin`], + /// [`SurfParamType::VMax`], [`SurfParamType::Both`] (a corner) or + /// [`SurfParamType::Neither`]. + /// + /// # Errors + /// + /// [`SolidError::Unported`] from the evaluation. + pub fn pushing_boundary(&mut self, tol: f64) -> Result<(bool, SurfParamType)> { + let mut ret_val = SurfParamType::Neither; + + self.eval()?; + + let b = &self.base; + let test_pt = b.test_point(); + let dir_vec = b.dir_vec(); + let directed = b.solver_op() == SolverOperation::Signed3dDirectedMaximize; + let domain = *b.domain(); + + let mut gap = test_pt - self.pos()?; + if directed { + let dot = gap.dot(dir_vec); + gap = dot * dir_vec; + } + + // Method: find the uv direction that the 3d step will take us, + // and see whether it's heading out of bounds. + let su = self.su()?; + let sv = self.sv()?; + let Ok(uv_step) = vec_lin_comb_two_vectors(su, sv, gap) else { + return Ok((false, SurfParamType::Neither)); + }; + + let bdry_tol = EFF_ZERO; + let pushing_tol = tol; + let uv = self.uv(); + + if uv.x < domain.u_interval().min() + bdry_tol && uv_step.x < -pushing_tol { + ret_val = SurfParamType::UMin; + } else if uv.x > domain.u_interval().max() - bdry_tol && uv_step.x > pushing_tol { + ret_val = SurfParamType::UMax; + } + + if uv.y < domain.v_interval().min() + bdry_tol && uv_step.y < -pushing_tol { + ret_val = if ret_val == SurfParamType::Neither { + SurfParamType::VMin + } else { + SurfParamType::Both + }; + } else if uv.y > domain.v_interval().max() - bdry_tol && uv_step.y > pushing_tol { + ret_val = if ret_val == SurfParamType::Neither { + SurfParamType::VMax + } else { + SurfParamType::Both + }; + } + + Ok((ret_val != SurfParamType::Neither, ret_val)) + } + + /// When the UV is within `tol` of a domain boundary but not on it, snaps to the + /// boundary (or boundaries) if that is no worse (`SnapToBoundaries`). Returns + /// whether the stepper moved. + /// + /// # Errors + /// + /// [`SolidError::Unported`] from the evaluations. + pub fn snap_to_boundaries(&mut self, tol: f64) -> Result { + let mut uv = self.uv(); + let domain = *self.base.domain(); + + if !domain.is_point_on_boundary(uv, tol).0 { + return Ok(false); + } + + if domain.is_point_on_boundary(uv, 0.0).0 { + return Ok(false); + } + + // Ok, we're within dTol of a boundary but not exactly on it. + let op = self.base.solver_op(); + let norm_type = if op == SolverOperation::Normalize + || op == SolverOperation::Signed3dDirectedMaximize + { + ErrorNormType::SumAbs + } else { + ErrorNormType::Gap + }; + let current_error = self.error_norm(norm_type)?; + + let save_uv = uv; + + if uv.x <= domain.min().x + tol { + uv.x = domain.min().x; + } else if uv.x >= domain.max().x - tol { + uv.x = domain.max().x; + } + + if uv.y <= domain.min().y + tol { + uv.y = domain.min().y; + } else if uv.y >= domain.max().y - tol { + uv.y = domain.max().y; + } + + self.set_uv(uv); + + let bdry_error = self.error_norm(norm_type)?; + + let mut ret_val = true; + + // If not better, restore original uv. + if bdry_error > current_error { + // If exactly equal, leave on bdry. + self.set_uv(save_uv); + ret_val = false; + } + + Ok(ret_val) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::geom::global_solver::SolutionRequested; + use crate::geom::solution_array::ObjectKind; + use crate::nurbs::b_spline_surface::BSplineSurface; + + #[derive(Debug, Clone, Copy, PartialEq)] + struct Srf; + + impl SolutionObject for Srf { + fn kind(&self) -> ObjectKind { + ObjectKind::Surface + } + fn natural_uv_domain(&self) -> Option { + Some(unit_dom()) + } + } + + fn bilinear(twist: f64) -> Surface { + Surface::from( + BSplineSurface::create_bilinear_surface( + Point3::new(0.0, 0.0, 0.0), + Point3::new(1.0, 0.0, 0.0), + Point3::new(0.0, 1.0, 0.0), + Point3::new(1.0, 1.0, twist), + ) + .unwrap(), + ) + } + + fn unit_dom() -> Extent2d { + Extent2d::from_bounds(0.0, 0.0, 1.0, 1.0) + } + + fn global<'s>( + k: &Kernel, + s: &'s Surface, + op: SolverOperation, + test: Point3, + ) -> SurfPtGlobalSolver<'s> { + SurfPtGlobalSolver::new( + k, + s, + &unit_dom(), + op, + test, + 1e-6, + None, + SolutionRequested::Single, + ) + .unwrap() + } + + fn stepper<'s>( + k: &Kernel, + s: &'s Surface, + test: Point3, + uv: Point2, + op: SolverOperation, + ) -> SurfPtSolveStepper<'s> { + SurfPtSolveStepper::new(k, s, &unit_dom(), test, uv, 1e-12, 1e-6, op).unwrap() + } + + #[test] + fn rough_and_polygon_guesses() { + let k = Kernel::default(); + let s = bilinear(0.0); + let mut g = global( + &k, + &s, + SolverOperation::Minimize, + Point3::new(0.25, 0.75, 2.0), + ); + let r = g.find_rough_guess().unwrap(); + assert!((r.uv.x - 0.25).abs() < 1e-12 && (r.uv.y - 0.75).abs() < 1e-12); + let p = g.find_polygon_guess((0, 0)).unwrap(); + assert!(!p.suspicious); + assert_eq!(p.uv, r.uv); + + // a twisted bilinear surface with a test point far off its net: the grid + // search still lands near the closest point + let s = bilinear(0.8); + let mut g = global( + &k, + &s, + SolverOperation::Minimize, + Point3::new(0.9, 0.9, 0.7), + ); + let r = g.find_rough_guess().unwrap(); + assert!(r.uv.x > 0.5 && r.uv.y > 0.5); + } + + #[test] + fn refinements() { + let k = Kernel::default(); + let s = bilinear(0.0); + let mut g = global( + &k, + &s, + SolverOperation::Minimize, + Point3::new(0.3, 0.6, 1.0), + ); + // the seam test of a B-spline is not ported yet + let mut uv = Point2::new(0.5, 0.5); + assert!(g.refine_guess(&mut uv).unwrap_err().is_unported()); + + let mut best = BIG_DOUBLE; + let mut guess = Point2::new(0.0, 0.0); + g.refine_patch_by_grid(&unit_dom(), 3, 1, &mut guess, &mut best) + .unwrap(); + assert!((best - 1.0).abs() < 1e-12); + + let uv = g.refine_guess_by_grid(Point2::new(0.9, 0.1), 1, 1).unwrap(); + assert!(g.calc_func_val(uv).unwrap() <= g.calc_func_val(Point2::new(0.9, 0.1)).unwrap()); + } + + #[test] + fn seams_and_solve() { + let k = Kernel::default(); + let s = bilinear(0.0); + let mut g = global( + &k, + &s, + SolverOperation::Minimize, + Point3::new(0.0, 0.5, 0.0), + ); + let mut sol: Solution = Solution::new(); + sol.num_variables = 2; + sol.num_objects = 1; + sol.objects[0] = Some(Srf); + sol.start.parameters[0] = 0.0; + sol.start.parameters[1] = 0.5; + // an open surface: just the solution + let sols = g.check_seams(&k, 1e-12, &sol, Some(Srf)).unwrap(); + assert_eq!(sols.len(), 1); + + // the local solve is not ported yet + let err = g.solve(&k, Some(Srf)).unwrap_err(); + assert!(err.is_unported()); + } + + #[test] + fn stepper_norms_and_tests() { + let k = Kernel::default(); + let s = bilinear(0.0); + let mut st = stepper( + &k, + &s, + Point3::new(0.5, 0.25, 2.0), + Point2::new(0.25, 0.5), + SolverOperation::Minimize, + ); + assert_eq!(st.error_norm(ErrorNormType::F).unwrap(), -0.25); + assert_eq!(st.error_norm(ErrorNormType::G).unwrap(), 0.25); + assert_eq!(st.error_norm(ErrorNormType::SumAbs).unwrap(), 0.5); + assert_eq!(st.error_norm(ErrorNormType::Sum).unwrap(), 0.0); + assert_eq!(st.error_norm(ErrorNormType::Max).unwrap(), -0.25); + assert_eq!(st.error_norm(ErrorNormType::MaxAbs).unwrap(), 0.25); + assert_eq!(st.is_singular().unwrap(), (false, SurfParamType::Neither)); + assert!(!st.converged(1e-6).unwrap()); + + // a Newton step on a plane lands on the foot point + let mut next = SurfPtSolveStepper::copy_of(&st); + assert_eq!( + next.step(&mut st, false).unwrap(), + TerminationReason::StillIterating + ); + assert!((next.uv().x - 0.5).abs() < 1e-12 && (next.uv().y - 0.25).abs() < 1e-12); + assert!(next.converged(1e-6).unwrap()); + assert!(next.converged_with(false).unwrap()); + + // pushing a boundary + let mut b = stepper( + &k, + &s, + Point3::new(2.0, 0.5, 0.0), + Point2::new(1.0, 0.5), + SolverOperation::Minimize, + ); + assert_eq!( + b.pushing_boundary(1e-9).unwrap(), + (true, SurfParamType::UMax) + ); + assert_eq!( + b.converged_on_boundary(1e-9).unwrap(), + (true, SurfParamType::UMax) + ); + + // snapping + let mut sn = stepper( + &k, + &s, + Point3::new(1.5, 0.5, 0.0), + Point2::new(1.0 - 1e-4, 0.5), + SolverOperation::Minimize, + ); + assert!(sn.snap_to_boundaries(1e-3).unwrap()); + assert_eq!(sn.uv(), Point2::new(1.0, 0.5)); + assert!(!sn.snap_to_boundaries(1e-3).unwrap()); + } + + #[test] + fn adjust_and_line_search() { + let k = Kernel::default(); + let s = bilinear(0.3); + let test = Point3::new(0.5, 0.5, 1.0); + let prev = stepper( + &k, + &s, + test, + Point2::new(0.5, 0.5), + SolverOperation::Minimize, + ); + let mut this = SurfPtSolveStepper::copy_of(&prev); + this.set_uv(Point2::new(1.5, 0.5)); + let (adjusted, scale) = this.adjust_step(&prev).unwrap(); + assert!(adjusted); + assert!((scale - 0.5).abs() < 1e-12); + assert_eq!(this.uv(), Point2::new(1.0, 0.5)); + assert_eq!(this.bound_hits_u(), 1); + + // no step, outside: clamped, SM_ERR + let mut out = SurfPtSolveStepper::copy_of(&prev); + out.set_uv(Point2::new(1.5, 0.5)); + let prev_out = SurfPtSolveStepper::copy_of(&out); + assert!(out.adjust_step(&prev_out).is_err()); + assert_eq!(out.uv(), Point2::new(1.0, 0.5)); + + let mut prev = stepper( + &k, + &s, + test, + Point2::new(0.1, 0.1), + SolverOperation::Minimize, + ); + let mut this = SurfPtSolveStepper::copy_of(&prev); + this.set_uv(Point2::new(0.9, 0.9)); + let before = this.error_norm(ErrorNormType::Gap).unwrap(); + this.line_search(&mut prev).unwrap(); + assert!(this.error_norm(ErrorNormType::Gap).unwrap() <= before); + + assert!(this + .line_search_bracketted(&mut prev, ErrorNormType::Gap) + .is_err()); + } +} diff --git a/harness/cpp/shim/surface_surface_point_solve_p2.cpp b/harness/cpp/shim/surface_surface_point_solve_p2.cpp new file mode 100644 index 0000000000..99be367cc9 --- /dev/null +++ b/harness/cpp/shim/surface_surface_point_solve_p2.cpp @@ -0,0 +1,516 @@ +// SPDX-License-Identifier: Apache-2.0 +// difftest shim for unit `surface.surface_point_solve_p2`: the global solver's guesses, +// refinements, seam solutions and Solve, and the local stepper's step adjustment, norms, +// convergence tests, steps and line searches, of the file-local classes of upstream +// source/SMLib/src/SmSurfacePointSolve.cpp. +// +// As the p1 shim (surface_surface_point_solve_p1.cpp), this translation unit compiles the +// upstream source itself with the classes renamed (sc_spsp2_*), SmSurface::GlobalPointSolve_1 +// / LocalPointSolve_1 renamed before any header is read, and `private` / `protected` made +// public for the class definitions only. +// +// One entry point, sc_surface_surface_point_solve_p2_eval(op, in, n_in, out, cap, out_len). +// Layouts are documented in crates/solidcore-ffi/src/units/surface_surface_point_solve_p2.rs, +// which validates every input. +#define GlobalPointSolve_1 sc_spsp2_GlobalPointSolve_1 +#define LocalPointSolve_1 sc_spsp2_LocalPointSolve_1 + +#include "sc_shim.h" + +#include "StdAfx.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#define SmSurfacePointSolver sc_spsp2_SmSurfacePointSolver +#define SmSurfPtGlobalSolver sc_spsp2_SmSurfPtGlobalSolver +#define SmSurfPtSolveStepper sc_spsp2_SmSurfPtSolveStepper +#define SmSrfPtErrorNormType sc_spsp2_SmSrfPtErrorNormType +#define GetTrianglesDist sc_spsp2_GetTrianglesDist +#define WellBehavedControlNet sc_spsp2_WellBehavedControlNet +#define private public +#define protected public +#include "../../../upstream/source/SMLib/src/SmSurfacePointSolve.cpp" +#undef private +#undef protected + +namespace sc_spsp2 +{ +enum +{ + OP_ROUGH = 0, + OP_POLY = 1, + OP_REFINE = 2, + OP_PATCH = 3, + OP_GRID = 4, + OP_SEAMS = 5, + OP_SOLVE = 6, + OP_ADJUST = 7, + OP_NORMS = 8, + OP_CONV = 9, + OP_PUSH = 10, + OP_SNAP = 11, + OP_STEP = 12, + OP_LINE = 13, + OP_LINE_B = 14, +}; + +struct Rd +{ + const double* in; + int64_t at; + double num() { return in[at++]; } + int64_t idx() { return (int64_t)in[at++]; } + bool flag() { return in[at++] != 0.0; } + SmVector3d vec() + { + double x = num(); + double y = num(); + double z = num(); + return SmVector3d(x, y, z); + } + SmPoint2d uv() + { + double x = num(); + double y = num(); + return SmPoint2d(x, y); + } + SmExtent2d domain() + { + double a = num(); + double b = num(); + double c = num(); + double d = num(); + return SmExtent2d(a, b, c, d); + } +}; + +double st(SmStatus s) +{ + return s == SM_SUCCESS ? 0.0 : (double)s; +} + +double flg(bool b) +{ + return b ? 1.0 : 0.0; +} + +void put_uv(std::vector& o, const SmVector2d& v) +{ + o.push_back(v.x); + o.push_back(v.y); +} + +gw_SURFACE* read_nurb(Rd& rd) +{ + gw_INDEX n = (gw_INDEX)rd.idx(); + gw_INDEX m = (gw_INDEX)rd.idx(); + gw_DEGREE p = (gw_DEGREE)rd.idx(); + gw_DEGREE q = (gw_DEGREE)rd.idx(); + gw_INDEX r = (gw_INDEX)rd.idx(); + gw_INDEX s = (gw_INDEX)rd.idx(); + gw_SURFACE* sur = sm_AllocateNurbSurface(n, m, p, q, r, s); + for (gw_INDEX i = 0; i <= n; i++) + for (gw_INDEX j = 0; j <= m; j++) + { + sur->net->Pw[i][j].x = rd.num(); + sur->net->Pw[i][j].y = rd.num(); + sur->net->Pw[i][j].z = rd.num(); + sur->net->Pw[i][j].w = rd.num(); + } + for (gw_INDEX i = 0; i <= r; i++) + sur->knu->U[i] = rd.num(); + for (gw_INDEX i = 0; i <= s; i++) + sur->knv->U[i] = rd.num(); + return sur; +} + +// `oob form nurb`. +SmBSplineSurface* read_bspline(Rd& rd, const SmContext& ctx) +{ + bool oob = rd.flag(); + SmBSplineSurfaceForm form = (SmBSplineSurfaceForm)rd.idx(); + SmBSplineSurface* s = new (ctx) SmBSplineSurface(read_nurb(rd), FALSE, NULL); + s->GetOutOfBoundsEnabled() = oob ? TRUE : FALSE; + s->SetBSplineSurfaceForm(form); + return s; +} + +// A surface `kind ...`: 0 B-spline `oob form nurb`; 1 bounded plane `o x y sx sy D`; +// 2 curve-bounded surface (implicit outer, no boundary curves) over a B-spline block. +SmSurface* read_surface(Rd& rd, const SmContext& ctx) +{ + int64_t kind = rd.idx(); + if (kind == 1) + { + SmPoint3d o = rd.vec(); + SmVector3d x = rd.vec(); + SmVector3d y = rd.vec(); + double sx = rd.num(); + double sy = rd.num(); + SmExtent2d d = rd.domain(); + return new (ctx) SmPlane(o, x, y, SmVector2d(sx, sy), d, &ctx); + } + SmBSplineSurface* b = read_bspline(rd, ctx); + if (kind == 0) + return b; + // the curve-bounded surface owns (deletes) its basis + SmTArray c3d(ctx), cuv(ctx); + return new (ctx) SmCurveBoundedSurface(*b, &c3d, &cuv, TRUE); +} + +// `W S`: the surface, wrapped in an SmSTEPSurface owning it when W = 1. +SmSurface* read_wrapped(Rd& rd, const SmContext& ctx) +{ + bool wrap = rd.flag(); + SmSurface* s = read_surface(rd, ctx); + if (!wrap) + return s; + return new (ctx) SmSTEPSurface(*s, TRUE); +} + +// The problem: `W S D op test tol`. +struct Problem +{ + SmSurface* srf; + SmExtent2d dom; + SmSolverOperationType op; + SmPoint3d test; + double tol; +}; + +Problem read_problem(Rd& rd, const SmContext& ctx) +{ + Problem p; + p.srf = read_wrapped(rd, ctx); + p.dom = rd.domain(); + p.op = (SmSolverOperationType)rd.idx(); + p.test = rd.vec(); + p.tol = rd.num(); + return p; +} + +// `(u v value nobj nvar obj0_is_srf)` of each solution, after their count. +void put_solutions(std::vector& o, SmSolutionArray& sols, const SmSurface* srf) +{ + o.push_back((double)sols.GetSize()); + for (ULONG i = 0; i < sols.GetSize(); i++) + { + const SmSolution& s = sols[i]; + o.push_back(s.m_vStart.m_adParameters[0]); + o.push_back(s.m_vStart.m_adParameters[1]); + o.push_back(s.m_vStart.m_dSolutionValue); + o.push_back((double)s.m_lNumObjects); + o.push_back((double)s.m_lNumVariables); + o.push_back(flg(s.m_apObjects[0] == (const SmObject*)srf)); + } +} + +// `P`: FindRoughGuess. +void op_rough(Rd& rd, const SmContext& ctx, std::vector& o) +{ + Problem p = read_problem(rd, ctx); + SmSurfPtGlobalSolver g(p.srf, p.dom, p.op, p.test, p.tol, NULL, SM_SR_SINGLE); + SmPoint2d uv(0.0, 0.0); + ULONG iu = 7, iv = 7; + SmStatus s = g.FindRoughGuess(uv, iu, iv); + o.push_back(st(s)); + if (s != SM_SUCCESS) + return; + put_uv(o, uv); + o.push_back((double)iu); + o.push_back((double)iv); +} + +// `P iu iv`: FindPolygonGuess from the given best indices. +void op_poly(Rd& rd, const SmContext& ctx, std::vector& o) +{ + Problem p = read_problem(rd, ctx); + ULONG iu = (ULONG)rd.idx(); + ULONG iv = (ULONG)rd.idx(); + SmSurfPtGlobalSolver g(p.srf, p.dom, p.op, p.test, p.tol, NULL, SM_SR_SINGLE); + SmPoint2d uv(0.0, 0.0); + SmBoolean susp = FALSE; + SmStatus s = g.FindPolygonGuess(uv, iu, iv, susp); + o.push_back(st(s)); + if (s != SM_SUCCESS) + return; + put_uv(o, uv); + o.push_back((double)iu); + o.push_back((double)iv); + o.push_back(flg(susp != FALSE)); +} + +// `P uv`: RefineGuess. +void op_refine(Rd& rd, const SmContext& ctx, std::vector& o) +{ + Problem p = read_problem(rd, ctx); + SmPoint2d uv = rd.uv(); + SmSurfPtGlobalSolver g(p.srf, p.dom, p.op, p.test, p.tol, NULL, SM_SR_SINGLE); + ULONG iu = 0, iv = 0; + SmStatus s = g.RefineGuess(uv, iu, iv); + o.push_back(st(s)); + put_uv(o, uv); +} + +// `P sub nu nv uv best`: RefinePatchByGrid. +void op_patch(Rd& rd, const SmContext& ctx, std::vector& o) +{ + Problem p = read_problem(rd, ctx); + SmExtent2d sub = rd.domain(); + int nu = (int)rd.idx(); + int nv = (int)rd.idx(); + SmPoint2d uv = rd.uv(); + double best = rd.num(); + SmSurfPtGlobalSolver g(p.srf, p.dom, p.op, p.test, p.tol, NULL, SM_SR_SINGLE); + SmStatus s = g.RefinePatchByGrid(sub, nu, nv, uv, best); + o.push_back(st(s)); + put_uv(o, uv); + o.push_back(best); +} + +// `P uv iu iv`: RefineGuessByGrid. +void op_grid(Rd& rd, const SmContext& ctx, std::vector& o) +{ + Problem p = read_problem(rd, ctx); + SmPoint2d uv = rd.uv(); + ULONG iu = (ULONG)rd.idx(); + ULONG iv = (ULONG)rd.idx(); + SmSurfPtGlobalSolver g(p.srf, p.dom, p.op, p.test, p.tol, NULL, SM_SR_SINGLE); + SmStatus s = g.RefineGuessByGrid(uv, iu, iv); + o.push_back(st(s)); + put_uv(o, uv); +} + +// `P uv value`: CheckSeams of the solution (u, v) on the surface with that value. +void op_seams(Rd& rd, const SmContext& ctx, std::vector& o) +{ + Problem p = read_problem(rd, ctx); + SmPoint2d uv = rd.uv(); + double value = rd.num(); + SmSurfPtGlobalSolver g(p.srf, p.dom, p.op, p.test, p.tol, NULL, SM_SR_SINGLE); + SmSolution sol; + sol.m_lNumVariables = 2; + sol.m_lNumObjects = 1; + sol.m_apObjects[0] = p.srf; + sol.m_vStart[0] = uv.x; + sol.m_vStart[1] = uv.y; + sol.m_vStart.m_dSolutionValue = value; + SmSolutionArray sols; + SmStatus s = g.CheckSeams(g.m_dTightTol, sol, sols); + o.push_back(st(s)); + put_solutions(o, sols, p.srf); +} + +// `P`: Solve. +void op_solve(Rd& rd, const SmContext& ctx, std::vector& o) +{ + Problem p = read_problem(rd, ctx); + SmSurfPtGlobalSolver g(p.srf, p.dom, p.op, p.test, p.tol, NULL, SM_SR_SINGLE); + SmSolutionArray sols; + SmStatus s = g.Solve(sols); + o.push_back(st(s)); + if (s != SM_SUCCESS) + return; + put_solutions(o, sols, p.srf); +} + +// `HitBound(which, sign)` |n| times. +void hits(SmSurfPtSolveStepper& s, SmSurfParamType which, int64_t n) +{ + int sign = n < 0 ? -1 : 1; + for (int64_t i = 0; i < (n < 0 ? -n : n); i++) + s.HitBound(which, sign); +} + +// The stepper of `P uv0 tight`. +SmSurfPtSolveStepper read_stepper(Rd& rd, const Problem& p) +{ + SmPoint2d uv0 = rd.uv(); + double tight = rd.num(); + return SmSurfPtSolveStepper(p.srf, p.dom, p.test, uv0, tight, p.tol, p.op); +} + +// `UV BoundHitsU BoundHitsV NRStepSuspect IterCount` of a stepper. +void put_state(std::vector& o, SmSurfPtSolveStepper& s) +{ + put_uv(o, s.UV()); + o.push_back((double)s.BoundHitsU()); + o.push_back((double)s.BoundHitsV()); + o.push_back(flg(s.NRStepSuspect() != FALSE)); + o.push_back((double)s.IterCount()); +} + +// `P uv0 tight uv1 hu hv`: AdjustStep of a copy of the stepper at uv0 (after hu / hv +// boundary hits) moved to uv1. +void op_adjust(Rd& rd, const SmContext& ctx, std::vector& o) +{ + Problem p = read_problem(rd, ctx); + SmSurfPtSolveStepper prev = read_stepper(rd, p); + SmPoint2d uv1 = rd.uv(); + int64_t hu = rd.idx(); + int64_t hv = rd.idx(); + hits(prev, SM_SP_U, hu); + hits(prev, SM_SP_V, hv); + SmSurfPtSolveStepper next(prev); + next.SetUV(uv1); + SmBoolean adjusted = FALSE; + double scale = -7.0; + SmStatus s = next.AdjustStep(prev, adjusted, scale); + o.push_back(st(s)); + o.push_back(flg(adjusted != FALSE)); + o.push_back(scale); + put_state(o, next); +} + +// `P uv0 tight`: the nine ErrorNorms, IsSingular. +void op_norms(Rd& rd, const SmContext& ctx, std::vector& o) +{ + Problem p = read_problem(rd, ctx); + SmSurfPtSolveStepper s = read_stepper(rd, p); + o.push_back(0.0); + for (int n = 0; n <= (int)SM_SPN_GAP; n++) + o.push_back(s.ErrorNorm((SmSrfPtErrorNormType)n)); + SmSurfParamType which = SM_SP_UNKNOWN; + SmBoolean b = s.IsSingular(&which); + o.push_back(flg(b != FALSE)); + o.push_back((double)which); +} + +// `P uv0 tight ctol`: Converged(ctol), Converged(TRUE), Converged(FALSE). +void op_conv(Rd& rd, const SmContext& ctx, std::vector& o) +{ + Problem p = read_problem(rd, ctx); + SmSurfPtSolveStepper s = read_stepper(rd, p); + double tol = rd.num(); + o.push_back(0.0); + o.push_back(flg(s.Converged(tol) != FALSE)); + o.push_back(flg(s.Converged((SmBoolean)TRUE) != FALSE)); + o.push_back(flg(s.Converged((SmBoolean)FALSE) != FALSE)); +} + +// `P uv0 tight ptol`: PushingBoundary(ptol), ConvergedOnBoundary(ptol) (fresh steppers). +void op_push(Rd& rd, const SmContext& ctx, std::vector& o) +{ + Problem p = read_problem(rd, ctx); + SmSurfPtSolveStepper s = read_stepper(rd, p); + double tol = rd.num(); + SmSurfPtSolveStepper c(s); + o.push_back(0.0); + SmSurfParamType which = SM_SP_NEITHER; + SmBoolean b = s.PushingBoundary(tol, &which); + o.push_back(flg(b != FALSE)); + o.push_back((double)which); + which = SM_SP_NEITHER; + b = c.ConvergedOnBoundary(tol, &which); + o.push_back(flg(b != FALSE)); + o.push_back((double)which); +} + +// `P uv0 tight stol`: SnapToBoundaries(stol). +void op_snap(Rd& rd, const SmContext& ctx, std::vector& o) +{ + Problem p = read_problem(rd, ctx); + SmSurfPtSolveStepper s = read_stepper(rd, p); + double tol = rd.num(); + SmBoolean b = s.SnapToBoundaries(tol); + o.push_back(0.0); + o.push_back(flg(b != FALSE)); + put_state(o, s); +} + +// `P uv0 tight plane`: Step of a copy of the stepper at uv0, from it. +void op_step(Rd& rd, const SmContext& ctx, std::vector& o) +{ + Problem p = read_problem(rd, ctx); + SmSurfPtSolveStepper prev = read_stepper(rd, p); + bool plane = rd.flag(); + SmSurfPtSolveStepper next(prev); + SmTerminationReasonType r = next.Step(prev, plane ? TRUE : FALSE); + o.push_back(0.0); + o.push_back((double)r); + put_state(o, next); +} + +// `P uv0 tight uv1`: LineSearch of a copy of the stepper at uv0 moved to uv1, from it. +void op_line(Rd& rd, const SmContext& ctx, std::vector& o) +{ + Problem p = read_problem(rd, ctx); + SmSurfPtSolveStepper prev = read_stepper(rd, p); + SmPoint2d uv1 = rd.uv(); + SmSurfPtSolveStepper next(prev); + next.SetUV(uv1); + SmStatus s = next.LineSearch(prev); + o.push_back(st(s)); + put_state(o, next); +} + +// `P uv0 tight uv1 norm`: LineSearchBracketted as op_line. +void op_line_b(Rd& rd, const SmContext& ctx, std::vector& o) +{ + Problem p = read_problem(rd, ctx); + SmSurfPtSolveStepper prev = read_stepper(rd, p); + SmPoint2d uv1 = rd.uv(); + SmSrfPtErrorNormType norm = (SmSrfPtErrorNormType)rd.idx(); + SmSurfPtSolveStepper next(prev); + next.SetUV(uv1); + SmStatus s = next.LineSearchBracketted(prev, norm); + o.push_back(st(s)); + put_state(o, next); +} +} // namespace sc_spsp2 + +SC_EXPORT int32_t sc_surface_surface_point_solve_p2_eval(int32_t op, const double* in, int64_t n_in, double* out, + int64_t cap, int64_t* out_len) +{ + (void)n_in; + SC_TRY + SmContext ctx; + std::vector o; + sc_spsp2::Rd rd{in, 0}; + switch (op) + { + case sc_spsp2::OP_ROUGH: sc_spsp2::op_rough(rd, ctx, o); break; + case sc_spsp2::OP_POLY: sc_spsp2::op_poly(rd, ctx, o); break; + case sc_spsp2::OP_REFINE: sc_spsp2::op_refine(rd, ctx, o); break; + case sc_spsp2::OP_PATCH: sc_spsp2::op_patch(rd, ctx, o); break; + case sc_spsp2::OP_GRID: sc_spsp2::op_grid(rd, ctx, o); break; + case sc_spsp2::OP_SEAMS: sc_spsp2::op_seams(rd, ctx, o); break; + case sc_spsp2::OP_SOLVE: sc_spsp2::op_solve(rd, ctx, o); break; + case sc_spsp2::OP_ADJUST: sc_spsp2::op_adjust(rd, ctx, o); break; + case sc_spsp2::OP_NORMS: sc_spsp2::op_norms(rd, ctx, o); break; + case sc_spsp2::OP_CONV: sc_spsp2::op_conv(rd, ctx, o); break; + case sc_spsp2::OP_PUSH: sc_spsp2::op_push(rd, ctx, o); break; + case sc_spsp2::OP_SNAP: sc_spsp2::op_snap(rd, ctx, o); break; + case sc_spsp2::OP_STEP: sc_spsp2::op_step(rd, ctx, o); break; + case sc_spsp2::OP_LINE: sc_spsp2::op_line(rd, ctx, o); break; + case sc_spsp2::OP_LINE_B: sc_spsp2::op_line_b(rd, ctx, o); break; + default: return SC_ERR_ARG; + } + return sc_copy_out(o.data(), (int64_t)o.size(), out, cap, out_len); + SC_CATCH +}