From 83c31d77c4c223aec7e7900c6528db94262c6e54 Mon Sep 17 00:00:00 2001 From: Irakli Kourtis Date: Sat, 17 Oct 2026 21:47:38 -0400 Subject: [PATCH] surface.boundary: surfaces through boundary curves and curve networks (#763) (#2129) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * feat(geom): surfaces through curve networks, and spline spaces to build them Adds spline_space (clamped B-spline spaces, Hermite interpolation, exact embedding by Bézier degree elevation and knot insertion, a monotone cubic for reparameterising) and network (loft, Coons and Gordon surfaces through curves in one or two directions, with G1/G2 cross-boundary conditions against a neighbouring face). Part of #763. Co-Authored-By: Claude Opus 5.5 * feat(kernel): boundary surfaces from sketch curves, 3D curves and edges kernel::boundary_network turns sketch chains, 3D sketch polylines and edge chains into network curves, gives G1 frames of the faces along open edges (EdgeFaces), builds the B-spline face (sides split at chain vertices so a join shares the body's edges) and joins it into a surface body by stitching. G2 is refused for now (not supported yet). Part of #763. Co-Authored-By: Claude Opus 5.5 * fix(kernel): second derivatives of B-spline surfaces in the vendored truck-geometry BSplineSurface::uuder and vvder multiplied the second differences by p instead of p·(p − 1), as BSplineCurve::der2 does, so every surface of degree p ≥ 3 (NURBS too, through its homogeneous surface) reported 1/(p − 1) of its second derivative: a cubic surface's curvature read as half. This is why a G2 boundary surface on a cylinder measured half the curvature it carries (#763). New tests hold both to finite differences for degrees 2 to 5. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01LGeCa5jWWdroMWvXxN3LLM * feat(kernel): G2 boundary surfaces against a body's face The network always carried the face's curvature across a G2 edge; the refusal came from measuring it with truck's halved surface second derivatives (fixed in the previous commit). EdgeFaces now reads the face's second fundamental form from the exact second derivatives instead of differencing the normal. A new test holds the curvature across a cylinder's edge to the cylinder's within 0.1%, also with a sheared far curve; G3 is refused (a network holds up to G2). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01LGeCa5jWWdroMWvXxN3LLM * test(kernel): a boundary patch with straight corners measures its own area A stadium's arcs meet its lines at 180°, so the Coons patch is singular at its corners. From such a corner the reference mesher's boundary walk ran onto the B-spline's extrapolation past its domain (a planar spline passes through every point of its plane there), and the face measured 289 mm² for its 158.5 (#763). The vendored mesher's bounded projection (nearest_projection, #2043) keeps the walk on the domain; this holds it there, by the faces, the measuring mesh and the display mesh. The display path (#1653) was never affected. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01LGeCa5jWWdroMWvXxN3LLM * feat(geom,kernel): three-sided boundary surfaces (one side collapsed to a point) A network curve may be a PointCurve: the apex where the two sides beside it meet. The curves through it cross it at the point itself, its row of the net collapses there, and the kernel's face leaves that side out of its loop (three edges, one vertex at the apex). One point per patch, as a side; refused elsewhere, when it moves, or with a continuity condition. Tests: a triangle and a quarter disc filled flat with their analytic areas (40 mm², 25π), by geom and by the kernel's faces and meshes. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01LGeCa5jWWdroMWvXxN3LLM * feat(doc): boundary surfaces through curve networks and loops of three or four FeatureKind::SurfaceNetwork saves a network as `u` and `v` curves, each a sketch chain (planar or 3D), a chain of a body's edges with an optional G1/G2 continuity, or a point (the apex of a three-sided patch), and an `operation` that can join the surface into the body of its edges. A SurfaceBoundary of three or four sketch curves, refused until now, builds as the flat Coons patch of its loop (three collapse the side opposite the first curve). The format change is one pending file; feature-walk goldens record the new variant's sketches, bodies and points. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01LGeCa5jWWdroMWvXxN3LLM * feat(engine): surface.boundary through loops of three or four and curve networks The command keeps its loop form (two to four curves of one sketch; three and four now build the flat patch of the loop, crown refused there) and gains a network form: `u` and `v`, each curve a sketch curve id, a {sketch?, curves} chain (planar or 3D), a {body, edges, continuity?} edge chain (g1/g2 with the face beside it) or a {point} for a three-sided patch, with `operation` new or join. Edits may move between the two forms (timeline.redefine accepts the command now that it has a dialog). Tests: analytic areas of loop patches and of a hyperbolic paraboloid bounded by four 3D lines (every mesh vertex on z = xy/10), a loft across two sketches, G1 and G2 measured across a cylinder's edge with inspect.continuity (angle < 0.1°; curvature jump 1e-16 with G2 against 0.027 with G1), 26 hostile networks refused without panicking, and u/v added to hostile_params_never_panic. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01LGeCa5jWWdroMWvXxN3LLM * feat(ui): Boundary Surface dialog Curves and Cross curves take sketch curves and edges (picks that meet end to end in one sketch or body are one chain); Continuity and Operation (New Body or Join). Two to four sketch curves of one sketch with nothing across are the loop patch; anything else is a network. Edit Feature brings either back (per-side orders, crown and stations the dialog does not show are kept). Strings in en, ja and zh. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01LGeCa5jWWdroMWvXxN3LLM * docs: ROADMAP M8 and the bridge patch's notes for boundary networks (#763) Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01LGeCa5jWWdroMWvXxN3LLM * test(oracle): a boundary surface through a loop, thickened (#763) Case 83: four sketch curves (a square's three sides and a half circle for its top) make one boundary surface, thickened 2 mm; analytic values from tools/oracle-analytic/cases_surfaces.py (V = 2A, area 2A + 2P, centroid by moments). Recipes gain boundary_surface, thicken and remove. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01LGeCa5jWWdroMWvXxN3LLM * refactor: boundary surfaces, review cleanups Project onto the face once per derivative evaluation (the frame is passed to natural() and reused by scale()), keep each track's finite-difference step instead of rescanning its table, export the kernel's edge-pick and chain caps and the default station count for the engine, doc and dialog to share, map TooLarge spline errors to KernelError::TooLarge, reuse the engine's operation parser, and skip the dialog's sketch lookup when the picks are no loop. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01LGeCa5jWWdroMWvXxN3LLM --------- Co-authored-by: jyviko <2501648+jyviko@users.noreply.github.com> Co-authored-by: Claude Opus 5.5 Co-authored-by: bherbruck --- ATTRIBUTION.md | 2 +- ROADMAP.md | 2 +- apps/solidcraft-cli/src/recipe.rs | 18 + assets/locales/en.json | 2 + assets/locales/ja.json | 2 + assets/locales/zh.json | 2 + crates/doc/src/document.rs | 43 +- crates/doc/src/eval.rs | 15 + crates/doc/src/eval_surface.rs | 170 +- crates/doc/src/feature_walk_tests/fixtures.rs | 4 + .../src/feature_walk_tests/golden_exprs.txt | 1 + .../feature_walk_tests/golden_geometry.txt | 2 + .../src/feature_walk_tests/golden_refs.txt | 4 + .../pending/boundary-surface-networks.rs | 7 + crates/engine/src/cmd/edit.rs | 6 +- crates/engine/src/cmd/surface_boundary.rs | 190 +- .../engine/src/cmd/surface_boundary_tests.rs | 267 +++ crates/engine/src/tests.rs | 2 + crates/geom/src/lib.rs | 7 +- crates/geom/src/network.rs | 1623 +++++++++++++++++ crates/geom/src/spline_space.rs | 584 ++++++ crates/kernel/src/boundary.rs | 13 +- crates/kernel/src/boundary_network.rs | 744 ++++++++ crates/kernel/src/lib.rs | 4 + .../src/surface_second_derivative_tests.rs | 96 + crates/kernel/src/surfaces.rs | 2 +- crates/ui-egui/src/boundary_dialog_tests.rs | 127 ++ crates/ui-egui/src/dialogs_surface.rs | 200 +- crates/ui-egui/src/lib.rs | 2 + crates/ui-egui/src/timeline.rs | 3 +- docs/oracle.md | 3 +- .../83-boundary-loop-thicken/expected.json | 25 + .../83-boundary-loop-thicken/recipe.json | 157 ++ tools/oracle-analytic/README.md | 1 + tools/oracle-analytic/cases_surfaces.py | 62 + vendor/README.md | 4 + vendor/truck-geometry/src/nurbs/bspsurface.rs | 6 + 37 files changed, 4358 insertions(+), 44 deletions(-) create mode 100644 crates/doc/src/format/pending/boundary-surface-networks.rs create mode 100644 crates/engine/src/cmd/surface_boundary_tests.rs create mode 100644 crates/geom/src/network.rs create mode 100644 crates/geom/src/spline_space.rs create mode 100644 crates/kernel/src/boundary_network.rs create mode 100644 crates/kernel/src/surface_second_derivative_tests.rs create mode 100644 crates/ui-egui/src/boundary_dialog_tests.rs create mode 100644 tests/oracle/83-boundary-loop-thicken/expected.json create mode 100644 tests/oracle/83-boundary-loop-thicken/recipe.json create mode 100644 tools/oracle-analytic/cases_surfaces.py 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`.