diff --git a/crates/doc/src/clipboard.rs b/crates/doc/src/clipboard.rs index b6dc2e17..ab98a7b6 100644 --- a/crates/doc/src/clipboard.rs +++ b/crates/doc/src/clipboard.rs @@ -194,6 +194,7 @@ impl crate::Feature { | K::SurfaceTrim { plane, .. } => checked_plane(plane, 0)?, K::Draft { neutral, .. } => checked_plane(neutral, 0)?, K::ReplaceFace { target, .. } => checked_plane(target, 0)?, + K::SurfaceExtend { up_to: Some(crate::ExtendTo::Plane(plane)), .. } => checked_plane(plane, 0)?, K::ConstructionAxis { def } => { for r in def.refs() { checked_geo(r, 0)?; @@ -248,6 +249,7 @@ impl crate::Feature { | K::SurfaceTrim { plane, .. } => remap_plane(plane, ids, names, bodies, 0)?, K::Draft { neutral, .. } => remap_plane(neutral, ids, names, bodies, 0)?, K::ReplaceFace { target, .. } => remap_plane(target, ids, names, bodies, 0)?, + K::SurfaceExtend { up_to: Some(crate::ExtendTo::Plane(plane)), .. } => remap_plane(plane, ids, names, bodies, 0)?, K::ConstructionAxis { def } => { for r in def.refs_mut() { remap_geo(r, ids, names, bodies, 0)?; diff --git a/crates/doc/src/document.rs b/crates/doc/src/document.rs index ea1b7399..076513c6 100644 --- a/crates/doc/src/document.rs +++ b/crates/doc/src/document.rs @@ -1137,11 +1137,20 @@ pub enum FeatureKind { #[serde(default, skip_serializing_if = "Option::is_none")] face: Option, }, - /// A planar surface body extended outward by `distance` past the edges through `edges`. + /// A surface body's boundary edges through `edges` extended outward: by `distance`, or up + /// to `up_to`; the surface continued (`mode` natural) or a run-off face along each edge + /// (tangent, curvature). `neighbours`: the edges next to an extended edge are lengthened to + /// meet it (else the extension's sides run square to the edge). SurfaceExtend { body: String, edges: Vec, distance: String, + #[serde(default, skip_serializing_if = "ExtendMode::is_natural")] + mode: ExtendMode, + #[serde(default, skip_serializing_if = "Option::is_none")] + up_to: Option, + #[serde(default = "yes", skip_serializing_if = "is_yes")] + neighbours: bool, }, /// A surface body trimmed by a plane, or by a face of a body (`tool`), keeping the side of /// `keep`. @@ -1189,6 +1198,37 @@ pub enum FeatureKind { }, } +/// How Extend continues a surface past its edges. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Serialize, Deserialize)] +#[serde(rename_all = "snake_case")] +pub enum ExtendMode { + /// The surface itself, continued. + #[default] + Natural, + /// A ruled run-off in the surface's tangent plane (G1). + Tangent, + /// A run-off that keeps the surface's curvature across the edge (G2). + Curvature, +} + +impl ExtendMode { + pub fn is_natural(&self) -> bool { + *self == ExtendMode::Natural + } +} + +/// What Extend extends up to instead of a distance. +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] +#[serde(rename_all = "snake_case")] +pub enum ExtendTo { + /// A plane. + Plane(PlaneRef), + /// The surface of a body's face (extended). + Face(FaceAt), + /// The surface's own bounds (natural only): a trim undone. + Bounds, +} + /// How a fillet's size is given. #[derive(Clone, Debug, Default, PartialEq, Serialize, Deserialize)] #[serde(tag = "type", rename_all = "snake_case")] diff --git a/crates/doc/src/eval_surface.rs b/crates/doc/src/eval_surface.rs index e0cc4fff..8ee61e94 100644 --- a/crates/doc/src/eval_surface.rs +++ b/crates/doc/src/eval_surface.rs @@ -78,11 +78,30 @@ pub(super) fn eval(doc: &Document, vals: &BTreeMap, f: &Feature, } Ok(()) } - FeatureKind::SurfaceExtend { body, edges, distance } => { - let d = val(vals, distance, Kind::Length)?; + FeatureKind::SurfaceExtend { body, edges, distance, mode, up_to, neighbours } => { let i = st.bodies.iter().position(|b| &b.name == body).ok_or_else(|| DocError::Unknown(format!("body `{body}`")))?; let Some(mb) = st.bodies.get(i).cloned() else { return Err(DocError::Invalid("body".into())) }; - let out = kernel::extend(&mb.body, edges, d)?; + let mode = match mode { + crate::ExtendMode::Natural => kernel::ExtendMode::Natural, + crate::ExtendMode::Tangent => kernel::ExtendMode::Tangent, + crate::ExtendMode::Curvature => kernel::ExtendMode::Curvature, + }; + let at = |limit| kernel::ExtendOptions { mode, limit, neighbours: *neighbours }; + let out = match up_to { + None => kernel::extend_surface(&mb.body, edges, &at(kernel::ExtendLimit::Distance(val(vals, distance, Kind::Length)?)))?, + Some(crate::ExtendTo::Bounds) => kernel::extend_surface(&mb.body, edges, &at(kernel::ExtendLimit::Bounds))?, + Some(crate::ExtendTo::Plane(plane)) => { + let tool = kernel::SplitTool::Plane(doc.resolve_plane_in(vals, plane, 0, Some(st))?); + kernel::extend_surface(&mb.body, edges, &at(kernel::ExtendLimit::UpTo(&tool)))? + } + Some(crate::ExtendTo::Face(t)) => { + let tb = st.body(&t.body).ok_or_else(|| DocError::Unknown(format!("body `{}`", t.body)))?; + let face = crate::appearance::face_index_at(tb, t.point) + .ok_or_else(|| DocError::Invalid(format!("no face of `{}` at the target point", t.body)))?; + let tool = kernel::SplitTool::Face { body: &tb.body, face }; + kernel::extend_surface(&mb.body, edges, &at(kernel::ExtendLimit::UpTo(&tool)))? + } + }; if let Some(slot) = st.bodies.get_mut(i) { *slot = ModelBody::new(mb.name, out, mb.feature); } diff --git a/crates/doc/src/feature_geometry.rs b/crates/doc/src/feature_geometry.rs index 37133bec..8d1aad92 100644 --- a/crates/doc/src/feature_geometry.rs +++ b/crates/doc/src/feature_geometry.rs @@ -139,11 +139,21 @@ macro_rules! geometry_walk { f(Pos(start)); } } - FeatureKind::Chamfer { edges, .. } | FeatureKind::SheetHem { edges, .. } | FeatureKind::SurfaceExtend { edges, .. } => { + FeatureKind::Chamfer { edges, .. } | FeatureKind::SheetHem { edges, .. } => { for e in edges { f(Pos(e)); } } + FeatureKind::SurfaceExtend { edges, up_to, .. } => { + for e in edges { + f(Pos(e)); + } + match up_to { + Some(crate::ExtendTo::Plane(plane)) => f(Plane(plane)), + Some(crate::ExtendTo::Face(t)) => f(Pos(& $($m)? t.point)), + Some(crate::ExtendTo::Bounds) | None => {} + } + } FeatureKind::Shell { faces, .. } | FeatureKind::OffsetFace { faces, .. } | FeatureKind::DeleteFace { faces, .. } diff --git a/crates/doc/src/feature_refs.rs b/crates/doc/src/feature_refs.rs index ecaa2879..9a8623a0 100644 --- a/crates/doc/src/feature_refs.rs +++ b/crates/doc/src/feature_refs.rs @@ -175,8 +175,15 @@ macro_rules! feature_walk { } } FeatureKind::SurfaceUntrim { face, .. } => f(B(& $($m)? face.body)), + FeatureKind::SurfaceExtend { body, up_to, .. } => { + f(B(body)); + match up_to { + Some(crate::ExtendTo::Plane(plane)) => $plane(plane, f), + Some(crate::ExtendTo::Face(t)) => f(B(& $($m)? t.body)), + Some(crate::ExtendTo::Bounds) | None => {} + } + } FeatureKind::SheetConvert { body, .. } - | FeatureKind::SurfaceExtend { body, .. } | FeatureKind::SurfaceRuled { body, .. } | FeatureKind::Unstitch { body, .. } | FeatureKind::ReverseNormal { body, .. } diff --git a/crates/doc/src/feature_walk_tests/fixtures.rs b/crates/doc/src/feature_walk_tests/fixtures.rs index 7aa70478..e6f79fc7 100644 --- a/crates/doc/src/feature_walk_tests/fixtures.rs +++ b/crates/doc/src/feature_walk_tests/fixtures.rs @@ -153,6 +153,14 @@ pub(super) fn fixtures() -> Vec<(&'static str, Value)> { ("reverse_normal", json!({"type": "reverse_normal", "body": "B"})), ("surface_untrim", json!({"type": "surface_untrim", "face": {"body": "B", "point": o}, "holes": [o], "outer": true})), ("surface_extend", json!({"type": "surface_extend", "body": "B", "edges": [o], "distance": "extend.distance"})), + ( + "surface_extend_to_plane", + json!({"type": "surface_extend", "body": "B", "edges": [o], "distance": "extend.distance", "mode": "curvature", "up_to": {"plane": plane("extend")}}), + ), + ( + "surface_extend_to_face", + json!({"type": "surface_extend", "body": "B", "edges": [o], "distance": "extend.distance", "up_to": {"face": {"body": "B", "point": o}}, "neighbours": false}), + ), ( "surface_ruled", json!({"type": "surface_ruled", "body": "B", "edges": [o], "distance": "ruled.distance", "mode": "angle", "angle": "ruled.angle", diff --git a/crates/doc/src/feature_walk_tests/golden_exprs.txt b/crates/doc/src/feature_walk_tests/golden_exprs.txt index 2fef4bc3..d4cb20e5 100644 --- a/crates/doc/src/feature_walk_tests/golden_exprs.txt +++ b/crates/doc/src/feature_walk_tests/golden_exprs.txt @@ -330,6 +330,16 @@ surface_extend inputs: Distance = extend.distance (Length) expressions: extend.distance +surface_extend_to_plane + inputs: + Distance = extend.distance (Length) + Angle = extend.angle (Angle) + Offset = extend.offset (Length) + expressions: extend.distance, extend.angle, extend.offset +surface_extend_to_face + inputs: + Distance = extend.distance (Length) + expressions: extend.distance surface_ruled inputs: Distance = ruled.distance (Length) diff --git a/crates/doc/src/feature_walk_tests/golden_geometry.txt b/crates/doc/src/feature_walk_tests/golden_geometry.txt index 98569900..c133c5c3 100644 --- a/crates/doc/src/feature_walk_tests/golden_geometry.txt +++ b/crates/doc/src/feature_walk_tests/golden_geometry.txt @@ -155,6 +155,10 @@ surface_untrim face points: [0, 0, 0] surface_extend moves: edges[0][0], edges[0][1], edges[0][2] +surface_extend_to_plane + moves: edges[0][0], edges[0][1], edges[0][2], up_to.plane.axis_dir[0], up_to.plane.axis_dir[1], up_to.plane.axis_origin[0], up_to.plane.axis_origin[1], up_to.plane.axis_origin[2] +surface_extend_to_face + moves: edges[0][0], edges[0][1], edges[0][2], up_to.face.point[0], up_to.face.point[1], up_to.face.point[2] surface_ruled moves: direction[0], direction[2], edges[0][0], edges[0][1], edges[0][2], face[1], face[2] face points: [1, 0, 0] diff --git a/crates/doc/src/feature_walk_tests/golden_refs.txt b/crates/doc/src/feature_walk_tests/golden_refs.txt index 39a698d5..de2068ef 100644 --- a/crates/doc/src/feature_walk_tests/golden_refs.txt +++ b/crates/doc/src/feature_walk_tests/golden_refs.txt @@ -172,6 +172,11 @@ surface_untrim remap_bodies: face.body surface_extend remap_bodies: body +surface_extend_to_plane + remap_refs: up_to.plane.base.base.name + remap_bodies: body +surface_extend_to_face + remap_bodies: body, up_to.face.body surface_ruled remap_bodies: body surface_trim diff --git a/crates/doc/src/format.rs b/crates/doc/src/format.rs index a0091aa9..eced0442 100644 --- a/crates/doc/src/format.rs +++ b/crates/doc/src/format.rs @@ -173,6 +173,12 @@ //! A `SurfaceRuled` feature (Ruled Surface: `body`, `edges` points, `distance`, optional //! `mode`, `angle`, `direction`, `flip`, `corners` and `face`) is also first saved in 21; //! 0.11.0 asks for an update instead of meeting a feature it can't read. +//! Under the same version, a `SurfaceExtend` can say how it continues the surface (`mode`: +//! `tangent` or `curvature`; `natural` when absent), extend up to something instead of its +//! `distance` (`up_to`: a `plane`, a `face` as `body` and `point`, or `bounds`) and keep the +//! neighbouring edges as they are (`neighbours: false`). An Extend without them builds the +//! planar extension it always did and saves as before; 0.11.0 asks for an update instead of +//! dropping them and extending a plane by the distance. use serde_json::Value; diff --git a/crates/doc/src/params.rs b/crates/doc/src/params.rs index f35812dd..5d28ba87 100644 --- a/crates/doc/src/params.rs +++ b/crates/doc/src/params.rs @@ -353,7 +353,12 @@ macro_rules! feature_walk { FeatureKind::MeshConvert { tolerance, .. } => f("Tolerance", tolerance, L, I), FeatureKind::Thicken { thickness, .. } => f("Thickness", thickness, L, I), FeatureKind::SurfaceTrim { plane, .. } => $plane(plane, f), - FeatureKind::SurfaceExtend { distance, .. } => f("Distance", distance, L, I), + FeatureKind::SurfaceExtend { distance, up_to, .. } => { + f("Distance", distance, L, I); + if let Some(crate::ExtendTo::Plane(plane)) = up_to { + $plane(plane, f); + } + } FeatureKind::SurfaceRuled { distance, angle, .. } => { f("Distance", distance, L, I); if let Some(a) = angle { diff --git a/crates/engine/src/cmd/surface.rs b/crates/engine/src/cmd/surface.rs index 6eeee282..2171b2b6 100644 --- a/crates/engine/src/cmd/surface.rs +++ b/crates/engine/src/cmd/surface.rs @@ -2,7 +2,7 @@ //! solids), and the planned surface toolset (`planned!`: registered, not built yet). use serde_json::Value; -use solidcraft_doc::{FaceAt, FeatureKind, PlaneRef, ProfileSel, RuledCorners, RuledMode, expr::Kind}; +use solidcraft_doc::{ExtendMode, ExtendTo, FaceAt, FeatureKind, PlaneRef, ProfileSel, RuledCorners, RuledMode, expr::Kind}; use super::CommandSpec; use super::features::{add_feature, axis_param, check_expr, extent_param, feature_sketch, operation, plane_param, profiles, sketch_id}; @@ -25,7 +25,7 @@ pub static COMMANDS: &[CommandSpec] = &[ CommandSpec::new("surface.extend", "Extend", extend) .at("SURFACE", "MODIFY") .icon("extend") - .params("body: planar surface name; edges: [[x,y,z]] (points on straight edges to push out); distance: expr"), + .params("body: surface name; edges: [[x,y,z]] (points on boundary edges to push out); distance: expr (unless up_to); mode?: natural (the surface itself continued: exact for planes and surfaces of revolution, a B-spline's end spans continued; the default) | tangent (a ruled run-off face along each edge, G1) | curvature (a run-off face keeping the surface's curvature across the edge, G2); up_to?: plane (XY|XZ|YZ|plane name|{origin, normal}) | {body, point} (a face, its surface extended) | \"bounds\" (natural: back to the surface's own bounds, undoing a trim); neighbours?: bool (default true: the edges next to an extended edge are lengthened to meet it; false: the extension's sides run square to the edge)"), CommandSpec::new("surface.trim", "Trim", trim) .at("SURFACE", "MODIFY") .icon("trim") @@ -329,9 +329,36 @@ fn extend(s: &mut Session, p: &Value) -> Result { if edges.is_empty() || edges.len() > 10_000 { return Err(bad(cmd, "`edges` must list points on the edges to extend")); } - let distance = expr(p, "distance").ok_or_else(|| bad(cmd, "`distance` is required"))?; + let mode = match str_(p, "mode").map(str::trim) { + None | Some("natural") => ExtendMode::Natural, + Some("tangent") => ExtendMode::Tangent, + Some("curvature") => ExtendMode::Curvature, + Some(m) => return Err(bad(cmd, format!("`mode` must be natural, tangent or curvature, not `{m}`"))), + }; + let up_to = match p.get("up_to").filter(|v| !v.is_null()) { + None => None, + Some(Value::String(b)) if b.trim() == "bounds" => Some(ExtendTo::Bounds), + Some(t @ Value::Object(o)) if o.contains_key("body") => { + let tb = t.get("body").and_then(Value::as_str).ok_or_else(|| bad(cmd, "`up_to.body` must name a body"))?.to_string(); + if s.model.state().body(&tb).is_none() { + return Err(bad(cmd, format!("no body `{tb}`"))); + } + let point = t.get("point").and_then(vec3).ok_or_else(|| bad(cmd, "`up_to.point` must be [x, y, z] on the target face"))?; + Some(ExtendTo::Face(FaceAt { body: tb, point })) + } + Some(v) => Some(ExtendTo::Plane(plane_param(s, Some(v), cmd)?)), + }; + if up_to == Some(ExtendTo::Bounds) && mode != ExtendMode::Natural { + return Err(bad(cmd, "only a natural extension goes back to the surface's bounds")); + } + let distance = match (expr(p, "distance"), &up_to) { + (Some(d), _) => d, + (None, Some(_)) => "0 mm".into(), + (None, None) => return Err(bad(cmd, "`distance` is required (or `up_to`)")), + }; check_expr(s, &distance, Kind::Length, cmd, "distance")?; - add_feature(s, p, FeatureKind::SurfaceExtend { body, edges, distance }) + let neighbours = bool_(p, "neighbours").unwrap_or(true); + add_feature(s, p, FeatureKind::SurfaceExtend { body, edges, distance, mode, up_to, neighbours }) } fn trim(s: &mut Session, p: &Value) -> Result { @@ -871,6 +898,102 @@ mod tests { assert!((area(&s, "Top") - 800.0).abs() < 1e-9, "{}", area(&s, "Top")); } + /// A quarter cylinder sheet (radius 10, heights 0…20) swept from a sketch arc. + fn quarter_cylinder(s: &mut Session, name: &str) { + run(s, "sketch.create", json!({"plane": "XY", "name": "Arc"})); + let r = 10.0 / 2f64.sqrt(); + let arc = run(s, "sketch.arc.three_point", json!({"start": [10, 0], "end": [0, 10], "through": [r, r]}))["created"]["curves"][0] + .as_str() + .unwrap() + .to_string(); + run(s, "sketch.finish", json!({})); + run(s, "surface.extrude", json!({"sketch": "Arc", "curves": [arc], "distance": 20, "body_name": name})); + } + + /// Extend a cylinder strip: naturally 5 along its axis (10 · π/2 · 25, one face, still a + /// cylinder), up to the plane z = 30 (10 · π/2 · 30), and by tangent and curvature run-offs + /// (a second face along a straight edge, 20 long: the tangent one a flat 5 × 20, the + /// curvature one the parabola s ↦ s·R + ½κs²·n with the cylinder's κ = 1/10, its section + /// ∫₀⁵ √(1 + κ²s²) ds long). Undo, redo and reopening; the default saves nothing new. + #[test] + fn extend_a_cylinder_strip_every_way() { + let pi = std::f64::consts::PI; + let mut s = Session::default(); + quarter_cylinder(&mut s, "Strip"); + let top = [10.0 * (pi / 4.0).cos(), 10.0 * (pi / 4.0).sin(), 20.0]; + run(&mut s, "surface.extend", json!({"body": "Strip", "edges": [top], "distance": "5 mm", "name": "Up"})); + let want = 10.0 * pi / 2.0 * 25.0; + assert!((area(&s, "Strip") - want).abs() < want * 1e-3, "{} vs {want}", area(&s, "Strip")); // (measured on meshes) + assert_eq!(s.world_state().body("Strip").unwrap().body.face_count(), 1); + let saved = s.doc.to_json(); + assert!(!saved.contains("\"mode\"") && !saved.contains("\"up_to\"") && !saved.contains("\"neighbours\""), "the default saves as before"); + run(&mut s, "edit.undo", json!({})); + assert!((area(&s, "Strip") - 200.0 * pi / 4.0 * 2.0).abs() < 1.0); + run(&mut s, "edit.redo", json!({})); + let cold = Session::new(solidcraft_doc::Document::from_json(&s.doc.to_json()).unwrap()); + assert!((area(&cold, "Strip") - area(&s, "Strip")).abs() < 1e-9); + run(&mut s, "edit.undo", json!({})); + + run(&mut s, "surface.extend", json!({"body": "Strip", "edges": [top], "up_to": {"origin": [0, 0, 30], "normal": [0, 0, 1]}})); + let want = 10.0 * pi / 2.0 * 30.0; + assert!((area(&s, "Strip") - want).abs() < want * 1e-3, "{} vs {want}", area(&s, "Strip")); // (measured on meshes) + run(&mut s, "edit.undo", json!({})); + + let side = [0.0, 10.0, 10.0]; + let k: f64 = 0.1; + let parabola = 2.5 * (1.0 + 25.0 * k * k).sqrt() + (5.0 * k).asinh() / (2.0 * k); + for (mode, run_off) in [("tangent", 5.0), ("curvature", parabola)] { + run(&mut s, "surface.extend", json!({"body": "Strip", "edges": [side], "distance": 5, "mode": mode})); + let st = s.world_state(); + assert_eq!(st.body("Strip").unwrap().body.face_count(), 2, "{mode}"); + let want = 10.0 * pi / 2.0 * 20.0 + run_off * 20.0; + assert!((area(&s, "Strip") - want).abs() < want * 1e-3, "{mode}: {} vs {want}", area(&s, "Strip")); // (measured on meshes) + let cold = Session::new(solidcraft_doc::Document::from_json(&s.doc.to_json()).unwrap()); + assert!((area(&cold, "Strip") - area(&s, "Strip")).abs() < 1e-9, "{mode}"); + run(&mut s, "edit.undo", json!({})); + } + } + + /// The half cylinder trimmed by a slanted plane goes back to its span (150π) when its top + /// edge is extended to the surface's bounds. + #[test] + fn extend_back_to_the_surfaces_bounds() { + let mut s = Session::default(); + run(&mut s, "sketch.create", json!({"plane": "XY", "name": "S"})); + let arc = run(&mut s, "sketch.arc.three_point", json!({"start": [10, 0], "end": [-10, 0], "through": [0, 10]}))["created"]["curves"][0] + .as_str() + .unwrap() + .to_string(); + run(&mut s, "sketch.finish", json!({})); + run(&mut s, "surface.extrude", json!({"sketch": "S", "curves": [arc], "distance": 20, "body_name": "Half"})); + run(&mut s, "surface.trim", json!({"body": "Half", "plane": {"origin": [0, 0, 15], "normal": [0, 0.5, 1]}, "keep": [0, 10, 1]})); + run(&mut s, "surface.extend", json!({"body": "Half", "edges": [[0, 10, 10]], "up_to": "bounds"})); + let want = 150.0 * std::f64::consts::PI; + assert!((area(&s, "Half") - want).abs() < want * 1e-3, "{} vs {want}", area(&s, "Half")); // (measured on meshes) + } + + /// Requests Extend can't take are refused and leave the design as it was. + #[test] + fn extend_refuses_bad_requests() { + let mut s = Session::default(); + quarter_cylinder(&mut s, "Strip"); + let before = s.doc.features.len(); + let top = [10.0 * (std::f64::consts::PI / 4.0).cos(), 10.0 * (std::f64::consts::PI / 4.0).sin(), 20.0]; + for p in [ + json!({"body": "Strip", "edges": [top]}), + json!({"body": "Strip", "edges": [top], "distance": 5, "mode": "sideways"}), + json!({"body": "Strip", "edges": [top], "up_to": "bounds", "mode": "tangent"}), + json!({"body": "Strip", "edges": [top], "up_to": {"body": "Nope", "point": [0, 0, 0]}}), + json!({"body": "Strip", "edges": [top], "up_to": {"body": "Strip"}}), + json!({"body": "Strip", "edges": [top], "up_to": "No plane"}), + json!({"body": "Strip", "edges": [], "distance": 5}), + json!({"body": "Nope", "edges": [top], "distance": 5}), + ] { + assert!(s.execute("surface.extend", &p).is_err(), "{p}"); + } + assert_eq!(s.doc.features.len(), before); + } + /// The planned surface commands are registered but refuse to run: disabled with the /// issue that tracks them, an error (never a panic) from every caller, and the design left /// as it was. diff --git a/crates/kernel/src/face_params.rs b/crates/kernel/src/face_params.rs new file mode 100644 index 00000000..d0951531 --- /dev/null +++ b/crates/kernel/src/face_params.rs @@ -0,0 +1,272 @@ +//! A face's boundary in its surface's parameters, and exact parameter lines on a surface: +//! shared by the operations that rework a face on the surface it lies on (Untrim, Extend). +//! +//! A face's loops are sampled into parameters with periodic directions unwrapped along each +//! loop; which loop is the outer one is worked out from those parameters (the loop enclosing +//! the most area), never from the order the loops are stored in: an imported face can list +//! its holes first. + +use mt::{BoundedCurve, Cut, InnerSpace, Invertible, MetricSpace, ParametricCurve, ParametricSurface, SearchNearestParameter, Transformed}; +use solidcraft_geom::Vec3; +use solidcraft_geom::tess::EvalSurface; +use truck_modeling as mt; + +use crate::body::from_p3; +use crate::tess::uv::{self, Domain}; +use crate::{KernelError, Result}; + +/// Samples along each boundary edge when working out the face's parameter loops. +pub(crate) const EDGE_SAMPLES: usize = 32; +/// Golden-section steps that place an extreme of the face's span on an edge. +const EXTREME_STEPS: usize = 60; + +/// A boundary loop of a face sampled in 3D and in the surface's parameters: for each of its +/// edges (as the loop runs) the curve parameters, points and parameters of its samples, ends +/// included. Periodic parameters are unwrapped along the loop. +pub(crate) struct LoopSamples { + pub(crate) edges: Vec<(mt::Curve, Vec<(f64, Vec3, [f64; 2])>)>, +} + +impl LoopSamples { + /// The loop as a polygon in parameters (each edge's samples but its last, which the next + /// edge starts on). + pub(crate) fn polygon(&self) -> Vec<[f64; 2]> { + self.edges.iter().flat_map(|(_, s)| s.iter().take(s.len().saturating_sub(1)).map(|x| x.2)).collect() + } + + /// Signed area of the loop in parameters (positive: counter-clockwise). + pub(crate) fn area(&self) -> f64 { + let poly = self.polygon(); + let n = poly.len(); + (0..n) + .filter_map(|i| { + let (a, b) = (poly.get(i)?, poly.get((i + 1) % n)?); + Some(0.5 * (a[0] * b[1] - b[0] * a[1])) + }) + .sum() + } + + /// How far `p` is from the loop (from its samples). + pub(crate) fn dist(&self, p: Vec3) -> f64 { + self.edges.iter().flat_map(|(_, s)| s.iter().map(|x| x.1.dist(p))).fold(f64::INFINITY, f64::min) + } +} + +/// A face's boundary loops (absolute: as they run on its surface) in the surface's parameters, +/// and which of them is the outer one. +pub(crate) struct FaceLoops { + pub(crate) eval: EvalSurface, + pub(crate) dom: Domain, + pub(crate) wires: Vec, + pub(crate) loops: Vec, + pub(crate) outer: usize, +} + +impl FaceLoops { + /// The loops of `face`, its boundary found on its surface within `accept`; `what` names the + /// operation for a surface no evaluator takes. + pub(crate) fn of(face: &mt::Face, accept: f64, what: &str) -> Result { + let eval = crate::tess::surface::eval_surface(&face.surface()) + .ok_or_else(|| KernelError::TooLarge(format!("{what} a face on this kind of surface")))?; + let dom = uv::domain(&eval); + let wires: Vec = face.absolute_boundaries().clone(); + let loops: Vec = wires.iter().map(|w| sample_loop(&eval, &dom, w, accept)).collect::>()?; + // The outer loop encloses the others: the most area in parameters. + let outer = loops + .iter() + .map(LoopSamples::area) + .enumerate() + .max_by(|a, b| a.1.abs().total_cmp(&b.1.abs())) + .map(|(i, _)| i) + .ok_or_else(|| KernelError::Invalid("the face has no boundary".into()))?; + Ok(FaceLoops { eval, dom, wires, loops, outer }) + } +} + +/// Sample a boundary loop (absolute: as it runs on the surface) in the surface's parameters. +pub(crate) fn sample_loop(s: &EvalSurface, dom: &Domain, wire: &mt::Wire, accept: f64) -> Result { + let mut edges = Vec::new(); + let mut prev: Option<[f64; 2]> = None; + for e in wire.edge_iter() { + let c = e.oriented_curve(); + let (t0, t1) = c.range_tuple(); + let mut samples = Vec::with_capacity(EDGE_SAMPLES + 1); + for k in 0..=EDGE_SAMPLES { + let t = t0 + (t1 - t0) * k as f64 / EDGE_SAMPLES as f64; + let p = from_p3(c.subs(t)); + let q = parameters(s, dom, p, prev, accept)?; + samples.push((t, p, q)); + prev = Some(q); + } + edges.push((c, samples)); + } + Ok(LoopSamples { edges }) +} + +/// The parameters of `p` on the surface, unwrapped to lie near `near` along periodic directions. +pub(crate) fn parameters(s: &EvalSurface, dom: &Domain, p: Vec3, near: Option<[f64; 2]>, accept: f64) -> Result<[f64; 2]> { + let mut q = uv::project(s, dom, p, near, accept).ok_or_else(|| KernelError::Failed("the face's boundary does not lie on its surface".into()))?; + if let Some(n) = near { + for d in 0..2 { + q[d] = uv::unwrap(q[d], n[d], dom.period[d]); + } + } + Ok(q) +} + +/// The parameter box the outer boundary is restored to: along each direction the surface's own +/// range where it has one of its own (a free-form surface's; not a periodic or unbounded +/// direction, and not an analytic surface's), else the face's own span. +pub(crate) fn natural_span(s: &EvalSurface, dom: &Domain, outer: &LoopSamples, analytic: bool, accept: f64) -> Result<[(f64, f64); 2]> { + let mut span = [(0.0, 0.0); 2]; + for (d, out) in span.iter_mut().enumerate() { + *out = match (dom.range.get(d).copied().flatten(), dom.period.get(d).copied().flatten()) { + (Some(range), None) if !analytic => range, + (_, period) => { + let lo = extreme(s, dom, outer, d, false, accept); + let hi = extreme(s, dom, outer, d, true, accept); + if let Some(p) = period + && hi - lo >= p * (1.0 - 1e-9) + { + return Err(KernelError::TooLarge("a face that goes all the way round its surface".into())); + } + (lo, hi) + } + }; + if !(out.1 - out.0 > 1e-12 * (1.0 + out.0.abs() + out.1.abs())) || !out.0.is_finite() || !out.1.is_finite() { + return Err(KernelError::Invalid("the face has no extent along its surface".into())); + } + } + Ok(span) +} + +/// The least (`max == false`) or greatest parameter along direction `d` the loop reaches: the +/// extreme sample, then placed on its edge between the neighbouring samples. Samples where the +/// surface does not move along `d` (a pole) say nothing about `d` and are passed over. +pub(crate) fn extreme(s: &EvalSurface, dom: &Domain, lp: &LoopSamples, d: usize, max: bool, accept: f64) -> f64 { + let sign = if max { 1.0 } else { -1.0 }; + let mut best: Option<(f64, usize, usize)> = None; + for (ei, (_, samples)) in lp.edges.iter().enumerate() { + for (k, (_, _, q)) in samples.iter().enumerate() { + if uv::degenerate(s, *q, d) { + continue; + } + let x = sign * q[d]; + if best.is_none_or(|(b, _, _)| x > b) { + best = Some((x, ei, k)); + } + } + } + let Some((x, ei, k)) = best else { return f64::NAN }; + let Some((curve, samples)) = lp.edges.get(ei) else { return sign * x }; + let (Some(at), Some(lo), Some(hi)) = (samples.get(k), samples.get(k.saturating_sub(1)), samples.get(k + 1).or(samples.get(k))) else { + return sign * x; + }; + // Golden-section search for the extreme between the neighbouring samples. + let value = |t: f64| -> Option { + let p = from_p3(curve.subs(t)); + let q = uv::project(s, dom, p, Some(at.2), accept)?; + Some(sign * uv::unwrap(q[d], at.2[d], dom.period[d])) + }; + let g = (5f64.sqrt() - 1.0) / 2.0; + let (mut a, mut b) = (lo.0, hi.0); + let mut best = x; + for _ in 0..EXTREME_STEPS { + let (c, e) = (b - g * (b - a), a + g * (b - a)); + let (Some(fc), Some(fe)) = (value(c), value(e)) else { break }; + best = best.max(fc).max(fe); + if fc > fe { + b = e; + } else { + a = c; + } + } + sign * best +} + +/// Whether the curve passes within `tol` of `p`. +pub(crate) fn passes(curve: &mt::Curve, p: mt::Point3, tol: f64) -> bool { + let (t0, t1) = curve.range_tuple(); + let start = (0..=64) + .map(|k| t0 + (t1 - t0) * k as f64 / 64.0) + .min_by(|a, b| curve.subs(*a).distance(p).total_cmp(&curve.subs(*b).distance(p))) + .unwrap_or(t0); + let t = curve.search_nearest_parameter(p, Some(start), 100).unwrap_or(start).clamp(t0.min(t1), t0.max(t1)); + curve.subs(t).distance(p) <= tol +} + +/// The parameter line of `surface` at `w` (v = w when `along_u`, else u = w) from parameter +/// `from` to `to` along it, exact: a line on a plane, a B-spline surface's iso-curve, and on a +/// revolved surface the profile turned to its angle or an arc about the axis. +pub(crate) fn iso_curve(surface: &mt::Surface, along_u: bool, w: f64, from: f64, to: f64) -> Result { + let at = |t: f64| if along_u { surface.subs(t, w) } else { surface.subs(w, t) }; + match surface { + mt::Surface::Plane(_) => Ok(mt::Curve::Line(mt::Line(at(from), at(to)))), + mt::Surface::BSplineSurface(_) | mt::Surface::NurbsSurface(_) => { + let hom: mt::BSplineSurface = match surface { + mt::Surface::BSplineSurface(b) => mt::BSplineSurface::new( + b.knot_vecs().clone(), + b.control_points().iter().map(|row| row.iter().map(|p| p.to_homogeneous()).collect()).collect(), + ), + mt::Surface::NurbsSurface(n) => n.non_rationalized().clone(), + _ => return Err(KernelError::Failed("not a B-spline surface".into())), + }; + let iso = crate::step_in::split::iso_bsp(&hom, w, along_u).map_err(KernelError::Failed)?; + Ok(piece(mt::Curve::NurbsCurve(iso), from, to)) + } + mt::Surface::RevolutedCurve(proc) => { + let rev = proc.entity(); + // Truck's revolved surface is (profile parameter, angle), swapped when the processor + // is turned over; the processor's matrix places it. + let along_profile = along_u == proc.orientation(); + let transform = *proc.transform(); + let (origin, axis) = (rev.origin(), rev.axis().normalize()); + if along_profile { + let profile = rev.entity_curve(); + if let mt::Curve::Line(_) = profile { + return Ok(mt::Curve::Line(mt::Line(at(from), at(to)))); + } + if let mt::Curve::IntersectionCurve(_) = profile { + return Err(KernelError::TooLarge("untrimming a surface revolved from an intersection curve".into())); + } + let turn = mt::Matrix4::from_translation(origin - mt::Point3::new(0.0, 0.0, 0.0)) + * mt::Matrix4::from_axis_angle(axis, mt::Rad(w)) + * mt::Matrix4::from_translation(mt::Point3::new(0.0, 0.0, 0.0) - origin); + Ok(piece(profile.transformed(transform * turn), from, to)) + } else { + // The profile point at `w` turns about the axis: an arc from angle `from` to `to`. + let p = rev.entity_curve().subs(w); + let centre = origin + axis * (p - origin).dot(axis); + let r = (p - centre).magnitude(); + if !(r > 1e-12) { + return Err(KernelError::Failed("an arc on the axis of a revolved surface".into())); + } + let x = (p - centre) / r; + let frame = crate::step_in::geom::Frame { o: centre, x, y: axis.cross(x), z: axis }; + let arc = if to > from { + crate::step_in::geom::conic_arc(&frame, r, r, from, to - from) + } else { + crate::step_in::geom::conic_arc(&frame, r, r, to, from - to).inverse() + }; + Ok(mt::Curve::NurbsCurve(arc).transformed(transform)) + } + } + } +} + +/// The piece of `curve` between its parameters `from` and `to`, running from `from`. +pub(crate) fn piece(mut curve: mt::Curve, from: f64, to: f64) -> mt::Curve { + let (lo, hi) = if from < to { (from, to) } else { (to, from) }; + let (t0, t1) = curve.range_tuple(); + let eps = (t1 - t0).abs() * 1e-12; + if lo > t0 + eps { + curve = curve.cut(lo); + } + if hi < t1 - eps { + let _ = curve.cut(hi); + } + if from > to { + curve.invert(); + } + curve +} diff --git a/crates/kernel/src/lib.rs b/crates/kernel/src/lib.rs index 7cb138fa..9a3c4e29 100644 --- a/crates/kernel/src/lib.rs +++ b/crates/kernel/src/lib.rs @@ -28,6 +28,7 @@ mod coplanar; mod curveblend; mod delete_face; mod emboss; +mod face_params; #[cfg(test)] mod face_source_tests; mod fill; @@ -92,6 +93,7 @@ mod splitface; mod step; mod step_in; mod step_out; +mod surface_extend; mod surface_kind; mod surface_offset; mod surface_ruled; @@ -137,6 +139,7 @@ pub use splitface::{SplitTool, split_body, split_faces, split_faces_with_map}; pub use step::step_export; pub use step_in::{ImportedBody, StepImport, StepNode, step_cache_alive, step_import, step_import_shared, step_orientation_errors, step_validate}; pub use step_out::{ExportBody, ExportProduct, StepHeader, step_export_bodies, step_export_products}; +pub use surface_extend::{ExtendLimit, ExtendMode, ExtendOptions, extend_surface}; pub use surface_kind::{FaceSurface, face_surface}; pub use surface_offset::offset_surface; pub use surface_ruled::{Ruled, RuledCorners, RuledMode, ruled_surface}; diff --git a/crates/kernel/src/splitface.rs b/crates/kernel/src/splitface.rs index 4e4c47a8..88038c6b 100644 --- a/crates/kernel/src/splitface.rs +++ b/crates/kernel/src/splitface.rs @@ -903,9 +903,25 @@ pub(crate) fn analytic_field(s: &mt::Surface) -> Option { /// Which side of a tool a point is on. pub(crate) fn side_of_point(tool: &SplitTool, p: Vec3) -> Result { - let f = fields(tool)?; - let [g] = f.as_slice() else { return Err(KernelError::Invalid("one tool at a time".into())) }; - Ok(g.eval(p)) + Ok(ToolSide::of(tool)?.at(p)) +} + +/// Which side of a tool points are on, the tool's field built once for many points. +pub(crate) struct ToolSide(Field); + +impl ToolSide { + pub(crate) fn of(tool: &SplitTool) -> Result { + let mut f = fields(tool)?; + match (f.pop(), f.is_empty()) { + (Some(g), true) => Ok(ToolSide(g)), + _ => Err(KernelError::Invalid("one tool at a time".into())), + } + } + + /// The side of `p`: its sign says which, zero on the tool (NaN where it can't be told). + pub(crate) fn at(&self, p: Vec3) -> f64 { + self.0.eval(p) + } } /// Which side of a tool a face lies on (its field well inside it). diff --git a/crates/kernel/src/surface_extend.rs b/crates/kernel/src/surface_extend.rs new file mode 100644 index 00000000..20d12c22 --- /dev/null +++ b/crates/kernel/src/surface_extend.rs @@ -0,0 +1,1693 @@ +//! Extend (Surface › Extend): boundary edges of a surface body pushed outward across their +//! face, by a distance, up to a plane or a face, or (natural only) back to the surface's own +//! bounds. Three ways to continue past an edge: +//! +//! - natural: the face's own surface continued and the face bounded again on it. Exact for a +//! plane and a surface of revolution (cylinder, cone, sphere, torus); a B-spline surface's end +//! spans continue as their own polynomials. Each point of the edge moves along the parameter +//! line that leaves the edge square to it (in space), by the arc length asked for. +//! - tangent: a ruled run-off face along the edge, in the surface's tangent plane at each point +//! (tangent-continuous, G1). +//! - curvature: a run-off face that also keeps the surface's normal curvature across the edge +//! (curvature-continuous, G2): each section `C + s·R + ½·κ·s²·n`. +//! +//! The corners: where two extended edges meet, both run out along the bisector of their +//! outward directions (a mitre); next to an edge that is not extended, the extension's side +//! continues that neighbour (it is lengthened to meet the extension), or with `neighbours` off +//! runs straight out square to the edge, as it does where the edges meet smoothly. Along the +//! edge the run-out direction turns evenly between its ends, keeping the distance measured +//! square to the edge. Faces of one body extended at a corner they share use one side edge +//! there when their sides coincide, so the body's boundary stays closed. +//! +//! Refused rather than built: an extension that folds over (a concave edge pushed past its +//! centre of curvature), crosses itself, wraps all the way round its surface or runs into a +//! pole of it; edges between two faces, and edges of holes. + +use std::collections::HashMap; + +use mt::{BoundedCurve, ParametricCurve, ParametricSurface}; +use solidcraft_geom::Vec3; +use solidcraft_geom::tess::EvalSurface; +use truck_modeling as mt; + +use crate::body::{Body, from_p3, p3}; +use crate::face_params::{FaceLoops, iso_curve, natural_span, parameters}; +use crate::record::{Rebuild, Was, rebuilt_body}; +use crate::splitface::ToolSide; +use crate::surfaces::{edge_near, surface_body}; +use crate::{KernelError, Result, SplitTool, guard}; + +/// Samples along an edge a fitted curve or run-off starts with, and the most it is refined to. +const START_SAMPLES: usize = 17; +const MAX_SAMPLES: usize = 257; +/// Most edge picks. +const MAX_PICKS: usize = 10_000; +/// The run-out direction may lean at most this far from square to its edge (its cosine; about +/// 84°): further, the extension would run nearly along the edge. +const MIN_COS: f64 = 0.1; +/// Edges whose directions agree this closely where they meet (about 2.6°) meet smoothly. +const SMOOTH_COS: f64 = 0.999; +/// Steps a march along the surface takes over the distance it is asked for, and the most it +/// takes. +const MARCH_STEPS: usize = 64; +const MAX_MARCH_STEPS: usize = 8192; +/// Bisection steps that place where a march meets its target. +const BISECT_STEPS: usize = 60; + +/// How a surface continues past its edge. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub enum ExtendMode { + /// The surface itself, continued. + #[default] + Natural, + /// A ruled run-off in the surface's tangent plane (G1). + Tangent, + /// A run-off following the surface's normal curvature across the edge (G2). + Curvature, +} + +/// How far an extension goes. +#[derive(Clone, Copy, Debug)] +pub enum ExtendLimit<'a> { + /// This far from the edge, measured across the surface square to it. + Distance(f64), + /// To where it meets a plane, or the surface of a face (extended). + UpTo(&'a SplitTool<'a>), + /// Back to the surface's own bounds (natural only), undoing a trim: a free-form surface's + /// parameter box, an analytic surface's the box of the face's own span. + Bounds, +} + +/// What [`extend_surface`] does. +#[derive(Clone, Copy, Debug)] +pub struct ExtendOptions<'a> { + pub mode: ExtendMode, + pub limit: ExtendLimit<'a>, + /// Lengthen the edges next to an extended edge to meet the extension (else its sides run + /// square to the edge). + pub neighbours: bool, +} + +/// The boundary edges of `body` picked by points on them, extended as `opts` says. A natural +/// extension rebuilds each face on its continued surface; a run-off adds a face along each +/// edge. The result is a surface body whose history is over `[body]`: an extended face is the +/// face it was (trimmed differently) and its moved edges the edges they were moved from; a +/// run-off face is made from its edge. +pub fn extend_surface(body: &Body, at: &[Vec3], opts: &ExtendOptions) -> Result { + body.require_brep("extend")?; + if at.is_empty() { + return Err(KernelError::Invalid("pick the edges to extend".into())); + } + if at.len() > MAX_PICKS { + return Err(KernelError::Invalid("too many edges picked".into())); + } + if at.iter().any(|p| !p.is_finite()) { + return Err(KernelError::Invalid("an edge point is not finite".into())); + } + match opts.limit { + ExtendLimit::Distance(d) if !(d.is_finite() && d > 0.0 && d <= 1e6) => { + return Err(KernelError::Invalid("the extension distance must be positive".into())); + } + ExtendLimit::Bounds if opts.mode != ExtendMode::Natural => { + return Err(KernelError::Invalid("only a natural extension goes back to the surface's bounds".into())); + } + _ => {} + } + // A plane's straight edges pushed out with their neighbours lengthened: the exact planar + // extension (what designs saved before the other modes build). + if let (ExtendMode::Natural, ExtendLimit::Distance(d), true) = (opts.mode, opts.limit, opts.neighbours) + && plain_planar(body) + { + return crate::surfaces::extend(body, at, d); + } + let side = match opts.limit { + ExtendLimit::UpTo(tool) => Some(ToolSide::of(tool)?), + _ => None, + }; + guard("extend", || Extension::new(body, opts, side).run(at)) +} + +/// One planar face, one loop of straight edges: what [`crate::surfaces::extend`] takes. +fn plain_planar(body: &Body) -> bool { + let faces: Vec<&mt::Face> = body.solid.face_iter().collect(); + let [face] = faces.as_slice() else { return false }; + let bounds = face.boundaries(); + let [outer] = bounds.as_slice() else { return false }; + matches!(crate::splitface::analytic_field(&face.oriented_surface()), Some(crate::splitface::Analytic::Plane { .. })) + && outer.edge_iter().all(|e| matches!(e.curve(), mt::Curve::Line(_))) +} + +/// A point of an edge with the surface's frame there. +#[derive(Clone, Copy, Debug)] +struct Local { + /// On the edge. + p: Vec3, + /// Its parameters on the surface (unwrapped along the loop). + q: [f64; 2], + /// The surface's unit normal (`S_u × S_v`, as the parameters run). + n: Vec3, + /// The edge's unit direction as the loop runs, in the tangent plane. + d: Vec3, + /// Unit, in the tangent plane, square to the edge, away from the face. + out: Vec3, + su: Vec3, + sv: Vec3, +} + +/// Where the extension of a point of an edge ends. +#[derive(Clone, Copy, Debug)] +struct Far { + /// The edge's point and frame. + at: Local, + /// The run-out direction (unit, in the tangent plane) and the run-off's length along it. + r: Vec3, + l: f64, + /// The surface's normal curvature along `r` (curvature mode). + kappa: f64, + /// The far point, and (natural) its parameters. + f: Vec3, + fq: [f64; 2], +} + +impl Far { + /// The run-off's section at `s` along it (tangent and curvature modes). + fn section(&self, s: f64) -> Vec3 { + self.at.p + self.r * s + self.at.n * (0.5 * self.kappa * s * s) + } +} + +/// A corner of the outer loop (where edge `i − 1` ends and edge `i` starts) next to an +/// extended edge: the direction the extension runs out along there, and the loop index of the +/// neighbouring edge it lengthens (when it does). +#[derive(Clone, Copy, Debug)] +struct Corner { + r: Vec3, + lengthens: Option, +} + +/// An extended edge: its loop index, curve range and the run-out angles (from square to the +/// edge, toward its direction) at its ends. +#[derive(Clone, Copy, Debug)] +struct Run { + i: usize, + t0: f64, + t1: f64, + phi0: f64, + phi1: f64, +} + +/// A face being extended: its loops in parameters and its outer loop's edges. +struct Frame { + face: mt::Face, + surface: mt::Surface, + loops: FaceLoops, + /// +1 when the outer loop runs counter-clockwise in the parameters, else −1. + ccw: f64, + edges: Vec, + chosen: Vec, +} + +impl Frame { + fn new(face: &mt::Face, ids: &[mt::EdgeID], accept: f64) -> Result { + let loops = FaceLoops::of(face, accept, "extending")?; + let outer = loops.wires.get(loops.outer).ok_or_else(|| KernelError::Invalid("the face has no boundary".into()))?; + let edges: Vec = outer.edge_iter().cloned().collect(); + let chosen: Vec = edges.iter().map(|e| ids.contains(&e.id())).collect(); + if ids.iter().any(|id| !edges.iter().any(|e| e.id() == *id)) { + return Err(KernelError::Invalid("not supported yet: extending an edge of a hole (Untrim takes holes away)".into())); + } + let area = loops.loops.get(loops.outer).map_or(0.0, |l| l.area()); + if !(area.abs() > 0.0) { + return Err(KernelError::Invalid("the face has no area on its surface".into())); + } + Ok(Frame { face: face.clone(), surface: face.surface(), ccw: area.signum(), loops, edges, chosen }) + } + + fn m(&self) -> usize { + self.edges.len() + } + + /// The outer loop's edge `i` (as the loop runs) and its samples. + fn samples(&self, i: usize) -> Result<&(mt::Curve, Vec<(f64, Vec3, [f64; 2])>)> { + self.loops.loops.get(self.loops.outer).and_then(|l| l.edges.get(i)).ok_or_else(|| KernelError::Failed("extend: a loop edge".into())) + } + + /// The curve range of edge `i` (as the loop runs). + fn range(&self, i: usize) -> Result<(f64, f64)> { + Ok(self.samples(i)?.0.range_tuple()) + } + + /// The point of edge `i` at curve parameter `t`, with the surface's frame there. + fn local(&self, i: usize, t: f64, accept: f64) -> Result { + let (curve, samples) = self.samples(i)?; + let p = from_p3(curve.subs(t)); + let hint = samples.iter().min_by(|a, b| (a.0 - t).abs().total_cmp(&(b.0 - t).abs())).map(|s| s.2); + let q = parameters(&self.loops.eval, &self.loops.dom, p, hint, accept)?; + let (_, su, sv) = self.loops.eval.derivs(q[0], q[1]).ok_or_else(|| KernelError::Failed("extend: the surface at an edge".into()))?; + let n = normal(su, sv)?; + let dc = curve.der(t); + let dc = Vec3::new(dc.x, dc.y, dc.z); + let d = (dc - n * dc.dot(n)).normalized().ok_or_else(|| KernelError::Invalid("an edge to extend has no direction".into()))?; + Ok(Local { p, q, n, d, out: d.cross(n) * self.ccw, su, sv }) + } + + /// The point of the surface at `q`. + fn at(&self, q: [f64; 2]) -> Result { + self.loops.eval.derivs(q[0], q[1]).map(|x| x.0).ok_or_else(|| KernelError::Failed("extend: the surface can't be evaluated there".into())) + } + + fn is_plane(&self) -> bool { + matches!(self.loops.eval, EvalSurface::Plane { .. }) + } +} + +/// The unit normal from the parameter derivatives; refused at a pole (where it has none). +fn normal(su: Vec3, sv: Vec3) -> Result { + let c = su.cross(sv); + if !(c.len() > 1e-12 * (su.len2() + sv.len2())) { + return Err(KernelError::Invalid("the extension runs into a pole of the surface (a point with no normal)".into())); + } + c.normalized().ok_or_else(|| KernelError::Invalid("the extension runs into a pole of the surface".into())) +} + +/// The parameter step `(a, b)` with `a·S_u + b·S_v = w` for `w` in the tangent plane, its +/// near-zero parts snapped to zero (a run along a parameter line stays on it). +fn tangent_coords(su: Vec3, sv: Vec3, w: Vec3) -> Result<[f64; 2]> { + let (e, f, g) = (su.dot(su), su.dot(sv), sv.dot(sv)); + let det = e * g - f * f; + if !(det.abs() > 1e-24 * (e * g).max(1e-300)) { + return Err(KernelError::Invalid("the extension runs into a pole of the surface".into())); + } + let (x, y) = (w.dot(su), w.dot(sv)); + let mut ab = [(g * x - f * y) / det, (e * y - f * x) / det]; + let big = ab[0].abs().max(ab[1].abs()); + for c in &mut ab { + if c.abs() <= 1e-12 * big { + *c = 0.0; + } + } + Ok(ab) +} + +/// An extension in progress: the body, what to do, and what faces extended at shared corners +/// share. +struct Extension<'a> { + body: &'a Body, + opts: &'a ExtendOptions<'a>, + side: Option, + size: f64, + /// How far a fitted curve or run-off may stray. + tol: f64, + /// How far a boundary point may be from its surface. + accept: f64, + /// Faces using each edge. + users: HashMap, + /// The far corners made at each original vertex, and the side edges run out from it. + corners: HashMap>, + sides: HashMap>, + rec: Rebuild, +} + +impl<'a> Extension<'a> { + fn new(body: &'a Body, opts: &'a ExtendOptions<'a>, side: Option) -> Extension<'a> { + let size = body.size().max(1e-6); + let mut users: HashMap = HashMap::new(); + for f in body.solid.face_iter() { + for w in f.boundaries().iter() { + for e in w.edge_iter() { + *users.entry(e.id()).or_default() += 1; + } + } + } + Extension { + body, + opts, + side, + size, + tol: (size * 1e-6).max(1e-7), + accept: (size * 1e-5).max(1e-6), + users, + corners: HashMap::new(), + sides: HashMap::new(), + rec: Rebuild::default(), + } + } + + fn run(mut self, at: &[Vec3]) -> Result { + let faces: Vec = self.body.solid.face_iter().cloned().collect(); + // The picked edges, by the face they bound. + let mut picks: Vec<(usize, Vec)> = Vec::new(); + for p in at { + let e = edge_near(self.body, *p).ok_or_else(|| KernelError::Invalid(format!("no edge at {p:?}")))?; + if self.users.get(&e.id()).copied().unwrap_or(0) != 1 { + return Err(KernelError::Invalid("an edge between two faces can't be extended: pick edges on the surface's boundary".into())); + } + let fi = faces + .iter() + .position(|f| f.boundaries().iter().any(|w| w.edge_iter().any(|x| x.id() == e.id()))) + .ok_or_else(|| KernelError::Failed("extend: the face of an edge".into()))?; + match picks.iter_mut().find(|(f, _)| *f == fi) { + Some((_, ids)) if !ids.contains(&e.id()) => ids.push(e.id()), + Some(_) => {} + None => picks.push((fi, vec![e.id()])), + } + } + picks.sort_by_key(|(f, _)| *f); + let mut out = faces.clone(); + let mut added = Vec::new(); + for (fi, ids) in &picks { + let face = faces.get(*fi).ok_or_else(|| KernelError::Failed("extend: a face".into()))?; + let frame = Frame::new(face, ids, self.accept)?; + match self.opts.mode { + ExtendMode::Natural => { + let extended = self.natural(&frame)?; + if let Some(slot) = out.get_mut(*fi) { + *slot = extended; + } + } + ExtendMode::Tangent | ExtendMode::Curvature => added.extend(self.run_offs(&frame)?), + } + } + out.extend(added); + Ok(rebuilt_body(surface_body(out), &[&self.body.solid], &self.rec)) + } + + /// The corners of the frame's outer loop next to an extended edge (`None` elsewhere). + fn corners_of(&self, fr: &Frame) -> Result>> { + let m = fr.m(); + let mut out = Vec::with_capacity(m); + for i in 0..m { + let pi = (i + m - 1) % m; + let (pc, nc) = (fr.chosen.get(pi).copied().unwrap_or(false), fr.chosen.get(i).copied().unwrap_or(false)); + if !pc && !nc { + out.push(None); + continue; + } + let prev = fr.local(pi, fr.range(pi)?.1, self.accept)?; + let next = fr.local(i, fr.range(i)?.0, self.accept)?; + let smooth = prev.d.dot(next.d) >= SMOOTH_COS; + let (r, lengthens) = if pc && nc { + let r = (prev.out + next.out) + .normalized() + .ok_or_else(|| KernelError::Invalid("two extended edges turn back on each other at a corner".into()))?; + (r, None) + } else if pc { + if !self.opts.neighbours || smooth { (prev.out, None) } else { (-next.d, Some(i)) } + } else if !self.opts.neighbours || smooth { + (next.out, None) + } else { + (prev.d, Some(pi)) + }; + for (chosen, l) in [(pc, &prev), (nc, &next)] { + if chosen && r.dot(l.out) < MIN_COS { + return Err(KernelError::Invalid(if lengthens.is_some() { + "an edge next to the extended edge meets it at too sharp an angle to be lengthened (extend without lengthening the neighbouring edges)".into() + } else { + "two extended edges meet at too sharp a corner".into() + })); + } + } + out.push(Some(Corner { r, lengthens })); + } + Ok(out) + } + + /// The extended edge `i` of the frame with its run-out angles from the corners at its ends. + fn run_of(&self, fr: &Frame, corners: &[Option], i: usize) -> Result { + let m = fr.m(); + let (t0, t1) = fr.range(i)?; + let (Some(Some(c0)), Some(Some(c1))) = (corners.get(i), corners.get((i + 1) % m)) else { + return Err(KernelError::Failed("extend: the corners of an edge".into())); + }; + let (a, b) = (fr.local(i, t0, self.accept)?, fr.local(i, t1, self.accept)?); + Ok(Run { i, t0, t1, phi0: c0.r.dot(a.d).atan2(c0.r.dot(a.out)), phi1: c1.r.dot(b.d).atan2(c1.r.dot(b.out)) }) + } + + /// Where the extension of edge `run` at curve parameter `t` ends. + fn far(&self, fr: &Frame, run: &Run, t: f64) -> Result { + let at = fr.local(run.i, t, self.accept)?; + let w = if run.t1 != run.t0 { (t - run.t0) / (run.t1 - run.t0) } else { 0.0 }; + let phi = run.phi0 + (run.phi1 - run.phi0) * w; + let r = at.out * phi.cos() + at.d * phi.sin(); + let reach = |cos: f64| -> f64 { + match self.opts.limit { + ExtendLimit::Distance(d) => d / cos.max(MIN_COS), + _ => 0.0, + } + }; + let ab = tangent_coords(at.su, at.sv, r)?; + match self.opts.mode { + ExtendMode::Tangent | ExtendMode::Curvature => { + let kappa = if self.opts.mode == ExtendMode::Curvature { normal_curvature(&fr.surface, at.q, ab, at.n) } else { 0.0 }; + let mut far = Far { at, r, l: 0.0, kappa, f: at.p, fq: at.q }; + far.l = match self.opts.limit { + ExtendLimit::Distance(_) => reach(phi.cos()), + ExtendLimit::UpTo(_) => self.meet(|s| Ok(far.section(s)), at.p)?, + ExtendLimit::Bounds => return Err(KernelError::Invalid("only a natural extension goes back to the surface's bounds".into())), + }; + far.f = far.section(far.l); + Ok(far) + } + ExtendMode::Natural => { + let line = |s: f64| [at.q[0] + ab[0] * s, at.q[1] + ab[1] * s]; + let s = match self.opts.limit { + ExtendLimit::Distance(_) => arc_param(fr, at.q, ab, reach(phi.cos()))?, + ExtendLimit::UpTo(_) => self.meet(|s| fr.at(line(s)), at.p)?, + ExtendLimit::Bounds => self.to_bounds(fr, at.q, ab)?, + }; + let fq = line(s); + let (f, su, sv) = + fr.loops.eval.derivs(fq[0], fq[1]).ok_or_else(|| KernelError::Failed("extend: the surface past its edge".into()))?; + normal(su, sv)?; + Ok(Far { at, r, l: s, kappa: 0.0, f, fq }) + } + } + } + + /// Where a run `point(s)` (s ≥ 0, from `start`) first crosses the target: stepped out until + /// it changes sides, then bisected; zero from a point on the target. Refused when the run + /// doesn't reach it within a generous length. + fn meet(&self, point: impl Fn(f64) -> Result, start: Vec3) -> Result { + let side = self.side.as_ref().ok_or_else(|| KernelError::Failed("extend: no target".into()))?; + let f0 = side.at(start); + if !f0.is_finite() { + return Err(KernelError::Invalid("the edge's side of the target can't be told".into())); + } + if f0.abs() <= self.tol { + return Ok(0.0); + } + let h = self.size / 16.0; + let reach = self.size * 50.0; + let (mut s, mut prev, mut walked) = (0.0, start, 0.0); + for _ in 0..MAX_MARCH_STEPS { + let next = point(s + h)?; + walked += next.dist(prev); + let f1 = side.at(next); + if f1.is_finite() && f1 * f0 <= 0.0 { + let (mut a, mut b) = (s, s + h); + for _ in 0..BISECT_STEPS { + let m = 0.5 * (a + b); + let fm = side.at(point(m)?); + if fm.is_finite() && fm * f0 > 0.0 { a = m } else { b = m } + } + return Ok(0.5 * (a + b)); + } + if walked > reach { + break; + } + (s, prev) = (s + h, next); + } + Err(KernelError::Invalid("the extension does not reach the target".into())) + } + + /// How far along the parameter line `q + s·ab` the surface's own bounds are (zero from a + /// point on them). + fn to_bounds(&self, fr: &Frame, q: [f64; 2], ab: [f64; 2]) -> Result { + let outer = fr.loops.loops.get(fr.loops.outer).ok_or_else(|| KernelError::Invalid("the face has no boundary".into()))?; + let analytic = crate::offset::surf_of(&fr.face, (self.size * 1e-7).max(1e-9) * 100.0).is_some(); + let span = natural_span(&fr.loops.eval, &fr.loops.dom, outer, analytic, self.accept)?; + let big = ab[0].abs().max(ab[1].abs()); + let mut exit = f64::INFINITY; + for (d, (lo, hi)) in span.iter().enumerate() { + let (Some(&x), Some(&dx)) = (q.get(d), ab.get(d)) else { continue }; + if dx.abs() <= 1e-9 * big { + continue; + } + let bound = if dx > 0.0 { *hi } else { *lo }; + exit = exit.min((bound - x) / dx); + } + if !exit.is_finite() { + return Err(KernelError::Invalid("the extension does not reach its surface's bounds".into())); + } + let scale = (span[0].1 - span[0].0).abs().max((span[1].1 - span[1].0).abs()) / big.max(1e-300); + Ok(if exit > 1e-9 * scale { exit } else { 0.0 }) + } + + /// The extension's far ends of edge `run` at `n` even curve parameters; refused when none + /// moves (the edge lies on its target already). + fn fars(&self, fr: &Frame, run: &Run, n: usize) -> Result> { + let fars: Vec = (0..n).map(|k| self.far(fr, run, run.t0 + (run.t1 - run.t0) * k as f64 / (n - 1) as f64)).collect::>()?; + if fars.iter().all(|f| !self.moved(f)) { + return Err(KernelError::Invalid(match self.opts.limit { + ExtendLimit::Bounds => "the edge already lies on its surface's bounds".into(), + _ => "the edge already lies on the target".into(), + })); + } + Ok(fars) + } + + /// Whether the extension moves the point at all. + fn moved(&self, f: &Far) -> bool { + f.f.dist(f.at.p) > self.tol * 10.0 + } + + /// The far corner vertex at `v` (an original vertex) for the point `p`: one made there + /// before (by this face or a neighbour) when it is at the same place. + fn corner_vertex(&mut self, v: &mt::Vertex, p: Vec3) -> mt::Vertex { + let near = self.tol * 10.0; + if from_p3(v.point()).dist(p) <= near { + return v.clone(); + } + let list = self.corners.entry(v.id()).or_default(); + if let Some(found) = list.iter().find(|x| from_p3(x.point()).dist(p) <= near) { + return found.clone(); + } + let made = mt::Vertex::new(p3(p)); + list.push(made.clone()); + made + } + + /// The side edge from the original vertex `v` to the far corner `f`, through `mid`: one made + /// before (by this face or a neighbour) along the same path, else `make`'s curve. + fn side_edge(&mut self, v: &mt::Vertex, f: &mt::Vertex, mid: Vec3, make: impl FnOnce() -> Result) -> Result { + let near = self.tol * 10.0; + if let Some(found) = self.sides.get(&v.id()).and_then(|list| { + list.iter().find(|e| { + let c = e.curve(); + let (t0, t1) = c.range_tuple(); + e.back().id() == f.id() && from_p3(c.subs(0.5 * (t0 + t1))).dist(mid) <= near * 10.0 + }) + }) { + return Ok(found.clone()); + } + let e = mt::Edge::try_new(v, f, make()?).map_err(|e| KernelError::Failed(format!("extend: a side edge: {e}")))?; + self.sides.entry(v.id()).or_default().push(e.clone()); + Ok(e) + } + + /// The face extended naturally: its surface continued and its outer loop bounded anew. + fn natural(&mut self, fr: &Frame) -> Result { + let m = fr.m(); + let corners = self.corners_of(fr)?; + let chosen = |i: usize| fr.chosen.get(i).copied().unwrap_or(false); + let runs: Vec> = + (0..m).map(|i| if chosen(i) { self.run_of(fr, &corners, i).map(Some) } else { Ok(None) }).collect::>()?; + // The far end at each corner next to an extended edge (from that edge's run). + let far_at: Vec> = (0..m) + .map(|i| match (runs.get(i).copied().flatten(), runs.get((i + m - 1) % m).copied().flatten()) { + (Some(r), _) => self.far(fr, &r, r.t0).map(Some), + (None, Some(r)) => self.far(fr, &r, r.t1).map(Some), + (None, None) => Ok(None), + }) + .collect::>()?; + // The new loop, piece by piece, in the old parameters; no side where a corner stays. + let moves = |i: usize| far_at.get(i).copied().flatten().is_some_and(|f| self.moved(&f)); + let mut pieces: Vec = Vec::new(); + for i in 0..m { + let ni = (i + 1) % m; + if chosen(i) { + if !chosen((i + m - 1) % m) && moves(i) { + pieces.push(Piece::Side { corner: i, out: true }); + } + pieces.push(Piece::Far(i)); + if !chosen(ni) && moves(ni) { + pieces.push(Piece::Side { corner: ni, out: false }); + } + } else { + pieces.push(Piece::Keep(i)); + } + } + // A kept edge on a parameter line, lengthened by the side next to it along that line, + // becomes one edge (unless another face shares it). + let mergeable = |i: usize| -> bool { + let Some(e) = fr.edges.get(i) else { return false }; + self.users.get(&e.id()).copied().unwrap_or(0) == 1 && straight(fr, i) + }; + let lengthens = |corner: usize, edge: usize| corners.get(corner).copied().flatten().and_then(|c| c.lengthens) == Some(edge); + let n = pieces.len(); + let mut dropped = vec![false; n]; + let mut merged: Vec> = vec![None; n]; + for (k, piece) in pieces.iter().enumerate() { + let Piece::Keep(i) = *piece else { continue }; + let (pk, nk) = ((k + n - 1) % n, (k + 1) % n); + if pk == k || !mergeable(i) { + continue; + } + let before = matches!(pieces.get(pk), Some(Piece::Side { corner, out: false }) if lengthens(*corner, i)); + let after = matches!(pieces.get(nk), Some(Piece::Side { corner, out: true }) if lengthens(*corner, i)); + if before || after { + for (at, yes) in [(pk, before), (nk, after)] { + if let Some(d) = dropped.get_mut(at) { + *d |= yes; + } + } + if let Some(slot) = merged.get_mut(k) { + *slot = Some(Piece::Merged { edge: i, from_far: before, to_far: after }); + } + } + } + let pieces: Vec = pieces + .iter() + .enumerate() + .filter(|(k, _)| !dropped.get(*k).copied().unwrap_or(false)) + .map(|(k, p)| merged.get(k).copied().flatten().unwrap_or(*p)) + .collect(); + // Far samples of each extended edge (refined later where a curve is fitted). + let samples: Vec>> = runs + .iter() + .map(|r| match r { + Some(r) => { + let fars = self.fars(fr, r, START_SAMPLES)?; + check_no_fold(&fars)?; + Ok(Some(fars)) + } + None => Ok(None), + }) + .collect::>()?; + // The new loop as a polygon in the old parameters: checked before anything is built. + let corner_q = |i: usize| -> Result<[f64; 2]> { + far_at.get(i).copied().flatten().map(|f| f.fq).ok_or_else(|| KernelError::Failed("extend: a corner".into())) + }; + let vertex_q = |i: usize| -> Result<[f64; 2]> { + fr.samples(i).and_then(|s| s.1.first().map(|x| x.2).ok_or_else(|| KernelError::Failed("extend: a vertex".into()))) + }; + let mut poly: Vec<[f64; 2]> = Vec::new(); + for p in &pieces { + match *p { + Piece::Keep(i) => { + let s = &fr.samples(i)?.1; + poly.extend(s.iter().take(s.len().saturating_sub(1)).map(|x| x.2)); + } + Piece::Side { corner, out } => poly.push(if out { vertex_q(corner)? } else { corner_q(corner)? }), + Piece::Far(i) => { + let s = samples.get(i).and_then(|s| s.as_ref()).ok_or_else(|| KernelError::Failed("extend: far samples".into()))?; + poly.extend(s.iter().take(s.len().saturating_sub(1)).map(|f| f.fq)); + } + Piece::Merged { edge, from_far, .. } => poly.push(if from_far { corner_q(edge)? } else { vertex_q(edge)? }), + } + } + check_loop(fr, &poly)?; + // The surface continued over the new loop. + let ext = continued(&fr.surface, &fr.loops.eval, &bounds_of(&poly))?; + // The edges. + let mut wire: Vec = Vec::with_capacity(pieces.len()); + for p in &pieces { + match *p { + Piece::Keep(i) => wire.push(fr.edges.get(i).cloned().ok_or_else(|| KernelError::Failed("extend: an edge".into()))?), + Piece::Side { corner, out } => { + let far = far_at.get(corner).copied().flatten().ok_or_else(|| KernelError::Failed("extend: a corner".into()))?; + let v = fr.edges.get(corner).map(|e| e.front().clone()).ok_or_else(|| KernelError::Failed("extend: a vertex".into()))?; + let fv = self.corner_vertex(&v, far.f); + let (a, b) = (far.at.q, far.fq); + let mid = fr.at([0.5 * (a[0] + b[0]), 0.5 * (a[1] + b[1])])?; + let tol = self.tol; + let e = self.side_edge(&v, &fv, mid, || uv_line(fr, &ext, a, b, tol))?; + let lengthened = corners.get(corner).copied().flatten().and_then(|c| c.lengthens).and_then(|i| fr.edges.get(i)); + if let Some(old) = lengthened { + self.rec.piece(&e, old); + } + wire.push(if out { e } else { e.inverse() }); + } + Piece::Far(i) => { + let run = runs.get(i).copied().flatten().ok_or_else(|| KernelError::Failed("extend: a run".into()))?; + let old = fr.edges.get(i).ok_or_else(|| KernelError::Failed("extend: an edge".into()))?; + let (Some(a), Some(b)) = (far_at.get(i).copied().flatten(), far_at.get((i + 1) % m).copied().flatten()) else { + return Err(KernelError::Failed("extend: the corners of an edge".into())); + }; + let va = self.corner_vertex(old.front(), a.f); + let vb = self.corner_vertex(old.back(), b.f); + let first = samples.get(i).and_then(|s| s.as_ref()).ok_or_else(|| KernelError::Failed("extend: far samples".into()))?; + let curve = self.far_curve(fr, &ext, &run, first)?; + let e = mt::Edge::try_new(&va, &vb, curve).map_err(|e| KernelError::Failed(format!("extend: a moved edge: {e}")))?; + self.rec.piece(&e, old); + wire.push(e); + } + Piece::Merged { edge, from_far, to_far } => { + let old = fr.edges.get(edge).ok_or_else(|| KernelError::Failed("extend: an edge".into()))?; + let (a, va) = if from_far { + let f = far_at.get(edge).copied().flatten().ok_or_else(|| KernelError::Failed("extend: a corner".into()))?; + (f.fq, self.corner_vertex(old.front(), f.f)) + } else { + (vertex_q(edge)?, old.front().clone()) + }; + let (b, vb) = if to_far { + let f = far_at.get((edge + 1) % m).copied().flatten().ok_or_else(|| KernelError::Failed("extend: a corner".into()))?; + (f.fq, self.corner_vertex(old.back(), f.f)) + } else { + (vertex_q((edge + 1) % m)?, old.back().clone()) + }; + let b = unwrap_near(fr, b, a); + let curve = uv_line(fr, &ext, a, b, self.tol)?; + let e = mt::Edge::try_new(&va, &vb, curve).map_err(|e| KernelError::Failed(format!("extend: a lengthened edge: {e}")))?; + self.rec.piece(&e, old); + wire.push(e); + } + } + } + let holes = fr.loops.wires.iter().enumerate().filter(|(k, _)| *k != fr.loops.outer).map(|(_, w)| w.clone()); + let boundaries: Vec = std::iter::once(mt::Wire::from(wire)).chain(holes).collect(); + let mut f = mt::Face::try_new(boundaries, ext.surface).map_err(|e| KernelError::Failed(format!("extended face: {e}")))?; + if !fr.face.orientation() { + f.invert(); + } + self.rec.face(&f, &fr.face); + Ok(f) + } + + /// The far edge of a natural extension: exact on a parameter line (and straight on a plane), + /// else a cubic through the far points, refined until it follows them within the tolerance. + fn far_curve(&self, fr: &Frame, ext: &Continued, run: &Run, first: &[Far]) -> Result { + let (Some(a), Some(b)) = (first.first(), first.last()) else { return Err(KernelError::Failed("extend: far samples".into())) }; + if let Some(along_u) = iso(first.iter().map(|f| f.fq)) { + let (w, from, to) = if along_u { (a.fq[1], a.fq[0], b.fq[0]) } else { (a.fq[0], a.fq[1], b.fq[1]) }; + return ext.iso(along_u, w, from, to); + } + if fr.is_plane() && first.iter().all(|f| f.f.dist_to_segment(a.f, b.f) <= self.tol) { + return Ok(mt::Curve::Line(mt::Line(p3(a.f), p3(b.f)))); + } + self.fit(|t| self.far(fr, run, t).map(|f| f.f), run.t0, run.t1) + } + + /// A cubic through `point(t)` at even parameters from `t0` to `t1`, refined until it passes + /// within the tolerance of the points between its samples. + fn fit(&self, point: impl Fn(f64) -> Result, t0: f64, t1: f64) -> Result { + let mut n = START_SAMPLES; + loop { + let at = |k: f64| t0 + (t1 - t0) * k / (n - 1) as f64; + let pts: Vec = (0..n).map(|k| point(at(k as f64))).collect::>()?; + let c = crate::build::interpolate_cubic(&pts).ok_or_else(|| KernelError::Failed("extend: fitting an edge".into()))?; + let mut worst = 0.0f64; + for k in 0..n - 1 { + let u = (k as f64 + 0.5) / (n - 1) as f64; + worst = worst.max(from_p3(c.subs(u)).dist(point(at(k as f64 + 0.5))?)); + } + if worst <= self.tol { + return Ok(mt::Curve::BSplineCurve(c)); + } + if n >= MAX_SAMPLES { + return Err(KernelError::Invalid(format!("the extended edge can't be fitted closely enough (off by up to {worst:.6} mm)"))); + } + n = 2 * n - 1; + } + } + + /// The run-off faces of a tangent or curvature extension, one along each extended edge. + fn run_offs(&mut self, fr: &Frame) -> Result> { + let m = fr.m(); + let corners = self.corners_of(fr)?; + let mut out = Vec::new(); + for i in 0..m { + if !fr.chosen.get(i).copied().unwrap_or(false) { + continue; + } + let run = self.run_of(fr, &corners, i)?; + let edge = fr.edges.get(i).cloned().ok_or_else(|| KernelError::Failed("extend: an edge".into()))?; + let (surface, fars, far_curve) = self.run_off_surface(fr, &run)?; + let (Some(a), Some(b)) = (fars.first().copied(), fars.last().copied()) else { + return Err(KernelError::Failed("extend: run-off samples".into())); + }; + if !self.moved(&a) || !self.moved(&b) { + return Err(KernelError::Invalid("not supported yet: a run-off that ends at a corner already on its target".into())); + } + let (fa, fb) = (self.corner_vertex(edge.front(), a.f), self.corner_vertex(edge.back(), b.f)); + let curvature = self.opts.mode == ExtendMode::Curvature; + let side = |f: Far| -> mt::Curve { + if curvature { + let pts = vec![p3(f.at.p), p3(f.at.p + f.r * (0.5 * f.l)), p3(f.f)]; + mt::Curve::BSplineCurve(mt::BSplineCurve::new(mt::KnotVec::bezier_knot(2), pts)) + } else { + mt::Curve::Line(mt::Line(p3(f.at.p), p3(f.f))) + } + }; + let start = self.side_edge(edge.front(), &fa, a.section(0.5 * a.l), || Ok(side(a)))?; + let end = self.side_edge(edge.back(), &fb, b.section(0.5 * b.l), || Ok(side(b)))?; + let far = mt::Edge::try_new(&fa, &fb, mt::Curve::BSplineCurve(far_curve)) + .map_err(|e| KernelError::Failed(format!("extend: a run-off edge: {e}")))?; + let wire: mt::Wire = vec![edge.clone(), end, far.inverse(), start.inverse()].into(); + let mut f = mt::Face::try_new(vec![wire], mt::Surface::BSplineSurface(surface.clone())) + .map_err(|e| KernelError::Failed(format!("run-off face: {e}")))?; + // Facing the way the face it runs off faces. + let mid = self.far(fr, &run, 0.5 * (run.t0 + run.t1))?; + let facing = mid.at.n * if fr.face.orientation() { 1.0 } else { -1.0 }; + let own = surface.uder(0.5, 0.0).cross(surface.vder(0.5, 0.0)); + if Vec3::new(own.x, own.y, own.z).dot(facing) < 0.0 { + f.invert(); + } + self.rec.made(&f, vec![Was::Edge(edge.id())]); + out.push(f); + } + Ok(out) + } + + /// The run-off along edge `run`: degree 1 (tangent) or 2 (curvature) across it, a cubic + /// through its sections along it, refined until it is within the tolerance of its exact + /// sections (its near side of the edge itself) between the samples; and its samples. + fn run_off_surface(&self, fr: &Frame, run: &Run) -> Result<(mt::BSplineSurface, Vec, mt::BSplineCurve)> { + let curvature = self.opts.mode == ExtendMode::Curvature; + let mut n = START_SAMPLES; + loop { + let fars = self.fars(fr, run, n)?; + check_no_fold(&fars)?; + let rows: Vec> = if curvature { + vec![ + fars.iter().map(|f| f.at.p).collect(), + fars.iter().map(|f| f.at.p + f.r * (0.5 * f.l)).collect(), + fars.iter().map(|f| f.f).collect(), + ] + } else { + vec![fars.iter().map(|f| f.at.p).collect(), fars.iter().map(|f| f.f).collect()] + }; + let curves: Vec> = rows + .iter() + .map(|r| crate::build::interpolate_cubic(r)) + .collect::>() + .ok_or_else(|| KernelError::Failed("extend: fitting a run-off".into()))?; + let ku = curves.first().map(|c| c.knot_vec().clone()).ok_or_else(|| KernelError::Failed("extend: a run-off".into()))?; + let count = curves.first().map_or(0, |c| c.control_points().len()); + let net: Vec> = (0..count).map(|k| curves.iter().filter_map(|c| c.control_points().get(k).copied()).collect()).collect(); + let surface = mt::BSplineSurface::try_new((ku, mt::KnotVec::bezier_knot(rows.len() - 1)), net) + .map_err(|e| KernelError::Failed(format!("extend: a run-off surface: {e}")))?; + let mut worst = 0.0f64; + for k in 0..n - 1 { + let u = (k as f64 + 0.5) / (n - 1) as f64; + let exact = self.far(fr, run, run.t0 + (run.t1 - run.t0) * u)?; + for sigma in [0.0, 0.5, 1.0] { + let want = exact.section(sigma * exact.l); + worst = worst.max(from_p3(surface.subs(u, sigma)).dist(want)); + } + } + if worst <= self.tol { + check_regular(&surface, n)?; + let far = curves.last().cloned().ok_or_else(|| KernelError::Failed("extend: a run-off".into()))?; + return Ok((surface, fars, far)); + } + if n >= MAX_SAMPLES { + return Err(KernelError::Invalid(format!("the run-off can't be fitted closely enough (off by up to {worst:.6} mm)"))); + } + n = 2 * n - 1; + } + } +} + +/// A piece of a naturally extended face's new outer loop. +#[derive(Clone, Copy, Debug)] +enum Piece { + /// Edge `i` of the old loop, kept. + Keep(usize), + /// The side from corner `corner`'s original vertex out to its far corner (`out`), or back. + Side { corner: usize, out: bool }, + /// Edge `i` moved out. + Far(usize), + /// Edge `edge` lengthened along its parameter line, from its start's far corner + /// (`from_far`) or vertex, to its end's far corner (`to_far`) or vertex. + Merged { edge: usize, from_far: bool, to_far: bool }, +} + +/// Whether edge `i` of the frame runs straight in the surface's parameters: along a parameter +/// line, or a straight line on a plane. +fn straight(fr: &Frame, i: usize) -> bool { + let Ok((curve, samples)) = fr.samples(i) else { return false }; + if fr.is_plane() { + return matches!(curve, mt::Curve::Line(_)); + } + iso(samples.iter().map(|s| s.2)).is_some() +} + +/// Whether the parameters all lie on one parameter line: `Some(true)` along u (v the same), +/// `Some(false)` along v. +fn iso(qs: impl Iterator + Clone) -> Option { + let span = |d: usize| -> (f64, f64) { qs.clone().fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), q| (lo.min(q[d]), hi.max(q[d]))) }; + let (su, sv) = (span(0), span(1)); + let (wu, wv) = (su.1 - su.0, sv.1 - sv.0); + let scale = wu.max(wv); + if !(scale > 0.0) { + return None; + } + let eps = 1e-9 * scale; + if wv <= eps && wu > eps { + Some(true) + } else if wu <= eps && wv > eps { + Some(false) + } else { + None + } +} + +/// `b` moved by whole periods to lie near `a`. +fn unwrap_near(fr: &Frame, b: [f64; 2], a: [f64; 2]) -> [f64; 2] { + [0, 1].map(|d| crate::tess::uv::unwrap(b[d], a[d], fr.loops.dom.period[d])) +} + +/// The curve along the straight parameter segment `a`–`b` (old parameters) on the continued +/// surface: exact along a parameter line or on a plane, else a fitted cubic. +fn uv_line(fr: &Frame, ext: &Continued, a: [f64; 2], b: [f64; 2], tol: f64) -> Result { + let dir = [b[0] - a[0], b[1] - a[1]]; + let big = dir[0].abs().max(dir[1].abs()); + if fr.is_plane() { + return Ok(mt::Curve::Line(mt::Line(p3(fr.at(a)?), p3(fr.at(b)?)))); + } + if dir[1].abs() <= 1e-9 * big { + return ext.iso(true, a[1], a[0], b[0]); + } + if dir[0].abs() <= 1e-9 * big { + return ext.iso(false, a[0], a[1], b[1]); + } + let point = |t: f64| fr.at([a[0] + dir[0] * t, a[1] + dir[1] * t]); + let mut n = START_SAMPLES; + loop { + let pts: Vec = (0..n).map(|k| point(k as f64 / (n - 1) as f64)).collect::>()?; + let c = crate::build::interpolate_cubic(&pts).ok_or_else(|| KernelError::Failed("extend: fitting a side edge".into()))?; + let mut worst = 0.0f64; + for k in 0..n - 1 { + let u = (k as f64 + 0.5) / (n - 1) as f64; + worst = worst.max(from_p3(c.subs(u)).dist(point(u)?)); + } + if worst <= tol { + return Ok(mt::Curve::BSplineCurve(c)); + } + if n >= MAX_SAMPLES { + return Err(KernelError::Invalid(format!("a side of the extension can't be fitted closely enough (off by up to {worst:.6} mm)"))); + } + n = 2 * n - 1; + } +} + +/// The parameter range `q + s·ab` covers from `q` for arc length `length` on the surface +/// (Simpson's rule over even steps, the last step bisected). +fn arc_param(fr: &Frame, q: [f64; 2], ab: [f64; 2], length: f64) -> Result { + let speed = |s: f64| -> Result { + let (_, su, sv) = fr + .loops + .eval + .derivs(q[0] + ab[0] * s, q[1] + ab[1] * s) + .ok_or_else(|| KernelError::Failed("extend: the surface past its edge".into()))?; + normal(su, sv)?; + Ok((su * ab[0] + sv * ab[1]).len()) + }; + // |S_u a + S_v b| = 1 at the start: steps of the length's share. + let h = length / MARCH_STEPS as f64; + let (mut s, mut walked, mut v0) = (0.0, 0.0, speed(0.0)?); + for _ in 0..MAX_MARCH_STEPS { + let (vm, v1) = (speed(s + 0.5 * h)?, speed(s + h)?); + let step = h / 6.0 * (v0 + 4.0 * vm + v1); + if walked + step >= length { + let part = |x: f64| -> Result { Ok(x / 6.0 * (v0 + 4.0 * speed(s + 0.5 * x)? + speed(s + x)?)) }; + let (mut a, mut b) = (0.0, h); + for _ in 0..BISECT_STEPS { + let mid = 0.5 * (a + b); + if walked + part(mid)? < length { a = mid } else { b = mid } + } + return Ok(s + 0.5 * (a + b)); + } + (s, walked, v0) = (s + h, walked + step, v1); + } + Err(KernelError::Invalid("the extension runs off its surface".into())) +} + +/// The surface's normal curvature along the tangent `a·S_u + b·S_v` (unit) at `q`, signed +/// toward `n`. +fn normal_curvature(s: &mt::Surface, q: [f64; 2], ab: [f64; 2], n: Vec3) -> f64 { + let (uu, uv, vv) = (s.uuder(q[0], q[1]), s.uvder(q[0], q[1]), s.vvder(q[0], q[1])); + let second = uu * (ab[0] * ab[0]) + uv * (2.0 * ab[0] * ab[1]) + vv * (ab[1] * ab[1]); + let k = Vec3::new(second.x, second.y, second.z).dot(n); + if k.is_finite() { k } else { 0.0 } +} + +/// Refused where the far ends run back against the edge: the extension folds over (a concave +/// edge pushed past its centre of curvature). +fn check_no_fold(fars: &[Far]) -> Result<()> { + for w in fars.windows(2) { + let (Some(a), Some(b)) = (w.first(), w.get(1)) else { continue }; + if (b.f - a.f).dot(b.at.p - a.at.p) <= 0.0 { + return Err(KernelError::Invalid( + "the extension folds over itself (a concave edge pushed past its centre of curvature): extend it less".into(), + )); + } + } + Ok(()) +} + +/// Refused where a run-off's normal turns over anywhere across it: it folds. +fn check_regular(s: &mt::BSplineSurface, n: usize) -> Result<()> { + for k in 0..n { + let u = k as f64 / (n - 1) as f64; + let base = s.uder(u, 0.0).cross(s.vder(u, 0.0)); + for j in 1..=4 { + let v = j as f64 / 4.0; + let here = s.uder(u, v).cross(s.vder(u, v)); + if !(here.x * base.x + here.y * base.y + here.z * base.z > 0.0) { + return Err(KernelError::Invalid("the extension folds over itself: extend it less".into())); + } + } + } + Ok(()) +} + +/// Refused where the new outer loop (a polygon in the old parameters) crosses itself or the +/// face's holes, or spans a whole turn of a periodic surface. +fn check_loop(fr: &Frame, poly: &[[f64; 2]]) -> Result<()> { + for d in 0..2 { + if let Some(p) = fr.loops.dom.period.get(d).copied().flatten() { + let (lo, hi) = poly.iter().fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), q| (lo.min(q[d]), hi.max(q[d]))); + if hi - lo >= p * (1.0 - 1e-9) { + return Err(KernelError::Invalid("the extension would wrap all the way round its surface onto itself".into())); + } + } + } + let segs = |pts: &[[f64; 2]]| -> Vec<([f64; 2], [f64; 2])> { + let n = pts.len(); + (0..n).filter_map(|i| Some((*pts.get(i)?, *pts.get((i + 1) % n)?))).collect() + }; + let outer = segs(poly); + let n = outer.len(); + for i in 0..n { + for j in i + 2..n { + if i == 0 && j == n - 1 { + continue; + } + if let (Some(a), Some(b)) = (outer.get(i), outer.get(j)) + && crosses(*a, *b) + { + return Err(KernelError::Invalid("the extension crosses itself".into())); + } + } + } + for (k, lp) in fr.loops.loops.iter().enumerate() { + if k == fr.loops.outer { + continue; + } + for h in segs(&lp.polygon()) { + if outer.iter().any(|o| crosses(*o, h)) { + return Err(KernelError::Invalid("the extension crosses a hole of its face".into())); + } + } + } + Ok(()) +} + +/// Whether two segments properly cross (touching ends don't). +fn crosses((a, b): ([f64; 2], [f64; 2]), (c, d): ([f64; 2], [f64; 2])) -> bool { + let orient = |p: [f64; 2], q: [f64; 2], r: [f64; 2]| (q[0] - p[0]) * (r[1] - p[1]) - (q[1] - p[1]) * (r[0] - p[0]); + let scale = [a, b, c, d].iter().map(|p| p[0].abs().max(p[1].abs())).fold(1e-300, f64::max); + let eps = 1e-12 * scale * scale; + let (o1, o2, o3, o4) = (orient(a, b, c), orient(a, b, d), orient(c, d, a), orient(c, d, b)); + ((o1 > eps && o2 < -eps) || (o1 < -eps && o2 > eps)) && ((o3 > eps && o4 < -eps) || (o3 < -eps && o4 > eps)) +} + +/// The parameter box of the points. +fn bounds_of(poly: &[[f64; 2]]) -> [(f64, f64); 2] { + let span = |d: usize| poly.iter().fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), q| (lo.min(q[d]), hi.max(q[d]))); + [span(0), span(1)] +} + +/// A face's surface continued to cover a parameter box, and how its parameters relate to the +/// old ones (`new = old · scale + offset` per direction). +struct Continued { + surface: mt::Surface, + map: [(f64, f64); 2], +} + +impl Continued { + fn to_new(&self, d: usize, x: f64) -> f64 { + let (scale, offset) = self.map.get(d).copied().unwrap_or((1.0, 0.0)); + x * scale + offset + } + + /// The exact parameter line at `w` (old parameters; v = w when `along_u`) from `from` to `to`. + fn iso(&self, along_u: bool, w: f64, from: f64, to: f64) -> Result { + let (dw, dt) = if along_u { (1, 0) } else { (0, 1) }; + iso_curve(&self.surface, along_u, self.to_new(dw, w), self.to_new(dt, from), self.to_new(dt, to)) + } +} + +/// `surface` continued to cover `need` (old parameters, with a margin): a B-spline surface's +/// end spans continued as their own polynomials, a surface of revolution's profile continued +/// (a line lengthened, a B-spline profile's end spans continued); a plane as it is. Exact on the +/// old surface. +fn continued(surface: &mt::Surface, eval: &EvalSurface, need: &[(f64, f64); 2]) -> Result { + let identity = [(1.0, 0.0), (1.0, 0.0)]; + let dom = crate::tess::uv::domain(eval); + // The range each direction needs: `None` where the surface's own covers it. + let wanted = |d: usize| -> Option<(f64, f64)> { + let (lo, hi) = need.get(d).copied()?; + let (a, b) = dom.range.get(d).copied().flatten()?; + if dom.period.get(d).copied().flatten().is_some() { + return None; + } + let margin = 0.01 * (hi - lo).max(b - a); + let eps = 1e-12 * (b - a).abs().max(1e-300); + (lo < a - eps || hi > b + eps) + .then_some((lo.min(a) - if lo < a - eps { margin } else { 0.0 }, hi.max(b) + if hi > b + eps { margin } else { 0.0 })) + }; + let unsupported = || KernelError::Invalid("not supported yet: continuing this kind of surface past its bounds".into()); + match surface { + mt::Surface::Plane(_) => Ok(Continued { surface: surface.clone(), map: identity }), + mt::Surface::BSplineSurface(_) | mt::Surface::NurbsSurface(_) => { + let (Some(EvalSurface::Nurbs { degree, knots, count, ctrl }), rational) = (Some(eval), matches!(surface, mt::Surface::NurbsSurface(_))) + else { + return Err(unsupported()); + }; + let [nu, nv] = *count; + // rows[i][j]: control point i along u, j along v. + let mut rows: Vec> = (0..nu).map(|i| (0..nv).filter_map(|j| ctrl.get(j * nu + i).copied()).collect()).collect(); + let [mut ku, mut kv] = knots.clone(); + if let Some(range) = wanted(0) { + (ku, rows) = continue_span(&ku, degree[0], &rows, range).ok_or_else(unsupported)?; + } + if let Some(range) = wanted(1) { + let cols = transpose(&rows); + let (k, cols) = continue_span(&kv, degree[1], &cols, range).ok_or_else(unsupported)?; + kv = k; + rows = transpose(&cols); + } + Ok(Continued { surface: nurbs_surface(ku, kv, rows, rational)?, map: identity }) + } + mt::Surface::RevolutedCurve(proc) => { + let pd = if proc.orientation() { 0 } else { 1 }; + let Some((lo, hi)) = wanted(pd) else { return Ok(Continued { surface: surface.clone(), map: identity }) }; + let rev = proc.entity(); + let (profile, map) = match rev.entity_curve() { + mt::Curve::Line(l) => { + let at = |t: f64| l.0 + (l.1 - l.0) * t; + let mut map = identity; + if let Some(slot) = map.get_mut(pd) { + *slot = (1.0 / (hi - lo), -lo / (hi - lo)); + } + (mt::Curve::Line(mt::Line(at(lo), at(hi))), map) + } + c @ (mt::Curve::BSplineCurve(_) | mt::Curve::NurbsCurve(_)) => { + let ec = crate::tess::surface::eval_curve(c).ok_or_else(unsupported)?; + let rows: Vec> = ec.ctrl.iter().map(|p| vec![*p]).collect(); + let (k, rows) = continue_span(&ec.knots, ec.degree, &rows, (lo, hi)).ok_or_else(unsupported)?; + (nurbs_curve(k, rows.iter().filter_map(|r| r.first().copied()).collect(), matches!(c, mt::Curve::NurbsCurve(_)))?, identity) + } + mt::Curve::IntersectionCurve(_) => return Err(unsupported()), + }; + let turned = proc.map_ref(|r| mt::RevolutedCurve::by_revolution(profile, r.origin(), r.axis())); + Ok(Continued { surface: mt::Surface::RevolutedCurve(turned), map }) + } + } +} + +fn transpose(rows: &[Vec<[f64; 4]>]) -> Vec> { + let width = rows.first().map_or(0, Vec::len); + (0..width).map(|j| rows.iter().filter_map(|r| r.get(j).copied()).collect()).collect() +} + +/// `a + (b − a)·t` on homogeneous points. +fn lerp4(a: [f64; 4], b: [f64; 4], t: f64) -> [f64; 4] { + [0, 1, 2, 3].map(|k| a[k] + (b[k] - a[k]) * t) +} + +/// A clamped B-spline (degree `p`, control `rows[i]` the points at index `i` along the +/// direction, across the other) continued to cover `(lo, hi)`: each end span's polynomial over +/// a longer span. `None` for one that is not clamped or malformed, or whose weights would not +/// stay positive. +fn continue_span(knots: &[f64], p: usize, rows: &[Vec<[f64; 4]>], (lo, hi): (f64, f64)) -> Option<(Vec, Vec>)> { + let (k, c) = continue_end(knots, p, rows, hi)?; + // The start: the same on the curve run backwards. + let k: Vec = k.iter().rev().map(|x| -x).collect(); + let c: Vec> = c.into_iter().rev().collect(); + let (k, c) = continue_end(&k, p, &c, -lo)?; + let k: Vec = k.iter().rev().map(|x| -x).collect(); + let c: Vec> = c.into_iter().rev().collect(); + c.iter().flatten().all(|x| x[3] > 0.0 && x.iter().all(|v| v.is_finite())).then_some((k, c)) +} + +/// [`continue_span`] at the end only: the last knot `a` below the end `b` repeated to the +/// degree (Boehm's knot insertion; the B-spline is unchanged), which leaves the last span's +/// own Bézier points; then de Casteljau's subdivision at `(to − a)/(b − a) > 1` gives the +/// Bézier points of the same polynomial over `a…to` (NURBS Book §5.2, A5.1; §1.4). +fn continue_end(knots: &[f64], p: usize, rows: &[Vec<[f64; 4]>], to: f64) -> Option<(Vec, Vec>)> { + let mut k = knots.to_vec(); + let mut c = rows.to_vec(); + if p == 0 || c.len() <= p || k.len() != c.len() + p + 1 || k.iter().any(|x| !x.is_finite()) { + return None; + } + let b = *k.last()?; + if !(to > b) { + return Some((k, c)); + } + if k.iter().rev().take(p + 1).any(|x| *x != b) || !(*k.get(c.len() - 1)? < b) { + return None; + } + let a = *k.get(c.len() - 1)?; + while k.iter().filter(|x| **x == a).count() < p { + let n = c.len(); + let span = n - 1; + let mut next = Vec::with_capacity(n + 1); + for i in 0..=n { + let row = if i + p <= span { + c.get(i)?.clone() + } else if i > span { + c.get(i - 1)?.clone() + } else { + let (ui, uip) = (*k.get(i)?, *k.get(i + p)?); + let alpha = if uip > ui { (a - ui) / (uip - ui) } else { 0.0 }; + let (pi, pm) = (c.get(i)?, c.get(i - 1)?); + pm.iter().zip(pi).map(|(y, x)| lerp4(*y, *x, alpha)).collect() + }; + next.push(row); + } + k.insert(span + 1, a); + c = next; + } + let n = c.len(); + let first = n - p - 1; + let tau = (to - a) / (b - a); + let width = c.first()?.len(); + for j in 0..width { + let mut pts: Vec<[f64; 4]> = (first..n).map(|i| c.get(i).and_then(|r| r.get(j)).copied()).collect::>()?; + let mut left = Vec::with_capacity(p + 1); + left.push(*pts.first()?); + for r in 1..=p { + for i in 0..=(p - r) { + let (x, y) = (*pts.get(i)?, *pts.get(i + 1)?); + *pts.get_mut(i)? = lerp4(x, y, tau); + } + left.push(*pts.first()?); + } + for (i, x) in left.into_iter().enumerate() { + *c.get_mut(first + i)?.get_mut(j)? = x; + } + } + let len = k.len(); + for x in k.iter_mut().skip(len - p - 1) { + *x = to; + } + Some((k, c)) +} + +/// A truck B-spline surface from knots and homogeneous control rows (rows[i][j], i along u): +/// rational when `rational`, else of plain points. +fn nurbs_surface(ku: Vec, kv: Vec, rows: Vec>, rational: bool) -> Result { + let kvs = (mt::KnotVec::from(ku), mt::KnotVec::from(kv)); + let bad = |e: String| KernelError::Failed(format!("extend: the continued surface: {e}")); + if rational { + let net: Vec> = rows.iter().map(|r| r.iter().map(|h| mt::Vector4::new(h[0], h[1], h[2], h[3])).collect()).collect(); + Ok(mt::Surface::NurbsSurface(mt::NurbsSurface::new(mt::BSplineSurface::try_new(kvs, net).map_err(|e| bad(e.to_string()))?))) + } else { + let net: Vec> = + rows.iter().map(|r| r.iter().map(|h| mt::Point3::new(h[0] / h[3], h[1] / h[3], h[2] / h[3])).collect()).collect(); + Ok(mt::Surface::BSplineSurface(mt::BSplineSurface::try_new(kvs, net).map_err(|e| bad(e.to_string()))?)) + } +} + +/// A truck B-spline curve from knots and homogeneous control points. +fn nurbs_curve(knots: Vec, ctrl: Vec<[f64; 4]>, rational: bool) -> Result { + let kv = mt::KnotVec::from(knots); + let bad = |e: String| KernelError::Failed(format!("extend: the continued profile: {e}")); + if rational { + let pts: Vec = ctrl.iter().map(|h| mt::Vector4::new(h[0], h[1], h[2], h[3])).collect(); + Ok(mt::Curve::NurbsCurve(mt::NurbsCurve::new(mt::BSplineCurve::try_new(kv, pts).map_err(|e| bad(e.to_string()))?))) + } else { + let pts: Vec = ctrl.iter().map(|h| mt::Point3::new(h[0] / h[3], h[1] / h[3], h[2] / h[3])).collect(); + Ok(mt::Curve::BSplineCurve(mt::BSplineCurve::try_new(kv, pts).map_err(|e| bad(e.to_string()))?)) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::{measure, patch_region}; + use mt::InnerSpace; + use solidcraft_geom::{Loop2, Plane, Region2, Seg2, Vec2}; + use std::collections::HashSet; + use std::f64::consts::{FRAC_PI_2, FRAC_PI_4, PI}; + + fn rel(a: f64, b: f64) -> f64 { + (a - b).abs() / b.abs().max(1e-9) + } + fn area(b: &Body) -> f64 { + measure(b).unwrap().area + } + fn with(mode: ExtendMode, limit: ExtendLimit, neighbours: bool) -> ExtendOptions { + ExtendOptions { mode, limit, neighbours } + } + fn natural(d: f64) -> ExtendOptions<'static> { + with(ExtendMode::Natural, ExtendLimit::Distance(d), true) + } + fn edges(b: &Body) -> usize { + b.solid.edge_iter().map(|e| e.id()).collect::>().len() + } + /// Edges bounding two faces. + fn shared(b: &Body) -> usize { + let mut count: HashMap = HashMap::new(); + for f in b.solid.face_iter() { + for w in f.boundaries().iter() { + for e in w.edge_iter() { + *count.entry(e.id()).or_default() += 1; + } + } + } + count.values().filter(|n| **n == 2).count() + } + fn square() -> Body { + let r = Region2 { + outer: Loop2::polygon(&[Vec2::new(0.0, 0.0), Vec2::new(10.0, 0.0), Vec2::new(10.0, 10.0), Vec2::new(0.0, 10.0)]), + holes: vec![], + }; + patch_region(&Plane::XY, &r).unwrap() + } + /// A strip of the cylinder of radius 10 about Z, heights 0…20, angles 0…`sweep`: one face on + /// an exact surface of revolution (a line turned about the axis), its edges two generators + /// and two arcs. + fn cylinder_strip(sweep: f64) -> Body { + let pt = mt::Point3::new; + let rev = mt::RevolutedCurve::by_revolution( + mt::Curve::Line(mt::Line(pt(10.0, 0.0, 0.0), pt(10.0, 0.0, 20.0))), + pt(0.0, 0.0, 0.0), + mt::Vector3::unit_z(), + ); + let s = mt::Surface::RevolutedCurve(mt::Processor::new(rev)); + // Counter-clockwise in (height, angle): up the first generator, along the top arc, down + // the last generator, back along the bottom arc. + let corners = [(0.0, 0.0), (1.0, 0.0), (1.0, sweep), (0.0, sweep)]; + let v: Vec = corners.iter().map(|(u, a)| mt::Vertex::new(s.subs(*u, *a))).collect(); + let edges: Vec = (0..4) + .map(|i| { + let ((u0, a0), (u1, a1)) = (corners[i], corners[(i + 1) % 4]); + let c = if a0 == a1 { iso_curve(&s, true, a0, u0, u1) } else { iso_curve(&s, false, u0, a0, a1) }; + mt::Edge::new(&v[i], &v[(i + 1) % 4], c.unwrap()) + }) + .collect(); + surface_body(vec![mt::Face::new(vec![edges.into()], s)]) + } + + /// Arc length of z = x²/20 from 0 to x. + fn parabola_arc(x: f64) -> f64 { + 0.5 * x * (1.0 + (x / 10.0).powi(2)).sqrt() + 5.0 * (x / 10.0).asinh() + } + /// The parabolic sheet z = x²/20 over 0 ≤ x, y ≤ 10: a B-spline quadratic along x (u), + /// straight along y (v), its four edges on its parameter lines. + fn parabola() -> Body { + let pt = mt::Point3::new; + let net = vec![ + vec![pt(0.0, 0.0, 0.0), pt(0.0, 10.0, 0.0)], + vec![pt(5.0, 0.0, 0.0), pt(5.0, 10.0, 0.0)], + vec![pt(10.0, 0.0, 5.0), pt(10.0, 10.0, 5.0)], + ]; + let s = mt::BSplineSurface::new((mt::KnotVec::bezier_knot(2), mt::KnotVec::bezier_knot(1)), net); + let v = [pt(0.0, 0.0, 0.0), pt(10.0, 0.0, 5.0), pt(10.0, 10.0, 5.0), pt(0.0, 10.0, 0.0)].map(mt::Vertex::new); + let bez = + |a: mt::Point3, m: mt::Point3, b: mt::Point3| mt::Curve::BSplineCurve(mt::BSplineCurve::new(mt::KnotVec::bezier_knot(2), vec![a, m, b])); + let edges = vec![ + mt::Edge::new(&v[0], &v[1], bez(pt(0.0, 0.0, 0.0), pt(5.0, 0.0, 0.0), pt(10.0, 0.0, 5.0))), + mt::Edge::new(&v[1], &v[2], mt::Curve::Line(mt::Line(v[1].point(), v[2].point()))), + mt::Edge::new(&v[2], &v[3], bez(pt(10.0, 10.0, 5.0), pt(5.0, 10.0, 0.0), pt(0.0, 10.0, 0.0))), + mt::Edge::new(&v[3], &v[0], mt::Curve::Line(mt::Line(v[3].point(), v[0].point()))), + ]; + surface_body(vec![mt::Face::new(vec![edges.into()], mt::Surface::BSplineSurface(s))]) + } + + /// The default (natural, by a distance, neighbours lengthened) on a plane's straight edges + /// is today's planar extension, exactly; the general path builds the same square. + #[test] + fn a_planar_extension_is_unchanged() { + let s = square(); + let old = crate::extend(&s, &[Vec3::new(10.0, 5.0, 0.0)], 5.0).unwrap(); + let new = extend_surface(&s, &[Vec3::new(10.0, 5.0, 0.0)], &natural(5.0)).unwrap(); + assert_eq!(measure(&old).unwrap().area, measure(&new).unwrap().area); + assert_eq!(edges(&old), edges(&new)); + // Up to the plane x = 15 takes the general path: the same 15 × 10, its lengthened sides + // one edge each. + let tool = SplitTool::Plane(Plane::from_normal(Vec3::new(15.0, 0.0, 0.0), Vec3::X).unwrap()); + let up = extend_surface(&s, &[Vec3::new(10.0, 5.0, 0.0)], &with(ExtendMode::Natural, ExtendLimit::UpTo(&tool), true)).unwrap(); + assert!(rel(area(&up), 150.0) < 1e-9, "{}", area(&up)); + assert_eq!((up.face_count(), edges(&up)), (1, 4)); + // Up to x = 17: 170. + let tool = SplitTool::Plane(Plane::from_normal(Vec3::new(17.0, 0.0, 0.0), Vec3::X).unwrap()); + let up = extend_surface(&s, &[Vec3::new(10.0, 5.0, 0.0)], &with(ExtendMode::Natural, ExtendLimit::UpTo(&tool), true)).unwrap(); + assert!(rel(area(&up), 170.0) < 1e-9, "{}", area(&up)); + let bb = measure(&up).unwrap().bbox; + assert!((bb.max.x - 17.0).abs() < 1e-9, "{bb:?}"); + } + + /// A quadrilateral (0,0), (10,0), (10,10), (0,20) (area 150) with its right edge pushed out + /// by 5. Lengthening the neighbours follows the slanted top down to (15, 5): + 37.5. Without + /// them the extension is a 5 × 10 strip: + 50. A tangent run-off is the same trapezoid as a + /// face of its own, sharing the edge. + #[test] + fn neighbours_lengthened_or_square() { + let r = Region2 { + outer: Loop2::polygon(&[Vec2::new(0.0, 0.0), Vec2::new(10.0, 0.0), Vec2::new(10.0, 10.0), Vec2::new(0.0, 20.0)]), + holes: vec![], + }; + let q = patch_region(&Plane::XY, &r).unwrap(); + let at = [Vec3::new(10.0, 5.0, 0.0)]; + let lengthened = extend_surface(&q, &at, &natural(5.0)).unwrap(); + assert!(rel(area(&lengthened), 187.5) < 1e-9, "{}", area(&lengthened)); + let square_ends = extend_surface(&q, &at, &with(ExtendMode::Natural, ExtendLimit::Distance(5.0), false)).unwrap(); + assert!(rel(area(&square_ends), 200.0) < 1e-9, "{}", area(&square_ends)); + let run_off = extend_surface(&q, &at, &with(ExtendMode::Tangent, ExtendLimit::Distance(5.0), true)).unwrap(); + assert_eq!((run_off.face_count(), shared(&run_off)), (2, 1)); + assert!(rel(area(&run_off), 187.5) < 1e-6, "{}", area(&run_off)); + let h = run_off.history().unwrap(); + assert_eq!(h.faces.iter().map(|f| f.change).collect::>(), vec![crate::history::Change::Kept, crate::history::Change::Generated]); + } + + /// A cylinder strip (radius 10, a quarter turn, heights 0…20), on an exact surface of + /// revolution: extended 5 along its axis (its top arc picked), 10 · π/2 · 25; 5 around its + /// circumference (a straight edge picked), 10 · (π/2 + ½) · 20; both, meeting at a mitre, + /// 10 · (π/2 + ½) · 25; up to the plane z = 30, 10 · π/2 · 30. Each stays one cylindrical + /// face; the sides along the extension are the old edges lengthened. + #[test] + fn a_cylinder_strip_extended_along_its_axis_and_around() { + let c = cylinder_strip(FRAC_PI_2); + assert!(rel(area(&c), 100.0 * PI) < 1e-3); + let top = Vec3::new(10.0 * FRAC_PI_4.cos(), 10.0 * FRAC_PI_4.sin(), 20.0); + let side = Vec3::new(0.0, 10.0, 10.0); + let cylinder = + |b: &Body| matches!(crate::face_surface(b, 0), Some(crate::FaceSurface::Cylinder { radius, .. }) if (radius - 10.0).abs() < 1e-6); + + let up = extend_surface(&c, &[top], &natural(5.0)).unwrap(); + let want = 10.0 * FRAC_PI_2 * 25.0; + assert!(rel(area(&up), want) < 1e-3, "{} vs {want}", area(&up)); // (measured on meshes) + assert!((measure(&up).unwrap().bbox.max.z - 25.0).abs() < 1e-6); + assert!(up.face_count() == 1 && cylinder(&up) && edges(&up) == 4, "{}", edges(&up)); + + let round = extend_surface(&c, &[side], &natural(5.0)).unwrap(); + let want = 10.0 * (FRAC_PI_2 + 0.5) * 20.0; + assert!(rel(area(&round), want) < 1e-3, "{} vs {want}", area(&round)); // (measured on meshes) + assert!(round.face_count() == 1 && cylinder(&round) && edges(&round) == 4); + // The moved edge is at angle π/2 + ½ (5 mm of arc): x = −10 sin ½. + assert!((measure(&round).unwrap().bbox.min.x + 10.0 * 0.5f64.sin()).abs() < 1e-6, "{:?}", measure(&round).unwrap().bbox); + + let both = extend_surface(&c, &[top, side], &natural(5.0)).unwrap(); + let want = 10.0 * (FRAC_PI_2 + 0.5) * 25.0; + assert!(rel(area(&both), want) < 1e-3, "{} vs {want}", area(&both)); // (measured on meshes) + + let tool = SplitTool::Plane(Plane::from_normal(Vec3::new(0.0, 0.0, 30.0), Vec3::Z).unwrap()); + let to = extend_surface(&c, &[top], &with(ExtendMode::Natural, ExtendLimit::UpTo(&tool), true)).unwrap(); + let want = 10.0 * FRAC_PI_2 * 30.0; + assert!(rel(area(&to), want) < 1e-3, "{} vs {want}", area(&to)); // (measured on meshes) + assert!((measure(&to).unwrap().bbox.max.z - 30.0).abs() < 1e-6); + + // The history: the face is the one it was, trimmed differently. + let h = up.history().unwrap(); + assert_eq!(h.faces.first().map(|f| f.change), Some(crate::history::Change::Modified)); + } + + /// The same strip swept as a B-spline sheet (a rational arc extruded): its end spans + /// continue as themselves, so the extension stays on the cylinder. + #[test] + fn an_extruded_cylinder_strip_continues_on_its_cylinder() { + let arc = vec![Seg2::Arc { center: Vec2::new(0.0, 0.0), radius: 10.0, start: 0.0, sweep: FRAC_PI_2 }]; + let sheet = crate::extrude_sheet(&Plane::XY, &[arc], 0.0, 20.0).unwrap(); + let top = Vec3::new(10.0 * FRAC_PI_4.cos(), 10.0 * FRAC_PI_4.sin(), 20.0); + let up = extend_surface(&sheet, &[top], &natural(5.0)).unwrap(); + let want = 10.0 * FRAC_PI_2 * 25.0; + assert!(rel(area(&up), want) < 1e-3, "{} vs {want}", area(&up)); // (measured on meshes) + let round = extend_surface(&sheet, &[Vec3::new(0.0, 10.0, 10.0)], &natural(5.0)).unwrap(); + let want = 10.0 * (FRAC_PI_2 + 0.5) * 20.0; + assert!(rel(area(&round), want) < 1e-3, "{} vs {want}", area(&round)); // (measured on meshes) + assert!(matches!(crate::face_surface(&round, 0), Some(crate::FaceSurface::Cylinder { radius, .. }) if (radius - 10.0).abs() < 1e-6)); + } + + /// A half cylinder pushed round by 40 mm would close on itself (π + 4 > 2π): refused. So + /// is a concave edge pushed past its centre. + #[test] + fn self_intersecting_extensions_are_refused() { + let c = cylinder_strip(PI); + let side = Vec3::new(-10.0, 0.0, 10.0); + let e = extend_surface(&c, &[side], &natural(40.0)).unwrap_err().to_string(); + assert!(e.contains("wrap"), "{e}"); + // A 20 × 10 plate with a half-round notch of radius 5 in its top: the notch's edge + // pushed out by 2 leaves a notch of radius 3 (area 200 − 4.5π, the top edges + // lengthened into it); by 6, past the notch's centre, it folds and is refused. + let outer = Loop2 { + segs: vec![ + Seg2::Line { a: Vec2::new(0.0, 0.0), b: Vec2::new(20.0, 0.0) }, + Seg2::Line { a: Vec2::new(20.0, 0.0), b: Vec2::new(20.0, 10.0) }, + Seg2::Line { a: Vec2::new(20.0, 10.0), b: Vec2::new(15.0, 10.0) }, + Seg2::Arc { center: Vec2::new(10.0, 10.0), radius: 5.0, start: 0.0, sweep: -PI }, + Seg2::Line { a: Vec2::new(5.0, 10.0), b: Vec2::new(0.0, 10.0) }, + Seg2::Line { a: Vec2::new(0.0, 10.0), b: Vec2::new(0.0, 0.0) }, + ], + }; + let plate = patch_region(&Plane::XY, &Region2 { outer, holes: vec![] }).unwrap(); + assert!(rel(area(&plate), 200.0 - 12.5 * PI) < 1e-4); // (measured on meshes) + let notch = [Vec3::new(10.0, 5.0, 0.0)]; + let less = extend_surface(&plate, ¬ch, &natural(2.0)).unwrap(); + assert!(rel(area(&less), 200.0 - 4.5 * PI) < 1e-4, "{} vs {}", area(&less), 200.0 - 4.5 * PI); // (measured on meshes) + for mode in [ExtendMode::Natural, ExtendMode::Tangent, ExtendMode::Curvature] { + let e = extend_surface(&plate, ¬ch, &with(mode, ExtendLimit::Distance(6.0), true)).unwrap_err().to_string(); + assert!(e.contains("folds") || e.contains("crosses"), "{mode:?}: {e}"); + } + } + + /// The parabolic B-spline sheet z = x²/20 (0 ≤ x, y ≤ 10) extended 5 past its edge x = 10. + /// Natural: the parabola itself continues, 5 mm of arc further (to x* where its arc length + /// grows by 5): area 10 · (s(10) + 5), reaching z = x*²/20. Curvature: a run-off face along + /// the edge, its normal and its normal curvature across the edge those of the parabola + /// there (G1, G2: κ = 0.1 / 2^1.5), area 10 · ∫₀⁵ √(1 + κ²s²) ds more. Tangent: a plane + /// run-off, 50 more. + #[test] + fn a_bspline_extended_naturally_and_by_curvature() { + let p = parabola(); + let base = 10.0 * parabola_arc(10.0); + assert!(rel(area(&p), base) < 1e-4, "{} vs {base}", area(&p)); // (measured on meshes) + let edge = [Vec3::new(10.0, 5.0, 5.0)]; + + let nat = extend_surface(&p, &edge, &natural(5.0)).unwrap(); + let (mut lo, mut hi) = (10.0, 20.0); + for _ in 0..100 { + let m = 0.5 * (lo + hi); + if parabola_arc(m) - parabola_arc(10.0) < 5.0 { lo = m } else { hi = m } + } + let x = 0.5 * (lo + hi); + let bb = measure(&nat).unwrap().bbox; + assert!((bb.max.x - x).abs() < 1e-6 && (bb.max.z - x * x / 20.0).abs() < 1e-6, "{bb:?} vs x* = {x}"); + assert!(rel(area(&nat), base + 50.0) < 1e-4, "{} vs {}", area(&nat), base + 50.0); // (measured on meshes) + assert_eq!(nat.face_count(), 1); + // On the extended surface the parabola goes on: its height at x = 12, past the old edge. + let s = nat.solid.face_iter().next().unwrap().surface(); + let (u, v) = mt::SearchNearestParameter::::search_nearest_parameter(&s, mt::Point3::new(12.0, 5.0, 7.2), None, 100).unwrap(); + let off = from_p3(s.subs(u, v)).dist(Vec3::new(12.0, 5.0, 7.2)); + assert!(off < 1e-7, "{off}"); + + let kappa = 0.1 / 2f64.powf(1.5); + let curv = extend_surface(&p, &edge, &with(ExtendMode::Curvature, ExtendLimit::Distance(5.0), true)).unwrap(); + assert_eq!((curv.face_count(), shared(&curv)), (2, 1)); + let arc = 2.5 * (1.0 + 25.0 * kappa * kappa).sqrt() + (5.0 * kappa).asinh() / (2.0 * kappa); + assert!(rel(area(&curv), base + 10.0 * arc) < 1e-4, "{} vs {}", area(&curv), base + 10.0 * arc); // (measured on meshes) + // G1 and G2 at the old edge, measured on the run-off along its length. + let run = curv.solid.face_iter().nth(1).unwrap().surface(); + let up = mt::Vector3::new(-1.0, 0.0, 1.0) / 2f64.sqrt(); + for k in 0..=8 { + let u = k as f64 / 8.0; + let (su, sv) = (run.uder(u, 0.0), run.vder(u, 0.0)); + let n = su.cross(sv).normalize(); + let n = if n.dot(up) < 0.0 { -n } else { n }; + assert!((n - up).magnitude() < 1e-6, "normal {n:?} at {u}"); + let k_across = run.vvder(u, 0.0).dot(up) / sv.magnitude2(); + assert!(rel(k_across, kappa) < 1e-6, "κ {k_across} vs {kappa} at {u}"); + } + + let flat = extend_surface(&p, &edge, &with(ExtendMode::Tangent, ExtendLimit::Distance(5.0), true)).unwrap(); + assert!(rel(area(&flat), base + 50.0) < 1e-4, "{} vs {}", area(&flat), base + 50.0); // (measured on meshes) + let run = flat.solid.face_iter().nth(1).unwrap().surface(); + let n = run.uder(0.5, 0.0).cross(run.vder(0.5, 0.0)).normalize(); + assert!(n.cross(up).magnitude() < 1e-6, "G1: {n:?}"); + } + + /// A disc of radius 5 (two half-circle edges meeting smoothly) with both edges pushed out + /// by 2: naturally, one disc of radius 7 (49π); as tangent run-offs, an annulus of two + /// faces sharing their sides where the half circles meet (25π + 24π). + #[test] + fn a_disc_grows_by_its_curved_edges() { + let disc = patch_region(&Plane::XY, &Region2 { outer: Loop2::circle(Vec2::ZERO, 5.0), holes: vec![] }).unwrap(); + let both = [Vec3::new(0.0, 5.0, 0.0), Vec3::new(0.0, -5.0, 0.0)]; + let grown = extend_surface(&disc, &both, &natural(2.0)).unwrap(); + assert_eq!(grown.face_count(), 1); + assert!(rel(area(&grown), 49.0 * PI) < 1e-4, "{} vs {}", area(&grown), 49.0 * PI); // (measured on meshes) + let ring = extend_surface(&disc, &both, &with(ExtendMode::Tangent, ExtendLimit::Distance(2.0), true)).unwrap(); + assert_eq!((ring.face_count(), shared(&ring)), (3, 4)); + assert!(rel(area(&ring), 49.0 * PI) < 1e-4, "{} vs {}", area(&ring), 49.0 * PI); // (measured on meshes) + } + + /// Two faces of a box (its top and front, sharing their edge along x) with their edges at + /// x = 10 pushed out by 5: the shared edge is lengthened in both, as one edge, so the two + /// extensions stay joined: 150 + 150, two shared edges. + #[test] + fn neighbouring_faces_extend_together() { + let b = crate::box_solid(Vec3::ZERO, Vec3::new(10.0, 10.0, 10.0)).unwrap(); + let two = crate::copy_faces(&b, &[Vec3::new(5.0, 5.0, 10.0), Vec3::new(5.0, 0.0, 5.0)]).unwrap(); + assert_eq!((two.face_count(), shared(&two)), (2, 1)); + let out = extend_surface(&two, &[Vec3::new(10.0, 5.0, 10.0), Vec3::new(10.0, 0.0, 5.0)], &natural(5.0)).unwrap(); + assert!(rel(area(&out), 300.0) < 1e-9, "{}", area(&out)); + assert_eq!((out.face_count(), shared(&out)), (2, 2)); + } + + /// Back to the surface's own bounds: the half saddle (a B-spline over 0…10 × 0…10, x from 0 + /// to 5) goes back to its whole patch; the half cylinder trimmed by a slanted plane (top + /// from z = 15 down to 10 and back) goes back to its span, heights 0…15: 150π. + #[test] + fn extended_back_to_the_surfaces_bounds() { + let pt = |x: f64, y: f64| mt::Point3::new(x, y, x * y / 10.0); + let knots = mt::KnotVec::bezier_knot(1); + let saddle = mt::BSplineSurface::new((knots.clone(), knots), vec![vec![pt(0.0, 0.0), pt(0.0, 10.0)], vec![pt(10.0, 0.0), pt(10.0, 10.0)]]); + let corners = [pt(0.0, 0.0), pt(5.0, 0.0), pt(5.0, 10.0), pt(0.0, 10.0)]; + let v: Vec = corners.iter().map(|p| mt::Vertex::new(*p)).collect(); + let e: Vec = + (0..4).map(|i| mt::Edge::new(&v[i], &v[(i + 1) % 4], mt::Curve::Line(mt::Line(corners[i], corners[(i + 1) % 4])))).collect(); + let half = surface_body(vec![mt::Face::new(vec![e.into()], mt::Surface::BSplineSurface(saddle))]); + let bounds = with(ExtendMode::Natural, ExtendLimit::Bounds, true); + let whole = extend_surface(&half, &[Vec3::new(5.0, 5.0, 2.5)], &bounds).unwrap(); + let bb = measure(&whole).unwrap().bbox; + assert!((bb.max.x - 10.0).abs() < 1e-9 && (bb.max.z - 10.0).abs() < 1e-9, "{bb:?}"); + // Already there: refused. + assert!(extend_surface(&whole, &[Vec3::new(10.0, 5.0, 5.0)], &bounds).is_err()); + + let arc = vec![Seg2::Arc { center: Vec2::new(0.0, 0.0), radius: 10.0, start: 0.0, sweep: PI }]; + let sheet = crate::extrude_sheet(&Plane::XY, &[arc], 0.0, 20.0).unwrap(); + let tool = SplitTool::Plane(Plane::from_normal(Vec3::new(0.0, 0.0, 15.0), Vec3::new(0.0, 0.5, 1.0)).unwrap()); + let cut = crate::trim(&sheet, &tool, Vec3::new(0.0, 10.0, 1.0)).unwrap(); + let back = extend_surface(&cut, &[Vec3::new(0.0, 10.0, 10.0)], &bounds).unwrap(); + let a = area(&back); + assert!(rel(a, 150.0 * PI) < 1e-3, "{a} vs {}", 150.0 * PI); // (measured on meshes) + assert!((measure(&back).unwrap().bbox.max.z - 15.0).abs() < 1e-6); + } + + /// Hostile and unsupported requests are refused, never panic. + #[test] + fn hostile_requests_are_refused() { + let s = square(); + let at = [Vec3::new(10.0, 5.0, 0.0)]; + for d in [0.0, -1.0, f64::NAN, f64::INFINITY, 1e9] { + assert!(extend_surface(&s, &at, &natural(d)).is_err(), "{d}"); + } + assert!(extend_surface(&s, &[], &natural(1.0)).is_err()); + assert!(extend_surface(&s, &[Vec3::new(f64::NAN, 0.0, 0.0)], &natural(1.0)).is_err()); + assert!(extend_surface(&s, &vec![Vec3::ZERO; MAX_PICKS + 1], &natural(1.0)).is_err()); + assert!(extend_surface(&s, &[Vec3::new(50.0, 50.0, 0.0)], &natural(1.0)).is_err()); + assert!(extend_surface(&s, &at, &with(ExtendMode::Tangent, ExtendLimit::Bounds, true)).is_err()); + // A solid's edges lie between two faces. + let b = crate::box_solid(Vec3::ZERO, Vec3::new(10.0, 10.0, 10.0)).unwrap(); + assert!(extend_surface(&b, &[Vec3::new(10.0, 5.0, 10.0)], &natural(1.0)).is_err()); + // A hole's edge. + let r = Region2 { + outer: Loop2::polygon(&[Vec2::new(0.0, 0.0), Vec2::new(20.0, 0.0), Vec2::new(20.0, 20.0), Vec2::new(0.0, 20.0)]), + holes: vec![Loop2::circle(Vec2::new(10.0, 10.0), 3.0).reversed()], + }; + let holed = patch_region(&Plane::XY, &r).unwrap(); + let e = extend_surface(&holed, &[Vec3::new(13.0, 10.0, 0.0)], &with(ExtendMode::Tangent, ExtendLimit::Distance(1.0), true)).unwrap_err(); + assert!(e.to_string().contains("hole"), "{e}"); + // A target the extension never reaches (behind the edge). + let behind = SplitTool::Plane(Plane::from_normal(Vec3::new(-5.0, 0.0, 0.0), Vec3::X).unwrap()); + assert!(extend_surface(&s, &at, &with(ExtendMode::Natural, ExtendLimit::UpTo(&behind), true)).is_err()); + } + + /// Continuing a B-spline's end span leaves it the same on its old range. + #[test] + fn a_continued_bspline_is_unchanged_on_its_range() { + // A cubic with interior knots, rational weights. + let knots = vec![0.0, 0.0, 0.0, 0.0, 0.3, 0.6, 1.0, 1.0, 1.0, 1.0]; + let ctrl: Vec> = + [[0.0, 0.0, 0.0, 1.0], [1.0, 2.0, 0.0, 0.8], [2.0, -1.0, 1.0, 1.2], [3.0, 1.0, 0.5, 1.0], [4.0, 0.0, 2.0, 0.9], [5.0, 1.0, 0.0, 1.0]] + .iter() + .map(|p| vec![[p[0] * p[3], p[1] * p[3], p[2] * p[3], p[3]]]) + .collect(); + let before = solidcraft_geom::tess::EvalCurve { degree: 3, knots: knots.clone(), ctrl: ctrl.iter().map(|r| r[0]).collect() }; + let (k, c) = continue_span(&knots, 3, &ctrl, (-0.2, 1.25)).unwrap(); + let after = solidcraft_geom::tess::EvalCurve { degree: 3, knots: k.clone(), ctrl: c.iter().map(|r| r[0]).collect() }; + assert_eq!((k.first().copied(), k.last().copied()), (Some(-0.2), Some(1.25))); + // The same curve on the old range, and past it the old end spans' own polynomials + // (which the evaluator also continues). + for i in -4..=25 { + let t = i as f64 / 20.0; + let (a, b) = (before.eval(t).unwrap().0, after.eval(t).unwrap().0); + assert!(a.dist(b) < 1e-9, "{t}: {a:?} vs {b:?}"); + } + } +} diff --git a/crates/kernel/src/surface_untrim.rs b/crates/kernel/src/surface_untrim.rs index 4e19c583..acc435bd 100644 --- a/crates/kernel/src/surface_untrim.rs +++ b/crates/kernel/src/surface_untrim.rs @@ -6,26 +6,18 @@ //! boundary runs along the surface's parameter lines, each edge exact: a line, an arc about a //! revolved surface's axis, a rotated profile, or a B-spline iso-curve. //! -//! Which loop is the outer one is worked out from the face's parameters (the loop enclosing the -//! most area), never from the order the loops are stored in: an imported face can list its holes -//! first. +//! Which loop is the outer one is worked out from the face's parameters ([`FaceLoops`]). -use mt::{BoundedCurve, Cut, InnerSpace, Invertible, MetricSpace, ParametricCurve, ParametricSurface, SearchNearestParameter, Transformed}; +use mt::{BoundedCurve, MetricSpace, ParametricCurve, ParametricSurface}; use solidcraft_geom::Vec3; -use solidcraft_geom::tess::EvalSurface; use truck_modeling as mt; -use crate::body::{Body, from_p3}; +use crate::body::Body; +use crate::face_params::{FaceLoops, iso_curve, natural_span, passes}; use crate::record::{Rebuild, rebuilt_body}; use crate::surfaces::surface_body; -use crate::tess::uv::{self, Domain}; use crate::{KernelError, Result, guard}; -/// Samples along each boundary edge when working out the face's parameter loops. -const EDGE_SAMPLES: usize = 32; -/// Golden-section steps that place an extreme of the face's span on an edge. -const EXTREME_STEPS: usize = 60; - /// The inner loops Untrim takes away. #[derive(Clone, Copy, Debug, PartialEq)] pub enum UntrimHoles<'a> { @@ -65,86 +57,11 @@ pub fn untrim(body: &Body, face: usize, holes: UntrimHoles, outer: bool, in_plac Ok(rebuilt_body(result, &[&body.solid], &rec)) } -/// A boundary loop of a face sampled in 3D and in the surface's parameters: for each of its -/// edges (as the loop runs) the curve parameters, points and parameters of its samples, ends -/// included. Periodic parameters are unwrapped along the loop. -struct LoopSamples { - edges: Vec<(mt::Curve, Vec<(f64, Vec3, [f64; 2])>)>, -} - -impl LoopSamples { - /// The loop as a polygon in parameters (each edge's samples but its last, which the next - /// edge starts on). - fn polygon(&self) -> Vec<[f64; 2]> { - self.edges.iter().flat_map(|(_, s)| s.iter().take(s.len().saturating_sub(1)).map(|x| x.2)).collect() - } - - /// Signed area of the loop in parameters (positive: counter-clockwise). - fn area(&self) -> f64 { - let poly = self.polygon(); - let n = poly.len(); - (0..n) - .filter_map(|i| { - let (a, b) = (poly.get(i)?, poly.get((i + 1) % n)?); - Some(0.5 * (a[0] * b[1] - b[0] * a[1])) - }) - .sum() - } - - /// How far `p` is from the loop (from its samples). - fn dist(&self, p: Vec3) -> f64 { - self.edges.iter().flat_map(|(_, s)| s.iter().map(|x| x.1.dist(p))).fold(f64::INFINITY, f64::min) - } -} - -/// Sample a boundary loop (absolute: as it runs on the surface) in the surface's parameters. -fn sample_loop(s: &EvalSurface, dom: &Domain, wire: &mt::Wire, accept: f64) -> Result { - let mut edges = Vec::new(); - let mut prev: Option<[f64; 2]> = None; - for e in wire.edge_iter() { - let c = e.oriented_curve(); - let (t0, t1) = c.range_tuple(); - let mut samples = Vec::with_capacity(EDGE_SAMPLES + 1); - for k in 0..=EDGE_SAMPLES { - let t = t0 + (t1 - t0) * k as f64 / EDGE_SAMPLES as f64; - let p = from_p3(c.subs(t)); - let q = parameters(s, dom, p, prev, accept)?; - samples.push((t, p, q)); - prev = Some(q); - } - edges.push((c, samples)); - } - Ok(LoopSamples { edges }) -} - -/// The parameters of `p` on the surface, unwrapped to lie near `near` along periodic directions. -fn parameters(s: &EvalSurface, dom: &Domain, p: Vec3, near: Option<[f64; 2]>, accept: f64) -> Result<[f64; 2]> { - let mut q = uv::project(s, dom, p, near, accept).ok_or_else(|| KernelError::Failed("the face's boundary does not lie on its surface".into()))?; - if let Some(n) = near { - for d in 0..2 { - q[d] = uv::unwrap(q[d], n[d], dom.period[d]); - } - } - Ok(q) -} - /// The face untrimmed (see [`untrim`]). fn untrimmed_face(face: &mt::Face, holes: UntrimHoles, outer: bool, size: f64) -> Result { let surface = face.surface(); - let eval = - crate::tess::surface::eval_surface(&surface).ok_or_else(|| KernelError::TooLarge("untrimming a face on this kind of surface".into()))?; - let dom = uv::domain(&eval); let accept = (size * 1e-5).max(1e-6); - let wires = face.absolute_boundaries(); - let loops: Vec = wires.iter().map(|w| sample_loop(&eval, &dom, w, accept)).collect::>()?; - // The outer loop encloses the others: the most area in parameters. - let areas: Vec = loops.iter().map(LoopSamples::area).collect(); - let outer_index = areas - .iter() - .enumerate() - .max_by(|a, b| a.1.abs().total_cmp(&b.1.abs())) - .map(|(i, _)| i) - .ok_or_else(|| KernelError::Invalid("the face has no boundary".into()))?; + let FaceLoops { eval, dom, wires, loops, outer: outer_index } = FaceLoops::of(face, accept, "untrimming")?; let inner: Vec = (0..wires.len()).filter(|i| *i != outer_index).collect(); let removed: Vec = match holes { UntrimHoles::All => inner.clone(), @@ -173,13 +90,13 @@ fn untrimmed_face(face: &mt::Face, holes: UntrimHoles, outer: bool, size: f64) - "nothing to untrim: pick holes to remove or restore the outer boundary".into() })); } - let (Some(outer_wire), Some(outer_loop), Some(outer_area)) = (wires.get(outer_index), loops.get(outer_index), areas.get(outer_index)) else { + let (Some(outer_wire), Some(outer_loop)) = (wires.get(outer_index), loops.get(outer_index)) else { return Err(KernelError::Invalid("the face has no boundary".into())); }; let first = if outer { let analytic = crate::offset::surf_of(face, (size * 1e-7).max(1e-9) * 100.0).is_some(); let span = natural_span(&eval, &dom, outer_loop, analytic, accept)?; - bounds_wire(&surface, span, *outer_area > 0.0, size)? + bounds_wire(&surface, span, outer_loop.area() > 0.0, size)? } else { outer_wire.clone() }; @@ -192,76 +109,6 @@ fn untrimmed_face(face: &mt::Face, holes: UntrimHoles, outer: bool, size: f64) - Ok(f) } -/// The parameter box the outer boundary is restored to: along each direction the surface's own -/// range where it has one of its own (a free-form surface's; not a periodic or unbounded -/// direction, and not an analytic surface's), else the face's own span. -fn natural_span(s: &EvalSurface, dom: &Domain, outer: &LoopSamples, analytic: bool, accept: f64) -> Result<[(f64, f64); 2]> { - let mut span = [(0.0, 0.0); 2]; - for (d, out) in span.iter_mut().enumerate() { - *out = match (dom.range.get(d).copied().flatten(), dom.period.get(d).copied().flatten()) { - (Some(range), None) if !analytic => range, - (_, period) => { - let lo = extreme(s, dom, outer, d, false, accept); - let hi = extreme(s, dom, outer, d, true, accept); - if let Some(p) = period - && hi - lo >= p * (1.0 - 1e-9) - { - return Err(KernelError::TooLarge("a face that goes all the way round its surface".into())); - } - (lo, hi) - } - }; - if !(out.1 - out.0 > 1e-12 * (1.0 + out.0.abs() + out.1.abs())) || !out.0.is_finite() || !out.1.is_finite() { - return Err(KernelError::Invalid("the face has no extent along its surface".into())); - } - } - Ok(span) -} - -/// The least (`max == false`) or greatest parameter along direction `d` the loop reaches: the -/// extreme sample, then placed on its edge between the neighbouring samples. Samples where the -/// surface does not move along `d` (a pole) say nothing about `d` and are passed over. -fn extreme(s: &EvalSurface, dom: &Domain, lp: &LoopSamples, d: usize, max: bool, accept: f64) -> f64 { - let sign = if max { 1.0 } else { -1.0 }; - let mut best: Option<(f64, usize, usize)> = None; - for (ei, (_, samples)) in lp.edges.iter().enumerate() { - for (k, (_, _, q)) in samples.iter().enumerate() { - if uv::degenerate(s, *q, d) { - continue; - } - let x = sign * q[d]; - if best.is_none_or(|(b, _, _)| x > b) { - best = Some((x, ei, k)); - } - } - } - let Some((x, ei, k)) = best else { return f64::NAN }; - let Some((curve, samples)) = lp.edges.get(ei) else { return sign * x }; - let (Some(at), Some(lo), Some(hi)) = (samples.get(k), samples.get(k.saturating_sub(1)), samples.get(k + 1).or(samples.get(k))) else { - return sign * x; - }; - // Golden-section search for the extreme between the neighbouring samples. - let value = |t: f64| -> Option { - let p = from_p3(curve.subs(t)); - let q = uv::project(s, dom, p, Some(at.2), accept)?; - Some(sign * uv::unwrap(q[d], at.2[d], dom.period[d])) - }; - let g = (5f64.sqrt() - 1.0) / 2.0; - let (mut a, mut b) = (lo.0, hi.0); - let mut best = x; - for _ in 0..EXTREME_STEPS { - let (c, e) = (b - g * (b - a), a + g * (b - a)); - let (Some(fc), Some(fe)) = (value(c), value(e)) else { break }; - best = best.max(fc).max(fe); - if fc > fe { - b = e; - } else { - a = c; - } - } - sign * best -} - /// The closed boundary of the parameter box `span` on `surface`, counter-clockwise in its /// parameters when `ccw`: four parameter-line edges, fewer where a side shrinks to a point (a /// pole: the corners at its ends are one vertex). @@ -313,93 +160,6 @@ fn bounds_wire(surface: &mt::Surface, span: [(f64, f64); 2], ccw: bool, size: f6 Ok(edges.into()) } -/// Whether the curve passes within `tol` of `p`. -fn passes(curve: &mt::Curve, p: mt::Point3, tol: f64) -> bool { - let (t0, t1) = curve.range_tuple(); - let start = (0..=64) - .map(|k| t0 + (t1 - t0) * k as f64 / 64.0) - .min_by(|a, b| curve.subs(*a).distance(p).total_cmp(&curve.subs(*b).distance(p))) - .unwrap_or(t0); - let t = curve.search_nearest_parameter(p, Some(start), 100).unwrap_or(start).clamp(t0.min(t1), t0.max(t1)); - curve.subs(t).distance(p) <= tol -} - -/// The parameter line of `surface` at `w` (v = w when `along_u`, else u = w) from parameter -/// `from` to `to` along it, exact: a line on a plane, a B-spline surface's iso-curve, and on a -/// revolved surface the profile turned to its angle or an arc about the axis. -fn iso_curve(surface: &mt::Surface, along_u: bool, w: f64, from: f64, to: f64) -> Result { - let at = |t: f64| if along_u { surface.subs(t, w) } else { surface.subs(w, t) }; - match surface { - mt::Surface::Plane(_) => Ok(mt::Curve::Line(mt::Line(at(from), at(to)))), - mt::Surface::BSplineSurface(_) | mt::Surface::NurbsSurface(_) => { - let hom: mt::BSplineSurface = match surface { - mt::Surface::BSplineSurface(b) => mt::BSplineSurface::new( - b.knot_vecs().clone(), - b.control_points().iter().map(|row| row.iter().map(|p| p.to_homogeneous()).collect()).collect(), - ), - mt::Surface::NurbsSurface(n) => n.non_rationalized().clone(), - _ => return Err(KernelError::Failed("not a B-spline surface".into())), - }; - let iso = crate::step_in::split::iso_bsp(&hom, w, along_u).map_err(KernelError::Failed)?; - Ok(piece(mt::Curve::NurbsCurve(iso), from, to)) - } - mt::Surface::RevolutedCurve(proc) => { - let rev = proc.entity(); - // Truck's revolved surface is (profile parameter, angle), swapped when the processor - // is turned over; the processor's matrix places it. - let along_profile = along_u == proc.orientation(); - let transform = *proc.transform(); - let (origin, axis) = (rev.origin(), rev.axis().normalize()); - if along_profile { - let profile = rev.entity_curve(); - if let mt::Curve::Line(_) = profile { - return Ok(mt::Curve::Line(mt::Line(at(from), at(to)))); - } - if let mt::Curve::IntersectionCurve(_) = profile { - return Err(KernelError::TooLarge("untrimming a surface revolved from an intersection curve".into())); - } - let turn = mt::Matrix4::from_translation(origin - mt::Point3::new(0.0, 0.0, 0.0)) - * mt::Matrix4::from_axis_angle(axis, mt::Rad(w)) - * mt::Matrix4::from_translation(mt::Point3::new(0.0, 0.0, 0.0) - origin); - Ok(piece(profile.transformed(transform * turn), from, to)) - } else { - // The profile point at `w` turns about the axis: an arc from angle `from` to `to`. - let p = rev.entity_curve().subs(w); - let centre = origin + axis * (p - origin).dot(axis); - let r = (p - centre).magnitude(); - if !(r > 1e-12) { - return Err(KernelError::Failed("an arc on the axis of a revolved surface".into())); - } - let x = (p - centre) / r; - let frame = crate::step_in::geom::Frame { o: centre, x, y: axis.cross(x), z: axis }; - let arc = if to > from { - crate::step_in::geom::conic_arc(&frame, r, r, from, to - from) - } else { - crate::step_in::geom::conic_arc(&frame, r, r, to, from - to).inverse() - }; - Ok(mt::Curve::NurbsCurve(arc).transformed(transform)) - } - } - } -} - -/// The piece of `curve` between its parameters `from` and `to`, running from `from`. -fn piece(mut curve: mt::Curve, from: f64, to: f64) -> mt::Curve { - let (lo, hi) = if from < to { (from, to) } else { (to, from) }; - let (t0, t1) = curve.range_tuple(); - let eps = (t1 - t0).abs() * 1e-12; - if lo > t0 + eps { - curve = curve.cut(lo); - } - if hi < t1 - eps { - let _ = curve.cut(hi); - } - if from > to { - curve.invert(); - } - curve -} - #[cfg(test)] mod tests { use super::*; diff --git a/crates/ui-egui/src/dialogs_surface.rs b/crates/ui-egui/src/dialogs_surface.rs index ce5cf23b..27d37890 100644 --- a/crates/ui-egui/src/dialogs_surface.rs +++ b/crates/ui-egui/src/dialogs_surface.rs @@ -661,8 +661,33 @@ pub fn for_feature(app: &SolveApp, kind: &FeatureKind, frame: &Mat) -> Option<(K Surf::Thicken { thickness: thickness.clone(), symmetric: *symmetric, operation: op_index(operation) }, vec![bodies.iter().map(|n| Sel::Body { name: n.clone() }).collect()], ), - FeatureKind::SurfaceExtend { edges, distance, .. } => { - (Surf::Extend { distance: distance.clone() }, vec![edges.iter().filter_map(edge).collect()]) + FeatureKind::SurfaceExtend { edges, distance, mode, up_to, neighbours, .. } => { + // The dialog edits the edges and the distance; the rest is kept as it was. + let mut extra = Map::new(); + if !mode.is_natural() { + extra.insert("mode".into(), json!(mode)); + } + if !*neighbours { + extra.insert("neighbours".into(), json!(false)); + } + match up_to { + Some(solidcraft_engine::doc::ExtendTo::Bounds) => { + extra.insert("up_to".into(), json!("bounds")); + } + Some(solidcraft_engine::doc::ExtendTo::Face(t)) => { + extra.insert("up_to".into(), json!({"body": t.body, "point": pt(t.point)})); + } + Some(solidcraft_engine::doc::ExtendTo::Plane(PlaneRef::Origin { name } | PlaneRef::Construction { name })) => { + extra.insert("up_to".into(), json!(name)); + } + _ => {} + } + let k = Surf::Extend { distance: distance.clone() }; + let mut ins = inputs(&k); + if let Some(inp) = ins.first_mut() { + inp.items = edges.iter().filter_map(edge).collect(); + } + return Some((Kind::Surface(k), ins, extra)); } FeatureKind::SurfaceTrim { body, plane, tool, keep } => ( Surf::Trim,