mirror of
https://github.com/bherbruck/solidcraft.git
synced 2026-11-10 18:55:14 +00:00
Compare commits
1 Commits
cancel-check
...
kernel
| Author | SHA1 | Date | |
|---|---|---|---|
| 6d01515911 |
+7
-11
@@ -7,8 +7,8 @@
|
||||
//! kind, by what it became; a number within [`ROUNDING`] of its size (relative, no floor) is
|
||||
//! counted as rounding, with its count and the largest.
|
||||
//!
|
||||
//! On the whole corpus (1,332 bodies, at #3187): 1,256 come back bit for bit, 36 within
|
||||
//! rounding (the largest 4.9e-16 relative), and 40 differ, each in exactly the classes
|
||||
//! On the whole corpus (1,333 bodies, at #3199 with kernel case 11991 (#3231) in): 1,259 come
|
||||
//! back bit for bit, 36 within rounding (the largest 4.9e-16 relative), and 38 differ, each in exactly the classes
|
||||
//! [`pinned`] lists for it.
|
||||
#![allow(
|
||||
clippy::unwrap_used,
|
||||
@@ -390,8 +390,8 @@ fn pinned() -> BTreeMap<String, BTreeSet<String>> {
|
||||
],
|
||||
),
|
||||
// a label its form is not proved on goes as the form alone (#2540): cylinders whose
|
||||
// frames are about 1e-14 rad off their forms (two of them at the rounding of a derived
|
||||
// bound, #3185's thread)
|
||||
// frames are off their forms past the rounding `evaluates_as` derives (margins 1.23 to
|
||||
// 2.03, tests/issue_3185_label_margins.rs)
|
||||
(
|
||||
label("cylinder -> rational_spline", &["radius"]),
|
||||
&[
|
||||
@@ -401,18 +401,14 @@ fn pinned() -> BTreeMap<String, BTreeSet<String>> {
|
||||
"fillets/2903-original-rim-round.json:b0",
|
||||
],
|
||||
),
|
||||
// tori whose frame and fitted sweep are about 1e-14 rad off their forms, rightly
|
||||
// refused (#3185)
|
||||
// a torus label whose frame is off its form's own first sweep row past the rounding
|
||||
// `evaluates_as` derives (margin 1.69, tests/issue_3185_label_margins.rs)
|
||||
(
|
||||
label(
|
||||
"torus -> rational_spline",
|
||||
&["major_radius", "minor_radius"],
|
||||
),
|
||||
&[
|
||||
"sweeps/3897-sweep-ring-along-a-short-arc-joined-to-a-collar.json:b6",
|
||||
"sweeps/4280-circular-pattern-of-sweep-cut-groove-on-hub.json:b46",
|
||||
"sweeps/4280-circular-pattern-of-sweep-cut-groove-on-hub.json:b85",
|
||||
],
|
||||
&["sweeps/4280-circular-pattern-of-sweep-cut-groove-on-hub.json:b46"],
|
||||
),
|
||||
// an extruded circle the flat hop left a spline, proved a cylinder
|
||||
(
|
||||
|
||||
@@ -0,0 +1,258 @@
|
||||
//! (#3185) An analytic label is kept on its NURBS form when the kernel's analytic surface,
|
||||
//! carrying the label's numbers, evaluates as the form within the rounding
|
||||
//! `AnalyticSupport::build` derives term by term (`surface::implicit`'s `evaluates_as`):
|
||||
//! [`AnalyticSupport::parametrisation_margin`] is the largest ratio of distance to bound over
|
||||
//! the form's Greville points, at most 1 where the label stands.
|
||||
//!
|
||||
//! - The corpus labels #3137's parity run refused under the old constant `gamma(16 (p + q +
|
||||
//! 2)) size`, each pinned to its margin as measured under the derived bound: the three tori
|
||||
//! of b6 and b85 (0.02 to 0.50) stand; b46's torus (1.69) and the four cylinders (b94 1.23,
|
||||
//! b43 1.40, the 2903 rims 1.52 and 2.03) do not: their frames are off their forms past
|
||||
//! the rounding at their size. Since #3199 the kernel's recognised sweep frame is the form's
|
||||
//! own (its end meridians), so a label is measured against the form, not against a fit
|
||||
//! that shared the label's error.
|
||||
//! - A torus the kernel makes at the groove's sizes stands on its own form, and the same label
|
||||
//! turned about its axis by four times what the bound allows is refused.
|
||||
#![allow(
|
||||
clippy::unwrap_used,
|
||||
clippy::expect_used,
|
||||
clippy::panic,
|
||||
clippy::indexing_slicing,
|
||||
clippy::print_stderr,
|
||||
clippy::disallowed_methods
|
||||
)] // test (the turn's cos and sin)
|
||||
|
||||
use std::fmt::Write as _;
|
||||
use std::io::Read as _;
|
||||
use std::path::{Path, PathBuf};
|
||||
|
||||
use serde_json::Value;
|
||||
use solidcore::blend::fillet_section::station::as_b_spline;
|
||||
use solidcore::geom::axis2_placement::Axis2Placement;
|
||||
use solidcore::geom::extent2d::Extent2d;
|
||||
use solidcore::geom::vector3d::{Point3, Vec3};
|
||||
use solidcore::kernel::{Kernel, KernelOptions, UpstreamFixes};
|
||||
use solidcore::surface::from_parameters::from_parameters;
|
||||
use solidcore::surface::implicit::AnalyticSupport;
|
||||
use solidcore::surface::surface::{Surface, SurfaceKind};
|
||||
use solidcore::surface::torus::{Torus, TorusOptions};
|
||||
use solidcore_record::brep::Brep as Record;
|
||||
|
||||
fn corpus() -> PathBuf {
|
||||
Path::new(env!("CARGO_MANIFEST_DIR")).join("../solidcore/tests/kernel_cases")
|
||||
}
|
||||
|
||||
/// The corpus's kernel, its fixes on (as the parity run and SolidCraft run it).
|
||||
fn kernel() -> Kernel {
|
||||
let v: Value =
|
||||
serde_json::from_str(&std::fs::read_to_string(corpus().join("fixes.json")).unwrap())
|
||||
.unwrap();
|
||||
let mut fixes = UpstreamFixes::default();
|
||||
for flag in v.as_array().unwrap() {
|
||||
*fixes.flag_mut(flag.as_str().unwrap()).unwrap() = true;
|
||||
}
|
||||
let mut k = Kernel::new(&KernelOptions::default());
|
||||
k.fixes = fixes;
|
||||
k
|
||||
}
|
||||
|
||||
fn body(case: &str, id: &str) -> Record {
|
||||
let gz = corpus().join(format!("{case}.bodies.json.gz"));
|
||||
let mut text = String::new();
|
||||
flate2::read::GzDecoder::new(std::fs::File::open(&gz).unwrap())
|
||||
.read_to_string(&mut text)
|
||||
.unwrap();
|
||||
let records: Value = serde_json::from_str(&text).unwrap();
|
||||
serde_json::from_value(records[id].clone()).unwrap()
|
||||
}
|
||||
|
||||
fn placement(f: &solidcore::eval::Frame) -> Axis2Placement {
|
||||
let v = |a: [f64; 3]| Vec3::new(a[0], a[1], a[2]);
|
||||
Axis2Placement::from_stored(
|
||||
Point3::new(f.origin[0], f.origin[1], f.origin[2]),
|
||||
v(f.x_axis),
|
||||
v(f.y_axis),
|
||||
)
|
||||
}
|
||||
|
||||
/// A corpus face's label as the kernel's support, and its NURBS form.
|
||||
fn label_and_form(
|
||||
k: &Kernel,
|
||||
case: &str,
|
||||
id: &str,
|
||||
face: usize,
|
||||
) -> (
|
||||
AnalyticSupport,
|
||||
solidcore::nurbs::b_spline_surface::BSplineSurface,
|
||||
) {
|
||||
let r = body(case, id);
|
||||
let f = &r.faces[face];
|
||||
let nurbs = solidcore_convert::surface(f.nurbs.as_ref().unwrap()).unwrap();
|
||||
let (s, _) = from_parameters(k, &nurbs, f.domain).unwrap();
|
||||
let SurfaceKind::BSpline(form) = s.kind() else {
|
||||
panic!("{case} {id} {face}: the form is not a B-spline")
|
||||
};
|
||||
let support = match solidcore_convert::surface(&f.surface).unwrap() {
|
||||
solidcore::eval::Surface::Torus {
|
||||
frame,
|
||||
major_radius,
|
||||
minor_radius,
|
||||
..
|
||||
} => AnalyticSupport::Torus {
|
||||
placement: placement(&frame),
|
||||
major_radius,
|
||||
minor_radius,
|
||||
},
|
||||
solidcore::eval::Surface::Cylinder { frame, radius } => AnalyticSupport::Cylinder {
|
||||
placement: placement(&frame),
|
||||
radius,
|
||||
},
|
||||
other => panic!("{case} {id} {face}: {other:?}"),
|
||||
};
|
||||
(support, form.copy())
|
||||
}
|
||||
|
||||
/// (case, body, face, the margin measured at #3185, whether the label stands)
|
||||
const CASES: &[(&str, &str, usize, f64, bool)] = &[
|
||||
(
|
||||
"sweeps/3897-sweep-ring-along-a-short-arc-joined-to-a-collar",
|
||||
"b6",
|
||||
1,
|
||||
0.0201,
|
||||
true,
|
||||
),
|
||||
(
|
||||
"sweeps/4280-circular-pattern-of-sweep-cut-groove-on-hub",
|
||||
"b46",
|
||||
1,
|
||||
1.6943,
|
||||
false,
|
||||
),
|
||||
(
|
||||
"sweeps/4280-circular-pattern-of-sweep-cut-groove-on-hub",
|
||||
"b85",
|
||||
15,
|
||||
0.4995,
|
||||
true,
|
||||
),
|
||||
(
|
||||
"sweeps/4280-circular-pattern-of-sweep-cut-groove-on-hub",
|
||||
"b85",
|
||||
21,
|
||||
0.4245,
|
||||
true,
|
||||
),
|
||||
(
|
||||
"booleans/3456-deep-cut-through-drilled-taper-leaves-base",
|
||||
"b94",
|
||||
3,
|
||||
1.2339,
|
||||
false,
|
||||
),
|
||||
(
|
||||
"crashes/2243-step-export-flanged-enclosure",
|
||||
"b43",
|
||||
1,
|
||||
1.4008,
|
||||
false,
|
||||
),
|
||||
("fillets/2903-original-rim-chamfer", "b0", 1, 1.5243, false),
|
||||
("fillets/2903-original-rim-chamfer", "b0", 3, 2.0275, false),
|
||||
("fillets/2903-original-rim-round", "b0", 1, 1.5243, false),
|
||||
("fillets/2903-original-rim-round", "b0", 3, 2.0275, false),
|
||||
];
|
||||
|
||||
/// Each corpus label's margin, pinned to 1 % of what #3185 measured, on its side of 1, and
|
||||
/// `build` agreeing with it.
|
||||
#[test]
|
||||
fn the_corpus_labels_stand_or_fall_by_the_derived_bound() {
|
||||
let k = kernel();
|
||||
let mut report = String::new();
|
||||
let mut wrong = Vec::new();
|
||||
for &(case, id, face, pinned, stands) in CASES {
|
||||
let (support, form) = label_and_form(&k, case, id, face);
|
||||
let m = support.parametrisation_margin(&k, &form).unwrap().unwrap();
|
||||
let _ = writeln!(report, "{case} {id} face {face}: margin {m:.4}");
|
||||
let built = support.build(&k, form).is_ok();
|
||||
if (m <= 1.0) != stands || built != stands || (m - pinned).abs() > 0.01 * pinned {
|
||||
wrong.push(format!(
|
||||
"{case} {id} face {face}: margin {m} (pinned {pinned}, stands {stands}), build {built}"
|
||||
));
|
||||
}
|
||||
}
|
||||
eprintln!("{report}");
|
||||
assert!(wrong.is_empty(), "{wrong:#?}");
|
||||
}
|
||||
|
||||
/// A torus at the groove's sizes (#3185's b46: centre (19.375, -33.558, 37.5), R 33.657, r 2,
|
||||
/// a 43.6 degree sweep of the whole tube), made by the kernel: its own label stands on its
|
||||
/// form; turned about its axis by four times the turn the bound allows (the margin is linear
|
||||
/// in a small turn: a turn of `t` reads `m(t)`, so `4 t / m(t)` reads 4), it is refused.
|
||||
#[test]
|
||||
fn a_torus_label_turned_past_its_bound_is_refused() {
|
||||
let k = kernel();
|
||||
let o = Point3::new(19.374_999_999_999_975, -33.558_484_396_647_01, 37.5);
|
||||
let x = Vec3::new(
|
||||
-0.427_099_277_807_218_7,
|
||||
0.739_757_649_038_078_7,
|
||||
-0.519_946_946_895_745_4,
|
||||
);
|
||||
let y = Vec3::new(
|
||||
0.259_973_473_447_872_36,
|
||||
-0.450_287_264_631_874_2,
|
||||
-0.854_198_555_614_438_5,
|
||||
);
|
||||
let domain = Extent2d::from_bounds(0.0, -180.0, 43.602_818_972_703_645, 180.0);
|
||||
let made = Torus::new(
|
||||
&k,
|
||||
o,
|
||||
x,
|
||||
y,
|
||||
domain,
|
||||
33.657_280_044_590_65,
|
||||
2.0,
|
||||
TorusOptions::default(),
|
||||
)
|
||||
.unwrap();
|
||||
let s = Surface::from(made);
|
||||
let form = as_b_spline(&k, &s).unwrap().unwrap().into_owned();
|
||||
let own = AnalyticSupport::of(&s).unwrap();
|
||||
let m0 = own.parametrisation_margin(&k, &form).unwrap().unwrap();
|
||||
assert!(m0 <= 1.0, "own label: {m0}");
|
||||
assert!(own.build(&k, form.copy()).is_ok());
|
||||
let turned = |t: f64| {
|
||||
let AnalyticSupport::Torus {
|
||||
placement,
|
||||
major_radius,
|
||||
minor_radius,
|
||||
} = own
|
||||
else {
|
||||
panic!()
|
||||
};
|
||||
let (c, sn) = (t.cos(), t.sin());
|
||||
let (px, py) = (placement.x_axis(), placement.y_axis());
|
||||
AnalyticSupport::Torus {
|
||||
placement: Axis2Placement::from_stored(
|
||||
placement.origin(),
|
||||
px * c + py * sn,
|
||||
py * c - px * sn,
|
||||
),
|
||||
major_radius,
|
||||
minor_radius,
|
||||
}
|
||||
};
|
||||
let probe = 1e-11;
|
||||
let m1 = turned(probe)
|
||||
.parametrisation_margin(&k, &form)
|
||||
.unwrap()
|
||||
.unwrap();
|
||||
assert!(m1 > 1.0, "a turn of {probe} rad reads {m1}");
|
||||
let four = 4.0 * probe / m1;
|
||||
let m4 = turned(four)
|
||||
.parametrisation_margin(&k, &form)
|
||||
.unwrap()
|
||||
.unwrap();
|
||||
eprintln!("own {m0:.4}; turn {probe:e}: {m1:.3}; turn {four:e}: {m4:.3}");
|
||||
assert!(m4 > 1.0, "a turn of {four:e} rad reads {m4}");
|
||||
assert!(turned(four).build(&k, form).is_err());
|
||||
}
|
||||
+333
-30
@@ -1123,12 +1123,37 @@ impl AnalyticSupport {
|
||||
// a surface of revolution: the one the kernel recognises in the form, which must be
|
||||
// this one, take the form as its NURB and evaluate as the form at the form's own
|
||||
// parameters
|
||||
let not_form = || {
|
||||
SolidError::invalid_input(
|
||||
OP_BUILD,
|
||||
"a support's analytic surface is not its NURBS form's parametrisation",
|
||||
)
|
||||
};
|
||||
let surface = self.on_form(k, &form)?;
|
||||
if !evaluates_as(k, &surface, &form)? {
|
||||
return Err(not_form());
|
||||
}
|
||||
Ok(surface)
|
||||
}
|
||||
|
||||
/// The largest ratio, over the form's Greville points, of this support's analytic
|
||||
/// surface's distance from `form` to the rounding [`evaluates_as`] derives there (#3185):
|
||||
/// at most 1 where [`Self::build`]'s parametrisation test passes. The surface is the one
|
||||
/// `build` makes (the kernel's recognition of `form`, carrying this support's numbers);
|
||||
/// `lies_on` and the placement's axes are not checked here. `None` for a plane, or a
|
||||
/// surface the bound does not cover.
|
||||
///
|
||||
/// # Errors
|
||||
///
|
||||
/// `InvalidInput` where the kernel recognises no surface of this kind in `form`; those of
|
||||
/// the recognisers, the conversions and the evaluations.
|
||||
pub fn parametrisation_margin(&self, k: &Kernel, form: &BSplineSurface) -> Result<Option<f64>> {
|
||||
if matches!(self, Self::Plane { .. }) {
|
||||
return Ok(None);
|
||||
}
|
||||
let surface = self.on_form(k, form)?;
|
||||
evaluation_margin(k, &surface, form)
|
||||
}
|
||||
|
||||
/// The surface of revolution of this support's kind the kernel recognises in `form`,
|
||||
/// carrying this support's own numbers and `form` as its NURB ([`Self::build`] after its
|
||||
/// checks of the point set, before the parametrisation's).
|
||||
fn on_form(&self, k: &Kernel, form: &BSplineSurface) -> Result<Surface> {
|
||||
let form = form.copy();
|
||||
let mut surface = recognise(k, self.kind(), &form)?.ok_or_else(not_form)?;
|
||||
// the recognised surface must be of this support's kind; the numbers recognition
|
||||
// recomputed are not compared with the support's, because it recomputes them from
|
||||
@@ -1178,13 +1203,18 @@ impl AnalyticSupport {
|
||||
}
|
||||
_ => return Err(not_form()),
|
||||
}
|
||||
if !evaluates_as(k, &surface, &form)? {
|
||||
return Err(not_form());
|
||||
}
|
||||
Ok(surface)
|
||||
}
|
||||
}
|
||||
|
||||
/// [`AnalyticSupport::build`]'s refusal of a surface that is not its form's parametrisation.
|
||||
fn not_form() -> SolidError {
|
||||
SolidError::invalid_input(
|
||||
OP_BUILD,
|
||||
"a support's analytic surface is not its NURBS form's parametrisation",
|
||||
)
|
||||
}
|
||||
|
||||
/// A cylinder's or cone's recognised parametrisation carried by the support's own
|
||||
/// `placement`: its heights are measured along the axis from the placement's origin, so the
|
||||
/// recognised heights move by how far along the axis the recognised origin lies from it.
|
||||
@@ -1259,41 +1289,314 @@ fn recognise(k: &Kernel, kind: SupportKind, form: &BSplineSurface) -> Result<Opt
|
||||
}
|
||||
|
||||
/// Whether `surface` evaluates as `form` at the form's Greville parameters (each control
|
||||
/// point's), within the rounding of the two evaluations: `gamma_K` times the form's size,
|
||||
/// `K = 16 (p + q + 2)` for the rational evaluation of degrees `(p, q)` and the analytic
|
||||
/// one's conversion of the parameters and its trigonometry, a few ulps per operation.
|
||||
/// point's): at each, the analytic surface's point and the form's agree within the bound
|
||||
/// [`Rounding::at`] derives for that point, term by term below. The analytic side is the
|
||||
/// **STEP** evaluation of the converted parameters ([`Surface::convert_uv_from_nurbs_to_step`]
|
||||
/// and [`Surface::evaluate_step_point`]), which a cone, sphere, torus and surface of
|
||||
/// revolution compute from their placement, radii, semi-angle and generator.
|
||||
/// [`Surface::evaluate_point`] would not do: a sphere, a torus and a surface of revolution
|
||||
/// take it from their NURBS base, and that base is the `form` the caller has just installed,
|
||||
/// so it would compare the form with itself and read nothing of what the caller imposed
|
||||
/// (#2588).
|
||||
///
|
||||
/// The analytic side is the **STEP** evaluation of the converted parameters
|
||||
/// ([`Surface::convert_uv_from_nurbs_to_step`] and [`Surface::evaluate_step_point`]),
|
||||
/// which a cone, sphere, torus and surface of revolution compute from their placement,
|
||||
/// radii, semi-angle and generator. [`Surface::evaluate_point`] would not do: a sphere, a
|
||||
/// torus and a surface of revolution take it from their NURBS base, and that base is the
|
||||
/// `form` the caller has just installed, so it would compare the form with itself and read
|
||||
/// nothing of what the caller imposed (#2588; the cone already converted and evaluated
|
||||
/// analytically, which is why only the other kinds went unchecked).
|
||||
/// What the bound forgives and what it tests (#3185): it forgives the rounding of both
|
||||
/// evaluations and of the two parameter conversions, nothing else. It does not forgive the
|
||||
/// label's frame turned about the axis against the converter's, its axis, origin or radii,
|
||||
/// a form whose parameter is not the analytic surface's (a reparametrised arc), nor the
|
||||
/// form's distance from the analytic surface: those are what the comparison proves.
|
||||
///
|
||||
/// `u` is the unit roundoff, `gamma_n = n u / (1 - n u)`. Stated assumption: the libm
|
||||
/// `sin`, `cos`, `tan` and `atan2` are within one ulp (relative `2 u`).
|
||||
///
|
||||
/// 1. **The form** (`Surface::eval`, Cox-de Boor then two blends): each basis value is a
|
||||
/// nonnegative combination over `p` levels of at most 6 operations (two knot differences,
|
||||
/// a sum, a division, a product and a sum), so relative `gamma_6p`; the blends sum
|
||||
/// `q + 1` then `p + 1` products. With positive weights the homogeneous point's error is
|
||||
/// at most `gamma_(7 (p + q) + 4) W size` and the weight's `gamma_(7 (p + q) + 4) W`, so
|
||||
/// the Euclidean point's, after the division, `E_form = (2 gamma_(7 (p + q) + 4) + u)
|
||||
/// size` (`gamma_(7 (p + q) + 4) size` when polynomial), `size` the net's largest point.
|
||||
/// 2. **The STEP point**, `O + m (cos u X + sin u Y) + h Z`, `m` and `h` the generator's
|
||||
/// radial and axial parts (the torus `R + r cos v`, `r sin v`; the sphere `r cos v`,
|
||||
/// `r sin v`; the cone `r0 + v tan a`, `v`). An angle `x` in degrees goes to radians in
|
||||
/// two operations of the rounded `PI`: `gamma_3 |x|`; with one ulp of the cosine and sine,
|
||||
/// each is off by `dc(x) = 2 u + gamma_3 |x|`. Per coordinate (unit axes) the error is at
|
||||
/// most `gamma_3 (|O| + M + H) + M (gamma_2 + 2 dc(u)) + dm + dh`, `M`, `H` the sizes of
|
||||
/// `m`, `h` and `dm`, `dh` their errors: the torus `gamma_3 (R + r) + 2 r dc(v)`, the
|
||||
/// sphere `gamma_1 r + 2 r dc(v)`, the cone `gamma_2 (|r0| + |v t|) + |v| dt` (`dt =
|
||||
/// (1 + t^2) gamma_3 |a| + 2 u |t|` for the tangent of the half angle `a`); a general
|
||||
/// surface of revolution rotates its generator's point `G` (a B-spline: item 1's bound of
|
||||
/// its own net, `E_gen`) about the axis, `gamma_8 (|O| + 3 D) + 2 D dc(u) + E_gen`, `D`
|
||||
/// the net's reach from `O`. `E_step` is `sqrt 3` times the coordinate bound.
|
||||
/// 3. **The generator's parameter**: a circle (torus, sphere) inverts its NURBS point by
|
||||
/// `atan2` about its centre `c`: at its radius `r` the error is that point's (`E_form`)
|
||||
/// plus `gamma_4 (size + |c|)` and `gamma_4 r |v|` of the angle's own rounding; a line
|
||||
/// (cone) maps it affinely in 6 operations, `gamma_6 2 (|v0| + |v1|)` along the axis and
|
||||
/// its slope, times `1 + |t|`; a B-spline generator's parameter is its NURBS's (none).
|
||||
/// 4. **The sweep's parameter**: the polar converter evaluates its circle's NURBS at the
|
||||
/// parameter (item 1's bound of that NURBS's own net, `E_polar`) and inverts the point by
|
||||
/// `atan2` about the circle's centre `c` in its frame: at the circle's radius `rc` the
|
||||
/// angle is off by `(E_polar + gamma_4 (rc + 2 |c|)) / rc`, plus `gamma_4 |s|` of the
|
||||
/// angle's own rounding and the degrees' period (`gamma_1 2 pi`), times the STEP point's
|
||||
/// distance from that axis (the lever). The converter's ends are its domain's, which
|
||||
/// recognition takes from the form's own end meridians (#3199), so they are the inverted
|
||||
/// ends to the same rounding.
|
||||
///
|
||||
/// A surface of another kind, or a general surface of revolution whose generator is not a
|
||||
/// B-spline, has no bound here and is not proved.
|
||||
///
|
||||
/// # Errors
|
||||
///
|
||||
/// Those of the parameter conversion and the evaluations and of the form's knots.
|
||||
fn evaluates_as(k: &Kernel, surface: &Surface, form: &BSplineSurface) -> Result<bool> {
|
||||
Ok(evaluation_margin(k, surface, form)?.is_some_and(|m| m <= 1.0))
|
||||
}
|
||||
|
||||
/// [`evaluates_as`]'s largest ratio of distance to bound over the Greville points (infinite
|
||||
/// for a distance that is not a number); `None` for a surface the bound does not cover.
|
||||
fn evaluation_margin(k: &Kernel, surface: &Surface, form: &BSplineSurface) -> Result<Option<f64>> {
|
||||
let (net, p, q) = euclidean_net(form)?;
|
||||
let size = net.iter().flatten().fold(0.0_f64, |a, x| a.max(x.length()));
|
||||
let tol = gamma(16 * (p + q + 2)) * size;
|
||||
let Some(r) = Rounding::of(surface, form.is_rational(), p, q, size, &net)? else {
|
||||
return Ok(None);
|
||||
};
|
||||
let mut worst = 0.0_f64;
|
||||
for u in greville(form, true)? {
|
||||
for v in greville(form, false)? {
|
||||
let uv = Vec2::new(u, v);
|
||||
let step_uv = surface.convert_uv_from_nurbs_to_step(k, uv)?;
|
||||
let off =
|
||||
(surface.evaluate_step_point(k, step_uv)? - form.evaluate_point(uv)?).length();
|
||||
if !matches!(
|
||||
off.partial_cmp(&tol),
|
||||
Some(core::cmp::Ordering::Less | core::cmp::Ordering::Equal)
|
||||
) {
|
||||
return Ok(false);
|
||||
}
|
||||
let s = surface.evaluate_step_point(k, step_uv)?;
|
||||
let f = form.evaluate_point(uv)?;
|
||||
let ratio = (s - f).length() / r.at(step_uv, s, f);
|
||||
worst = if ratio.is_nan() {
|
||||
f64::INFINITY
|
||||
} else {
|
||||
worst.max(ratio)
|
||||
};
|
||||
}
|
||||
}
|
||||
Ok(true)
|
||||
Ok(Some(worst))
|
||||
}
|
||||
|
||||
/// A surface of revolution's generator, for [`Rounding`]'s item 2 and 3.
|
||||
#[derive(Debug, Clone, Copy)]
|
||||
enum Generator {
|
||||
/// `R + r cos v`, `r sin v` about the circle of centre `c`.
|
||||
Torus { major: f64, minor: f64, c: f64 },
|
||||
/// `r cos v`, `r sin v` about the circle of centre `c`.
|
||||
Sphere { radius: f64, c: f64 },
|
||||
/// `r0 + v t`, `v` (the half angle `a` in radians), its STEP `v` over `[v0, v1]`.
|
||||
Cone {
|
||||
r0: f64,
|
||||
t: f64,
|
||||
a: f64,
|
||||
v0: f64,
|
||||
v1: f64,
|
||||
},
|
||||
/// A B-spline generator, its own evaluation's bound.
|
||||
Spline { e_gen: f64 },
|
||||
}
|
||||
|
||||
/// The pieces of [`evaluates_as`]'s bound that do not depend on the point.
|
||||
#[derive(Debug, Clone, Copy)]
|
||||
struct Rounding {
|
||||
e_form: f64,
|
||||
origin: Vec3,
|
||||
/// The net's reach from the origin.
|
||||
reach: f64,
|
||||
generator: Generator,
|
||||
/// The sweep converter's circle: centre, axes.
|
||||
polar: (Vec3, Vec3, Vec3, f64, f64),
|
||||
}
|
||||
|
||||
/// `|v|`'s largest coordinate.
|
||||
fn max_coord(v: Vec3) -> f64 {
|
||||
v.x.abs().max(v.y.abs()).max(v.z.abs())
|
||||
}
|
||||
|
||||
/// A rational (or polynomial) evaluation's bound of degrees summing to `n` over points of
|
||||
/// size `size` ([`evaluates_as`] item 1).
|
||||
fn eval_bound(rational: bool, n: usize, size: f64) -> f64 {
|
||||
let g = gamma(7 * n + 4);
|
||||
if rational {
|
||||
(2.0 * g + crate::nurbs::math::UNIT) * size
|
||||
} else {
|
||||
g * size
|
||||
}
|
||||
}
|
||||
|
||||
impl Rounding {
|
||||
/// The bound's data for `surface` (a cone, sphere, torus or surface of revolution with a
|
||||
/// B-spline generator) on its form of degrees `(p, q)`, net `net`; `None` for any other.
|
||||
fn of(
|
||||
surface: &Surface,
|
||||
rational: bool,
|
||||
p: usize,
|
||||
q: usize,
|
||||
size: f64,
|
||||
net: &[Vec<Vec3>],
|
||||
) -> Result<Option<Self>> {
|
||||
use crate::curve::curve::CurveKind;
|
||||
let e_form = eval_bound(rational, p + q, size);
|
||||
let circle_centre = |rev: &crate::surface::surf_of_revolution::SurfOfRevolution| match rev
|
||||
.gen_curve
|
||||
.as_ref()
|
||||
.map(crate::curve::curve::Curve::kind)
|
||||
{
|
||||
Some(CurveKind::Ellipse(e)) => Some(e.position.origin().length()),
|
||||
_ => None,
|
||||
};
|
||||
let (rev, generator) = match surface.kind() {
|
||||
SurfaceKind::Torus(t) => {
|
||||
let Some(c) = circle_centre(&t.base) else {
|
||||
return Ok(None);
|
||||
};
|
||||
(
|
||||
&t.base,
|
||||
Generator::Torus {
|
||||
major: t.major_radius,
|
||||
minor: t.minor_radius,
|
||||
c,
|
||||
},
|
||||
)
|
||||
}
|
||||
SurfaceKind::Sphere(sp) => {
|
||||
let Some(c) = circle_centre(&sp.base) else {
|
||||
return Ok(None);
|
||||
};
|
||||
(
|
||||
&sp.base,
|
||||
Generator::Sphere {
|
||||
radius: sp.radius,
|
||||
c,
|
||||
},
|
||||
)
|
||||
}
|
||||
SurfaceKind::Cone(c) => {
|
||||
let a = c.semi_angle_deg * PI / 180.0;
|
||||
let d = c.rev.anal_uv_domain;
|
||||
(
|
||||
&c.rev,
|
||||
Generator::Cone {
|
||||
r0: c.base_radius,
|
||||
t: c.tan_semi_angle(),
|
||||
a,
|
||||
v0: d.v_min(),
|
||||
v1: d.v_max(),
|
||||
},
|
||||
)
|
||||
}
|
||||
SurfaceKind::SurfOfRevolution(rev) => {
|
||||
let Some(CurveKind::BSpline(g)) =
|
||||
rev.gen_curve.as_ref().map(crate::curve::curve::Curve::kind)
|
||||
else {
|
||||
return Ok(None);
|
||||
};
|
||||
let pts = g.control_points()?;
|
||||
let gsize = pts.iter().fold(0.0_f64, |m, c| {
|
||||
let w = if c.w == NOW { 1.0 } else { c.w };
|
||||
m.max((Vec3::new(c.x, c.y, c.z) * (1.0 / w)).length())
|
||||
});
|
||||
let rational = pts.iter().any(|c| c.w != NOW);
|
||||
(
|
||||
&**rev,
|
||||
Generator::Spline {
|
||||
e_gen: eval_bound(rational, g.degree()?, gsize),
|
||||
},
|
||||
)
|
||||
}
|
||||
_ => return Ok(None),
|
||||
};
|
||||
let Some(pc) = rev.polar_curve() else {
|
||||
return Ok(None);
|
||||
};
|
||||
let e_polar = {
|
||||
let pts = pc.base().control_points()?;
|
||||
let psize = pts.iter().fold(0.0_f64, |m, c| {
|
||||
let w = if c.w == NOW { 1.0 } else { c.w };
|
||||
m.max((Vec3::new(c.x, c.y, c.z) * (1.0 / w)).length())
|
||||
});
|
||||
eval_bound(pts.iter().any(|c| c.w != NOW), pc.base().degree()?, psize)
|
||||
};
|
||||
let origin = rev.position.origin();
|
||||
let reach = net
|
||||
.iter()
|
||||
.flatten()
|
||||
.fold(0.0_f64, |m, x| m.max((*x - origin).length()));
|
||||
Ok(Some(Self {
|
||||
e_form,
|
||||
origin,
|
||||
reach,
|
||||
generator,
|
||||
polar: (
|
||||
pc.position.origin(),
|
||||
pc.position.x_axis(),
|
||||
pc.position.y_axis(),
|
||||
pc.radius_x,
|
||||
e_polar,
|
||||
),
|
||||
}))
|
||||
}
|
||||
|
||||
/// The bound at the STEP parameters `step` (degrees, and the generator's), the STEP
|
||||
/// point `s` and the form's point `f` ([`evaluates_as`] items 1 to 4).
|
||||
fn at(&self, step: Vec2, s: Vec3, f: Vec3) -> f64 {
|
||||
let u = crate::nurbs::math::UNIT;
|
||||
let rad = |deg: f64| deg * PI / 180.0;
|
||||
let dc = |deg: f64| 2.0 * u + gamma(3) * rad(deg).abs();
|
||||
let o = max_coord(self.origin);
|
||||
let size_f = f.length();
|
||||
// item 2 (per coordinate) and item 3 (as a length)
|
||||
let (coord, generator) = match self.generator {
|
||||
Generator::Torus { major, minor, c } => {
|
||||
let (m, h) = (major + minor, minor);
|
||||
let dm = gamma(3) * (major + minor) + 2.0 * minor * dc(step.y);
|
||||
let dh = gamma(1) * minor + 2.0 * minor * dc(step.y);
|
||||
(
|
||||
gamma(3) * (o + m + h) + m * (gamma(2) + 2.0 * dc(step.x)) + dm + dh,
|
||||
self.e_form + gamma(4) * (size_f + c) + gamma(4) * minor * rad(step.y).abs(),
|
||||
)
|
||||
}
|
||||
Generator::Sphere { radius, c } => {
|
||||
let dm = gamma(1) * radius + 2.0 * radius * dc(step.y);
|
||||
(
|
||||
gamma(3) * (o + 2.0 * radius)
|
||||
+ radius * (gamma(2) + 2.0 * dc(step.x))
|
||||
+ 2.0 * dm,
|
||||
self.e_form + gamma(4) * (size_f + c) + gamma(4) * radius * rad(step.y).abs(),
|
||||
)
|
||||
}
|
||||
Generator::Cone { r0, t, a, v0, v1 } => {
|
||||
let v = step.y;
|
||||
let m = r0.abs() + (v * t).abs();
|
||||
let dt = (1.0 + t * t) * gamma(3) * a.abs() + 2.0 * u * t.abs();
|
||||
let dm = gamma(2) * m + v.abs() * dt;
|
||||
(
|
||||
gamma(3) * (o + m + v.abs()) + m * (gamma(2) + 2.0 * dc(step.x)) + dm,
|
||||
gamma(6) * 2.0 * (v0.abs() + v1.abs()) * (1.0 + t.abs()),
|
||||
)
|
||||
}
|
||||
Generator::Spline { e_gen } => (
|
||||
gamma(8) * (o + 3.0 * self.reach) + 2.0 * self.reach * dc(step.x) + e_gen,
|
||||
0.0,
|
||||
),
|
||||
};
|
||||
let e_step = 3.0_f64.sqrt() * coord;
|
||||
// item 4: the polar converter's inversion of its own circle's point
|
||||
let (c, x, y, rc, e_polar) = self.polar;
|
||||
let z = x.cross(y);
|
||||
let along = |v: Vec3| v - z * v.dot(z);
|
||||
let lever = along(s - c).length();
|
||||
let sweep = if rc > 0.0 {
|
||||
lever
|
||||
* ((e_polar + gamma(4) * (rc + 2.0 * c.length())) / rc
|
||||
+ gamma(4) * rad(step.x).abs()
|
||||
+ gamma(1) * 2.0 * PI)
|
||||
} else {
|
||||
// a circle of no radius has no angle: any, at the STEP point's lever
|
||||
lever * PI
|
||||
};
|
||||
self.e_form + e_step + generator + sweep
|
||||
}
|
||||
}
|
||||
|
||||
const OP_DERIVE: &str = "surface::implicit::AnalyticSupport::derive";
|
||||
|
||||
+99
-6
@@ -1061,7 +1061,6 @@ impl Torus {
|
||||
let ref_frame = ref_iso.reference_frame.ok_or_else(missing)?;
|
||||
let other_frame = other_iso.reference_frame.ok_or_else(missing)?;
|
||||
let min_rad = other_iso.radius.ok_or_else(missing)?;
|
||||
let sweep_ivl = ref_iso.angles.ok_or_else(missing)?;
|
||||
let common_point = other_iso.common_line_point.ok_or_else(missing)?;
|
||||
let common_vec = other_iso.common_line_vec.ok_or_else(missing)?;
|
||||
|
||||
@@ -1069,6 +1068,17 @@ impl Torus {
|
||||
let torus_origin =
|
||||
common_point + common_vec.dot(other_frame.origin() - common_point) * common_vec;
|
||||
|
||||
// not upstream (#3199): the x axis and the sweep interval from the form's own first
|
||||
// sweep row's ends about the axis (`sweep_frame`), not from `is_arc`'s fit of it
|
||||
let Some((torus_x_axis, torus_y_axis, sweep_ivl)) = sweep_frame(
|
||||
torus_origin,
|
||||
ref_frame.z_axis(),
|
||||
fitted_sides(&ref_frame, ref_iso.angles.ok_or_else(missing)?),
|
||||
meridians(test, swap_uv)?,
|
||||
) else {
|
||||
return Ok(None);
|
||||
};
|
||||
|
||||
// get torus orientation and size
|
||||
// not upstream (WR-501, #1470, opt-in; tangency-labs/smlib@6136a04): the major
|
||||
// radius runs from the torus origin, not from the first sweep circle's centre,
|
||||
@@ -1079,9 +1089,7 @@ impl Torus {
|
||||
other_frame.origin() - ref_frame.origin()
|
||||
};
|
||||
let maj_rad = major_vec.length();
|
||||
let torus_x_axis = ref_frame.x_axis(); // the frame where all circles are parallel
|
||||
let torus_y_axis = ref_frame.y_axis();
|
||||
let torus_z_axis = ref_frame.z_axis();
|
||||
let torus_z_axis = ref_frame.z_axis(); // the frame where all circles are parallel
|
||||
|
||||
// torus must have positive radii
|
||||
if min_rad <= EFF_ZERO || maj_rad <= EFF_ZERO {
|
||||
@@ -1294,7 +1302,6 @@ impl Torus {
|
||||
let mut ref_frame = ref_iso.reference_frame.ok_or_else(missing)?;
|
||||
let other_frame = other_iso.reference_frame.ok_or_else(missing)?;
|
||||
let min_rad = other_iso.radius.ok_or_else(missing)?;
|
||||
let sweep_ivl = ref_iso.angles.ok_or_else(missing)?;
|
||||
let common_point = other_iso.common_line_point.ok_or_else(missing)?;
|
||||
let common_vec = other_iso.common_line_vec.ok_or_else(missing)?;
|
||||
|
||||
@@ -1302,8 +1309,18 @@ impl Torus {
|
||||
let torus_origin =
|
||||
common_point + common_vec.dot(other_frame.origin() - common_point) * common_vec;
|
||||
|
||||
// not upstream (#3199): the x axis and the sweep interval from the form's own first
|
||||
// sweep row's ends about the axis (`sweep_frame`), not from `is_arc`'s fit of it
|
||||
let (x_axis, y_axis, sweep_ivl) = sweep_frame(
|
||||
torus_origin,
|
||||
ref_frame.z_axis(),
|
||||
fitted_sides(&ref_frame, ref_iso.angles.ok_or_else(missing)?),
|
||||
meridians(srf, swap_uv)?,
|
||||
)
|
||||
.ok_or_else(|| SolidError::sm_err(OP_REBUILD))?;
|
||||
|
||||
// Install the origin into pRefFrame. // [B328]
|
||||
let _ = ref_frame.set_canonical(torus_origin, ref_frame.x_axis(), ref_frame.y_axis());
|
||||
let _ = ref_frame.set_canonical(torus_origin, x_axis, y_axis);
|
||||
|
||||
// get torus orientation and size
|
||||
let major_vec = other_frame.origin() - ref_frame.origin();
|
||||
@@ -1587,6 +1604,82 @@ impl Revolved for Torus {
|
||||
|
||||
/// `SmTorus::CreateIsoParametricCurve` of the surface of revolution `rev` of a torus
|
||||
/// of `shape` (also called from the base class's rebuilds, where it is virtual).
|
||||
/// The control points, Euclidean, of the form's first and last meridians: the columns of its
|
||||
/// net at the sweep's two ends (the sweep is along u, or along v with `swap_uv`). A clamped
|
||||
/// end's column is that boundary curve's own control polygon, and a meridian lies in a plane
|
||||
/// through the axis, so every one of its control points has the meridian's angle about the
|
||||
/// axis.
|
||||
fn meridians(srf: &BSplineSurface, swap_uv: bool) -> Result<[Vec<Point3>; 2]> {
|
||||
let nurb = srf
|
||||
.gw_nurb()
|
||||
.ok_or_else(|| SolidError::null_pointer(OP_NURB_TORUS))?;
|
||||
let (n, m) = (nurb.n(), nurb.m());
|
||||
let column = |end: bool| -> Result<Vec<Point3>> {
|
||||
if swap_uv {
|
||||
let j = if end { m } else { 0 };
|
||||
(0..=n).map(|i| Ok(nurb.cpt(i, j)?.to_euclid())).collect()
|
||||
} else {
|
||||
let i = if end { n } else { 0 };
|
||||
(0..=m).map(|j| Ok(nurb.cpt(i, j)?.to_euclid())).collect()
|
||||
}
|
||||
};
|
||||
Ok([column(false)?, column(true)?])
|
||||
}
|
||||
|
||||
/// The directions, about the axis, of the fitted first sweep row's ends (`is_arc`'s frame and
|
||||
/// angles): [`sweep_frame`] reads only their sides from them.
|
||||
fn fitted_sides(frame: &Axis2Placement, angles: Extent1d) -> [Vec3; 2] {
|
||||
let at = |deg: f64| {
|
||||
let a = deg * PI / 180.0;
|
||||
a.cos() * frame.x_axis() + a.sin() * frame.y_axis()
|
||||
};
|
||||
[at(angles.min()), at(angles.max())]
|
||||
}
|
||||
|
||||
/// Not upstream (#3199): a recognised torus's x axis, y axis and sweep interval (degrees) from
|
||||
/// its form's first and last meridians (`meridians`: their control points, exact data), about
|
||||
/// the axis through `origin` along the unit `z`. A meridian lies in a plane through the axis,
|
||||
/// so its control points share its angle on the side the tube's centre is on, and the
|
||||
/// opposite angle past the axis (a spindle torus's tube crosses it); `sides` (the fitted
|
||||
/// row's end directions, [`fitted_sides`]) says which side that is, and nothing else. Each
|
||||
/// meridian's angle is read at its control point farthest from the axis on that side: x is
|
||||
/// the first meridian's direction square to the axis (so the interval starts at 0), y `z x
|
||||
/// x`, and the interval's end the last meridian's angle from x about z (`atan2` in that
|
||||
/// frame), a closed sweep (within `SM_EFF_ZERO_SQRT` of 0) at 360, as `is_arc` maps its end.
|
||||
/// `is_arc` fitted the first sweep row's circle through three samples, so on a short arc its
|
||||
/// centre, and with it the x axis and the end angle, was off the axis by its conditioning
|
||||
/// (1e-14 rad on #3185's grooves); the converter centres the sweep circle on the axis, so its
|
||||
/// ends and its inside now agree to rounding. `None` when a meridian has no control point off
|
||||
/// the axis on its side.
|
||||
fn sweep_frame(
|
||||
origin: Point3,
|
||||
z: Vec3,
|
||||
sides: [Vec3; 2],
|
||||
[first, last]: [Vec<Point3>; 2],
|
||||
) -> Option<(Vec3, Vec3, Extent1d)> {
|
||||
let radial = |p: &Point3| {
|
||||
let s = *p - origin;
|
||||
s - s.dot(z) * z
|
||||
};
|
||||
let farthest = |ps: &[Point3], side: Vec3| {
|
||||
ps.iter()
|
||||
.map(radial)
|
||||
.filter(|r| r.dot(side) > 0.0)
|
||||
.max_by(|a, b| a.length().total_cmp(&b.length()))
|
||||
};
|
||||
let [first_side, last_side] = sides;
|
||||
let x = farthest(&first, first_side)?.unitize().ok()?;
|
||||
let y = z.cross(x);
|
||||
let e = farthest(&last, last_side)?;
|
||||
let mut end_deg = e.dot(y).atan2(e.dot(x)) * 180.0 / PI;
|
||||
if end_deg.abs() < EFF_ZERO_SQRT {
|
||||
end_deg = 360.0;
|
||||
} else if end_deg < 0.0 {
|
||||
end_deg += 360.0;
|
||||
}
|
||||
Some((x, y, Extent1d::from_min_max(0.0, end_deg)))
|
||||
}
|
||||
|
||||
fn iso_curve(
|
||||
rev: &SurfOfRevolution,
|
||||
shape: Shape,
|
||||
|
||||
Reference in New Issue
Block a user