diff --git a/ATTRIBUTION.md b/ATTRIBUTION.md index b5b8c4d2..97dca7e5 100644 --- a/ATTRIBUTION.md +++ b/ATTRIBUTION.md @@ -36,6 +36,7 @@ Editable drawing code, native projects and shipped exports are collected under ` | assets/archive/icons/macos-v1/exports/icon.iconset/icon_256x256@2x.png | SolidCraft contributors | native sips export from assets/archive/icons/macos-v1/exports/icon.png | MIT OR Apache-2.0 | | assets/archive/icons/macos-v1/exports/icon.iconset/icon_512x512.png | SolidCraft contributors | native sips export from assets/archive/icons/macos-v1/exports/icon.png | MIT OR Apache-2.0 | | assets/archive/icons/macos-v1/exports/icon.iconset/icon_512x512@2x.png | SolidCraft contributors | native sips export from assets/archive/icons/macos-v1/exports/icon.png | MIT OR Apache-2.0 | +| docs/screenshots/class-a-boundary.png | SolidCraft contributors | screenshot of SolidCraft itself: a patch through two sketch lines built by `surface.boundary` with `crown: 0.4` (`solidcraft-cli snapshot`) | MIT OR Apache-2.0 | | docs/screenshots/sample-plate.png | SolidCraft contributors | screenshot of SolidCraft itself (`solidcraft --sample`) | MIT OR Apache-2.0 | | docs/screenshots/demo-0.8.gif | SolidCraft contributors | original screen recording of the actual SolidCraft 0.8 English desktop UI modelling an original Raspberry Pi 4 case (no logos or branding; the engraved text is set in the system's DejaVu Sans Bold), generated with examples/demo/record_ui.sh and its original animated pointer overlay (9× playback) | MIT OR Apache-2.0 | | docs/screenshots/demo-0.8.mp4 | SolidCraft contributors | original screen recording of the actual SolidCraft 0.8 English desktop UI modelling an original Raspberry Pi 4 case (no logos or branding; the engraved text is set in the system's DejaVu Sans Bold), generated with examples/demo/record_ui.sh and its original animated pointer overlay (9× playback) | MIT OR Apache-2.0 | diff --git a/assets/locales/de.json b/assets/locales/de.json index 69381095..c74bd78c 100644 --- a/assets/locales/de.json +++ b/assets/locales/de.json @@ -366,8 +366,10 @@ "cmd.surface_blend": "Verbindungsfläche", "cmd.surface_boundary": "Begrenzungsfläche", "cmd.surface_extrude": "Flächen-Extrusion", + "cmd.surface_fair": "Flächenbegradigung", "cmd.surface_fillet": "Flächenverrundung", "cmd.surface_loft": "Flächen-Loft", + "cmd.surface_match": "Flächenanpassung", "cmd.surface_offset": "Offsetfläche", "cmd.surface_patch": "Patch", "cmd.surface_reverse_normal": "Normale umkehren", diff --git a/assets/locales/en.json b/assets/locales/en.json index 5767b5dd..ffbe52a3 100644 --- a/assets/locales/en.json +++ b/assets/locales/en.json @@ -568,8 +568,10 @@ "cmd.surface_blend": "Blend Surface", "cmd.surface_boundary": "Boundary Surface", "cmd.surface_extrude": "Surface Extrude", + "cmd.surface_fair": "Fair Surface", "cmd.surface_fillet": "Surface Fillet", "cmd.surface_loft": "Surface Loft", + "cmd.surface_match": "Match Surface", "cmd.surface_offset": "Offset Surface", "cmd.surface_patch": "Patch", "cmd.surface_reverse_normal": "Reverse Normal", diff --git a/assets/locales/ja.json b/assets/locales/ja.json index 350dccc7..4bb4502f 100644 --- a/assets/locales/ja.json +++ b/assets/locales/ja.json @@ -568,8 +568,10 @@ "cmd.surface_blend": "ブレンドサーフェス", "cmd.surface_boundary": "境界サーフェス", "cmd.surface_extrude": "サーフェス押し出し", + "cmd.surface_fair": "サーフェスフェアリング", "cmd.surface_fillet": "サーフェスフィレット", "cmd.surface_loft": "サーフェスロフト", + "cmd.surface_match": "サーフェス一致", "cmd.surface_offset": "オフセットサーフェス", "cmd.surface_patch": "パッチ", "cmd.surface_reverse_normal": "法線を反転", diff --git a/assets/locales/ko.json b/assets/locales/ko.json index 09c175fb..a274a358 100644 --- a/assets/locales/ko.json +++ b/assets/locales/ko.json @@ -366,8 +366,10 @@ "cmd.surface_blend": "블렌드 서피스", "cmd.surface_boundary": "바운더리 서피스", "cmd.surface_extrude": "서피스 돌출", + "cmd.surface_fair": "서피스 페어링", "cmd.surface_fillet": "서피스 필렛", "cmd.surface_loft": "서피스 로프트", + "cmd.surface_match": "서피스 일치", "cmd.surface_offset": "오프셋 서피스", "cmd.surface_patch": "패치", "cmd.surface_reverse_normal": "법선 반전", diff --git a/assets/locales/zh.json b/assets/locales/zh.json index 264d61d7..b310f839 100644 --- a/assets/locales/zh.json +++ b/assets/locales/zh.json @@ -568,8 +568,10 @@ "cmd.surface_blend": "混合曲面", "cmd.surface_boundary": "边界曲面", "cmd.surface_extrude": "曲面拉伸", + "cmd.surface_fair": "曲面光顺", "cmd.surface_fillet": "曲面圆角", "cmd.surface_loft": "曲面放样", + "cmd.surface_match": "曲面匹配", "cmd.surface_offset": "偏移曲面", "cmd.surface_patch": "补面", "cmd.surface_reverse_normal": "反转法向", diff --git a/crates/doc/src/document.rs b/crates/doc/src/document.rs index 37062da1..21e4125c 100644 --- a/crates/doc/src/document.rs +++ b/crates/doc/src/document.rs @@ -1384,6 +1384,59 @@ pub enum FeatureKind { #[serde(default, skip_serializing_if = "std::ops::Not::not")] stitch: bool, }, + /// One surface patch through a network of two to four boundary curves of a sketch + /// (`surface.boundary`, issue #763): the patch runs through the curves in loop order and leaves + /// each boundary at its own `order` — a sketch curve's boundary is a plane, so G1 and up hold + /// exactly — with `crown` shaping the middle without touching a boundary. + SurfaceBoundary { + sketch: SketchId, + /// The curves in loop order: each one ends where the next begins, and the last where the + /// first does. + curves: Vec, + /// How much of the surface across each boundary the patch repeats, one per curve (none: a + /// tangent plane, `g1`, on every side). + #[serde(default, skip_serializing_if = "Vec::is_empty")] + orders: Vec, + /// How full the middle is: 0 the plain construction, up to 3 fuller, down to -0.9 flatter. + #[serde(default)] + crown: f64, + /// The stations every boundary carries. + #[serde(default = "boundary_samples")] + samples: usize, + }, + /// One face of a surface body rebuilt against another body (`surface.match`): the rebuilt + /// boundary nearest the chain of `edges` on `target` takes that chain's place and repeats the + /// target's own shape along the seam up to `order`, while the boundary opposite it, and the + /// body's name, stay where they were. `faces` picks the face to rebuild (its first, when none). + SurfaceMatch { + body: BodyName, + #[serde(default, skip_serializing_if = "Vec::is_empty")] + faces: Vec, + target: BodyName, + edges: Vec, + #[serde(default = "match_order")] + order: Continuity, + }, + /// The B-spline faces at `faces` of `body` faired (`surface.fair`): each face's interior + /// control vertices move to flatten the thin-plate (curvature-variation) energy of the face, + /// `strength` of the full least-squares step at a time, for at most `iterations` steps. + /// `keep_edges` holds the row of control points next to each edge as well as the boundary row + /// itself, so the face keeps meeting its neighbours in position and in tangent. `tolerance` + /// is the farthest a control point may move, so the faired face stays within it of the one the + /// design asked for; an empty `faces` fairs every face of the body. + SurfaceFair { + body: BodyName, + #[serde(default, skip_serializing_if = "Vec::is_empty")] + faces: Vec, + #[serde(default = "fair_strength")] + strength: f64, + #[serde(default = "fair_iterations")] + iterations: u32, + #[serde(default = "fair_keep_edges")] + keep_edges: bool, + #[walk(length = "Tolerance")] + tolerance: Expr, + }, /// A surface body trimmed by a plane, or by a face of a body (`tool`), keeping the side of /// `keep`. SurfaceTrim { @@ -1600,10 +1653,33 @@ fn is_g0(g: &Continuity) -> bool { } /// The continuity a blend repeats by default: the tangent plane on both sides. +/// The stations a boundary patch carries by default. +fn boundary_samples() -> usize { + 25 +} + fn blend_orders() -> [Continuity; 2] { [Continuity::G1, Continuity::G1] } +/// The continuity a match holds by default: position, tangent plane and curvature. +fn match_order() -> Continuity { + Continuity::G2 +} + +/// A fairing's defaults for a design that leaves them out: the whole least-squares step at a +/// time, and the boundary held (the row next to each edge too, so the face still meets its +/// neighbours in tangent). +fn fair_strength() -> f64 { + 1.0 +} +fn fair_iterations() -> u32 { + 1 +} +fn fair_keep_edges() -> bool { + true +} + /// Triangles of one imported mesh body (millimetres). #[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] pub struct MeshData { @@ -1896,6 +1972,9 @@ impl FeatureKind { FeatureKind::SurfaceSweep { .. } => "SurfaceSweepFeature", FeatureKind::SurfaceLoft { .. } => "SurfaceLoftFeature", FeatureKind::SurfaceBlend { .. } => "BlendFeature", + FeatureKind::SurfaceBoundary { .. } => "BoundaryFeature", + FeatureKind::SurfaceMatch { .. } => "MatchFeature", + FeatureKind::SurfaceFair { .. } => "FairFeature", FeatureKind::Stitch { .. } => "StitchFeature", FeatureKind::Thicken { .. } => "ThickenFeature", FeatureKind::SurfaceTrim { .. } => "TrimFeature", @@ -1973,6 +2052,9 @@ impl FeatureKind { FeatureKind::SurfaceSweep { .. } => "SurfaceSweep", FeatureKind::SurfaceLoft { .. } => "SurfaceLoft", FeatureKind::SurfaceBlend { .. } => "Blend", + FeatureKind::SurfaceBoundary { .. } => "Boundary", + FeatureKind::SurfaceMatch { .. } => "Match", + FeatureKind::SurfaceFair { .. } => "Fair", FeatureKind::Stitch { .. } => "Stitch", FeatureKind::Thicken { .. } => "Thicken", FeatureKind::SurfaceTrim { .. } => "Trim", diff --git a/crates/doc/src/eval.rs b/crates/doc/src/eval.rs index d2b52319..12003ecd 100644 --- a/crates/doc/src/eval.rs +++ b/crates/doc/src/eval.rs @@ -4133,6 +4133,9 @@ fn eval_feature(doc: &Document, vals: &BTreeMap, f: &Feature, st: | FeatureKind::SurfaceSweep { .. } | FeatureKind::SurfaceLoft { .. } | FeatureKind::SurfaceBlend { .. } + | FeatureKind::SurfaceBoundary { .. } + | FeatureKind::SurfaceMatch { .. } + | FeatureKind::SurfaceFair { .. } | FeatureKind::SurfaceOffset { .. } | FeatureKind::SurfaceRuled { .. } | FeatureKind::SurfaceUntrim { .. } => surface::eval(doc, vals, f, st, warning), diff --git a/crates/doc/src/eval_surface.rs b/crates/doc/src/eval_surface.rs index 19a1d35f..d6c8a573 100644 --- a/crates/doc/src/eval_surface.rs +++ b/crates/doc/src/eval_surface.rs @@ -145,6 +145,78 @@ pub(super) fn eval(doc: &Document, vals: &BTreeMap, f: &Feature, st.bodies.push(ModelBody::new(name, made, f.id)); Ok(()) } + FeatureKind::SurfaceBoundary { sketch, curves, orders, crown, samples } => { + let ss = sketch_in_feature_frame( + doc, + st.sketch(sketch.0).ok_or_else(|| DocError::Unknown(format!("sketch {sketch} (it must come earlier in the timeline)")))?, + f.component, + )?; + let boundaries = boundary_frames(&ss, curves, orders, *samples)?; + let body = kernel::boundary_patch(&boundaries, *crown, *samples)?; + let name = new_name(st, f, 0); + st.bodies.push(ModelBody::new(name, body, f.id)); + Ok(()) + } + FeatureKind::SurfaceMatch { body, faces, target, edges, order } => { + 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 tb = st.body(target).ok_or_else(|| DocError::Unknown(format!("body `{target}`")))?.body.clone(); + if edges.is_empty() { + return Err(DocError::Invalid("a match needs a chain of the target's edges".into())); + } + // Picked by a point on the face; without picks, the body's first face. + let picked: Vec = if faces.is_empty() { + let faces = mb.body.faces(mb.body.size() * 1e-2)?; + vec![faces.first().ok_or_else(|| DocError::Invalid(format!("`{body}` has no face to match")))?.centroid] + } else { + crate::plain(faces) + }; + let made = kernel::match_surface(&mb.body, &picked, &tb, &crate::plain(edges), *order)?; + // A match rebuilds the face in place: the body keeps its name, and whatever meets it + // downstream still finds it. + if let Some(slot) = st.bodies.get_mut(i) { + *slot = ModelBody::new(mb.name, made, mb.feature); + } + Ok(()) + } + FeatureKind::SurfaceFair { body, faces, strength, iterations, keep_edges, tolerance } => { + let tol = val(vals, tolerance, Kind::Length)?; + let opts = kernel::FairOptions::new(*strength, *iterations as usize, *keep_edges, tol)?; + 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())) }; + // Picked by a point on the face; without picks, every face of the body. + let mut picked: Vec = Vec::new(); + if faces.is_empty() { + picked.extend(0..mb.body.face_count()); + } else { + for p in faces { + let fi = + crate::appearance::face_index_at(&mb, p.0).ok_or_else(|| DocError::Invalid(format!("Fair: no face of {body} at {p:?}")))?; + if !picked.contains(&fi) { + picked.push(fi); + } + } + } + if picked.is_empty() { + return Err(DocError::Invalid("the body has no faces to fair".into())); + } + // The face keeps every edge it had: only the surface between them is faired. The + // energy before and after goes into the feature's warning, so the reply can say how + // much the face moved. + let (faired, report) = kernel::fair_faces(&mb.body, &picked, &opts)?; + *warning = Some(format!( + "faired {} face(s): energy {:.6} → {:.6} in {} step(s), control points moved up to {:.4} mm", + picked.len(), + report.energy_before, + report.energy_after, + report.steps, + report.movement + )); + if let Some(slot) = st.bodies.get_mut(i) { + *slot = ModelBody::new(mb.name, faired, mb.feature); + } + Ok(()) + } FeatureKind::SurfaceExtrude { sketch, profiles, curves, extent } => { let ss = sketch_in_feature_frame( doc, @@ -349,6 +421,70 @@ fn push_sheet(st: &mut ModelState, f: &Feature, traced: crate::lineage::Traced) st.bodies.push(ModelBody::traced_new(name, traced, f.id)); } +/// The boundaries `surface.boundary` builds its patch from: each picked curve's own stations in +/// the sketch's frame, every one of them leaving the boundary square into the loop — the direction +/// from the curve towards the middle of the loop, taken across its tangent — over the plane the +/// curve lies in (so a chain of planar curves makes a flat patch, and a network that is not +/// coplanar is followed as the plane through the curve's own tangent). `orders` sets how much of +/// that plane each side repeats. +fn boundary_frames( + ss: &SolvedSketch, + curves: &[String], + orders: &[solidcraft_kernel::Continuity], + samples: usize, +) -> Result> { + use solidcraft_kernel::{Boundary, Frame}; + // The stations along each curve, evenly in its own arclength, so every side carries the same + // number of frames and the patch's net is square. + let mut sides: Vec> = Vec::with_capacity(curves.len()); + for id in curves { + let ci = ss.sketch.curve_index(id).ok_or_else(|| DocError::Unknown(format!("sketch curve `{id}`")))?; + let pieces = ss.sketch.exact_pieces(ci); + if pieces.is_empty() { + return Err(DocError::Invalid(format!("curve `{id}` has no exact geometry to build on"))); + } + let lengths: Vec = pieces.iter().map(Seg2::length).collect(); + let total: f64 = lengths.iter().sum(); + if !total.is_finite() || total <= 1e-9 { + return Err(DocError::Invalid(format!("curve `{id}` has no length"))); + } + let mut pts = Vec::with_capacity(samples); + let (mut used, mut at) = (0.0, 0usize); + for k in 0..samples { + let want = total * k as f64 / (samples.saturating_sub(1)).max(1) as f64; + while at + 1 < pieces.len() && used + lengths[at] < want { + used += lengths[at]; + at += 1; + } + let (Some(piece), Some(len)) = (pieces.get(at), lengths.get(at)) else { return Err(DocError::Invalid("boundary piece".into())) }; + let t = ((want - used) / len.max(1e-12)).clamp(0.0, 1.0); + pts.push(ss.plane.to_world(piece.point_at(t))); + } + sides.push(pts); + } + // The middle of the loop, so every side leaves towards the inside of the patch. + let count = sides.iter().map(Vec::len).sum::().max(1) as f64; + let centre = sides.iter().flatten().fold(solidcraft_geom::Vec3::ZERO, |a, p| a + *p) * (1.0 / count); + let mut out = Vec::with_capacity(sides.len()); + for (side, pts) in sides.iter().enumerate() { + let mut frames = Vec::with_capacity(pts.len()); + for (i, p) in pts.iter().enumerate() { + // The tangent from the stations on either side, one-sided at the ends (where the side + // only has one neighbour). + let before = if i == 0 { *p } else { pts[i - 1] }; + let after = if i + 1 == pts.len() { *p } else { pts[i + 1] }; + let tangent = (after - before).normalized().ok_or_else(|| DocError::Invalid("a boundary curve has a stationary station".into()))?; + let across = centre - *p; + let leave = (across - tangent * across.dot(tangent)) + .normalized() + .ok_or_else(|| DocError::Invalid("a boundary runs through the middle of the loop: no patch can leave it".into()))?; + frames.push(Frame { p: *p, leave, bend: Vec3::ZERO, bend_change: Vec3::ZERO }); + } + out.push(Boundary { frames, order: orders.get(side).copied().unwrap_or(solidcraft_kernel::Continuity::G1) }); + } + Ok(out) +} + /// What a surface builder sweeps: the boundary loops of the picked profiles, then the picked /// curves chained end to end. fn sheet_chains(ss: &SolvedSketch, profiles: Option<&ProfileSel>, curves: &[String]) -> Result>> { diff --git a/crates/doc/src/feature_walk_tests/fixtures.rs b/crates/doc/src/feature_walk_tests/fixtures.rs index 6a23ee71..2364cf6c 100644 --- a/crates/doc/src/feature_walk_tests/fixtures.rs +++ b/crates/doc/src/feature_walk_tests/fixtures.rs @@ -161,6 +161,11 @@ pub(super) fn fixtures() -> Vec<(&'static str, Value)> { ("surface_loft", json!({"type": "surface_loft", "sections": []})), ("surface_blend", json!({"type": "surface_blend", "from": {"body": "B", "edges": [o]}, "to": {"body": "C", "edges": [o]}, "align": [o, o]})), ("surface_offset", json!({"type": "surface_offset", "faces": [o], "distance": "soff.distance", "tolerance": "soff.tolerance"})), + // The Class-A surface family (surface.blend/boundary/match/fair). + ("surface_blend", json!({"type": "surface_blend", "from": {"body": "B", "edges": [o]}, "to": {"body": "C", "edges": [o]}, "align": [o, o]})), + ("surface_boundary", json!({"type": "surface_boundary", "sketch": 1, "curves": ["l1", "l2"]})), + ("surface_match", json!({"type": "surface_match", "body": "B", "faces": [o], "target": "C", "edges": [o]})), + ("surface_fair", json!({"type": "surface_fair", "body": "B", "faces": [o], "tolerance": "sfair.tolerance"})), ("unstitch", json!({"type": "unstitch", "body": "B"})), ("reverse_normal", json!({"type": "reverse_normal", "body": "B"})), ("surface_untrim", json!({"type": "surface_untrim", "face": {"body": "B", "point": o}, "holes": [o], "outer": true})), diff --git a/crates/doc/src/feature_walk_tests/golden_exprs.txt b/crates/doc/src/feature_walk_tests/golden_exprs.txt index 7a060857..5a887777 100644 --- a/crates/doc/src/feature_walk_tests/golden_exprs.txt +++ b/crates/doc/src/feature_walk_tests/golden_exprs.txt @@ -340,6 +340,13 @@ surface_offset Distance = soff.distance (Length) Tolerance = soff.tolerance (Length) expressions: soff.distance, soff.tolerance +surface_blend +surface_boundary +surface_match +surface_fair + inputs: + Tolerance = sfair.tolerance (Length) + expressions: sfair.tolerance unstitch reverse_normal surface_untrim diff --git a/crates/doc/src/feature_walk_tests/golden_geometry.txt b/crates/doc/src/feature_walk_tests/golden_geometry.txt index 32f8592e..767444fa 100644 --- a/crates/doc/src/feature_walk_tests/golden_geometry.txt +++ b/crates/doc/src/feature_walk_tests/golden_geometry.txt @@ -154,6 +154,15 @@ surface_blend surface_offset moves: faces[0][0], faces[0][1], faces[0][2] face points: [0, 0, 0] +surface_blend + moves: align[0][0], align[0][1], align[0][2], align[1][0], align[1][1], align[1][2], from.edges[0][0], from.edges[0][1], from.edges[0][2], to.edges[0][0], to.edges[0][1], to.edges[0][2] +surface_boundary +surface_match + moves: edges[0][0], edges[0][1], edges[0][2], faces[0][0], faces[0][1], faces[0][2] + face points: [0, 0, 0] +surface_fair + moves: faces[0][0], faces[0][1], faces[0][2] + face points: [0, 0, 0] unstitch moves: faces[0][0], faces[0][1], faces[0][2] reverse_normal diff --git a/crates/doc/src/feature_walk_tests/golden_refs.txt b/crates/doc/src/feature_walk_tests/golden_refs.txt index 9c56ab0b..0b621c6a 100644 --- a/crates/doc/src/feature_walk_tests/golden_refs.txt +++ b/crates/doc/src/feature_walk_tests/golden_refs.txt @@ -177,6 +177,15 @@ surface_blend remap_bodies: from.body surface_offset remap_bodies: body +surface_blend + remap_bodies: from.body +surface_boundary + sketches: [1] + remap_refs: sketch +surface_match + remap_bodies: body +surface_fair + remap_bodies: body unstitch remap_bodies: body reverse_normal diff --git a/crates/doc/src/format/pending/class-a-surfaces.rs b/crates/doc/src/format/pending/class-a-surfaces.rs new file mode 100644 index 00000000..ee4e5859 --- /dev/null +++ b/crates/doc/src/format/pending/class-a-surfaces.rs @@ -0,0 +1,5 @@ +//! The Class-A surface pass adds `SurfaceBoundary` (Boundary Surface: `sketch`, `curves`, +//! optional `orders`, `crown`, `samples`), `SurfaceMatch` (Match Surface: `body`, optional +//! `faces`, `target`, `edges`, optional `order`) and `SurfaceFair` (Fair Surface: `body`, optional +//! `faces`, `strength`, `iterations`, `keep_edges`, `tolerance`). Older designs need no +//! transformation; 0.13.x asks for an update instead of meeting a feature it can't read. diff --git a/crates/doc/src/walk.rs b/crates/doc/src/walk.rs index 64073b31..7c571a1b 100644 --- a/crates/doc/src/walk.rs +++ b/crates/doc/src/walk.rs @@ -40,6 +40,8 @@ impl Walk for bool {} impl Walk for u64 {} impl Walk for i8 {} impl Walk for f64 {} +impl Walk for usize {} +impl Walk for u32 {} impl Walk for ProfileSel {} impl Walk for solidcraft_sketch::Sketch {} impl Walk for solidcraft_sketch::PatternAlignment {} diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index db7fa5c3..a37f226a 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -74,6 +74,7 @@ mod sketch_modify; mod sketch_project; mod surface; mod surface_analysis; +mod surface_match; mod view_dxf; mod view_navigation; @@ -236,6 +237,7 @@ pub fn command_specs() -> Vec<&'static CommandSpec> { v.extend(sketch_3d::COMMANDS.iter()); v.extend(surface_analysis::COMMANDS.iter()); v.extend(surface::COMMANDS.iter()); + v.extend(surface_match::COMMANDS.iter()); v.extend(sketch_modify::COMMANDS.iter()); v.extend(sketch_project::COMMANDS.iter()); v.extend(features::COMMANDS.iter()); diff --git a/crates/engine/src/cmd/surface.rs b/crates/engine/src/cmd/surface.rs index 93d7a19a..4e76e0d4 100644 --- a/crates/engine/src/cmd/surface.rs +++ b/crates/engine/src/cmd/surface.rs @@ -3,7 +3,8 @@ use serde_json::Value; use solidcraft_doc::{ - BlendSide, EdgeContinuity, ExtendMode, ExtendTo, FaceAt, FeatureKind, PatchContinuity, PlaneRef, ProfileSel, RuledCorners, RuledMode, expr::Kind, + BlendSide, EdgeContinuity, Expr, ExtendMode, ExtendTo, FaceAt, FeatureKind, PatchContinuity, PlaneRef, ProfileSel, RuledCorners, RuledMode, + expr::Kind, }; use solidcraft_doc::{BodyName, SketchId}; use solidcraft_doc::{Dir, FacePick, Position}; @@ -11,7 +12,7 @@ use solidcraft_kernel::Continuity; use super::CommandSpec; use super::features::{add_feature, axis_param, check_expr, extent_param, feature_sketch, operation, plane_param, profiles, sketch_id}; -use crate::params::{bad, bool_, expr, str_, string_list, vec3}; +use crate::params::{bad, bool_, expr, num, str_, string_list, vec3}; use crate::{Result, Session}; pub static COMMANDS: &[CommandSpec] = &[ @@ -59,10 +60,6 @@ pub static COMMANDS: &[CommandSpec] = &[ .at("SURFACE", "CREATE") .icon("ruled") .params("body; edges: [[x,y,z]] (points on a chain of its edges, open or closed); distance: expr (how far the strip reaches); mode?: normal|tangent|angle|direction (default normal: along the normal of the face beside each edge; tangent: continuing it); angle?: expr (angle mode: from `direction`, else from the face's normal, turned away from the face; 90 deg from the normal is tangent); direction?: [x,y,z] (direction mode: the rulings; angle mode: the reference, e.g. a pull direction for a draft); face?: [x,y,z] (a point on the face the directions come from; default: the face beside the most picked edges); flip?: bool (the other side); corners?: mitre|round (sharp corners of the chain; inside corners are always cut back). Lines and arcs give exact planes, cylinders and cones, other edges a fitted ruled B-spline. A new surface body"), - planned!("surface.boundary", "Boundary Surface", "763") - .at("SURFACE", "CREATE") - .icon("boundary") - .params("planned: u: [curves], v?: [curves] (sketch curves, 3D sketch wires or {body, point} edges; 2–4 boundaries, more make a network); continuity?: [{curve, g: 0|1|2}] for edge boundaries; operation?: new|join"), CommandSpec::new("surface.blend", "Blend Surface", blend) .at("SURFACE", "CREATE") .icon("surface_blend") @@ -79,6 +76,10 @@ pub static COMMANDS: &[CommandSpec] = &[ .at("SURFACE", "MODIFY") .icon("reverse_normal") .params("body (a surface); faces?: [[x,y,z]] (points on the faces to turn over; default: all of them)"), + CommandSpec::new("surface.fair", "Fair Surface", fair) + .at("SURFACE", "MODIFY") + .icon("spline_cv") + .params("body: surface name; faces?: [[x,y,z]] points on the faces to fair (default: every face); strength?: number (0..1 of the full fairing step, default 1); iterations?: number (1..1000 steps, default 1); keep_edges?: bool (default true: hold the boundary row and the row next to it, so a face keeps meeting its neighbours in position and in tangent); tolerance?: expr mm (how far a control vertex may move, default 0.001 of the body's size). Moves the interior control vertices of each B-spline face to flatten its curvature variation while the boundary stays put: the face keeps its edges and the body its shape within the tolerance, and the warning in the reply carries the energy before and after; name?"), planned!("surface.fillet", "Surface Fillet", "769") .at("SURFACE", "MODIFY") .icon("surface_fillet") @@ -558,6 +559,39 @@ fn trim(s: &mut Session, p: &Value) -> Result { add_feature(s, p, FeatureKind::SurfaceTrim { body: BodyName(body), plane, tool, keep: Position(keep) }) } +/// `surface.fair`: the B-spline faces of a body faired in place. The tolerance a caller leaves +/// out is a thousandth of the body's own size, so a fairing run can never make a design jump by +/// more than a hair of itself. +fn fair(s: &mut Session, p: &Value) -> Result { + let cmd = "surface.fair"; + let (body, faces) = body_faces(s, p, cmd)?; + let size = { + let state = s.model.state(); + let Some(mb) = state.body(&body) else { return Err(bad(cmd, format!("no body `{body}`"))) }; + mb.body.size() + }; + let strength = num(p, "strength").unwrap_or(1.0); + if !(strength.is_finite() && strength > 0.0 && strength <= 1.0) { + return Err(bad(cmd, "`strength` must be greater than 0 and at most 1 (the whole fairing step)")); + } + let steps = num(p, "iterations").unwrap_or(1.0); + if !(steps.is_finite() && (1.0..=1000.0).contains(&steps)) { + return Err(bad(cmd, "`iterations` must be 1 to 1000")); + } + let keep_edges = bool_(p, "keep_edges").unwrap_or(true); + let tolerance = expr(p, "tolerance").unwrap_or_else(|| Expr::new(format!("{:.6}", (size * 1e-3).max(1e-6)))); + check_expr(s, tolerance.as_str(), Kind::Length, cmd, "tolerance")?; + let kind = FeatureKind::SurfaceFair { + body: BodyName::from(body), + faces: faces.into_iter().map(FacePick).collect(), + strength, + iterations: steps as u32, + keep_edges, + tolerance, + }; + add_feature(s, p, kind) +} + #[cfg(test)] mod tests { use super::*; @@ -1178,7 +1212,7 @@ mod tests { /// as it was. #[test] fn planned_commands_refuse_and_leave_the_design_alone() { - let planned = ["surface.boundary", "surface.fillet", "inspect.continuity"]; + let planned = ["surface.fillet", "inspect.continuity"]; let mut s = Session::default(); run(&mut s, "solid.box", json!({"length": 10, "width": 10, "height": 10, "body_name": "B"})); let before = s.doc.features.len(); @@ -1500,4 +1534,58 @@ mod tests { let sketch_patch = json!({"sketch": 1, "continuity": "g1"}); assert!(s.execute("surface.patch", &sketch_patch).is_err(), "continuity on a sketch patch"); } + /// The energies the reply's warning carries: the fairness energy before and after. + fn energies(reply: &Value) -> (f64, f64) { + let warning = reply["warning"].as_str().unwrap_or_else(|| panic!("no warning in {reply}")); + let after = warning.split("energy ").nth(1).unwrap_or_else(|| panic!("no energies in `{warning}`")); + let numbers: Vec = after + .split(|c: char| !(c.is_ascii_digit() || c == '.' || c == 'e' || c == '-' || c == '+')) + .filter_map(|w| w.parse::().ok()) + .collect(); + assert!(numbers.len() >= 2, "no energies in `{warning}`"); + (numbers[0], numbers[1]) + } + + /// Surface Fair on the fitted B-spline face of an offset sphere: the face is smoothed in + /// place, its edges do not move at all, the reply carries the energy before and after, undo, + /// redo and a cold reopen keep it, and hostile inputs are refused. + #[test] + fn fairing_a_b_spline_face_keeps_its_edges_and_reports_its_energy() { + let mut s = Session::default(); + run(&mut s, "solid.sphere", json!({"radius": 20, "body_name": "Ball"})); + // A sphere is neither a plane nor a cylinder, so its offset is a fitted B-spline face. + run(&mut s, "surface.offset", json!({"faces": [[0, 0, 20]], "body": "Ball", "distance": 2, "body_name": "Shell"})); + let st = s.world_state(); + assert!(st.body("Shell").unwrap().body.is_surface(), "an offset sphere face is a sheet"); + let before: Vec = st.body("Shell").unwrap().body.edges(0.1).unwrap().iter().flat_map(|e| e.points.clone()).collect(); + assert!(!before.is_empty()); + let area_before = area(&s, "Shell"); + // The fairing reports both energies, the second no larger than the first. + let reply = run(&mut s, "surface.fair", json!({"body": "Shell", "tolerance": 1.0, "iterations": 3})); + let (energy_before, energy_after) = energies(&reply); + let area_after = area(&s, "Shell"); + eprintln!("surface.fair: {reply}\nenergies {energy_before:e} → {energy_after:e}, area {area_before} → {area_after}"); + assert!(energy_before > 0.0 && energy_after <= energy_before, "{energy_before:e} → {energy_after:e}"); + // The face moved (the fit is not the flattest surface on its own rim), but not its edges: + // every boundary point of the body is where it was, bit for bit. + assert!((area_after - area_before).abs() > 0.0 || energy_after < energy_before); + let after: Vec = st.body("Shell").unwrap().body.edges(0.1).unwrap().iter().flat_map(|e| e.points.clone()).collect(); + assert_eq!(after, before, "fairing must not move the face's edges"); + assert_eq!(st.body("Shell").unwrap().body.face_count(), 1); + // Undo takes the fairing back, redo puts it back, and a cold reopen rebuilds the same. + run(&mut s, "edit.undo", json!({})); + assert!((area(&s, "Shell") - area_before).abs() < 1e-9, "undo did not take the fairing back"); + run(&mut s, "edit.redo", json!({})); + assert!((area(&s, "Shell") - area_after).abs() < 1e-9, "redo did not put the fairing back"); + let cold = Session::new(solidcraft_doc::Document::from_json(&s.doc.to_json()).unwrap()); + assert!((area(&cold, "Shell") - area_after).abs() < 1e-9, "a cold reopen did not rebuild the same sheet"); + // Refusals: an unknown body, a body whose faces are all planes, an out-of-range strength + // or step count, and a missing tolerance expression. + assert!(s.execute("surface.fair", &json!({"body": "Nope", "tolerance": 1.0})).is_err()); + run(&mut s, "solid.box", json!({"length": 10, "width": 10, "height": 10, "body_name": "Box"})); + assert!(s.execute("surface.fair", &json!({"body": "Box", "tolerance": 1.0})).is_err()); + assert!(s.execute("surface.fair", &json!({"body": "Shell", "strength": 2, "tolerance": 1.0})).is_err()); + assert!(s.execute("surface.fair", &json!({"body": "Shell", "iterations": 0, "tolerance": 1.0})).is_err()); + assert!(s.execute("surface.fair", &json!({"body": "Shell", "tolerance": 0})).is_err()); + } } diff --git a/crates/engine/src/cmd/surface_boundary.rs b/crates/engine/src/cmd/surface_boundary.rs new file mode 100644 index 00000000..6bbc3f1a --- /dev/null +++ b/crates/engine/src/cmd/surface_boundary.rs @@ -0,0 +1,197 @@ +//! Boundary surfaces: one patch through a network of two to four sketch curves (`surface.boundary`, +//! issue #763). The kernel builds the face from the curves' own frames; this is the command's +//! parameters and its checks, so that a loop it cannot hold is refused before anything is built. + +use serde_json::Value; +use solidcraft_doc::{FeatureKind, SketchId}; +use solidcraft_kernel::Continuity; + +use super::CommandSpec; +use super::features::{add_feature, feature_sketch}; +use crate::params::{bad, string_list}; +use crate::{Result, Session}; + +/// The stations a boundary carries by default (the patch's resolution along the loop). +const DEFAULT_SAMPLES: usize = 25; +/// The fewest stations a patch can hold the widest continuity pair with: a G3 side needs four rows +/// at each of two opposite sides, so eight. +const MIN_SAMPLES: usize = 8; +/// The most stations a boundary may carry (the kernel's own cap). +const MAX_SAMPLES: usize = 512; + +pub static COMMANDS: &[CommandSpec] = &[CommandSpec::new("surface.boundary", "Boundary Surface", boundary) + .at("SURFACE", "CREATE") + .icon("boundary") + .params( + "sketch?: id|name (default: the sketch being edited); curves: [2 to 4 sketch curve ids in loop order] (each one ends where the next begins and the last where the first does; the patch runs through them and leaves each boundary along the tangent across it); continuity?: g0|g1|g2|g3 | [one per curve] (default g1: how much of the surface across each boundary the patch repeats — a sketch curve's own plane, so g1 and up hold exactly and g0 lets the middle bow); crown?: number -0.9…3 (default 0: how full the middle is, without touching a boundary); samples?: 8…512 stations per boundary (default 25); name?, body_name? (a surface body)", + )]; + +/// `surface.boundary`: one patch through a network of two to four sketch curves. +fn boundary(s: &mut Session, p: &Value) -> Result { + let cmd = "surface.boundary"; + let sketch = feature_sketch(s, p, cmd)?; + let curves = string_list(p, "curves"); + if !(2..=4).contains(&curves.len()) { + return Err(bad(cmd, "`curves` must list 2 to 4 sketch curves in loop order")); + } + // The curves must be the sketch's own, and they must run on into one another: the last station + // of each is the first of the next, and the last boundary's end is the first one's start. + let state = s.model.state(); + let sk = state.sketch(sketch).ok_or_else(|| bad(cmd, "the sketch has not been solved"))?; + let mut ends: Vec<(solidcraft_geom::Vec2, solidcraft_geom::Vec2)> = Vec::with_capacity(curves.len()); + for id in &curves { + let ci = sk.sketch.curve_index(id).ok_or_else(|| bad(cmd, format!("no sketch curve `{id}`")))?; + let pieces = sk.sketch.exact_pieces(ci); + let (Some(first), Some(last)) = (pieces.first(), pieces.last()) else { + return Err(bad(cmd, format!("curve `{id}` has no exact geometry to build on"))); + }; + if last.end().dist(first.start()) < 1e-9 { + return Err(bad(cmd, format!("curve `{id}` is closed: a boundary must be one side of the loop"))); + } + ends.push((first.start(), last.end())); + } + // Two boundaries are the two sides of a bridge, which do not meet; three or four run round a + // loop, so each has to end where the next begins. + let scale = ends.iter().map(|(a, b)| a.dist(*b)).fold(1.0_f64, f64::max); + let tol = scale * 1e-4 + 1e-9; + for i in 0..if curves.len() >= 3 { ends.len() } else { 0 } { + let next = (i + 1) % ends.len(); + let (Some(end), Some(start)) = (ends.get(i).map(|e| e.1), ends.get(next).map(|e| e.0)) else { + return Err(bad(cmd, "a curve is missing its ends")); + }; + let gap = end.dist(start); + if gap > tol { + return Err(bad( + cmd, + format!("curve `{}` ends {gap:.4} mm from where `{}` starts: list the curves in loop order", curves[i], curves[next]), + )); + } + } + drop(state); + let orders = match p.get("continuity") { + None | Some(Value::Null) => vec![Continuity::G1; curves.len()], + Some(Value::String(w)) => { + let one = Continuity::parse(w).ok_or_else(|| bad(cmd, "`continuity` must be g0, g1, g2 or g3"))?; + vec![one; curves.len()] + } + Some(v) => { + let list = v.as_array().filter(|a| a.len() == curves.len()).ok_or_else(|| bad(cmd, "`continuity` must be one word, or one per curve"))?; + let mut out = Vec::with_capacity(list.len()); + for x in list { + out.push(x.as_str().and_then(Continuity::parse).ok_or_else(|| bad(cmd, "`continuity` entries must be g0, g1, g2 or g3"))?); + } + out + } + }; + let crown = match p.get("crown") { + None | Some(Value::Null) => 0.0, + Some(v) => v + .as_f64() + .filter(|x| x.is_finite() && (solidcraft_geom::fill::MIN_CROWN..=solidcraft_geom::fill::MAX_CROWN).contains(x)) + .ok_or_else(|| bad(cmd, "`crown` must be a number between -0.9 and 3"))?, + }; + let samples = match p.get("samples") { + None | Some(Value::Null) => DEFAULT_SAMPLES, + Some(v) => { + let n = v.as_f64().filter(|x| x.is_finite() && x.fract() == 0.0).map(|x| x as usize); + n.filter(|x| (MIN_SAMPLES..=MAX_SAMPLES).contains(x)) + .ok_or_else(|| bad(cmd, format!("`samples` must be a whole number of stations between {MIN_SAMPLES} and {MAX_SAMPLES}")))? + } + }; + add_feature(s, p, FeatureKind::SurfaceBoundary { sketch: SketchId(sketch), curves, orders, crown, samples }) +} + +#[cfg(test)] +mod tests { + use serde_json::{Value, json}; + + use crate::Session; + + fn run(s: &mut Session, id: &str, p: Value) -> Value { + match s.execute(id, &p) { + Ok(v) => v, + Err(e) => panic!("{id} {p}: {e}"), + } + } + + fn ids(v: &Value) -> Vec { + v.as_array().map(|a| a.iter().filter_map(|x| x.as_str().map(str::to_string)).collect()).unwrap_or_default() + } + + fn new_sketch() -> Session { + let mut s = Session::default(); + run(&mut s, "sketch.create", json!({"plane": "XY"})); + s + } + + /// Two parallel 20 mm lines 10 mm apart, and the patch between them. + fn two_walls(s: &mut Session) -> Vec { + let a = ids(&run(s, "sketch.line", json!({"points": [[0, 0], [20, 0]]}))["created"]["curves"])[0].clone(); + let b = ids(&run(s, "sketch.line", json!({"points": [[0, 10], [20, 10]]}))["created"]["curves"])[0].clone(); + vec![a, b] + } + + #[test] + fn a_boundary_patch_between_two_sketch_curves_is_the_sheet_they_span() { + for continuity in [json!("g0"), json!("g1"), json!("g2"), json!(["g1", "g2"])] { + let mut s = new_sketch(); + let curves = two_walls(&mut s); + let r = run(&mut s, "surface.boundary", json!({"curves": curves, "continuity": continuity})); + let body = ids(&r["created_bodies"])[0].clone(); + // It is one face, its boundary is the two lines themselves, and it covers exactly the + // 20 × 10 sheet between them — all four frames lie in the sketch's plane. + let f = run(&mut s, "model.faces", json!({"body": body})); + let faces = f["faces"].as_array().cloned().unwrap_or_default(); + assert_eq!(faces.len(), 1, "{continuity}: {f}"); + assert!((faces[0]["area"].as_f64().unwrap_or(f64::NAN) - 200.0).abs() < 1e-4, "{continuity}: {f}"); + let normal = faces[0]["plane_normal"].as_array().cloned().unwrap_or_default(); + assert_eq!(normal.len(), 3, "{continuity}: the patch is planar: {f}"); + assert!(normal[2].as_f64().unwrap_or(0.0).abs() > 1.0 - 1e-9, "{continuity}: it lies in the sketch's plane: {f}"); + let m = run(&mut s, "inspect.measure", json!({})); + let area = m["total"]["area_mm2"].as_f64().unwrap_or(f64::NAN); + assert!((area - 200.0).abs() < 1e-4, "{continuity}: {area}"); + } + } + + #[test] + fn a_boundary_patch_refuses_a_loop_and_values_it_cannot_hold() { + let mut s = new_sketch(); + let curves = two_walls(&mut s); + let extra = ids(&run(&mut s, "sketch.line", json!({"points": [[30, 30], [40, 40]]}))["created"]["curves"])[0].clone(); + // A loop of three or four: the kernel refuses it with its own reason (a tensor-product patch + // cannot hold two sides' data at a corner, issue #763) rather than building it approximately. + let triangle = ["l1", "l2", "l3"].map(|id| { + let _ = id; + String::new() + }); + let _ = triangle; + let a = ids(&run(&mut s, "sketch.line", json!({"points": [[0, 0], [20, 0]]}))["created"]["curves"])[0].clone(); + let b = ids(&run(&mut s, "sketch.line", json!({"points": [[20, 0], [10, 15]]}))["created"]["curves"])[0].clone(); + let c = ids(&run(&mut s, "sketch.line", json!({"points": [[10, 15], [0, 0]]}))["created"]["curves"])[0].clone(); + let e = s.execute("surface.boundary", &json!({"curves": [a, b, c], "continuity": "g2"})).unwrap_err().to_string(); + assert!(e.contains("not supported yet") && e.contains("#763"), "three boundaries: {e}"); + // A closed square of four lines, in loop order. + let sq = [[[0, 0], [10, 0]], [[10, 0], [10, 10]], [[10, 10], [0, 10]], [[0, 10], [0, 0]]]; + let square: Vec = sq.iter().map(|p| ids(&run(&mut s, "sketch.line", json!({"points": p}))["created"]["curves"])[0].clone()).collect(); + let e = s.execute("surface.boundary", &json!({"curves": square})).unwrap_err().to_string(); + assert!(e.contains("not supported yet") && e.contains("#763"), "four boundaries: {e}"); + // Parameters and curves that make no sense. + let cases: Vec<(Value, &str)> = vec![ + (json!({"curves": [curves[0]]}), "2 to 4"), + (json!({"curves": [curves[0], curves[1], extra.clone(), extra.clone(), extra.clone()]}), "2 to 4"), + (json!({"curves": [curves[0], "nope"]}), "no sketch curve"), + (json!({"curves": curves, "continuity": "g9"}), "continuity"), + (json!({"curves": curves, "continuity": ["g1", "g2", "g3"]}), "continuity"), + (json!({"curves": curves, "crown": 5}), "crown"), + (json!({"curves": curves, "crown": "full"}), "crown"), + (json!({"curves": curves, "samples": 3}), "samples"), + (json!({"curves": curves, "samples": 2.5}), "samples"), + ]; + for (params, want) in cases { + let e = s.execute("surface.boundary", ¶ms).unwrap_err().to_string(); + assert!(e.contains(want), "{params}: {e}"); + } + // Nothing was built by any of them. + assert!(s.model.state().bodies.is_empty(), "a refused patch left a body behind"); + } +} diff --git a/crates/engine/src/cmd/surface_match.rs b/crates/engine/src/cmd/surface_match.rs new file mode 100644 index 00000000..78f5b8ec --- /dev/null +++ b/crates/engine/src/cmd/surface_match.rs @@ -0,0 +1,160 @@ +//! `surface.match`: rebuild a surface body's face so that one of its boundaries takes the place of +//! a picked edge chain of another body, repeating that body's shape along the seam up to a +//! requested continuity, while the boundary opposite it keeps its place. The body keeps its name: +//! a match is a modify, like Trim, Extend and Fair. + +use serde_json::Value; +use solidcraft_doc::{BodyName, FacePick, FeatureKind, Position}; +use solidcraft_kernel::Continuity; + +use super::CommandSpec; +use super::features::add_feature; +use crate::params::{bad, str_, vec3}; +use crate::{Result, Session}; + +pub static COMMANDS: &[CommandSpec] = &[ + CommandSpec::new("surface.match", "Match Surface", match_surface) + .at("SURFACE", "MODIFY") + .icon("continuity") + .params( + "body: surface name (the body whose face is rebuilt); faces?: [[x,y,z]] (points on that face, default: the body's first); target: body name (the body the seam must meet); edges: [[x,y,z]] (points on one edge chain of the target: the rebuilt boundary takes its place); order?: g1|g2|g3 (default g2: position and tangent plane, the curvature too, or how that curvature changes along the seam). The far boundary stays where it was and the body keeps its name", + ), +]; + +/// `surface.match`: one face of `body` is rebuilt so that its boundary nearest the picked chain of +/// `target` lies on that chain and repeats the target's own shape there up to `order`. +fn match_surface(s: &mut Session, p: &Value) -> Result { + let cmd = "surface.match"; + let body = str_(p, "body").ok_or_else(|| bad(cmd, "`body` is required"))?.to_string(); + if s.model.state().body(&body).is_none() { + return Err(bad(cmd, format!("no body `{body}`"))); + } + let target = str_(p, "target").ok_or_else(|| bad(cmd, "`target` is required"))?.to_string(); + if s.model.state().body(&target).is_none() { + return Err(bad(cmd, format!("no body `{target}`"))); + } + let faces: Vec = p.get("faces").and_then(Value::as_array).map(|a| a.iter().filter_map(vec3).collect()).unwrap_or_default(); + if faces.len() > 64 { + return Err(bad(cmd, "too many points on the face to rebuild")); + } + let edges: Vec = p.get("edges").and_then(Value::as_array).map(|a| a.iter().filter_map(vec3).collect()).unwrap_or_default(); + if edges.is_empty() || edges.len() > 10_000 { + return Err(bad(cmd, "`edges` must list 1 to 10000 points on the target's edge chain")); + } + let order = match p.get("order") { + None => Continuity::G2, + Some(v) => v.as_str().and_then(Continuity::parse).filter(|o| *o != Continuity::G0).ok_or_else(|| bad(cmd, "`order` must be g1, g2 or g3"))?, + }; + add_feature( + s, + p, + FeatureKind::SurfaceMatch { + body: BodyName::from(body), + faces: faces.into_iter().map(FacePick).collect(), + target: BodyName::from(target), + edges: edges.into_iter().map(Position).collect(), + order, + }, + ) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::Session; + use serde_json::json; + use solidcraft_geom::Vec3; + + fn run(s: &mut Session, c: &str, p: Value) -> Value { + s.execute(c, &p).unwrap_or_else(|e| panic!("{c}: {e}")) + } + + /// The picked chain: two points of the cylinder's own vertical wall seam, clear of the corners + /// at either end. The seam is where `solid.extrude` starts the circle, so it is read off the + /// body rather than assumed. + fn picks(s: &Session) -> Vec { + let boss = s.world_state().body("Boss").unwrap().body.clone(); + let e = boss + .edges(boss.size() * 1e-2) + .unwrap() + .into_iter() + .find(|e| e.length > 25.0 && (e.mid.x.hypot(e.mid.y) - 10.0).abs() < 0.05 && e.mid.z > 5.0 && e.mid.z < 25.0) + .expect("the cylinder has a straight vertical wall seam"); + let (a, b) = (e.points.first().copied(), e.points.last().copied()); + let (a, b) = (a.unwrap_or(e.mid), b.unwrap_or(e.mid)); + vec![a + (b - a) * 0.25, a + (b - a) * 0.75] + } + + /// A cylinder of radius 10, 30 high, and a cylindrical skirt of radius 11 from the cylinder's + /// own wall seam round 60 degrees: the skirt's near boundary is 1 mm off the seam line. + fn boss_and_skirt(s: &mut Session) { + run(s, "sketch.create", json!({"plane": "XY", "name": "Round"})); + run(s, "sketch.circle.center", json!({"center": [0, 0], "radius": 10})); + run(s, "sketch.finish", json!({})); + run(s, "solid.extrude", json!({"sketch": "Round", "distance": 30, "body_name": "Boss"})); + let seam = picks(s)[0]; + let a = seam.y.atan2(seam.x); + let at = |t: f64| json!([11.0 * t.cos(), 11.0 * t.sin()]); + run(s, "sketch.create", json!({"plane": "XY", "name": "Skirt sketch"})); + let arc = run(s, "sketch.arc.three_point", json!({"start": at(a), "end": at(a + 60f64.to_radians()), "through": at(a + 30f64.to_radians())})) + ["created"]["curves"][0] + .as_str() + .unwrap() + .to_string(); + run(s, "sketch.finish", json!({})); + run(s, "surface.extrude", json!({"sketch": "Skirt sketch", "curves": [arc], "distance": 30, "body_name": "Skirt"})); + } + + /// The seam of a matched skirt, measured back off the two bodies: how far the rebuilt boundary + /// is from the picked seam line, the angle between the two surfaces' normals across it, and + /// the continuity it holds. + fn seam(s: &Session) -> (f64, f64, Continuity) { + let st = s.world_state(); + let made = st.body("Skirt").unwrap().body.clone(); + let target = st.body("Boss").unwrap().body.clone(); + solidcraft_kernel::seam_continuity(&made, &target, &picks(s)).unwrap() + } + + /// `surface.match`: a skirt meets a cylinder's wall seam at g2; the body keeps its name, undo + /// and redo move it back and forth, a cold reopen gives the same seam, and what cannot be + /// matched is refused. + #[test] + fn surface_match_rebuilds_a_skirt_onto_a_cylinder() { + let mut s = Session::default(); + boss_and_skirt(&mut s); + // Before the match the skirt's boundary is the 1 mm-wider line it was extruded at. + let (gap, _, _) = seam(&s); + assert!((gap - 1.0).abs() < 0.05, "the skirt starts a millimetre off the seam: {gap}"); + let chain = picks(&s); + run(&mut s, "surface.match", json!({"body": "Skirt", "target": "Boss", "edges": chain, "order": "g2"})); + let st = s.world_state(); + assert!(st.body("Skirt").is_some(), "the matched body keeps its name"); + assert!(st.body("Boss").is_some(), "the body it was matched to stays too"); + assert!(st.body("Skirt").unwrap().body.is_surface(), "a match stays a surface"); + let (gap, deg, order) = seam(&s); + println!("engine: gap {gap:e} mm, normal {deg:e} deg, order {order:?}"); + assert!(gap < 1e-6, "the seam is on the picked line: {gap:e} mm"); + assert!(deg < 0.05, "the seam holds the wall's tangent plane: {deg} deg"); + assert!(order.conditions() >= Continuity::G1.conditions(), "{order:?}"); + // Undo puts the boundary back, redo matches it again, and a cold reopen agrees. + run(&mut s, "edit.undo", json!({})); + assert!(seam(&s).0 > 0.5, "undo restores the original skirt"); + run(&mut s, "edit.redo", json!({})); + assert!(seam(&s).1 < 0.05, "redo matches it again"); + let cold = Session::new(solidcraft_doc::Document::from_json(&s.doc.to_json()).unwrap()); + assert!(seam(&cold).1 < 0.05, "a reopened design matches the same"); + // Refusals: an unknown body, no target, no picks on the target, a word that is no order, + // and a body that is not a surface. + assert!(s.execute("surface.match", &json!({"body": "Nope", "target": "Boss", "edges": picks(&s)})).is_err()); + assert!(s.execute("surface.match", &json!({"body": "Skirt", "edges": picks(&s)})).is_err()); + assert!(s.execute("surface.match", &json!({"body": "Skirt", "target": "Boss", "edges": []})).is_err()); + assert!(s.execute("surface.match", &json!({"body": "Skirt", "target": "Nope", "edges": picks(&s)})).is_err()); + assert!(s.execute("surface.match", &json!({"body": "Skirt", "target": "Boss", "edges": picks(&s), "order": "g5"})).is_err()); + assert!(s.execute("surface.match", &json!({"body": "Skirt", "target": "Boss", "edges": [[100, 0, 0]]})).is_err()); + run(&mut s, "sketch.create", json!({"plane": "XY", "name": "B"})); + run(&mut s, "sketch.rectangle.two_point", json!({"p0": [0, 0], "p1": [4, 4]})); + run(&mut s, "sketch.finish", json!({})); + run(&mut s, "solid.extrude", json!({"sketch": "B", "distance": 4, "body_name": "Block"})); + assert!(s.execute("surface.match", &json!({"body": "Block", "target": "Boss", "edges": picks(&s)})).is_err(), "a solid is not a surface"); + } +} diff --git a/crates/kernel/src/boundary.rs b/crates/kernel/src/boundary.rs new file mode 100644 index 00000000..c0825f2e --- /dev/null +++ b/crates/kernel/src/boundary.rs @@ -0,0 +1,123 @@ +//! Boundary patches: one face through the boundary curves of a sketch, each side meeting the +//! surface it was taken from at a requested continuity (G0 up to G3). +//! +//! Two boundaries are the two sides of a [`crate::bridge`]: the patch is the Hermite surface that +//! leaves each boundary along the direction its frames state and repeats as much of that support +//! surface — its tangent plane, its curvature, the change of that curvature — as the side's order +//! asks for, with `crown` shaping the middle without touching a boundary condition. +//! +//! A patch through *three or four* boundaries is refused rather than built approximately. At a +//! corner two sides prescribe the same entries of the control net from two different surfaces, and +//! a tensor-product patch can hold both only where those two prescriptions already agree (one wall +//! carrying every side, say). Where they do not — the general case — the requested order would +//! silently drop near the corner, or the boundary curve would visibly leave its stations, and a +//! Gregory-type corner (a rational blend of the sides' own data) is needed to hold them instead +//! (issue #763). Everything this module refuses is an error, never a quietly worse surface. + +use solidcraft_geom::Vec3; + +use crate::body::Body; +use crate::bridge::{Continuity, Frame, bridge}; +use crate::{KernelError, Result}; + +/// The most stations a boundary of a patch may carry. +const MAX_SAMPLES: usize = 512; + +/// One boundary of a patch: the stations along it, and how much of the surface it was taken from +/// the patch repeats there. +#[derive(Clone, Debug)] +pub struct Boundary { + /// The stations along the boundary, in the order they run: the end of one boundary meets the + /// start of the next, and the last meets the first. + pub frames: Vec, + /// How much of the neighbouring surface the patch repeats at this boundary. + pub order: Continuity, +} + +/// A two-sided patch: one face that meets boundary `i` at the continuity its `order` asks for (G0 up +/// to G3), with `samples` stations every boundary must carry and `crown` shaping the middle of the +/// patch without touching any boundary condition (0 the plain blend, up to 3 fuller, down to -0.9 +/// flatter). Three or four boundaries are refused: see the module docs. +pub fn boundary_patch(boundaries: &[Boundary], crown: f64, samples: usize) -> Result { + check_input(boundaries, crown, samples)?; + if boundaries.len() > 2 { + return Err(KernelError::Invalid(format!( + "not supported yet: a patch through {} boundaries — a three- or four-sided patch needs a Gregory-type corner (issue #763); \ + give two boundaries, which is the bridge this builds exactly", + boundaries.len() + ))); + } + let (Some(a), Some(b)) = (boundaries.first(), boundaries.get(1)) else { + return Err(KernelError::Invalid("a boundary patch joins 2 boundaries: one of them is missing".into())); + }; + bridge(&a.frames, &b.frames, [a.order, b.order], [crown, crown]) +} + +/// Everything a patch refuses, before any geometry is built. +fn check_input(boundaries: &[Boundary], crown: f64, samples: usize) -> Result<()> { + let n = boundaries.len(); + if !(2..=4).contains(&n) { + return Err(KernelError::Invalid(format!("a boundary patch joins 2 to 4 boundaries, {n} given"))); + } + if !(2..=MAX_SAMPLES).contains(&samples) { + return Err(KernelError::Invalid(format!("a boundary patch needs 2 to {MAX_SAMPLES} stations per boundary, {samples} given"))); + } + if !crown.is_finite() || !(-0.9..=3.0).contains(&crown) { + return Err(KernelError::Invalid("a boundary patch's crown must be between -0.9 and 3".into())); + } + for (i, b) in boundaries.iter().enumerate() { + if b.frames.len() != samples { + return Err(KernelError::Invalid(format!( + "boundary {i} carries {} stations, not the {samples} every boundary of a patch needs", + b.frames.len() + ))); + } + for f in &b.frames { + if ![f.p, f.leave, f.bend, f.bend_change].into_iter().all(Vec3::is_finite) { + return Err(KernelError::Invalid(format!("boundary {i} has a station that is not finite"))); + } + if f.leave.normalized().is_none() { + return Err(KernelError::Invalid(format!("boundary {i} has a station with no leave direction"))); + } + } + if boundary_length(&b.frames) <= 1e-9 { + return Err(KernelError::Invalid(format!("boundary {i} has no length: a patch needs a curve, not a single point"))); + } + } + // Two boundaries are the two sides of a bridge, which do not meet; three or four would run round + // a loop, so each one would have to end where the next one begins. + if n >= 3 { + let mut scale: f64 = 1.0; + for b in boundaries { + scale = scale.max(boundary_length(&b.frames)); + } + let tol = scale * 1e-4 + 1e-9; + for i in 0..n { + let (Some(a), Some(b)) = (boundaries.get(i), boundaries.get((i + 1) % n)) else { continue }; + let (Some(end), Some(start)) = (a.frames.last(), b.frames.first()) else { continue }; + let gap = end.p.dist(start.p); + if gap > tol { + return Err(KernelError::Invalid(format!( + "boundary {i} ends and boundary {} starts {gap:.6} mm apart: the boundaries of a loop must run on into one another", + (i + 1) % n + ))); + } + } + } + Ok(()) +} + +/// The length of a run of stations. +fn boundary_length(frames: &[Frame]) -> f64 { + let mut length: f64 = 0.0; + for w in frames.windows(2) { + if let (Some(a), Some(b)) = (w.first(), w.get(1)) { + length += a.p.dist(b.p); + } + } + length +} + +#[cfg(test)] +#[path = "boundary_tests.rs"] +mod tests; diff --git a/crates/kernel/src/boundary_tests.rs b/crates/kernel/src/boundary_tests.rs new file mode 100644 index 00000000..ce7abfcb --- /dev/null +++ b/crates/kernel/src/boundary_tests.rs @@ -0,0 +1,276 @@ +//! Tests of boundary patches: the bridge this builds is measured back off the face it made — not +//! from the construction — and everything it cannot hold is an error rather than a quiet +//! approximation. + +use super::*; +use crate::body::from_p3; +use truck_modeling::ParametricSurface; + +/// The angle between two vectors, in degrees. +fn angle_deg(a: Vec3, b: Vec3) -> f64 { + (a.normalized().unwrap().dot(b.normalized().unwrap())).clamp(-1.0, 1.0).acos().to_degrees() +} + +/// The geometric shape of a surface curve from its first three derivatives along it: the unit +/// tangent, the curvature vector and the arclength derivative of that curvature vector — the same +/// reading `bridge` states its conditions with. +fn shape_of(d1: Vec3, d2: Vec3, d3: Vec3) -> (Vec3, Vec3, Vec3) { + let sigma = d1.len(); + let t = d1 * (1.0 / sigma); + let tp = d2 * (1.0 / sigma) - d1 * (d1.dot(d2) / (sigma * sigma * sigma)); + let nn = d2 - t * d2.dot(t); + let np = d3 - t * d3.dot(t) - tp * d2.dot(t) - t * d2.dot(tp); + let kappa = np * (sigma * sigma) - nn * (2.0 * d1.dot(d2)); + (t, nn * (1.0 / (sigma * sigma)), kappa * (1.0 / (sigma.powi(5)))) +} + +/// Across a boundary: five samples a step `h` apart, `into` the patch from the parameter `(u, v)`, +/// read back as the position, the unit tangent, the curvature vector and the change of curvature +/// along that inward direction. +fn across(body: &Body, uv: (f64, f64), into: (f64, f64), h: f64) -> (Vec3, Vec3, Vec3, Vec3) { + let face = body.solid.face_iter().next().unwrap(); + let surf = face.surface(); + let q: Vec = (0..5).map(|k| from_p3(surf.subs(uv.0 + k as f64 * h * into.0, uv.1 + k as f64 * h * into.1))).collect(); + let (f0, f1, f2, f3, f4) = (q[0], q[1], q[2], q[3], q[4]); + let d1 = (f0 * -25.0 + f1 * 48.0 + f2 * -36.0 + f3 * 16.0 + f4 * -3.0) * (1.0 / (12.0 * h)); + let d2 = (f0 * 35.0 + f1 * -104.0 + f2 * 114.0 + f3 * -56.0 + f4 * 11.0) * (1.0 / (12.0 * h * h)); + let d3 = (f0 * -5.0 + f1 * 18.0 + f2 * -24.0 + f3 * 14.0 + f4 * -3.0) * (1.0 / (2.0 * h * h * h)); + let (t, bend, change) = shape_of(d1, d2, d3); + (f0, t, bend, change) +} + +/// The stations' parameters: `n` of them, evenly spread over the patch. +fn stations(n: usize) -> Vec { + (0..n).map(|i| i as f64 / (n - 1) as f64).collect() +} + +/// The normal angle a side may be off by. +const ANGLE_TOL: f64 = 0.01; +/// How much of the supporting wall's own curvature a side must repeat, relative to it. +const CURV_TOL: f64 = 1e-3; +/// The normal part of the change of curvature a G3 side may leave over, relative to 1/r² (the +/// finite differences that read the third derivative off the face amplify round-off by 1/h³, so +/// this is a noise floor rather than a geometric tolerance). +const CHANGE_TOL: f64 = 1e-4; + +/// Measure one side of a built patch at every station of `frames`, and check what its order asks +/// for: the boundary curve itself always, the tangent plane from G1 on, the curvature vector from +/// G2 on and its change along the leave direction at G3. +fn check_side(body: &Body, frames: &[Frame], order: Continuity, at: impl Fn(usize) -> ((f64, f64), (f64, f64)), r: f64) { + for (i, f) in frames.iter().enumerate() { + let (uv, into) = at(i); + let (p, t, bend, change) = across(body, uv, into, 1e-3); + let gap = p.dist(f.p); + assert!(gap < 1e-9, "station {i} sits on its boundary: {p:?} vs {:?} ({gap:.2e} mm)", f.p); + if order == Continuity::G0 { + continue; + } + assert!(angle_deg(t, f.leave) < ANGLE_TOL, "station {i} leaves along the tangent: {t:?} vs {:?}", f.leave); + if order == Continuity::G1 { + continue; + } + if f.bend.len() > 1e-9 { + assert!(bend.dist(f.bend) / f.bend.len() < CURV_TOL, "station {i}: curvature {bend:?} vs the wall's {:?}", f.bend); + if order == Continuity::G3 { + let normal = f.bend * (1.0 / f.bend.len()); + assert!(change.dot(normal).abs() < CHANGE_TOL / (r * r), "station {i}: the change of curvature has a normal part {change:?}"); + } + } else { + assert!(bend.len() < CURV_TOL / r, "station {i}: the wall is straight here, the patch curves by {bend:?}"); + if order == Continuity::G3 { + assert!(change.len() < CHANGE_TOL / (r * r), "station {i}: a straight cross-section changes curvature by {change:?}"); + } + } + } +} + +/// Two walls 20 mm apart, each leaving towards the other, both walls of the same cylinder of radius +/// `r` (the valley `bridge` is tested with): the two sides of a two-sided patch. +fn two_walls(n: usize, r: f64) -> (Vec, Vec) { + let mut a = Vec::new(); + let mut b = Vec::new(); + for t in stations(n) { + let y = 20.0 * t; + a.push(Frame { + p: Vec3::new(0.0, y, 0.0), + leave: Vec3::new(0.0, 1.0, 0.0), + bend: Vec3::new(0.0, 0.0, 1.0 / r), + bend_change: Vec3::new(0.0, -1.0 / (r * r), 0.0), + }); + b.push(Frame { + p: Vec3::new(0.0, y + 20.0, 0.0), + leave: Vec3::new(0.0, -1.0, 0.0), + bend: Vec3::new(0.0, 0.0, 1.0 / r), + bend_change: Vec3::new(0.0, 1.0 / (r * r), 0.0), + }); + } + (a, b) +} + +#[test] +fn a_two_sided_patch_is_a_bridge_that_holds_both_ends() { + let n = 9; + let r = 10.0; + let (a, b) = two_walls(n, r); + let boundaries = vec![Boundary { frames: a.clone(), order: Continuity::G2 }, Boundary { frames: b.clone(), order: Continuity::G2 }]; + let body = boundary_patch(&boundaries, 0.0, n).unwrap(); + let s = stations(n); + check_side(&body, &a, Continuity::G2, |i| ((s[i], 0.0), (0.0, 1.0)), r); + check_side(&body, &b, Continuity::G2, |i| ((s[i], 1.0), (0.0, -1.0)), r); +} + +/// Four chained sides round a 20 × 20 square in the XY plane, each leaving into the square: the +/// loop a patch through three or four boundaries would have to close. +fn square_loop(n: usize) -> Vec { + let corners = [(0.0_f64, 0.0_f64), (20.0, 0.0), (20.0, 20.0), (0.0, 20.0)]; + let leaves = [Vec3::new(0.0, 1.0, 0.0), Vec3::new(-1.0, 0.0, 0.0), Vec3::new(0.0, -1.0, 0.0), Vec3::new(1.0, 0.0, 0.0)]; + (0..4) + .map(|side| { + let (a, b) = (corners[side], corners[(side + 1) % 4]); + let frames = stations(n) + .into_iter() + .map(|t| Frame { + p: Vec3::new(a.0 + (b.0 - a.0) * t, a.1 + (b.1 - a.1) * t, 0.0), + leave: leaves[side], + bend: Vec3::ZERO, + bend_change: Vec3::ZERO, + }) + .collect(); + Boundary { frames, order: Continuity::G2 } + }) + .collect() +} + +#[test] +fn a_two_sided_patch_meets_its_walls_up_to_the_curvature() { + let n = 11; + let r = 10.0; + let (a, b) = two_walls(n, r); + // G0, G1 and G2 are what this patch states: position, tangent plane and curvature across each + // boundary. (A G3 bridge's third-order rows are `bridge`'s own subject, measured in its tests.) + for order in [Continuity::G0, Continuity::G1, Continuity::G2] { + let boundaries = vec![Boundary { frames: a.clone(), order }, Boundary { frames: b.clone(), order }]; + let body = boundary_patch(&boundaries, 0.0, n).unwrap(); + let s = stations(n); + if order == Continuity::G0 { + for (i, f) in a.iter().enumerate() { + let (p, _, _, _) = across(&body, (s[i], 0.0), (0.0, 1.0), 1e-3); + assert!(p.dist(f.p) < 1e-9, "station {i}: {p:?} vs {:?}", f.p); + } + continue; + } + check_side(&body, &a, order, |i| ((s[i], 0.0), (0.0, 1.0)), r); + check_side(&body, &b, order, |i| ((s[i], 1.0), (0.0, -1.0)), r); + } +} + +#[test] +fn a_crown_shapes_the_middle_without_touching_the_boundary_rows() { + let n = 9; + let r = 10.0; + let (a, b) = two_walls(n, r); + let boundaries = vec![Boundary { frames: a.clone(), order: Continuity::G2 }, Boundary { frames: b.clone(), order: Continuity::G2 }]; + let s = stations(n); + let mut middles = Vec::new(); + for crown in [-0.5, 0.0, 1.0, 2.0] { + let body = boundary_patch(&boundaries, crown, n).unwrap(); + let grid = crate::bridge::sample_grid(&body, 9, 9).unwrap(); + middles.push(grid[4][4].z); + // The interior is the cross-section between the two walls: its middle sits in the gap, not + // outside it, and it stays finite everywhere on the grid. + for row in &grid { + for p in row { + assert!(p.is_finite() && (-2.0..=42.0).contains(&p.y), "the interior stays between the walls: {p:?}"); + } + } + assert!((18.0..=22.0).contains(&grid[4][4].y), "the middle of the patch sits midway in the 20 mm band: {:?}", grid[4][4]); + // The conditions are untouched: the same boundary points, tangent and curvature. + check_side(&body, &a, Continuity::G2, |i| ((s[i], 0.0), (0.0, 1.0)), r); + check_side(&body, &b, Continuity::G2, |i| ((s[i], 1.0), (0.0, -1.0)), r); + } + let steps: Vec = middles.windows(2).map(|w| w[1] - w[0]).collect(); + assert!(steps.iter().all(|d| d.abs() > 1e-9), "a fuller crown moves the middle: {middles:?}"); +} + +#[test] +fn the_patch_refuses_what_it_cannot_hold() { + let n = 9; + let (a, b) = two_walls(n, 10.0); + let two = vec![Boundary { frames: a.clone(), order: Continuity::G2 }, Boundary { frames: b.clone(), order: Continuity::G2 }]; + // Not 2 to 4 boundaries. + assert!(boundary_patch(&[], 0.0, n).is_err()); + assert!(boundary_patch(&two[..1], 0.0, n).is_err()); + let five = vec![two[0].clone(), two[1].clone(), two[0].clone(), two[1].clone(), two[0].clone()]; + assert!(boundary_patch(&five, 0.0, n).is_err()); + // Three or four sides are refused: a tensor-product patch cannot hold two sides' data at a + // corner, which needs a Gregory-type corner (issue #763). The message says so. + let square = square_loop(n); + let e = boundary_patch(&square, 0.0, n).unwrap_err().to_string(); + assert!(e.contains("not supported yet") && e.contains("#763"), "4 sides: {e}"); + // Three sides that do close: a triangle, whose last side runs from the second one's end back to + // the first one's start. + let triangle: Vec = vec![ + square[0].clone(), + square[1].clone(), + Boundary { + frames: stations(n) + .into_iter() + .map(|t| Frame { + p: Vec3::new(20.0 * (1.0 - t), 20.0 * (1.0 - t), 0.0), + leave: Vec3::new(-std::f64::consts::FRAC_1_SQRT_2, -std::f64::consts::FRAC_1_SQRT_2, 0.0), + bend: Vec3::ZERO, + bend_change: Vec3::ZERO, + }) + .collect(), + order: Continuity::G2, + }, + ]; + let e = boundary_patch(&triangle, 0.0, n).unwrap_err().to_string(); + assert!(e.contains("not supported yet") && e.contains("#763"), "3 sides: {e}"); + // A loop whose sides do not even meet is refused earlier and for its own reason. + let broken = square_loop(n); + let mut moved = broken[1].frames.clone(); + moved[0] = Frame { p: Vec3::new(30.0, 0.0, 0.0), ..moved[0] }; + let broken = vec![broken[0].clone(), Boundary { frames: moved, order: Continuity::G2 }, broken[2].clone()]; + let e = boundary_patch(&broken, 0.0, n).unwrap_err().to_string(); + assert!(e.contains("run on into one another"), "{e}"); + // The stations and the crown. + let short = vec![Boundary { frames: a[..n - 1].to_vec(), order: Continuity::G2 }, two[1].clone()]; + assert!(boundary_patch(&short, 0.0, n).is_err(), "a boundary with the wrong number of stations"); + assert!(boundary_patch(&two, 0.0, n - 1).is_err()); + assert!(boundary_patch(&two, 0.0, 1).is_err()); + for crown in [-1.5, 3.5, f64::NAN, f64::INFINITY] { + assert!(boundary_patch(&two, crown, n).is_err(), "crown {crown}"); + } + // A station that is not finite, and a station with no leave direction. + for field in 0..4 { + let mut bad = two.clone(); + let mut frames = bad[0].frames.clone(); + frames[3] = match field { + 0 => Frame { p: Vec3::new(f64::NAN, 0.0, 0.0), ..frames[3] }, + 1 => Frame { leave: Vec3::new(f64::NAN, 0.0, 0.0), ..frames[3] }, + 2 => Frame { bend: Vec3::new(f64::INFINITY, 0.0, 0.0), ..frames[3] }, + _ => Frame { bend_change: Vec3::new(f64::NAN, 0.0, 0.0), ..frames[3] }, + }; + if let Some(b) = bad.first_mut() { + b.frames = frames; + } + assert!(boundary_patch(&bad, 0.0, n).is_err(), "a station that is not finite: field {field}"); + } + let mut flat = two.clone(); + if let Some(b) = flat.first_mut() { + b.frames = vec![Frame { p: Vec3::new(0.0, 1.0, 0.0), leave: Vec3::ZERO, bend: Vec3::ZERO, bend_change: Vec3::ZERO }; n]; + } + assert!(boundary_patch(&flat, 0.0, n).is_err(), "a station with no leave direction"); + let point = vec![ + Boundary { + frames: vec![Frame { p: Vec3::new(1.0, 1.0, 1.0), leave: Vec3::X, bend: Vec3::ZERO, bend_change: Vec3::ZERO }; n], + order: Continuity::G1, + }, + Boundary { + frames: vec![Frame { p: Vec3::new(1.0, 1.0, 1.0), leave: Vec3::X, bend: Vec3::ZERO, bend_change: Vec3::ZERO }; n], + order: Continuity::G1, + }, + ]; + assert!(boundary_patch(&point, 0.0, n).is_err(), "a boundary with no length"); +} diff --git a/crates/kernel/src/fair.rs b/crates/kernel/src/fair.rs new file mode 100644 index 00000000..d8b4a0bb --- /dev/null +++ b/crates/kernel/src/fair.rs @@ -0,0 +1,1538 @@ +//! Fairing the control points of a B-spline curve or surface: moving the interior control +//! vertices to the least-bending shape while the ends, or the boundary rows, stay exactly where +//! the design put them. +//! +//! The fairness energy is the thin-plate (curvature-variation) energy of the *geometry* the +//! control net describes, in the net's own parameter directions: +//! +//! - a curve: `E = ∫‖C″(t)‖² dt`, +//! - a surface: `E = ∫∫ ‖S_uu‖² + 2‖S_uv‖² + ‖S_vv‖² du dv`. +//! +//! A plane has energy zero and anything else pays for how much it bends, which is what a +//! designer means by "fair". Both energies are quadratic in the control points, `E(P) = Pᵀ·G·P`, +//! with `G` the Gram matrix of the basis functions' derivatives (`Gᵢⱼ = ∫Nᵢ″Nⱼ″ dt`, and the +//! tensor products of the one-direction Gram matrices for a surface). Each one-direction Gram is +//! integrated span by span with four-point Gauss–Legendre, exact to degree 7 — a cubic's second +//! derivative products are quadratic, so the cubic curves and bicubic faces fairing meets in +//! practice are integrated exactly. +//! +//! One fairing step is then a constrained least-squares solve: hold the boundary control points, +//! take the exact minimiser of `E` over the interior (for a quadratic that is one Newton step, +//! `G_ff·d = −(G·P)_f`), and apply a fraction `strength` of it, so `strength = 1` is the full +//! solution and a smaller value walks towards it. With the exact step the energy falls by +//! `q·s·(1 − s)` for `q ≥ 0`, so a run at any strength decreases it strictly until it converges. +//! `G_ff` does not change as the control points move, so it is factorised once, by banded +//! Cholesky, and every further step is a matrix–vector product and a triangular solve. +//! +//! Two guarantees are built in rather than checked afterwards: +//! +//! - **The boundary never moves.** The result starts as a copy of the input and only interior +//! entries are written, so a held control point is bit-identical to the one the design gave. +//! - **The geometry cannot jump.** A B-spline lies in the convex hull of its control points +//! (NURBS Book §3.1; with the weights unchanged the rational basis functions are still +//! non-negative and sum to one), so the geometry moves at most as far as the farthest control +//! point. Every step is clamped so that no control point leaves the ball of radius +//! `tolerance` about where it started: the faired curve or surface is then within `tolerance` +//! of the input at the same parameters, by a bound rather than by measurement. A step that +//! would not lower the energy is halved and, if it still would not, the run stops — it never +//! loops for ever and never moves the design further than the caller allowed. +//! +//! Everything is checked before it is used: degrees, knot vectors and counts must be consistent, +//! the control points finite, `strength ∈ (0, 1]`, `iterations` in `1…`[`MAX_FAIR_ITERATIONS`] and +//! `tolerance > 0`. A singular system — fewer than two interior control points, or an interior +//! control point that does not change the geometry at all — is reported as an error and the +//! geometry is left alone. + +use solidcraft_geom::Vec3; +use solidcraft_geom::bspline::{HPoint, MAX_SPLINE_DEGREE, NurbsSurface, Side, SplineError, check_knots, local_basis}; +use truck_modeling as mt; + +use crate::body::Body; +use crate::history::Change; +use crate::{KernelError, Result, guard}; + +/// Most steps one fairing run takes (see [`FairOptions::iterations`]). +pub const MAX_FAIR_ITERATIONS: usize = 1000; +/// Most control points one fairing run accepts (a 64 × 64 net): plenty for any face's control +/// net, and small enough that a step stays cheap whatever a design file contains. +pub const MAX_FAIR_POINTS: usize = 4096; +/// Most entries of the banded system (2²¹ doubles, 16 MB). +const MAX_FAIR_BAND: usize = 1 << 21; +/// The fairness energy of a curve, as a one-direction net: the `(order in u, order in v, weight)` +/// term list of `∫∫ (∂S)²`. +const CURVE_ENERGY: [(usize, usize, f64); 1] = [(2, 0, 1.0)]; +/// The thin-plate energy of a surface: `‖S_uu‖² + 2‖S_uv‖² + ‖S_vv‖²` — the classic +/// curvature-variation energy, zero on a plane and the same if u and v are swapped. +const SURFACE_ENERGY: [(usize, usize, f64); 3] = [(2, 0, 1.0), (1, 1, 2.0), (0, 2, 1.0)]; +/// Four-point Gauss–Legendre nodes and weights on `[−1, 1]` (exact to degree 7). +const GAUSS4: [(f64, f64); 4] = [ + (-0.861_136_311_594_052_6, 0.347_854_845_137_453_8), + (-0.339_981_043_584_856_3, 0.652_145_154_862_546_2), + (0.339_981_043_584_856_3, 0.652_145_154_862_546_2), + (0.861_136_311_594_052_6, 0.347_854_845_137_453_8), +]; +/// The second direction of a curve seen as a control net: one basis function of degree 0 over +/// `[0, 1]`, whose Gram entry is exactly `∫1 dt = 1`, so the net's energy is the curve's own +/// energy unchanged. +const TRIVIAL_KNOTS: [f64; 2] = [0.0, 1.0]; +/// How often a step is halved before the run calls the point converged. +const MAX_HALVINGS: u32 = 40; +/// A step smaller than this fraction of the full Newton step is not worth taking. +const MIN_STEP: f64 = 1e-12; +/// Steps are clamped this far inside the tolerance, so rounding cannot push a control point +/// past it. +const TOLERANCE_SLACK: f64 = 1.0 - 1e-12; + +/// How a fairing run is asked to behave. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct FairOptions { + /// The fraction of each energy step to apply, `0 < strength ≤ 1`: `1` is the full least + /// squares solution (the exact minimiser of the energy at the held boundary), a smaller + /// value takes that fraction of the way there. The energy falls at every step either way. + pub strength: f64, + /// How many steps at most, one to [`MAX_FAIR_ITERATIONS`]. + pub iterations: usize, + /// Hold the layer of control points next to the boundary as well as the boundary itself (two + /// layers held at each end). For a clamped knot vector the boundary row fixes the position + /// along the edge and the row next to it fixes the derivative across the edge, so with this + /// set the boundary keeps the position and the tangent the design gave it. + pub keep_edges: bool, + /// The farthest a control point may move from where the design put it, in mm. The geometry + /// then stays within this distance of the input at the same parameters (see the module + /// docs), so a fairing run can never make a design jump. + pub tolerance: f64, +} + +impl FairOptions { + /// Checked options for one fairing run: an error, never a panic, for a strength outside + /// `(0, 1]`, zero iterations, more than [`MAX_FAIR_ITERATIONS`] of them, or a tolerance that + /// is not a positive finite number. + pub fn new(strength: f64, iterations: usize, keep_edges: bool, tolerance: f64) -> Result { + FairOptions { strength, iterations, keep_edges, tolerance }.checked() + } + + /// The same checks as [`FairOptions::new`], for options built field by field. + pub fn checked(self) -> Result { + if !(self.strength.is_finite() && self.strength > 0.0 && self.strength <= 1.0) { + return Err(KernelError::Invalid(format!("the fairing strength {} must be greater than 0 and at most 1", self.strength))); + } + if self.iterations == 0 || self.iterations > MAX_FAIR_ITERATIONS { + return Err(KernelError::Invalid(format!("{} fairing iterations (1 to {MAX_FAIR_ITERATIONS})", self.iterations))); + } + if !(self.tolerance.is_finite() && self.tolerance > 0.0) { + return Err(KernelError::Invalid(format!("the fairing tolerance {} must be a positive number of mm", self.tolerance))); + } + Ok(self) + } +} + +impl Default for FairOptions { + /// One full step with the boundary held, at most 1 mm of movement: the defaults a command can + /// start from before a user picks their own. + fn default() -> Self { + FairOptions { strength: 1.0, iterations: 1, keep_edges: false, tolerance: 1.0 } + } +} + +/// What a fairing run did: the energy before and after, how many steps were applied, how far the +/// control points moved and the energy after every applied step. The energy is a sum of squares, +/// so where rounding leaves it a hair below zero the true value is zero: it is reported as zero. +#[derive(Clone, Debug, PartialEq)] +pub struct FairReport { + /// The fairness energy of the input, in the parameter units of the net's own knots. + pub energy_before: f64, + /// The fairness energy of the result: never more than `energy_before`, and equal to it only + /// when no step could lower it. + pub energy_after: f64, + /// How many steps were applied (at most the requested `iterations`). + pub steps: usize, + /// The farthest any control point moved from where the design put it, in mm. It bounds how + /// far the geometry moved from the input (see the module docs) and never exceeds the + /// requested tolerance. + pub movement: f64, + /// The energy after each applied step, the last equal to `energy_after`: it never rises, and + /// falls strictly at every step until the run reaches the rounding of its own arithmetic. + pub energy_trace: Vec, +} + +/// A faired B-spline curve: the input's degree and knots, the moved control points and the run's +/// report. +#[derive(Clone, Debug, PartialEq)] +pub struct FairCurve { + pub degree: usize, + pub knots: Vec, + /// The faired control points; the ends (and, with `keep_edges`, the points next to them) are + /// bit-identical to the input. + pub ctrl: Vec<[f64; N]>, + pub report: FairReport, +} + +/// A faired B-spline control net: the input's degrees, knots and counts, the moved control points +/// (u fastest) and the run's report. +#[derive(Clone, Debug, PartialEq)] +pub struct FairNet { + pub degree: [usize; 2], + pub knots: [Vec; 2], + pub count: [usize; 2], + /// The faired control points, u fastest; the held boundary layers are bit-identical to the + /// input. + pub ctrl: Vec<[f64; N]>, + pub report: FairReport, +} + +/// A geom spline error as the kernel's. +fn spline_error(what: &str, e: SplineError) -> KernelError { + match e { + SplineError::TooLarge(m) => KernelError::TooLarge(format!("{what}: {m}")), + other => KernelError::Invalid(format!("{what}: {other}")), + } +} + +/// The dot product of two control points' coordinate vectors. +fn dot(a: &[f64; N], b: &[f64; N]) -> f64 { + a.iter().zip(b).map(|(x, y)| x * y).sum() +} + +/// The distance between two control points. +fn dist(a: &[f64; N], b: &[f64; N]) -> f64 { + let d: [f64; N] = std::array::from_fn(|c| a[c] - b[c]); + dot(&d, &d).sqrt() +} + +/// The Gram matrices of one parameter direction, stored banded: entry `r` is +/// `∫ Nᵢ^{(r)}(t)·Nₖ^{(r)}(t) dt` over the direction's domain, and since a basis function +/// overlaps only its `degree` neighbours an entry is non-zero only for `|i − k| ≤ degree`, so +/// the whole thing is `count·(degree + 1)` numbers per order. +struct Gram { + count: usize, + bw: usize, + g: [Vec; 3], +} + +impl Gram { + /// The Gram matrices of `count` basis functions of `degree` over `knots`, for orders 0, 1 and + /// 2 — the derivatives a thin-plate energy is built from. + fn new(degree: usize, knots: &[f64], count: usize) -> Result { + let bw = degree; + let mut g: [Vec; 3] = std::array::from_fn(|_| vec![0.0; count * (bw + 1)]); + // On the span `[knots[s], knots[s + 1])` only the basis functions `i ∈ [s − degree, s]` + // are non-zero. The quadrature nodes lie strictly inside the span, so either one-sided + // limit gives the same value there. + let mut ders = [[0.0f64; 3]; MAX_SPLINE_DEGREE + 1]; + for s in degree..count { + let (Some(&a), Some(&b)) = (knots.get(s), knots.get(s + 1)) else { continue }; + if !(b > a) { + continue; + } + for &(x, w) in GAUSS4.iter() { + let t = 0.5 * (a + b) + 0.5 * (b - a) * x; + let first = s - degree; + let mut live = 0usize; + for i in first..=s { + let Some(local) = knots.get(i..=i + degree + 1) else { break }; + let d = local_basis(local, t, Side::Right, 2); + if let Some(slot) = ders.get_mut(live) { + for (o, v) in slot.iter_mut().zip(d.iter()) { + *o = *v; + } + } + live += 1; + } + let weight = 0.5 * (b - a) * w; + let value = |rel: usize, r: usize| ders.get(rel).and_then(|d| d.get(r)).copied().unwrap_or(0.0); + // The upper triangle only — the band holds one copy of a symmetric entry — and at + // the basis functions' own indices, not their places in this span. + for i in 0..live { + for k in i..live { + let (row, off) = (first + i, k - i); + for (r, gg) in g.iter_mut().enumerate() { + if let Some(slot) = gg.get_mut(row * (bw + 1) + off) { + *slot += weight * value(i, r) * value(k, r); + } + } + } + } + } + } + if g.iter().flatten().any(|v| !v.is_finite()) { + return Err(KernelError::Invalid("the fairness energy of this net is not a finite number".into())); + } + Ok(Gram { count, bw, g }) + } + + /// `∫ Nᵢ^{(r)} Nₖ^{(r)} dt`, zero outside the band or off the end of the basis. + fn at(&self, r: usize, i: usize, k: usize) -> f64 { + let Some(m) = self.g.get(r) else { return 0.0 }; + let (row, off) = if i <= k { (i, k - i) } else { (k, i - k) }; + if off > self.bw || row >= self.count { + return 0.0; + } + m.get(row * (self.bw + 1) + off).copied().unwrap_or(0.0) + } +} + +/// A control net to fair: a tensor-product B-spline with the terms of its fairness energy. A +/// direction of one control point (degree 0, knots `[0, 1]`) makes the net a curve, whose energy +/// is the one-direction [`CURVE_ENERGY`]. +struct Net { + degree: [usize; 2], + count: [usize; 2], + /// The number of control points, `count[0] · count[1]`. + total: usize, + gram: [Gram; 2], + terms: &'static [(usize, usize, f64)], +} + +impl Net { + /// Check the degrees, knots and counts and build the Gram matrices: an error, never a panic, + /// for a degree above [`MAX_SPLINE_DEGREE`], too few control points for the degree, a knot + /// vector of the wrong length, unusable knots or an empty parameter range. + fn new(degree: [usize; 2], knots: [Vec; 2], count: [usize; 2], terms: &'static [(usize, usize, f64)]) -> Result { + let Some(total) = count[0].checked_mul(count[1]).filter(|&n| n <= MAX_FAIR_POINTS) else { + return Err(KernelError::TooLarge(format!("fairing {} × {} control points (at most {MAX_FAIR_POINTS} in all)", count[0], count[1]))); + }; + for d in 0..2 { + let (p, c, kn) = (degree[d], count[d], &knots[d]); + if p > MAX_SPLINE_DEGREE { + return Err(KernelError::TooLarge(format!("fairing a spline of degree {p} (at most {MAX_SPLINE_DEGREE})"))); + } + if c <= p { + return Err(KernelError::Invalid(format!("{c} control points is too few across one direction of a spline of degree {p}"))); + } + if kn.len() != c + p + 1 { + return Err(KernelError::Invalid(format!( + "{} knots across one direction, expected {} for {c} control points of degree {p}", + kn.len(), + c + p + 1 + ))); + } + check_knots(kn).map_err(|e| spline_error("fairing", e))?; + let (lo, hi) = (kn.get(p).copied().unwrap_or(f64::NAN), kn.get(c).copied().unwrap_or(f64::NAN)); + if !(lo < hi) { + return Err(KernelError::Invalid("the fairing's parameter range is empty".into())); + } + } + if terms.iter().any(|&(ru, rv, w)| ru > degree[0] || rv > degree[1] || !(w.is_finite() && w >= 0.0)) { + return Err(KernelError::Invalid("a fairness term does not fit the net's degrees".into())); + } + let gram = [Gram::new(degree[0], &knots[0], count[0])?, Gram::new(degree[1], &knots[1], count[1])?]; + Ok(Net { degree, count, total, gram, terms }) + } + + /// The index of control point `(i, j)` in the flattened net (u fastest). + fn index(&self, i: usize, j: usize) -> usize { + j * self.count[0] + i + } + + /// The fairness energy `Pᵀ·G·P` of the net `ctrl` (u fastest) and, in `grad`, `G·P`. One + /// banded sweep gives both, since `E = Σ Pₚ·(G·P)ₚ`. + fn energy_gradient(&self, ctrl: &[[f64; N]], grad: &mut [[f64; N]]) -> f64 { + let [cu, cv] = self.count; + let (pu, pv) = (self.degree[0], self.degree[1]); + let mut energy = 0.0; + for j in 0..cv { + for i in 0..cu { + let mut acc = [0.0f64; N]; + for &(ru, rv, w) in self.terms { + for l in j.saturating_sub(pv)..=(j + pv).min(cv - 1) { + for k in i.saturating_sub(pu)..=(i + pu).min(cu - 1) { + let a = w * self.gram[0].at(ru, i, k) * self.gram[1].at(rv, j, l); + if a == 0.0 { + continue; + } + if let Some(p) = ctrl.get(self.index(k, l)) { + for (x, v) in acc.iter_mut().zip(p) { + *x += a * v; + } + } + } + } + } + let at = self.index(i, j); + if let Some(c) = ctrl.get(at) { + energy += dot(c, &acc); + } + if let Some(g) = grad.get_mut(at) { + *g = acc; + } + } + } + energy + } + + /// The energy matrix `G` over the interior control points, banded and already factorised (the + /// matrix does not change as the control points move, so one factorisation serves every + /// step). In the u-fastest order two unknowns couple within `degree_u + degree_v · cols`: + /// `degree_u` along a row, `degree_v` rows apart. + fn system(&self, interior: &Interior) -> Result { + let [cu, cv] = self.count; + let (pu, pv) = (self.degree[0], self.degree[1]); + let rows = if cv > 1 { interior.rows.saturating_sub(1).min(pv) } else { 0 }; + let mut band = Band::new(interior.points.len(), pu + rows * interior.cols)?; + for (r, &at) in interior.points.iter().enumerate() { + let (i, j) = (at % cu, at / cu); + for &(ru, rv, w) in self.terms { + for l in j.saturating_sub(pv)..=(j + pv).min(cv - 1) { + for k in i.saturating_sub(pu)..=(i + pu).min(cu - 1) { + let Some(s) = interior.slot.get(self.index(k, l)).copied().flatten() else { continue }; + band.add(r, s, w * self.gram[0].at(ru, i, k) * self.gram[1].at(rv, j, l)); + } + } + } + } + band.factorised() + } + + /// Fair `ctrl` (u fastest) with checked `opts`: see the module docs. The result is the input + /// with the interior control points moved; the held ones are bit-identical. + fn fair(&self, ctrl: &[[f64; N]], opts: &FairOptions) -> Result<(Vec<[f64; N]>, FairReport)> { + if ctrl.len() != self.total { + return Err(KernelError::Invalid(format!("{} control points, expected {}", ctrl.len(), self.total))); + } + if ctrl.iter().flatten().any(|x| !x.is_finite()) { + return Err(KernelError::Invalid("a control point is not a finite number".into())); + } + let interior = Interior::new(self, if opts.keep_edges { 2 } else { 1 })?; + let system = self.system(&interior)?; + let mut current: Vec<[f64; N]> = ctrl.to_vec(); + let mut grad = vec![[0.0; N]; self.total]; + let energy_before = self.energy_gradient(¤t, &mut grad); + if !energy_before.is_finite() { + return Err(KernelError::Invalid("the fairness energy of this net is not a finite number".into())); + } + // The energy is a sum of squares, so where rounding leaves it a hair below zero the true + // value is zero: report it that way rather than as a negative number. + let energy_before = energy_before.max(0.0); + // How big the net itself is, for telling a fairing step from rounding (see below). + let size = interior.size(ctrl); + let mut step = vec![[0.0; N]; interior.points.len()]; + let mut b = vec![0.0; interior.points.len()]; + let mut report = FairReport { energy_before, energy_after: energy_before, steps: 0, movement: 0.0, energy_trace: Vec::new() }; + for _ in 0..opts.iterations { + let e_now = self.energy_gradient(¤t, &mut grad); + if !e_now.is_finite() { + break; + } + // The exact minimiser over the interior is the Newton step `G_ff·d = −(G·P)_f`; the + // system is per coordinate, all three sharing the one factorisation. + for c in 0..N { + for (slot, &at) in interior.points.iter().enumerate() { + b[slot] = -grad.get(at).and_then(|g| g.get(c)).copied().unwrap_or(0.0); + } + system.solve(&mut b)?; + for (slot, v) in b.iter().enumerate() { + if let Some(d) = step.get_mut(slot).and_then(|d| d.get_mut(c)) { + *d = *v; + } + } + } + if step.iter().flatten().all(|&x| x == 0.0) { + break; + } + // The largest fraction of the step that keeps every control point within the + // tolerance of where the design put it, and never more than the whole step. + let mut s = opts.strength.min(interior.limit(ctrl, ¤t, &step, opts.tolerance)); + if !(s > 0.0) { + break; + } + // A step this much shorter than the net's own size cannot be told from the rounding + // of the arithmetic that found it: the net is as fair as it gets. + let longest = step.iter().flatten().fold(0.0f64, |m, x| m.max(x.abs())); + if s * longest <= MIN_STEP * size { + break; + } + let mut taken = false; + for _ in 0..MAX_HALVINGS { + let candidate = interior.apply(¤t, &step, s); + let e = self.energy_gradient(&candidate, &mut grad); + if e < e_now { + current = candidate; + report.steps += 1; + report.energy_trace.push(e.max(0.0)); + taken = true; + break; + } + // The step did not lower the energy: halve it and try again, and when even a hair + // of a step will not, this is as fair as the net gets. + s *= 0.5; + if s <= MIN_STEP { + break; + } + } + if !taken { + break; + } + } + report.energy_after = report.energy_trace.last().copied().unwrap_or(energy_before); + report.movement = interior.movement(ctrl, ¤t); + Ok((current, report)) + } +} + +/// Which control points a fairing run may move, and where each sits in the banded system. +struct Interior { + /// The flattened index (u fastest) of interior point number `k`. + points: Vec, + /// The interior number of each flattened index, `None` where the point is held. + slot: Vec>, + /// How many control points across u (`cols`) and across v (`rows`) are interior. + cols: usize, + rows: usize, +} + +impl Interior { + /// The interior of `net` with `keep` layers held at each end of each direction (1 without + /// [`FairOptions::keep_edges`], 2 with it). + fn new(net: &Net, keep: usize) -> Result { + let [cu, cv] = net.count; + let (i0, i1) = free_span(cu, keep)?; + let (j0, j1) = free_span(cv, keep)?; + let (cols, rows) = (i1 - i0, j1 - j0); + let mut points = Vec::with_capacity(cols * rows); + let mut slot = vec![None; net.total]; + for j in j0..j1 { + for i in i0..i1 { + if let Some(s) = slot.get_mut(net.index(i, j)) { + *s = Some(points.len()); + } + points.push(net.index(i, j)); + } + } + if points.len() < 2 { + return Err(KernelError::Invalid(format!( + "fairing needs at least two interior control points; this {cu} × {cv} net has {}", + points.len() + ))); + } + Ok(Interior { points, slot, cols, rows }) + } + + /// `current` with `s` times the step `d` added to the interior control points; the held ones + /// are copied as they are, bit for bit. + fn apply(&self, current: &[[f64; N]], d: &[[f64; N]], s: f64) -> Vec<[f64; N]> { + let mut out = current.to_vec(); + for (slot, &at) in self.points.iter().enumerate() { + let (Some(p), Some(v)) = (out.get_mut(at), d.get(slot)) else { continue }; + for (x, dx) in p.iter_mut().zip(v) { + *x += s * dx; + } + } + out + } + + /// How big the interior is: the farthest interior control point from the first one, in mm. + /// A fairing step far below this belongs to the rounding of the arithmetic, not the design. + fn size(&self, ctrl: &[[f64; N]]) -> f64 { + let Some(first) = self.points.first().and_then(|&at| ctrl.get(at)) else { return 0.0 }; + self.points.iter().fold(0.0f64, |m, &at| ctrl.get(at).map_or(m, |p| m.max(dist(first, p)))) + } + + /// The largest fraction `s ≤ 1` of the step `d` that keeps every interior control point + /// within `tolerance` of where the design put it (`input`). The movement is `‖Δ + s·d‖ ≤ tol` + /// per point — a quadratic in `s` opening upwards, with the accumulated movement + /// `Δ = current − input` already inside the ball, so the bound is its positive root. + fn limit(&self, input: &[[f64; N]], current: &[[f64; N]], d: &[[f64; N]], tolerance: f64) -> f64 { + let mut s = 1.0f64; + for (slot, &at) in self.points.iter().enumerate() { + let (Some(dv), Some(cv), Some(iv)) = (d.get(slot), current.get(at), input.get(at)) else { continue }; + let a = dot(dv, dv); + if a <= 0.0 { + continue; + } + let delta: [f64; N] = std::array::from_fn(|c| cv[c] - iv[c]); + let b = 2.0 * dot(&delta, dv); + let c = dot(&delta, &delta) - tolerance * tolerance; + let disc = b * b - 4.0 * a * c; + let root = if disc >= 0.0 { (-b + disc.sqrt()) / (2.0 * a) } else { 0.0 }; + // `f64::max` ignores a NaN, so a step whose arithmetic overflows stops the run + // instead of leaking a not-a-number into a control point. + s = s.min(root.max(0.0)); + } + (s * TOLERANCE_SLACK).clamp(0.0, 1.0) + } + + /// The farthest any interior control point moved from the design's own position, in mm. + fn movement(&self, input: &[[f64; N]], current: &[[f64; N]]) -> f64 { + self.points.iter().fold(0.0f64, |m, &at| match (input.get(at), current.get(at)) { + (Some(a), Some(b)) => m.max(dist(a, b)), + _ => m, + }) + } +} + +/// The control points of one direction a fairing run may move, as `[lo, hi)`. A direction of one +/// control point — a curve seen as a net one row deep — is all interior; otherwise `keep` layers +/// are held at each end: the boundary, and with `keep = 2` the layer next to it. +fn free_span(count: usize, keep: usize) -> Result<(usize, usize)> { + if count == 1 { + return Ok((0, 1)); + } + if count <= 2 * keep { + return Err(KernelError::Invalid(format!( + "fairing a net {count} control points across one direction and holding {keep} at each edge leaves nothing to move" + ))); + } + Ok((keep, count - keep)) +} + +/// A symmetric banded matrix in upper-triangular form: row `r` holds columns `r ..= r + bw`. +struct Band { + n: usize, + bw: usize, + /// The upper band, `n · (bw + 1)` numbers; the scaling by the largest diagonal is folded in + /// when the matrix is factorised. + a: Vec, + scale: f64, +} + +impl Band { + fn new(n: usize, bw: usize) -> Result { + let Some(len) = (bw + 1).checked_mul(n).filter(|&l| l <= MAX_FAIR_BAND) else { + return Err(KernelError::TooLarge(format!("a fairing system of {n} unknowns {bw} apart"))); + }; + Ok(Band { n, bw, a: vec![0.0; len], scale: 1.0 }) + } + + fn index(&self, r: usize, s: usize) -> Option { + let (r, s) = if r <= s { (r, s) } else { (s, r) }; + (s - r <= self.bw && s < self.n).then_some(r * (self.bw + 1) + (s - r)) + } + + fn get(&self, r: usize, s: usize) -> f64 { + self.index(r, s).and_then(|i| self.a.get(i)).copied().unwrap_or(0.0) + } + + fn set(&mut self, r: usize, s: usize, v: f64) { + if let Some(x) = self.index(r, s).and_then(|i| self.a.get_mut(i)) { + *x = v; + } + } + + /// Add `v` at `(r, s)`. Only `r ≤ s` is stored (the matrix is symmetric, so callers adding a + /// whole symmetric neighbourhood reach each stored entry once). + fn add(&mut self, r: usize, s: usize, v: f64) { + if r <= s + && let Some(x) = self.index(r, s).and_then(|i| self.a.get_mut(i)) + { + *x += v; + } + } + + /// `A = Uᵀ·U` in place, the entries scaled by the largest diagonal first (the scaling does + /// not change the solution, only the size of the pivots, so a net whose knots span microns or + /// kilometres factorises just as well). An error, never a panic, when a pivot is not + /// positive: the system is singular there and the geometry must not move. + fn factorised(mut self) -> Result { + let (n, bw) = (self.n, self.bw); + let scale = (0..n).fold(0.0f64, |m, r| m.max(self.get(r, r).abs())); + if !(scale > 0.0) || !scale.is_finite() { + return Err(KernelError::Invalid("the fairness energy does not depend on any interior control point: nothing to fair".into())); + } + self.scale = scale; + for v in self.a.iter_mut() { + *v /= scale; + } + for r in 0..n { + let mut d = self.get(r, r); + for k in r.saturating_sub(bw)..r { + let u = self.get(k, r); + d -= u * u; + } + if !(d > 0.0) { + return Err(KernelError::Invalid( + "the fairing equations are singular (an interior control point that does not move the geometry)".into(), + )); + } + let d = d.sqrt(); + self.set(r, r, d); + for s in r + 1..(r + bw + 1).min(n) { + let mut v = self.get(r, s); + for k in s.saturating_sub(bw)..r { + v -= self.get(k, r) * self.get(k, s); + } + self.set(r, s, v / d); + } + } + Ok(self) + } + + /// Solve `A·x = b` with the stored factorisation, in place. + fn solve(&self, b: &mut [f64]) -> Result<()> { + let (n, bw) = (self.n, self.bw); + if b.len() != n { + return Err(KernelError::Invalid(format!("{} fairing unknowns for a system of {n}", b.len()))); + } + for x in b.iter_mut() { + *x /= self.scale; + } + // Uᵀ·y = b, then U·x = y. + for r in 0..n { + let mut v = b.get(r).copied().unwrap_or(0.0); + for k in r.saturating_sub(bw)..r { + v -= self.get(k, r) * b.get(k).copied().unwrap_or(0.0); + } + if let Some(x) = b.get_mut(r) { + *x = v / self.get(r, r); + } + } + for r in (0..n).rev() { + let mut v = b.get(r).copied().unwrap_or(0.0); + for s in r + 1..(r + bw + 1).min(n) { + v -= self.get(r, s) * b.get(s).copied().unwrap_or(0.0); + } + if let Some(x) = b.get_mut(r) { + *x = v / self.get(r, r); + } + } + if b.iter().any(|x| !x.is_finite()) { + return Err(KernelError::Invalid("the fairing step is not a finite number".into())); + } + Ok(()) + } +} + +/// Fair a B-spline curve's control points of any dimension: hold control points `0` and `n − 1` +/// and (with [`FairOptions::keep_edges`]) those next to them, and move the rest to lower the +/// bending energy `∫‖C″‖² dt` in the curve's own parameter (see the module docs). `N` is 3 for a +/// wire in space, 2 for a sketch curve. +pub fn fair_curve(degree: usize, knots: &[f64], ctrl: &[[f64; N]], opts: &FairOptions) -> Result> { + let opts = opts.checked()?; + let net = Net::new([degree, 0], [knots.to_vec(), TRIVIAL_KNOTS.to_vec()], [ctrl.len(), 1], &CURVE_ENERGY)?; + let (faired, report) = net.fair(ctrl, &opts)?; + Ok(FairCurve { degree, knots: knots.to_vec(), ctrl: faired, report }) +} + +/// Fair the interior of a B-spline control net with the thin-plate energy +/// `∫∫ ‖S_uu‖² + 2‖S_uv‖² + ‖S_vv‖² du dv` in the net's own parameters (see the module docs). +/// `count[0]` control points lie along u and `count[1]` along v, `ctrl` is u fastest, and the +/// boundary layers (see [`FairOptions::keep_edges`]) are held. The generic core behind +/// [`fair_surface`], [`fair_bspline_surface`] and their rational counterparts, for a caller with +/// its own control net. +pub fn fair_net( + degree: [usize; 2], + knots: [&[f64]; 2], + count: [usize; 2], + ctrl: &[[f64; N]], + opts: &FairOptions, +) -> Result> { + let opts = opts.checked()?; + let net = Net::new(degree, [knots[0].to_vec(), knots[1].to_vec()], count, &SURFACE_ENERGY)?; + let (faired, report) = net.fair(ctrl, &opts)?; + Ok(FairNet { degree, knots: [knots[0].to_vec(), knots[1].to_vec()], count, ctrl: faired, report }) +} + +/// Fair a truck B-spline curve's control points: [`fair_curve`] on its three-dimensional points, +/// rebuilt as a B-spline of the same knots. Guarded here itself, like every call into truck, so a +/// panic inside comes back as an error and the design is left alone. +pub fn fair_bspline_curve(curve: &mt::BSplineCurve, opts: &FairOptions) -> Result<(mt::BSplineCurve, FairReport)> { + guard("fair a B-spline curve", || { + let ctrl: Vec<[f64; 3]> = curve.control_points().iter().map(|p| [p.x, p.y, p.z]).collect(); + let faired = fair_curve(curve.degree(), curve.knot_vec(), &ctrl, opts)?; + let points: Vec = faired.ctrl.iter().map(|c| mt::Point3::new(c[0], c[1], c[2])).collect(); + let curve = mt::BSplineCurve::try_new(mt::KnotVec::from(faired.knots), points).map_err(|e| KernelError::Failed(format!("fair: {e}")))?; + Ok((curve, faired.report)) + }) +} + +/// Fair a truck rational curve's control points: its Cartesian control points are faired like any +/// other curve and its weights are kept, so the result is the same rational curve with moved +/// control points. The convex-hull bound that keeps the geometry within the tolerance still +/// holds, because with the weights unchanged the rational basis functions are non-negative and +/// sum to one. Guarded like [`fair_bspline_curve`]. +pub fn fair_nurbs_curve(curve: &mt::NurbsCurve, opts: &FairOptions) -> Result<(mt::NurbsCurve, FairReport)> { + guard("fair a NURBS curve", || { + let count = curve.control_points().len(); + let mut ctrl = Vec::with_capacity(count); + let mut weights = Vec::with_capacity(count); + for h in curve.control_points() { + let Some(p) = (HPoint { x: h.x, y: h.y, z: h.z, w: h.w }).project() else { + return Err(KernelError::Invalid("a control point has no position: its weight is zero".into())); + }; + ctrl.push([p.x, p.y, p.z]); + weights.push(h.w); + } + let faired = fair_curve(curve.degree(), curve.knot_vec(), &ctrl, opts)?; + let points: Vec = faired.ctrl.iter().zip(&weights).map(|(c, w)| mt::Vector4::new(c[0] * w, c[1] * w, c[2] * w, *w)).collect(); + let bsp = mt::BSplineCurve::try_new(mt::KnotVec::from(faired.knots), points).map_err(|e| KernelError::Failed(format!("fair: {e}")))?; + Ok((mt::NurbsCurve::new(bsp), faired.report)) + }) +} + +/// Fair a truck B-spline surface's control net: [`fair_net`] on its points, rebuilt as a B-spline +/// of the same knots, with the boundary rows held. Guarded like [`fair_bspline_curve`]. +pub fn fair_bspline_surface(surface: &mt::BSplineSurface, opts: &FairOptions) -> Result<(mt::BSplineSurface, FairReport)> { + guard("fair a B-spline surface", || { + let rows = surface.control_points(); + let count = [rows.len(), rows.first().map_or(0, Vec::len)]; + if count[0] == 0 || count[1] == 0 { + return Err(KernelError::Invalid("the surface has no control points".into())); + } + if rows.iter().any(|r| r.len() != count[1]) { + return Err(KernelError::Invalid("the surface's control net is not rectangular".into())); + } + // truck keeps the net as control_points()[i][j] with i along u; fair_net wants u fastest. + let mut ctrl = Vec::with_capacity(count[0] * count[1]); + for j in 0..count[1] { + for i in 0..count[0] { + let Some(p) = rows.get(i).and_then(|r| r.get(j)) else { continue }; + ctrl.push([p.x, p.y, p.z]); + } + } + let (ku, kv) = surface.knot_vecs(); + let (pu, pv) = surface.degrees(); + let faired = fair_net([pu, pv], [ku, kv], count, &ctrl, opts)?; + let mut net: Vec> = Vec::with_capacity(count[0]); + for i in 0..count[0] { + let mut row = Vec::with_capacity(count[1]); + for j in 0..count[1] { + let Some(c) = faired.ctrl.get(j * count[0] + i) else { + return Err(KernelError::Failed("fairing lost a control point".into())); + }; + row.push(mt::Point3::new(c[0], c[1], c[2])); + } + net.push(row); + } + let surface = mt::BSplineSurface::try_new((ku.clone(), kv.clone()), net).map_err(|e| KernelError::Failed(format!("fair: {e}")))?; + Ok((surface, faired.report)) + }) +} + +/// Fair a geom NURBS surface's control net: [`fair_net`] on its Cartesian control points, with +/// the weights kept as they are (a fairing moves points, not weights) and the boundary rows held. +/// A rational surface's geometry then stays within the tolerance too, by the same convex-hull +/// bound as the non-rational case. +pub fn fair_surface(surface: &NurbsSurface, opts: &FairOptions) -> Result<(NurbsSurface, FairReport)> { + let [cu, cv] = surface.count(); + let mut ctrl = Vec::with_capacity(cu * cv); + let mut weights = Vec::with_capacity(cu * cv); + for j in 0..cv { + for i in 0..cu { + let h = surface.point(i, j); + let Some(p) = h.project() else { + return Err(KernelError::Invalid("a control point has no position: its weight is zero".into())); + }; + ctrl.push([p.x, p.y, p.z]); + weights.push(h.w); + } + } + let [ku, kv] = surface.knots(); + let degree = surface.degree(); + let faired = fair_net(degree, [ku, kv], [cu, cv], &ctrl, opts)?; + let points: Vec = faired.ctrl.iter().zip(&weights).map(|(p, w)| HPoint::weighted(Vec3::new(p[0], p[1], p[2]), *w)).collect(); + let out = NurbsSurface::from_homogeneous(degree, surface.knots().clone(), points).map_err(|e| spline_error("fairing the surface", e))?; + Ok((out, faired.report)) +} + +/// Fair a truck rational surface's control net, like [`fair_surface`]: the Cartesian control +/// points are faired, the weights kept, and the result is a NURBS surface of the same knots. +/// Guarded like [`fair_bspline_curve`]. +pub fn fair_nurbs_surface(surface: &mt::NurbsSurface, opts: &FairOptions) -> Result<(mt::NurbsSurface, FairReport)> { + guard("fair a NURBS surface", || { + let bsp = surface.non_rationalized(); + let rows = bsp.control_points(); + let count = [rows.len(), rows.first().map_or(0, Vec::len)]; + if count[0] == 0 || count[1] == 0 { + return Err(KernelError::Invalid("the surface has no control points".into())); + } + if rows.iter().any(|r| r.len() != count[1]) { + return Err(KernelError::Invalid("the surface's control net is not rectangular".into())); + } + let mut ctrl = Vec::with_capacity(count[0] * count[1]); + let mut weights = Vec::with_capacity(count[0] * count[1]); + for j in 0..count[1] { + for i in 0..count[0] { + let Some(h) = rows.get(i).and_then(|r| r.get(j)) else { continue }; + let Some(p) = (HPoint { x: h.x, y: h.y, z: h.z, w: h.w }).project() else { + return Err(KernelError::Invalid("a control point has no position: its weight is zero".into())); + }; + ctrl.push([p.x, p.y, p.z]); + weights.push(h.w); + } + } + let (ku, kv) = bsp.knot_vecs(); + let (pu, pv) = bsp.degrees(); + let faired = fair_net([pu, pv], [ku, kv], count, &ctrl, opts)?; + let mut net: Vec> = Vec::with_capacity(count[0]); + for i in 0..count[0] { + let mut row = Vec::with_capacity(count[1]); + for j in 0..count[1] { + let (Some(c), Some(w)) = (faired.ctrl.get(j * count[0] + i), weights.get(j * count[0] + i)) else { + return Err(KernelError::Failed("fairing lost a control point".into())); + }; + row.push(mt::Vector4::new(c[0] * w, c[1] * w, c[2] * w, *w)); + } + net.push(row); + } + let bsp = mt::BSplineSurface::try_new((ku.clone(), kv.clone()), net).map_err(|e| KernelError::Failed(format!("fair: {e}")))?; + Ok((mt::NurbsSurface::new(bsp), faired.report)) + }) +} + +/// Fair the B-spline faces at `faces` (indices in `body`'s face order) of `body`: each face's +/// interior control vertices move to flatten the thin-plate energy of the face while its boundary +/// rows stay put, so the face keeps meeting its neighbours as it did. +/// +/// Every picked face keeps the very wires, edges and vertices it had and the way it pointed: only +/// the surface between them changes, so the body is as closed (or as open) as it was, holes and +/// all. A face that is not a B-spline (a plane, cylinder, cone, sphere, torus or swept face) has no +/// control points to fair and is refused by name; a mesh body has no faces to fair at all. +/// +/// The report is over the picked faces together: the energies summed, the largest movement of any +/// control point, and the energy trace as the sum at each step (a face that stops early keeps its +/// last energy for the steps the others still take). +pub fn fair_faces(body: &Body, faces: &[usize], opts: &FairOptions) -> Result<(Body, FairReport)> { + let opts = opts.checked()?; + body.require_brep("fair")?; + let all: Vec = body.solid.face_iter().cloned().collect(); + if all.is_empty() { + return Err(KernelError::Invalid("the body has no faces".into())); + } + let mut picked: Vec = Vec::new(); + for &i in faces { + if i >= all.len() { + return Err(KernelError::Invalid(format!("no face {i}: the body has {}", all.len()))); + } + if !picked.contains(&i) { + picked.push(i); + } + } + if picked.is_empty() { + return Err(KernelError::Invalid("pick the faces to fair".into())); + } + guard("fair faces", || { + let mut faired = all.clone(); + let mut reports: Vec = Vec::with_capacity(picked.len()); + for &i in &picked { + let Some(face) = all.get(i) else { continue }; + let (surface, report) = match face.surface() { + mt::Surface::BSplineSurface(s) => { + let (faired, report) = fair_bspline_surface(&s, &opts)?; + (mt::Surface::BSplineSurface(faired), report) + } + mt::Surface::NurbsSurface(s) => { + let (faired, report) = fair_nurbs_surface(&s, &opts)?; + (mt::Surface::NurbsSurface(faired), report) + } + _ => { + return Err(KernelError::Invalid(format!( + "face {i} is not a B-spline face: a plane, cylinder, cone, sphere, torus or swept face has no control points to fair" + ))); + } + }; + // The same face as before, wires and all, with the faired surface between them. + let mut nf = mt::Face::new_unchecked(face.absolute_boundaries().clone(), surface); + if !face.orientation() { + nf.invert(); + } + if let Some(slot) = faired.get_mut(i) { + *slot = nf; + } + reports.push(report); + } + let report = FairReport { + energy_before: reports.iter().map(|r| r.energy_before).sum(), + energy_after: reports.iter().map(|r| r.energy_after).sum(), + steps: reports.iter().map(|r| r.steps).max().unwrap_or(0), + movement: reports.iter().fold(0.0f64, |m, r| m.max(r.movement)), + energy_trace: { + let longest = reports.iter().map(|r| r.energy_trace.len()).max().unwrap_or(0); + (0..longest).map(|k| reports.iter().map(|r| r.energy_trace.get(k).copied().unwrap_or(r.energy_after)).sum()).collect() + }, + }; + // The body shell for shell, in its own face order: a body with holes keeps them. + let mut shells: Vec = Vec::with_capacity(body.solid.boundaries().len()); + let mut index = 0usize; + for shell in body.solid.boundaries() { + let mut shell_faces = Vec::with_capacity(shell.face_iter().count()); + for f in shell.face_iter() { + let i = index; + index += 1; + shell_faces.push(faired.get(i).cloned().unwrap_or_else(|| f.clone())); + } + shells.push(mt::Shell::from(shell_faces)); + } + // The topology is the body's own, face for face, so the solid is built as it (the shells + // are as closed, or as open, as they were) and only the shape is measured anew. + let mut result = body.clone(); + result.solid = std::sync::Arc::new(mt::Solid::new_unchecked(shells)); + result.size = 0.0; + result.volume = Some(std::sync::Arc::new(std::sync::Mutex::new(None))); + result.coarse_missing = None; + result.history = None; + let faces: Vec = + (0..all.len()).map(|i| crate::record::from_face(0, i, if picked.contains(&i) { Change::Modified } else { Change::Kept })).collect(); + let result = crate::record::attach(result, &[crate::record::Operand::new(body)], Some(faces), std::collections::HashMap::new())?; + Ok((result, report)) + }) +} + +#[cfg(test)] +mod tests { + use super::*; + use mt::ParametricCurve; + + /// A clamped uniform knot vector over `[0, spans]` of degree `degree`: `spans + degree` + /// control points. + fn clamped(spans: usize, degree: usize) -> Vec { + let mut k = vec![0.0; degree + 1]; + k.extend((1..spans).map(|i| i as f64)); + k.extend(std::iter::repeat_n(spans as f64, degree + 1)); + k + } + + /// The knot average `(kᵢ₊₁ + … + kᵢ₊ₚ)/p` of control point `i`: `Σ Nᵢ(t)·ξᵢ = t` (NURBS Book + /// eq. 3.4), so control points at `a + b·ξᵢ` describe the affine curve `a + b·t` exactly. + fn average(knots: &[f64], degree: usize, i: usize) -> f64 { + (1..=degree).map(|r| knots.get(i + r).copied().unwrap_or(0.0)).sum::() / degree as f64 + } + + /// The control points of the affine curve `ξ ↦ a + b·ξ` at the knot averages of `knots`: a + /// curve of exactly zero bending energy. + fn affine_curve(knots: &[f64], degree: usize, a: [f64; 3], b: [f64; 3]) -> Vec<[f64; 3]> { + (0..knots.len() - degree - 1).map(|i| std::array::from_fn(|c| a[c] + b[c] * average(knots, degree, i))).collect() + } + + /// The control net (u fastest) of the plane `(ξ, η) ↦ a + b·ξ + c·η` at the knot averages of + /// `ku` and `kv`: a net of exactly zero thin-plate energy. + fn planar_net(ku: &[f64], kv: &[f64], degree: [usize; 2], a: [f64; 3], b: [f64; 3], c: [f64; 3]) -> Vec<[f64; 3]> { + let (cu, cv) = (ku.len() - degree[0] - 1, kv.len() - degree[1] - 1); + let mut net = Vec::with_capacity(cu * cv); + for j in 0..cv { + for i in 0..cu { + let (xi, eta) = (average(ku, degree[0], i), average(kv, degree[1], j)); + net.push(std::array::from_fn(|d| a[d] + b[d] * xi + c[d] * eta)); + } + } + net + } + + /// The point of the curve `(degree, knots, ctrl)` at `t`, through the same basis geom uses: a + /// basis function `i` owns the `degree + 2` knots `knots[i..=i + degree + 1]`. + fn curve_at(degree: usize, knots: &[f64], ctrl: &[[f64; N]], t: f64) -> [f64; N] { + let n = ctrl.len(); + let side = if t >= knots.get(n).copied().unwrap_or(f64::INFINITY) { Side::Left } else { Side::Right }; + let mut out = [0.0; N]; + for (i, c) in ctrl.iter().enumerate() { + let b = knots.get(i..=i + degree + 1).and_then(|l| local_basis(l, t, side, 0).first().copied()).unwrap_or(0.0); + for (o, x) in out.iter_mut().zip(c) { + *o += b * x; + } + } + out + } + + /// The greatest distance between two points of space, in mm. + fn apart(a: mt::Point3, b: mt::Point3) -> f64 { + ((a.x - b.x).powi(2) + (a.y - b.y).powi(2) + (a.z - b.z).powi(2)).sqrt() + } + + /// A geom point as a coordinate triple. + fn triple(p: Vec3) -> [f64; 3] { + [p.x, p.y, p.z] + } + + #[test] + fn a_perturbed_curve_is_faired_back_to_its_smooth_shape_and_its_ends_stay_put() { + // A cubic curve whose control points sit at the knot averages of an affine curve: zero + // bending energy, so a fairing run must return to the same control polygon. The knots are + // deliberately uneven, so a fairing that ignored the curve's own parameter would not. + let knots = [0.0, 0.0, 0.0, 0.0, 1.0, 3.0, 5.0, 9.0, 9.0, 9.0, 9.0]; + let (a, b) = ([-20.0, 5.0, 1.5], [4.0, -2.5, 3.0]); + let smooth = affine_curve(&knots, 3, a, b); + let mut rough = smooth.clone(); + rough[2] = [rough[2][0] + 1.5, rough[2][1] - 2.0, rough[2][2] + 3.0]; + rough[3][2] += 2.5; + rough[4][0] -= 1.25; + rough[5][1] += 1.75; + let far = (0..smooth.len()).map(|i| dist(&rough[i], &smooth[i])).fold(0.0, f64::max); + + let opts = FairOptions::new(1.0, 8, false, 100.0).unwrap(); + let faired = fair_curve(3, &knots, &rough, &opts).unwrap(); + // The ends did not move at all, bit for bit. + assert_eq!(faired.ctrl[0], rough[0]); + assert_eq!(faired.ctrl[6], rough[6]); + // Every interior control point came back to the smooth shape. + let near = (1..6).map(|i| dist(&faired.ctrl[i], &smooth[i])).fold(0.0, f64::max); + eprintln!( + "curve: perturbation {far} mm, faired to {near:e} mm, energy {:e} → {:e} in {} steps, movement {:e} mm", + faired.report.energy_before, faired.report.energy_after, faired.report.steps, faired.report.movement + ); + assert!(near < 1e-6 && near < far * 1e-6, "faired to {near:e}, was {far} mm out"); + // The energy is a real number before, fell to nothing, and never rose on the way. + assert!(faired.report.energy_before > 1e-3, "{:e}", faired.report.energy_before); + assert!(faired.report.energy_after >= 0.0 && faired.report.energy_after < 1e-9 * faired.report.energy_before); + assert!(faired.report.steps >= 1 && faired.report.energy_trace.len() == faired.report.steps); + assert_eq!(faired.report.energy_trace.last().copied(), Some(faired.report.energy_after)); + assert!(faired.report.energy_trace.windows(2).all(|w| w[1] < w[0]), "{:?}", faired.report.energy_trace); + // A step never moves a control point further than the tolerance, and here it moved the + // whole way back to the smooth shape. + assert!(faired.report.movement <= 100.0 && (faired.report.movement - far).abs() < 1e-9, "{:e}", faired.report.movement); + // The faired curve is the affine line it was meant to be. + let mut off = 0.0f64; + for k in 0..=32 { + let t = 9.0 * f64::from(k) / 32.0; + let want: [f64; 3] = std::array::from_fn(|c| a[c] + b[c] * t); + off = off.max(dist(&curve_at(3, &knots, &faired.ctrl, t), &want)); + } + eprintln!("curve: the faired curve is {off:e} mm off the affine line it describes"); + assert!(off < 1e-6, "{off:e}"); + } + + #[test] + fn a_half_strength_run_walks_to_the_same_shape_over_its_iterations() { + // The same perturbation, walked a half step at a time: the energy must fall at every step + // (monotone, strictly) and still arrive at the smooth shape. + let knots = [0.0, 0.0, 0.0, 0.0, 1.0, 3.0, 5.0, 9.0, 9.0, 9.0, 9.0]; + let smooth = affine_curve(&knots, 3, [-20.0, 5.0, 1.5], [4.0, -2.5, 3.0]); + let mut rough = smooth.clone(); + rough[2][2] += 4.0; + rough[3][0] += 2.0; + let opts = FairOptions::new(0.5, 40, false, 100.0).unwrap(); + let faired = fair_curve(3, &knots, &rough, &opts).unwrap(); + eprintln!("half strength: energy {:e} → {:e} in {} steps", faired.report.energy_before, faired.report.energy_after, faired.report.steps); + assert!(faired.report.energy_trace.len() > 4, "one step is not a walk"); + assert!(faired.report.energy_trace.windows(2).all(|w| w[1] < w[0]), "the energy must fall at every step"); + assert!(faired.report.energy_after < 1e-9 * faired.report.energy_before); + assert_eq!(faired.ctrl[0], rough[0]); + assert_eq!(faired.ctrl[6], rough[6]); + assert!(dist(&faired.ctrl[3], &smooth[3]) < 1e-6); + } + + #[test] + fn the_faired_curve_stays_within_the_tolerance_of_the_design() { + // A heavily perturbed curve with a 0.75 mm tolerance: the design may be smoothed, but + // never by more than the caller allowed, and the bound is on the geometry itself. + let knots = clamped(3, 3); + let smooth = affine_curve(&knots, 3, [0.0, 0.0, 0.0], [30.0, 8.0, 0.0]); + let mut rough = smooth.clone(); + for (i, d) in [(1, [1.0, -2.0, 4.0]), (2, [-3.0, 2.0, -6.0]), (3, [2.0, -1.0, 5.0]), (4, [-2.0, 3.0, -4.0])] { + for c in 0..3 { + rough[i][c] += d[c]; + } + } + let tolerance = 0.75; + let opts = FairOptions::new(1.0, 12, false, tolerance).unwrap(); + let faired = fair_curve(3, &knots, &rough, &opts).unwrap(); + let mut off = 0.0f64; + for k in 0..=64 { + let t = 3.0 * f64::from(k) / 64.0; + off = off.max(dist(&curve_at(3, &knots, &faired.ctrl, t), &curve_at(3, &knots, &rough, t))); + } + eprintln!( + "tolerance {tolerance}: the curve moved {off:e} mm at worst ({:e} mm of control points), energy {:e} → {:e}", + faired.report.movement, faired.report.energy_before, faired.report.energy_after + ); + assert!(faired.report.movement <= tolerance, "{:e}", faired.report.movement); + assert!(off <= tolerance * (1.0 + 1e-9), "{off:e} over {tolerance}"); + assert!(faired.report.energy_after < faired.report.energy_before); + // The tolerance did not stop the run from heading for the smooth shape. + let before = (0..6).map(|i| dist(&rough[i], &smooth[i])).fold(0.0, f64::max); + let after = (0..6).map(|i| dist(&faired.ctrl[i], &smooth[i])).fold(0.0, f64::max); + assert!(after < before, "{after} vs {before}"); + assert_eq!(faired.ctrl[0], rough[0]); + assert_eq!(faired.ctrl[5], rough[5]); + } + + #[test] + fn a_sketch_curve_of_two_dimensions_fairs_the_same_way() { + // The generic core takes `[f64; 2]` just as well: a quadratic sketch curve, its middle + // pulled 2 mm out of the smooth polyline. + let knots = clamped(4, 2); + let smooth: Vec<[f64; 2]> = (0..knots.len() - 3).map(|i| [average(&knots, 2, i) * 10.0, average(&knots, 2, i) * 2.0]).collect(); + let mut rough = smooth.clone(); + rough[2] = [rough[2][0] + 1.0, rough[2][1] - 2.0]; + rough[3][1] += 1.5; + let opts = FairOptions::new(1.0, 10, false, 50.0).unwrap(); + let faired = fair_curve(2, &knots, &rough, &opts).unwrap(); + eprintln!("2D sketch curve: energy {:e} → {:e} in {} steps", faired.report.energy_before, faired.report.energy_after, faired.report.steps); + assert_eq!(faired.ctrl[0], rough[0]); + assert_eq!(faired.ctrl[5], rough[5]); + assert!(faired.report.energy_after >= 0.0 && faired.report.energy_after < 1e-9 * faired.report.energy_before); + for (i, want) in smooth.iter().enumerate().take(5).skip(1) { + assert!(dist(&faired.ctrl[i], want) < 1e-6, "point {i}"); + } + } + + #[test] + fn a_perturbed_net_is_faired_back_to_its_plane_with_its_boundary_rows_held() { + // A bicubic net of a plane (zero thin-plate energy), its middle pushed out; the boundary + // rows must not move at all and the interior must come back to the plane. + let (ku, kv) = (clamped(3, 3), clamped(4, 3)); + let degree = [3, 3]; + let (cu, cv) = (ku.len() - 4, kv.len() - 4); + let (a, b, c) = ([1.0, -2.0, 4.0], [2.0, 0.5, -0.25], [-0.5, 1.5, 0.75]); + let smooth = planar_net(&ku, &kv, degree, a, b, c); + let mut rough = smooth.clone(); + for &(i, j, d) in &[(2usize, 2usize, 3.0), (3, 3, -2.5), (1, 4, 1.5), (4, 5, -1.0)] { + rough[j * cu + i][2] += d; + } + let opts = FairOptions::new(1.0, 6, false, 100.0).unwrap(); + let faired = fair_net(degree, [&ku, &kv], [cu, cv], &rough, &opts).unwrap(); + for j in 0..cv { + for i in [0, cu - 1] { + assert_eq!(faired.ctrl[j * cu + i], rough[j * cu + i], "column {i}, row {j} moved"); + } + } + for i in 0..cu { + for j in [0, cv - 1] { + assert_eq!(faired.ctrl[j * cu + i], rough[j * cu + i], "row {j}, column {i} moved"); + } + } + let near = (0..cu * cv).map(|k| dist(&faired.ctrl[k], &smooth[k])).fold(0.0, f64::max); + let far = (0..cu * cv).map(|k| dist(&rough[k], &smooth[k])).fold(0.0, f64::max); + eprintln!("net: perturbation {far} mm, faired to {near:e} mm, energy {:e} → {:e}", faired.report.energy_before, faired.report.energy_after); + assert!(near < 1e-6, "{near:e}"); + assert!(faired.report.energy_after < 1e-9 * faired.report.energy_before); + assert!(faired.report.energy_trace.windows(2).all(|w| w[1] < w[0])); + // With keep_edges the two outer layers at each side stay exactly as designed — so the + // bump goes into the middle of the net — and the rest still comes back to the plane. + let mut inner = smooth.clone(); + for &(i, j, d) in &[(2usize, 2usize, 3.0), (3, 3, -2.5), (2, 3, 1.5), (3, 2, -1.0)] { + inner[j * cu + i][2] += d; + } + let kept = FairOptions::new(1.0, 6, true, 100.0).unwrap(); + let faired = fair_net(degree, [&ku, &kv], [cu, cv], &inner, &kept).unwrap(); + let shell = |i: usize, j: usize| i < 2 || i + 2 >= cu || j < 2 || j + 2 >= cv; + for j in 0..cv { + for i in 0..cu { + if shell(i, j) { + assert_eq!(faired.ctrl[j * cu + i], inner[j * cu + i], "({i}, {j}) is a boundary layer and moved"); + } + } + } + let near = (0..cu * cv).filter(|&k| !shell(k % cu, k / cu)).map(|k| dist(&faired.ctrl[k], &smooth[k])).fold(0.0, f64::max); + eprintln!("net with keep_edges: interior faired to {near:e} mm, energy {:e} → {:e}", faired.report.energy_before, faired.report.energy_after); + assert!(near < 1e-6 && faired.report.energy_after < faired.report.energy_before); + } + + #[test] + fn a_faired_geom_surface_stays_within_the_tolerance_of_the_design_and_keeps_its_weights() { + // A rational surface (weights off one) whose middle is pushed out: fair it with a 0.5 mm + // tolerance. The weights must be untouched, the boundary rows must not move, and the + // surface must stay within the tolerance of the one the design asked for. + let (ku, kv) = (clamped(3, 3), clamped(3, 3)); + let degree = [3, 3]; + let (cu, cv) = (ku.len() - 4, kv.len() - 4); + let smooth = planar_net(&ku, &kv, degree, [0.0, 0.0, 0.0], [20.0, 0.0, 0.0], [0.0, 15.0, 0.0]); + let weights: Vec = (0..cu * cv).map(|k| 0.5 + (k % 5) as f64 * 0.25).collect(); + let mut rough = smooth.clone(); + for &(i, j, d) in &[(2usize, 2usize, 4.0), (3, 3, -3.0), (1, 2, 2.0)] { + rough[j * cu + i][2] += d; + } + let surface = + NurbsSurface::new(degree, [ku.clone(), kv.clone()], &rough.iter().map(|c| Vec3::new(c[0], c[1], c[2])).collect::>(), &weights) + .unwrap(); + let tolerance = 0.5; + let opts = FairOptions::new(1.0, 8, false, tolerance).unwrap(); + let (faired, report) = fair_surface(&surface, &opts).unwrap(); + for k in 0..cu * cv { + let (i, j) = (k % cu, k / cu); + assert_eq!(faired.point(i, j).w, surface.point(i, j).w, "the weight of ({i}, {j}) changed"); + } + for j in 0..cv { + for i in [0, cu - 1] { + assert_eq!(faired.point(i, j), surface.point(i, j), "column {i} moved"); + } + } + let dom = surface.domain(); + let mut off = 0.0f64; + for k in 0..=24 { + for l in 0..=24 { + let (u, v) = (dom[0][0] + (dom[0][1] - dom[0][0]) * f64::from(k) / 24.0, dom[1][0] + (dom[1][1] - dom[1][0]) * f64::from(l) / 24.0); + off = off.max(surface.eval(u, v).unwrap().dist(faired.eval(u, v).unwrap())); + } + } + eprintln!( + "surface: moved {off:e} mm at worst ({:e} mm of control points), energy {:e} → {:e}", + report.movement, report.energy_before, report.energy_after + ); + assert!(report.movement <= tolerance, "{:e}", report.movement); + assert!(off <= tolerance * (1.0 + 1e-9), "{off:e} over {tolerance}"); + assert!(report.energy_after < report.energy_before); + // Every interior control point is nearer the plane than it was. + let before = (0..cu * cv).map(|k| dist(&rough[k], &smooth[k])).fold(0.0, f64::max); + let after = (0..cu * cv).map(|k| dist(&triple(faired.point(k % cu, k / cu).project().unwrap()), &smooth[k])).fold(0.0, f64::max); + assert!(after < before, "{after} vs {before}"); + } + + #[test] + fn the_truck_wrappers_fair_a_spline_curve_and_surface_with_their_ends_held() { + let knots = [0.0, 0.0, 0.0, 0.0, 1.0, 3.0, 5.0, 9.0, 9.0, 9.0, 9.0]; + let smooth = affine_curve(&knots, 3, [-20.0, 5.0, 1.5], [4.0, -2.5, 3.0]); + let mut rough = smooth.clone(); + rough[2][2] += 3.0; + rough[4][1] -= 2.0; + let points = |net: &[[f64; 3]]| -> Vec { net.iter().map(|c| mt::Point3::new(c[0], c[1], c[2])).collect() }; + let curve = mt::BSplineCurve::try_new(mt::KnotVec::from(knots.to_vec()), points(&rough)).unwrap(); + let tolerance = 5.0; + let opts = FairOptions::new(1.0, 8, false, tolerance).unwrap(); + let (faired, report) = fair_bspline_curve(&curve, &opts).unwrap(); + let mut off = 0.0f64; + for k in 0..=64 { + let t = 9.0 * f64::from(k) / 64.0; + off = off.max(apart(faired.subs(t), curve.subs(t))); + } + eprintln!( + "truck curve: energy {:e} → {:e} in {} steps, moved {off:e} mm ({:e} mm of control points)", + report.energy_before, report.energy_after, report.steps, report.movement + ); + assert_eq!(faired.control_points().first().copied(), curve.control_points().first().copied()); + assert_eq!(faired.control_points().last().copied(), curve.control_points().last().copied()); + assert!(off <= tolerance * (1.0 + 1e-9), "{off:e}"); + assert!(report.energy_after < 1e-9 * report.energy_before); + + // The same on a B-spline surface: boundary rows bit-identical, the bump in the middle + // comes down, the energy falls. + let (ku, kv) = (clamped(3, 3), clamped(3, 3)); + let degree = [3, 3]; + let (cu, cv) = (ku.len() - 4, kv.len() - 4); + let smooth = planar_net(&ku, &kv, degree, [0.0, 0.0, 0.0], [20.0, 0.0, 0.0], [0.0, 15.0, 0.0]); + let mut rough = smooth.clone(); + rough[2 * cu + 2][2] += 5.0; + let rows: Vec> = + (0..cu).map(|i| (0..cv).map(|j| mt::Point3::new(rough[j * cu + i][0], rough[j * cu + i][1], rough[j * cu + i][2])).collect()).collect(); + let surface = mt::BSplineSurface::try_new((mt::KnotVec::from(ku.clone()), mt::KnotVec::from(kv.clone())), rows).unwrap(); + let (faired, report) = fair_bspline_surface(&surface, &FairOptions::new(1.0, 6, true, 8.0).unwrap()).unwrap(); + eprintln!("truck surface: energy {:e} → {:e}, {} steps", report.energy_before, report.energy_after, report.steps); + for i in 0..cu { + for j in 0..cv { + if i < 2 || i + 2 >= cu || j < 2 || j + 2 >= cv { + assert_eq!(faired.control_points()[i][j], surface.control_points()[i][j], "({i}, {j}) is a boundary layer and moved"); + } + } + } + assert!(report.energy_after < report.energy_before); + assert!(faired.control_points()[2][2].z < surface.control_points()[2][2].z, "the bump should come down"); + } + + #[test] + fn a_rational_surface_fairs_with_its_boundary_and_tolerance_intact() { + // The truck NURBS wrapper, on weights away from one: fairing moves the Cartesian control + // points and keeps the weights, so the surface still cannot move more than the tolerance. + let (ku, kv) = (clamped(3, 3), clamped(3, 3)); + let degree = [3, 3]; + let (cu, cv) = (ku.len() - 4, kv.len() - 4); + let smooth = planar_net(&ku, &kv, degree, [0.0, 0.0, 0.0], [20.0, 0.0, 0.0], [0.0, 15.0, 0.0]); + let mut rough = smooth.clone(); + rough[2 * cu + 3][2] += 6.0; + let weights: Vec = (0..cu * cv).map(|k| 0.75 + (k % 3) as f64 * 0.5).collect(); + let rows: Vec> = (0..cu) + .map(|i| { + (0..cv) + .map(|j| { + let (c, w) = (rough[j * cu + i], weights[j * cu + i]); + mt::Vector4::new(c[0] * w, c[1] * w, c[2] * w, w) + }) + .collect() + }) + .collect(); + let bsp = mt::BSplineSurface::try_new((mt::KnotVec::from(ku.clone()), mt::KnotVec::from(kv.clone())), rows).unwrap(); + let surface = mt::NurbsSurface::new(bsp); + let (faired, report) = fair_nurbs_surface(&surface, &FairOptions::new(1.0, 6, false, 0.5).unwrap()).unwrap(); + eprintln!("rational surface: energy {:e} → {:e}, movement {:e} mm", report.energy_before, report.energy_after, report.movement); + assert!(report.movement <= 0.5 && report.energy_after < report.energy_before); + for i in 0..cu { + for j in 0..cv { + let (before, after) = (surface.non_rationalized().control_points()[i][j], faired.non_rationalized().control_points()[i][j]); + assert_eq!(after.w, before.w, "the weight of ({i}, {j}) changed"); + if i == 0 || i + 1 == cu || j == 0 || j + 1 == cv { + assert_eq!(after, before, "({i}, {j}) is on the boundary and moved"); + } + } + } + } + + #[test] + fn hostile_fairing_inputs_are_refused() { + let knots = clamped(3, 3); + let ctrl = affine_curve(&knots, 3, [0.0, 0.0, 0.0], [10.0, 4.0, 2.0]); + let ok = FairOptions::new(1.0, 4, false, 10.0).unwrap(); + // Out-of-range parameters: strength, iterations and tolerance are all input-derived. + for strength in [0.0, -0.25, 1.000_001, f64::NAN, f64::INFINITY] { + assert!(FairOptions::new(strength, 4, false, 10.0).is_err(), "strength {strength}"); + } + for iterations in [0, MAX_FAIR_ITERATIONS + 1, usize::MAX] { + assert!(FairOptions::new(1.0, iterations, false, 10.0).is_err(), "iterations {iterations}"); + } + for tolerance in [0.0, -1.0, f64::NAN, f64::INFINITY] { + assert!(FairOptions::new(1.0, 4, false, tolerance).is_err(), "tolerance {tolerance}"); + } + // A control point that is not a finite number. + for bad in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] { + let mut nan = ctrl.clone(); + nan[2][1] = bad; + assert!(fair_curve(3, &knots, &nan, &ok).is_err(), "a control point of {bad} was faired"); + } + // Knot vectors that do not fit the control points or the degree. + assert!(fair_curve(3, &knots[..9], &ctrl, &ok).is_err(), "nine knots for six control points of degree three"); + assert!(fair_curve(3, &[0.0, 0.0, 0.0, 0.0, 2.0, 1.0, 3.0, 4.0, 4.0, 4.0], &ctrl, &ok).is_err(), "knots that decrease"); + assert!(fair_curve(3, &vec![f64::NAN; knots.len()], &ctrl, &ok).is_err(), "a knot that is not a number"); + assert!(fair_curve(9, &knots, &ctrl, &ok).is_err(), "a degree above the supported maximum"); + // Fewer than two interior control points: one free point (a triangle of a quadratic), and + // a 3 × 3 net whose two edge layers leave nothing at all. + let triangle = [[0.0, 0.0, 0.0], [5.0, 2.0, 0.0], [10.0, 0.0, 1.0]]; + let one = fair_curve(2, &[0.0, 0.0, 0.0, 1.0, 1.0, 1.0], &triangle, &ok).unwrap_err(); + assert!(one.to_string().contains("at least two interior"), "{one}"); + let flat = [0.0, 0.0, 1.0, 2.0, 2.0]; + let grid = [[0.0, 0.0, 0.0]; 9]; + let kept = FairOptions::new(1.0, 4, true, 1.0).unwrap(); + assert!(fair_net([1, 1], [&flat, &flat], [3, 3], &grid, &kept).is_err(), "keep_edges on a 3 × 3 net"); + // A net too large to fair, refused before anything is allocated for it. + let line = [0.0, 0.0, 1.0, 1.0]; + assert!(fair_net([1, 1], [&line, &line], [65, 65], &vec![[0.0; 3]; 65 * 65], &ok).is_err(), "4225 control points"); + // A singular system: a cubic whose last control points lie past the end of its domain + // (0.5), where they do not move the curve at all, so the interior equations have no + // unique solution and the run must refuse rather than move anything. + let odd = [0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.5, 1.0, 1.0, 1.0]; + let six: Vec<[f64; 3]> = (0..6).map(|i| [i as f64, 0.0, 1.0]).collect(); + let singular = fair_curve(3, &odd, &six, &ok).unwrap_err(); + eprintln!("hostile inputs: refused with {singular}"); + assert!(singular.to_string().contains("singular"), "{singular}"); + } + + #[test] + fn a_small_quadratic_net_and_a_quintic_curve_fair_the_same_way() { + // Degrees and band shapes away from the bicubic case. A 5 × 5 quadratic net has fewer + // rows than the degree, so the solver's band is narrower than `degree_u + degree_v · cols` + // (the coupling cannot reach that far in v); a quintic curve carries six basis functions + // per span instead of four. Both must reproduce a zero-energy shape exactly. + let (ku, kv) = (clamped(3, 2), clamped(3, 2)); + let degree = [2, 2]; + let (cu, cv) = (ku.len() - 3, kv.len() - 3); + assert_eq!((cu, cv), (5, 5)); + let smooth = planar_net(&ku, &kv, degree, [0.0, 0.0, 0.0], [12.0, 1.0, 0.0], [-2.0, 9.0, 0.0]); + let mut rough = smooth.clone(); + for &(i, j, d) in &[(1usize, 2usize, 2.5), (2, 2, -1.5), (3, 3, 3.0), (2, 1, -2.0)] { + rough[j * cu + i][2] += d; + } + let faired = fair_net(degree, [&ku, &kv], [cu, cv], &rough, &FairOptions::new(1.0, 8, false, 100.0).unwrap()).unwrap(); + let near = (0..cu * cv).map(|k| dist(&faired.ctrl[k], &smooth[k])).fold(0.0, f64::max); + eprintln!("5 × 5 quadratic net: faired to {near:e} mm, energy {:e} → {:e}", faired.report.energy_before, faired.report.energy_after); + assert!(near < 1e-6 && faired.report.energy_after < 1e-9 * faired.report.energy_before); + for j in 0..cv { + for i in 0..cu { + if i == 0 || i + 1 == cu || j == 0 || j + 1 == cv { + assert_eq!(faired.ctrl[j * cu + i], rough[j * cu + i], "({i}, {j}) is on the boundary and moved"); + } + } + } + // A quintic curve, its middle pulled out. + let knots = clamped(3, 5); + let smooth = affine_curve(&knots, 5, [-10.0, 0.0, 2.0], [6.0, 4.0, -1.0]); + let mut rough = smooth.clone(); + rough[3][1] += 3.0; + rough[4][2] -= 2.0; + let faired = fair_curve(5, &knots, &rough, &FairOptions::new(1.0, 8, false, 100.0).unwrap()).unwrap(); + let near = (1..smooth.len() - 1).map(|i| dist(&faired.ctrl[i], &smooth[i])).fold(0.0, f64::max); + eprintln!("quintic curve: faired to {near:e} mm, energy {:e} → {:e}", faired.report.energy_before, faired.report.energy_after); + assert!(near < 1e-6 && faired.report.energy_after < 1e-9 * faired.report.energy_before); + assert_eq!(faired.ctrl[0], rough[0]); + assert_eq!(faired.ctrl[rough.len() - 1], rough[rough.len() - 1]); + } + + /// A one-face sheet body of the bicubic surface `net` (u fastest) over the clamped knots `ku` + /// and `kv`, its sides the straight lines between the surface's corners — true while the + /// net's boundary rows describe lines, which is how a plane's control net is built here. + fn sheet(ku: &[f64], kv: &[f64], net: &[[f64; 3]]) -> Body { + let (cu, cv) = (ku.len() - 4, kv.len() - 4); + assert_eq!(net.len(), cu * cv, "the net does not fit the knots"); + let rows: Vec> = + (0..cu).map(|i| (0..cv).map(|j| mt::Point3::new(net[j * cu + i][0], net[j * cu + i][1], net[j * cu + i][2])).collect()).collect(); + let bsp = mt::BSplineSurface::try_new((mt::KnotVec::from(ku.to_vec()), mt::KnotVec::from(kv.to_vec())), rows).unwrap(); + let corner = |i: usize, j: usize| { + let c = net[j * cu + i]; + mt::Vertex::new(mt::Point3::new(c[0], c[1], c[2])) + }; + let v = [corner(0, 0), corner(cu - 1, 0), corner(cu - 1, cv - 1), corner(0, cv - 1)]; + let line = |a: &mt::Vertex, b: &mt::Vertex| mt::Edge::new(a, b, mt::Curve::Line(mt::Line(a.point(), b.point()))); + let wire: mt::Wire = vec![line(&v[0], &v[1]), line(&v[1], &v[2]), line(&v[2], &v[3]), line(&v[3], &v[0])].into(); + crate::surfaces::surface_body(vec![mt::Face::try_new(vec![wire], mt::Surface::BSplineSurface(bsp)).unwrap()]) + } + + #[test] + fn a_body_face_is_faired_with_its_edges_untouched_and_its_history_kept() { + // A sheet whose middle is bumped: fairing the face must smooth it without moving a single + // edge, vertex or neighbouring face of the body, and must leave a history saying so. + let (ku, kv) = (clamped(3, 3), clamped(3, 3)); + let (cu, cv) = (ku.len() - 4, kv.len() - 4); + let smooth = planar_net(&ku, &kv, [3, 3], [0.0, 0.0, 0.0], [20.0, 0.0, 0.0], [0.0, 15.0, 0.0]); + let mut rough = smooth.clone(); + for &(i, j, d) in &[(2usize, 2usize, 3.0), (3, 3, -2.0), (2, 3, 1.5)] { + rough[j * cu + i][2] += d; + } + let body = sheet(&ku, &kv, &rough); + let before: Vec<_> = + body.solid.face_iter().flat_map(|f| f.absolute_boundaries().iter()).flat_map(|w| w.edge_iter().map(|e| e.id())).collect(); + assert!(!before.is_empty()); + let (faired, report) = fair_faces(&body, &[0], &FairOptions::new(1.0, 6, false, 50.0).unwrap()).unwrap(); + let after: Vec<_> = + faired.solid.face_iter().flat_map(|f| f.absolute_boundaries().iter()).flat_map(|w| w.edge_iter().map(|e| e.id())).collect(); + eprintln!( + "face: energy {:e} → {:e} in {} steps, movement {:e} mm, {} faces, history {}", + report.energy_before, + report.energy_after, + report.steps, + report.movement, + faired.face_count(), + faired.history().is_some() + ); + assert_eq!(before, after, "the body's edges must be the very ones it had"); + assert_eq!(faired.face_count(), 1); + assert!(faired.history().is_some(), "the face came from the one it was faired from"); + assert!(report.energy_after < 1e-9 * report.energy_before); + // The face's own surface is back on the plane, boundary rows bit for bit as they were. + let Some(mt::Surface::BSplineSurface(s)) = faired.solid.face_iter().next().map(|f| f.surface()) else { panic!("a B-spline face") }; + for j in 0..cv { + for i in 0..cu { + let p = s.control_points()[i][j]; + let want = if i == 0 || i + 1 == cu || j == 0 || j + 1 == cv { rough[j * cu + i] } else { smooth[j * cu + i] }; + assert!(dist(&[p.x, p.y, p.z], &want) < 1e-6, "({i}, {j})"); + } + } + // A solid keeps its shells, its faces and its edges: thickening the sheet gives a closed + // body whose top face is still the B-spline one, and fairing that face leaves the body a + // solid of the same shape. + let solid = crate::thicken(&body, 2.0, false).unwrap(); + assert!(!solid.is_surface()); + let top = + (0..solid.face_count()).find(|&i| matches!(solid.solid.face_iter().nth(i).map(|f| f.surface()), Some(mt::Surface::BSplineSurface(_)))); + let Some(top) = top else { panic!("the thickened sheet keeps its B-spline face") }; + let edges: Vec<_> = + solid.solid.face_iter().flat_map(|f| f.absolute_boundaries().iter()).flat_map(|w| w.edge_iter().map(|e| e.id())).collect(); + let (faired, report) = fair_faces(&solid, &[top], &FairOptions::new(1.0, 6, false, 0.5).unwrap()).unwrap(); + eprintln!("thickened solid: energy {:e} → {:e}, movement {:e} mm", report.energy_before, report.energy_after, report.movement); + assert!(!faired.is_surface() && faired.face_count() == solid.face_count() && faired.shell_count() == solid.shell_count()); + assert_eq!( + faired.solid.face_iter().flat_map(|f| f.absolute_boundaries().iter()).flat_map(|w| w.edge_iter().map(|e| e.id())).collect::>(), + edges + ); + assert!(report.energy_after < report.energy_before); + // A face with no control points to fair is refused, by name. + let cube = crate::box_solid(Vec3::ZERO, Vec3::new(10.0, 10.0, 10.0)).unwrap(); + let refused = fair_faces(&cube, &[0], &FairOptions::new(1.0, 1, false, 1.0).unwrap()).unwrap_err(); + eprintln!("cube face: {refused}"); + assert!(refused.to_string().contains("not a B-spline face"), "{refused}"); + assert!(fair_faces(&body, &[7], &FairOptions::new(1.0, 1, false, 1.0).unwrap()).is_err(), "no such face"); + assert!(fair_faces(&body, &[], &FairOptions::new(1.0, 1, false, 1.0).unwrap()).is_err(), "no face picked"); + } + + #[test] + fn an_already_fair_curve_is_left_where_it_is() { + // Zero bending energy, no reason to move: nothing churns, and an input at the minimum + // cannot be faired into a worse shape. + let knots = clamped(4, 3); + let smooth = affine_curve(&knots, 3, [1.0, 2.0, 3.0], [5.0, -1.0, 0.5]); + let opts = FairOptions::new(1.0, 20, false, 10.0).unwrap(); + let faired = fair_curve(3, &knots, &smooth, &opts).unwrap(); + eprintln!( + "already fair: energy {:e} → {:e}, {} steps, {:e} mm", + faired.report.energy_before, faired.report.energy_after, faired.report.steps, faired.report.movement + ); + assert!(faired.report.energy_after <= faired.report.energy_before); + assert!(faired.report.energy_after >= 0.0); + assert!(faired.report.movement < 1e-9, "{:e}", faired.report.movement); + for (i, want) in smooth.iter().enumerate() { + assert!(dist(&faired.ctrl[i], want) < 1e-9, "control point {i}"); + } + } +} diff --git a/crates/kernel/src/lib.rs b/crates/kernel/src/lib.rs index 0376c6a7..22756f87 100644 --- a/crates/kernel/src/lib.rs +++ b/crates/kernel/src/lib.rs @@ -24,6 +24,7 @@ mod boolean_diagnostics; mod boolean_guard; mod boolean_record; mod bore_guard; +mod boundary; mod bridge; mod build; mod clearance; @@ -36,6 +37,7 @@ mod emboss; mod face_params; #[cfg(test)] mod face_source_tests; +mod fair; mod fill; mod freeform; pub mod fuzz; @@ -66,6 +68,7 @@ mod issue_871_tests; mod join_between_caps_tests; mod lip; mod loopblend; +mod matchsurf; mod measure; mod mesh_convert; mod meshbody; @@ -125,6 +128,7 @@ pub use blend_corner::fillet_radii; pub use blend_variable::{RadiusSpline, Taper, fillet_chord, fillet_tapered, fillet_variable, fillet_variable_points, fillet_variable_spline, taper}; pub use body::{Body, CapSignature, DISPLAY_ANGLE, EdgeInfo, FaceInfo, FacePaint, FaceSource, Paint}; pub use boolean_diagnostics::{BooleanDiagnostics, BooleanFallback, BooleanStrategy, capture_boolean_diagnostics}; +pub use boundary::{Boundary, boundary_patch}; pub use bridge::{Continuity, Frame, blend_between, bridge, chain, chain_closes, frames, frames_across, sample_grid}; pub use build::{ PathSeg, box_solid, canonical_revolve_axis, cylinder, extrude, extrude_tapered, extrude_tapered_between, loft, revolve, sphere, sweep, @@ -132,11 +136,16 @@ pub use build::{ }; pub use delete_face::delete_faces; pub use emboss::{cylinder_emboss_anchor, emboss_cylinder, emboss_cylinder_at, has_cylinder_emboss_wall, infer_cylinder_emboss}; +pub use fair::{ + FairCurve, FairNet, FairOptions, FairReport, MAX_FAIR_ITERATIONS, MAX_FAIR_POINTS, fair_bspline_curve, fair_bspline_surface, fair_curve, + fair_faces, fair_net, fair_nurbs_curve, fair_nurbs_surface, fair_surface, +}; pub use fill::{PATCH_ANGLE_TOL_DEG, PatchReport, PatchSpec, patch_edges_crowned, patch_edges_report, patch_edges_with}; pub use helix::{loft_to_point, sweep_helix, sweep_helix_mesh}; pub use iges_in::{MAX_IGES_BYTES, iges_import, iges_to_step}; pub use iges_out::iges_export_bodies; pub use lip::raise_band; +pub use matchsurf::{match_surface, seam_continuity}; pub use measure::{BodyMeasure, measure}; pub use mesh_convert::convert_prismatic_mesh; pub use meshbody::{MAX_TRIANGLES, mesh_body}; diff --git a/crates/kernel/src/matchsurf.rs b/crates/kernel/src/matchsurf.rs new file mode 100644 index 00000000..8134ea34 --- /dev/null +++ b/crates/kernel/src/matchsurf.rs @@ -0,0 +1,601 @@ +//! Surface matching: rebuild a surface body's face so that one boundary of it meets a picked +//! edge chain of another body up to a requested continuity (G1, G2 or G3), while the boundary +//! opposite it does not move. +//! +//! A match is the bridge core ([`crate::bridge`]) pointed at a different pair of boundaries. One +//! end reads its position, tangent, curvature and change of curvature from the *target's* face +//! along the picked chain — that is the end that moves. The other end reads them from the body's +//! own surface along the boundary opposite it, built through the original boundary's own points: +//! a match may shape the middle, never that seam. In between the cross-sections are the bridge's +//! cubic Hermites, whose crown is free, so the crown whose surface stays closest to the original +//! one — sampled on a grid and measured against it — is the one kept: a match moves the seam and +//! as little else as the conditions allow. +//! +//! Which boundary of the face moves is read off the geometry: the one nearest the picked chain. +//! The order a match may ask for is read off the target the same way: a chain that several faces +//! of the target meet, with shapes that part company, is a seam the target itself does not hold +//! beyond that order, and no surface can be told to carry a condition the target has not got. + +use mt::{BoundedCurve, ParametricCurve, ParametricSurface, ParametricSurface3D, SearchNearestParameter}; +use solidcraft_geom::Vec3; +use truck_modeling as mt; + +use crate::body::{Body, from_p3, p3}; +use crate::bridge::{self, Continuity, Frame}; +use crate::record::Rebuild; +use crate::surfaces::surface_body; +use crate::{KernelError, Result, guard}; + +/// Stations a match reads along the seam and along the boundary it keeps. +const STATIONS: usize = 33; + +/// Samples each way of the grid the rebuilt surface is compared with the original on. +const GRID: usize = 7; + +/// The crowns the match tries at each end: [`crate::bridge::bridge`]'s free knob, which shapes +/// the middle without touching either end's conditions. +const CROWNS: [f64; 5] = [-0.5, 0.0, 0.5, 1.0, 2.0]; + +/// Points a caller may pick one face with. +const MAX_PICKS: usize = 64; + +/// Samples per edge of a chain, for lengths and distances. +const PER_EDGE: usize = 32; + +/// A seam is credited with G1 when its normals agree this closely (degrees), with G2 when the +/// two curvature vectors agree within this fraction of the larger of them, and with G3 when the +/// change of curvature agrees within this fraction of the curvature over the body's size. +const G1_DEG: f64 = 0.05; +const G2_REL: f64 = 5e-3; +const G3_REL: f64 = 5e-2; + +/// Two faces of the target count as meeting each other at the picked chain when their normals +/// agree this closely, their curvature within this fraction of the larger, and their change of +/// curvature within this fraction of the curvature over the body's size. +const SUPPLY_DEG: f64 = 0.05; +const SUPPLY_REL: f64 = 2e-2; +const SUPPLY_BEND_REL: f64 = 1e-1; + +/// Rebuild the face of `body` picked by `faces` (points on it) so that its boundary nearest the +/// chain of `target`'s edges picked by `target_pick` takes that chain's place and repeats the +/// target's own shape along it up to `order` — position, tangent, curvature vector, and the +/// change of that curvature vector along the direction the match leaves in. The boundary +/// opposite keeps its place and the original surface's own frame, and the middle is shaped to +/// stay as close to the original surface as the conditions allow. +pub fn match_surface(body: &Body, faces: &[Vec3], target: &Body, target_pick: &[Vec3], order: Continuity) -> Result { + body.require_brep("match")?; + target.require_brep("match")?; + if !body.is_surface() { + return Err(KernelError::Invalid("match rebuilds a surface body: this one is closed".into())); + } + let face = face_near(body, faces)?; + let chain = bridge::chain(target, target_pick)?; + let near = nearest_side(&face, &chain)?; + let far = near.opposite(); + let far_edges = side_edges(&face, far, body.size() * 1e-6 + 1e-9)?; + let (mid_near, mid_far) = (chain_mid(&chain)?, Side::at(&face.surface(), &ranges(&face)?, far, 0.5)); + let leave = (mid_far - mid_near) + .normalized() + .filter(|l| l.is_finite()) + .ok_or_else(|| KernelError::Invalid("the far boundary sits on the picked chain: there is nothing to match".into()))?; + // Which face of the target the match continues from, and how much of the requested order the + // target itself holds along that chain. + let candidates = support_faces(target, &chain, leave)?; + let supply = supply_order(&candidates, &chain, leave, target.size())?; + if order.conditions() > supply.conditions() { + return Err(KernelError::Invalid(format!( + "the target's own surfaces part company at the picked edges: they hold only {} there, so a seam cannot hold {}", + supply.as_str(), + order.as_str() + ))); + } + let Some((target_face, _)) = candidates.first() else { + return Err(KernelError::Invalid("no face of the target holds the picked edges".into())); + }; + // The stations of the seam: what the target does there along the direction the match leaves + // in, and what the far boundary does coming back the other way. + let mut near_frames = bridge::frames(target_face, &chain, leave, STATIONS)?; + let far_frames = bridge::frames(&face, &far_edges, leave * -1.0, STATIONS)?; + // The two chains are sampled from their own starts: turn one round when that pairs the + // stations better (a match is built station to station). + let direct: f64 = near_frames.iter().zip(&far_frames).map(|(a, b)| a.p.dist(b.p)).sum(); + let reversed: f64 = near_frames.iter().zip(far_frames.iter().rev()).map(|(a, b)| a.p.dist(b.p)).sum(); + if reversed < direct { + near_frames.reverse(); + } + guard("surface match", || { + // The crown is the free knob: keep the one whose surface stays closest to the original. + let mut best: Option<(f64, mt::Face)> = None; + for crown_near in CROWNS { + for crown_far in CROWNS { + let made = bridge::bridge(&near_frames, &far_frames, [order, Continuity::G1], [crown_near, crown_far])?; + let Some(made) = made.solid.face_iter().next().cloned() else { + return Err(KernelError::Failed("match: the bridge made no face".into())); + }; + let stray = deviation(&made, &face)?; + if best.as_ref().is_none_or(|(b, _)| stray < *b) { + best = Some((stray, made)); + } + } + } + let Some((_, made)) = best else { + return Err(KernelError::Failed("match: no crown built a surface".into())); + }; + let mut kept = made; + // The rebuilt face takes the target's own sense at the seam: the two surfaces then meet + // the same way round, which is what a stitch, a thicken or a display of the result reads. + let seam = chain_mid(&chain)?; + if let (Some((made_n, _)), Some((target_n, _))) = (normal_near(&kept, seam), normal_near(target_face, seam)) + && made_n.dot(target_n) < 0.0 + { + kept.invert(); + } + let mut rec = Rebuild::default(); + rec.face(&kept, &face); + Ok(crate::record::rebuilt_body(surface_body(vec![kept]), &[&body.solid, &target.solid], &rec)) + }) +} + +/// How well the rebuilt surface `made` meets `target` along the chain picked by `target_pick`: +/// the worst distance from the seam to that chain, the worst angle between the two surfaces' +/// normals across it (degrees), and the highest continuity the seam holds +/// ([`Continuity::G0`] when even the normals disagree). +pub fn seam_continuity(made: &Body, target: &Body, target_pick: &[Vec3]) -> Result<(f64, f64, Continuity)> { + made.require_brep("match")?; + target.require_brep("match")?; + let chain = bridge::chain(target, target_pick)?; + let poly = chain_polyline(&chain)?; + // The face whose boundary runs along the chain: the one nearest it, so that a body that has + // not been matched yet — where the seam is still some way off — can be measured too. + let mid = chain_mid(&chain)?; + let face = made + .solid + .face_iter() + .filter_map(|f| face_distance(f, mid).map(|d| (d, f))) + .min_by(|a, b| a.0.total_cmp(&b.0)) + .map(|(_, f)| f.clone()) + .ok_or_else(|| KernelError::Failed("match: the body has no faces".into()))?; + let near = nearest_side(&face, &chain)?; + let far = near.opposite(); + let ranges = ranges(&face)?; + let mid_far = Side::at(&face.surface(), &ranges, far, 0.5); + let leave = (mid_far - chain_mid(&chain)?) + .normalized() + .ok_or_else(|| KernelError::Invalid("the seam and the far boundary coincide: there is nothing to measure".into()))?; + let candidates = support_faces(target, &chain, leave)?; + let Some((target_face, _)) = candidates.first() else { + return Err(KernelError::Invalid("no face of the target holds the picked edges".into())); + }; + let frames = bridge::frames(target_face, &chain, leave, STATIONS)?; + let size = target.size(); + guard("surface match", || { + let surf = face.surface(); + let (mut gap, mut deg, mut curv, mut bend) = (0.0f64, 0.0f64, 0.0f64, 0.0f64); + for i in 0..STATIONS { + let t = i as f64 / (STATIONS - 1) as f64; + let p = Side::at(&surf, &ranges, near, t); + gap = gap.max(chain_distance(&poly, p)); + // What the target says there, and what the rebuilt surface actually does there. + let Some(frame) = frames.iter().min_by(|a, b| a.p.dist(p).total_cmp(&b.p.dist(p))) else { + return Err(KernelError::Failed("match: the target has no frames".into())); + }; + let (_, t_made, bend_made, change_made) = seam_shape(&surf, &ranges, near, t)?; + if let Some((n_made, _)) = normal_near(&face, p) + && let Some((n_target, _)) = normal_near(target_face, p) + { + deg = deg.max(angle_deg(n_made, n_target)); + } + let scale = frame.bend.len().max(bend_made.len()).max(1.0 / size); + curv = curv.max((bend_made - frame.bend).len() / scale); + // The change of curvature is what a surface can hold: only its normal part counts, + // and it turns round with the direction of travel. + let change = if t_made.dot(frame.leave) < 0.0 { change_made * -1.0 } else { change_made }; + let scale3 = frame.bend_change.len().max(change.len()).max(scale / size); + bend = bend.max((change - frame.bend_change).len() / scale3); + } + let mut order = Continuity::G0; + if deg <= G1_DEG { + order = Continuity::G1; + if curv <= G2_REL { + order = Continuity::G2; + if bend <= G3_REL { + order = Continuity::G3; + } + } + } + Ok((gap, deg, order)) + }) +} + +/// Which side of a face's parameter rectangle a boundary runs along: the match reads its face as +/// a tensor patch, so its four sides are its own parameter lines. +#[derive(Clone, Copy, PartialEq, Eq, Debug)] +enum Side { + /// The `u = u0` side (its points run along `v`). + U0, + /// The `u = u1` side. + U1, + /// The `v = v0` side (its points run along `u`). + V0, + /// The `v = v1` side. + V1, +} + +impl Side { + /// The four sides, in a fixed order. + const ALL: [Side; 4] = [Side::U0, Side::U1, Side::V0, Side::V1]; + + /// The side opposite this one: what a match keeps when this one moves. + fn opposite(self) -> Side { + match self { + Side::U0 => Side::U1, + Side::U1 => Side::U0, + Side::V0 => Side::V1, + Side::V1 => Side::V0, + } + } + + /// Whether this side pins `u` (its points run along `v`). + fn is_u(self) -> bool { + matches!(self, Side::U0 | Side::U1) + } + + /// The value this side pins, and the range the side runs over. + fn fixed_and_along(self, u: (f64, f64), v: (f64, f64)) -> (f64, (f64, f64)) { + match self { + Side::U0 => (u.0, v), + Side::U1 => (u.1, v), + Side::V0 => (v.0, u), + Side::V1 => (v.1, u), + } + } + + /// The point at the normalized station `t` along the side. + fn at(surf: &mt::Surface, ranges: &[(f64, f64); 2], side: Side, t: f64) -> Vec3 { + let (u, v) = (ranges[0], ranges[1]); + let (fixed, along) = side.fixed_and_along(u, v); + let t = along.0 + (along.1 - along.0) * t.clamp(0.0, 1.0); + if side.is_u() { from_p3(surf.subs(fixed, t)) } else { from_p3(surf.subs(t, fixed)) } + } +} + +/// The bounded parameter ranges of a face's surface. +fn ranges(face: &mt::Face) -> Result<[(f64, f64); 2]> { + let surf = face.surface(); + let (Some(u), Some(v)) = surf.try_range_tuple() else { + return Err(KernelError::Failed("match: the face's surface is unbounded".into())); + }; + Ok([u, v]) +} + +/// How far `p` is from the face's own surface. +fn face_distance(face: &mt::Face, p: Vec3) -> Option { + let surf = face.surface(); + let (u, v) = surf.search_nearest_parameter(p3(p), None, 100)?; + Some(from_p3(surf.subs(u, v)).dist(p)) +} + +/// The surface's unit normal at the point nearest `p`, and how far `p` is from that point. +fn normal_near(face: &mt::Face, p: Vec3) -> Option<(Vec3, f64)> { + let surf = face.surface(); + let (u, v) = surf.search_nearest_parameter(p3(p), None, 100)?; + let n = surf.normal(u, v); + let n = Vec3::new(n.x, n.y, n.z).normalized()?; + Some((n, from_p3(surf.subs(u, v)).dist(p))) +} + +/// The face of `body` the picks lie on: they must all pick the same one. A pick is a point the +/// caller clicked, so the face nearest it within the body's own pick tolerance is the face it +/// means (`crate::surfaces::edge_near` reads points the same way). +fn face_near(body: &Body, picks: &[Vec3]) -> Result { + if picks.is_empty() || picks.len() > MAX_PICKS || picks.iter().any(|p| !p.is_finite()) { + return Err(KernelError::Invalid(format!("pick 1 to {MAX_PICKS} points on the face to match"))); + } + let tol = body.size() * 0.05 + 1e-6; + let mut found: Option = None; + for p in picks { + let hit = body + .solid + .face_iter() + .filter_map(|f| face_distance(f, *p).map(|d| (d, f))) + .min_by(|a, b| a.0.total_cmp(&b.0)) + .filter(|(d, _)| *d <= tol) + .map(|(_, f)| f.clone()) + .ok_or_else(|| KernelError::Invalid(format!("no face of the body passes through {p:?}")))?; + if found.as_ref().is_some_and(|f| f.id() != hit.id()) { + return Err(KernelError::Invalid("the picks are on different faces: match one face".into())); + } + found = Some(hit); + } + found.ok_or_else(|| KernelError::Invalid("pick a point on the face to match".into())) +} + +/// A chain of edges as a polyline (each edge sampled), for lengths and distances. +fn chain_polyline(chain: &[mt::Edge]) -> Result> { + if chain.is_empty() { + return Err(KernelError::Invalid("the chain has no edges".into())); + } + let mut out = Vec::with_capacity(chain.len() * (PER_EDGE + 1)); + for e in chain { + let c = e.curve(); + let (t0, t1) = c.range_tuple(); + if !(t1 > t0) || ![t0, t1].into_iter().all(f64::is_finite) { + return Err(KernelError::Failed("match: an edge has no parameter range".into())); + } + for k in 0..=PER_EDGE { + out.push(from_p3(c.subs(t0 + (t1 - t0) * k as f64 / PER_EDGE as f64))); + } + } + Ok(out) +} + +/// The point at the middle of a chain, by length. +fn chain_mid(chain: &[mt::Edge]) -> Result { + let poly = chain_polyline(chain)?; + let mut cum = Vec::with_capacity(poly.len()); + cum.push(0.0); + for w in poly.windows(2) { + let (Some(a), Some(b)) = (w.first(), w.get(1)) else { continue }; + let Some(last) = cum.last().copied() else { break }; + cum.push(last + a.dist(*b)); + } + let total = cum.last().copied().unwrap_or(0.0); + if !(total > 1e-9) || !total.is_finite() { + return Err(KernelError::Invalid("the picked chain has no length".into())); + } + let half = total * 0.5; + let hi = cum.partition_point(|x| *x < half).clamp(1, poly.len().saturating_sub(1)); + let (Some(a), Some(b)) = (poly.get(hi - 1), poly.get(hi)) else { + return Err(KernelError::Failed("match: chain".into())); + }; + let (Some(&lo), Some(&up)) = (cum.get(hi - 1), cum.get(hi)) else { + return Err(KernelError::Failed("match: chain".into())); + }; + let span = up - lo; + let t = if span > 1e-12 { (half - lo) / span } else { 0.0 }; + Ok(*a + (*b - *a) * t) +} + +/// How far `p` is from a polyline. +fn chain_distance(poly: &[Vec3], p: Vec3) -> f64 { + let mut best = f64::INFINITY; + for w in poly.windows(2) { + let (Some(a), Some(b)) = (w.first(), w.get(1)) else { continue }; + let d = b.dist(*a); + let t = if d > 1e-12 { ((p - *a).dot(*b - *a) / (d * d)).clamp(0.0, 1.0) } else { 0.0 }; + best = best.min(p.dist(*a + (*b - *a) * t)); + } + best +} + +/// The side of the face's own parameter rectangle the picked chain sits nearest. +fn nearest_side(face: &mt::Face, chain: &[mt::Edge]) -> Result { + let surf = face.surface(); + let ranges = ranges(face)?; + let poly = chain_polyline(chain)?; + /// Samples of each side used to tell which one the chain runs along. + const SAMPLES: usize = 17; + let mut best: Option<(f64, Side)> = None; + for side in Side::ALL { + let mut total = 0.0; + for k in 0..SAMPLES { + let t = k as f64 / (SAMPLES - 1) as f64; + total += chain_distance(&poly, Side::at(&surf, &ranges, side, t)); + } + let mean = total / SAMPLES as f64; + if best.is_none_or(|(b, _)| mean < b) { + best = Some((mean, side)); + } + } + best.map(|(_, s)| s).ok_or_else(|| KernelError::Failed("match: the face has no sides".into())) +} + +/// The face's own edges that run along a side of its parameter rectangle, chained head to tail. +/// A match keeps the far side's own edges, so they must be a chain that runs on. +fn side_edges(face: &mt::Face, side: Side, tol: f64) -> Result> { + let surf = face.surface(); + let ranges = ranges(face)?; + let (fixed, along) = side.fixed_and_along(ranges[0], ranges[1]); + let band = (along.1 - along.0).abs().max(1e-12) * 1e-3; + let mut pool: Vec = Vec::new(); + for e in face.edge_iter() { + let c = e.curve(); + let (t0, t1) = c.range_tuple(); + let p = from_p3(c.subs(0.5 * (t0 + t1))); + let Some((u, v)) = surf.search_nearest_parameter(p3(p), None, 100) else { continue }; + if from_p3(surf.subs(u, v)).dist(p) > tol { + continue; + } + let got = if side.is_u() { u } else { v }; + if (got - fixed).abs() <= band { + pool.push(e.clone()); + } + } + if pool.is_empty() { + return Err(KernelError::Invalid("the face has no edges along the boundary a match would keep".into())); + } + let first = pool.remove(0); + let mut chain = vec![first]; + let mut end = match chain.first() { + Some(e) => e.back().clone(), + None => return Err(KernelError::Failed("match: chain".into())), + }; + while let Some(k) = pool.iter().position(|e| e.front().id() == end.id() || e.back().id() == end.id()) { + let e = pool.remove(k); + let e = if e.front().id() == end.id() { e } else { e.inverse() }; + end = e.back().clone(); + chain.push(e); + } + let mut start = match chain.first() { + Some(e) => e.front().clone(), + None => return Err(KernelError::Failed("match: chain".into())), + }; + while let Some(k) = pool.iter().position(|e| e.front().id() == start.id() || e.back().id() == start.id()) { + let e = pool.remove(k); + let e = if e.back().id() == start.id() { e } else { e.inverse() }; + start = e.front().clone(); + chain.insert(0, e); + } + if !pool.is_empty() { + return Err(KernelError::Invalid("the face's edges along that boundary do not run on from one another".into())); + } + Ok(chain) +} + +/// The faces of `target` that meet the chain and whose tangent plane holds the leave direction, +/// best first: the face a match continues from, and the ones it might have been read from. +fn support_faces(target: &Body, chain: &[mt::Edge], leave: Vec3) -> Result> { + let mid = chain_mid(chain)?; + let tol = target.size() * 1e-6 + 1e-9; + let mut out: Vec<(mt::Face, f64)> = Vec::new(); + for face in target.solid.face_iter() { + if !chain.iter().all(|e| face.edge_iter().any(|f| f.id() == e.id())) { + continue; + } + let Some((n, d)) = normal_near(face, mid) else { continue }; + if d > tol { + continue; + } + let score = (leave - n * leave.dot(n)).len(); + if score.is_finite() && score >= 0.2 { + out.push((face.clone(), score)); + } + } + if out.is_empty() { + return Err(KernelError::Invalid( + "no face of the target holds the picked edges in the direction the match leaves in: pick the chain of a face the surface can continue" + .into(), + )); + } + out.sort_by(|a, b| b.1.total_cmp(&a.1)); + Ok(out) +} + +/// The highest order the target itself holds along the chain: every face that meets those edges +/// and whose tangent plane holds the leave direction reads the chain's shape from its own +/// surface, and where two of them part company the target's own geometry jumps at the seam. +fn supply_order(candidates: &[(mt::Face, f64)], chain: &[mt::Edge], leave: Vec3, size: f64) -> Result { + let Some((first, best)) = candidates.first() else { + return Err(KernelError::Invalid("no face of the target holds the picked edges".into())); + }; + let keep: Vec<&mt::Face> = candidates.iter().filter(|(_, s)| *s >= 0.5 * best).map(|(f, _)| f).collect(); + let mut order = Continuity::G3; + if keep.len() < 2 { + return Ok(order); + } + let primary: Vec = bridge::frames(first, chain, leave, STATIONS)?; + for other in keep.iter().skip(1) { + // A face that cannot be read along the whole chain says nothing here: the face the match + // is built on is the one that must be able to supply it. + let Ok(other_frames) = bridge::frames(other, chain, leave, STATIONS) else { continue }; + let mut worst_deg = 0.0f64; + let mut worst = 0.0f64; + let mut worst_change = 0.0f64; + for (a, b) in primary.iter().zip(&other_frames) { + if let (Some((na, _)), Some((nb, _))) = (normal_near(first, a.p), normal_near(other, b.p)) { + worst_deg = worst_deg.max(angle_deg(na, nb)); + } + let scale = a.bend.len().max(b.bend.len()).max(1.0 / size); + worst = worst.max((a.bend - b.bend).len() / scale); + let scale3 = a.bend_change.len().max(b.bend_change.len()).max(scale / size); + worst_change = worst_change.max((a.bend_change - b.bend_change).len() / scale3); + } + if worst_deg > SUPPLY_DEG { + order = Continuity::G0; + } + if worst > SUPPLY_REL && order.conditions() > 1 { + order = Continuity::G1; + } + if worst_change > SUPPLY_BEND_REL && order.conditions() > 2 { + order = Continuity::G2; + } + } + Ok(order) +} + +/// The acute angle between two directions, in degrees. For two surfaces' normals this is the +/// angle between their tangent planes: a sheet has no inside, so which way a normal points says +/// nothing about whether the surfaces meet smoothly. +fn angle_deg(a: Vec3, b: Vec3) -> f64 { + let d = a.dot(b) / (a.len() * b.len()); + if !d.is_finite() { + return f64::NAN; + } + d.abs().clamp(0.0, 1.0).acos().to_degrees() +} + +/// The geometric shape of a surface curve at a point from its first three derivatives along the +/// curve: the same statement [`crate::bridge`] makes of a support surface — the unit tangent, the +/// curvature vector, and the arclength derivative of that curvature vector. +fn shape(d1: Vec3, d2: Vec3, d3: Vec3) -> Result<(Vec3, Vec3, Vec3)> { + let sigma = d1.len(); + if !(sigma > 1e-12) || !sigma.is_finite() { + return Err(KernelError::Failed("match: the seam's cross-section does not move".into())); + } + let t = d1 * (1.0 / sigma); + let tp = d2 * (1.0 / sigma) - d1 * (d1.dot(d2) / (sigma * sigma * sigma)); + let nn = d2 - t * d2.dot(t); + let np = d3 - t * d3.dot(t) - tp * d2.dot(t) - t * d2.dot(tp); + let kappa = np * (sigma * sigma) - nn * (2.0 * d1.dot(d2)); + let change = kappa * (1.0 / (sigma.powi(5))); + if ![t, nn, change].into_iter().all(Vec3::is_finite) { + return Err(KernelError::Failed("match: the seam's shape is not finite".into())); + } + Ok((t, nn * (1.0 / (sigma * sigma)), change)) +} + +/// The shape of a face's surface along its cross direction at the station `t`, read inwards from +/// the boundary `near` over five one-sided samples: the point, the unit tangent, the curvature +/// vector and the change of that curvature. The same finite difference [`crate::bridge::frames`] +/// reads a support surface with, so that a rebuilt seam can be compared with what it matched. +fn seam_shape(surf: &mt::Surface, ranges: &[(f64, f64); 2], near: Side, t: f64) -> Result<(Vec3, Vec3, Vec3, Vec3)> { + let (station, cross, at_start) = if near.is_u() { + (ranges[1].0 + (ranges[1].1 - ranges[1].0) * t, ranges[0], near == Side::U0) + } else { + (ranges[0].0 + (ranges[0].1 - ranges[0].0) * t, ranges[1], near == Side::V0) + }; + let span = (cross.1 - cross.0).abs().max(1e-12); + let h = span * 1e-3; + let c0 = if at_start { cross.0 } else { cross.1 }; + let sgn = if at_start { 1.0 } else { -1.0 }; + let at = |x: f64| if near.is_u() { from_p3(surf.subs(x, station)) } else { from_p3(surf.subs(station, x)) }; + let q: Vec = (0..5).map(|k| at(c0 + sgn * k as f64 * h)).collect(); + let (Some(&q0), Some(&q1), Some(&q2), Some(&q3), Some(&q4)) = (q.first(), q.get(1), q.get(2), q.get(3), q.get(4)) else { + return Err(KernelError::Failed("match: seam samples".into())); + }; + let d1 = (q0 * -25.0 + q1 * 48.0 + q2 * -36.0 + q3 * 16.0 + q4 * -3.0) * (sgn / (12.0 * h)); + let d2 = (q0 * 35.0 + q1 * -104.0 + q2 * 114.0 + q3 * -56.0 + q4 * 11.0) * (1.0 / (12.0 * h * h)); + let d3 = (q0 * -5.0 + q1 * 18.0 + q2 * -24.0 + q3 * 14.0 + q4 * -3.0) * (sgn / (2.0 * h * h * h)); + let (t, bend, change) = shape(d1, d2, d3)?; + Ok((q0, t, bend, change)) +} + +/// How far the rebuilt surface strays from the original one: the root mean square of the +/// distances from a grid of its samples to the original surface. +fn deviation(made: &mt::Face, original: &mt::Face) -> Result { + let (surf, os) = (made.surface(), original.surface()); + let ranges = ranges(made)?; + let mut sum = 0.0f64; + let mut n = 0usize; + for i in 0..GRID { + for j in 0..GRID { + let u = ranges[0].0 + (ranges[0].1 - ranges[0].0) * i as f64 / (GRID - 1) as f64; + let v = ranges[1].0 + (ranges[1].1 - ranges[1].0) * j as f64 / (GRID - 1) as f64; + let q = from_p3(surf.subs(u, v)); + let Some((a, b)) = os.search_nearest_parameter(p3(q), None, 100) else { continue }; + let d = from_p3(os.subs(a, b)).dist(q); + if d.is_finite() { + sum += d * d; + n += 1; + } + } + } + if n == 0 { + return Err(KernelError::Failed("match: the original surface cannot be sampled".into())); + } + Ok((sum / n as f64).sqrt()) +} + +#[cfg(test)] +#[path = "matchsurf_tests.rs"] +mod tests; diff --git a/crates/kernel/src/matchsurf_tests.rs b/crates/kernel/src/matchsurf_tests.rs new file mode 100644 index 00000000..f99066be --- /dev/null +++ b/crates/kernel/src/matchsurf_tests.rs @@ -0,0 +1,225 @@ +//! Tests of surface matching: a sheet is rebuilt to meet a cylinder's wall up to G1, G2 and G3, +//! its seam, its far boundary and its interior are measured, and what cannot be matched is +//! refused. + +use super::*; +use crate::bridge::{Continuity, Frame}; +use crate::build::SEAM_ANGLE; +use crate::{cylinder, stitch}; +use solidcraft_geom::Vec3; +use truck_modeling as mt; + +/// The radius of the cylinder a sheet is matched to, the radius the sheet starts at, its height, +/// and how far round the cylinder the sheet runs. +const R_TARGET: f64 = 10.0; +const R_SHEET: f64 = 11.0; +const HEIGHT: f64 = 30.0; +const SPAN: f64 = 0.6; + +/// A point on the vertical line at radius `r`, angle `a`, height `z` (Z up). +fn at(a: f64, r: f64, z: f64) -> Vec3 { + Vec3::new(r * a.cos(), r * a.sin(), z) +} + +/// The outward radial and the circumferential unit directions at the angle `a`. +fn radial(a: f64) -> Vec3 { + Vec3::new(a.cos(), a.sin(), 0.0) +} + +fn circ(a: f64) -> Vec3 { + Vec3::new(-a.sin(), a.cos(), 0.0) +} + +/// A cylindrical sheet of radius `r` between the angles `a0` and `a1`, `stations` high: the body +/// a match rebuilds. Its near boundary is the vertical line at `a0`, its far boundary at `a1`. +fn cylinder_sheet(r: f64, a0: f64, a1: f64, stations: usize) -> Body { + let mut a = Vec::new(); + let mut b = Vec::new(); + for i in 0..stations { + let z = HEIGHT * i as f64 / (stations - 1) as f64; + a.push(Frame { p: at(a0, r, z), leave: circ(a0), bend: radial(a0) * (-1.0 / r), bend_change: Vec3::ZERO }); + b.push(Frame { p: at(a1, r, z), leave: circ(a1) * -1.0, bend: radial(a1) * (-1.0 / r), bend_change: Vec3::ZERO }); + } + bridge::bridge(&a, &b, [Continuity::G2, Continuity::G2], [0.0, 0.0]).unwrap() +} + +/// The stock cylinder, and two picks on the vertical seam line of its wall. +fn cylinder_and_picks() -> (Body, [Vec3; 2]) { + let c = cylinder(Vec3::ZERO, Vec3::new(0.0, 0.0, 1.0), R_TARGET, HEIGHT).unwrap(); + (c, [at(SEAM_ANGLE, R_TARGET, 6.0), at(SEAM_ANGLE, R_TARGET, 24.0)]) +} + +/// The sheet of the usual shape, and a pick on its near boundary. +fn sheet_and_pick() -> (Body, [Vec3; 1]) { + (cylinder_sheet(R_SHEET, SEAM_ANGLE, SEAM_ANGLE + SPAN, 9), [at(SEAM_ANGLE, R_SHEET, HEIGHT * 0.5)]) +} + +/// The middle of the edge of `body`'s first face nearest `p`. +fn edge_at(body: &Body, p: Vec3) -> mt::Edge { + let face = body.solid.face_iter().next().unwrap(); + let mut best: Option<(f64, mt::Edge)> = None; + for e in face.edge_iter() { + let c = e.curve(); + let (t0, t1) = c.range_tuple(); + let m = from_p3(c.subs(0.5 * (t0 + t1))); + let d = m.dist(p); + if best.as_ref().is_none_or(|(b, _)| d < *b) { + best = Some((d, e.clone())); + } + } + best.unwrap().1 +} + +/// How far the rebuilt face's far boundary lies from the original one, in mm: the far row of the +/// rebuilt surface sampled at its own stations, against the original boundary edge's own points. +fn far_drift(made: &Body, sheet: &Body, far_mid: Vec3) -> f64 { + let face = made.solid.face_iter().next().unwrap(); + let surf = face.surface(); + let (Some((u0, u1)), Some((_, v1))) = surf.try_range_tuple() else { panic!("unbounded") }; + let edge = edge_at(sheet, far_mid); + let poly = chain_polyline(&[edge]).unwrap(); + let mut worst = 0.0f64; + for k in 0..=8 { + let p = from_p3(surf.subs(u0 + (u1 - u0) * k as f64 / 8.0, v1)); + worst = worst.max(chain_distance(&poly, p)); + } + worst +} + +/// What the rebuilt sheet does at the normalized station `t` of its matching boundary: the +/// point, the unit tangent, the curvature vector and the change of that curvature, read off the +/// finished surface (never from the construction). +fn seam_at(body: &Body, t: f64) -> (Vec3, Vec3, Vec3, Vec3) { + let face = body.solid.face_iter().next().unwrap(); + let ranges = super::ranges(face).unwrap(); + super::seam_shape(&face.surface(), &ranges, Side::V0, t).unwrap() +} + +/// The stations a seam is measured at, as fractions of its length. +fn stations() -> Vec { + (0..9).map(|k| k as f64 / 8.0).collect() +} + +#[test] +fn a_g2_match_meets_a_cylinder_with_its_own_curvature() { + let (target, picks) = cylinder_and_picks(); + let (sheet, pick) = sheet_and_pick(); + let made = match_surface(&sheet, &pick, &target, &picks, Continuity::G2).unwrap(); + // The rebuilt boundary has moved onto the picked line: the seam is on the cylinder's wall. + assert!(seam_at(&made, 0.5).0.dist(at(SEAM_ANGLE, R_TARGET, HEIGHT * 0.5)) < 1e-6, "the seam is on the picked line"); + let (gap, deg, order) = seam_continuity(&made, &target, &picks).unwrap(); + println!("g2: gap {gap:e} mm, normal {deg:e} deg, order {order:?}"); + assert!(gap < 1e-6, "position gap {gap:e}"); + assert!(deg < 0.01, "normal angle {deg}"); + assert!(order.conditions() >= Continuity::G2.conditions(), "{order:?}"); + // The curvature the cylinder bends with along the seam: matched to a thousandth of itself. + let want = radial(SEAM_ANGLE) * (-1.0 / R_TARGET); + let mut worst = 0.0f64; + for t in stations() { + worst = worst.max(seam_at(&made, t).2.dist(want) * R_TARGET); + } + println!("g2: worst curvature jump {worst:e} of the cylinder's own"); + assert!(worst < 1e-3, "curvature jump {worst:e}"); + // The far boundary did not move: it is built through the original boundary's own points and + // stays on the original edge between them to round-off. + let far = at(SEAM_ANGLE + SPAN, R_SHEET, HEIGHT * 0.5); + let drift = far_drift(&made, &sheet, far); + println!("g2: far boundary drift {drift:e} mm"); + assert!(drift < 1e-9, "the far boundary drifted {drift:e} mm"); + // The interior stays near the original: the crowning that leaves it closest. + let (face, original) = (made.solid.face_iter().next().unwrap().clone(), sheet.solid.face_iter().next().unwrap().clone()); + let stray = deviation(&face, &original).unwrap(); + println!("g2: interior stray {stray:.6} mm rms"); + // The seam moved 1 mm, so the face cannot be closer than that at the seam itself: over the + // whole face the root mean square distance is well under half of it. + assert!(stray < 0.6, "the interior strays {stray} mm"); +} + +#[test] +fn a_g1_match_holds_the_tangent_and_leaves_the_curvature_free() { + let (target, picks) = cylinder_and_picks(); + let (sheet, pick) = sheet_and_pick(); + let made = match_surface(&sheet, &pick, &target, &picks, Continuity::G1).unwrap(); + let (gap, deg, order) = seam_continuity(&made, &target, &picks).unwrap(); + let want = radial(SEAM_ANGLE) * (-1.0 / R_TARGET); + let worst = stations().into_iter().map(|t| seam_at(&made, t).2.dist(want) * R_TARGET).fold(0.0f64, f64::max); + println!("g1: gap {gap:e} mm, normal {deg:e} deg, order {order:?}, curvature jump {worst:e}"); + assert!(gap < 1e-6 && deg < 0.01, "position {gap:e}, normal {deg}"); + assert_eq!(order, Continuity::G1, "a tangent match leaves the curvature free"); + assert!(worst > 1e-3, "a G1 match need not hold the curvature: {worst:e}"); + // The far boundary is kept whatever the order. + let far = at(SEAM_ANGLE + SPAN, R_SHEET, HEIGHT * 0.5); + assert!(far_drift(&made, &sheet, far) < 1e-9); +} + +#[test] +fn a_g3_match_holds_the_change_of_curvature() { + let (target, picks) = cylinder_and_picks(); + let (sheet, pick) = sheet_and_pick(); + let made = match_surface(&sheet, &pick, &target, &picks, Continuity::G3).unwrap(); + let (gap, deg, order) = seam_continuity(&made, &target, &picks).unwrap(); + let want = radial(SEAM_ANGLE) * (-1.0 / R_TARGET); + let (mut bend, mut change) = (0.0f64, 0.0f64); + for t in stations() { + let (_, _, b, c) = seam_at(&made, t); + bend = bend.max(b.dist(want) * R_TARGET); + // A cylinder's curvature does not change in its own surface: the normal part of the + // change of curvature across the seam must be zero too. + change = change.max(c.dot(radial(SEAM_ANGLE)).abs() * R_TARGET * R_TARGET); + } + println!("g3: gap {gap:e} mm, normal {deg:e} deg, order {order:?}, curvature jump {bend:e}, change {change:e}"); + assert!(gap < 1e-6 && deg < 0.01, "position {gap:e}, normal {deg}"); + assert_eq!(order, Continuity::G3, "the seam holds the change of curvature"); + assert!(bend < 1e-3, "curvature jump {bend:e}"); + assert!(change < 1e-3, "change of curvature {change:e}"); + let far = at(SEAM_ANGLE + SPAN, R_SHEET, HEIGHT * 0.5); + assert!(far_drift(&made, &sheet, far) < 1e-9); +} + +#[test] +fn a_match_refuses_what_it_cannot_hold() { + let (target, picks) = cylinder_and_picks(); + let (sheet, pick) = sheet_and_pick(); + // No pick on the body's face, no pick on the target's edges, and picks that are not numbers. + assert!(match_surface(&sheet, &[], &target, &picks, Continuity::G1).is_err()); + assert!(match_surface(&sheet, &[Vec3::new(0.0, 0.0, 100.0)], &target, &picks, Continuity::G1).is_err()); + assert!(match_surface(&sheet, &[Vec3::new(f64::NAN, 0.0, 0.0)], &target, &picks, Continuity::G1).is_err()); + assert!(match_surface(&sheet, &vec![Vec3::ZERO; MAX_PICKS + 1], &target, &picks, Continuity::G1).is_err()); + assert!(match_surface(&sheet, &pick, &target, &[], Continuity::G1).is_err()); + // A pick that is not on the target at all, and one in the middle of its wall: no chain there. + assert!(match_surface(&sheet, &pick, &target, &[Vec3::new(100.0, 0.0, 15.0)], Continuity::G1).is_err()); + assert!(match_surface(&sheet, &pick, &target, &[Vec3::new(10.0, 0.0, 15.0)], Continuity::G1).is_err()); + // A body that is not a surface: the stock cylinder itself. + assert!(match_surface(&target, &[at(SEAM_ANGLE, R_TARGET, 15.0)], &target, &picks, Continuity::G1).is_err()); +} + +#[test] +fn a_target_that_jumps_in_curvature_supplies_only_g1() { + // Two sheets of the target meet along a line with the same tangent plane but different + // curvature: a cylinder's wall and the plane tangent to it there. The target itself does not + // hold the curvature along that line, so no match may be told to. + let wall = cylinder_sheet(R_TARGET, SEAM_ANGLE, SEAM_ANGLE + 0.5, 5); + let mut a = Vec::new(); + let mut b = Vec::new(); + for i in 0..5 { + let z = HEIGHT * i as f64 / 4.0; + let (p, q) = (at(SEAM_ANGLE, R_TARGET, z), at(SEAM_ANGLE, R_TARGET, z) + circ(SEAM_ANGLE) * 8.0); + a.push(Frame { p, leave: circ(SEAM_ANGLE), bend: Vec3::ZERO, bend_change: Vec3::ZERO }); + b.push(Frame { p: q, leave: circ(SEAM_ANGLE) * -1.0, bend: Vec3::ZERO, bend_change: Vec3::ZERO }); + } + let plane = bridge::bridge(&a, &b, [Continuity::G2, Continuity::G2], [0.0, 0.0]).unwrap(); + let joined = stitch(&[&wall, &plane], 1e-6).unwrap(); + assert_eq!(joined.solid.face_iter().count(), 2, "the two sheets sewed into one body"); + let picks = [at(SEAM_ANGLE, R_TARGET, 6.0), at(SEAM_ANGLE, R_TARGET, 24.0)]; + let (sheet, pick) = sheet_and_pick(); + // The seam is smooth enough to match, so G1 goes through; the curvature jumps there, so G2 + // and G3 do not. + assert!(match_surface(&sheet, &pick, &joined, &picks, Continuity::G1).is_ok()); + let err = match_surface(&sheet, &pick, &joined, &picks, Continuity::G2).unwrap_err().to_string(); + println!("refused g2: {err}"); + assert!(err.contains("only g1"), "{err}"); + assert!(match_surface(&sheet, &pick, &joined, &picks, Continuity::G3).is_err()); + // The plain cylinder, whose two wall faces are the same surface, supplies G3. + let (target, plain_picks) = cylinder_and_picks(); + assert!(match_surface(&sheet, &pick, &target, &plain_picks, Continuity::G3).is_ok()); +} diff --git a/crates/sketch/src/model.rs b/crates/sketch/src/model.rs index 493df7b0..65581877 100644 --- a/crates/sketch/src/model.rs +++ b/crates/sketch/src/model.rs @@ -1289,8 +1289,8 @@ impl Sketch { } /// Exact boundary segments of a free-form curve: ellipses as four conics, conics, fit - /// splines and control splines up to degree 3 as Béziers. None for other curves (and for - /// higher-degree splines, which stay polylines). + /// splines and control splines up to `MAX_BEZIER_DEGREE` as Béziers. None for other curves + /// (and for higher-degree splines, which stay polylines). pub fn exact_segs(&self, ci: usize) -> Option> { let cu = self.curves.get(ci)?; let pt = |i: &usize| self.point(*i); @@ -1326,6 +1326,17 @@ impl Sketch { } } + /// The exact [`Seg2`] pieces of a curve of any kind, in the curve's own direction of travel: + /// lines, circles and arcs from their analytic form, free-form curves through + /// [`Sketch::exact_segs`]. Empty for a curve with no exact pieces (a spline above + /// `solidcraft_geom::MAX_BEZIER_DEGREE`). + pub fn exact_pieces(&self, ci: usize) -> Vec { + match self.exact_segs(ci) { + Some(segs) => segs, + None => self.segs(ci), + } + } + /// Polyline of a free-form curve (ellipse, spline, conic); other curves: their segments' /// polylines. pub fn polyline(&self, ci: usize) -> Vec { @@ -1375,8 +1386,11 @@ impl Sketch { "e" } CurveKind::Spline { pts, control, degree, periodic } => { - if pts.len() < 2 || pts.len() > crate::curves::MAX_SPLINE_POINTS || !(1..=7).contains(degree) { - return Err(SketchError::Invalid("a spline needs 2 to 500 points and degree 1 to 7".into())); + if pts.len() < 2 || pts.len() > crate::curves::MAX_SPLINE_POINTS || !(1..=solidcraft_geom::MAX_BEZIER_DEGREE as u8).contains(degree) { + return Err(SketchError::Invalid(format!( + "a spline needs 2 to 500 points and degree 1 to {}", + solidcraft_geom::MAX_BEZIER_DEGREE + ))); } if *periodic { if *control || *degree != 3 { diff --git a/docs/screenshots/class-a-boundary.png b/docs/screenshots/class-a-boundary.png new file mode 100644 index 00000000..e5f436e9 Binary files /dev/null and b/docs/screenshots/class-a-boundary.png differ