diff --git a/ATTRIBUTION.md b/ATTRIBUTION.md index 97a67f23..17be64a9 100644 --- a/ATTRIBUTION.md +++ b/ATTRIBUTION.md @@ -462,7 +462,7 @@ copied into SolidCraft. | vendor/truck-meshalgo | truck-meshalgo 0.4.0, https://github.com/ricosjp/truck (RICOS Co. Ltd.) | Apache-2.0 | meshing: bounded parameter searches, collapsed-boundary projection, subdivision and provisional fixed-grid previews (vendor/README.md) | | vendor/truck-stepio | truck-stepio 0.3.0, https://github.com/ricosjp/truck (RICOS Co. Ltd.) | Apache-2.0 | STEP writer: intersection-curve surfaces written at unique ids (vendor/README.md) | | vendor/truck-shapeops | truck-shapeops 0.4.0, https://github.com/ricosjp/truck (RICOS Co. Ltd.) | Apache-2.0 | robust meshing and nearest-parameter fallback in booleans; degenerate intersection-curve rejection; face and edge history of booleans; trace macro; corner choice at vertices several intersection curves end on (vendor/README.md) | -| vendor/truck-geometry | truck-geometry 0.5.0, https://github.com/ricosjp/truck (RICOS Co. Ltd.) | Apache-2.0 | intersection curves divided from their neighbours' surface parameters (vendor/README.md) | +| vendor/truck-geometry | truck-geometry 0.5.0, https://github.com/ricosjp/truck (RICOS Co. Ltd.) | Apache-2.0 | intersection curves divided from their neighbours' surface parameters; surface second derivatives (vendor/README.md) | | vendor/matext4cgmath | matext4cgmath 0.1.0, https://github.com/ricosjp/matext4cgmath (RICOS Co. Ltd.) | Apache-2.0 | root refinement stops after 100 Newton steps (vendor/README.md) | ## Vendored renderer compatibility fix diff --git a/ROADMAP.md b/ROADMAP.md index 58799605..852c3348 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -53,7 +53,7 @@ PCB / Package, Render, Animation, Simulation, Manufacture/CAM, Drawings. | M5 | Persistent naming of faces and edges (docs/book/src/dev/naming.md): faces named from their exact kernel-reported history (sketch curve sides, caps, revolve sections, primitive sides, blends, shell inner walls, pattern/mirror copies; split pieces told apart by the faces around them, lumps of a cut body named after their faces), edges by their faces; fillets/chamfers, sketches on faces, shell/draft/offset/replace faces, hole placements, joint origins and face appearances resolve by name first and warn when the entity split or vanished: done. sheet metal edges by base sketch line and flange side, shell inner faces, sweep and loft roles: done. kernel-reported face history: done (merged faces answer to every name; operations without history yet fall back to geometric names, flagged unverified; format 14 keeps older designs' geometric names). Left: history through sheet metal, plastic, surface and emboss features; FaceSource certificates through history | in progress | 20 | | M6 | UI depth: selection feedback and origin (done), timeline editing (done: hover and selection highlight modifier faces, including fillets and chamfers; edit feature, history marker drag, dependency-checked single and multi-item reorder, suppress, rename, delete with dependents, re-resolved references; redefining selections refreshes persistent edge/face references while retaining feature identity and dimension names; optional compact named left feature history with rollback bar, multi-selection and persistent folders above the parts browser), right-docked dialogs (done), live previews for every feature dialog (done), keyboard-first commands with on-canvas value boxes and sketch inline dimensions (done), drag manipulators with automatic join/cut (done: extrude incl. two sides, fillet, shell, hole depth, move, pattern spacing, revolve rotator; taper handle open), dark theme (done), right-click menu with Repeat (done), measure tool (done), section analysis (done), selection filter (done), drag-solve of sketch points (done); appearance controls now preserve body/face/component targeting, named transparent looks, reuse/search and override removal; preferences | in progress | 6 | | M7 | Thin extrude (walls along open or closed curves: centre, inside, outside or two sides; mitred corners, crossing walls joined; #721), sweep, loft (incl. to a point), rib and web (to next / depth), emboss (planar faces and exact radial emboss/deboss of bounded orthogonal profiles on complete exterior cylinders/sleeves), thread (cosmetic, and modelled M3–M20 outside and inside: the ISO profile screwed in as a faceted cut), coil (smooth B-spline helical walls), pipe, pattern on path, patterns of bodies and components, replace face (parallel), align, remove, feature copy/paste: done. Unique eligible cylindrical targets are inferred for new Emboss features. Left: general curved-face emboss and exact (unfaceted) modelled threads | in progress | 12 | -| M8 | Surface workspace (#772). Done: Patch (sketch regions; edge loops filled G0, plus `crown`), Stitch, Thicken (planes and cylinders), Extend (straight edges of planar surfaces), Trim (by a plane or a face), Delete Face (healing), Boundary Fill, zebra/draft/curvature-map analysis, Surface Extrude/Revolve/Sweep/Loft from open and closed curves, Offset, Unstitch and Reverse Normal; the continuity-driven core (`kernel/src/bridge.rs`, G0–G3 Hermite rows) with **Blend** (`surface.blend`, #764), **Match** (`surface.match`), **Boundary** (`surface.boundary`, #763 — two boundaries exact, three/four refused pending a Gregory-type corner), **Fair** (`surface.fair`, CV fairing), **continuity analysis** (`inspect.continuity`, #770; a toolbar tool with edge markers and a legend, and curvature combs along face isocurves) and a multi-section **loft** with per-end continuity and crown. Also in the Class-A pass: G1…G4 blend curves in the sketch, chord and variable-radius fillets on curved edges (with honest refusals where the fit cannot hold), and a G2 corner patch where filleted edges meet (kernel; the command refuses until the patch's sections can be sewn to the blends' contacts). Left: bias/form-factor fillet sections (design recorded), a Gregory-type corner, surface Fillet between separate surface bodies (#769), Ruled (#762), Untrim (#761), curved-face Extend (#760), loft guide rails (#201; sweep guide rails and two-rail sweeps are done, #771), and every remaining #757–#772 item | deferred | — | +| M8 | Surface workspace (#772). Done: Patch (sketch regions; edge loops filled G0, plus `crown`), Stitch, Thicken (planes and cylinders), Extend (straight edges of planar surfaces), Trim (by a plane or a face), Delete Face (healing), Boundary Fill, zebra/draft/curvature-map analysis, Surface Extrude/Revolve/Sweep/Loft from open and closed curves, Offset, Unstitch and Reverse Normal; the continuity-driven core (`kernel/src/bridge.rs`, G0–G3 Hermite rows) with **Blend** (`surface.blend`, #764), **Match** (`surface.match`), **Boundary** (`surface.boundary`, #763 — two sketch curves exact with G0–G3, three or four as the flat patch of their loop, and curve networks: Coons and Gordon surfaces through sketch, 3D-sketch or edge curves, a point for a three-sided side, G1/G2 against the faces along one direction's edges, joined into their body on request; G1/G2 on two sides that meet at a corner waits for the Gregory-type corner), **Fair** (`surface.fair`, CV fairing), **continuity analysis** (`inspect.continuity`, #770; a toolbar tool with edge markers and a legend, and curvature combs along face isocurves) and a multi-section **loft** with per-end continuity and crown. Also in the Class-A pass: G1…G4 blend curves in the sketch, chord and variable-radius fillets on curved edges (with honest refusals where the fit cannot hold), and a G2 corner patch where filleted edges meet (kernel; the command refuses until the patch's sections can be sewn to the blends' contacts). Left: bias/form-factor fillet sections (design recorded), a Gregory-type corner, surface Fillet between separate surface bodies (#769), Ruled (#762), Untrim (#761), curved-face Extend (#760), loft guide rails (#201; sweep guide rails and two-rail sweeps are done, #771), and every remaining #757–#772 item | deferred | — | | M9 | Components and assemblies: New Component / Create Components from Bodies, active component, per-component origin/bodies/sketches/timeline entries, occurrences with transforms, Copy/Paste vs Paste New, nested Browser tree with visibility and isolate; Ground, Move/Copy of occurrences, joints (rigid, revolute, slider, cylindrical, pin-slot, planar, ball) by joint origins with limits and Drive Joints, closed kinematic loops (four-bar, slider-crank: the loop's other joints follow the driven one on its branch, and driving stops where a loop can no longer close), cams (`joint.cam`: a revolute cam drives a slider or revolute follower along a closed sketch profile or the cam's edges, knife, roller or flat followers, worked out again when the profile changes; drives, studies, closed loops, contact sets), As-built Joint, Rigid Group, Contact Sets (driving and motion studies stop where bodies in a set would pass through each other), Motion Study (joint values keyed along a timeline, playback, export of occurrence positions as JSON/CSV), motion analysis (joint values and a point's path charted over a study as value, velocity and acceleration; the parts that collide along a study, with the steps where they do), Exploded views (automatic or by hand, shown without changing placements); interference, per-component physical properties; STEP with assembly structure; URDF export for robot simulators (links with meshes, mass and inertia; joint chains whose values match Drive Joints; docs/urdf.md). Physical properties now measure material-aware component/occurrence placements, descendants and repeated instances, including pending moves and body reassignment offsets; contact geometry and automatic exploded centres respect those offsets. Groundwork done: component tree, occurrences with transforms (instances, Paste New), world placement, active component as authoring frame, ground, Capture/Revert Position, moving bodies and sketches. UI: ASSEMBLE panel dialogs for New Component, Joint (snaps on faces, edges, circles and vertices; animated motion preview; offsets, flip, limits), As-built Joint, Joint Origin, Rigid Group, Drive Joints (live slider), Motion Link, Cam and Interference (overlaps shown red); joints in the timeline and the Browser | in progress | 32 | | M10 | Sheet metal: flange (base, edge, contour), hem, unfold/refold, convert, flat pattern + DXF, rules: done. PLASTIC enclosure subset: Boss, Rib, Web, Lip/Groove, Snap Fit, Rest, Plastic Rules (library, edit, assign; rule draft and clearance as feature defaults): done, tested against hand-computed geometry. UI: SHEET METAL and PLASTIC tabs in Fusion's layout with dialogs (live preview, Edit Feature, pre-selection) for Flange, Hem, Unfold/Refold, Convert, Sheet Metal Rules, Create Flat Pattern (flat view with bend lines, bend table CSV/JSON and DXF export with exact circle/arc/line cutouts), Boss, Lip/Groove, Snap Fit, Rest, Manage and Assign Plastic Rules; Fold (bend along a sketch line); flanges meeting at a corner closed at any angle with the rule's gap; edge flanges with their own angle and height per edge in one feature; corner reliefs (notch, round, rectangular, tear) chosen by the rule; Convert to Sheet Metal from a solid or shelled body by its bend edges (flat faces, sharp bends) | in progress | 36 | | M11 | 2D drawings | deferred | — | diff --git a/apps/solidcraft-cli/src/recipe.rs b/apps/solidcraft-cli/src/recipe.rs index 68d5bfdd..d5a810b0 100644 --- a/apps/solidcraft-cli/src/recipe.rs +++ b/apps/solidcraft-cli/src/recipe.rs @@ -423,6 +423,24 @@ pub fn to_script(recipe: &Value) -> Result { "sketch": profile_sketch(f), "profiles": profiles_of(f), "path_sketch": path.get("sketch"), "path": path.get("curves"), "operation": op(f), "name": name, "body_names": bodies}})); } + // A boundary surface: a loop of one sketch's curves (`sketch`, `curves`), or a network + // (`u`, `v`), as `surface.boundary` takes them. + Some("boundary_surface") => { + let mut p = json!({"name": name}); + for k in ["sketch", "curves", "continuity", "u", "v", "operation"] { + if let Some(v) = f.get(k) { + p[k] = v.clone(); + } + } + if let Some(b) = bodies.as_array().and_then(|a| a.first()) { + p["body_name"] = b.clone(); + } + out.push(json!({"command": "surface.boundary", "params": p})); + } + Some("thicken") => out.push(json!({"command": "surface.thicken", "params": { + "bodies": f.get("bodies"), "thickness": f.get("thickness"), "symmetric": f.get("symmetric").cloned().unwrap_or(json!(false)), + "operation": op(f), "name": name}})), + Some("remove") => out.push(json!({"command": "solid.remove", "params": {"bodies": f.get("bodies"), "name": name}})), Some("shell") => { let faces: Vec = f.get("faces_removed").and_then(Value::as_array).into_iter().flatten().filter_map(|x| x.get("point").cloned()).collect(); out.push(json!({"command": "solid.shell", "params": {"faces": faces, "thickness": f.get("inside_thickness"), "name": name}})); diff --git a/assets/locales/en.json b/assets/locales/en.json index 623a0f9e..2bb067ed 100644 --- a/assets/locales/en.json +++ b/assets/locales/en.json @@ -31,6 +31,7 @@ "Arc length": "Arc length", "Area": "Area", "Attachment": "Attachment", + "BOUNDARY SURFACE": "BOUNDARY SURFACE", "BROWSER": "BROWSER", "Bend Line Position": "Bend Line Position", "Bend Position": "Bend Position", @@ -73,6 +74,7 @@ "Counterbore depth": "Counterbore depth", "Countersink angle": "Countersink angle", "Create Issue — review draft": "Create Issue — review draft", + "Cross curves": "Cross curves", "Cube": "Cube", "Current orthographic view (1:1 mm)": "Current orthographic view (1:1 mm)", "Current perspective view": "Current perspective view", diff --git a/assets/locales/ja.json b/assets/locales/ja.json index 34c12b36..62a43a25 100644 --- a/assets/locales/ja.json +++ b/assets/locales/ja.json @@ -31,6 +31,7 @@ "Arc length": "円弧の長さ", "Area": "表面積", "Attachment": "アタッチメント", + "BOUNDARY SURFACE": "境界サーフェス", "BROWSER": "ブラウザ", "Bend Line Position": "曲げ線位置", "Bend Position": "曲げ位置", @@ -73,6 +74,7 @@ "Counterbore depth": "座ぐり深さ", "Countersink angle": "皿穴角度", "Create Issue — review draft": "Issueを作成 — 下書きの確認", + "Cross curves": "交差する曲線", "Cube": "キューブ", "Current orthographic view (1:1 mm)": "現在の平行投影（実寸・mm）", "Current perspective view": "現在の透視投影", diff --git a/assets/locales/zh.json b/assets/locales/zh.json index 8028d7ef..2a19478f 100644 --- a/assets/locales/zh.json +++ b/assets/locales/zh.json @@ -31,6 +31,7 @@ "Arc length": "弧长", "Area": "面积", "Attachment": "附着", + "BOUNDARY SURFACE": "边界曲面", "BROWSER": "浏览器", "Bend Line Position": "弯折线位置", "Bend Position": "折弯位置", @@ -73,6 +74,7 @@ "Counterbore depth": "沉孔深度", "Countersink angle": "沉头角度", "Create Issue — review draft": "创建 Issue——审阅草稿", + "Cross curves": "交叉曲线", "Cube": "立方体", "Current orthographic view (1:1 mm)": "当前正交视图（1:1毫米）", "Current perspective view": "当前透视视图", diff --git a/crates/doc/src/document.rs b/crates/doc/src/document.rs index 426a9c20..2b2f6753 100644 --- a/crates/doc/src/document.rs +++ b/crates/doc/src/document.rs @@ -1465,6 +1465,19 @@ pub enum FeatureKind { #[serde(default = "boundary_samples")] samples: usize, }, + /// One surface through a curve network (`surface.boundary` with `u` and `v`, issue #763): the + /// curves `u` run one way (two or more: the surface lofts through them) and `v` cross them the + /// other way (two each way: a Coons patch; more: a Gordon surface through every curve). A side + /// may be a point (a three-sided patch). A curve that is a chain of a surface's open edges may + /// carry G1 or G2 with the face beside it. `operation` join sews the surface into the body of + /// its first edge curve (a solid when it closes the gap there). + SurfaceNetwork { + u: Vec, + #[serde(default, skip_serializing_if = "Vec::is_empty")] + v: Vec, + #[serde(default)] + operation: Operation, + }, /// 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 @@ -1692,6 +1705,27 @@ pub struct BlendSide { pub edges: Vec, } +/// One curve of a boundary surface's network ([`FeatureKind::SurfaceNetwork`]): a chain of a +/// sketch's `curves` (planar or 3D, meeting end to end), a chain of `body`'s edges picked by points +/// on them (`edges`), or a side collapsed to the point `point` (the apex of a three-sided patch). +/// A chain of a surface's open edges may carry `continuity` with the face beside it (G0, the +/// default, G1 or G2). +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize, solidcraft_derive::Walk)] +pub struct NetworkCurve { + #[serde(default, skip_serializing_if = "Option::is_none")] + pub sketch: Option, + #[serde(default, skip_serializing_if = "Vec::is_empty")] + pub curves: Vec, + #[serde(default, skip_serializing_if = "Option::is_none")] + pub body: Option, + #[serde(default, skip_serializing_if = "Vec::is_empty")] + pub edges: Vec, + #[serde(default, skip_serializing_if = "Option::is_none")] + pub point: Option, + #[serde(default = "patch_g0", skip_serializing_if = "is_g0")] + pub continuity: Continuity, +} + /// How a patch on a loop of edges meets the faces beside the loop: every edge with `all`, but the /// edges picked in `edges` with their own. #[derive(Clone, Debug, PartialEq, Serialize, Deserialize, solidcraft_derive::Walk)] @@ -1733,9 +1767,12 @@ 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 + BOUNDARY_SAMPLES } +/// The stations a boundary patch carries when a design does not say. +pub const BOUNDARY_SAMPLES: usize = 25; + fn blend_orders() -> [Continuity; 2] { [Continuity::G1, Continuity::G1] } @@ -2088,7 +2125,7 @@ impl FeatureKind { FeatureKind::SurfaceSweep { .. } => "SurfaceSweepFeature", FeatureKind::SurfaceLoft { .. } => "SurfaceLoftFeature", FeatureKind::SurfaceBlend { .. } => "BlendFeature", - FeatureKind::SurfaceBoundary { .. } => "BoundaryFeature", + FeatureKind::SurfaceBoundary { .. } | FeatureKind::SurfaceNetwork { .. } => "BoundaryFeature", FeatureKind::SurfaceMatch { .. } => "MatchFeature", FeatureKind::SurfaceFair { .. } => "FairFeature", FeatureKind::Stitch { .. } => "StitchFeature", @@ -2168,7 +2205,7 @@ impl FeatureKind { FeatureKind::SurfaceSweep { .. } => "SurfaceSweep", FeatureKind::SurfaceLoft { .. } => "SurfaceLoft", FeatureKind::SurfaceBlend { .. } => "Blend", - FeatureKind::SurfaceBoundary { .. } => "Boundary", + FeatureKind::SurfaceBoundary { .. } | FeatureKind::SurfaceNetwork { .. } => "Boundary", FeatureKind::SurfaceMatch { .. } => "Match", FeatureKind::SurfaceFair { .. } => "Fair", FeatureKind::Stitch { .. } => "Stitch", diff --git a/crates/doc/src/eval.rs b/crates/doc/src/eval.rs index bae42e3b..61f8df0c 100644 --- a/crates/doc/src/eval.rs +++ b/crates/doc/src/eval.rs @@ -1035,6 +1035,20 @@ fn placement_digest(doc: &Document, f: &Feature) -> u64 { } return hash.finish(); } + // A network through 3D sketch curves reads their points through the components the + // feature and those sketches are placed in, as a sweep's wired path does. + if let FeatureKind::SurfaceNetwork { u, v, .. } = &f.kind { + let sketches: Vec = u.iter().chain(v).filter_map(|c| c.sketch).collect(); + if sketches.iter().any(wired) { + let mut hash = std::collections::hash_map::DefaultHasher::new(); + doc.component_transform(f.component).iter().flatten().for_each(|value| value.to_bits().hash(&mut hash)); + for sketch in &sketches { + let owner = doc.feature(sketch.0).map(|feature| feature.component).unwrap_or(0); + doc.component_transform(owner).iter().flatten().for_each(|value| value.to_bits().hash(&mut hash)); + } + return hash.finish(); + } + } // To / From Object finds its target where that target's component is placed (#928). let to_object = match &f.kind { FeatureKind::Extrude { extent, .. } => extent.to.is_some() || extent.from.is_some(), @@ -4691,6 +4705,7 @@ fn eval_feature(doc: &Document, vals: &BTreeMap, f: &Feature, st: | FeatureKind::SurfaceLoft { .. } | FeatureKind::SurfaceBlend { .. } | FeatureKind::SurfaceBoundary { .. } + | FeatureKind::SurfaceNetwork { .. } | FeatureKind::SurfaceMatch { .. } | FeatureKind::SurfaceFair { .. } | FeatureKind::SurfaceOffset { .. } diff --git a/crates/doc/src/eval_surface.rs b/crates/doc/src/eval_surface.rs index 1464ecda..227a126a 100644 --- a/crates/doc/src/eval_surface.rs +++ b/crates/doc/src/eval_surface.rs @@ -151,12 +151,49 @@ pub(super) fn eval(doc: &Document, vals: &BTreeMap, f: &Feature, 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)?; + // Two curves are the two sides of a bridge. Three or four run round a loop: the Coons + // patch of the loop, which lies in the sketch's plane with them, so every order holds. + let body = if curves.len() >= 3 { + if *crown != 0.0 { + return Err(DocError::Invalid( + "not supported yet: crown on a patch through three or four curves (it is the flat patch of its sides)".into(), + )); + } + loop_patch(&ss, curves)? + } else { + let boundaries = boundary_frames(&ss, curves, orders, *samples)?; + kernel::boundary_patch(&boundaries, *crown, *samples)? + }; let name = new_name(st, f, 0); st.bodies.push(ModelBody::new(name, body, f.id)); Ok(()) } + FeatureKind::SurfaceNetwork { u, v, operation } => { + if !matches!(operation, crate::Operation::NewBody | crate::Operation::Join) { + return Err(DocError::Invalid("a boundary surface is a new body or joins the body of its edges".into())); + } + let us = network_sources(doc, st, f, u, "U")?; + let vs = network_sources(doc, st, f, v, "V")?; + let made = kernel::boundary_surface(&network_curves(&us), &network_curves(&vs))?; + if *operation == crate::Operation::Join { + // Into the body of the first curve picked on a body's edges. + let target = u + .iter() + .chain(v) + .find_map(|c| c.body.clone()) + .ok_or_else(|| DocError::Invalid("join needs a curve on a body's edges: the surface is sewn into that body".into()))?; + let i = st.bodies.iter().position(|b| b.name == target).ok_or_else(|| DocError::Unknown(format!("body `{target}`")))?; + let Some(mb) = st.bodies.get(i).cloned() else { return Err(DocError::Invalid("body".into())) }; + let joined = kernel::boundary_join(&mb.body, &made)?; + if let Some(slot) = st.bodies.get_mut(i) { + *slot = ModelBody::new(mb.name, joined, mb.feature); + } + return Ok(()); + } + let name = new_name(st, f, 0); + st.bodies.push(ModelBody::new(name, made, 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())) }; @@ -439,6 +476,135 @@ fn push_sheet(st: &mut ModelState, f: &Feature, traced: crate::lineage::Traced) st.bodies.push(ModelBody::traced_new(name, traced, f.id)); } +/// The most curves one direction of a network may list (the kernel's own cap). +const MAX_NETWORK_CURVES: usize = solidcraft_geom::network::MAX_NETWORK_CURVES; + +/// A loop of three or four sketch curves as one patch: four give the Coons patch of opposite pairs +/// (the first and third one way, the second and fourth the other); three collapse the side +/// opposite the first curve to the corner where the other two meet. +fn loop_patch(ss: &SolvedSketch, curves: &[String]) -> Result { + let mut chains: 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 segs = ss.sketch.exact_segs(ci).unwrap_or_else(|| ss.sketch.segs(ci)); + if segs.is_empty() { + return Err(DocError::Invalid(format!("curve `{id}` has no geometry to build on"))); + } + chains.push(segs); + } + let side = |k: usize| -> Result> { + let (Some(id), Some(segs)) = (curves.get(k), chains.get(k)) else { return Err(DocError::Invalid("a boundary curve is missing".into())) }; + Ok(kernel::BoundaryCurve { + label: format!("curve `{id}`"), + source: kernel::CurveSource::Sketch { plane: ss.plane, segs: segs.clone() }, + continuity: kernel::Continuity::G0, + }) + }; + match curves.len() { + 4 => Ok(kernel::boundary_surface(&[side(0)?, side(2)?], &[side(1)?, side(3)?])?), + 3 => { + // The corner the second and third curves share, whichever way each was drawn. + let ends = |k: usize| -> Option<[Vec2; 2]> { + let c = chains.get(k)?; + Some([c.first()?.start(), c.last()?.end()]) + }; + let (Some(b), Some(c)) = (ends(1), ends(2)) else { return Err(DocError::Invalid("a boundary curve is missing".into())) }; + let mut best = (f64::INFINITY, Vec2::new(0.0, 0.0)); + for p in b { + for q in c { + if p.dist(q) < best.0 { + best = (p.dist(q), (p + q) * 0.5); + } + } + } + let apex = kernel::BoundaryCurve { + label: "the corner opposite the first curve".into(), + source: kernel::CurveSource::Point(ss.plane.to_world(best.1)), + continuity: kernel::Continuity::G0, + }; + Ok(kernel::boundary_surface(&[side(0)?, apex], &[side(1)?, side(2)?])?) + } + _ => Err(DocError::Invalid("a loop patch takes three or four curves".into())), + } +} + +/// A network curve resolved in the feature's frame, owning what it is built from. +enum NetSource { + Sketch(solidcraft_geom::Plane, Vec), + Points(Vec), + Edges(kernel::Body, Vec), + Point(Vec3), +} + +/// The curves of one direction of a network, resolved: each one sketch chain, 3D sketch chain, +/// edge chain or point, labelled `U1`, `U2`, … (with `dir` "U") for messages. +fn network_sources( + doc: &Document, + st: &ModelState, + f: &Feature, + list: &[crate::NetworkCurve], + dir: &str, +) -> Result> { + if list.len() > MAX_NETWORK_CURVES { + return Err(DocError::Invalid(format!("at most {MAX_NETWORK_CURVES} curves each way"))); + } + let mut out = Vec::with_capacity(list.len()); + for (k, c) in list.iter().enumerate() { + let label = format!("{dir}{}", k + 1); + let given = [c.sketch.is_some() || !c.curves.is_empty(), c.body.is_some() || !c.edges.is_empty(), c.point.is_some()]; + if given.iter().filter(|g| **g).count() != 1 { + return Err(DocError::Invalid(format!("{label} must be sketch curves, a body's edges or a point"))); + } + let source = if let Some(p) = c.point { + NetSource::Point(p.0) + } else if let Some(name) = &c.body { + if c.edges.is_empty() { + return Err(DocError::Invalid(format!("{label} needs points on the edges of {name}"))); + } + let b = st.body(name).ok_or_else(|| DocError::Unknown(format!("body `{name}`")))?; + NetSource::Edges(b.body.clone(), crate::plain(&c.edges)) + } else { + let sk = c.sketch.ok_or_else(|| DocError::Invalid(format!("{label} needs the sketch its curves are in")))?; + if c.curves.is_empty() || c.curves.len() > MAX_CHAIN_CURVES { + return Err(DocError::Invalid(format!("{label} needs 1 to {MAX_CHAIN_CURVES} curves of its sketch"))); + } + let raw = st.sketch(sk.0).ok_or_else(|| DocError::Unknown(format!("sketch {sk} (it must come earlier in the timeline)")))?; + // 3D sketch curves are followed through their points; planar ones exactly. + let mapped = path_in_frame(doc, raw, f.component)?; + match path_points(&mapped, &c.curves, Vec3::new(f64::NAN, f64::NAN, f64::NAN))? { + Some(pts) => NetSource::Points(pts), + None => { + let ss = sketch_in_feature_frame(doc, raw, f.component)?; + let mut chains = curve_chains(&ss, &c.curves)?.into_iter(); + let (Some(chain), None) = (chains.next(), chains.next()) else { + return Err(DocError::Invalid(format!("{label}: its curves must meet end to end in one chain"))); + }; + NetSource::Sketch(ss.plane, chain) + } + } + }; + out.push((label, source, c.continuity)); + } + Ok(out) +} + +/// The kernel's curves over resolved sources (borrowing their bodies). +fn network_curves(sources: &[(String, NetSource, kernel::Continuity)]) -> Vec> { + sources + .iter() + .map(|(label, src, g)| kernel::BoundaryCurve { + label: label.clone(), + source: match src { + NetSource::Sketch(plane, segs) => kernel::CurveSource::Sketch { plane: *plane, segs: segs.clone() }, + NetSource::Points(pts) => kernel::CurveSource::Points(pts.clone()), + NetSource::Edges(body, at) => kernel::CurveSource::Edges { body, at: at.clone() }, + NetSource::Point(p) => kernel::CurveSource::Point(*p), + }, + continuity: *g, + }) + .collect() +} + /// 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 diff --git a/crates/doc/src/feature_walk_tests/fixtures.rs b/crates/doc/src/feature_walk_tests/fixtures.rs index 2364cf6c..71cecf9b 100644 --- a/crates/doc/src/feature_walk_tests/fixtures.rs +++ b/crates/doc/src/feature_walk_tests/fixtures.rs @@ -164,6 +164,10 @@ pub(super) fn fixtures() -> Vec<(&'static str, Value)> { // 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_network", + json!({"type": "surface_network", "u": [{"sketch": 1, "curves": ["l1"]}, {"body": "B", "edges": [o], "continuity": "g1"}], "v": [{"sketch": 2, "curves": ["l1"]}, {"point": o}], "operation": "join"}), + ), ("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"})), diff --git a/crates/doc/src/feature_walk_tests/golden_exprs.txt b/crates/doc/src/feature_walk_tests/golden_exprs.txt index 5a887777..d8987147 100644 --- a/crates/doc/src/feature_walk_tests/golden_exprs.txt +++ b/crates/doc/src/feature_walk_tests/golden_exprs.txt @@ -342,6 +342,7 @@ surface_offset expressions: soff.distance, soff.tolerance surface_blend surface_boundary +surface_network surface_match surface_fair inputs: diff --git a/crates/doc/src/feature_walk_tests/golden_geometry.txt b/crates/doc/src/feature_walk_tests/golden_geometry.txt index 767444fa..314c4607 100644 --- a/crates/doc/src/feature_walk_tests/golden_geometry.txt +++ b/crates/doc/src/feature_walk_tests/golden_geometry.txt @@ -157,6 +157,8 @@ surface_offset 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_network + moves: u[1].edges[0][0], u[1].edges[0][1], u[1].edges[0][2], v[1].point[0], v[1].point[1], v[1].point[2] 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] diff --git a/crates/doc/src/feature_walk_tests/golden_refs.txt b/crates/doc/src/feature_walk_tests/golden_refs.txt index 0b621c6a..386158c4 100644 --- a/crates/doc/src/feature_walk_tests/golden_refs.txt +++ b/crates/doc/src/feature_walk_tests/golden_refs.txt @@ -182,6 +182,10 @@ surface_blend surface_boundary sketches: [1] remap_refs: sketch +surface_network + sketches: [1, 2] + remap_refs: u[0].sketch, v[0].sketch + remap_bodies: u[1].body surface_match remap_bodies: body surface_fair diff --git a/crates/doc/src/format/pending/boundary-surface-networks.rs b/crates/doc/src/format/pending/boundary-surface-networks.rs new file mode 100644 index 00000000..8a28c820 --- /dev/null +++ b/crates/doc/src/format/pending/boundary-surface-networks.rs @@ -0,0 +1,7 @@ +//! A `SurfaceNetwork` (Boundary Surface through a curve network, #763) saves `u` and `v`, each a +//! list of curves: a `sketch` and its `curves`, a `body` and points on its `edges`, or a `point` +//! (the apex of a three-sided patch), with an optional `continuity` (`g1` or `g2` with the face +//! beside an edge chain), and an `operation` (`new_body`, or `join` into the body of its edges). +//! A `SurfaceBoundary` of three or four curves, which older builds refused, now builds and saves +//! as it always did. Older designs need no transformation; the last release asks for an update +//! instead of meeting a feature it can't read. diff --git a/crates/engine/src/cmd/edit.rs b/crates/engine/src/cmd/edit.rs index 58e33c55..3d596cbe 100644 --- a/crates/engine/src/cmd/edit.rs +++ b/crates/engine/src/cmd/edit.rs @@ -978,6 +978,7 @@ fn redefine(s: &mut Session, p: &Value) -> Result { .chain(super::features_more::COMMANDS) .chain(super::construct::COMMANDS) .chain(super::surface::COMMANDS) + .chain(super::surface_boundary::COMMANDS) .any(|candidate| candidate.id == spec.id) || matches!( spec.id, @@ -1049,7 +1050,10 @@ fn redefine(s: &mut Session, p: &Value) -> Result { { *after = before.clone(); } - if std::mem::discriminant(&new.kind) != std::mem::discriminant(&old.kind) { + // A boundary surface's loop of sketch curves and its curve network are one command's two + // forms: an edit may turn one into the other. + let boundary = |k: &FeatureKind| matches!(k, FeatureKind::SurfaceBoundary { .. } | FeatureKind::SurfaceNetwork { .. }); + if std::mem::discriminant(&new.kind) != std::mem::discriminant(&old.kind) && !(boundary(&new.kind) && boundary(&old.kind)) { return Err(bad(cmd, format!("{} is not a {} feature", old.name, cid))); } // Optional inputs can move other inputs' positions (e.g. Distance 2 goes before Taper). diff --git a/crates/engine/src/cmd/surface_boundary.rs b/crates/engine/src/cmd/surface_boundary.rs index 79532302..985ab6c6 100644 --- a/crates/engine/src/cmd/surface_boundary.rs +++ b/crates/engine/src/cmd/surface_boundary.rs @@ -1,38 +1,56 @@ -//! 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. +//! Boundary surfaces (`surface.boundary`, issue #763): one patch through two to four sketch +//! curves in loop order, or one surface through a network of curves in two directions (sketch +//! curves, 3D sketch curves, a body's edges, or a point for a three-sided patch). The kernel +//! builds the face; this is the command's parameters and its checks, so that what it cannot hold +//! is refused before anything is built. use serde_json::Value; -use solidcraft_doc::{FeatureKind, SketchId}; +use solidcraft_doc::{BodyName, FeatureKind, NetworkCurve, Operation, Position, SketchId}; use solidcraft_kernel::Continuity; use super::CommandSpec; -use super::features::{add_feature, feature_sketch}; -use crate::params::{bad, string_list}; +use super::features::{add_feature, feature_sketch, sketch_id}; +use crate::params::{bad, str_, string_list, vec3}; use crate::{Result, Session}; /// The stations a boundary carries by default (the patch's resolution along the loop). -const DEFAULT_SAMPLES: usize = 25; +const DEFAULT_SAMPLES: usize = solidcraft_doc::BOUNDARY_SAMPLES; /// 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; +/// The most curves one direction of a network may list (the kernel's own cap). +const MAX_NETWORK_CURVES: usize = solidcraft_geom::network::MAX_NETWORK_CURVES; +/// The most sketch curves one network curve may chain, and points picking one chain of edges +/// (the kernel's own caps). +use solidcraft_kernel::{MAX_CHAIN_CURVES as MAX_CHAIN, MAX_EDGE_PICKS as MAX_PICKS}; pub static COMMANDS: &[CommandSpec] = &[CommandSpec::new("surface.boundary", "Boundary Surface", boundary) .at("SURFACE", "CREATE") .icon("surface_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)", + "Either a loop of one sketch's curves — 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; two: the patch spans between them and leaves each boundary along the tangent across it; three or four: the flat patch of the loop, three collapsing the side opposite the first curve to the corner there); 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; two curves only); samples?: 8…512 stations per boundary (default 25) — or a curve network — u: [curve] (two or more one way: the surface lofts through them), v?: [curve] (crossing every u-curve once: two each way give a Coons patch, more a Gordon surface through every curve; a curve that misses another by more than 0.01 mm is refused, naming the pair), each curve \"id\" (a curve of `sketch`) | {sketch?: id|name, curves: [ids meeting end to end]} (planar or 3D sketch curves) | {body: name, edges: [[x,y,z]] (points on a chain of its edges), continuity?: g0|g1|g2 (default g0; g1 the surface leaves tangent to the face beside a surface's open edges, g2 with its curvature too; on curves of one direction only)} | {point: [x,y,z]} (a side collapsed to a point: a three-sided patch); operation?: new|join (default new; join sews the surface into the body of its first edge curve, a solid when it closes the gap). name?, body_name? (a surface body)", )]; -/// `surface.boundary`: one patch through a network of two to four sketch curves. +/// `surface.boundary`: a loop of sketch curves, or a network when `u` is given. fn boundary(s: &mut Session, p: &Value) -> Result { + if p.get("u").is_some_and(|v| !v.is_null()) { + return network(s, p); + } + loop_patch(s, p) +} + +/// `surface.boundary` with `curves`: one patch through two to four sketch curves in loop order. +fn loop_patch(s: &mut Session, p: &Value) -> Result { let cmd = "surface.boundary"; + if p.get("v").is_some_and(|v| !v.is_null()) || p.get("operation").is_some_and(|v| !v.is_null()) { + return Err(bad(cmd, "`v` and `operation` belong to a curve network: give its curves as `u` and `v`")); + } 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")); + return Err(bad(cmd, "`curves` must list 2 to 4 sketch curves in loop order (or give a network as `u` and `v`)")); } // 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. @@ -98,9 +116,134 @@ fn boundary(s: &mut Session, p: &Value) -> Result { .ok_or_else(|| bad(cmd, format!("`samples` must be a whole number of stations between {MIN_SAMPLES} and {MAX_SAMPLES}")))? } }; + if crown != 0.0 && curves.len() >= 3 { + return Err(bad(cmd, "not supported yet: `crown` on a patch through three or four curves (it is the flat patch of its sides)")); + } add_feature(s, p, FeatureKind::SurfaceBoundary { sketch: SketchId(sketch), curves, orders, crown, samples }) } +/// `surface.boundary` with `u` (and `v`): one surface through a curve network. +fn network(s: &mut Session, p: &Value) -> Result { + let cmd = "surface.boundary"; + for key in ["curves", "continuity", "crown", "samples"] { + if p.get(key).is_some_and(|v| !v.is_null()) { + return Err(bad(cmd, format!("`{key}` belongs to a loop of sketch curves: a network gives each curve its own continuity"))); + } + } + let operation = super::features::operation(p, cmd)?; + if !matches!(operation, Operation::NewBody | Operation::Join) { + return Err(bad(cmd, "`operation` must be new or join")); + } + // The sketch a bare curve id belongs to, found only when one is used. + let mut default_sketch: Option = None; + let u = network_curves(s, p, "u", &mut default_sketch)?; + let v = network_curves(s, p, "v", &mut default_sketch)?; + if u.is_empty() { + return Err(bad(cmd, "`u` must list the curves of one direction")); + } + if u.len() + v.len() < 2 { + return Err(bad(cmd, "a boundary surface needs two curves in one direction, or crossing curves")); + } + if operation == Operation::Join && !u.iter().chain(&v).any(|c| c.body.is_some()) { + return Err(bad(cmd, "`operation` join needs a curve on a body's edges: the surface is sewn into that body")); + } + add_feature(s, p, FeatureKind::SurfaceNetwork { u, v, operation }) +} + +/// One direction of a network (`key` is `u` or `v`), checked: each curve sketch curves that exist +/// (planar or 3D), a body's edges with points on them, or a point; a continuity only on edges. +fn network_curves(s: &Session, p: &Value, key: &str, default_sketch: &mut Option) -> Result> { + let cmd = "surface.boundary"; + let list = match p.get(key) { + None | Some(Value::Null) => return Ok(Vec::new()), + Some(Value::Array(a)) => a, + Some(_) => return Err(bad(cmd, format!("`{key}` must list curves"))), + }; + if list.len() > MAX_NETWORK_CURVES { + return Err(bad(cmd, format!("`{key}` may list at most {MAX_NETWORK_CURVES} curves"))); + } + let mut out = Vec::with_capacity(list.len()); + for (k, item) in list.iter().enumerate() { + let name = format!("{key}[{k}]"); + let mut curve = NetworkCurve { sketch: None, curves: Vec::new(), body: None, edges: Vec::new(), point: None, continuity: Continuity::G0 }; + let obj = match item { + Value::String(id) => { + curve.curves = vec![id.clone()]; + None + } + Value::Object(_) => Some(item), + _ => return Err(bad(cmd, format!("{name} must be a curve id or {{sketch?, curves}}, {{body, edges}} or {{point}}"))), + }; + if let Some(o) = obj { + let has = |k: &str| o.get(k).is_some_and(|v| !v.is_null()); + let kinds = [has("curves"), has("body") || has("edges"), has("point")]; + if kinds.iter().filter(|x| **x).count() != 1 { + return Err(bad(cmd, format!("{name} must be one of {{sketch?, curves}}, {{body, edges}} or {{point}}"))); + } + if has("point") { + curve.point = Some(Position(o.get("point").and_then(vec3).ok_or_else(|| bad(cmd, format!("{name}.point must be [x, y, z]")))?)); + } else if kinds[1] { + let body = str_(o, "body").ok_or_else(|| bad(cmd, format!("{name}.body must name a body")))?; + if s.model.state().body(body).is_none() { + return Err(bad(cmd, format!("no body `{body}`"))); + } + let edges = + o.get("edges").and_then(Value::as_array).ok_or_else(|| bad(cmd, format!("{name}.edges must list points on the edges")))?; + if edges.is_empty() || edges.len() > MAX_PICKS { + return Err(bad(cmd, format!("{name}.edges must list 1 to {MAX_PICKS} points on a chain of edges"))); + } + let pts: Vec = edges + .iter() + .map(|e| vec3(e).map(Position).ok_or_else(|| bad(cmd, format!("{name}.edges: a point must be [x, y, z]")))) + .collect::>()?; + curve.body = Some(BodyName::from(body.to_string())); + curve.edges = pts; + } else { + curve.curves = string_list(o, "curves"); + if o.get("sketch").is_some_and(|v| !v.is_null()) { + curve.sketch = Some(SketchId(sketch_id(s, o.get("sketch"), cmd, &format!("{name}.sketch"))?)); + } + } + if let Some(g) = o.get("continuity").filter(|v| !v.is_null()) { + let g = g.as_str().and_then(Continuity::parse).ok_or_else(|| bad(cmd, format!("{name}.continuity must be g0, g1 or g2")))?; + if g == Continuity::G3 { + return Err(bad(cmd, format!("not supported yet: G3 on {name} (a curve network holds up to g2)"))); + } + if g != Continuity::G0 && curve.body.is_none() { + return Err(bad(cmd, format!("{name}: g1 and g2 apply to a curve on a body's edges (a surface's open edges)"))); + } + curve.continuity = g; + } + } + if curve.point.is_none() && curve.body.is_none() { + if curve.curves.is_empty() || curve.curves.len() > MAX_CHAIN { + return Err(bad(cmd, format!("{name} must list 1 to {MAX_CHAIN} sketch curves"))); + } + let sketch = match curve.sketch { + Some(id) => id.0, + None => match *default_sketch { + Some(id) => id, + None => { + let id = feature_sketch(s, p, cmd)?; + *default_sketch = Some(id); + id + } + }, + }; + let state = s.model.state(); + let sk = state.sketch(sketch).ok_or_else(|| bad(cmd, "the sketch has not been solved"))?; + for id in &curve.curves { + if sk.sketch.curve_index(id).is_none() && sk.sketch.wire_index(id).is_none() { + return Err(bad(cmd, format!("{name}: no sketch curve `{id}`"))); + } + } + curve.sketch = Some(SketchId(sketch)); + } + out.push(curve); + } + Ok(out) +} + #[cfg(test)] mod tests { use serde_json::{Value, json}; @@ -154,27 +297,19 @@ mod tests { } #[test] - fn a_boundary_patch_refuses_a_loop_and_values_it_cannot_hold() { + fn a_boundary_patch_refuses_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; + // A loop of three: the crown is the bridge's, not the flat loop patch's. 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}"); + let e = s.execute("surface.boundary", &json!({"curves": [a, b, c], "crown": 1})).unwrap_err().to_string(); + assert!(e.contains("not supported yet") && e.contains("crown"), "three boundaries with a crown: {e}"); + // Curves out of loop order. + let e = s.execute("surface.boundary", &json!({"curves": [curves[0], curves[1], extra]})).unwrap_err().to_string(); + assert!(e.contains("loop order"), "{e}"); // Parameters and curves that make no sense. let cases: Vec<(Value, &str)> = vec![ (json!({"curves": [curves[0]]}), "2 to 4"), @@ -186,6 +321,7 @@ mod tests { (json!({"curves": curves, "crown": "full"}), "crown"), (json!({"curves": curves, "samples": 3}), "samples"), (json!({"curves": curves, "samples": 2.5}), "samples"), + (json!({"curves": curves, "operation": "join"}), "network"), ]; for (params, want) in cases { let e = s.execute("surface.boundary", ¶ms).unwrap_err().to_string(); @@ -195,3 +331,7 @@ mod tests { assert!(s.model.state().bodies.is_empty(), "a refused patch left a body behind"); } } + +#[cfg(test)] +#[path = "surface_boundary_tests.rs"] +mod network_tests; diff --git a/crates/engine/src/cmd/surface_boundary_tests.rs b/crates/engine/src/cmd/surface_boundary_tests.rs new file mode 100644 index 00000000..fd9f0f81 --- /dev/null +++ b/crates/engine/src/cmd/surface_boundary_tests.rs @@ -0,0 +1,267 @@ +//! `surface.boundary` through loops of three or four sketch curves and through curve networks +//! (#763): analytic areas and shapes, G1 and G2 measured across the edge a surface shares with a +//! body (`inspect.continuity`), and refusals of hostile networks that never panic. + +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 first(v: &Value) -> String { + v.as_array().and_then(|a| a.first()).and_then(Value::as_str).map(str::to_string).unwrap_or_default() +} + +/// A sketch line's id. +fn line(s: &mut Session, a: [f64; 2], b: [f64; 2]) -> String { + first(&run(s, "sketch.line", json!({"points": [a, b]}))["created"]["curves"]) +} + +/// A 3D sketch line's id (the active sketch must have 3D Sketch on). +fn line3d(s: &mut Session, a: [f64; 3], b: [f64; 3]) -> String { + first(&run(s, "sketch.line3d", json!({"points": [a, b]}))["curves"]) +} + +/// The area of the one face of `body`, and the area `inspect.measure` gives for the whole model. +fn areas(s: &mut Session, body: &str) -> (f64, f64) { + let f = run(s, "model.faces", json!({"body": body})); + let faces = f["faces"].as_array().cloned().unwrap_or_default(); + assert_eq!(faces.len(), 1, "one face: {f}"); + let face = faces[0]["area"].as_f64().unwrap_or(f64::NAN); + let m = run(s, "inspect.measure", json!({})); + (face, m["total"]["area_mm2"].as_f64().unwrap_or(f64::NAN)) +} + +/// Three and four sketch curves in loop order: the flat patch of the loop. A triangle and a square +/// of lines, and a square whose top is a half circle bulging out, each measure their own area. +#[test] +fn three_and_four_curve_loops_are_flat_patches() { + use std::f64::consts::PI; + let cases: Vec<(Vec<[[f64; 2]; 2]>, f64, f64)> = vec![ + (vec![[[0.0, 0.0], [20.0, 0.0]], [[20.0, 0.0], [10.0, 15.0]], [[10.0, 15.0], [0.0, 0.0]]], 150.0, 1e-6), + (vec![[[0.0, 0.0], [10.0, 0.0]], [[10.0, 0.0], [10.0, 10.0]], [[10.0, 10.0], [0.0, 10.0]], [[0.0, 10.0], [0.0, 0.0]]], 100.0, 1e-6), + ]; + for (lines, want, tol) in cases { + for continuity in [json!("g0"), json!("g1"), json!("g3")] { + let mut s = Session::default(); + run(&mut s, "sketch.create", json!({"plane": "XY"})); + let curves: Vec = lines.iter().map(|[a, b]| line(&mut s, *a, *b)).collect(); + let r = run(&mut s, "surface.boundary", json!({"curves": curves, "continuity": continuity})); + let body = first(&r["created_bodies"]); + let (face, total) = areas(&mut s, &body); + assert!( + (face - want).abs() < tol * want && (total - want).abs() < tol * want, + "{} sides, {continuity}: {face} / {total} vs {want}", + lines.len() + ); + } + } + // A square with a half circle for its top, bulging out: 100 + 12.5π. + let mut s = Session::default(); + run(&mut s, "sketch.create", json!({"plane": "XY"})); + let bottom = line(&mut s, [0.0, 0.0], [10.0, 0.0]); + let right = line(&mut s, [10.0, 0.0], [10.0, 10.0]); + let arc = first(&run(&mut s, "sketch.arc.three_point", json!({"start": [10, 10], "end": [0, 10], "through": [5, 15]}))["created"]["curves"]); + let left = line(&mut s, [0.0, 10.0], [0.0, 0.0]); + let r = run(&mut s, "surface.boundary", json!({"curves": [bottom, right, arc, left]})); + let (face, total) = areas(&mut s, &first(&r["created_bodies"])); + let want = 100.0 + 12.5 * PI; + assert!((face - want).abs() < 2e-3 * want && (total - want).abs() < 2e-3 * want, "{face} / {total} vs {want}"); +} + +/// Simpson's rule over [0, a] × [0, a] with n (even) intervals each way. +fn simpson2(f: impl Fn(f64, f64) -> f64, a: f64, n: usize) -> f64 { + let h = a / n as f64; + let w = |k: usize| { + if k == 0 || k == n { + 1.0 + } else if k % 2 == 1 { + 4.0 + } else { + 2.0 + } + }; + let mut sum = 0.0; + for i in 0..=n { + for j in 0..=n { + sum += w(i) * w(j) * f(i as f64 * h, j as f64 * h); + } + } + sum * h * h / 9.0 +} + +/// Four straight 3D sketch lines bounding the hyperbolic paraboloid z = x·y/10 over the 10 × 10 +/// square: the network's Coons patch is that surface (every vertex of its mesh on it), and its +/// area is the paraboloid's, ∫∫ √(1 + (x/10)² + (y/10)²). +#[test] +fn a_network_of_3d_lines_is_a_hyperbolic_paraboloid() { + let mut s = Session::default(); + run(&mut s, "sketch.create", json!({"plane": "XY"})); + run(&mut s, "sketch.options", json!({"sketch_3d": true})); + let a = line3d(&mut s, [0.0, 0.0, 0.0], [10.0, 0.0, 0.0]); + let c = line3d(&mut s, [0.0, 10.0, 0.0], [10.0, 10.0, 10.0]); + let d = line3d(&mut s, [0.0, 0.0, 0.0], [0.0, 10.0, 0.0]); + let b = line3d(&mut s, [10.0, 0.0, 0.0], [10.0, 10.0, 10.0]); + let r = run(&mut s, "surface.boundary", json!({"u": [a, c], "v": [d, b]})); + let body = first(&r["created_bodies"]); + let want = simpson2(|x, y| (1.0 + (x / 10.0).powi(2) + (y / 10.0).powi(2)).sqrt(), 10.0, 200); + let (face, total) = areas(&mut s, &body); + // The face's own area comes from the coarse reference mesh (a chord of a fiftieth of a per + // cent of the body); the measuring mesh is refined onto the surface. + assert!((face - want).abs() < 1e-2 * want && (total - want).abs() < 1e-3 * want, "{face} / {total} vs {want}"); + let st = s.world_state(); + let mesh = st.body(&body).map(|b| b.body.tessellate(0.01)).and_then(Result::ok).unwrap_or_default(); + assert!(!mesh.positions.is_empty()); + for p in &mesh.positions { + assert!((p.z - p.x * p.y / 10.0).abs() < 1e-6, "{p:?} is off the paraboloid"); + } +} + +/// Curves of two planar sketches, one picked by a bare id in the sketch being edited: the loft +/// between a line on the ground and one 10 mm across and 5 mm up is the flat strip between them. +#[test] +fn a_network_takes_curves_from_several_sketches() { + let mut s = Session::default(); + run(&mut s, "sketch.create", json!({"plane": "XY", "name": "Low"})); + let low = line(&mut s, [0.0, 0.0], [10.0, 0.0]); + run(&mut s, "sketch.finish", json!({})); + run(&mut s, "sketch.create", json!({"plane": "XY", "offset": 5, "name": "High"})); + let high = line(&mut s, [0.0, 10.0], [10.0, 10.0]); + let r = run(&mut s, "surface.boundary", json!({"u": [{"sketch": "Low", "curves": [low]}, high]})); + let (face, _) = areas(&mut s, &first(&r["created_bodies"])); + let want = 10.0 * 125f64.sqrt(); + assert!((face - want).abs() < 1e-6 * want, "{face} vs {want}"); +} + +/// A quarter cylinder sheet (R = 10, axis along x, its straight top edge at y = 10, z = 0 open) +/// and a line 3D-sketched above it; returns the sheet's name and the line's id. +fn cylinder_and_line(s: &mut Session) -> (String, String) { + run(s, "sketch.create", json!({"plane": "YZ", "name": "Arc"})); + let h = 10.0 * std::f64::consts::FRAC_1_SQRT_2; + let arc = first(&run(s, "sketch.arc.three_point", json!({"start": [0, -10], "end": [10, 0], "through": [h, -h]}))["created"]["curves"]); + run(s, "sketch.finish", json!({})); + let sheet = first(&run(s, "surface.extrude", json!({"sketch": "Arc", "curves": [arc], "distance": 10}))["created_bodies"]); + run(s, "sketch.create", json!({"plane": "XY", "name": "Top"})); + run(s, "sketch.options", json!({"sketch_3d": true})); + let top = line3d(s, [0.0, 3.0, 12.0], [10.0, 3.0, 12.0]); + run(s, "sketch.finish", json!({})); + (sheet, top) +} + +/// The edge a surface built from the cylinder's top edge shares with it, measured with +/// `inspect.continuity`: (normal angle in degrees, curvature jump per mm). +fn seam(s: &mut Session, body: &str) -> (f64, f64) { + let v = run(s, "inspect.continuity", json!({"body": body, "edges": [[5.0, 10.0, 0.0]]})); + let e = v["bodies"][0]["edges"].as_array().and_then(|a| a.first()).cloned().unwrap_or_default(); + (e["normal_angle_deg"]["value"].as_f64().unwrap_or(f64::NAN), e["curvature_jump_per_mm"]["value"].as_f64().unwrap_or(f64::NAN)) +} + +/// G1 and G2 against a body: a network from the cylinder's open top edge to a line above, joined +/// into the cylinder. Across the shared edge the normals agree within 0.1° for G1 and G2; the +/// curvature jumps by much of the cylinder's 1/R with G1 and by under 1% of it with G2. +#[test] +fn g1_and_g2_networks_meet_a_body_across_its_edge() { + let mut jumps = Vec::new(); + for g in ["g1", "g2"] { + let mut s = Session::default(); + let (sheet, top) = cylinder_and_line(&mut s); + let r = run( + &mut s, + "surface.boundary", + json!({"u": [{"body": sheet, "edges": [[5.0, 10.0, 0.0]], "continuity": g}, {"sketch": "Top", "curves": [top]}], "operation": "join"}), + ); + assert!(r["created_bodies"].as_array().is_none_or(Vec::is_empty), "a join makes no new body: {r}"); + let f = run(&mut s, "model.faces", json!({"body": sheet})); + assert_eq!(f["faces"].as_array().map(Vec::len), Some(2), "the surface is sewn into the sheet: {f}"); + let (angle, jump) = seam(&mut s, &sheet); + assert!(angle < 0.1, "{g}: normals {angle}° apart across the edge"); + jumps.push(jump); + } + let r = 10.0; + assert!(jumps[0] > 0.2 / r, "G1 leaves the curvature behind: jump {}", jumps[0]); + assert!(jumps[1] < 0.01 / r, "G2 carries the cylinder's curvature across: jump {}", jumps[1]); +} + +/// Networks that cannot be built are refused with a reason — by the command's own checks or by the +/// kernel — and leave nothing behind; none of them panics. +#[test] +fn hostile_networks_are_refused() { + let mut s = Session::default(); + let (sheet, top) = cylinder_and_line(&mut s); + run(&mut s, "solid.box", json!({"length": 10, "width": 10, "height": 10, "corner": [40, 40, 0]})); + let boxed = s.world_state().bodies.last().map(|b| b.name.clone()).unwrap_or_default(); + run(&mut s, "sketch.create", json!({"plane": "XY", "name": "Flat"})); + let a = line(&mut s, [0.0, 0.0], [10.0, 0.0]); + let b = line(&mut s, [0.0, 5.0], [10.0, 5.0]); + let far = line(&mut s, [0.0, 20.0], [10.0, 20.0]); + let ring = first(&run(&mut s, "sketch.circle.center", json!({"center": [50, 50], "radius": 3}))["created"]["curves"]); + let bodies = s.world_state().bodies.len(); + let edge = |g: &str| json!({"body": sheet, "edges": [[5.0, 10.0, 0.0]], "continuity": g}); + let many: Vec = (0..33).map(|_| json!(a)).collect(); + let cases: Vec<(Value, &str)> = vec![ + (json!({"u": "l1"}), "list"), + (json!({"u": [42, a]}), "curve id"), + (json!({"u": [{}, a]}), "one of"), + (json!({"u": [{"curves": [], "sketch": "Flat"}, a]}), "sketch curves"), + (json!({"u": [{"curves": ["nope"]}, a]}), "no sketch curve"), + (json!({"u": [{"curves": [a], "sketch": "Nope"}, b]}), "no sketch"), + (json!({"u": [{"body": "Nope", "edges": [[0, 0, 0]]}, a]}), "no body"), + (json!({"u": [{"body": sheet, "edges": []}, a]}), "edges"), + (json!({"u": [{"body": sheet, "edges": [["x", 0, 0]]}, a]}), "[x, y, z]"), + (json!({"u": [{"body": sheet, "edges": [[1e300, 0, 0]]}, a]}), "[x, y, z]"), + (json!({"u": [{"point": [0, 0]}, a]}), "point"), + (json!({"u": [{"curves": [a], "body": sheet, "edges": [[5, 10, 0]]}, b]}), "one of"), + (json!({"u": [{"curves": [a], "continuity": "g1"}, b]}), "body's edges"), + (json!({"u": [edge("g3"), {"sketch": "Top", "curves": [top]}]}), "not supported yet"), + (json!({"u": [edge("g9"), {"sketch": "Top", "curves": [top]}]}), "continuity"), + (json!({"u": many}), "at most"), + (json!({"u": [a, b], "curves": [a, b]}), "loop"), + (json!({"u": [a, b], "crown": 1}), "loop"), + (json!({"u": [a, b], "operation": "cut"}), "operation"), + (json!({"u": [a, b], "operation": "join"}), "join"), + (json!({"u": [a]}), "two curves"), + (json!({"u": [], "v": [a, b]}), "`u`"), + // The kernel's refusals: curves that miss, a closed curve, G1 on an edge between two faces. + (json!({"u": [a, b], "v": [{"curves": [far]}, {"curves": [a]}]}), "don't meet"), + (json!({"u": [{"curves": [ring]}, b]}), "closed"), + (json!({"u": [{"body": boxed, "edges": [[45, 40, 0]], "continuity": "g1"}, b]}), "open edge"), + (json!({"u": [a, {"point": [5, 2, 0]}, b]}), "point"), + ]; + for (params, want) in cases { + let e = match s.execute("surface.boundary", ¶ms) { + Ok(v) => panic!("{params} was built: {v}"), + Err(e) => e, + }; + assert!(!matches!(e, crate::EngineError::Internal(..)), "{params}: {e}"); + assert!(e.to_string().contains(want), "{params}: {e}"); + } + assert_eq!(s.world_state().bodies.len(), bodies, "a refused network left a body behind"); +} + +/// An edit (`timeline.redefine`) may turn a loop of sketch curves into the network of the same +/// curves and back: the feature keeps its place and its name, and the surface its area. +#[test] +fn an_edit_turns_a_loop_into_a_network_and_back() { + let mut s = Session::default(); + run(&mut s, "sketch.create", json!({"plane": "XY", "name": "Square"})); + let sides: Vec = [[[0.0, 0.0], [10.0, 0.0]], [[10.0, 0.0], [10.0, 10.0]], [[10.0, 10.0], [0.0, 10.0]], [[0.0, 10.0], [0.0, 0.0]]] + .iter() + .map(|[a, b]| line(&mut s, *a, *b)) + .collect(); + run(&mut s, "sketch.finish", json!({})); + run(&mut s, "surface.boundary", json!({"sketch": "Square", "curves": sides, "name": "Patch"})); + let network = json!({"sketch": "Square", "u": [sides[0], sides[2]], "v": [sides[1], sides[3]]}); + run(&mut s, "timeline.redefine", json!({"feature": "Patch", "command": "surface.boundary", "params": network})); + let kind = |s: &Session| s.doc.find_feature("Patch").map(|f| f.kind.clone()); + assert!(matches!(kind(&s), Some(solidcraft_doc::FeatureKind::SurfaceNetwork { .. })), "{:?}", kind(&s)); + let body = s.world_state().bodies.first().map(|b| b.name.clone()).unwrap_or_default(); + assert!((areas(&mut s, &body).0 - 100.0).abs() < 1e-6); + run(&mut s, "timeline.redefine", json!({"feature": "Patch", "command": "surface.boundary", "params": {"sketch": "Square", "curves": sides}})); + assert!(matches!(kind(&s), Some(solidcraft_doc::FeatureKind::SurfaceBoundary { .. })), "{:?}", kind(&s)); + assert!((areas(&mut s, &body).0 - 100.0).abs() < 1e-6); +} diff --git a/crates/engine/src/tests.rs b/crates/engine/src/tests.rs index 41505eb2..fe739c29 100644 --- a/crates/engine/src/tests.rs +++ b/crates/engine/src/tests.rs @@ -807,6 +807,8 @@ fn hostile_params_never_panic() { "continuity", "crown", "faces", + "u", + "v", ]; // Commands are independent; worker threads take them one at a time (their costs vary a lot). let ids: Vec<&str> = command_specs() diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index ab75d895..9b968092 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -1,7 +1,8 @@ //! Geometry foundation for SolidCraft: 2D/3D vectors, sketch planes, axis-aligned boxes, //! triangle meshes with their measures (volume, area, centroid), and freeform surface maths -//! (B-spline basis and NURBS surfaces in [`bspline`], one-surface fills of closed loops in -//! [`fill`], T-splines in [`tspline`]). Units are millimetres. +//! (B-spline basis and NURBS surfaces in [`bspline`], spline spaces and the operators that make +//! splines compatible in [`spline_space`], one-surface fills of closed loops in [`fill`], +//! surfaces through curve networks in [`network`], T-splines in [`tspline`]). Units are millimetres. //! //! This crate has no knowledge of sketches, documents or a modelling backend; every other crate //! builds on it. @@ -19,8 +20,10 @@ pub mod fill; pub mod fill_general; mod mesh; pub mod mesh_distance; +pub mod network; mod plane; mod profile; +pub mod spline_space; pub mod tess; pub mod tspline; mod vec; diff --git a/crates/geom/src/network.rs b/crates/geom/src/network.rs new file mode 100644 index 00000000..06ecd688 --- /dev/null +++ b/crates/geom/src/network.rs @@ -0,0 +1,1623 @@ +//! Surfaces through curve networks (Surface › Boundary): one non-rational B-spline surface +//! through curves in one direction (a loft that skins them), or through curves in two +//! directions that cross each other (a Gordon surface; with two curves each way, the bilinearly +//! blended Coons patch). A boundary curve may carry a continuity condition with a neighbouring +//! face ([`CrossFace`]): G1, the surface leaves that face tangent to it, or G2, with its normal +//! curvature as well. +//! +//! Method: +//! - Each curve is measured by (normalised) arc length. Where curves of the two directions cross +//! ([`closest`]: refused when they miss by more than the crossing tolerance, or cross twice), +//! the crossing positions give the network's parameters: `ū_j` (where v-curve `j` crosses the +//! u-curves) and `v̄_i` the other way, averages of the relative positions over the curves. +//! With curves in one direction only, `v̄_i` is the mean distance between neighbours (chord +//! length, NURBS Book §9.2.1 applied across). +//! - Each curve is reparameterised onto the network parameter by a monotone cubic through its +//! crossings, so that every u-curve passes crossing `j` at `ū_j`; the small gap at a crossing +//! (within tolerance) is closed by moving each curve onto the crossing's mid-point, the move +//! spread linearly between crossings. +//! - All curves of one direction, and the cross-boundary derivative functions, are fitted in one +//! common cubic spline space whose breaks include every crossing parameter (refined until each +//! is within tolerance): so they are compatible from the start and pass exactly through the +//! crossing points. +//! - The surface is Gordon's `S = P_u + P_v − P_uv` (W. J. Gordon, "Spline-blended surface +//! interpolation through curve networks", 1969; the NURBS Book, ch. 10, interpolation of a +//! bidirectional curve network): `P_u` interpolates the u-curves across v (skinning: the same +//! interpolation applied to every column of control points), `P_v` the v-curves across u, and +//! `P_uv` the crossing points both ways, by the same two interpolation operators. `S` then +//! reproduces every curve. Two curves each way and linear operators give the bilinearly blended +//! Coons patch (S. A. Coons, 1967); curves in one direction only give `S = P_u`. +//! - A continuity condition enters the across-v operator as Hermite end conditions (NURBS Book +//! §9.2.2): the surface's first (G1) and second (G2) v-derivatives along the boundary curve. +//! The first lies in the face's tangent plane, pointing away from the face (or into it at the +//! far end); the second has the normal component `II(T, T)` that matches the face's normal +//! curvature in that direction (the face's second fundamental form). A curve that crosses +//! such a boundary is part of the surface, so it must already leave the face tangent to it +//! (G1) and with its curvature (G2), within tolerance; otherwise the network is refused, +//! naming the curve. +//! - The three surfaces are brought into one space ([`crate::spline_space::embed`]: Bézier degree +//! elevation and knot insertion) and summed. +//! +//! - A side may collapse to a point ([`PointCurve`]: the apex of a three-sided patch, where the +//! two sides beside it meet). The curves through it cross it at the point itself, so the side +//! stays one point; the surface is singular there, as a three-sided tensor patch must be. +//! +//! Continuity conditions may sit on curves of one direction only (on two boundaries that meet +//! at a corner the cross derivatives would also have to agree on the twist there: not supported +//! yet), and only on the first and last curves of that direction. Closed curves are not +//! supported yet. +//! +//! Input is hostile: sizes are capped, every evaluation is checked, and failures are +//! [`SplineError`]s naming the curves; nothing here panics. + +use crate::Vec3; +use crate::bspline::{NurbsSurface, SplineError, SplineResult, invalid}; +use crate::spline_space::{Cond, Interpolant, Linear, Monotone, Space, embed}; + +/// Most curves in one direction of a network. +pub const MAX_NETWORK_CURVES: usize = 32; +/// Most spans of the common fitting space in one direction. +const MAX_SPANS: usize = 512; +/// Fewest spans the fitting starts with in one direction. +const MIN_SPANS: usize = 8; +/// Samples of a curve's arc-length table. +const TABLE: usize = 256; +/// How far a curve that crosses a G1 or G2 boundary may lean out of the face's tangent plane +/// there (radians, about 1°). +const TANGENT_TOL: f64 = 0.0175; +/// How much a curve that crosses a G2 boundary may differ from the face's normal curvature in +/// its direction there, relative to that curvature (or to one over the network's size, when +/// larger). +const CURVATURE_TOL: f64 = 0.02; + +/// A curve of a network: its points along a parameter running over `[0, 1]` (any regular +/// parameterisation; the network measures its arc length). +pub trait NetCurve { + fn point(&mut self, t: f64) -> Option; + + /// A side collapsed to a point (a three-sided patch: the apex where two sides meet). Any + /// other curve of no length is refused. + fn is_point(&self) -> bool { + false + } +} + +/// A side of a network collapsed to one point: the apex of a three-sided patch, where the two +/// sides beside it meet. +pub struct PointCurve(pub Vec3); + +impl NetCurve for PointCurve { + fn point(&mut self, _: f64) -> Option { + Some(self.0) + } + + fn is_point(&self) -> bool { + true + } +} + +/// A face along a boundary curve, seen from a point of the curve. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct FaceFrame { + /// The face's unit normal there. + pub normal: Vec3, + /// A unit vector in the face's tangent plane, square to the curve, pointing away from the + /// face (out across its edge). + pub outward: Vec3, + /// The face's second fundamental form against `normal` (`II(w, w)`: the normal component of + /// the second derivative of a curve on the face with first derivative `w`), in the + /// orthonormal basis of the tangent asked about `t` (made unit in the tangent plane, running + /// the way it was given) and `outward` `b`: `[II(t, t), II(t, b), II(b, b)]`. + pub form: [f64; 3], +} + +/// The face a boundary curve runs along: its frame at a point of the curve (`tangent`: the +/// curve's direction there, either way). +pub trait CrossFace { + fn frame(&mut self, at: Vec3, tangent: Vec3) -> Option; +} + +/// How smoothly a surface meets a neighbouring face across a boundary curve. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub enum Continuity { + /// Position only. + #[default] + G0, + /// Tangent: the same tangent plane along the curve. + G1, + /// Curvature: the same tangent plane and normal curvature along the curve. + G2, +} + +impl Continuity { + /// The number of cross-boundary derivatives it fixes. + pub fn order(self) -> usize { + match self { + Continuity::G0 => 0, + Continuity::G1 => 1, + Continuity::G2 => 2, + } + } +} + +/// One curve of a network: a name for messages, the curve, and its continuity with a face. +pub struct NetworkCurve<'a> { + pub label: String, + pub curve: Box, + pub condition: Option<(Continuity, Box)>, +} + +/// The tolerances of a network (mm): how close curves must come to cross, and how far the fitted +/// curves may stray from the given ones. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Tolerances { + pub cross: f64, + pub fit: f64, +} + +/// A curve measured by arc length: `point(s)` with `s` in `[0, 1]` (reversed when `flipped`). +struct Track<'a> { + label: String, + curve: Box, + /// Normalised arc length and curve parameter at each table sample. + s: Vec, + t: Vec, + pts: Vec, + length: f64, + flipped: bool, + /// A side collapsed to a point ([`PointCurve`]). + point: bool, + /// The finite-difference step of [`Track::ders`]: a quarter of the smallest table cell. + step: f64, + condition: Option<(Continuity, Box)>, +} + +impl<'a> Track<'a> { + fn new(c: NetworkCurve<'a>) -> SplineResult> { + let NetworkCurve { label, mut curve, condition } = c; + let mut pts = Vec::with_capacity(TABLE + 1); + let mut t = Vec::with_capacity(TABLE + 1); + for k in 0..=TABLE { + let tk = k as f64 / TABLE as f64; + let p = curve.point(tk).filter(|p| p.is_finite()).ok_or_else(|| SplineError::Invalid(format!("{label} can't be evaluated")))?; + pts.push(p); + t.push(tk); + } + let mut s = vec![0.0]; + let mut acc = 0.0; + for w in pts.windows(2) { + acc += w[0].dist(w[1]); + s.push(acc); + } + let point = curve.is_point(); + if point { + // Its parameter runs evenly: every value is the same point. + if acc > 1e-9 || !acc.is_finite() { + return invalid(format!("{label} is meant to be a point but moves")); + } + if condition.as_ref().is_some_and(|(g, _)| *g != Continuity::G0) { + return invalid(format!("{label} is a point: a continuity condition needs a curve")); + } + let step = step(&t); + return Ok(Track { label, curve, s: t.clone(), t, pts, length: 0.0, flipped: false, point, step, condition: None }); + } + if !(acc > 1e-9 && acc.is_finite()) { + return invalid(format!("{label} has no length")); + } + for x in &mut s { + *x /= acc; + } + let condition = condition.filter(|(g, _)| *g != Continuity::G0); + let step = step(&s); + Ok(Track { label, curve, s, t, pts, length: acc, flipped: false, point, step, condition }) + } + + fn continuity(&self) -> Continuity { + self.condition.as_ref().map_or(Continuity::G0, |(g, _)| *g) + } + + /// The point at normalised arc length `s` (clamped to `[0, 1]`). + fn point(&mut self, s: f64) -> SplineResult { + let s = if self.flipped { 1.0 - s } else { s }.clamp(0.0, 1.0); + let k = self.s.partition_point(|x| *x <= s).saturating_sub(1).min(TABLE - 1); + let (s0, s1, t0, t1) = (self.s[k], self.s[k + 1], self.t[k], self.t[k + 1]); + let t = if s1 > s0 { t0 + (t1 - t0) * (s - s0) / (s1 - s0) } else { t0 }; + self.curve.point(t).filter(|p| p.is_finite()).ok_or_else(|| SplineError::Invalid(format!("{} can't be evaluated", self.label))) + } + + /// The table's points in the curve's current direction. + fn polyline(&self) -> Vec { + if self.flipped { self.pts.iter().rev().copied().collect() } else { self.pts.clone() } + } + + /// The first two derivatives by `s` at `s` (one-sided at the ends), by finite differences + /// within one table cell. + fn ders(&mut self, s: f64) -> SplineResult<(Vec3, Vec3)> { + let h = self.step; + if s - 2.0 * h < 0.0 { + let f: Vec = (0..4).map(|k| self.point(s + k as f64 * h)).collect::>()?; + Ok(((f[0] * -3.0 + f[1] * 4.0 - f[2]) / (2.0 * h), (f[0] * 2.0 - f[1] * 5.0 + f[2] * 4.0 - f[3]) / (h * h))) + } else if s + 2.0 * h > 1.0 { + let f: Vec = (0..4).map(|k| self.point(s - k as f64 * h)).collect::>()?; + Ok(((f[0] * 3.0 - f[1] * 4.0 + f[2]) / (2.0 * h), (f[0] * 2.0 - f[1] * 5.0 + f[2] * 4.0 - f[3]) / (h * h))) + } else { + let (a, b, c) = (self.point(s - h)?, self.point(s)?, self.point(s + h)?); + Ok(((c - a) / (2.0 * h), (a - b * 2.0 + c) / (h * h))) + } + } +} + +/// The finite-difference step along a table of parameters: a quarter of its smallest cell, at +/// most 1e-3. +fn step(s: &[f64]) -> f64 { + let cell = s.windows(2).map(|w| w[1] - w[0]).filter(|d| *d > 0.0).fold(f64::INFINITY, f64::min); + (cell / 4.0).min(1e-3) +} + +/// Where two curves cross: the arc-length positions on each and the mid-point between them. +#[derive(Clone, Copy, Debug)] +struct Crossing { + s: f64, + r: f64, + at: Vec3, +} + +/// The one place where `a` and `b` come within `tol` of each other. Refused when they don't, or +/// do at two places. +fn closest(a: &mut Track, b: &mut Track, tol: f64) -> SplineResult { + let (pa, pb) = (a.polyline(), b.polyline()); + let longest = |p: &[Vec3]| p.windows(2).map(|w| w[0].dist(w[1])).fold(0.0, f64::max); + let near = tol + longest(&pa) + longest(&pb); + // Distance from each sample of `a` to the polyline of `b`, and the nearest segment. + let mut d = Vec::with_capacity(pa.len()); + for p in &pa { + let mut best = (f64::INFINITY, 0usize); + for (k, w) in pb.windows(2).enumerate() { + let x = p.dist_to_segment(w[0], w[1]); + if x < best.0 { + best = (x, k); + } + } + d.push(best); + } + // Runs of samples near `b`: one candidate crossing each, refined from its nearest sample. + let n = d.len(); + let mut starts = Vec::new(); + let mut k = 0; + while k < n { + if d[k].0 > near { + k += 1; + continue; + } + let mut m = k; + while k < n && d[k].0 <= near { + if d[k].0 < d[m].0 { + m = k; + } + k += 1; + } + starts.push(m); + } + if starts.is_empty() { + starts.push((0..n).min_by(|x, y| d[*x].0.total_cmp(&d[*y].0)).unwrap_or(0)); + } + let mut found: Vec = Vec::new(); + let mut miss: Option<(f64, Vec3)> = None; + for m in starts { + let c = refine(a, b, m as f64 / TABLE as f64, (d[m].1 as f64 + 0.5) / TABLE as f64)?; + let gap = a.point(c.s)?.dist(b.point(c.r)?); + if gap <= tol { + if !found.iter().any(|f| (f.s - c.s).abs() < 1e-6 && (f.r - c.r).abs() < 1e-6) { + found.push(c); + } + } else if miss.is_none_or(|(g, _)| gap < g) { + miss = Some((gap, c.at)); + } + } + match found.as_slice() { + [c] => Ok(*c), + [] => { + let (gap, at) = miss.unwrap_or((f64::INFINITY, Vec3::ZERO)); + invalid(format!( + "{} and {} don't meet: they pass {gap:.4} mm apart near ({:.3}, {:.3}, {:.3}) (the crossing tolerance is {tol} mm)", + a.label, b.label, at.x, at.y, at.z + )) + } + _ => invalid(format!("{} and {} cross more than once", a.label, b.label)), + } +} + +/// The closest points of `a` and `b` near `(s, r)`: Gauss–Newton on `|a(s) − b(r)|²`, the +/// parameters kept in `[0, 1]`. +fn refine(a: &mut Track, b: &mut Track, mut s: f64, mut r: f64) -> SplineResult { + let h = 1e-7; + let slope = |t: &mut Track, x: f64| -> SplineResult { + let (x0, x1) = ((x - h).max(0.0), (x + h).min(1.0)); + Ok((t.point(x1)? - t.point(x0)?) / (x1 - x0)) + }; + for _ in 0..40 { + let (pa, pb) = (a.point(s)?, b.point(r)?); + let (da, db) = (slope(a, s)?, slope(b, r)?); + let e = pa - pb; + // Normal equations of the 3×2 system [da, −db]·(Δs, Δr) = −e. + let (m11, m12, m22) = (da.dot(da), -da.dot(db), db.dot(db)); + let (g1, g2) = (-da.dot(e), db.dot(e)); + let det = m11 * m22 - m12 * m12; + let (ds, dr) = if det.abs() > 1e-14 * (m11 * m22).max(1e-300) { + ((g1 * m22 - m12 * g2) / det, (m11 * g2 - m12 * g1) / det) + } else { + // Parallel there: a plain step along the gradient. + (g1 / m11.max(1e-300), g2 / m22.max(1e-300)) + }; + let (ns, nr) = ((s + ds).clamp(0.0, 1.0), (r + dr).clamp(0.0, 1.0)); + let done = (ns - s).abs() < 1e-13 && (nr - r).abs() < 1e-13; + (s, r) = (ns, nr); + if done { + break; + } + } + let (pa, pb) = (a.point(s)?, b.point(r)?); + Ok(Crossing { s, r, at: (pa + pb) * 0.5 }) +} + +/// What an interpolation condition across the network is matched with: curve `i`'s value, the +/// `order`-th cross derivative at the first (`end` 0) or last (`end` 1) curve, or zero (the +/// slope that holds a single curve's value constant). +#[derive(Clone, Copy, Debug, PartialEq)] +enum Datum { + Curve(usize), + Der { end: usize, order: usize }, + Zero, +} + +/// An interpolation operator across a direction: values at `params`, cross derivatives of the +/// orders in `ders` at the first and last, on `[0, 1]`. +struct Across { + interp: Interpolant, + data: Vec, +} + +impl Across { + fn new(params: &[f64], ders: [usize; 2]) -> SplineResult { + let mut conds = Vec::new(); + let mut data = Vec::new(); + for (i, &t) in params.iter().enumerate() { + conds.push(Cond { t, der: 0 }); + data.push(Datum::Curve(i)); + } + for (end, &count) in ders.iter().enumerate() { + let t = if end == 0 { params.first() } else { params.last() }; + let Some(&t) = t else { return invalid("no curves") }; + for order in 1..=count { + conds.push(Cond { t, der: order }); + data.push(Datum::Der { end, order }); + } + } + let interp = Interpolant::lowest(conds, 0.0, 1.0)?; + // A lone value gets a zero slope to hold it constant. + data.resize(interp.conds().len(), Datum::Zero); + Ok(Across { interp, data }) + } + + /// Interpolation of values given by `datum` (one per condition). + fn fit(&self, mut datum: impl FnMut(Datum) -> SplineResult) -> SplineResult> { + let values: Vec = self.data.iter().map(|d| if *d == Datum::Zero { Ok(T::zero()) } else { datum(*d) }).collect::>()?; + self.interp.fit(&values) + } + + fn space(&self) -> &Space { + self.interp.space() + } + + /// The interpolant's interior knots (the fitting space must contain them). + fn knots(&self) -> Vec { + let [a, b] = self.space().domain(); + self.space().knots().iter().copied().filter(|k| *k > a && *k < b).collect() + } +} + +/// A function along one direction that is fitted: its values (and, at the ends, derivatives) and +/// how close the fit must come. +struct Target<'f> { + label: String, + value: Box SplineResult + 'f>, + /// Derivatives (orders 1, 2, …) at the first and last parameter. + ders: [Vec; 2], + tol: f64, +} + +/// The common cubic space of a direction and each target's control points in it: breaks at +/// `base` (increasing, from 0 to 1) each split into equal parts, refined until every target is +/// within its tolerance at the quarter points of every span. A G2 end gets one more knot in its +/// end span, room for its second derivative. +fn fit_direction(base: &[f64], ders: [usize; 2], targets: &mut [Target]) -> SplineResult<(Space, Vec>)> { + let intervals = base.len().saturating_sub(1).max(1); + let mut m = MIN_SPANS.div_ceil(intervals); + loop { + let mut breaks = Vec::with_capacity(intervals * m + 1); + for w in base.windows(2) { + for k in 0..m { + breaks.push(if k == 0 { w[0] } else { w[0] + (w[1] - w[0]) * k as f64 / m as f64 }); + } + } + if let Some(&last) = base.last() { + breaks.push(last); + } + let spans = breaks.len() - 1; + let mut space = Space::with_breaks(3, &breaks)?; + let mid = |k: usize| (breaks[k] + breaks[k + 1]) / 2.0; + let mut extra = Vec::new(); + if ders[0] >= 2 { + extra.push(mid(0)); + } + if ders[1] >= 2 { + extra.push(mid(spans - 1)); + } + space = space.refined(&extra)?; + // Values at the breaks, the end derivatives, and a value mid-span at an end without one. + let mut conds: Vec = breaks.iter().map(|&t| Cond { t, der: 0 }).collect(); + let mut kinds: Vec<(usize, usize)> = breaks.iter().map(|_| (0, 0)).collect(); + for (end, &count) in ders.iter().enumerate() { + let t = if end == 0 { 0.0 } else { 1.0 }; + for order in 1..=count { + conds.push(Cond { t, der: order }); + kinds.push((end, order)); + } + if count == 0 { + conds.push(Cond { t: mid(if end == 0 { 0 } else { spans - 1 }), der: 0 }); + kinds.push((0, 0)); + } + } + let interp = Interpolant::new(space.clone(), conds.clone())?; + let mut ctrls = Vec::with_capacity(targets.len()); + let mut worst: Option<(f64, String)> = None; + for tg in targets.iter_mut() { + let data: Vec = conds + .iter() + .zip(&kinds) + .map(|(c, (end, order))| match order { + 0 => (tg.value)(c.t), + o => tg.ders.get(*end).and_then(|d| d.get(o - 1)).copied().ok_or_else(|| SplineError::Invalid("a derivative is missing".into())), + }) + .collect::>()?; + let ctrl = interp.fit(&data)?; + for k in 0..spans { + for q in [0.25, 0.75] { + let t = breaks[k] + (breaks[k + 1] - breaks[k]) * q; + let off = space.eval(&ctrl, t, 0)?.dist((tg.value)(t)?); + if off > tg.tol && worst.as_ref().is_none_or(|(w, _)| off / tg.tol > *w) { + worst = Some((off / tg.tol, format!("{} (off by {off:.2e} at the best fit)", tg.label))); + } + } + } + ctrls.push(ctrl); + } + match worst { + None => return Ok((space, ctrls)), + Some((_, what)) if spans * 2 > MAX_SPANS => { + return Err(SplineError::Unsupported(format!( + "{what} can't be fitted closely enough: a curve with a sharp corner, or one that turns very sharply, can't be part of a boundary surface" + ))); + } + Some(_) => m *= 2, + } + } +} + +/// Linear interpolation at `x` of `values` given at the increasing parameters `xs` (held +/// constant beyond the first and last; a lone value everywhere). +fn lerp_between(xs: &[f64], values: &[T], x: f64) -> T { + match (xs.len(), values.len()) { + (0, _) | (_, 0) => T::zero(), + (1, _) | (_, 1) => values[0], + (n, _) => { + let k = xs.partition_point(|t| *t <= x).saturating_sub(1).min(n - 2); + let (x0, x1) = (xs[k], xs[k + 1]); + let w = if x1 > x0 { ((x - x0) / (x1 - x0)).clamp(0.0, 1.0) } else { 0.0 }; + values[k].times(1.0 - w).plus(values.get(k + 1).copied().unwrap_or(values[k]).times(w)) + } + } +} + +/// How a curve runs over the network parameter: its arc-length position at network parameter +/// `x`, and the shift that puts it through the crossing points exactly. +struct Placement { + map: Monotone, + /// The network parameters of the crossings and each crossing's shift. + at: Vec, + shift: Vec, +} + +impl Placement { + fn point(&self, track: &mut Track, x: f64) -> SplineResult { + Ok(track.point(self.map.eval(x)[0])? + lerp_between(&self.at, &self.shift, x)) + } + + /// The first two derivatives by the network parameter at `x`. + fn ders(&self, track: &mut Track, x: f64) -> SplineResult<(Vec3, Vec3)> { + let [s, ds, dds] = self.map.eval(x); + let (p1, p2) = track.ders(s)?; + // The shift is linear between crossings: its slope on the span at x. + let slope = match self.at.len() { + 0 | 1 => Vec3::ZERO, + n => { + let k = self.at.partition_point(|t| *t <= x).saturating_sub(1).min(n - 2); + let (x0, x1) = (self.at[k], self.at[k + 1]); + if x1 > x0 { (self.shift[k + 1] - self.shift[k]) / (x1 - x0) } else { Vec3::ZERO } + } + }; + Ok((p1 * ds + slope, p2 * (ds * ds) + p1 * dds)) + } +} + +/// A placement through `nodes` (network parameter, arc-length position), with the shifts that +/// take the curve onto `targets` at the crossings `(x, s)` it was measured at. +fn placement(track: &mut Track, nodes: &[(f64, f64)], crossings: &[(f64, f64, Vec3)]) -> SplineResult { + let map = Monotone::new(nodes.iter().map(|n| n.0).collect(), nodes.iter().map(|n| n.1).collect()) + .map_err(|_| SplineError::Invalid(format!("{} is crossed out of order", track.label)))?; + let mut at = Vec::with_capacity(crossings.len()); + let mut shift = Vec::with_capacity(crossings.len()); + for &(x, s, q) in crossings { + at.push(x); + shift.push(q - track.point(s)?); + } + Ok(Placement { map, at, shift }) +} + +/// The frame of the face along conditioned curve `e`, with the tangent and outward directions +/// made orthonormal in its tangent plane. +struct Frame { + n: Vec3, + t: Vec3, + b: Vec3, + form: [f64; 3], +} + +impl Frame { + /// `II(w, w)` for `w = a·b + c·t`. + fn second(&self, a: f64, c: f64) -> f64 { + let [tt, tb, bb] = self.form; + c * c * tt + 2.0 * a * c * tb + a * a * bb + } +} + +/// A surface through a curve network (see the module documentation): curves `u` in one +/// direction, `v` (possibly none) in the other. +pub fn network_surface(u: Vec, v: Vec, tol: Tolerances) -> SplineResult { + if !(tol.cross > 0.0 && tol.cross.is_finite() && tol.fit > 0.0 && tol.fit.is_finite()) { + return invalid("the tolerances must be positive"); + } + if u.len() > MAX_NETWORK_CURVES || v.len() > MAX_NETWORK_CURVES { + return Err(SplineError::TooLarge(format!("more than {MAX_NETWORK_CURVES} curves in one direction"))); + } + let mut us: Vec = u.into_iter().map(Track::new).collect::>()?; + let mut vs: Vec = v.into_iter().map(Track::new).collect::>()?; + for t in us.iter().chain(&vs).filter(|t| !t.point) { + if t.pts.first().zip(t.pts.last()).is_some_and(|(a, b)| a.dist(*b) <= tol.cross) { + return Err(SplineError::Unsupported(format!("{} is closed: closed curves in a boundary surface", t.label))); + } + } + let conditioned = |ts: &[Track]| ts.iter().any(|t| t.condition.is_some()); + if conditioned(&us) && conditioned(&vs) { + return Err(SplineError::Unsupported("continuity conditions on curves of both directions (on two boundaries that meet at a corner)".into())); + } + // Conditions go on the u-curves, and there are two u-curves where possible. + if conditioned(&vs) || (us.len() < 2 && vs.len() >= 2 && !conditioned(&us)) { + std::mem::swap(&mut us, &mut vs); + } + let size = us + .iter() + .chain(&vs) + .flat_map(|t| t.pts.iter()) + .fold(None, |acc: Option<(Vec3, Vec3)>, p| Some(acc.map_or((*p, *p), |(lo, hi)| (lo.min(*p), hi.max(*p))))); + let size = size.map_or(0.0, |(lo, hi)| lo.dist(hi)); + match (us.len(), vs.len()) { + (0, _) | (1, 0) => return invalid("a boundary surface needs two curves in one direction, or crossing curves"), + (1, _) if !conditioned(&us) => return invalid("a boundary surface needs two curves in one direction"), + _ => {} + } + let (nu, nv) = (us.len(), vs.len()); + // A side collapsed to a point: one at most, with another curve beside it in its own direction, + // and the curves the other way (if any) running from it to the far side. + let points = |ts: &[Track]| ts.iter().filter(|t| t.point).count(); + match (points(&us), points(&vs)) { + (0, 0) => {} + (1, 0) if nu >= 2 && nv != 1 => {} + (0, 1) if nv >= 2 && nu >= 2 => {} + _ => return invalid("a point can stand for one side of the patch only, with curves running from it to the side across"), + } + + // Crossings: x[i][j] where u-curve i meets v-curve j (s along u-curve i, r along v-curve j). + // Where one of them is a point, the crossing is that point itself, so that every curve + // through it is moved onto it and the side stays collapsed. + let mut x: Vec> = Vec::with_capacity(nu); + for a in us.iter_mut() { + let mut row = Vec::with_capacity(nv); + for b in vs.iter_mut() { + let mut c = closest(a, b, tol.cross)?; + if a.point { + c.at = a.point(0.0)?; + } else if b.point { + c.at = b.point(0.0)?; + } + row.push(c); + } + x.push(row); + } + // The first curve of a direction that is not a point (orders are read along it). + let first_curve = |ts: &[Track]| ts.iter().position(|t| !t.point).unwrap_or(0); + + // Order and orient: the u-curves in their order along the first v-curve, each v-curve running + // from the first u-curve to the last; likewise the v-curves along the first u-curve. + if nv >= 1 && nu >= 2 { + let mut order: Vec = (0..nu).collect(); + let j0 = first_curve(&vs); + order.sort_by(|a, b| x[*a][j0].r.total_cmp(&x[*b][j0].r)); + us = reorder(us, &order); + x = reorder(x, &order); + for j in 0..nv { + // A point runs no way: nothing to orient. + if vs[j].point { + continue; + } + let rs: Vec = x.iter().map(|row| row[j].r).collect(); + match direction(&rs) { + Some(true) => {} + Some(false) => { + vs[j].flipped = !vs[j].flipped; + for row in x.iter_mut() { + row[j].r = 1.0 - row[j].r; + } + } + None => return invalid(format!("{} crosses the other curves out of order (or two of them at one point)", vs[j].label)), + } + } + } else if nv >= 1 { + // One conditioned u-curve: each v-curve starts on it. + for j in 0..nv { + if x[0][j].r > 0.5 { + vs[j].flipped = !vs[j].flipped; + x[0][j].r = 1.0 - x[0][j].r; + } + } + } + if nv >= 2 { + let mut order: Vec = (0..nv).collect(); + let i0 = first_curve(&us); + order.sort_by(|a, b| x[i0][*a].s.total_cmp(&x[i0][*b].s)); + vs = reorder(vs, &order); + for row in x.iter_mut() { + *row = reorder(std::mem::take(row), &order); + } + for (i, row) in x.iter_mut().enumerate() { + if us[i].point { + continue; + } + let ss: Vec = row.iter().map(|c| c.s).collect(); + match direction(&ss) { + Some(true) => {} + Some(false) => { + us[i].flipped = !us[i].flipped; + for c in row.iter_mut() { + c.s = 1.0 - c.s; + } + } + None => return invalid(format!("{} crosses the other curves out of order (or two of them at one point)", us[i].label)), + } + } + } else { + // No order along the u-curves: each runs the way the one before it does. + for i in 1..nu { + if us[i].point || us[i - 1].point { + continue; + } + let (ps, pe) = (us[i - 1].point(0.0)?, us[i - 1].point(1.0)?); + let (s, e) = (us[i].point(0.0)?, us[i].point(1.0)?); + if s.dist(ps) + e.dist(pe) > s.dist(pe) + e.dist(ps) { + us[i].flipped = !us[i].flipped; + for c in x[i].iter_mut() { + c.s = 1.0 - c.s; + } + } + } + } + + // A point is a side: the first or last curve of its direction. + for ts in [&us, &vs] { + for (i, t) in ts.iter().enumerate() { + if t.point && i != 0 && i + 1 != ts.len() { + return invalid(format!("{} lies inside the network: a point can only stand for a side", t.label)); + } + } + } + // Conditions only on the first and last u-curve. + for (i, t) in us.iter().enumerate() { + if t.condition.is_some() && i != 0 && i + 1 != nu { + return invalid(format!("{} lies inside the network: continuity applies only to a boundary curve", t.label)); + } + } + let end_order = [us.first().map_or(0, |t| t.continuity().order()), if nu >= 2 { us.last().map_or(0, |t| t.continuity().order()) } else { 0 }]; + + // Network parameters (a point has no positions along it to average). + let rows: Vec> = x.iter().zip(&us).filter(|(_, t)| !t.point).map(|(row, _)| row.iter().map(|c| c.s).collect()).collect(); + let ubar = network_params(&rows, nv)?; + let vbar: Vec = if nv == 0 { + chord_params(&mut us)? + } else if nu >= 2 { + let cols: Vec> = (0..nv).filter(|j| !vs[*j].point).map(|j| x.iter().map(|row| row[j].r).collect()).collect(); + network_params(&cols, nu)? + } else { + // One conditioned curve, crossed at the start of every v-curve. + for j in 0..nv { + if x[0][j].r * vs[j].length > tol.cross { + return invalid(format!("{} must start on {}: it crosses it partway along", vs[j].label, us[0].label)); + } + } + vec![0.0] + }; + + // A point is placed evenly: its crossings sit at the network parameters themselves. + for (i, row) in x.iter_mut().enumerate() { + for (j, c) in row.iter_mut().enumerate() { + if us.get(i).is_some_and(|t| t.point) { + c.s = ubar.get(j).copied().unwrap_or(c.s); + } + if vs.get(j).is_some_and(|t| t.point) { + c.r = vbar.get(i).copied().unwrap_or(c.r); + } + } + } + + // Placements: each curve through its crossings at the network parameters. + let mut up = Vec::with_capacity(nu); + for (i, t) in us.iter_mut().enumerate() { + let cross: Vec<(f64, f64, Vec3)> = (0..nv).map(|j| (ubar[j], x[i][j].s, x[i][j].at)).collect(); + let nodes = place_nodes(&cross, nv); + up.push(placement(t, &nodes, &cross)?); + } + let mut vp = Vec::with_capacity(nv); + for (j, t) in vs.iter_mut().enumerate() { + let cross: Vec<(f64, f64, Vec3)> = (0..nu).map(|i| (vbar[i], x[i][j].r, x[i][j].at)).collect(); + let nodes = place_nodes(&cross, nu); + vp.push(placement(t, &nodes, &cross)?); + } + + // The operators across: Π_v over the u-curves (with the end conditions), Π_u over the v-curves. + let across_v = Across::new(&vbar, end_order)?; + let across_u = if nv >= 1 { Some(Across::new(&ubar, [0, 0])?) } else { None }; + + // Cross-boundary conditions: for each conditioned end, the frame along it, the natural + // derivatives, and the corrections that make them agree with the v-curves at the crossings. + let mut ends: Vec> = vec![None, None]; + for end in 0..2 { + if end_order[end] == 0 { + continue; + } + let i = if end == 0 { 0 } else { nu - 1 }; + let mut cond = EndCondition { end, curve: i, order: end_order[end], neighbour: None, at: Vec::new(), corr: [None, None, None, None] }; + if nu >= 2 { + let k = if end == 0 { 1 } else { nu - 2 }; + cond.neighbour = Some((k, (vbar[k] - vbar[i]).abs())); + } + // The v-curves where they meet this boundary: their derivatives, checked against the face. + let mut values: Vec<[f64; 4]> = Vec::with_capacity(nv); + for j in 0..nv { + let xj = ubar[j]; + let frame = end_frame(&mut us, &up, i, xj)?; + let (d1, d2) = vp[j].ders(&mut vs[j], vbar[i])?; + let (vl, ul) = (vs[j].label.clone(), us[i].label.clone()); + let len = d1.len(); + if !(len > 0.0) { + return invalid(format!("{vl} has no direction where it meets {ul}")); + } + if d1.dot(frame.n).abs() > TANGENT_TOL.sin() * len { + let deg = (d1.dot(frame.n).abs() / len).clamp(0.0, 1.0).asin().to_degrees(); + return invalid(format!( + "{vl} meets the {} boundary {ul} at {deg:.2}° to its face: it must leave the face tangent to it there", + if end_order[end] == 1 { "G1" } else { "G2" } + )); + } + let a = d1.dot(frame.b); + let outward = if end == 0 { a } else { -a }; + if outward <= TANGENT_TOL.sin() * len { + return invalid(format!("{vl} turns back over the face of {ul} where it meets it")); + } + let c = d1.dot(frame.t); + let (mut kb, mut kt) = (0.0, 0.0); + if end_order[end] >= 2 { + let tproj = frame.b * a + frame.t * c; + let want = frame.second(a, c) / tproj.len2(); + let got = d2.dot(frame.n) / d1.len2(); + let scale = want.abs().max(1.0 / size.max(1e-9)); + if (got - want).abs() > CURVATURE_TOL * scale { + return invalid(format!( + "{vl} meets the G2 boundary {ul} with a normal curvature of {got:.5} /mm where the face has {want:.5} /mm: it must follow the face's curvature there" + )); + } + kb = d2.dot(frame.b); + kt = d2.dot(frame.t); + } + values.push([a, c, kb, kt]); + } + cond.at = values; + ends[end] = Some(cond); + } + // Corrections: Π_u through (crossing value − natural value), so each function hits the + // v-curves' own derivatives at the crossings. + if let Some(au) = &across_u { + for slot in ends.iter_mut() { + let Some(mut cond) = slot.take() else { continue }; + let mut diffs: [Vec; 4] = [Vec::new(), Vec::new(), Vec::new(), Vec::new()]; + for (&xj, at) in ubar.iter().zip(&cond.at) { + let frame = end_frame(&mut us, &up, cond.curve, xj)?; + let natural = cond.natural(&mut us, &up, &frame, xj)?; + for (k, d) in diffs.iter_mut().enumerate() { + d.push(at[k] - natural[k]); + } + } + for (k, d) in diffs.iter().enumerate() { + let ctrl = au.fit(|datum| match datum { + Datum::Curve(j) => d.get(j).copied().ok_or_else(|| SplineError::Invalid("crossing".into())), + _ => Ok(0.0), + })?; + cond.corr[k] = Some(ctrl); + } + *slot = Some(cond); + } + } + + // Fit the u-direction: the u-curves, then each condition's first and second derivatives. + let mut base_u: Vec = vec![0.0, 1.0]; + base_u.extend(ubar.iter().copied()); + if let Some(au) = &across_u { + base_u.extend(au.knots()); + } + let base_u = sorted_breaks(base_u); + let (u_space, u_ctrl) = { + let us_cell = std::cell::RefCell::new(&mut us); + let up_ref = &up; + let mut targets: Vec = Vec::new(); + for i in 0..nu { + let label = us_cell.borrow()[i].label.clone(); + let cell = &us_cell; + targets.push(Target { + label, + value: Box::new(move |t| { + let mut tracks = cell.borrow_mut(); + up_ref[i].point(&mut tracks[i], t) + }), + ders: [Vec::new(), Vec::new()], + tol: tol.fit, + }); + } + for cond in ends.iter().flatten() { + let label = format!("the cross-boundary tangent along {}", us_cell.borrow()[cond.curve].label); + let scale = cond.scale(&mut us_cell.borrow_mut(), up_ref, across_u.as_ref())?; + for order in 1..=cond.order { + let cell = &us_cell; + let au = across_u.as_ref(); + targets.push(Target { + label: label.clone(), + value: Box::new(move |t| { + let mut tracks = cell.borrow_mut(); + let (d1, d2) = cond.derivatives(&mut tracks, up_ref, au, t)?; + Ok(if order == 1 { d1 } else { d2 }) + }), + ders: [Vec::new(), Vec::new()], + tol: scale[order - 1] * 1e-4, + }); + } + } + fit_direction(&base_u, [0, 0], &mut targets)? + }; + // Which fitted function is which: the u-curves, then per conditioned end its derivatives. + let mut der_index: [[Option; 2]; 2] = [[None; 2]; 2]; + let mut next = nu; + for cond in ends.iter().flatten() { + for order in 1..=cond.order { + der_index[cond.end][order - 1] = Some(next); + next += 1; + } + } + let u_fn = |d: Datum, a: usize| -> SplineResult { + let k = match d { + Datum::Curve(i) => Some(i), + Datum::Der { end, order } => der_index.get(end).and_then(|e| e.get(order.wrapping_sub(1))).copied().flatten(), + Datum::Zero => None, + }; + Ok(k.and_then(|k| u_ctrl.get(k)).and_then(|c| c.get(a)).copied().unwrap_or(Vec3::ZERO)) + }; + // P_u: across v, column by column of the u-direction control points. + let nu_ctrl = u_space.count(); + let mut pu: Vec> = Vec::with_capacity(nu_ctrl); + for a in 0..nu_ctrl { + pu.push(across_v.fit(|d| u_fn(d, a))?); + } + let v_op_space = across_v.space().clone(); + let Some(across_u) = across_u else { + // Curves in one direction only: the skinned surface. + return surface(&u_space, &v_op_space, &pu); + }; + + // The derivative data at each crossing of a conditioned boundary, exactly as the fitted + // functions take them there. + let mut cross_ders: Vec<[[Vec3; 2]; 2]> = vec![[[Vec3::ZERO; 2]; 2]; nv]; + for cond in ends.iter().flatten() { + for (j, cd) in cross_ders.iter_mut().enumerate() { + let (d1, d2) = cond.derivatives(&mut us, &up, Some(&across_u), ubar[j])?; + cd[cond.end] = [d1, d2]; + } + } + + // Fit the v-direction: the v-curves, with the conditioned ends' derivatives. + let mut base_v: Vec = vec![0.0, 1.0]; + base_v.extend(vbar.iter().copied()); + base_v.extend(across_v.knots()); + let base_v = sorted_breaks(base_v); + let (v_space, v_ctrl) = { + let vs_cell = std::cell::RefCell::new(&mut vs); + let vp_ref = &vp; + let mut targets: Vec = Vec::new(); + for (j, cd) in cross_ders.iter().enumerate() { + let label = vs_cell.borrow()[j].label.clone(); + let cell = &vs_cell; + let ders = [cd[0][..end_order[0]].to_vec(), cd[1][..end_order[1]].to_vec()]; + targets.push(Target { + label, + value: Box::new(move |t| { + let mut tracks = cell.borrow_mut(); + vp_ref[j].point(&mut tracks[j], t) + }), + ders, + tol: tol.fit, + }); + } + fit_direction(&base_v, end_order, &mut targets)? + }; + let nv_ctrl = v_space.count(); + // P_v: across u, row by row of the v-direction control points. + let mut pv_rows: Vec> = Vec::with_capacity(nv_ctrl); + for b in 0..nv_ctrl { + pv_rows.push(across_u.fit(|d| match d { + Datum::Curve(j) => v_ctrl.get(j).and_then(|c| c.get(b)).copied().ok_or_else(|| SplineError::Invalid("v-curve".into())), + _ => Ok(Vec3::ZERO), + })?); + } + // P_uv: the crossing points (and cross derivatives) across v, then across u. + let mut cols: Vec> = Vec::with_capacity(nv); + for (j, cd) in cross_ders.iter().enumerate() { + cols.push(across_v.fit(|d| match d { + Datum::Curve(i) => x.get(i).and_then(|row| row.get(j)).map(|c| c.at).ok_or_else(|| SplineError::Invalid("crossing".into())), + Datum::Der { end, order } => Ok(cd.get(end).and_then(|e| e.get(order.wrapping_sub(1))).copied().unwrap_or(Vec3::ZERO)), + Datum::Zero => Ok(Vec3::ZERO), + })?); + } + let nvop = v_op_space.count(); + let mut puv_rows: Vec> = Vec::with_capacity(nvop); + for b in 0..nvop { + puv_rows.push(across_u.fit(|d| match d { + Datum::Curve(j) => cols.get(j).and_then(|c| c.get(b)).copied().ok_or_else(|| SplineError::Invalid("crossing".into())), + _ => Ok(Vec3::ZERO), + })?); + } + // Into the common space KU × KV and summed. + let u_op_space = across_u.space().clone(); + let pu_full: Vec> = pu.iter().map(|col| embed(&v_op_space, col, &v_space)).collect::>()?; + let pv_cols = transpose(&pv_rows.iter().map(|row| embed(&u_op_space, row, &u_space)).collect::>>()?); + let puv_rows_u: Vec> = puv_rows.iter().map(|row| embed(&u_op_space, row, &u_space)).collect::>()?; + let puv_cols: Vec> = transpose(&puv_rows_u).iter().map(|col| embed(&v_op_space, col, &v_space)).collect::>()?; + let mut net = Vec::with_capacity(nu_ctrl); + for a in 0..nu_ctrl { + let col: Vec = (0..nv_ctrl) + .map(|b| { + let get = |m: &Vec>| m.get(a).and_then(|c| c.get(b)).copied().unwrap_or(Vec3::ZERO); + get(&pu_full) + get(&pv_cols) - get(&puv_cols) + }) + .collect(); + net.push(col); + } + surface(&u_space, &v_space, &net) +} + +/// A conditioned end of the u-direction: which curve, the derivative orders, its neighbouring +/// curve (and the parameter gap to it) for the natural tangent, the scalar values the v-curves +/// give at each crossing (`[a, c, kb, kt]`: the tangent's outward and along-curve components, the +/// second derivative's), and the corrections to the natural values (control values across u). +struct EndCondition { + end: usize, + curve: usize, + order: usize, + neighbour: Option<(usize, f64)>, + at: Vec<[f64; 4]>, + corr: [Option>; 4], +} + +impl EndCondition { + /// The natural `[a, c, kb, kt]` at `x`: the tangent as long as the step to the neighbouring + /// curve, outward from the face (into it at the far end), sheared along the curve as that + /// step is; no tangential second derivative. + fn natural(&self, us: &mut [Track], up: &[Placement], frame: &Frame, x: f64) -> SplineResult<[f64; 4]> { + let Some((k, gap)) = self.neighbour else { return Ok([0.0; 4]) }; + let p = up.get(self.curve).ok_or_else(|| SplineError::Invalid("curve".into()))?.point(track(us, self.curve)?, x)?; + let q = up.get(k).ok_or_else(|| SplineError::Invalid("curve".into()))?.point(track(us, k)?, x)?; + let step = if self.end == 0 { q - p } else { p - q }; + let gap = gap.max(1e-12); + let a = step.len() / gap * if self.end == 0 { 1.0 } else { -1.0 }; + Ok([a, step.dot(frame.t) / gap, 0.0, 0.0]) + } + + /// The first and second v-derivatives of the surface along this boundary at `x`. + fn derivatives(&self, us: &mut [Track], up: &[Placement], across_u: Option<&Across>, x: f64) -> SplineResult<(Vec3, Vec3)> { + self.derivatives_framed(us, up, across_u, x).map(|(d1, d2, _)| (d1, d2)) + } + + /// [`EndCondition::derivatives`] with the face's frame there (projecting onto the face is + /// the costly part: it is done once). + fn derivatives_framed(&self, us: &mut [Track], up: &[Placement], across_u: Option<&Across>, x: f64) -> SplineResult<(Vec3, Vec3, Frame)> { + let frame = end_frame(us, up, self.curve, x)?; + let mut v = self.natural(us, up, &frame, x)?; + if let Some(au) = across_u { + for (k, c) in self.corr.iter().enumerate() { + if let Some(c) = c { + v[k] += au.space().eval(c, x, 0)?; + } + } + } + let [a, c, kb, kt] = v; + let d1 = frame.b * a + frame.t * c; + let d2 = frame.n * frame.second(a, c) + frame.b * kb + frame.t * kt; + Ok((d1, d2, frame)) + } + + /// Typical sizes of the first and second derivatives along the boundary (for fitting + /// tolerances), checking that the tangent never turns back along it. + fn scale(&self, us: &mut [Track], up: &[Placement], across_u: Option<&Across>) -> SplineResult<[f64; 2]> { + let mut big = [0.0f64; 2]; + for k in 0..=32 { + let x = k as f64 / 32.0; + let (d1, d2, frame) = self.derivatives_framed(us, up, across_u, x)?; + let out = d1.dot(frame.b) * if self.end == 0 { 1.0 } else { -1.0 }; + if !(out > 0.0) { + let label = us.get(self.curve).map(|t| t.label.clone()).unwrap_or_default(); + return invalid(format!("the surface's tangent across {label} would turn back over its face")); + } + big = [big[0].max(d1.len()), big[1].max(d2.len())]; + } + Ok([big[0].max(1e-9), big[1].max(big[0] * 1e-3).max(1e-9)]) + } +} + +fn track<'t, 'a>(us: &'t mut [Track<'a>], i: usize) -> SplineResult<&'t mut Track<'a>> { + us.get_mut(i).ok_or_else(|| SplineError::Invalid("curve".into())) +} + +/// The face's frame along conditioned u-curve `i` at network parameter `x`. +fn end_frame(us: &mut [Track], up: &[Placement], i: usize, x: f64) -> SplineResult { + let place = up.get(i).ok_or_else(|| SplineError::Invalid("curve".into()))?; + let t = track(us, i)?; + let p = place.point(t, x)?; + let (d1, _) = place.ders(t, x)?; + let label = t.label.clone(); + let face = t.condition.as_mut().map(|(_, f)| f).ok_or_else(|| SplineError::Invalid(format!("{label} has no face")))?; + let f = face + .frame(p, d1) + .ok_or_else(|| SplineError::Invalid(format!("the face along {label} can't be evaluated at ({:.3}, {:.3}, {:.3})", p.x, p.y, p.z)))?; + let n = f.normal.normalized().ok_or_else(|| SplineError::Invalid(format!("the face along {label} has no normal")))?; + let t = (d1 - n * d1.dot(n)).normalized().ok_or_else(|| SplineError::Invalid(format!("{label} runs along its face's normal")))?; + let mut b = f.outward - n * f.outward.dot(n); + b = (b - t * b.dot(t)).normalized().ok_or_else(|| SplineError::Invalid(format!("the face along {label} has no outward direction")))?; + Ok(Frame { n, t, b, form: f.form }) +} + +/// `true` when `xs` increase strictly, `false` when they decrease strictly, `None` otherwise. +fn direction(xs: &[f64]) -> Option { + let gap = 1e-9; + if xs.windows(2).all(|w| w[1] - w[0] > gap) { + Some(true) + } else if xs.windows(2).all(|w| w[0] - w[1] > gap) { + Some(false) + } else { + None + } +} + +fn reorder(items: Vec, order: &[usize]) -> Vec { + let mut slots: Vec> = items.into_iter().map(Some).collect(); + order.iter().filter_map(|&k| slots.get_mut(k).and_then(Option::take)).collect() +} + +fn transpose(m: &[Vec]) -> Vec> { + let cols = m.first().map_or(0, Vec::len); + (0..cols).map(|c| m.iter().map(|row| row.get(c).copied().unwrap_or(Vec3::ZERO)).collect()).collect() +} + +fn sorted_breaks(mut b: Vec) -> Vec { + b.retain(|x| x.is_finite() && (0.0..=1.0).contains(x)); + b.sort_by(f64::total_cmp); + b.dedup_by(|a, b| (*a - *b).abs() < 1e-12); + b +} + +/// The network parameters of the crossings along one direction: `rows[i][j]` is where curve `j` +/// of the other direction crosses curve `i` (arc length). With two crossings or more each curve +/// runs from its first to its last crossing; the parameters are the relative positions, +/// averaged over the curves. With one, it is its mean position along the whole curves. +fn network_params(rows: &[Vec], count: usize) -> SplineResult> { + let n = rows.len().max(1) as f64; + match count { + 0 => Ok(Vec::new()), + 1 => { + let firsts: Vec = rows.iter().filter_map(|r| r.first().copied()).collect(); + let at_end = |x: f64| x <= 1e-9 || x >= 1.0 - 1e-9; + if firsts.iter().any(|x| at_end(*x)) && !firsts.iter().all(|x| at_end(*x)) { + return invalid("a curve crosses the others at the end of some and partway along others"); + } + Ok(vec![firsts.iter().sum::() / n]) + } + _ => { + let mut out = vec![0.0; count]; + for r in rows { + let (Some(&a), Some(&b)) = (r.first(), r.last()) else { return invalid("crossings") }; + if !(b - a > 1e-12) { + return invalid("the curves cross out of order"); + } + for (o, x) in out.iter_mut().zip(r) { + *o += (x - a) / (b - a) / n; + } + } + out[0] = 0.0; + out[count - 1] = 1.0; + if out.windows(2).any(|w| !(w[1] - w[0] > 1e-9)) { + return invalid("two curves cross the others at the same place"); + } + Ok(out) + } + } +} + +/// Lofting parameters of curves in one direction only: cumulative mean distances between +/// neighbouring curves, from 0 to 1. +fn chord_params(us: &mut [Track]) -> SplineResult> { + let mut out = vec![0.0]; + let mut acc = 0.0; + for i in 1..us.len() { + let mut d = 0.0; + for k in 0..=32 { + let s = k as f64 / 32.0; + let p = track(us, i - 1)?.point(s)?; + let q = track(us, i)?.point(s)?; + d += p.dist(q) / 33.0; + } + if !(d > 1e-9) { + let (a, b) = (us[i - 1].label.clone(), us[i].label.clone()); + return invalid(format!("{a} and {b} lie on each other")); + } + acc += d; + out.push(acc); + } + Ok(out.into_iter().map(|x| x / acc).collect()) +} + +/// The nodes `(network parameter, arc length)` a curve is placed through: its crossings, and +/// with fewer than two the curve's ends at 0 and 1 (a crossing at an end standing for that end). +fn place_nodes(cross: &[(f64, f64, Vec3)], count: usize) -> Vec<(f64, f64)> { + if count >= 2 { + return cross.iter().map(|c| (c.0, c.1)).collect(); + } + let mut nodes = vec![(0.0, 0.0)]; + if let Some(&(x, s, _)) = cross.first() + && x > 1e-9 + && x < 1.0 - 1e-9 + && s > 1e-9 + && s < 1.0 - 1e-9 + { + nodes.push((x, s)); + } + nodes.push((1.0, 1.0)); + nodes +} + +/// The tensor surface of a net of control points (`net[a][b]`: a along u). +fn surface(su: &Space, sv: &Space, net: &[Vec]) -> SplineResult { + let (nu, nv) = (su.count(), sv.count()); + if net.len() != nu || net.iter().any(|c| c.len() != nv) { + return invalid("control net size"); + } + let mut pts = Vec::with_capacity(nu * nv); + for b in 0..nv { + for col in net { + pts.push(col[b]); + } + } + let weights = vec![1.0; pts.len()]; + NurbsSurface::new([su.degree(), sv.degree()], [su.knots().to_vec(), sv.knots().to_vec()], &pts, &weights) +} + +#[cfg(test)] +mod tests { + use super::*; + + /// A curve given by a function of its parameter. + struct F Vec3>(G); + impl Vec3> NetCurve for F { + fn point(&mut self, t: f64) -> Option { + Some((self.0)(t)) + } + } + + fn curve<'a>(label: &str, f: impl FnMut(f64) -> Vec3 + 'a) -> NetworkCurve<'a> { + NetworkCurve { label: label.into(), curve: Box::new(F(f)), condition: None } + } + + fn line<'a>(label: &str, a: Vec3, b: Vec3) -> NetworkCurve<'a> { + curve(label, move |t| a.lerp(b, t)) + } + + const TOL: Tolerances = Tolerances { cross: 1e-3, fit: 1e-5 }; + + fn grid(s: &NurbsSurface, n: usize) -> Vec { + let [du, dv] = s.domain(); + let mut out = Vec::new(); + for i in 0..=n { + for j in 0..=n { + let (u, v) = (du[0] + (du[1] - du[0]) * i as f64 / n as f64, dv[0] + (dv[1] - dv[0]) * j as f64 / n as f64); + out.push(s.eval(u, v).unwrap()); + } + } + out + } + + /// Distance from `p` to the polyline of `f` sampled finely. + fn dist_to(p: Vec3, f: &dyn Fn(f64) -> Vec3) -> f64 { + (0..2000).map(|k| p.dist_to_segment(f(k as f64 / 2000.0), f((k + 1) as f64 / 2000.0))).fold(f64::INFINITY, f64::min) + } + + /// The bilinearly blended Coons patch of four straight lines reproduces the hyperbolic + /// paraboloid z = x·y they bound, exactly; four coplanar curves (two arcs, two lines) give a + /// patch in their plane. + #[test] + fn coons_patch_reproduces_a_hyperbolic_paraboloid_and_a_plane() { + let p = |x: f64, y: f64| Vec3::new(x, y, x * y); + let s = network_surface( + vec![line("u0", p(0.0, 0.0), p(4.0, 0.0)), line("u1", p(0.0, 3.0), p(4.0, 3.0))], + vec![line("v0", p(0.0, 0.0), p(0.0, 3.0)), line("v1", p(4.0, 3.0), p(4.0, 0.0))], + TOL, + ) + .unwrap(); + for q in grid(&s, 20) { + assert!((q.z - q.x * q.y).abs() < 1e-9, "{q:?} off z = xy by {}", q.z - q.x * q.y); + } + assert!(s.eval(0.5, 0.5).unwrap().dist(p(2.0, 1.5)) < 1e-9, "the middle of the patch"); + // A plane z = 0.3x + 0.2y + 1 bounded by two arcs bulging out and two lines. + let on = |x: f64, y: f64| Vec3::new(x, y, 0.3 * x + 0.2 * y + 1.0); + let arc = |y: f64, sign: f64| move |t: f64| on(10.0 * t, y + sign * 2.0 * (std::f64::consts::PI * t).sin()); + let s = network_surface( + vec![curve("bottom", arc(0.0, -1.0)), curve("top", arc(8.0, 1.0))], + vec![line("left", on(0.0, 0.0), on(0.0, 8.0)), line("right", on(10.0, 0.0), on(10.0, 8.0))], + TOL, + ) + .unwrap(); + for q in grid(&s, 20) { + assert!((q.z - (0.3 * q.x + 0.2 * q.y + 1.0)).abs() < 1e-9, "{q:?} off the plane"); + } + } + + /// Area of a surface over its whole domain: |S_u × S_v| by central differences on an n × n grid + /// of cell centres. + fn area(s: &NurbsSurface, n: usize) -> f64 { + let [du, dv] = s.domain(); + let (hu, hv) = ((du[1] - du[0]) / n as f64, (dv[1] - dv[0]) / n as f64); + let mut a = 0.0; + for i in 0..n { + for j in 0..n { + let (u, v) = (du[0] + hu * (i as f64 + 0.5), dv[0] + hv * (j as f64 + 0.5)); + let su = (s.eval(u + hu * 0.25, v).unwrap() - s.eval(u - hu * 0.25, v).unwrap()) / (hu * 0.5); + let sv = (s.eval(u, v + hv * 0.25).unwrap() - s.eval(u, v - hv * 0.25).unwrap()) / (hv * 0.5); + a += su.cross(sv).len() * hu * hv; + } + } + a + } + + fn apex<'a>(label: &str, p: Vec3) -> NetworkCurve<'a> { + NetworkCurve { label: label.into(), curve: Box::new(PointCurve(p)), condition: None } + } + + /// Three sides: one of them a point. A triangle of lines and a quarter disc (two radii and the + /// arc, the apex at the centre) are filled flat, the whole patch on their plane, every point of + /// the collapsed side at the apex, with the area of the triangle and of the quarter disc. + #[test] + fn three_sides_collapse_one_to_a_point() { + let (a, b, c) = (Vec3::new(0.0, 0.0, 2.0), Vec3::new(10.0, 0.0, 2.0), Vec3::new(4.0, 8.0, 2.0)); + let s = network_surface(vec![line("ab", a, b), apex("c", c)], vec![line("ac", a, c), line("cb", c, b)], TOL).unwrap(); + for q in grid(&s, 20) { + assert!((q.z - 2.0).abs() < 1e-9, "{q:?} off the plane"); + } + let [du, dv] = s.domain(); + let collapsed = [s.eval(du[0], dv[0]).unwrap(), s.eval(du[1], dv[0]).unwrap(), s.eval(du[0], dv[1]).unwrap(), s.eval(du[1], dv[1]).unwrap()]; + assert!(collapsed.iter().filter(|p| p.dist(c) < 1e-12).count() == 2, "two corners at the apex: {collapsed:?}"); + let tri = area(&s, 200); + assert!((tri - 40.0).abs() < 1e-3, "triangle area {tri}"); + // The quarter disc of radius 10 about the origin. + let r = 10.0; + let arc = move |t: f64| Vec3::new(r * (std::f64::consts::FRAC_PI_2 * t).cos(), r * (std::f64::consts::FRAC_PI_2 * t).sin(), 0.0); + let s = network_surface( + vec![curve("arc", arc), apex("centre", Vec3::ZERO)], + vec![line("x", Vec3::ZERO, Vec3::X * r), line("y", Vec3::Y * r, Vec3::ZERO)], + TOL, + ) + .unwrap(); + for q in grid(&s, 20) { + assert!(q.z.abs() < 1e-9 && q.len() <= r + 1e-4, "{q:?} off the quarter disc"); + } + let quarter = area(&s, 200); + let want = std::f64::consts::PI * r * r / 4.0; + assert!((quarter - want).abs() < 2e-3 * want, "quarter disc area {quarter} vs {want}"); + // A point given as a v-curve is the same patch. + let s2 = network_surface( + vec![line("x", Vec3::ZERO, Vec3::X * r), line("y", Vec3::Y * r, Vec3::ZERO)], + vec![curve("arc", arc), apex("centre", Vec3::ZERO)], + TOL, + ) + .unwrap(); + assert!((area(&s2, 200) - quarter).abs() < 1e-6); + } + + /// Points where a patch has no side for them, a point that moves, and a point with a + /// continuity condition are refused, never a panic. + #[test] + fn misplaced_points_are_refused() { + let (a, b, c) = (Vec3::ZERO, Vec3::new(10.0, 0.0, 0.0), Vec3::new(4.0, 8.0, 0.0)); + let two = network_surface(vec![apex("p", a), apex("q", c)], vec![], TOL).unwrap_err(); + assert!(two.to_string().contains("point"), "{two}"); + let both = network_surface(vec![line("ab", a, b), apex("c", c)], vec![apex("a", a), line("cb", c, b)], TOL).unwrap_err(); + assert!(both.to_string().contains("point"), "{both}"); + let lone = network_surface(vec![apex("c", c)], vec![line("ac", a, c), line("cb", c, b)], TOL); + assert!(lone.is_err()); + let one_across = network_surface(vec![line("ab", a, b), apex("c", c)], vec![line("ac", a, c)], TOL).unwrap_err(); + assert!(one_across.to_string().contains("point"), "{one_across}"); + // A point between two curves of its direction. + let mid = network_surface( + vec![line("low", a, b), apex("mid", Vec3::new(5.0, 4.0, 0.0)), line("high", Vec3::new(0.0, 8.0, 0.0), Vec3::new(10.0, 8.0, 0.0))], + vec![], + TOL, + ) + .unwrap_err(); + assert!(mid.to_string().contains("inside"), "{mid}"); + struct Moving; + impl NetCurve for Moving { + fn point(&mut self, t: f64) -> Option { + Some(Vec3::new(t, 0.0, 0.0)) + } + fn is_point(&self) -> bool { + true + } + } + let moving = NetworkCurve { label: "m".into(), curve: Box::new(Moving), condition: None }; + let e = network_surface(vec![line("ab", a, b), moving], vec![], TOL).unwrap_err(); + assert!(e.to_string().contains("moves"), "{e}"); + let with_g1 = NetworkCurve { label: "g".into(), curve: Box::new(PointCurve(c)), condition: Some((Continuity::G1, Box::new(Plane))) }; + let e = network_surface(vec![line("ab", a, b), with_g1], vec![], TOL).unwrap_err(); + assert!(e.to_string().contains("point"), "{e}"); + } + + /// Two curves loft straight between them; three or more pass smoothly through all; curves + /// given running opposite ways are turned to agree. + #[test] + fn curves_one_way_are_lofted() { + let ruled = + network_surface(vec![line("a", Vec3::ZERO, Vec3::X * 10.0), line("b", Vec3::new(10.0, 5.0, 2.0), Vec3::new(0.0, 5.0, 2.0))], vec![], TOL) + .unwrap(); + assert_eq!(ruled.degree()[1], 1, "ruled between two curves"); + for q in grid(&ruled, 10) { + assert!((q.z - 0.4 * q.y).abs() < 1e-9, "{q:?}: a flat strip, not twisted"); + } + let arc = |z: f64, r: f64| move |t: f64| Vec3::new(r * (std::f64::consts::PI * t).cos(), r * (std::f64::consts::PI * t).sin(), z); + let fs: Vec<(f64, f64)> = vec![(0.0, 5.0), (4.0, 6.0), (8.0, 5.5), (12.0, 7.0)]; + let s = network_surface(fs.iter().enumerate().map(|(k, &(z, r))| curve(&format!("c{k}"), arc(z, r))).collect(), vec![], TOL).unwrap(); + // Every section lies on the surface: its iso-curve at the section's parameter is close. + let pts = grid(&s, 60); + for &(z, r) in &fs { + let f = arc(z, r); + for k in 0..=10 { + let p = f(k as f64 / 10.0); + let d = pts.iter().map(|q| q.dist(p)).fold(f64::INFINITY, f64::min); + assert!(d < 0.3, "section at z = {z} is {d} from the surface's samples"); + } + } + } + + /// A 3 × 3 network: the Gordon surface passes through all six curves (within the fit + /// tolerance), the interior ones included. + #[test] + fn gordon_surface_passes_through_every_curve() { + let z = |x: f64, y: f64| 0.05 * x * x - 0.03 * y * y + 0.02 * x * y; + let ux: Vec = vec![0.0, 4.0, 10.0]; + let vy: Vec = vec![0.0, 3.0, 8.0]; + let ucurve = |y: f64| move |t: f64| Vec3::new(10.0 * t, y, z(10.0 * t, y)); + let vcurve = |x: f64| move |t: f64| Vec3::new(x, 8.0 * t, z(x, 8.0 * t)); + let s = network_surface( + vy.iter().map(|&y| curve(&format!("y={y}"), ucurve(y))).collect(), + ux.iter().rev().map(|&x| curve(&format!("x={x}"), vcurve(x))).collect(), + TOL, + ) + .unwrap(); + // The surface's iso-curves at the crossing parameters are the given curves. + let [du, dv] = s.domain(); + for &y in &vy { + let f = ucurve(y); + // Find the v of this curve by its start point. + let start = f(0.0); + let off = |v: f64| s.eval(du[0], v).unwrap().dist(start); + let mut v = (0..=1000).map(|k| dv[0] + (dv[1] - dv[0]) * k as f64 / 1000.0).min_by(|a, b| off(*a).total_cmp(&off(*b))).unwrap(); + let (mut lo, mut hi) = ((v - 0.002).max(dv[0]), (v + 0.002).min(dv[1])); + for _ in 0..100 { + let (a, b) = (lo + (hi - lo) / 3.0, hi - (hi - lo) / 3.0); + if off(a) < off(b) { hi = b } else { lo = a } + } + v = (lo + hi) / 2.0; + for k in 0..=20 { + let q = s.eval(du[0] + (du[1] - du[0]) * k as f64 / 20.0, v).unwrap(); + assert!(dist_to(q, &f) < 1e-4, "y = {y}: {q:?} is {} off", dist_to(q, &f)); + } + } + for q in grid(&s, 10) { + assert!((q.z - z(q.x, q.y)).abs() < 0.05, "{q:?} strays from the network's surface"); + } + } + + /// A plane face along the x axis (normal +z, the face on the −y side): a frame for G1 and G2. + struct Plane; + impl CrossFace for Plane { + fn frame(&mut self, _: Vec3, _: Vec3) -> Option { + Some(FaceFrame { normal: Vec3::Z, outward: Vec3::Y, form: [0.0; 3] }) + } + } + + /// A cylinder of radius R about the x axis, the face below its edge at (y = R, z = 0): it + /// leaves the edge straight up (+z), bending toward the axis. + struct Cylinder(f64); + impl CrossFace for Cylinder { + fn frame(&mut self, at: Vec3, tangent: Vec3) -> Option { + let n = Vec3::new(0.0, at.y, at.z).normalized()?; + let _ = tangent; + // Outward normal n; curvature −1/R in the direction around the axis, none along it. + Some(FaceFrame { normal: n, outward: Vec3::Z, form: [0.0, 0.0, -1.0 / self.0] }) + } + } + + fn deriv_v(s: &NurbsSurface, u: f64, h: f64) -> (Vec3, Vec3) { + let f: Vec = (0..4).map(|k| s.eval(u, k as f64 * h).unwrap()).collect(); + ((f[0] * -3.0 + f[1] * 4.0 - f[2]) / (2.0 * h), (f[0] * 2.0 - f[1] * 5.0 + f[2] * 4.0 - f[3]) / (h * h)) + } + + /// G1: the surface leaves a plane edge tangent to the plane. G2 against a cylinder: the + /// surface's normal curvature across the edge is the cylinder's, 1/R; with G1 only it is + /// not. + #[test] + fn g1_and_g2_conditions_meet_the_face() { + let edge = |c: Box, g| NetworkCurve { + label: "edge".into(), + curve: Box::new(F(|t: f64| Vec3::new(10.0 * t, 0.0, 0.0))), + condition: Some((g, c)), + }; + let s = network_surface( + vec![edge(Box::new(Plane), Continuity::G1), line("far", Vec3::new(0.0, 10.0, 5.0), Vec3::new(10.0, 10.0, 5.0))], + vec![], + TOL, + ) + .unwrap(); + for k in 0..=10 { + let u = k as f64 / 10.0; + let (d1, _) = deriv_v(&s, u, 1e-5); + let angle = (d1.z / d1.len()).asin().to_degrees().abs(); + assert!(angle < 1e-3 && d1.y > 0.0, "u = {u}: leaves the plane at {angle}°"); + } + let r = 10.0; + let cyl = |g| { + let edge = NetworkCurve { + label: "rim".into(), + curve: Box::new(F(move |t: f64| Vec3::new(10.0 * t, r, 0.0))), + condition: Some((g, Box::new(Cylinder(r)) as Box)), + }; + network_surface(vec![edge, line("top", Vec3::new(0.0, 3.0, 12.0), Vec3::new(10.0, 3.0, 12.0))], vec![], TOL).unwrap() + }; + let curvature = |s: &NurbsSurface, u: f64| { + let (d1, d2) = deriv_v(s, u, 1e-4); + // Curvature of the v-curve: |d1 × d2| / |d1|³, toward the axis (−y). + (d1.cross(d2).len() / d1.len().powi(3), d2.y) + }; + let g2 = cyl(Continuity::G2); + for k in 0..=4 { + let (kappa, _) = curvature(&g2, k as f64 / 4.0); + assert!((kappa - 1.0 / r).abs() < 1e-3 / r, "G2: curvature {kappa} vs {}", 1.0 / r); + } + let g1 = cyl(Continuity::G1); + let (kappa, _) = curvature(&g1, 0.5); + assert!((kappa - 1.0 / r).abs() > 0.1 / r, "G1 alone does not match the curvature: {kappa}"); + } + + /// A Coons patch with a G1 edge: the side curves leave the plane tangent to it (arcs), so the + /// patch does too; a side curve that leaves it at an angle is refused, named. + #[test] + fn crossing_curves_must_agree_with_a_g1_face() { + let edge = || NetworkCurve { + label: "edge".into(), + curve: Box::new(F(|t: f64| Vec3::new(10.0 * t, 0.0, 0.0))), + condition: Some((Continuity::G1, Box::new(Plane) as Box)), + }; + // Arcs in x = const from (y, z) = (0, 0), tangent to +y, to (10, 5): centre (0, 12.5). + let side = |x: f64| { + move |t: f64| { + let a = -std::f64::consts::FRAC_PI_2 + t * 0.927_295_218_001_612_2; + Vec3::new(x, 12.5 * a.cos(), 12.5 + 12.5 * a.sin()) + } + }; + let far = line("far", Vec3::new(0.0, 10.0, 5.0), Vec3::new(10.0, 10.0, 5.0)); + let s = network_surface(vec![edge(), far], vec![curve("left", side(0.0)), curve("right", side(10.0))], TOL).unwrap(); + for k in 0..=10 { + let (d1, _) = deriv_v(&s, k as f64 / 10.0, 1e-5); + let angle = (d1.z / d1.len()).asin().to_degrees().abs(); + assert!(angle < 0.01, "leaves the plane at {angle}°"); + } + let far = line("far", Vec3::new(0.0, 10.0, 5.0), Vec3::new(10.0, 10.0, 5.0)); + let steep = line("steep", Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 10.0, 5.0)); + let err = network_surface(vec![edge(), far], vec![steep, curve("right", side(10.0))], TOL).unwrap_err(); + assert!(err.to_string().contains("steep") && err.to_string().contains("tangent"), "{err}"); + } + + /// Networks that can't be built are refused, naming the curves; none panics. + #[test] + fn bad_networks_are_refused() { + let l = |n: &str, a: [f64; 3], b: [f64; 3]| line(n, Vec3::from(a), Vec3::from(b)); + let miss = network_surface( + vec![l("a", [0.0, 0.0, 0.0], [10.0, 0.0, 0.0]), l("b", [0.0, 5.0, 0.0], [10.0, 5.0, 0.0])], + vec![l("c", [0.0, 0.0, 0.0], [0.0, 5.0, 0.0]), l("d", [10.0, 0.0, 1.0], [10.0, 5.0, 1.0])], + TOL, + ) + .unwrap_err(); + assert!(miss.to_string().contains("a and d") || miss.to_string().contains("b and d"), "{miss}"); + assert!(miss.to_string().contains("don't meet"), "{miss}"); + let one = network_surface(vec![l("a", [0.0, 0.0, 0.0], [1.0, 0.0, 0.0])], vec![], TOL); + assert!(one.is_err()); + assert!(network_surface(vec![], vec![], TOL).is_err()); + let circle = curve("ring", |t: f64| Vec3::new((6.3 * t).cos(), (6.3 * t).sin(), 0.0).lerp(Vec3::X, if t >= 1.0 { 1.0 } else { 0.0 })); + let closed = network_surface(vec![circle, l("b", [0.0, 0.0, 1.0], [1.0, 0.0, 1.0])], vec![], TOL).unwrap_err(); + assert!(closed.to_string().contains("closed"), "{closed}"); + let zero = network_surface(vec![l("a", [0.0; 3], [0.0; 3]), l("b", [0.0, 0.0, 1.0], [1.0, 0.0, 1.0])], vec![], TOL).unwrap_err(); + assert!(zero.to_string().contains("no length"), "{zero}"); + let nan = network_surface(vec![l("a", [f64::NAN, 0.0, 0.0], [1.0, 0.0, 0.0]), l("b", [0.0, 0.0, 1.0], [1.0, 0.0, 1.0])], vec![], TOL); + assert!(nan.is_err()); + // A parabola crossing a line twice. + let wave = curve("wave", |t: f64| Vec3::new(10.0 * t, 4.0 * (2.0 * t - 1.0).powi(2) - 1.0, 0.0)); + let twice = network_surface( + vec![l("a", [0.0, 0.0, 0.0], [10.0, 0.0, 0.0]), l("b", [0.0, 5.0, 0.0], [10.0, 5.0, 0.0])], + vec![l("c", [0.0, 0.0, 0.0], [0.0, 5.0, 0.0]), wave], + TOL, + ) + .unwrap_err(); + assert!(twice.to_string().contains("more than once"), "{twice}"); + // Conditions on both directions. + let cond = |n: &str, a: [f64; 3], b: [f64; 3]| { + let (a, b) = (Vec3::from(a), Vec3::from(b)); + NetworkCurve { + label: n.into(), + curve: Box::new(F(move |t: f64| a.lerp(b, t))), + condition: Some((Continuity::G1, Box::new(Plane) as Box)), + } + }; + let both = network_surface( + vec![cond("a", [0.0, 0.0, 0.0], [10.0, 0.0, 0.0]), l("b", [0.0, 5.0, 0.0], [10.0, 5.0, 0.0])], + vec![cond("c", [0.0, 0.0, 0.0], [0.0, 5.0, 0.0]), l("d", [10.0, 0.0, 0.0], [10.0, 5.0, 0.0])], + TOL, + ) + .unwrap_err(); + assert!(matches!(both, SplineError::Unsupported(_)), "{both}"); + let bad_tol = network_surface( + vec![l("a", [0.0; 3], [1.0, 0.0, 0.0]), l("b", [0.0, 1.0, 0.0], [1.0, 1.0, 0.0])], + vec![], + Tolerances { cross: f64::NAN, fit: 1.0 }, + ); + assert!(bad_tol.is_err()); + } +} diff --git a/crates/geom/src/spline_space.rs b/crates/geom/src/spline_space.rs new file mode 100644 index 00000000..6456ee8c --- /dev/null +++ b/crates/geom/src/spline_space.rs @@ -0,0 +1,584 @@ +//! Spline function spaces and the linear operators surface construction is built from: a clamped +//! B-spline space ([`Space`]), interpolation in it with derivative (Hermite) conditions +//! ([`Interpolant`]), carrying a spline into a larger space that contains it ([`embed`]: Bézier +//! degree elevation and Boehm knot insertion, exact), and a monotone piecewise cubic through +//! increasing data ([`Monotone`], for reparameterising curves). +//! +//! These are the pieces that make several splines *compatible* (one degree, one knot vector) so +//! that their control points can be combined linearly: lofting, Coons and Gordon surfaces, and +//! cross-boundary derivative conditions all reduce to them. +//! +//! References: Piegl & Tiller, *The NURBS Book* (2nd ed.) — §2.2–2.3 (basis functions and +//! their derivatives), §5.2 (knot insertion), §5.5 (degree elevation; the Bézier case), §9.2.1 +//! (global interpolation, knots by averaging the parameters) and §9.2.2 (interpolation with +//! end derivatives); Fritsch & Carlson 1980, "Monotone piecewise cubic interpolation" (SIAM J. +//! Numer. Anal. 17), with the weighted harmonic-mean slopes of Fritsch & Butland 1984. +//! +//! Input is hostile: constructors check sizes and finiteness and return [`SplineError`]; no +//! function here panics. + +use crate::Vec3; +use crate::bspline::{HPoint, MAX_SPLINE_DEGREE, Side, SplineError, SplineResult, check_knots, find_span, insert_knot, invalid, local_basis}; + +/// Most basis functions a space used for interpolation may have (its systems are solved densely). +pub const MAX_SPACE_DIMENSION: usize = 2048; + +/// Values a spline can carry: numbers and points (vectors). +pub trait Linear: Copy { + fn zero() -> Self; + fn plus(self, o: Self) -> Self; + fn times(self, s: f64) -> Self; + /// The value as a homogeneous point of weight one (knot insertion works on those). + fn to_h(self) -> HPoint; + fn from_h(h: HPoint) -> Self; + fn is_finite(self) -> bool; +} + +impl Linear for f64 { + fn zero() -> Self { + 0.0 + } + fn plus(self, o: Self) -> Self { + self + o + } + fn times(self, s: f64) -> Self { + self * s + } + fn to_h(self) -> HPoint { + HPoint { x: self, y: 0.0, z: 0.0, w: 1.0 } + } + fn from_h(h: HPoint) -> Self { + h.x + } + fn is_finite(self) -> bool { + f64::is_finite(self) + } +} + +impl Linear for Vec3 { + fn zero() -> Self { + Vec3::ZERO + } + fn plus(self, o: Self) -> Self { + self + o + } + fn times(self, s: f64) -> Self { + self * s + } + fn to_h(self) -> HPoint { + HPoint::weighted(self, 1.0) + } + fn from_h(h: HPoint) -> Self { + h.xyz() + } + fn is_finite(self) -> bool { + Vec3::is_finite(self) + } +} + +/// A clamped B-spline space: a degree and a knot vector whose first and last `degree + 1` knots +/// are the domain's ends, interior knots repeated at most `degree` times (the space's splines +/// are continuous). +#[derive(Clone, Debug, PartialEq)] +pub struct Space { + degree: usize, + knots: Vec, +} + +impl Space { + pub fn new(degree: usize, knots: Vec) -> SplineResult { + if degree == 0 || degree > MAX_SPLINE_DEGREE { + return Err(SplineError::Unsupported(format!("degree {degree} (1 to {MAX_SPLINE_DEGREE})"))); + } + check_knots(&knots)?; + let n = knots.len().checked_sub(degree + 1).filter(|n| *n > degree).ok_or_else(|| SplineError::Invalid("too few knots".into()))?; + if n > MAX_SPACE_DIMENSION { + return Err(SplineError::TooLarge(format!("{n} basis functions (at most {MAX_SPACE_DIMENSION})"))); + } + let (a, b) = (knots[0], knots[knots.len() - 1]); + if !(a < b) || knots[..=degree].iter().any(|k| *k != a) || knots[n..].iter().any(|k| *k != b) { + return invalid("the knot vector must be clamped on a non-empty domain"); + } + let inner = &knots[degree + 1..n]; + let mut run = 0; + for w in inner.windows(2) { + run = if w[0] == w[1] { run + 1 } else { 0 }; + if run + 1 > degree { + return invalid("an interior knot is repeated more often than the degree"); + } + } + if inner.first().is_some_and(|k| *k == a) || inner.last().is_some_and(|k| *k == b) { + return invalid("an interior knot lies on a domain end"); + } + Ok(Space { degree, knots }) + } + + /// Polynomials of `degree` on `[a, b]` (one Bézier piece). + pub fn polynomial(degree: usize, a: f64, b: f64) -> SplineResult { + let mut knots = vec![a; degree + 1]; + knots.extend(std::iter::repeat_n(b, degree + 1)); + Space::new(degree, knots) + } + + /// Splines of `degree` with a simple knot at every interior break: `breaks` increase + /// strictly, the first and last are the domain's ends. + pub fn with_breaks(degree: usize, breaks: &[f64]) -> SplineResult { + let (Some(&a), Some(&b)) = (breaks.first(), breaks.last()) else { return invalid("no breaks") }; + if breaks.len() > MAX_SPACE_DIMENSION || breaks.windows(2).any(|w| !(w[0] < w[1])) { + return invalid("the breaks must increase"); + } + let mut knots = vec![a; degree + 1]; + knots.extend(breaks.iter().skip(1).take(breaks.len().saturating_sub(2)).copied()); + knots.extend(std::iter::repeat_n(b, degree + 1)); + Space::new(degree, knots) + } + + pub fn degree(&self) -> usize { + self.degree + } + + pub fn knots(&self) -> &[f64] { + &self.knots + } + + /// The number of basis functions (control points of a spline in the space). + pub fn count(&self) -> usize { + self.knots.len() - self.degree - 1 + } + + /// The domain `[a, b]`. + pub fn domain(&self) -> [f64; 2] { + [self.knots[0], self.knots[self.knots.len() - 1]] + } + + /// The `der`-th derivatives at `t` (clamped to the domain) of the basis functions that do not + /// vanish there: the index of the first one and their values. + pub fn basis(&self, t: f64, der: usize) -> SplineResult<(usize, Vec)> { + let [a, b] = self.domain(); + if !t.is_finite() { + return invalid("a parameter is not a finite number"); + } + let t = t.clamp(a, b); + let (p, n) = (self.degree, self.count()); + let span = find_span(&self.knots, p, n, t).ok_or_else(|| SplineError::Invalid(format!("{t} is outside the domain")))?; + let side = if t >= b { Side::Left } else { Side::Right }; + let first = span - p; + let values = + (first..=span).map(|i| self.knots.get(i..i + p + 2).map_or(0.0, |k| local_basis(k, t, side, der).get(der).copied().unwrap_or(0.0))); + Ok((first, values.collect())) + } + + /// The `der`-th derivative at `t` of the spline with control values `ctrl`. + pub fn eval(&self, ctrl: &[T], t: f64, der: usize) -> SplineResult { + if ctrl.len() != self.count() { + return invalid(format!("{} control values for {} basis functions", ctrl.len(), self.count())); + } + let (first, values) = self.basis(t, der)?; + Ok(values.iter().zip(ctrl.iter().skip(first)).fold(T::zero(), |acc, (n, c)| acc.plus(c.times(*n)))) + } + + /// The splines of this space with `breaks` added as simple knots (those already knots are + /// left as they are). The result contains this space. + pub fn refined(&self, breaks: &[f64]) -> SplineResult { + let [a, b] = self.domain(); + let mut knots = self.knots.clone(); + for &x in breaks { + if x > a && x < b && !knots.contains(&x) { + knots.push(x); + } + } + knots.sort_by(f64::total_cmp); + Space::new(self.degree, knots) + } +} + +/// A condition on a spline: the value (`der` = 0) or a derivative at `t`. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Cond { + pub t: f64, + pub der: usize, +} + +/// Interpolation in a space: the spline whose values and derivatives at the conditions are given +/// (as many conditions as the space has basis functions). The system is factored once; each +/// [`Interpolant::fit`] solves it for new data, so one interpolant serves many curves and the +/// rows of a surface alike (the operator is linear). +#[derive(Clone, Debug)] +pub struct Interpolant { + space: Space, + conds: Vec, + lu: Lu, +} + +impl Interpolant { + pub fn new(space: Space, conds: Vec) -> SplineResult { + let n = space.count(); + if conds.len() != n { + return invalid(format!("{} conditions for {n} basis functions", conds.len())); + } + let mut a = vec![0.0; n * n]; + for (r, c) in conds.iter().enumerate() { + let (first, values) = space.basis(c.t, c.der)?; + for (k, v) in values.iter().enumerate() { + if let Some(x) = a.get_mut(r * n + first + k) { + *x = *v; + } + } + } + let lu = Lu::new(n, a).map_err(|_| SplineError::Invalid("the interpolation conditions do not determine a spline".into()))?; + Ok(Interpolant { space, conds, lu }) + } + + /// Interpolation by a single polynomial piece of the lowest degree that meets `conds` (degree + /// one less than their number, at least one: a lone value is held constant), on `[a, b]`; + /// with more than four conditions, a cubic spline whose interior knots average the + /// conditions' parameters (taken as a sequence with one entry per condition, NURBS Book eq. + /// 9.8 applied to that confluent sequence, which keeps every basis function's support over + /// the condition it is matched with). + pub fn lowest(conds: Vec, a: f64, b: f64) -> SplineResult { + let mut conds = conds; + if conds.iter().any(|c| !c.t.is_finite() || c.t < a || c.t > b) { + return invalid("an interpolation parameter is outside the domain"); + } + if conds.len() == 1 + && let Some(c) = conds.first().copied() + { + // A constant: the value and a zero slope (the caller supplies zero data for it). + conds.push(Cond { t: c.t, der: 1 }); + } + let k = conds.len(); + if k <= 4 { + return Interpolant::new(Space::polynomial(k.saturating_sub(1).max(1), a, b)?, conds); + } + let mut tau: Vec = conds.iter().map(|c| c.t).collect(); + tau.sort_by(f64::total_cmp); + let mut breaks = vec![a]; + for j in 1..k - 3 { + breaks.push((tau[j] + tau[j + 1] + tau[j + 2]) / 3.0); + } + breaks.push(b); + if breaks.windows(2).any(|w| !(w[0] < w[1])) { + return invalid("the interpolation parameters repeat too often"); + } + Interpolant::new(Space::with_breaks(3, &breaks)?, conds) + } + + pub fn space(&self) -> &Space { + &self.space + } + + pub fn conds(&self) -> &[Cond] { + &self.conds + } + + /// The control values of the spline meeting the conditions with `data` (one per condition). + pub fn fit(&self, data: &[T]) -> SplineResult> { + if data.len() != self.conds.len() { + return invalid(format!("{} values for {} conditions", data.len(), self.conds.len())); + } + let out = self.lu.solve(data); + if out.iter().any(|x| !x.is_finite()) { + return invalid("interpolation did not give finite control values"); + } + Ok(out) + } +} + +/// An LU factorisation with partial pivoting of a dense square matrix (row-major). +#[derive(Clone, Debug)] +struct Lu { + n: usize, + a: Vec, + piv: Vec, +} + +impl Lu { + fn new(n: usize, mut a: Vec) -> SplineResult { + if a.len() != n * n { + return invalid("matrix size"); + } + let scale = a.iter().fold(0.0f64, |m, x| m.max(x.abs())); + if !(scale > 0.0 && scale.is_finite()) { + return invalid("singular system"); + } + let mut piv = vec![0; n]; + for k in 0..n { + let (mut best, mut big) = (k, 0.0); + for r in k..n { + let x = a[r * n + k].abs(); + if x > big { + (best, big) = (r, x); + } + } + if !(big > scale * 1e-13) { + return invalid("singular system"); + } + piv[k] = best; + if best != k { + for c in 0..n { + a.swap(k * n + c, best * n + c); + } + } + let d = a[k * n + k]; + for r in k + 1..n { + let f = a[r * n + k] / d; + a[r * n + k] = f; + if f != 0.0 { + for c in k + 1..n { + a[r * n + c] -= f * a[k * n + c]; + } + } + } + } + Ok(Lu { n, a, piv }) + } + + fn solve(&self, b: &[T]) -> Vec { + let n = self.n; + let mut x = b.to_vec(); + if x.len() != n { + return x; + } + for (k, &p) in self.piv.iter().enumerate() { + x.swap(k, p); + } + for r in 0..n { + let row = &self.a[r * n..r * n + r]; + let acc = row.iter().zip(&x[..r]).fold(x[r], |acc, (l, xc)| acc.plus(xc.times(-l))); + x[r] = acc; + } + for r in (0..n).rev() { + let row = &self.a[r * n + r + 1..(r + 1) * n]; + let acc = row.iter().zip(&x[r + 1..]).fold(x[r], |acc, (u, xc)| acc.plus(xc.times(-u))); + x[r] = acc.times(1.0 / self.a[r * n + r]); + } + x + } +} + +/// The control values in `to` of the spline with control values `ctrl` in `from`, which `to` +/// contains: the same domain, a degree at least as high, and every knot of `from` among those of +/// `to` (with a polynomial `from` when the degree rises). Exact: Bézier degree elevation (NURBS +/// Book §5.5) to the target degree, then Boehm knot insertion (§5.2) of the knots `to` adds. +pub fn embed(from: &Space, ctrl: &[T], to: &Space) -> SplineResult> { + if ctrl.len() != from.count() { + return invalid("control values do not match the space"); + } + let ([a0, b0], [a1, b1]) = (from.domain(), to.domain()); + if a0 != a1 || b0 != b1 { + return invalid("the spaces' domains differ"); + } + let mut p = from.degree(); + let mut knots = from.knots().to_vec(); + let mut pts: Vec = ctrl.to_vec(); + if p > to.degree() { + return invalid("the target space's degree is lower"); + } + if p < to.degree() { + if from.count() != p + 1 { + return Err(SplineError::Unsupported("raising the degree of a spline with interior knots".into())); + } + while p < to.degree() { + pts = elevate_bezier(&pts); + p += 1; + } + knots = vec![a0; p + 1]; + knots.extend(std::iter::repeat_n(b0, p + 1)); + } + // The knots `to` adds: a merge of the two sorted vectors. + let target = to.knots(); + let mut missing = Vec::new(); + let mut j = 0; + for &x in target { + match knots.get(j) { + Some(&k) if k == x => j += 1, + Some(&k) if k < x => return invalid("a knot of the spline is not in the target space"), + _ => missing.push(x), + } + } + if j != knots.len() { + return invalid("a knot of the spline is not in the target space"); + } + let mut h: Vec = pts.iter().map(|x| x.to_h()).collect(); + for x in missing { + (knots, h) = insert_knot(p, &knots, &h, x)?; + } + if knots.as_slice() != target { + return invalid("knot insertion did not reach the target space"); + } + Ok(h.into_iter().map(T::from_h).collect()) +} + +/// A Bézier piece's control values one degree higher, the same polynomial (NURBS Book eq. 5.36): +/// `Qᵢ = (i/(p+1))·Pᵢ₋₁ + (1 − i/(p+1))·Pᵢ`. +fn elevate_bezier(pts: &[T]) -> Vec { + let p = pts.len().saturating_sub(1); + let at = |i: usize| pts.get(i).copied().unwrap_or_else(T::zero); + (0..=p + 1) + .map(|i| { + let a = i as f64 / (p + 1) as f64; + let prev = if i == 0 { T::zero() } else { at(i - 1).times(a) }; + prev.plus(at(i).times(1.0 - a)) + }) + .collect() +} + +/// A monotone piecewise cubic Hermite function through `(xₖ, yₖ)`, both strictly increasing +/// (Fritsch & Carlson): slopes are the weighted harmonic means of the neighbouring secants +/// (Fritsch & Butland), the secant at the ends. Its slope is positive everywhere, and through two +/// points it is the straight line. +#[derive(Clone, Debug, PartialEq)] +pub struct Monotone { + x: Vec, + y: Vec, + m: Vec, +} + +impl Monotone { + pub fn new(x: Vec, y: Vec) -> SplineResult { + let n = x.len(); + if n < 2 || y.len() != n { + return invalid("a monotone cubic needs two points or more"); + } + let increasing = |v: &[f64]| v.iter().all(|t| t.is_finite()) && v.windows(2).all(|w| w[0] < w[1]); + if !increasing(&x) || !increasing(&y) { + return invalid("the points of a monotone cubic must increase"); + } + let h: Vec = x.windows(2).map(|w| w[1] - w[0]).collect(); + let d: Vec = y.windows(2).zip(&h).map(|(w, h)| (w[1] - w[0]) / h).collect(); + let mut m = Vec::with_capacity(n); + m.push(d[0]); + for k in 1..n - 1 { + let (w1, w2) = (2.0 * h[k] + h[k - 1], h[k] + 2.0 * h[k - 1]); + m.push((w1 + w2) / (w1 / d[k - 1] + w2 / d[k])); + } + m.push(d[n - 2]); + Ok(Monotone { x, y, m }) + } + + /// The value and the first two derivatives at `t` (clamped to the points' range). + pub fn eval(&self, t: f64) -> [f64; 3] { + let n = self.x.len(); + let t = t.clamp(self.x[0], self.x[n - 1]); + let k = self.x.partition_point(|x| *x <= t).saturating_sub(1).min(n - 2); + let (x0, x1, y0, y1, m0, m1) = (self.x[k], self.x[k + 1], self.y[k], self.y[k + 1], self.m[k], self.m[k + 1]); + let h = x1 - x0; + let s = (t - x0) / h; + let (s2, s3) = (s * s, s * s * s); + let value = (2.0 * s3 - 3.0 * s2 + 1.0) * y0 + (s3 - 2.0 * s2 + s) * h * m0 + (-2.0 * s3 + 3.0 * s2) * y1 + (s3 - s2) * h * m1; + let d1 = ((6.0 * s2 - 6.0 * s) * y0 + (3.0 * s2 - 4.0 * s + 1.0) * h * m0 + (-6.0 * s2 + 6.0 * s) * y1 + (3.0 * s2 - 2.0 * s) * h * m1) / h; + let d2 = ((12.0 * s - 6.0) * y0 + (6.0 * s - 4.0) * h * m0 + (-12.0 * s + 6.0) * y1 + (6.0 * s - 2.0) * h * m1) / (h * h); + [value, d1, d2] + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn close(a: f64, b: f64) -> bool { + (a - b).abs() < 1e-10 + } + + /// Hermite interpolation reproduces a cubic from its values and end derivatives, and the + /// lowest-degree interpolant through two values is the straight line. + #[test] + fn interpolation_meets_values_and_derivatives() { + let f = |t: f64| 1.0 + 2.0 * t - 3.0 * t * t + 0.5 * t * t * t; + let df = |t: f64| 2.0 - 6.0 * t + 1.5 * t * t; + let d2f = |t: f64| -6.0 + 3.0 * t; + let conds = vec![Cond { t: 0.0, der: 0 }, Cond { t: 0.0, der: 1 }, Cond { t: 0.0, der: 2 }, Cond { t: 2.0, der: 0 }]; + let it = Interpolant::lowest(conds, 0.0, 2.0).unwrap(); + let c = it.fit(&[f(0.0), df(0.0), d2f(0.0), f(2.0)]).unwrap(); + for k in 0..=10 { + let t = k as f64 * 0.2; + assert!(close(it.space().eval(&c, t, 0).unwrap(), f(t)), "{t}"); + } + let line = Interpolant::lowest(vec![Cond { t: 0.0, der: 0 }, Cond { t: 1.0, der: 0 }], 0.0, 1.0).unwrap(); + let c = line.fit(&[2.0, 5.0]).unwrap(); + assert!(close(line.space().eval(&c, 0.25, 0).unwrap(), 2.75)); + let constant = Interpolant::lowest(vec![Cond { t: 0.3, der: 0 }], 0.0, 1.0).unwrap(); + let c = constant.fit(&[4.0, 0.0]).unwrap(); + assert!(close(constant.space().eval(&c, 0.9, 0).unwrap(), 4.0)); + } + + /// Many conditions: a cubic spline with averaged knots through samples of a smooth function + /// with derivative conditions at both ends (values exact at the samples). + #[test] + fn spline_interpolation_with_end_derivatives() { + let f = |t: f64| (2.0 * t).sin(); + let df = |t: f64| 2.0 * (2.0 * t).cos(); + let ts: Vec = (0..=8).map(|k| k as f64 / 8.0).collect(); + let mut conds: Vec = ts.iter().map(|&t| Cond { t, der: 0 }).collect(); + conds.push(Cond { t: 0.0, der: 1 }); + conds.push(Cond { t: 1.0, der: 1 }); + let it = Interpolant::lowest(conds, 0.0, 1.0).unwrap(); + let mut data: Vec = ts.iter().map(|&t| f(t)).collect(); + data.push(df(0.0)); + data.push(df(1.0)); + let c = it.fit(&data).unwrap(); + for &t in &ts { + assert!(close(it.space().eval(&c, t, 0).unwrap(), f(t))); + } + assert!(close(it.space().eval(&c, 0.0, 1).unwrap(), df(0.0))); + assert!(close(it.space().eval(&c, 1.0, 1).unwrap(), df(1.0))); + assert!((it.space().eval(&c, 0.53, 0).unwrap() - f(0.53)).abs() < 1e-4); + } + + /// Embedding is exact: a quadratic polynomial raised to a cubic with knots inserted, and a + /// cubic spline refined, evaluate the same everywhere; a space that does not contain the + /// spline is refused. + #[test] + fn embedding_keeps_the_spline() { + let from = Space::polynomial(2, 0.0, 1.0).unwrap(); + let ctrl = [Vec3::new(0.0, 0.0, 0.0), Vec3::new(1.0, 2.0, 0.0), Vec3::new(3.0, -1.0, 1.0)]; + let to = Space::with_breaks(3, &[0.0, 0.2, 0.5, 0.7, 1.0]).unwrap(); + let out = embed(&from, &ctrl, &to).unwrap(); + let spline = Space::with_breaks(3, &[0.0, 0.5, 1.0]).unwrap(); + let sctrl: Vec = (0..spline.count()).map(|i| (i as f64 * 1.7).sin()).collect(); + let finer = spline.refined(&[0.25, 0.5, 0.8]).unwrap(); + let sout = embed(&spline, &sctrl, &finer).unwrap(); + for k in 0..=20 { + let t = k as f64 / 20.0; + assert!(from.eval(&ctrl, t, 0).unwrap().dist(to.eval(&out, t, 0).unwrap()) < 1e-12); + assert!(close(spline.eval(&sctrl, t, 0).unwrap(), finer.eval(&sout, t, 0).unwrap())); + } + assert!(embed(&finer, &sout, &spline).is_err(), "a coarser space"); + assert!(embed(&spline, &sctrl, &Space::with_breaks(3, &[0.0, 0.4, 1.0]).unwrap()).is_err(), "a knot missing"); + assert!(embed(&spline, &sctrl, &Space::with_breaks(3, &[0.0, 0.5, 2.0]).unwrap()).is_err(), "another domain"); + } + + #[test] + fn monotone_cubic_is_monotone_and_straight_through_two_points() { + let m = Monotone::new(vec![0.0, 0.3, 0.5, 1.0], vec![0.0, 0.05, 0.9, 1.0]).unwrap(); + let mut last = -1.0; + for k in 0..=200 { + let [y, d, _] = m.eval(k as f64 / 200.0); + assert!(y >= last && d >= 0.0, "{k}"); + last = y; + } + assert!(close(m.eval(0.3)[0], 0.05) && close(m.eval(0.5)[0], 0.9)); + let line = Monotone::new(vec![0.0, 1.0], vec![2.0, 4.0]).unwrap(); + assert!(close(line.eval(0.25)[0], 2.5) && close(line.eval(0.7)[1], 2.0) && close(line.eval(0.7)[2], 0.0)); + assert!(Monotone::new(vec![0.0, 0.0], vec![0.0, 1.0]).is_err()); + assert!(Monotone::new(vec![0.0, 1.0], vec![1.0, 0.0]).is_err()); + assert!(Monotone::new(vec![0.0, f64::NAN], vec![0.0, 1.0]).is_err()); + } + + #[test] + fn hostile_spaces_are_refused() { + assert!(Space::new(3, vec![0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]).is_ok()); + assert!(Space::new(3, vec![0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]).is_err(), "too few knots"); + assert!(Space::new(3, vec![0.0, 0.0, 0.0, 0.1, 1.0, 1.0, 1.0, 1.0]).is_err(), "not clamped"); + assert!(Space::new(1, vec![0.0, 0.0, 0.5, 0.5, 1.0, 1.0]).is_err(), "a discontinuous space"); + assert!(Space::new(0, vec![0.0, 1.0]).is_err()); + assert!(Space::with_breaks(3, &[0.0, 0.5, 0.5, 1.0]).is_err()); + assert!(Space::with_breaks(3, &[0.0, f64::NAN, 1.0]).is_err()); + let s = Space::with_breaks(3, &[0.0, 1.0]).unwrap(); + assert!(s.eval(&[1.0, 2.0], 0.5, 0).is_err(), "wrong control count"); + assert!(s.basis(f64::NAN, 0).is_err()); + // Repeated conditions do not determine a spline. + assert!(Interpolant::new(s.clone(), vec![Cond { t: 0.5, der: 0 }; 4]).is_err()); + assert!(Interpolant::lowest(vec![Cond { t: 2.0, der: 0 }], 0.0, 1.0).is_err()); + } +} diff --git a/crates/kernel/src/boundary.rs b/crates/kernel/src/boundary.rs index c0825f2e..0a63df1b 100644 --- a/crates/kernel/src/boundary.rs +++ b/crates/kernel/src/boundary.rs @@ -6,13 +6,16 @@ //! 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 +//! A patch through *three or four* boundaries is refused here 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. +//! (issue #763). A loop of one sketch's curves is the case where they agree: it lies in the +//! sketch's plane, and its flat Coons patch ([`crate::boundary_network`]) holds every order +//! exactly; the network there refuses G1/G2 on two sides that meet at a corner for the same +//! reason. Everything this module refuses is an error, never a quietly worse surface. use solidcraft_geom::Vec3; diff --git a/crates/kernel/src/boundary_network.rs b/crates/kernel/src/boundary_network.rs new file mode 100644 index 00000000..ea109e4a --- /dev/null +++ b/crates/kernel/src/boundary_network.rs @@ -0,0 +1,744 @@ +//! Boundary surfaces (Surface › Boundary): one B-spline face through curves in one or two +//! directions (`solidcraft_geom::network`: a loft between curves, a Coons patch inside four, a +//! Gordon surface through a network), made from sketch curves, 3D sketch curves or a body's +//! edges. A boundary curve that is a chain of open edges of a surface may carry a continuity +//! condition with the face beside it: G1 (tangent) or G2 (curvature). +//! +//! The reusable pieces: [`CurveSource`] turns each kind of curve into a curve the network can +//! measure, [`EdgeFaces`] gives the frame and second fundamental form of the faces along a chain +//! of edges (what a G1 or G2 condition needs), and [`boundary_surface`] builds the face. +//! +//! The two-curve bridge with G0–G3 against a sketch's plane is [`crate::boundary`]; this is the +//! general case: three or four sides, interior curves, 3D curves and edges. + +use std::collections::HashSet; + +use mt::{BoundedCurve, Cut, InnerSpace, ParametricCurve, ParametricSurface, SearchNearestParameter}; +use solidcraft_geom::bspline::{CurveN, NurbsSurface, SplineError}; +use solidcraft_geom::network::{self, CrossFace, FaceFrame, NetCurve, NetworkCurve, PointCurve, Tolerances}; +use solidcraft_geom::{Plane, Seg2, Vec3}; +use truck_modeling as mt; + +use crate::body::{Body, from_p3, p3}; +use crate::bridge::Continuity; +use crate::plane_fit::EdgeSampler; +use crate::record::Rebuild; +use crate::{KernelError, Result, guard}; + +/// How far apart curves of the two directions may pass and still cross (mm). +pub const CROSS_TOLERANCE: f64 = 0.01; +/// How far the surface may stray from its curves (mm). +pub const FIT_TOLERANCE: f64 = 1e-4; +/// Most points of a 3D sketch curve a boundary surface follows (longer ones are resampled). +const MAX_WIRE_POINTS: usize = 400; +/// Most points picking the edges of one curve. +pub const MAX_EDGE_PICKS: usize = 10_000; +/// Most sketch curves one curve of a boundary surface chains. +pub const MAX_CHAIN_CURVES: usize = 1000; + +/// The network's continuity for a curve's: G3 is refused (a network holds up to G2). +fn net_order(c: &BoundaryCurve) -> Result { + match c.continuity { + Continuity::G0 => Ok(network::Continuity::G0), + Continuity::G1 => Ok(network::Continuity::G1), + Continuity::G2 => Ok(network::Continuity::G2), + Continuity::G3 => Err(KernelError::Invalid(format!("not supported yet: G3 continuity on {} (a curve network holds up to G2)", c.label))), + } +} + +/// Where a curve of a boundary surface comes from. +pub enum CurveSource<'a> { + /// Sketch segments on a plane, each starting where the one before it ends. + Sketch { plane: Plane, segs: Vec }, + /// A 3D sketch curve's points: followed smoothly through them (a cubic spline). + Points(Vec), + /// A chain of a body's edges, each picked by a point on it. + Edges { body: &'a Body, at: Vec }, + /// A side collapsed to a point: the apex of a three-sided patch, where the sides beside it + /// meet. + Point(Vec3), +} + +/// A curve of a boundary surface: a name for messages, where it comes from, and (for edges of a +/// surface) its continuity with the face beside it. +pub struct BoundaryCurve<'a> { + pub label: String, + pub source: CurveSource<'a>, + pub continuity: Continuity, +} + +/// One B-spline face through curves `u` (two or more; or one with a continuity condition, crossed +/// by `v`) and `v` (none: a loft through `u`; two each way and no more: a Coons patch; more: a +/// Gordon surface through the network). A surface body; its history has no operands (the face +/// and its edges are new). Refused, naming the curves: curves that miss each other (by more than +/// [`CROSS_TOLERANCE`]) or cross twice, curves out of order, closed curves, a continuity +/// condition on anything but an open edge of a surface, or on curves of both directions. +pub fn boundary_surface(u: &[BoundaryCurve], v: &[BoundaryCurve]) -> Result { + for c in u.iter().chain(v) { + if let CurveSource::Edges { body, .. } = &c.source { + body.require_brep("boundary surface")?; + } + } + guard("boundary surface", || { + let mut corners: Vec = Vec::new(); + let mut to_net = |c: &BoundaryCurve| -> Result> { + let (curve, face): (Box, Option) = match &c.source { + CurveSource::Sketch { plane, segs } => (Box::new(SketchChain::new(&c.label, plane, segs)?), None), + CurveSource::Points(pts) => (Box::new(Polyline::new(&c.label, pts)?), None), + CurveSource::Point(p) => { + if !p.is_finite() { + return Err(KernelError::Invalid(format!("{} is not a finite point", c.label))); + } + (Box::new(PointCurve(*p)), None) + } + CurveSource::Edges { body, at } => { + if at.is_empty() || at.len() > MAX_EDGE_PICKS || at.iter().any(|p| !p.is_finite()) { + return Err(KernelError::Invalid(format!("{}: pick 1 to {MAX_EDGE_PICKS} points on its edges", c.label))); + } + let (chain, closed) = crate::surfaces::picked_chain(body, at)?; + if closed { + return Err(KernelError::Invalid(format!("not supported yet: {} is a closed chain of edges", c.label))); + } + // The chain's inner vertices: the face's sides are split there, so that a + // join shares the body's edges one for one. + corners.extend(chain.iter().skip(1).map(|e| from_p3(e.front().point()))); + let face = match c.continuity { + Continuity::G0 => None, + _ => Some(EdgeFaces::new(&c.label, body, &chain)?), + }; + (Box::new(EdgeChain::new(&c.label, &chain)?), face) + } + }; + if c.continuity != Continuity::G0 && face.is_none() { + return Err(KernelError::Invalid(format!("{}: G1 and G2 apply only to a body's edges", c.label))); + } + let order = net_order(c)?; + Ok(NetworkCurve { label: c.label.clone(), curve, condition: face.map(|f| (order, Box::new(f) as Box)) }) + }; + let nu: Vec = u.iter().map(&mut to_net).collect::>()?; + let nv: Vec = v.iter().map(&mut to_net).collect::>()?; + let surface = network::network_surface(nu, nv, Tolerances { cross: CROSS_TOLERANCE, fit: FIT_TOLERANCE }).map_err(|e| match e { + SplineError::Unsupported(m) => KernelError::Invalid(format!("not supported yet: {m}")), + SplineError::TooLarge(m) => KernelError::TooLarge(m), + SplineError::Invalid(m) => KernelError::Invalid(m), + })?; + let face = face_of(&surface, &corners)?; + let mut rec = Rebuild::default(); + rec.made(&face, Vec::new()); + crate::surfaces::body_of(vec![face], &[], &rec) + }) +} + +/// `made` (a boundary surface) sewn into `target` along the edges they share: one body, a solid +/// when it closes the gap. Refused when they share no edge. +pub fn boundary_join(target: &Body, made: &Body) -> Result { + target.require_brep("join")?; + let edges = |b: &Body| b.solid.edge_iter().map(|e| e.id()).collect::>().len(); + let joined = crate::surfaces::stitch(&[target, made], CROSS_TOLERANCE)?; + if edges(&joined) == edges(target) + edges(made) { + return Err(KernelError::Invalid("the surface shares no edge with the body it joins".into())); + } + Ok(joined) +} + +/// Sketch segments on a plane, run end to end by arc length. +struct SketchChain { + plane: Plane, + segs: Vec, + /// Cumulative length at the start of each segment, and the total. + starts: Vec, + total: f64, +} + +impl SketchChain { + fn new(label: &str, plane: &Plane, segs: &[Seg2]) -> Result { + if segs.is_empty() || segs.len() > MAX_CHAIN_CURVES { + return Err(KernelError::Invalid(format!("{label} has no curves, or too many"))); + } + if segs.windows(2).any(|w| w[0].end().dist(w[1].start()) > 1e-6) { + return Err(KernelError::Invalid(format!("{label}: its curves don't meet end to end"))); + } + let mut starts = Vec::with_capacity(segs.len()); + let mut total = 0.0; + for s in segs { + starts.push(total); + total += s.length(); + } + if !(total > 0.0 && total.is_finite()) { + return Err(KernelError::Invalid(format!("{label} has no length"))); + } + Ok(SketchChain { plane: *plane, segs: segs.to_vec(), starts, total }) + } +} + +impl NetCurve for SketchChain { + fn point(&mut self, t: f64) -> Option { + let at = t.clamp(0.0, 1.0) * self.total; + let k = self.starts.partition_point(|s| *s <= at).saturating_sub(1); + let (seg, start) = (self.segs.get(k)?, self.starts.get(k)?); + let len = seg.length(); + let local = if len > 0.0 { ((at - start) / len).clamp(0.0, 1.0) } else { 0.0 }; + let p = self.plane.to_world(seg.point_at(local)); + p.is_finite().then_some(p) + } +} + +/// A 3D polyline followed smoothly: a cubic through its points by chord length (fewer points, a +/// lower degree; a long one resampled first). +struct Polyline { + curve: CurveN<3>, + end: f64, +} + +impl Polyline { + fn new(label: &str, pts: &[Vec3]) -> Result { + let mut clean: Vec = Vec::with_capacity(pts.len()); + for p in pts { + if !p.is_finite() { + return Err(KernelError::Invalid(format!("{label} has a point that is not a finite number"))); + } + if clean.last().is_none_or(|q: &Vec3| q.dist(*p) > 1e-9) { + clean.push(*p); + } + } + let pts = resample(&clean, MAX_WIRE_POINTS); + if pts.len() < 2 { + return Err(KernelError::Invalid(format!("{label} has no length"))); + } + let mut params = vec![0.0]; + for w in pts.windows(2) { + params.push(params.last().copied().unwrap_or(0.0) + w[0].dist(w[1])); + } + let end = params.last().copied().unwrap_or(0.0); + let data: Vec<[f64; 3]> = pts.iter().map(|p| [p.x, p.y, p.z]).collect(); + let curve = CurveN::interpolate(¶ms, &data, (pts.len() - 1).min(3)) + .map_err(|e| KernelError::Invalid(format!("{label} can't be followed smoothly: {e}")))?; + Ok(Polyline { curve, end }) + } +} + +impl NetCurve for Polyline { + fn point(&mut self, t: f64) -> Option { + let ([x, y, z], _) = self.curve.eval(t.clamp(0.0, 1.0) * self.end)?; + Some(Vec3::new(x, y, z)) + } +} + +/// At most `max` points evenly along a polyline by arc length, its ends kept. +fn resample(pts: &[Vec3], max: usize) -> Vec { + if pts.len() <= max || max < 2 { + return pts.to_vec(); + } + let mut cum = vec![0.0]; + for w in pts.windows(2) { + cum.push(cum.last().copied().unwrap_or(0.0) + w[0].dist(w[1])); + } + let total = cum.last().copied().unwrap_or(0.0); + (0..max) + .filter_map(|k| { + let at = total * k as f64 / (max - 1) as f64; + let i = cum.partition_point(|c| *c <= at).saturating_sub(1).min(pts.len().saturating_sub(2)); + let (c0, c1) = (cum.get(i)?, cum.get(i + 1)?); + let w = if c1 > c0 { ((at - c0) / (c1 - c0)).clamp(0.0, 1.0) } else { 0.0 }; + Some(pts.get(i)?.lerp(*pts.get(i + 1)?, w)) + }) + .collect() +} + +/// A chain of edges (each turned to run along it), by their lengths. +struct EdgeChain { + edges: Vec, + starts: Vec, + total: f64, +} + +impl EdgeChain { + fn new(label: &str, chain: &[mt::Edge]) -> Result { + let mut edges = Vec::with_capacity(chain.len()); + let mut starts = Vec::with_capacity(chain.len()); + let mut total = 0.0; + for e in chain { + let mut s = EdgeSampler::new(e).ok_or_else(|| KernelError::Invalid(format!("{label}: an edge has no parameter range")))?; + let mut len = 0.0; + let mut last: Option = None; + for k in 0..=16 { + let p = s.at(k as f64 / 16.0).ok_or_else(|| KernelError::Invalid(format!("{label}: an edge can't be evaluated")))?; + if let Some(q) = last { + len += p.dist(q); + } + last = Some(p); + } + starts.push(total); + total += len; + edges.push(s); + } + if !(total > 0.0 && total.is_finite()) { + return Err(KernelError::Invalid(format!("{label} has no length"))); + } + Ok(EdgeChain { edges, starts, total }) + } +} + +impl NetCurve for EdgeChain { + fn point(&mut self, t: f64) -> Option { + let at = t.clamp(0.0, 1.0) * self.total; + let k = self.starts.partition_point(|s| *s <= at).saturating_sub(1); + let start = *self.starts.get(k)?; + let end = self.starts.get(k + 1).copied().unwrap_or(self.total); + let local = if end > start { ((at - start) / (end - start)).clamp(0.0, 1.0) } else { 0.0 }; + self.edges.get_mut(k)?.at(local) + } +} + +/// The faces along a chain of edges, each edge an open edge of one face of its body: their +/// frames and second fundamental forms where a boundary surface meets them (G1 and G2 +/// conditions). +pub(crate) struct EdgeFaces { + sides: Vec, +} + +/// One edge of the chain: points along it (to find the edge nearest a point), the surface of +/// its face, whether the face's boundary runs along the chain (`along`), and the surface +/// parameters found last (the next search starts there). +struct EdgeSide { + pts: Vec, + surface: mt::Surface, + along: bool, + hint: Option<(f64, f64)>, +} + +impl EdgeFaces { + pub(crate) fn new(label: &str, body: &Body, chain: &[mt::Edge]) -> Result { + let mut sides = Vec::with_capacity(chain.len()); + for e in chain { + let mut faces = body.solid.face_iter().filter(|f| f.boundaries().iter().flatten().any(|x| x.id() == e.id())); + let (Some(face), None) = (faces.next(), faces.next()) else { + return Err(KernelError::Invalid(format!( + "{label}: G1 and G2 need an open edge of a surface (an edge with a face on one side only)" + ))); + }; + // How the face's own (absolute) boundary runs this edge. + let in_wire = face.absolute_boundaries().iter().flatten().find(|x| x.id() == e.id()).cloned(); + let in_wire = in_wire.ok_or_else(|| KernelError::Failed("boundary surface: an edge's face".into()))?; + let along = in_wire.front().id() == e.front().id(); + let mut sampler = EdgeSampler::new(e).ok_or_else(|| KernelError::Invalid(format!("{label}: an edge has no parameter range")))?; + let pts = (0..=16).map(|k| sampler.at(k as f64 / 16.0)).collect::>>(); + let pts = pts.ok_or_else(|| KernelError::Invalid(format!("{label}: an edge can't be evaluated")))?; + sides.push(EdgeSide { pts, surface: face.surface(), along, hint: None }); + } + Ok(EdgeFaces { sides }) + } +} + +impl CrossFace for EdgeFaces { + fn frame(&mut self, at: Vec3, tangent: Vec3) -> Option { + // The edge nearest the point, and the chain's direction along it there. + let near = |s: &EdgeSide| s.pts.windows(2).map(|w| at.dist_to_segment(w[0], w[1])).fold(f64::INFINITY, f64::min); + let k = (0..self.sides.len()).min_by(|a, b| near(&self.sides[*a]).total_cmp(&near(&self.sides[*b])))?; + let side = self.sides.get_mut(k)?; + let seg = side.pts.windows(2).min_by(|a, b| at.dist_to_segment(a[0], a[1]).total_cmp(&at.dist_to_segment(b[0], b[1])))?; + let chain_dir = seg[1] - seg[0]; + let s = &side.surface; + let (u, v) = s.search_nearest_parameter(p3(at), side.hint, 100).or_else(|| s.search_nearest_parameter(p3(at), None, 100))?; + side.hint = Some((u, v)); + let (su, sv) = (s.uder(u, v), s.vder(u, v)); + let n = su.cross(sv); + if !(n.magnitude2() > 0.0) { + return None; + } + let n = n.normalize(); + let normal = Vec3::new(n.x, n.y, n.z); + // The face's boundary runs counter-clockwise about its surface's normal, so the face lies + // to the left of it: outward is the boundary's direction × the normal. + let wire = if side.along { chain_dir } else { -chain_dir }; + let outward = wire.cross(normal).normalized()?; + // The second fundamental form, in the basis of the tangent asked about and outward. + let t = (tangent - normal * tangent.dot(normal)).normalized()?; + let b = (outward - t * outward.dot(t)).normalized()?; + let (e, f, g) = (su.dot(su), su.dot(sv), sv.dot(sv)); + let det = e * g - f * f; + if !(det.abs() > 1e-24) { + return None; + } + // II(w, w) = L·α² + 2M·αβ + N·β² for w = α·S_u + β·S_v, with L, M, N the surface's + // second derivatives along its unit normal (do Carmo, §3-3): exact, rational surfaces + // included, and with nothing sampled beyond the edge, where the face's domain ends. + let nv = mt::Vector3::new(normal.x, normal.y, normal.z); + let (l, m, nn) = (s.uuder(u, v).dot(nv), s.uvder(u, v).dot(nv), s.vvder(u, v).dot(nv)); + let ii = |w: Vec3| -> Option { + let wv = mt::Vector3::new(w.x, w.y, w.z); + let (p, q) = (wv.dot(su), wv.dot(sv)); + let (alpha, beta) = ((p * g - q * f) / det, (q * e - p * f) / det); + let out = l * alpha * alpha + 2.0 * m * alpha * beta + nn * beta * beta; + out.is_finite().then_some(out) + }; + let (tt, bb) = (ii(t)?, ii(b)?); + let tb = (ii(t + b)? - tt - bb) / 2.0; + let form = [tt, tb, bb]; + form.iter().all(|x| x.is_finite()).then_some(FaceFrame { normal, outward: b, form }) + } +} + +/// The face of a network surface: its four sides from the surface's boundary rows, each split +/// at the `corners` that lie on it (inner vertices of edge chains). Refused when a side shrinks to +/// a point. +fn face_of(s: &NurbsSurface, corners: &[Vec3]) -> Result { + let fail = |m: &str| KernelError::Failed(format!("boundary surface: {m}")); + let [nu, nv] = s.count(); + let [ku, kv] = s.knots(); + let rows: Vec> = (0..nu).map(|i| (0..nv).map(|j| p3(s.point(i, j).xyz())).collect()).collect(); + let (uk, vk) = (mt::KnotVec::from(ku.clone()), mt::KnotVec::from(kv.clone())); + let surface = mt::BSplineSurface::try_new((uk.clone(), vk.clone()), rows.clone()).map_err(|e| fail(&e.to_string()))?; + let row = |pts: Vec, k: &mt::KnotVec| mt::BSplineCurve::try_new(k.clone(), pts).map_err(|e| fail(&e.to_string())); + let (Some(first), Some(last)) = (rows.first(), rows.last()) else { return Err(fail("control net")) }; + let bottom = row(rows.iter().filter_map(|r| r.first().copied()).collect(), &uk)?; + let top = row(rows.iter().filter_map(|r| r.last().copied()).collect(), &uk)?; + let left = row(first.clone(), &vk)?; + let right = row(last.clone(), &vk)?; + let corner = |c: &mt::BSplineCurve, end: bool| { + let (t0, t1) = c.range_tuple(); + c.subs(if end { t1 } else { t0 }) + }; + let (c00, c10, c11, c01) = (corner(&bottom, false), corner(&bottom, true), corner(&top, true), corner(&top, false)); + let size = [c00, c10, c11, c01].iter().map(|p| from_p3(*p).dist(from_p3(c00))).fold(0.0, f64::max).max(1.0); + // A side collapsed to a point (a three-sided patch) is left out of the loop: the two sides + // beside it meet at one vertex there. + let sides = [(c00, c10), (c10, c11), (c01, c11), (c00, c01)]; + let collapsed: Vec = sides.iter().map(|(a, b)| from_p3(*a).dist(from_p3(*b)) <= size * 1e-9).collect(); + match collapsed.iter().filter(|c| **c).count() { + 0 | 1 => {} + _ => return Err(KernelError::Invalid("not supported yet: a boundary surface with more than one side shrunk to a point".into())), + } + let mut verts: Vec = Vec::with_capacity(4); + for c in [c00, c10, c11, c01] { + let same = verts.iter().find(|v| from_p3(v.point()).dist(from_p3(c)) <= size * 1e-9).cloned(); + verts.push(same.unwrap_or_else(|| mt::Vertex::new(c))); + } + let (Some(v00), Some(v10), Some(v11), Some(v01)) = (verts.first(), verts.get(1), verts.get(2), verts.get(3)) else { + return Err(fail("corners")); + }; + let is = |k: usize| collapsed.get(k).copied().unwrap_or(false); + let mut wire: Vec = Vec::new(); + if !is(0) { + wire.extend(side_edges(bottom, v00, v10, corners)); + } + if !is(1) { + wire.extend(side_edges(right, v10, v11, corners)); + } + if !is(2) { + wire.extend(side_edges(top, v01, v11, corners).into_iter().rev().map(|e| e.inverse())); + } + if !is(3) { + wire.extend(side_edges(left, v00, v01, corners).into_iter().rev().map(|e| e.inverse())); + } + // u along the bottom, v up the left: counter-clockwise in (u, v). + mt::Face::try_new(vec![wire.into()], mt::Surface::BSplineSurface(surface)).map_err(|e| fail(&e.to_string())) +} + +/// The edges of one side from `a` to `b`, split at the points of `corners` that lie on it +/// (within [`CROSS_TOLERANCE`]) strictly between its ends. +fn side_edges(mut curve: mt::BSplineCurve, a: &mt::Vertex, b: &mt::Vertex, corners: &[Vec3]) -> Vec { + let (t0, t1) = curve.range_tuple(); + let margin = (t1 - t0) * 1e-6; + let mut cuts: Vec = corners + .iter() + .filter_map(|c| { + let t = curve.search_nearest_parameter(p3(*c), None, 100)?; + (t > t0 + margin && t < t1 - margin && from_p3(curve.subs(t)).dist(*c) <= CROSS_TOLERANCE).then_some(t) + }) + .collect(); + cuts.sort_by(f64::total_cmp); + cuts.dedup_by(|x, y| (*x - *y).abs() <= margin); + let mut edges = Vec::with_capacity(cuts.len() + 1); + let mut from = a.clone(); + // Cut from the far end back, so the earlier parameters stay on the part kept. + let mut tails = Vec::with_capacity(cuts.len()); + for &t in cuts.iter().rev() { + tails.push(curve.cut(t)); + } + let mut pieces = vec![curve]; + pieces.extend(tails.into_iter().rev()); + let count = pieces.len(); + for (k, piece) in pieces.into_iter().enumerate() { + let to = if k + 1 == count { b.clone() } else { mt::Vertex::new(piece.subs(piece.range_tuple().1)) }; + edges.push(mt::Edge::new(&from, &to, mt::Curve::BSplineCurve(piece))); + from = to; + } + edges +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::patch_region; + use mt::ParametricSurface3D; + use solidcraft_geom::{Loop2, Region2, Vec2}; + + fn sketch(label: &str, plane: Plane, segs: Vec) -> BoundaryCurve<'static> { + BoundaryCurve { label: label.into(), source: CurveSource::Sketch { plane, segs }, continuity: Continuity::G0 } + } + + fn line3(label: &str, a: Vec3, b: Vec3) -> BoundaryCurve<'static> { + BoundaryCurve { label: label.into(), source: CurveSource::Points(vec![a, b]), continuity: Continuity::G0 } + } + + fn face_surface(b: &Body) -> mt::Surface { + b.solid.face_iter().next().unwrap().surface() + } + + /// The flat sheet y ∈ [−10, 0], x ∈ [0, 10] in z = 0, its edge along y = 0 open. + fn sheet() -> Body { + let sq = Loop2::polygon(&[Vec2::new(0.0, -10.0), Vec2::new(10.0, -10.0), Vec2::new(10.0, 0.0), Vec2::new(0.0, 0.0)]); + patch_region(&Plane::XY, &Region2 { outer: sq, holes: vec![] }).unwrap() + } + + /// A Coons patch of four lines bounding the hyperbolic paraboloid z = x·y/10: every point of + /// the face lies on it, and its area is the paraboloid's over the square. + #[test] + fn coons_patch_of_lines_is_a_hyperbolic_paraboloid() { + let p = |x: f64, y: f64| Vec3::new(x, y, x * y / 10.0); + let body = boundary_surface( + &[line3("u0", p(0.0, 0.0), p(10.0, 0.0)), line3("u1", p(0.0, 10.0), p(10.0, 10.0))], + &[line3("v0", p(0.0, 0.0), p(0.0, 10.0)), line3("v1", p(10.0, 0.0), p(10.0, 10.0))], + ) + .unwrap(); + assert!(body.is_surface() && body.face_count() == 1); + let s = face_surface(&body); + for i in 0..=10 { + for j in 0..=10 { + let q = from_p3(s.subs(i as f64 / 10.0, j as f64 / 10.0)); + assert!((q.z - q.x * q.y / 10.0).abs() < 1e-9, "{q:?}"); + } + } + } + + /// Sketch curves: a patch inside two arcs and two lines on one plane lies in that plane. + #[test] + fn coplanar_sketch_curves_give_a_plane() { + use std::f64::consts::PI; + // Two half-discs' arcs over [0, 10] (radius 5) facing out, joined by the lines x = 0, 10. + let low = Seg2::Arc { center: Vec2::new(5.0, 0.0), radius: 5.0, start: PI, sweep: PI }; + let high = Seg2::Arc { center: Vec2::new(5.0, 8.0), radius: 5.0, start: PI, sweep: -PI }; + let left = Seg2::Line { a: Vec2::new(0.0, 0.0), b: Vec2::new(0.0, 8.0) }; + let right = Seg2::Line { a: Vec2::new(10.0, 0.0), b: Vec2::new(10.0, 8.0) }; + let pl = Plane::XY.offset(3.0); + let body = boundary_surface( + &[sketch("low", pl, vec![low]), sketch("high", pl, vec![high])], + &[sketch("left", pl, vec![left]), sketch("right", pl, vec![right])], + ) + .unwrap(); + let s = face_surface(&body); + for i in 0..=8 { + for j in 0..=8 { + assert!((s.subs(i as f64 / 8.0, j as f64 / 8.0).z - 3.0).abs() < 1e-9); + } + } + } + + /// A stadium: two half-circle arcs (radius 5) bulging out of the 10 × 8 rectangle between two + /// lines. Each arc meets its lines at 180°, so the patch is singular at its corners, and the + /// reference mesher's walk round the boundary used to run from such a corner onto the + /// surface's extrapolation past its domain (a planar spline passes through every point of its + /// plane there): the face measured 289 mm² for its 158.5. The vendored mesher's bounded + /// projection (`nearest_projection`, vendor/README.md) keeps the walk on the domain. Its area + /// is the stadium's, 80 + 25π, by the faces, the measuring mesh and the display mesh. + #[test] + fn a_patch_with_straight_corners_measures_its_own_area() { + use std::f64::consts::PI; + let low = Seg2::Arc { center: Vec2::new(5.0, 0.0), radius: 5.0, start: PI, sweep: PI }; + let high = Seg2::Arc { center: Vec2::new(5.0, 8.0), radius: 5.0, start: PI, sweep: -PI }; + let left = Seg2::Line { a: Vec2::new(0.0, 0.0), b: Vec2::new(0.0, 8.0) }; + let right = Seg2::Line { a: Vec2::new(10.0, 0.0), b: Vec2::new(10.0, 8.0) }; + let pl = Plane::XY.offset(3.0); + let body = boundary_surface( + &[sketch("low", pl, vec![low]), sketch("high", pl, vec![high])], + &[sketch("left", pl, vec![left]), sketch("right", pl, vec![right])], + ) + .unwrap(); + let want = 80.0 + 25.0 * PI; + // Chords of the arcs at 0.01 mm sag cut off up to (2/3)·sag·length ≈ 0.21 mm². + let near = |a: f64| (a - want).abs() < 0.25; + let faces = body.faces(0.01).unwrap(); + assert!(faces.len() == 1 && near(faces[0].area), "faces: {faces:?}"); + let tess = body.tessellate(0.01).unwrap().measure().area; + assert!(near(tess), "measuring mesh: {tess}"); + let display = body.display_mesh(0.01).unwrap().measure().area; + assert!(near(display), "display mesh: {display}"); + } + + /// Three sides: a quarter disc of radius 10 (two radii and the arc) is one face with three + /// edges, the collapsed side left out of its loop, and measures 25π by its faces and both + /// meshes. + #[test] + fn a_three_sided_patch_is_one_face_with_three_edges() { + use std::f64::consts::{FRAC_PI_2, PI}; + let r = 10.0; + let arc = Seg2::Arc { center: Vec2::new(0.0, 0.0), radius: r, start: 0.0, sweep: FRAC_PI_2 }; + let x = Seg2::Line { a: Vec2::new(0.0, 0.0), b: Vec2::new(r, 0.0) }; + let y = Seg2::Line { a: Vec2::new(0.0, r), b: Vec2::new(0.0, 0.0) }; + let pl = Plane::XY; + let apex = BoundaryCurve { label: "centre".into(), source: CurveSource::Point(Vec3::ZERO), continuity: Continuity::G0 }; + let body = boundary_surface(&[sketch("arc", pl, vec![arc]), apex], &[sketch("x", pl, vec![x]), sketch("y", pl, vec![y])]).unwrap(); + assert!(body.is_surface() && body.face_count() == 1); + let edges: std::collections::HashSet<_> = body.solid.edge_iter().map(|e| e.id()).collect(); + assert_eq!(edges.len(), 3, "the collapsed side is no edge"); + let want = PI * r * r / 4.0; + let near = |a: f64| (a - want).abs() < 0.2; + let faces = body.faces(0.01).unwrap(); + assert!(near(faces[0].area), "faces: {faces:?}"); + let tess = body.tessellate(0.01).unwrap().measure().area; + assert!(near(tess), "measuring mesh: {tess}"); + let display = body.display_mesh(0.01).unwrap().measure().area; + assert!(near(display), "display mesh: {display}"); + // Two collapsed sides are no patch. + let p = |l: &str, q: Vec3| BoundaryCurve { label: l.into(), source: CurveSource::Point(q), continuity: Continuity::G0 }; + assert!( + boundary_surface(&[sketch("arc", pl, vec![arc]), p("o", Vec3::ZERO)], &[p("a", Vec3::new(r, 0.0, 0.0)), sketch("y", pl, vec![y])]) + .is_err() + ); + assert!(boundary_surface(&[sketch("arc", pl, vec![arc]), p("nan", Vec3::new(f64::NAN, 0.0, 0.0))], &[]).is_err()); + } + + /// The angle between the two faces' normals across the edge at `p`. + fn normal_angle(a: &mt::Surface, b: &mt::Surface, p: Vec3) -> f64 { + let n = |s: &mt::Surface| { + let (u, v) = s.search_nearest_parameter(p3(p), None, 100).unwrap(); + let n = s.normal(u, v); + Vec3::new(n.x, n.y, n.z) + }; + n(a).dot(n(b)).abs().clamp(0.0, 1.0).acos().to_degrees() + } + + /// G1 against a flat sheet's open edge: a loft from the edge to a line above, and a patch + /// whose sides are arcs leaving the sheet tangent to it. The normals across the shared edge + /// agree within 0.1°; the patch joins the sheet into one surface body. + #[test] + fn g1_meets_the_sheet_tangent() { + let base = sheet(); + let edge = + |g| BoundaryCurve { label: "edge".into(), source: CurveSource::Edges { body: &base, at: vec![Vec3::new(5.0, 0.0, 0.0)] }, continuity: g }; + let far = || line3("far", Vec3::new(0.0, 10.0, 5.0), Vec3::new(10.0, 10.0, 5.0)); + let loft = boundary_surface(&[edge(Continuity::G1), far()], &[]).unwrap(); + let plain = boundary_surface(&[edge(Continuity::G0), far()], &[]).unwrap(); + let (sb, sl, sp) = (face_surface(&base), face_surface(&loft), face_surface(&plain)); + for k in 0..=10 { + let p = Vec3::new(k as f64, 0.0, 0.0); + let a = normal_angle(&sb, &sl, p); + assert!(a < 0.1, "G1 at x = {k}: {a}°"); + } + assert!(normal_angle(&sb, &sp, Vec3::new(5.0, 0.0, 0.0)) > 10.0, "G0 has a crease"); + // Arcs in x = 0 and x = 10 from the edge, tangent to the sheet, up to the far line. + let arc_plane = |x: f64| Plane::new(Vec3::new(x, 0.0, 0.0), Vec3::Y, Vec3::Z).unwrap(); + let arc = Seg2::Arc { center: Vec2::new(0.0, 12.5), radius: 12.5, start: -std::f64::consts::FRAC_PI_2, sweep: 0.927_295_218_001_612_2 }; + let patch = boundary_surface( + &[edge(Continuity::G1), far()], + &[sketch("left", arc_plane(0.0), vec![arc]), sketch("right", arc_plane(10.0), vec![arc])], + ) + .unwrap(); + let spa = face_surface(&patch); + for k in 0..=10 { + let a = normal_angle(&sb, &spa, Vec3::new(k as f64, 0.0, 0.0)); + assert!(a < 0.1, "G1 patch at x = {k}: {a}°"); + } + let joined = boundary_join(&base, &patch).unwrap(); + assert_eq!(joined.face_count(), 2); + assert!( + boundary_join(&base, &boundary_surface(&[far(), line3("other", Vec3::new(0.0, 20.0, 5.0), Vec3::new(10.0, 20.0, 5.0))], &[]).unwrap()) + .is_err() + ); + } + + /// Normal curvature of `s` at the point nearest `p` in the direction `w` (projected into the + /// tangent plane), signed against the unit normal `n`. + fn normal_curvature(s: &mt::Surface, p: Vec3, w: Vec3, n: Vec3) -> f64 { + let (u, v) = s.search_nearest_parameter(p3(p), None, 100).unwrap(); + let v3 = |x: mt::Vector3| Vec3::new(x.x, x.y, x.z); + let (su, sv) = (v3(s.uder(u, v)), v3(s.vder(u, v))); + let (suu, suv, svv) = (v3(s.uuder(u, v)), v3(s.uvder(u, v)), v3(s.vvder(u, v))); + let (e, f, g) = (su.dot(su), su.dot(sv), sv.dot(sv)); + let (l, m, nn) = (suu.dot(n), suv.dot(n), svv.dot(n)); + let (p1, q1) = (w.dot(su), w.dot(sv)); + let det = e * g - f * f; + let (a, b) = ((p1 * g - q1 * f) / det, (q1 * e - p1 * f) / det); + (l * a * a + 2.0 * m * a * b + nn * b * b) / (e * a * a + 2.0 * f * a * b + g * b * b) + } + + /// G2 against a real cylinder's open edge (R = 10, the face below the edge): across the edge + /// the surface's normal curvature is the cylinder's, −1/R about its outward normal, within + /// 0.1 % — also when the far curve is sheared along the edge, so that the surface's own + /// cross direction leans. With G1 alone it is not. + #[test] + fn g2_matches_a_cylinder_across_its_edge() { + let r = 10.0; + let pl = Plane::new(Vec3::ZERO, Vec3::Y, Vec3::Z).unwrap(); + let arc = Seg2::Arc { center: Vec2::new(0.0, 0.0), radius: r, start: -std::f64::consts::FRAC_PI_2, sweep: std::f64::consts::FRAC_PI_2 }; + let cyl = crate::extrude_sheet(&pl, &[vec![arc]], 0.0, 10.0).unwrap(); + let sc = face_surface(&cyl); + let made = |g, shear: f64| { + let edge = + BoundaryCurve { label: "rim".into(), source: CurveSource::Edges { body: &cyl, at: vec![Vec3::new(5.0, r, 0.0)] }, continuity: g }; + let top = line3("top", Vec3::new(shear, 3.0, 12.0), Vec3::new(10.0 + shear, 3.0, 12.0)); + face_surface(&boundary_surface(&[edge, top], &[]).unwrap()) + }; + for shear in [0.0, 5.0] { + let sm = made(Continuity::G2, shear); + for k in 1..=9 { + let p = Vec3::new(k as f64, r, 0.0); + let kc = normal_curvature(&sc, p, Vec3::Z, Vec3::Y); + assert!((kc + 1.0 / r).abs() < 1e-9, "the cylinder's own curvature: {kc}"); + let km = normal_curvature(&sm, p, Vec3::Z, Vec3::Y); + assert!((km - kc).abs() < 1e-3 / r, "shear {shear}, x = {k}: {km} across the edge where the cylinder has {kc}"); + // Along the edge both are straight. + assert!(normal_curvature(&sm, p, Vec3::X, Vec3::Y).abs() < 1e-6); + assert!(normal_angle(&sc, &sm, p) < 0.1); + } + } + let g1 = made(Continuity::G1, 0.0); + let k1 = normal_curvature(&g1, Vec3::new(5.0, r, 0.0), Vec3::Z, Vec3::Y); + assert!((k1 + 1.0 / r).abs() > 0.1 / r, "G1 alone does not carry the curvature: {k1}"); + } + + /// G3 is more than a network holds: refused, never a panic. + #[test] + fn g3_is_not_supported() { + let base = sheet(); + let edge = BoundaryCurve { + label: "edge".into(), + source: CurveSource::Edges { body: &base, at: vec![Vec3::new(5.0, 0.0, 0.0)] }, + continuity: Continuity::G3, + }; + let e = boundary_surface(&[edge, line3("far", Vec3::new(0.0, 10.0, 5.0), Vec3::new(10.0, 10.0, 5.0))], &[]).unwrap_err(); + assert!(e.to_string().contains("not supported yet"), "{e}"); + } + + /// Refusals: G1 on a sketch curve, on an edge between two faces, curves that miss. + #[test] + fn bad_boundaries_are_refused() { + let cube = crate::box_solid(Vec3::ZERO, Vec3::new(10.0, 10.0, 10.0)).unwrap(); + let edge = BoundaryCurve { + label: "edge".into(), + source: CurveSource::Edges { body: &cube, at: vec![Vec3::new(5.0, 0.0, 0.0)] }, + continuity: Continuity::G1, + }; + let far = line3("far", Vec3::new(0.0, 10.0, 15.0), Vec3::new(10.0, 10.0, 15.0)); + let e = boundary_surface(&[edge, far], &[]).unwrap_err(); + assert!(e.to_string().contains("open edge"), "{e}"); + let mut tangent = line3("t", Vec3::ZERO, Vec3::X); + tangent.continuity = Continuity::G1; + assert!(boundary_surface(&[tangent, line3("b", Vec3::Y, Vec3::new(1.0, 1.0, 0.0))], &[]).is_err()); + let miss = boundary_surface( + &[line3("a", Vec3::ZERO, Vec3::X * 10.0), line3("b", Vec3::Y * 10.0, Vec3::new(10.0, 10.0, 0.0))], + &[line3("c", Vec3::ZERO, Vec3::Y * 10.0), line3("d", Vec3::new(10.0, 0.0, 2.0), Vec3::new(10.0, 10.0, 2.0))], + ) + .unwrap_err(); + assert!(miss.to_string().contains("don't meet"), "{miss}"); + let pts = BoundaryCurve { label: "p".into(), source: CurveSource::Points(vec![Vec3::ZERO; 3]), continuity: Continuity::G0 }; + assert!(boundary_surface(&[pts, line3("b", Vec3::Y, Vec3::new(1.0, 1.0, 0.0))], &[]).is_err()); + let nan = BoundaryCurve { + label: "n".into(), + source: CurveSource::Points(vec![Vec3::ZERO, Vec3::new(f64::NAN, 0.0, 0.0)]), + continuity: Continuity::G0, + }; + assert!(boundary_surface(&[nan, line3("b", Vec3::Y, Vec3::new(1.0, 1.0, 0.0))], &[]).is_err()); + } +} diff --git a/crates/kernel/src/lib.rs b/crates/kernel/src/lib.rs index 35eb423e..f29b7dd9 100644 --- a/crates/kernel/src/lib.rs +++ b/crates/kernel/src/lib.rs @@ -25,6 +25,7 @@ mod boolean_guard; mod boolean_record; mod bore_guard; mod boundary; +mod boundary_network; mod bridge; mod build; mod clearance; @@ -121,6 +122,8 @@ mod surface_extend; mod surface_kind; mod surface_offset; mod surface_ruled; +#[cfg(test)] +mod surface_second_derivative_tests; mod surface_untrim; mod surfaces; mod sweep_path; @@ -145,6 +148,7 @@ pub use blend_variable::{RadiusSpline, Taper, fillet_chord, fillet_tapered, fill 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 boundary_network::{BoundaryCurve, CurveSource, MAX_CHAIN_CURVES, MAX_EDGE_PICKS, boundary_join, boundary_surface}; 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, path_polyline, revolve, diff --git a/crates/kernel/src/surface_second_derivative_tests.rs b/crates/kernel/src/surface_second_derivative_tests.rs new file mode 100644 index 00000000..76265437 --- /dev/null +++ b/crates/kernel/src/surface_second_derivative_tests.rs @@ -0,0 +1,96 @@ +//! The vendored truck-geometry's second derivatives of B-spline and NURBS surfaces (`uuder`, +//! `vvder`) against finite differences of the surface itself, for degrees 2 to 5. Upstream's +//! were short of the factor p − 1, so every curvature read from a cubic surface came out half its +//! value (#763: a G2 boundary surface measured half the cylinder's curvature it does carry). +//! These pin the vendored fix, so a refresh of the crate can't quietly bring the error back. + +use solidcraft_geom::Vec3; +use truck_modeling as mt; +use truck_modeling::ParametricSurface; + +use crate::body::from_p3; + +/// A degree-`p` × degree-`q` surface on a 5 × 6 net with an interior knot each way, bent in z. +fn surface(p: usize, q: usize) -> mt::BSplineSurface { + let knots = |d: usize, n: usize| { + let inner: Vec = (1..n - d).map(|k| k as f64 / (n - d) as f64).collect(); + let mut k = vec![0.0; d + 1]; + k.extend(inner); + k.extend(vec![1.0; d + 1]); + mt::KnotVec::from(k) + }; + let (nu, nv) = (p + 3, q + 3); + let rows: Vec> = (0..nu) + .map(|i| { + (0..nv) + .map(|j| { + let (x, y) = (i as f64 * 2.0, j as f64 * 1.5); + mt::Point3::new(x + 0.3 * y, y, 0.4 * x * x - 0.2 * y * y + 0.1 * x * y) + }) + .collect() + }) + .collect(); + mt::BSplineSurface::new((knots(p, nu), knots(q, nv)), rows) +} + +/// Second differences of `f` along the parameter line, central inside and one-sided at the ends. +fn second(f: impl Fn(f64) -> Vec3, t: f64) -> Vec3 { + let h = 1e-4; + if t < 2.0 * h { + (f(t) * 2.0 - f(t + h) * 5.0 + f(t + 2.0 * h) * 4.0 - f(t + 3.0 * h)) / (h * h) + } else if t > 1.0 - 2.0 * h { + (f(t) * 2.0 - f(t - h) * 5.0 + f(t - 2.0 * h) * 4.0 - f(t - 3.0 * h)) / (h * h) + } else { + (f(t - h) - f(t) * 2.0 + f(t + h)) / (h * h) + } +} + +fn v3(x: mt::Vector3) -> Vec3 { + Vec3::new(x.x, x.y, x.z) +} + +#[test] +fn bspline_surface_second_derivatives_match_finite_differences() { + for p in 2..=5 { + for q in 2..=5 { + let s = surface(p, q); + for (u, v) in [(0.0, 0.0), (0.37, 0.0), (0.5, 0.61), (1.0, 0.2), (0.13, 1.0)] { + let uu = second(|t| from_p3(s.subs(t, v)), u); + let vv = second(|t| from_p3(s.subs(u, t)), v); + let scale = uu.len().max(vv.len()).max(1.0); + assert!(v3(s.uuder(u, v)).dist(uu) < 1e-4 * scale, "degree ({p}, {q}) uuder at ({u}, {v}): {:?} vs {uu:?}", s.uuder(u, v)); + assert!(v3(s.vvder(u, v)).dist(vv) < 1e-4 * scale, "degree ({p}, {q}) vvder at ({u}, {v}): {:?} vs {vv:?}", s.vvder(u, v)); + } + } + } +} + +/// A rational surface's second derivatives come from its homogeneous B-spline's: they are right +/// too. +#[test] +fn nurbs_surface_second_derivatives_match_finite_differences() { + let b = surface(3, 3); + let (ku, kv) = b.knot_vecs().clone(); + let rows: Vec> = b + .control_points() + .iter() + .enumerate() + .map(|(i, row)| { + row.iter() + .enumerate() + .map(|(j, p)| { + let w = 1.0 + 0.3 * ((i + 2 * j) % 3) as f64; + mt::Vector4::new(p.x * w, p.y * w, p.z * w, w) + }) + .collect() + }) + .collect(); + let s = mt::NurbsSurface::new(mt::BSplineSurface::new((ku, kv), rows)); + for (u, v) in [(0.0, 0.0), (0.37, 0.0), (0.5, 0.61), (1.0, 0.2)] { + let uu = second(|t| from_p3(s.subs(t, v)), u); + let vv = second(|t| from_p3(s.subs(u, t)), v); + let scale = uu.len().max(vv.len()).max(1.0); + assert!(v3(s.uuder(u, v)).dist(uu) < 1e-4 * scale, "uuder at ({u}, {v}): {:?} vs {uu:?}", s.uuder(u, v)); + assert!(v3(s.vvder(u, v)).dist(vv) < 1e-4 * scale, "vvder at ({u}, {v}): {:?} vs {vv:?}", s.vvder(u, v)); + } +} diff --git a/crates/kernel/src/surfaces.rs b/crates/kernel/src/surfaces.rs index 82416f0c..70f86215 100644 --- a/crates/kernel/src/surfaces.rs +++ b/crates/kernel/src/surfaces.rs @@ -48,7 +48,7 @@ fn surface_of(shell: mt::Shell) -> Body { /// A body from faces: a solid when they close, else a surface; with the history of the /// operation that made the faces from `operands` as `rec` records. -fn body_of(faces: Vec, operands: &[&mt::Solid], rec: &Rebuild) -> Result { +pub(crate) fn body_of(faces: Vec, operands: &[&mt::Solid], rec: &Rebuild) -> Result { let shell = mt::Shell::from(faces); if shell.shell_condition() == mt::ShellCondition::Closed { crate::record::rebuilt(mt::Solid::new_unchecked(vec![shell]), operands, rec) diff --git a/crates/ui-egui/src/boundary_dialog_tests.rs b/crates/ui-egui/src/boundary_dialog_tests.rs new file mode 100644 index 00000000..561470e7 --- /dev/null +++ b/crates/ui-egui/src/boundary_dialog_tests.rs @@ -0,0 +1,127 @@ +//! The Boundary Surface dialog (#763): a loop of sketch curves becomes the loop patch, edges and +//! curves across become a network with its continuity and operation, curves that meet end to end +//! are one chain, and Edit Feature brings a network back as it was. +use crate::dialogs::Kind; +use crate::dialogs_surface::Surf; +use crate::scenario::Harness; +use serde_json::{Map, json}; +use solidcraft_engine::Sel; +use solidcraft_engine::doc::{FeatureKind, Operation, kernel::Continuity}; +use solidcraft_engine::geom::Vec3; + +fn ok(h: &mut Harness) { + h.frames(2); + h.step(&json!({"widget":"OK"})).unwrap(); + // An edit is applied off the UI thread: give it the frames it takes. + for _ in 0..200 { + if h.app.dialog.is_none() { + break; + } + h.frames(1); + } + assert!(h.app.dialog.is_none(), "{:?} {:?} {:?}", h.app.preview.error, h.app.dialog.as_ref().map(|d| (&d.error, &d.inputs)), h.app.session.log); +} + +#[test] +fn a_loop_of_sketch_curves_is_the_loop_patch() { + let mut h = Harness::new(); + h.app.run("sketch.create", json!({"plane":"XY","name":"Loop"})).unwrap(); + for (id, a, b) in [("s1", [0, 0], [10, 0]), ("s2", [10, 0], [10, 10]), ("s3", [10, 10], [0, 10]), ("s4", [0, 10], [0, 0])] { + h.app.run("sketch.line", json!({"points":[a, b],"ids":[id]})).unwrap(); + } + h.app.run("sketch.finish", json!({})).unwrap(); + h.step(&json!({"call":"ui.start","params":{"command":"surface.boundary"}})).unwrap(); + h.frames(3); + assert!(h.widget("Cross curves").is_some() && h.widget("Continuity").is_some() && h.widget("Operation").is_some()); + let d = h.app.dialog.as_mut().unwrap(); + d.inputs[0].items = ["s1", "s2", "s3", "s4"].iter().map(|id| Sel::SketchCurve { id: (*id).into() }).collect(); + ok(&mut h); + let f = h.app.session.doc.features.last().unwrap(); + let FeatureKind::SurfaceBoundary { curves, orders, .. } = &f.kind else { panic!("{:?}", f.kind) }; + assert_eq!(curves, &["s1", "s2", "s3", "s4"]); + assert!(orders.iter().all(|o| *o == Continuity::G1), "{orders:?}"); + assert_eq!(h.app.session.world_state().bodies.len(), 1); +} + +/// A quarter cylinder sheet (its top edge open at y = 10, z = 0) and a 3D line above it. +fn cylinder_and_line(h: &mut Harness) -> String { + h.app.run("sketch.create", json!({"plane":"YZ","name":"Arc"})).unwrap(); + let c = 10.0 * std::f64::consts::FRAC_1_SQRT_2; + let r = h.app.run("sketch.arc.three_point", json!({"start":[0,-10],"end":[10,0],"through":[c,-c]})).unwrap(); + let arc = r["created"]["curves"][0].as_str().unwrap().to_string(); + h.app.run("sketch.finish", json!({})).unwrap(); + h.app.run("surface.extrude", json!({"sketch":"Arc","curves":[arc],"distance":10})).unwrap(); + h.app.run("sketch.create", json!({"plane":"XY","name":"Top"})).unwrap(); + h.app.run("sketch.options", json!({"sketch_3d": true})).unwrap(); + let r = h.app.run("sketch.line3d", json!({"points":[[0,3,12],[10,3,12]]})).unwrap(); + h.app.run("sketch.finish", json!({})).unwrap(); + r["curves"][0].as_str().unwrap().to_string() +} + +#[test] +fn edges_and_curves_make_a_network_and_come_back_in_an_edit() { + let mut h = Harness::new(); + let top = cylinder_and_line(&mut h); + let rim = crate::dialogs::edge_sel(&h.app.session, Vec3::new(5.0, 10.0, 0.0)).unwrap(); + let Sel::Edge { body, .. } = &rim else { panic!("{rim:?}") }; + let body = body.clone(); + h.step(&json!({"call":"ui.start","params":{"command":"surface.boundary"}})).unwrap(); + h.frames(3); + let d = h.app.dialog.as_mut().unwrap(); + d.inputs[0].items = vec![rim.clone(), Sel::SketchCurve { id: top.clone() }]; + let Kind::Surface(Surf::Boundary { continuity, join }) = &mut d.kind else { panic!("{:?}", d.kind) }; + *continuity = 2; + *join = true; + ok(&mut h); + let name = h.app.session.doc.features.last().unwrap().name.clone(); + let network = |h: &Harness| match &h.app.session.doc.find_feature(&name).unwrap().kind { + FeatureKind::SurfaceNetwork { u, v, operation } => (u.clone(), v.clone(), *operation), + k => panic!("{k:?}"), + }; + let (u, v, operation) = network(&h); + assert_eq!((u.len(), v.len(), operation), (2, 0, Operation::Join)); + assert_eq!(u[0].body.as_ref().map(|b| b.to_string()), Some(body.clone())); + assert_eq!(u[0].continuity, Continuity::G2); + assert_eq!(u[1].curves, vec![top.clone()]); + // Joined into the sheet: one body, two faces. + let st = h.app.session.world_state(); + assert_eq!(st.bodies.len(), 1); + assert_eq!(st.body(&body).map(|b| b.body.face_count()), Some(2)); + drop(st); + h.step(&json!({"call":"ui.editFeature","params":{"feature":name}})).unwrap(); + h.frames(8); + let d = h.app.dialog.as_ref().unwrap(); + assert!(matches!(d.kind, Kind::Surface(Surf::Boundary { continuity: 2, join: true })), "{:?}", d.kind); + assert_eq!(d.inputs[0].items.len(), 2, "{:?}", d.inputs[0].items); + ok(&mut h); + let (u2, v2, op2) = network(&h); + assert_eq!((u2.len(), v2.len(), op2), (2, 0, Operation::Join)); + assert_eq!((u2[0].continuity, &u2[1].curves), (Continuity::G2, &vec![top])); +} + +/// Picks that meet end to end in one sketch are one curve of the network; a pick elsewhere is +/// another. +#[test] +fn curves_that_meet_are_one_chain() { + let mut h = Harness::new(); + h.app.run("sketch.create", json!({"plane":"XY","name":"Net"})).unwrap(); + h.app.run("sketch.line", json!({"points":[[0,0],[5,0]],"ids":["a1"]})).unwrap(); + h.app.run("sketch.line", json!({"points":[[5,0],[10,0]],"ids":["a2"]})).unwrap(); + h.app.run("sketch.line", json!({"points":[[0,8],[10,8]],"ids":["b"]})).unwrap(); + h.app.run("sketch.line", json!({"points":[[0,0],[0,8]],"ids":["c"]})).unwrap(); + h.app.run("sketch.line", json!({"points":[[10,0],[10,8]],"ids":["d"]})).unwrap(); + h.app.run("sketch.finish", json!({})).unwrap(); + let sketch = h.app.session.doc.find_feature("Net").unwrap().id; + let k = Surf::Boundary { continuity: 0, join: false }; + let (_, mut inputs) = crate::dialogs_surface::start(&h.app, "surface.boundary").unwrap(); + inputs[0].items = ["a2", "b", "a1"].iter().map(|id| Sel::SketchCurve { id: (*id).into() }).collect(); + inputs[1].items = ["c", "d"].iter().map(|id| Sel::SketchCurve { id: (*id).into() }).collect(); + let cmds = crate::dialogs_surface::commands(&h.app, &k, &inputs, &Map::new()).unwrap(); + let (id, p) = &cmds[0]; + assert_eq!(id, "surface.boundary"); + assert_eq!(p["u"], json!([{"sketch": sketch, "curves": ["a2", "a1"]}, {"sketch": sketch, "curves": ["b"]}]), "{p}"); + assert_eq!(p["v"], json!([{"sketch": sketch, "curves": ["c"]}, {"sketch": sketch, "curves": ["d"]}]), "{p}"); + // And it builds: the chain a1–a2 is one boundary of the patch. + h.app.run("surface.boundary", p.clone()).unwrap(); + assert_eq!(h.app.session.world_state().bodies.len(), 1); +} diff --git a/crates/ui-egui/src/dialogs_surface.rs b/crates/ui-egui/src/dialogs_surface.rs index 2c5375fc..a984cf3e 100644 --- a/crates/ui-egui/src/dialogs_surface.rs +++ b/crates/ui-egui/src/dialogs_surface.rs @@ -1,5 +1,6 @@ //! Dialogs for the surface tools and face splitting (#19): Surface Extrude, Revolve, Sweep and -//! Loft, Offset and Ruled Surface, Patch, Stitch, Thicken, Extend, Trim, Untrim and Split Face. Each has a live preview and Edit Feature. +//! Loft, Offset and Ruled Surface, Boundary Surface, Patch, Stitch, Thicken, Extend, Trim, Untrim +//! and Split Face. Each has a live preview and Edit Feature. use serde_json::{Map, Value, json}; use solidcraft_engine::Sel; @@ -27,6 +28,13 @@ pub enum Surf { }, Sweep(crate::dialogs_sweep::SweepOpts), Loft, + /// `continuity`: an index into [`PATCH_CONTINUITY`] (for a loop of sketch curves, how much of + /// the sketch's plane each side repeats; for picked edges, how the surface meets the faces + /// beside them); `join`: sew the surface into the body of its edges. + Boundary { + continuity: usize, + join: bool, + }, /// `continuity`: an index into [`PATCH_CONTINUITY`] (how a patch on edges meets the faces /// beside its loop). Patch { @@ -74,6 +82,7 @@ impl Surf { Surf::Revolve { .. } => "SURFACE REVOLVE", Surf::Sweep(_) => "SURFACE SWEEP", Surf::Loft => "SURFACE LOFT", + Surf::Boundary { .. } => "BOUNDARY SURFACE", Surf::Patch { .. } => "PATCH", Surf::Unstitch => "UNSTITCH", Surf::ReverseNormal => "REVERSE NORMAL", @@ -104,6 +113,9 @@ impl Surf { /// How a patch on edges meets the faces beside its loop (labels, and the `continuity` parameter). const PATCH_CONTINUITY: [(&str, &str); 3] = [("Position (G0)", "g0"), ("Tangent (G1)", "g1"), ("Curvature (G2)", "g2")]; +/// What a boundary surface becomes: a body of its own, or part of the body of its edges. +const BOUNDARY_OPS: [&str; 2] = ["New Body", "Join"]; + /// Which holes Untrim takes away: every one, the picked ones, or none (the outer boundary only). const UNTRIM_HOLES: [&str; 3] = ["All", "Picked", "None"]; @@ -128,6 +140,7 @@ fn inputs(k: &Surf) -> Vec { ] } Surf::Loft => vec![SelInput::new("Sections (one sketch each)", PROFILES | CURVES, true)], + Surf::Boundary { .. } => vec![SelInput::new("Curves", CURVES | EDGES, true), SelInput::new("Cross curves", CURVES | EDGES, true)], Surf::Patch { .. } => vec![SelInput::new("Boundary", PROFILES | EDGES, true), SelInput::new("Pass through", VERTICES | FACES, true)], Surf::Unstitch => vec![SelInput::new("Body, or faces to take out", BODIES | FACES, true)], Surf::ReverseNormal => vec![SelInput::new("Surface, or faces to turn over", BODIES | FACES, true)], @@ -155,6 +168,7 @@ pub fn start(_app: &SolveApp, id: &str) -> Option<(Kind, Vec)> { "surface.revolve" => Surf::Revolve { angle: "360 deg".into(), symmetric: false, axis: None }, "surface.sweep" => Surf::Sweep(crate::dialogs_sweep::SweepOpts::default()), "surface.loft" => Surf::Loft, + "surface.boundary" => Surf::Boundary { continuity: 1, join: false }, "surface.patch" => Surf::Patch { continuity: 0 }, "surface.unstitch" => Surf::Unstitch, "surface.reverse_normal" => Surf::ReverseNormal, @@ -264,6 +278,16 @@ pub fn rows(ui: &mut egui::Ui, k: &mut Surf) -> bool { combo(ui, "pa_cont", &PATCH_CONTINUITY.map(|(label, _)| label), continuity); ui.end_row(); } + Surf::Boundary { continuity, join } => { + row_label(ui, "Continuity"); + combo(ui, "bo_cont", &PATCH_CONTINUITY.map(|(label, _)| label), continuity); + ui.end_row(); + row_label(ui, "Operation"); + let mut op = usize::from(*join); + combo(ui, "bo_op", &BOUNDARY_OPS, &mut op); + *join = op == 1; + ui.end_row(); + } Surf::Sweep(opts) => enter |= crate::dialogs_sweep::rows(ui, opts, "ss"), Surf::Trim | Surf::SplitFace | Surf::Unstitch | Surf::ReverseNormal | Surf::Loft => {} } @@ -342,6 +366,98 @@ fn tool(app: &SolveApp, sel: Option<&Sel>) -> Result<(&'static str, Value), Stri } } +/// One curve of a boundary surface as the dialog picks it: sketch curves of one sketch, or edges +/// of one body, that meet end to end. +#[derive(Clone, Debug, PartialEq)] +enum Chain { + Sketch { sketch: u64, ids: Vec }, + Edges { body: String, points: Vec }, +} + +/// The picks grouped into curves: picks of one sketch (or one body) whose ends meet run on into +/// one chain, in the order they were first picked; anything else is a curve of its own. +fn chains(app: &SolveApp, picks: &[Sel]) -> Result, String> { + let st = app.session.world_state(); + // Each pick: what it belongs to (a sketch or a body), its ends, and how it joins a chain. + let mut parts: Vec<(Chain, [Vec3; 2])> = Vec::with_capacity(picks.len()); + for x in picks { + match x { + Sel::SketchCurve { id } => { + let sketch = crate::dialogs::curves_sketch(app, std::slice::from_ref(id), None).ok_or("a picked curve is no longer in a sketch")?; + let ss = st.sketch(sketch).ok_or("a picked curve's sketch has not been solved")?; + let ends = if let Some(ci) = ss.sketch.curve_index(id) { + let segs = ss.sketch.segs(ci); + match (segs.first(), segs.last()) { + (Some(a), Some(b)) => [ss.plane.to_world(a.start()), ss.plane.to_world(b.end())], + _ => return Err(format!("curve `{id}` has no geometry")), + } + } else { + let w = ss.sketch.wire_index(id).and_then(|wi| ss.sketch.wires.get(wi)).ok_or("a picked curve is no longer in a sketch")?; + match (w.pts.first(), w.pts.last()) { + (Some(a), Some(b)) => [*a, *b], + _ => return Err(format!("curve `{id}` has no geometry")), + } + }; + parts.push((Chain::Sketch { sketch, ids: vec![id.clone()] }, ends)); + } + Sel::Edge { body, index, point } => { + let b = st.body(body).ok_or("a picked edge's body is gone")?; + let m = b.mesh(); + let e = m.edges.get(*index).ok_or("a picked edge is gone")?; + let (Some(a), Some(z)) = (e.first(), e.last()) else { return Err("a picked edge has no points".into()) }; + parts.push((Chain::Edges { body: body.clone(), points: vec![*point] }, [*a, *z])); + } + _ => return Err("pick sketch curves or edges".into()), + } + } + let same_source = |a: &Chain, b: &Chain| match (a, b) { + (Chain::Sketch { sketch: x, .. }, Chain::Sketch { sketch: y, .. }) => x == y, + (Chain::Edges { body: x, .. }, Chain::Edges { body: y, .. }) => x == y, + _ => false, + }; + let meet = |a: &[Vec3; 2], b: &[Vec3; 2]| a.iter().any(|p| b.iter().any(|q| p.dist(*q) <= 1e-6 * (1.0 + p.len()))); + // Union the picks that meet, keeping the first pick of each chain in front. + let mut group: Vec = (0..parts.len()).collect(); + fn root(group: &[usize], mut i: usize) -> usize { + while let Some(&g) = group.get(i) { + if g == i { + break; + } + i = g; + } + i + } + for i in 0..parts.len() { + for j in 0..i { + let (Some((ci, ei)), Some((cj, ej))) = (parts.get(i), parts.get(j)) else { continue }; + if same_source(ci, cj) && meet(ei, ej) { + let (ri, rj) = (root(&group, i), root(&group, j)); + if let Some(g) = group.get_mut(ri.max(rj)) { + *g = ri.min(rj); + } + } + } + } + let mut out: Vec<(usize, Chain)> = Vec::new(); + for (i, (part, _)) in parts.into_iter().enumerate() { + let r = root(&group, i); + match out.iter_mut().find(|(k, _)| *k == r) { + Some((_, Chain::Sketch { ids, .. })) => { + if let Chain::Sketch { ids: more, .. } = part { + ids.extend(more); + } + } + Some((_, Chain::Edges { points, .. })) => { + if let Chain::Edges { points: more, .. } = part { + points.extend(more); + } + } + None => out.push((r, part)), + } + } + Ok(out.into_iter().map(|(_, c)| c).collect()) +} + /// The command a surface dialog makes. pub fn commands(app: &SolveApp, k: &Surf, inputs: &[SelInput], extra: &Map) -> Result, String> { let need = |i: usize, what: &str| if items(inputs, i).is_empty() { Err(format!("select {what} first")) } else { Ok(()) }; @@ -407,6 +523,48 @@ pub fn commands(app: &SolveApp, k: &Surf, inputs: &[SelInput], extra: &Map>()?; ("surface.loft", json!({ "sections": sections })) } + Surf::Boundary { continuity, join } => { + need(0, "the curves (sketch curves or edges)")?; + let word = PATCH_CONTINUITY.get(*continuity).map_or("g0", |(_, w)| *w); + let across = items(inputs, 1); + let picks = items(inputs, 0); + let sketch_ids: Vec = + picks.iter().filter_map(|x| if let Sel::SketchCurve { id } = x { Some(id.clone()) } else { None }).collect(); + // Two to four sketch curves of one sketch and nothing across: the loop patch. + let looped = across.is_empty() && !*join && sketch_ids.len() == picks.len() && (2..=4).contains(&picks.len()); + match looped.then(|| crate::dialogs::curves_sketch(app, &sketch_ids, None)).flatten() { + Some(sk) => { + let mut p = json!({"sketch": sk, "curves": sketch_ids}); + // Orders that differ from side to side (kept from the feature) stand as they are. + if !extra.contains_key("continuity") { + p["continuity"] = json!(word); + } + ("surface.boundary", p) + } + None => { + let curve = |c: &Chain| match c { + Chain::Sketch { sketch, ids } => json!({"sketch": sketch, "curves": ids}), + Chain::Edges { body, points } => { + let mut v = json!({"body": body, "edges": points.iter().map(|p| pt(*p)).collect::>()}); + if *continuity > 0 { + v["continuity"] = json!(word); + } + v + } + }; + let u: Vec = chains(app, picks)?.iter().map(curve).collect(); + let v: Vec = chains(app, across)?.iter().map(curve).collect(); + let mut p = json!({"u": u}); + if !v.is_empty() { + p["v"] = json!(v); + } + if *join { + p["operation"] = json!("join"); + } + ("surface.boundary", p) + } + } + } Surf::Patch { continuity } => { need(0, "a boundary (sketch profiles or a loop of edges)")?; let profiles: Vec<(u64, usize)> = @@ -629,6 +787,46 @@ pub fn for_feature(app: &SolveApp, kind: &FeatureKind, frame: &Mat) -> Option<(K ], ) } + FeatureKind::SurfaceBoundary { curves, orders, crown, samples, .. } => { + let first = orders.first().copied().unwrap_or(solidcraft_engine::doc::kernel::Continuity::G1); + let index = PATCH_CONTINUITY.iter().position(|(_, w)| *w == first.as_str()); + // What the dialog does not show is kept as it was: a crown, the stations, and orders + // that differ from side to side or that the list does not offer (G3). + let mut extra = Map::new(); + if *crown != 0.0 { + extra.insert("crown".into(), json!(crown)); + } + if *samples != solidcraft_engine::doc::BOUNDARY_SAMPLES { + extra.insert("samples".into(), json!(samples)); + } + if index.is_none() || orders.iter().any(|o| *o != first) { + extra.insert("continuity".into(), json!(orders.iter().map(|o| o.as_str()).collect::>())); + } + let k = Surf::Boundary { continuity: index.unwrap_or(1), join: false }; + let mut ins = inputs(&k); + if let Some(inp) = ins.first_mut() { + inp.items = curves.iter().map(|id| Sel::SketchCurve { id: id.clone() }).collect(); + } + return Some((Kind::Surface(k), ins, extra)); + } + FeatureKind::SurfaceNetwork { u, v, operation } => { + // A network with a point for a side has nothing to pick it with here. + if u.iter().chain(v).any(|c| c.point.is_some()) { + return None; + } + let sels = |list: &[solidcraft_engine::doc::NetworkCurve]| -> Vec { + list.iter() + .flat_map(|c| { + let mut out: Vec = c.curves.iter().map(|id| Sel::SketchCurve { id: id.clone() }).collect(); + out.extend(c.edges.iter().filter_map(edge)); + out + }) + .collect() + }; + let g = u.iter().chain(v).find(|c| c.body.is_some()).map_or(solidcraft_engine::doc::kernel::Continuity::G0, |c| c.continuity); + let index = PATCH_CONTINUITY.iter().position(|(_, w)| *w == g.as_str()).unwrap_or(0); + (Surf::Boundary { continuity: index, join: *operation == solidcraft_engine::doc::Operation::Join }, vec![sels(u), sels(v)]) + } FeatureKind::SurfaceLoft { sections, .. } => { (Surf::Loft, vec![sections.iter().flat_map(|sec| sheet_sels(app, sec.sketch.0, sec.profiles.as_ref(), &sec.curves)).collect()]) } diff --git a/crates/ui-egui/src/lib.rs b/crates/ui-egui/src/lib.rs index 638583fa..2b46b969 100644 --- a/crates/ui-egui/src/lib.rs +++ b/crates/ui-egui/src/lib.rs @@ -106,6 +106,8 @@ mod history_off_thread_tests; pub use file_picker::FilePickerFuture; #[cfg(test)] mod blend_arrow_tests; +#[cfg(test)] +mod boundary_dialog_tests; pub(crate) mod display_lod; #[cfg(test)] mod file_load_ui_tests; diff --git a/crates/ui-egui/src/timeline.rs b/crates/ui-egui/src/timeline.rs index 3904b026..a29e7dce 100644 --- a/crates/ui-egui/src/timeline.rs +++ b/crates/ui-egui/src/timeline.rs @@ -25,7 +25,7 @@ pub(crate) fn icon_of(k: &FeatureKind) -> &'static str { FeatureKind::SurfaceRuled { .. } => "surface_ruled", FeatureKind::SurfaceUntrim { .. } => "surface_untrim", FeatureKind::SurfaceBlend { .. } => "surface_blend", - FeatureKind::SurfaceBoundary { .. } => "surface_boundary", + FeatureKind::SurfaceBoundary { .. } | FeatureKind::SurfaceNetwork { .. } => "surface_boundary", FeatureKind::SurfaceMatch { .. } => "continuity", FeatureKind::SurfaceFair { .. } => "surface_fair", FeatureKind::Revolve { .. } => "revolve", @@ -549,6 +549,7 @@ mod tests { ("surface.reverse_normal", json!({"type":"reverse_normal", "body":"A"})), ("surface.blend", json!({"type":"surface_blend", "from":{"body":"A", "edges":[]}, "to":{"body":"B", "edges":[]}})), ("surface.boundary", json!({"type":"surface_boundary", "sketch":0, "curves":[]})), + ("surface.boundary", json!({"type":"surface_network", "u":[]})), ("surface.match", json!({"type":"surface_match", "body":"A", "target":"B", "edges":[]})), ("surface.fair", json!({"type":"surface_fair", "body":"A", "tolerance":"0.001"})), ("solid.rib", json!({"type":"rib", "sketch":0, "curves":[], "thickness":"2"})), diff --git a/docs/oracle.md b/docs/oracle.md index 88241c35..cf154d37 100644 --- a/docs/oracle.md +++ b/docs/oracle.md @@ -9,7 +9,7 @@ Where the expected values come from (the **Source** column): No expected value comes from a commercial CAD product. -**82 / 82 cases pass**; sources: 61 analytic, 21 freecad. 0 known failure(s) tracked in issues; 0 case(s) without expected values yet. +**83 / 83 cases pass**; sources: 62 analytic, 21 freecad. 0 known failure(s) tracked in issues; 0 case(s) without expected values yet. | Case | Source | Result | Surface distance mm (max / mean) | Method | Notes | |---|---|---|---|---|---| @@ -95,6 +95,7 @@ No expected value comes from a commercial CAD product. | 80-fillet-taper-to-fit | analytic | pass | | plate (top 20 -> 2 wide over x 0..40, 15 high) with its top front edge filleted at 3, fit: taper. The radius is the taper law: a clamped cubic B-spline in u = x/40 on 8 equal spans, control radii min(3, 0.9 room) at the Greville abscissae, room = min(20 - 18u, 15) tan 45 deg (top face width, front face height): 3.0000, 3.0000, 3.0000, 3.0000, 3.0000, 3.0000, 3.0000, 3.0000, 3.0000, 2.4750, 1.8000. Sections are quarter rounds; V, moments and areas integrated along the edge | | | 81-fillet-spline-radius | analytic | pass | | box 40x30x20 with its top front edge (y=0, z=20) rounded at a radius that is a clamped cubic B-spline in u = x/40, knots 0,0,0,0,0.25,0.6,1,1,1,1, control radii 2, 5, 1.5, 4, 1, 3 (r(0) = 2, r(40) = 3). Sections are quarter rounds; V, moments and areas integrated along the edge, exactly per knot span | | | 82-thin-extrude-frame | analytic | pass | | square frame: the lines' 40 mm square walled 3 mm inside and 1 mm outside (mitred corners) is 42 x 42 less 34 x 34, 10 mm tall: V = (42^2 - 34^2) * 10, A = 2 (42^2 - 34^2) + 4 * 42 * 10 + 4 * 34 * 10 | | +| 83-boundary-loop-thicken | analytic | pass | | region = 10 x 10 square + half disc r 5 on its top side: A = 100 + 12.5 pi, centroid by moments (half disc's at 4r / (3 pi) above its diameter); slab of thickness 2 over it (symmetric about z = 0): V = 2A, area = 2A + 2 (30 + 5 pi) | | ## Not covered by a reference case diff --git a/tests/oracle/83-boundary-loop-thicken/expected.json b/tests/oracle/83-boundary-loop-thicken/expected.json new file mode 100644 index 00000000..a2b050e3 --- /dev/null +++ b/tests/oracle/83-boundary-loop-thicken/expected.json @@ -0,0 +1,25 @@ +{ + "source": "analytic", + "method": "region = 10 x 10 square + half disc r 5 on its top side: A = 100 + 12.5 pi, centroid by moments (half disc's at 4r / (3 pi) above its diameter); slab of thickness 2 over it (symmetric about z = 0): V = 2A, area = 2A + 2 (30 + 5 pi)", + "body_count": 1, + "volume_mm3": 278.53981633974485, + "area_mm2": 369.9557428756428, + "center_of_mass": [ + 4.999999999999999, + 7.008207500514442, + 0.0 + ], + "bbox": { + "min": [ + 0.0, + 0.0, + -1.0 + ], + "max": [ + 10.0, + 15.0, + 1.0 + ] + }, + "faces": 6 +} diff --git a/tests/oracle/83-boundary-loop-thicken/recipe.json b/tests/oracle/83-boundary-loop-thicken/recipe.json new file mode 100644 index 00000000..78276e06 --- /dev/null +++ b/tests/oracle/83-boundary-loop-thicken/recipe.json @@ -0,0 +1,157 @@ +{ + "id": "83-boundary-loop-thicken", + "name": "Boundary surface through a loop, thickened", + "description": "Four sketch curves on XY (a 10x10 square's bottom, right and left sides, and a half circle of radius 5 bulging up for its top) make one boundary surface; thickened 2 mm symmetrically, the surface then removed.", + "exercises": [ + "boundary_surface", + "thicken" + ], + "units": { + "length": "mm", + "angle": "deg", + "area": "mm2", + "volume": "mm3" + }, + "conventions": { + "world": "right-handed; Z up. XY/XZ/YZ are the origin planes.", + "sketch_frame": "sketch-local 2D coords; world = plane.origin + x*plane.x_dir + y*plane.y_dir (origin in mm). Profile 'loops' list sketch entity ids.", + "point_refs": "'.start' / '.end' / '.center' reuse (share) that existing sketch point; 'at' is its position when referenced.", + "extrude_direction": "positive = along sketch plane normal unless stated", + "selections": "selected faces/edges carry 'point': a world-space point on the face / the edge's parametric midpoint; pick the nearest face/edge to it." + }, + "parameters": [], + "features": [ + { + "op": "sketch", + "name": "Loop", + "plane": { + "ref": "XY", + "origin": [ + 0.0, + 0.0, + 0.0 + ], + "x_dir": [ + 1.0, + 0.0, + 0.0 + ], + "y_dir": [ + 0.0, + 1.0, + 0.0 + ], + "normal": [ + 0.0, + 0.0, + 1.0 + ] + }, + "entities": [ + { + "type": "line", + "start": { + "at": [ + 0, + 0 + ] + }, + "end": { + "at": [ + 10, + 0 + ] + }, + "id": "s1" + }, + { + "type": "line", + "start": { + "ref": "s1.end", + "at": [ + 10, + 0 + ] + }, + "end": { + "at": [ + 10, + 10 + ] + }, + "id": "s2" + }, + { + "type": "arc_three_point", + "start": { + "ref": "s2.end", + "at": [ + 10, + 10 + ] + }, + "end": { + "at": [ + 0, + 10 + ] + }, + "mid": [ + 5, + 15 + ], + "id": "a3" + }, + { + "type": "line", + "start": { + "ref": "a3.end", + "at": [ + 0, + 10 + ] + }, + "end": { + "ref": "s1.start", + "at": [ + 0, + 0 + ] + }, + "id": "s4" + } + ] + }, + { + "op": "boundary_surface", + "name": "Patch", + "sketch": "Loop", + "curves": [ + "s1", + "s2", + "a3", + "s4" + ], + "result_bodies": [ + "Sheet" + ] + }, + { + "op": "thicken", + "name": "Slab", + "bodies": [ + "Sheet" + ], + "thickness": 2, + "symmetric": true, + "operation": "new" + }, + { + "op": "remove", + "name": "Clear", + "bodies": [ + "Sheet" + ] + } + ] +} diff --git a/tools/oracle-analytic/README.md b/tools/oracle-analytic/README.md index b4d1526a..4a6789a2 100644 --- a/tools/oracle-analytic/README.md +++ b/tools/oracle-analytic/README.md @@ -12,6 +12,7 @@ python3 tools/oracle-analytic/cases_basic.py # 01-29, 82: boxes, holes, r python3 tools/oracle-analytic/cases_features.py # 32-66: fillets, shells, sweeps, patterns python3 tools/oracle-analytic/cases_sheet_metal.py # 67-76: sheet metal python3 tools/oracle-analytic/cases_fillets.py # 77-81: variable, chord-length, tapered and spline fillets +python3 tools/oracle-analytic/cases_surfaces.py # 83: boundary surfaces ``` Where a value depends on a SolidCraft convention rather than on geometry alone (a shell's sharp diff --git a/tools/oracle-analytic/cases_surfaces.py b/tools/oracle-analytic/cases_surfaces.py new file mode 100644 index 00000000..cb84dc73 --- /dev/null +++ b/tools/oracle-analytic/cases_surfaces.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python3 +"""Closed-form expected values for the surface oracle cases. + +Run from the worktree root: python3 tools/oracle-analytic/cases_surfaces.py [case ...] +Writes tests/oracle//expected.json with "source": "analytic". + +Every value is derived from first principles (areas and centroids of plane regions bounded by +lines and circular arcs, a slab of constant thickness over such a region). Nothing here is +measured from a CAD kernel. Stdlib only. +""" +import json +import math +import os +import sys + +PI = math.pi +ROOT = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "..", "tests", "oracle") + + +def c83_boundary_loop_thicken(): + """A boundary surface through four sketch curves on XY (the 10 x 10 square's bottom, right and + left sides and, for its top, a half circle of radius 5 bulging out to y = 15) is the flat + region they bound; thickened 2 mm symmetrically it is a slab from z = -1 to z = 1.""" + side, r, t = 10.0, 5.0, 2.0 + sq_a, sq_c = side * side, (5.0, 5.0) + half_a, half_c = PI * r * r / 2, (5.0, side + 4 * r / (3 * PI)) + area = sq_a + half_a + cx = (sq_a * sq_c[0] + half_a * half_c[0]) / area + cy = (sq_a * sq_c[1] + half_a * half_c[1]) / area + perimeter = 3 * side + PI * r + return { + "source": "analytic", + "method": "region = 10 x 10 square + half disc r 5 on its top side: A = 100 + 12.5 pi, centroid by moments " + "(half disc's at 4r / (3 pi) above its diameter); slab of thickness 2 over it (symmetric about z = 0): " + "V = 2A, area = 2A + 2 (30 + 5 pi)", + "body_count": 1, + "volume_mm3": t * area, + "area_mm2": 2 * area + t * perimeter, + "center_of_mass": [cx, cy, 0.0], + "bbox": {"min": [0.0, 0.0, -t / 2], "max": [side, side + r, t / 2]}, + "faces": 6, + } + + +CASES = { + "83-boundary-loop-thicken": c83_boundary_loop_thicken, +} + + +def main(argv): + names = argv[1:] or list(CASES) + for n in names: + d = CASES[n]() + path = os.path.join(ROOT, n, "expected.json") + with open(path, "w") as f: + json.dump(d, f, indent=2) + f.write("\n") + print(f"{n}: V={d['volume_mm3']!r} A={d['area_mm2']!r} com={d['center_of_mass']}") + + +if __name__ == "__main__": + main(sys.argv) diff --git a/vendor/README.md b/vendor/README.md index 41521026..bc864466 100644 --- a/vendor/README.md +++ b/vendor/README.md @@ -59,6 +59,10 @@ Changes: back to the seam at the other end of the period and fold its parameter polygon over itself (no intersection loop then lay inside the face). `crates/kernel/src/seam_faces_tests.rs` holds it to a cone read from STEP. +- `nurbs/bspsurface.rs`: `uuder` and `vvder` multiply by p·(p − 1), as `BSplineCurve::der2` + does. They multiplied by p only, so the second derivatives of every B-spline and NURBS surface + of degree p ≥ 3 were 1/(p − 1) of their value (a cubic surface's curvature read as half, #763). + `crates/kernel/src/surface_second_derivative_tests.rs` holds them to finite differences. ## truck-stepio 0.3.0 (Apache-2.0, © RICOS Co. Ltd., https://github.com/ricosjp/truck) diff --git a/vendor/truck-geometry/src/nurbs/bspsurface.rs b/vendor/truck-geometry/src/nurbs/bspsurface.rs index fcedac20..4294d855 100644 --- a/vendor/truck-geometry/src/nurbs/bspsurface.rs +++ b/vendor/truck-geometry/src/nurbs/bspsurface.rs @@ -729,11 +729,15 @@ impl> ParametricSurface for BSplineSurface

{ }; basis1.iter().enumerate().fold(sum, closure) }; + // SolidCraft: the second derivative of a degree-p spline is p·(p − 1) times the sum of the + // second differences against the degree p − 2 basis (as `BSplineCurve::der2` has it). The + // factor p − 1 was missing, which halved the second derivative of a cubic surface. basis0 .into_iter() .enumerate() .fold(P::Diff::zero(), closure) * degree0 as f64 + * (degree0 - 1) as f64 } /// Substitutes 2nd-ord derived B-spline surface by the second parameter `v`. @@ -778,11 +782,13 @@ impl> ParametricSurface for BSplineSurface

{ }; basis0.iter().enumerate().fold(sum, closure) }; + // SolidCraft: times p·(p − 1), as in `uuder`. basis1 .into_iter() .enumerate() .fold(P::Diff::zero(), closure) * degree1 as f64 + * (degree1 - 1) as f64 } /// Substitutes 2nd-ord derived B-spline surface by the both parameters `u, v`.