From fac3f96bdd82a146c952d8a35d7e50bb85aa995c Mon Sep 17 00:00:00 2001 From: Irakli Kourtis Date: Thu, 15 Oct 2026 23:07:24 -0400 Subject: [PATCH] Patch: G1/G2 to the faces beside the loop, general and N-sided loops, points (#759) (#1521) Thanks @jyviko. --- ROADMAP.md | 2 + assets/locales/de.json | 6 + assets/locales/en.json | 6 + assets/locales/fr.json | 6 + assets/locales/ja.json | 6 + assets/locales/ko.json | 6 + assets/locales/ru.json | 6 + assets/locales/uk.json | 6 + assets/locales/zh.json | 6 + crates/doc/src/document.rs | 47 + crates/doc/src/eval_surface.rs | 10 +- .../pending/patch-continuity-and-points.rs | 6 + crates/engine/src/cmd/surface.rs | 142 +- crates/engine/src/tests.rs | 2 + crates/geom/src/fill.rs | 132 +- crates/geom/src/fill_general.rs | 1150 +++++++++++++++++ crates/geom/src/fill_general_tests.rs | 191 +++ crates/geom/src/lib.rs | 1 + crates/kernel/src/bridge.rs | 122 +- crates/kernel/src/fill.rs | 286 +++- crates/kernel/src/fill_general_tests.rs | 506 ++++++++ crates/kernel/src/lib.rs | 4 +- crates/ui-egui/src/dialogs_surface.rs | 65 +- crates/ui-egui/src/preselect_tests.rs | 51 + 24 files changed, 2688 insertions(+), 77 deletions(-) create mode 100644 crates/doc/src/format/pending/patch-continuity-and-points.rs create mode 100644 crates/geom/src/fill_general.rs create mode 100644 crates/geom/src/fill_general_tests.rs create mode 100644 crates/kernel/src/fill_general_tests.rs diff --git a/ROADMAP.md b/ROADMAP.md index 82459115..7cd14739 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -53,6 +53,8 @@ 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 of any shape and number of sides as one face, G0/G1/G2 to the faces beside each edge, through points, #759), 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. Proposed, in order: surface Extrude/Revolve/Sweep/Loft from open and closed curves (#757), Offset surface (#758), Unstitch and Reverse Normal (#766), Delete Face without healing (#768); Extend curved faces and edges (#760), mutual trim and trim by projected curves (#767), Untrim (#761), Thicken every surface kind (#765); Ruled (#762), Boundary surface and curve networks with G0/G1/G2 (#763), Blend surface (#764), sweep guide rails and two-rail sweeps (#771; orientation, twist and end scale done; loft rails #201), fillet between surface bodies (#769); continuity analysis, surface curvature combs and isocurves (#770) | 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), 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), 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 and Interference (overlaps shown red); joints in the timeline and the Browser | in progress | 32 | | M8 | Surface workspace (#772). Done: Patch (sketch regions; edge loops filled G0), 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. Proposed, in order: surface Extrude/Revolve/Sweep/Loft from open and closed curves (#757), Offset surface (#758), Unstitch and Reverse Normal (#766), Delete Face without healing (#768); Extend curved faces and edges (#760), mutual trim and trim by projected curves (#767), Untrim (#761), Thicken every surface kind (#765); Ruled (#762), Boundary surface and curve networks with G0/G1/G2 (#763), Patch with G1/G2 and general loops (#759), Blend surface (#764), sweep guide rails and two-rail sweeps (#771; orientation, twist and end scale done; loft rails #201), fillet between surface bodies (#769); continuity analysis, surface curvature combs and isocurves (#770) | 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), gears, rack and pinion and screws (motion links whose ratio follows picked radii, parameters or a pitch, also inside closed loops), 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), 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, Gear, Rack and Pinion, Screw 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 | diff --git a/assets/locales/de.json b/assets/locales/de.json index bf66c50d..e87a5cd3 100644 --- a/assets/locales/de.json +++ b/assets/locales/de.json @@ -8,14 +8,17 @@ "Area": "Fläche", "BROWSER": "BROWSER", "Bends": "Biegungen", + "Boundary": "Begrenzung", "Centre to centre": "Mitte zu Mitte", "Circumference": "Umfang", "Construction": "Konstruktionslinie", + "Continuity": "Stetigkeit", "Copy Error Message": "Fehlermeldung kopieren", "Copy Warning Message": "Warnmeldung kopieren", "Corner Relief": "Eckenentlastung", "Current orthographic view (1:1 mm)": "Aktuelle orthogonale Ansicht (1:1 mm)", "Current perspective view": "Aktuelle perspektivische Ansicht", + "Curvature (G2)": "Krümmung (G2)", "Do not export hidden edges": "Verdeckte Kanten nicht exportieren", "Drive Joint": "Gelenk antreiben", "Edges": "Kanten", @@ -42,7 +45,9 @@ "Notch": "Ausklinkung", "Outside": "Außen", "Overall length": "Gesamtlänge", + "Pass through": "Durch Punkte", "Picked": "Ausgewählte", + "Position (G0)": "Position (G0)", "Profiles/Curves": "Profile/Kurven", "Projection": "Projektion", "RULED SURFACE": "REGELFLÄCHE", @@ -62,6 +67,7 @@ "Sketch geometry and profiles": "Skizzengeometrie und Profile", "Sketch points": "Skizzenpunkte", "Slot width": "Langlochbreite", + "Tangent (G1)": "Tangente (G1)", "Tear": "Reißen", "Thickness Outside": "Dicke außen", "Thin": "Dünnwandig", diff --git a/assets/locales/en.json b/assets/locales/en.json index c4f790a1..1e00199b 100644 --- a/assets/locales/en.json +++ b/assets/locales/en.json @@ -34,6 +34,7 @@ "Bend Radius": "Bend Radius", "Bends": "Bends", "Body": "Body", + "Boundary": "Boundary", "Bounding box": "Bounding box", "CSV or JSON: updates existing parameters, adds new ones": "CSV or JSON: updates existing parameters, adds new ones", "Capture All Positions": "Capture All Positions", @@ -56,6 +57,7 @@ "Configuration rules": "Configuration rules", "Construction": "Construction", "Contact": "Contact", + "Continuity": "Continuity", "Copy Error Message": "Copy Error Message", "Copy Warning Message": "Copy Warning Message", "Copy all licences": "Copy all licences", @@ -70,6 +72,7 @@ "Cube": "Cube", "Current orthographic view (1:1 mm)": "Current orthographic view (1:1 mm)", "Current perspective view": "Current perspective view", + "Curvature (G2)": "Curvature (G2)", "Customize Toolbar": "Customize Toolbar", "Delete all": "Delete all", "Delete rule": "Delete rule", @@ -159,6 +162,7 @@ "Outside": "Outside", "Overall length": "Overall length", "Parameter values": "Parameter values", + "Pass through": "Pass through", "Path": "Path", "Per Edge": "Per Edge", "Pick the tool bodies": "Pick the tool bodies", @@ -167,6 +171,7 @@ "Playback": "Playback", "Position": "Position", "Position (0 to 1)": "Position (0 to 1)", + "Position (G0)": "Position (G0)", "Position along the edge, 0 at the start end, 1 at the other": "Position along the edge, 0 at the start end, 1 at the other", "Preview Row": "Preview Row", "Preview leaves the current design unchanged.": "Preview leaves the current design unchanged.", @@ -239,6 +244,7 @@ "Style": "Style", "Suppressed": "Suppressed", "Swap Ends": "Swap Ends", + "Tangent (G1)": "Tangent (G1)", "Tangent Chain": "Tangent Chain", "Tear": "Tear", "The dashed sketch is a saved plane guide. Select a new face to repair the attachment.": "The dashed sketch is a saved plane guide. Select a new face to repair the attachment.", diff --git a/assets/locales/fr.json b/assets/locales/fr.json index df64badb..bff40076 100644 --- a/assets/locales/fr.json +++ b/assets/locales/fr.json @@ -32,6 +32,7 @@ "Bend Radius": "Rayon de pliage", "Bends": "Plis", "Body": "Corps", + "Boundary": "Contour", "Bounding box": "Boîte englobante", "CSV or JSON: updates existing parameters, adds new ones": "CSV ou JSON : met à jour les paramètres existants, ajoute les nouveaux", "Capture All Positions": "Capturer toutes les positions", @@ -52,6 +53,7 @@ "Configuration rules": "Règles de configuration", "Construction": "Construction", "Contact": "Contact", + "Continuity": "Continuité", "Copy Error Message": "Copier le message d'erreur", "Copy Warning Message": "Copier le message d'avertissement", "Copy all licences": "Copier toutes les licences", @@ -64,6 +66,7 @@ "Countersink angle": "Angle de la fraisure", "Create Issue — review draft": "Créer un ticket — relire le brouillon", "Cube": "Cube", + "Curvature (G2)": "Courbure (G2)", "Customize Toolbar": "Personnaliser la barre d'outils", "Delete all": "Tout supprimer", "Delete rule": "Supprimer la règle", @@ -143,6 +146,7 @@ "Other side": "Autre côté", "Outside": "Extérieur", "Parameter values": "Valeurs des paramètres", + "Pass through": "Passer par", "Path": "Trajectoire", "Per Edge": "Par arête", "Pick the tool bodies": "Choisir les corps outils", @@ -151,6 +155,7 @@ "Playback": "Lecture", "Position": "Position", "Position (0 to 1)": "Position (0 à 1)", + "Position (G0)": "Position (G0)", "Position along the edge, 0 at the start end, 1 at the other": "Position le long de l'arête : 0 à l'extrémité de départ, 1 à l'autre", "Preview Row": "Aperçu de la ligne", "Preview leaves the current design unchanged.": "L'aperçu ne modifie pas la conception actuelle.", @@ -221,6 +226,7 @@ "Style": "Style", "Suppressed": "Désactivé", "Swap Ends": "Permuter les extrémités", + "Tangent (G1)": "Tangence (G1)", "Tangent Chain": "Chaîne tangente", "Tear": "Déchirure", "The dashed sketch is a saved plane guide. Select a new face to repair the attachment.": "L'esquisse en pointillés est un guide de plan enregistré. Sélectionnez une nouvelle face pour réparer le rattachement.", diff --git a/assets/locales/ja.json b/assets/locales/ja.json index b1ef7be0..df92e435 100644 --- a/assets/locales/ja.json +++ b/assets/locales/ja.json @@ -34,6 +34,7 @@ "Bend Radius": "曲げ半径", "Bends": "曲げ", "Body": "ボディ", + "Boundary": "境界", "Bounding box": "境界ボックス", "CSV or JSON: updates existing parameters, adds new ones": "CSVまたはJSON：既存のパラメータを更新し、新しいパラメータを追加します", "Capture All Positions": "すべての位置をキャプチャ", @@ -56,6 +57,7 @@ "Configuration rules": "構成ルール", "Construction": "コンストラクション", "Contact": "接触", + "Continuity": "連続性", "Copy Error Message": "エラーメッセージをコピー", "Copy Warning Message": "警告メッセージをコピー", "Copy all licences": "すべてのライセンスをコピー", @@ -70,6 +72,7 @@ "Cube": "キューブ", "Current orthographic view (1:1 mm)": "現在の平行投影（実寸・mm）", "Current perspective view": "現在の透視投影", + "Curvature (G2)": "曲率 (G2)", "Customize Toolbar": "ツールバーのカスタマイズ", "Delete all": "すべて削除", "Delete rule": "ルールを削除", @@ -159,6 +162,7 @@ "Outside": "外側", "Overall length": "全長", "Parameter values": "パラメータ値", + "Pass through": "通過点", "Path": "パス", "Per Edge": "エッジごと", "Pick the tool bodies": "ツールボディを選択", @@ -167,6 +171,7 @@ "Playback": "再生", "Position": "位置", "Position (0 to 1)": "位置 (0〜1)", + "Position (G0)": "位置 (G0)", "Position along the edge, 0 at the start end, 1 at the other": "エッジ上の位置（始点が0、終点が1）", "Preview Row": "行をプレビュー", "Preview leaves the current design unchanged.": "プレビューは現在のデザインを変更しません。", @@ -239,6 +244,7 @@ "Style": "スタイル", "Suppressed": "抑制", "Swap Ends": "端点を入れ替え", + "Tangent (G1)": "接線 (G1)", "Tangent Chain": "正接チェーン", "Tear": "引き裂き", "The dashed sketch is a saved plane guide. Select a new face to repair the attachment.": "破線のスケッチは保存された平面ガイドです。新しい面を選択してアタッチメントを修復してください。", diff --git a/assets/locales/ko.json b/assets/locales/ko.json index 611179b5..1ac034ef 100644 --- a/assets/locales/ko.json +++ b/assets/locales/ko.json @@ -8,14 +8,17 @@ "Area": "면적", "BROWSER": "브라우저", "Bends": "벤딩", + "Boundary": "경계", "Centre to centre": "중심 간", "Circumference": "원주", "Construction": "구성선", + "Continuity": "연속성", "Copy Error Message": "오류 메시지 복사", "Copy Warning Message": "경고 메시지 복사", "Corner Relief": "모서리 릴리프", "Current orthographic view (1:1 mm)": "현재 직교 투영 뷰 (1:1 mm)", "Current perspective view": "현재 원근 투영 뷰", + "Curvature (G2)": "곡률 (G2)", "Do not export hidden edges": "숨겨진 모서리 내보내지 않기", "Drive Joint": "조인트 구동", "Edges": "모서리", @@ -42,7 +45,9 @@ "Notch": "노치", "Outside": "바깥쪽", "Overall length": "전체 길이", + "Pass through": "통과점", "Picked": "선택한 항목", + "Position (G0)": "위치 (G0)", "Profiles/Curves": "프로파일/곡선", "Projection": "투영", "RULED SURFACE": "룰드 서피스", @@ -62,6 +67,7 @@ "Sketch geometry and profiles": "스케치 기하 및 프로파일", "Sketch points": "스케치 점", "Slot width": "슬롯 폭", + "Tangent (G1)": "접선 (G1)", "Tear": "찢김", "Thickness Outside": "바깥쪽 두께", "Thin": "얇은 벽", diff --git a/assets/locales/ru.json b/assets/locales/ru.json index e7d0e2b4..52fcd977 100644 --- a/assets/locales/ru.json +++ b/assets/locales/ru.json @@ -28,6 +28,7 @@ "Bend Radius": "Радиус гиба", "Bends": "Гибы", "Body": "Тело", + "Boundary": "Граница", "Bounding box": "Габаритная рамка", "CSV or JSON: updates existing parameters, adds new ones": "CSV или JSON: обновляет существующие параметры, добавляет новые", "Capture All Positions": "Зафиксировать все положения", @@ -48,6 +49,7 @@ "Configuration rules": "Правила конфигурации", "Construction": "Вспомогательная", "Contact": "Контакт", + "Continuity": "Непрерывность", "Copy Error Message": "Скопировать сообщение об ошибке", "Copy Warning Message": "Скопировать предупреждение", "Copy all licences": "Копировать все лицензии", @@ -60,6 +62,7 @@ "Countersink angle": "Угол зенковки", "Create Issue — review draft": "Создать Issue — проверить черновик", "Cube": "Куб", + "Curvature (G2)": "Кривизна (G2)", "Customize Toolbar": "Настроить панель инструментов", "Delete all": "Удалить все", "Delete rule": "Удалить правило", @@ -131,6 +134,7 @@ "Open-source parametric 3D CAD in Rust. Licensed MIT OR Apache-2.0.": "Открытая параметрическая 3D CAD на Rust. Лицензия MIT OR Apache-2.0.", "Other side": "Другая сторона", "Parameter values": "Значения параметров", + "Pass through": "Через точки", "Path": "Траектория", "Per Edge": "По каждому ребру", "Pick the tool bodies": "Выберите тела-инструменты", @@ -138,6 +142,7 @@ "Playback": "Воспроизведение", "Position": "Положение", "Position (0 to 1)": "Положение (от 0 до 1)", + "Position (G0)": "Положение (G0)", "Position along the edge, 0 at the start end, 1 at the other": "Положение вдоль ребра: 0 — начало, 1 — конец", "Preview Row": "Строка предпросмотра", "Preview leaves the current design unchanged.": "Предпросмотр не меняет текущий дизайн.", @@ -203,6 +208,7 @@ "Style": "Стиль", "Suppressed": "Подавлено", "Swap Ends": "Поменять концы", + "Tangent (G1)": "Касание (G1)", "Tangent Chain": "Касательная цепочка", "Tear": "Разрыв", "The dashed sketch is a saved plane guide. Select a new face to repair the attachment.": "Пунктирный эскиз — сохранённая подсказка плоскости. Выберите новую грань, чтобы восстановить привязку.", diff --git a/assets/locales/uk.json b/assets/locales/uk.json index 15b0e1b4..f5c5d015 100644 --- a/assets/locales/uk.json +++ b/assets/locales/uk.json @@ -32,6 +32,7 @@ "Bend Radius": "Радіус згину", "Bends": "Згини", "Body": "Тіло", + "Boundary": "Межа", "Bounding box": "Обмежувальна рамка", "CSV or JSON: updates existing parameters, adds new ones": "CSV або JSON: оновлює наявні параметри, додає нові", "Capture All Positions": "Зафіксувати всі положення", @@ -52,6 +53,7 @@ "Configuration rules": "Правила конфігурації", "Construction": "Конструкційна", "Contact": "Контакт", + "Continuity": "Неперервність", "Copy all licences": "Копіювати всі ліцензії", "Copy diagnostic information": "Копіювати діагностичну інформацію", "Copy report": "Копіювати звіт", @@ -62,6 +64,7 @@ "Countersink angle": "Кут зенкування", "Create Issue — review draft": "Створити повідомлення про проблему — перегляд чернетки", "Cube": "Куб", + "Curvature (G2)": "Кривина (G2)", "Customize Toolbar": "Налаштувати панель інструментів", "Delete all": "Видалити все", "Delete rule": "Видалити правило", @@ -135,6 +138,7 @@ "Other side": "Інший бік", "Outside": "Назовні", "Parameter values": "Значення параметрів", + "Pass through": "Через точки", "Path": "Траєкторія", "Per Edge": "На кожне ребро", "Pick the tool bodies": "Виберіть тіла інструмента", @@ -142,6 +146,7 @@ "Playback": "Відтворення", "Position": "Положення", "Position (0 to 1)": "Положення (від 0 до 1)", + "Position (G0)": "Положення (G0)", "Position along the edge, 0 at the start end, 1 at the other": "Положення вздовж ребра: 0 на початку, 1 на іншому кінці", "Preview Row": "Рядок перегляду", "Preview leaves the current design unchanged.": "Попередній перегляд не змінює поточний проєкт.", @@ -209,6 +214,7 @@ "Style": "Стиль", "Suppressed": "Придушено", "Swap Ends": "Поміняти кінці місцями", + "Tangent (G1)": "Дотик (G1)", "Tangent Chain": "Дотичний ланцюг", "Tear": "Розрив", "The dashed sketch is a saved plane guide. Select a new face to repair the attachment.": "Штриховий ескіз — це збережений орієнтир площини. Виберіть нову грань, щоб відновити прив'язку.", diff --git a/assets/locales/zh.json b/assets/locales/zh.json index ed6d0c8d..1481b970 100644 --- a/assets/locales/zh.json +++ b/assets/locales/zh.json @@ -34,6 +34,7 @@ "Bend Radius": "折弯半径", "Bends": "折弯", "Body": "实体", + "Boundary": "边界", "Bounding box": "边界框", "CSV or JSON: updates existing parameters, adds new ones": "CSV 或 JSON：更新现有参数、添加新参数", "Capture All Positions": "捕获所有位置", @@ -56,6 +57,7 @@ "Configuration rules": "配置规则", "Construction": "构造", "Contact": "接触", + "Continuity": "连续性", "Copy Error Message": "复制错误消息", "Copy Warning Message": "复制警告消息", "Copy all licences": "复制所有许可证", @@ -70,6 +72,7 @@ "Cube": "立方体", "Current orthographic view (1:1 mm)": "当前正交视图（1:1毫米）", "Current perspective view": "当前透视视图", + "Curvature (G2)": "曲率 (G2)", "Customize Toolbar": "自定义工具栏", "Delete all": "全部删除", "Delete rule": "删除规则", @@ -159,6 +162,7 @@ "Outside": "外侧", "Overall length": "总长", "Parameter values": "参数值", + "Pass through": "通过点", "Path": "路径", "Per Edge": "每条边", "Pick the tool bodies": "选择工具实体", @@ -167,6 +171,7 @@ "Playback": "播放", "Position": "位置", "Position (0 to 1)": "位置（0 到 1）", + "Position (G0)": "位置 (G0)", "Position along the edge, 0 at the start end, 1 at the other": "沿边的位置，起点端为 0，另一端为 1", "Preview Row": "预览行", "Preview leaves the current design unchanged.": "预览不会更改当前设计。", @@ -239,6 +244,7 @@ "Style": "样式", "Suppressed": "已抑制", "Swap Ends": "交换端点", + "Tangent (G1)": "相切 (G1)", "Tangent Chain": "相切链", "Tear": "撕裂", "The dashed sketch is a saved plane guide. Select a new face to repair the attachment.": "虚线草图是已保存的平面引导。选择一个新面以修复附着。", diff --git a/crates/doc/src/document.rs b/crates/doc/src/document.rs index 233c0a2e..3aeb554c 100644 --- a/crates/doc/src/document.rs +++ b/crates/doc/src/document.rs @@ -1276,6 +1276,14 @@ pub enum FeatureKind { /// exact plane whatever this says). #[serde(default, skip_serializing_if = "is_zero")] crown: f64, + /// How a patch on edges meets the faces beside its loop (G0 everywhere by default, and + /// then not saved). + #[serde(default, skip_serializing_if = "PatchContinuity::is_g0")] + continuity: PatchContinuity, + /// Points a patch on edges passes through. + #[serde(default, skip_serializing_if = "Vec::is_empty")] + #[walk(pos)] + points: Vec, }, /// Surface bodies sewn along shared edges into the first one (a solid when they close). Stitch { @@ -1666,6 +1674,45 @@ pub struct BlendSide { pub edges: Vec, } +/// 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)] +pub struct PatchContinuity { + #[serde(default = "patch_g0", skip_serializing_if = "is_g0")] + pub all: Continuity, + #[serde(default, skip_serializing_if = "Vec::is_empty")] + pub edges: Vec, +} + +/// One edge of a patch's loop (a point on it) and how the patch meets the face beside it. +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize, solidcraft_derive::Walk)] +pub struct EdgeContinuity { + #[walk(pos)] + pub edge: Vec3, + pub g: Continuity, +} + +impl Default for PatchContinuity { + fn default() -> PatchContinuity { + PatchContinuity { all: Continuity::G0, edges: Vec::new() } + } +} + +impl PatchContinuity { + /// G0 everywhere: the patch only meets its edges. + pub fn is_g0(&self) -> bool { + self.all == Continuity::G0 && self.edges.iter().all(|e| e.g == Continuity::G0) + } +} + +fn patch_g0() -> Continuity { + Continuity::G0 +} + +fn is_g0(g: &Continuity) -> bool { + *g == Continuity::G0 +} + /// The continuity a blend repeats by default: the tangent plane on both sides. fn blend_orders() -> [Continuity; 2] { [Continuity::G1, Continuity::G1] diff --git a/crates/doc/src/eval_surface.rs b/crates/doc/src/eval_surface.rs index 41d0aa7b..09969eef 100644 --- a/crates/doc/src/eval_surface.rs +++ b/crates/doc/src/eval_surface.rs @@ -8,7 +8,7 @@ use solidcraft_geom::Seg2; pub(super) fn eval(doc: &Document, vals: &BTreeMap, f: &Feature, st: &mut ModelState, warning: &mut Option) -> Result<()> { match &f.kind { - FeatureKind::Patch { sketch, profiles, body, edges, crown } => { + FeatureKind::Patch { sketch, profiles, body, edges, crown, continuity, points } => { let mut made = Vec::new(); if let Some(sk) = sketch { let ss = st.sketch(*sk).ok_or_else(|| DocError::Unknown(format!("sketch {sk} (it must come earlier in the timeline)")))?; @@ -20,7 +20,13 @@ pub(super) fn eval(doc: &Document, vals: &BTreeMap, f: &Feature, if !edges.is_empty() { let i = body_at(st, body, edges)?; let b = st.bodies.get(i).ok_or_else(|| DocError::Invalid("body".into()))?; - made.push(kernel::patch_edges_crowned(&b.body, edges, *crown)?); + let spec = kernel::PatchSpec { + continuity: continuity.all, + edges: continuity.edges.iter().map(|e| (e.edge, e.g)).collect(), + points: points.clone(), + crown: *crown, + }; + made.push(kernel::patch_edges_with(&b.body, edges, &spec)?); } if made.is_empty() { return Err(DocError::Invalid("a patch needs a sketch region or a loop of edges".into())); diff --git a/crates/doc/src/format/pending/patch-continuity-and-points.rs b/crates/doc/src/format/pending/patch-continuity-and-points.rs new file mode 100644 index 00000000..c9384646 --- /dev/null +++ b/crates/doc/src/format/pending/patch-continuity-and-points.rs @@ -0,0 +1,6 @@ +//! A `Patch` on edges can carry `continuity` (`all`: `g0`, `g1` or `g2`, and `edges`: +//! `[{edge: [x, y, z], g}]` for edges met otherwise; G0 everywhere when absent, and then not +//! saved) and `points` it passes through (none when absent). Older designs need no +//! transformation: a patch without them fills as it always did; a loop that is not a height field +//! over its plane, which 0.13.0 refused, now fills with one face. 0.13.x asks for an update +//! instead of dropping them. diff --git a/crates/engine/src/cmd/surface.rs b/crates/engine/src/cmd/surface.rs index 33c5aa2d..a2e9deac 100644 --- a/crates/engine/src/cmd/surface.rs +++ b/crates/engine/src/cmd/surface.rs @@ -2,7 +2,9 @@ //! solids), and the planned surface toolset (`planned!`: registered, not built yet). use serde_json::Value; -use solidcraft_doc::{BlendSide, ExtendMode, ExtendTo, FaceAt, FeatureKind, PlaneRef, ProfileSel, RuledCorners, RuledMode, expr::Kind}; +use solidcraft_doc::{ + BlendSide, EdgeContinuity, ExtendMode, ExtendTo, FaceAt, FeatureKind, PatchContinuity, PlaneRef, ProfileSel, RuledCorners, RuledMode, expr::Kind, +}; use solidcraft_kernel::Continuity; use super::CommandSpec; @@ -14,7 +16,7 @@ pub static COMMANDS: &[CommandSpec] = &[ CommandSpec::new("surface.patch", "Patch", patch) .at("SURFACE", "CREATE") .icon("patch") - .params("sketch?: id|name and profiles? (planar patches of sketch regions) | edges: [[x,y,z]] (points on a closed loop of a body's edges: one face on them, planar when the loop is, else a smooth B-spline surface) and body?; crown?: number −0.9…3 (default 0: only a fitted patch, whose middle becomes fuller by this many times the loop's size — the bump vanishes on the loop, so the face still meets the picked edges; a loop that is planar is the exact plane whatever the crown); name?, body_name?"), + .params("sketch?: id|name and profiles? (planar patches of sketch regions) | edges: [[x,y,z]] (points on a closed loop of a body's edges: one face on them, planar when the loop is, else a smooth B-spline surface) and body?; continuity?: g0|g1|g2 (edges only, default g0: how the patch meets the faces beside the loop — g1 their tangent planes, g2 also their curvature across the edges) or [{edge: [x,y,z], g: g0|g1|g2}] (those edges so, the rest g0) or {all, edges}; points?: [[x,y,z]] (edges only, up to 64: points the patch passes through); crown?: number −0.9…3 (default 0: only a fitted patch, whose middle becomes fuller by this many times the loop's size — the bump vanishes on the loop, so the face still meets the picked edges and the continuity across them; a loop that is planar and g0 is the exact plane whatever the crown; with points, crown needs a g1 or g2 edge); name?, body_name?"), CommandSpec::new("surface.stitch", "Stitch", stitch) .at("SURFACE", "MODIFY") .icon("stitch") @@ -84,24 +86,91 @@ pub static COMMANDS: &[CommandSpec] = &[ fn patch(s: &mut Session, p: &Value) -> Result { let cmd = "surface.patch"; let crown = crown_param(p, cmd)?; + let continuity = patch_continuity(p, cmd)?; + let points = patch_points(p, cmd)?; let edges: Vec = p.get("edges").and_then(Value::as_array).map(|a| a.iter().filter_map(vec3).collect()).unwrap_or_default(); let mut v = if !edges.is_empty() { if edges.len() > 10_000 { return Err(bad(cmd, "too many edges")); } let body = str_(p, "body").map(str::to_string); - add_feature(s, p, FeatureKind::Patch { sketch: None, profiles: ProfileSel::All, body, edges, crown })? + add_feature(s, p, FeatureKind::Patch { sketch: None, profiles: ProfileSel::All, body, edges, crown, continuity: continuity.clone(), points })? } else { + if !continuity.is_g0() || !points.is_empty() { + return Err(bad(cmd, "`continuity` and `points` are for a patch on a loop of edges")); + } let sketch = feature_sketch(s, p, cmd)?; let profiles = profiles(p, cmd)?; - add_feature(s, p, FeatureKind::Patch { sketch: Some(sketch), profiles, body: None, edges: Vec::new(), crown })? + let kind = FeatureKind::Patch { + sketch: Some(sketch), + profiles, + body: None, + edges: Vec::new(), + crown, + continuity: PatchContinuity::default(), + points, + }; + add_feature(s, p, kind)? }; if let Some(o) = v.as_object_mut() { o.insert("crown".into(), serde_json::json!(crown)); + o.insert("continuity".into(), serde_json::to_value(&continuity).unwrap_or(Value::Null)); } Ok(v) } +/// A continuity word for a patch: g0, g1 or g2 (a patch does not hold g3). +fn patch_order(v: Option<&Value>, cmd: &str) -> Result { + match v.and_then(Value::as_str).and_then(Continuity::parse) { + Some(g) if g != Continuity::G3 => Ok(g), + _ => Err(bad(cmd, "a patch's continuity is g0, g1 or g2")), + } +} + +/// The patch's `continuity`: one word for every edge, a list of `{edge, g}` for some edges (the +/// rest g0), or `{all, edges}`. +fn patch_continuity(p: &Value, cmd: &str) -> Result { + let edges = |list: &Value| -> Result> { + let list = list.as_array().ok_or_else(|| bad(cmd, "`continuity` edges must be a list of {edge: [x,y,z], g}"))?; + if list.len() > 10_000 { + return Err(bad(cmd, "`continuity` lists at most 10000 edges")); + } + list.iter() + .map(|e| { + let edge = e.get("edge").and_then(vec3).ok_or_else(|| bad(cmd, "a `continuity` entry's `edge` must be a point [x,y,z]"))?; + Ok(EdgeContinuity { edge, g: patch_order(e.get("g"), cmd)? }) + }) + .collect() + }; + match p.get("continuity") { + None => Ok(PatchContinuity::default()), + Some(v) if v.is_string() => Ok(PatchContinuity { all: patch_order(Some(v), cmd)?, edges: Vec::new() }), + Some(v) if v.is_array() => Ok(PatchContinuity { all: Continuity::G0, edges: edges(v)? }), + Some(v) if v.is_object() => { + let all = match v.get("all") { + None => Continuity::G0, + some => patch_order(some, cmd)?, + }; + let edges = match v.get("edges") { + None => Vec::new(), + Some(list) => edges(list)?, + }; + Ok(PatchContinuity { all, edges }) + } + Some(_) => Err(bad(cmd, "`continuity` must be g0, g1 or g2, or a list of {edge, g}")), + } +} + +/// The points a patch passes through (at most 64, each a finite point). +fn patch_points(p: &Value, cmd: &str) -> Result> { + let Some(v) = p.get("points") else { return Ok(Vec::new()) }; + let list = v + .as_array() + .filter(|a| a.len() <= solidcraft_geom::fill_general::MAX_FILL_POINTS) + .ok_or_else(|| bad(cmd, "`points` must list at most 64 points [x,y,z]"))?; + list.iter().map(|q| vec3(q).filter(|q| q.is_finite()).ok_or_else(|| bad(cmd, "a patch point must be a finite [x,y,z]"))).collect() +} + /// The patch's `crown`: how full a fitted patch's middle is made, -0.9…3 (default 0). A hostile /// value (NaN, infinity, out of range, a string) is refused. fn crown_param(p: &Value, cmd: &str) -> Result { @@ -1329,4 +1398,69 @@ mod tests { .is_err() ); } + + /// A tube whose top face is taken off (`surface.unstitch`), and points round its rim. + fn open_tube() -> (Session, Value) { + let mut s = Session::default(); + run(&mut s, "solid.cylinder", json!({"radius": 10, "height": 30, "body_name": "Tube"})); + run(&mut s, "surface.unstitch", json!({"body": "Tube", "faces": [[0, 0, 30]]})); + let rim: Vec = (0..8).map(|k| f64::from(k) * 0.785).map(|a| json!([10.0 * a.cos(), 10.0 * a.sin(), 30.0])).collect(); + (s, json!(rim)) + } + + /// `surface.patch continuity` and `points`: the tube's open end capped G1, per-edge + /// continuity, a point the cap passes through, all saved with the design and read back; + /// hostile values refused. + #[test] + fn patch_continuity_and_points_cap_a_tube() { + // A coarse mesh is enough to tell the patches apart (measuring a fitted face's area + // exactly tessellates it finely). + let mesh = |s: &Session, name: &str| s.world_state().body(name).unwrap().body.tessellate(0.05).unwrap(); + let (mut s, rim) = open_tube(); + let reply = run(&mut s, "surface.patch", json!({"edges": rim, "body": "Tube", "continuity": "g1", "body_name": "Cap"})); + assert_eq!(reply["continuity"], json!({"all": "g1"}), "{reply}"); + // A dome over the rim, not the flat disc a G0 patch of a planar rim is. + let disc = std::f64::consts::PI * 100.0; + let dome = mesh(&s, "Cap").measure().area; + assert!(dome > disc * 1.2, "a dome, not the disc: {dome}"); + // Through a point 8 above the rim, G1 to the wall: the cap reaches it. + run(&mut s, "surface.patch", json!({"edges": rim, "body": "Tube", "continuity": {"all": "g1"}, "points": [[0, 0, 38]], "body_name": "Tall"})); + let top = mesh(&s, "Tall").bounds().max.z; + assert!(top >= 38.0 - 1e-2, "the cap passes through the point: top {top}"); + // One edge G2, the rest G0 (the list form). + let first = rim[0].clone(); + let reply = + run(&mut s, "surface.patch", json!({"edges": rim, "body": "Tube", "continuity": [{"edge": first, "g": "g2"}], "body_name": "Mixed"})); + assert_eq!(reply["continuity"]["edges"][0]["g"], json!("g2"), "{reply}"); + // Saved and read back: the same patches. + let back = solidcraft_doc::Document::from_json(&s.doc.to_json()).unwrap(); + let reopened = Session::new(back); + for name in ["Cap", "Tall", "Mixed"] { + let (a, b) = (mesh(&s, name).positions, mesh(&reopened, name).positions); + assert_eq!(a, b, "{name}: the same surface read back"); + } + assert!(s.doc.to_json().contains("\"continuity\""), "the continuity is saved"); + // Hostile and misplaced values. + for (key, v) in [ + ("continuity", json!("g3")), + ("continuity", json!("smooth")), + ("continuity", json!(5)), + ("continuity", json!([{"edge": "here", "g": "g1"}])), + ("continuity", json!([{"edge": [f64::INFINITY, 0, 0], "g": "g1"}])), + ("continuity", json!({"all": "g9"})), + ("continuity", json!([{"edge": [0, 0, 0], "g": "g1"}])), + ("points", json!([[0, 0, f64::NAN]])), + ("points", json!("up")), + ("points", json!(vec![[0, 0, 35]; 65])), + ("points", json!([[100, 0, 30]])), + ] { + let mut p = json!({"edges": rim, "body": "Tube"}); + p[key] = v.clone(); + assert!(s.execute("surface.patch", &p).is_err(), "accepted {key}: {v}"); + } + run(&mut s, "sketch.create", json!({"plane": "XY"})); + run(&mut s, "sketch.rectangle.two_point", json!({"p0": [30, 0], "p1": [40, 10]})); + let sketch_patch = json!({"sketch": 1, "continuity": "g1"}); + assert!(s.execute("surface.patch", &sketch_patch).is_err(), "continuity on a sketch patch"); + } } diff --git a/crates/engine/src/tests.rs b/crates/engine/src/tests.rs index 000b8b58..02ff3555 100644 --- a/crates/engine/src/tests.rs +++ b/crates/engine/src/tests.rs @@ -779,6 +779,8 @@ fn hostile_params_never_panic() { "orientation", "twist", "scale", + "continuity", + "crown", ]; // 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/fill.rs b/crates/geom/src/fill.rs index 3122982b..7169a185 100644 --- a/crates/geom/src/fill.rs +++ b/crates/geom/src/fill.rs @@ -33,7 +33,7 @@ const GAUSS4: [(f64, f64); 4] = [ (0.861_136_311_594_052_6, 0.347_854_845_137_453_8), ]; /// The smoothness term's weight per fitted point (see [`fill_loop`]). -const SMOOTHNESS: f64 = 1e-7; +pub(crate) const SMOOTHNESS: f64 = 1e-7; /// A loop step shorter than this fraction of its length in space, seen along the normal, runs /// along the normal (steeper than about 87° from the plane). const MIN_STEP_SLOPE: f64 = 0.05; @@ -289,13 +289,32 @@ fn fit_once(fit: &[Vec3], all: &[Vec3], tol: f64, crown: f64) -> SplineResult SplineResult<(NurbsSurface, usize, f64)> { + if crown == 0.0 { + return Ok((surface, 0, 0.0)); + } + if !normal.is_finite() { + return invalid("the loop's extent is not finite"); + } + crown_fill_along(surface, uv, crown, &|_| Some(normal), MAX_FILL_SPANS) +} + +/// [`crown_fill`] moving each control point along `along(ξ, η)` (at its Greville point, the mean +/// of its knots each way) rather than one normal: what a fill that is not a height field takes, +/// with its own normal there. Its net may be refined up to `max_spans` each way. +pub(crate) fn crown_fill_along( + surface: NurbsSurface, + uv: &[[f64; 2]], + crown: f64, + along: &dyn Fn([f64; 2]) -> Option, + max_spans: usize, +) -> SplineResult<(NurbsSurface, usize, f64)> { if crown == 0.0 { return Ok((surface, 0, 0.0)); } let gap = loop_gap(uv); let mut base = surface; loop { - if let Some((crowned, points)) = crown_net(&base, uv, gap, crown, normal)? { + if let Some((crowned, points)) = crown_net(&base, uv, gap, crown, along)? { let mut drift = 0.0f64; for q in uv { if let (Ok(a), Ok(b)) = (base.eval(q[0], q[1]), crowned.eval(q[0], q[1])) { @@ -309,11 +328,11 @@ fn crown_fill(surface: NurbsSurface, uv: &[[f64; 2]], crown: f64, normal: Vec3) } // No room inside the loop yet: halve the cells (the surface itself does not change). let spans = [0, 1].map(|d| base.count()[d].saturating_sub(3)); - if spans[0] >= MAX_FILL_SPANS && spans[1] >= MAX_FILL_SPANS { + if spans[0] >= max_spans && spans[1] >= max_spans { return Err(SplineError::Unsupported("crowning a patch whose loop leaves no room inside it".into())); } for (d, &s) in spans.iter().enumerate() { - if s < MAX_FILL_SPANS { + if s < max_spans { base = subdivide(&base, d)?; } } @@ -324,13 +343,19 @@ fn crown_fill(surface: NurbsSurface, uv: &[[f64; 2]], crown: f64, normal: Vec3) /// loop takes `crown · size · B(ξ, η)`, `B` the bubble's blossom at the point's knots (so the /// surface between the perturbed points is exactly the bubble, not a fit of it) and `size` the /// box's diagonal. `None` when the net has no such point. -fn crown_net(base: &NurbsSurface, uv: &[[f64; 2]], gap: [f64; 2], crown: f64, normal: Vec3) -> SplineResult> { +fn crown_net( + base: &NurbsSurface, + uv: &[[f64; 2]], + gap: [f64; 2], + crown: f64, + along: &dyn Fn([f64; 2]) -> Option, +) -> SplineResult> { let [cu, cv] = base.count(); let degree = base.degree(); let knots = base.knots(); let domain = base.domain(); let size = (domain[0][1] - domain[0][0]).hypot(domain[1][1] - domain[1][0]); - if !(size > 0.0 && size.is_finite() && normal.is_finite()) { + if !(size > 0.0 && size.is_finite()) { return invalid("the loop's extent is not finite"); } // The bubble `4·s·(1−s)` per direction, as a quadratic in the plane coordinate. @@ -351,6 +376,8 @@ fn crown_net(base: &NurbsSurface, uv: &[[f64; 2]], gap: [f64; 2], crown: f64, no continue; } let bump = crown * size * blossom(qu, u_knots) * blossom(qv, v_knots); + let mean = |t: [f64; 3]| (t[0] + t[1] + t[2]) / 3.0; + let Some(normal) = along([mean(u_knots), mean(v_knots)]).filter(|n| n.is_finite()) else { continue }; let Some(slot) = ctrl.get_mut(j * cu + i) else { continue }; let h = *slot; let Some(p) = h.project() else { continue }; @@ -431,7 +458,7 @@ fn loop_gap(uv: &[[f64; 2]]) -> [f64; 2] { } /// Whether `p` is inside the closed loop `uv` (ray casting). -fn inside_loop(p: [f64; 2], uv: &[[f64; 2]]) -> bool { +pub(crate) fn inside_loop(p: [f64; 2], uv: &[[f64; 2]]) -> bool { let m = uv.len(); let mut inside = false; for k in 0..m { @@ -512,7 +539,7 @@ fn check_projectable(points: &[Vec3], uv: &[[f64; 2]]) -> SplineResult<()> { } /// Whether the closed segments `ab` and `cd` share a point. -fn segments_meet(a: [f64; 2], b: [f64; 2], c: [f64; 2], d: [f64; 2]) -> bool { +pub(crate) fn segments_meet(a: [f64; 2], b: [f64; 2], c: [f64; 2], d: [f64; 2]) -> bool { let orient = |p: [f64; 2], q: [f64; 2], r: [f64; 2]| (q[0] - p[0]) * (r[1] - p[1]) - (q[1] - p[1]) * (r[0] - p[0]); let within = |p: [f64; 2], q: [f64; 2], r: [f64; 2]| r[0] >= p[0].min(q[0]) && r[0] <= p[0].max(q[0]) && r[1] >= p[1].min(q[1]) && r[1] <= p[1].max(q[1]); @@ -525,9 +552,9 @@ fn segments_meet(a: [f64; 2], b: [f64; 2], c: [f64; 2], d: [f64; 2]) -> bool { /// A bicubic height field's grid: uniform clamped knots over `[lo, hi]` in each direction, /// `count` control points along each. -struct Grid { - knots: [Vec; 2], - count: [usize; 2], +pub(crate) struct Grid { + pub(crate) knots: [Vec; 2], + pub(crate) count: [usize; 2], } impl LoopFill { @@ -540,7 +567,7 @@ impl LoopFill { } impl Grid { - fn new(spans: [usize; 2], lo: [f64; 2], hi: [f64; 2]) -> Grid { + pub(crate) fn new(spans: [usize; 2], lo: [f64; 2], hi: [f64; 2]) -> Grid { let knots = [0, 1].map(|d| { let (a, b, s) = (lo[d], hi[d], spans[d]); let mut k = vec![a; 4]; @@ -552,7 +579,7 @@ impl Grid { } /// The length of one span in direction `d`. - fn cell(&self, d: usize) -> f64 { + pub(crate) fn cell(&self, d: usize) -> f64 { let (Some(k), Some(&c)) = (self.knots.get(d), self.count.get(d)) else { return 0.0 }; (k.last().copied().unwrap_or(0.0) - k.first().copied().unwrap_or(0.0)) / (c - 3) as f64 } @@ -564,7 +591,7 @@ impl Grid { /// The four non-zero basis functions' `r`-th derivatives and the first index, at `x` in /// direction `d`. - fn basis_derivative(&self, d: usize, x: f64, r: usize) -> Option<(usize, [f64; 4])> { + pub(crate) fn basis_derivative(&self, d: usize, x: f64, r: usize) -> Option<(usize, [f64; 4])> { let kn = self.knots.get(d)?; let x = x.clamp(*kn.first()?, *kn.last()?); let span = find_span(kn, 3, *self.count.get(d)?, x)?; @@ -622,6 +649,22 @@ impl Grid { // against the data term (squared heights) alike at any size. Each term is a product of // one-direction Gram matrices of the basis functions' derivatives. let lambda = SMOOTHNESS * fitted as f64 / n as f64 * self.cell(0) * self.cell(1); + self.add_energy(&mut band, lambda); + band.solve(rhs) + } + + /// Add `lambda` times the thin-plate energy `∫∫ h_uu² + 2·h_uv² + h_vv²` over the grid's box + /// to the normal equations `band` (unknowns i fastest): each term is a product of + /// one-direction Gram matrices of the basis functions' derivatives. + pub(crate) fn add_energy(&self, band: &mut Band, lambda: f64) { + let cu = self.count[0]; + self.energy(lambda, &mut |(i, j), (k, l), e| band.add(j * cu + i, l * cu + k, e)); + } + + /// The entries of `lambda` times the thin-plate energy's matrix (see [`Grid::add_energy`]), + /// control point `(i, j)` against `(k, l)`, for every pair whose supports overlap. + pub(crate) fn energy(&self, lambda: f64, add: &mut dyn FnMut((usize, usize), (usize, usize), f64)) { + let [cu, cv] = self.count; let gram = |d: usize, r: usize| -> Vec { let c = self.count.get(d).copied().unwrap_or(0); let mut g = vec![0.0; c * c]; @@ -652,16 +695,14 @@ impl Grid { }; for j in 0..cv { for i in 0..cu { - let r = j * cu + i; for l in j.saturating_sub(3)..(j + 4).min(cv) { for k in i.saturating_sub(3)..(i + 4).min(cu) { let e: f64 = ENERGY.iter().map(|&(ru, rv, w)| w * g(0, ru, i, k) * g(1, rv, j, l)).sum(); - band.add(r, l * cu + k, lambda * e); + add((i, j), (k, l), lambda * e); } } } } - band.solve(rhs) } /// The surface of the net `heights` over the plane `o`, `x`, `y`, `n`: control point `(i, j)` @@ -680,14 +721,15 @@ impl Grid { } /// A symmetric banded matrix: the upper band, row `r` holding columns `r ..= r + bw`. -struct Band { +#[derive(Clone)] +pub(crate) struct Band { n: usize, bw: usize, a: Vec, } impl Band { - fn new(n: usize, bw: usize) -> Band { + pub(crate) fn new(n: usize, bw: usize) -> Band { Band { n, bw, a: vec![0.0; n * (bw + 1)] } } fn index(&self, r: usize, s: usize) -> Option { @@ -704,7 +746,7 @@ impl Band { } /// Add `v` at `(r, s)` with `r ≤ s` (the lower half is the same by symmetry: callers adding /// a full outer product reach each stored entry once through `r ≤ s`). - fn add(&mut self, r: usize, s: usize, v: f64) { + pub(crate) fn add(&mut self, r: usize, s: usize, v: f64) { if r <= s && let Some(x) = self.index(r, s).and_then(|i| self.a.get_mut(i)) { @@ -713,7 +755,13 @@ impl Band { } /// Solve `A·x = b` by banded Cholesky factorisation `A = Uᵀ·U`; an error when `A` is not /// positive definite. - fn solve(mut self, mut b: Vec) -> SplineResult> { + fn solve(self, b: Vec) -> SplineResult> { + self.factor()?.substitute(b) + } + + /// The Cholesky factor `U` of `A = Uᵀ·U`, in place; an error when `A` is not positive + /// definite. + pub(crate) fn factor(mut self) -> SplineResult { let (n, bw) = (self.n, self.bw); for r in 0..n { let mut d = self.get(r, r); @@ -734,6 +782,48 @@ impl Band { self.set(r, s, v / d); } } + Ok(self) + } + + /// [`Band::factor`] by rows: once row `r` of `U` is known, the rows below it within the band + /// are updated from it at once (right-looking), which walks memory row by row instead of down + /// columns. The same factor, to rounding; what the general fill, whose bands are wide, uses. + pub(crate) fn factor_rows(mut self) -> SplineResult { + let (n, bw) = (self.n, self.bw); + let w = bw + 1; + let mut row = vec![0.0; w]; + for r in 0..n { + let Some(here) = self.a.get_mut(r * w..r * w + w) else { return invalid("the fit's equations are out of shape") }; + let d = here.first().copied().unwrap_or(0.0); + if !(d > 0.0) { + return invalid("the fit's equations are singular"); + } + let d = d.sqrt(); + for (k, x) in here.iter_mut().enumerate() { + *x = if k == 0 { d } else { *x / d }; + } + row.copy_from_slice(here); + let end = (r + bw + 1).min(n); + for s in r + 1..end { + let us = row.get(s - r).copied().unwrap_or(0.0); + if us == 0.0 { + continue; + } + let len = end - s; + let (Some(dst), Some(src)) = (self.a.get_mut(s * w..s * w + len), row.get(s - r..s - r + len)) else { + return invalid("the fit's equations are out of shape"); + }; + for (x, y) in dst.iter_mut().zip(src) { + *x -= us * y; + } + } + } + Ok(self) + } + + /// Solve `Uᵀ·U·x = b` with a factor from [`Band::factor`]. + pub(crate) fn substitute(&self, mut b: Vec) -> SplineResult> { + let (n, bw) = (self.n, self.bw); // Uᵀ·y = b, then U·x = y. for r in 0..n { let mut v = b.get(r).copied().unwrap_or(0.0); diff --git a/crates/geom/src/fill_general.rs b/crates/geom/src/fill_general.rs new file mode 100644 index 00000000..e6f5da58 --- /dev/null +++ b/crates/geom/src/fill_general.rs @@ -0,0 +1,1150 @@ +//! Filling a closed loop with one surface that meets the faces beside the loop with a chosen +//! continuity (G0, G1 or G2 per station) and passes through given points, for any loop that +//! turns round once: loops that are not a height field over any plane, and loops of any number of +//! sides, are filled whole. +//! +//! The formulation: +//! +//! - **The parameter domain is developed from the loop itself**, not projected. The loop is +//! walked as a polygon of its samples; each step keeps its length, and each turn between steps +//! is split, against the patch's expected normal there (`up`), into its part about that normal +//! (the turn *in* the patch: the boundary's geodesic turn, signed, so a reflex corner stays +//! reflex) and the rest (the turn the boundary makes *with* the patch bending under it, which a +//! plane cannot show). By Gauss–Bonnet the in-patch turns of a developable patch add up to a +//! full turn; a curved patch makes up the difference with its curvature, so a shortfall is given +//! back in proportion to the bending turns (a cap over a round rim, whose rim does not turn in +//! its neighbours' tangent planes at all, develops to the circle), and a saddle's excess comes +//! off the smooth stretches, so corners keep their angles. The headings then close exactly, and +//! the small gap left in position is spread along the steps by length. A domain that crosses +//! itself is refused. The domain's steps are as long as the loop's, so the parameters are +//! millimetres along the boundary. (`up` is the neighbour's normal at a G1/G2 station; at a G0 +//! station, the loop's tangent × the neighbour's leave direction tilted towards the loop's mean +//! plane.) +//! - **The surface is a bicubic B-spline over the domain's box** whose control points minimise, +//! by least squares, the misfit to the boundary conditions plus a vanishing multiple of the +//! thin-plate energy `∫∫ |S_uu|² + 2|S_uv|² + |S_vv|²` (which only decides what the conditions +//! leave free): +//! - G0, every station: the position. +//! - G1: `n·∂S/∂ν = 0` (`ν` the domain's inward normal, `n` the neighbour's normal): the +//! cross-boundary derivative lies in the neighbour's tangent plane, which with the boundary +//! tangent makes the normals agree; plus weak Hermite rows `∂S/∂ν = s·w` (`w` the direction +//! the neighbour leaves in across the edge, `s` the boundary's own speed), which say that the +//! patch leaves into itself rather than folding back over the neighbour, and about how fast, +//! but give way where full Hermite data cannot all hold (it is consistent only with a map +//! exactly conformal at the boundary, which a curved patch's is not). +//! - G2: `n·∂²S/∂ν² = (∂S/∂ν·w)²·κ` (`κ` the neighbour's normal curvature along `w`, read off +//! the bridge's frames): the patch's normal curvature across the edge is the neighbour's. The +//! speed is the fit's own, so the row is solved again with the last solution's speed. +//! +//! Without G1/G2 the three coordinates are fitted apart on one shared matrix; with them, together +//! (the geometric rows couple them), numbered so the band stays narrow. At corners the two +//! sides' conditions overdetermine the surface: see [`Near`] and [`fades`]. The face is bounded +//! by the body's own edges, so the surface only has to pass within the tolerance of them. +//! - **Points the patch must pass through** are located on the fitted surface (nearest point, +//! inside the domain), added as position data there, and the fit and the location are repeated +//! — Hoschek's parameter correction — so each point ends on the surface. +//! - **Grids** of 4 to [`MAX_GENERAL_SPANS`] spans across the domain's longer side are tried in +//! turn; the first whose surface is within the tolerance of every station, whose normals are +//! within the angle tolerance of the neighbours' at every G1 or G2 station (away from corners +//! where the neighbours' tangent planes disagree: no surface can meet both there), that passes +//! through the points, and that neither folds nor pinches inside the loop, is the fill. When a +//! grid needs more stations than it has, the sampling doubles. +//! - **A crown** fills the middle out along the surface's own normals, on control points clear of +//! the loop, so it moves neither the loop nor the conditions across it. + +use crate::Vec3; +use crate::bspline::{HPoint, NurbsSurface, SplineError, SplineResult, invalid}; +use crate::fill::{Band, Grid, MAX_LOOP_POINTS, SMOOTHNESS, inside_loop, segments_meet}; + +/// Most spans of the surface along the domain's longer side. +pub const MAX_GENERAL_SPANS: usize = 64; +/// Most points a patch can be asked to pass through. +pub const MAX_FILL_POINTS: usize = 64; +/// Fewest spans along the domain's longer side (the first grid tried). +const MIN_SPANS: usize = 4; +/// A turn between two loop steps sharper than this (radians) folds the loop back on itself. +const FOLD_TURN: f64 = 3.05; +/// How many points of the domain's box each way the fold check samples. +const FOLD_SAMPLES: usize = 32; +/// Rounds of locating the points on the surface and fitting again. +const POINT_ROUNDS: usize = 3; +/// A domain is rectilinear when this fraction of its length (weighted by how squarely it runs) +/// lies along two directions at right angles. +const RECTILINEAR: f64 = 0.5; +/// How far (radians) the domain's angle at a corner may be from the loop's for the corner to be +/// met conformally (see [`Near`]). +const CORNER_ANGLE_TOL: f64 = 0.5 * std::f64::consts::PI / 180.0; +/// Near a corner whose neighbours disagree, the G1 rows are left out over this share of the run to +/// the next corner (see [`fades`]). +const BLEND_SHARE: f64 = 0.4; +/// How close, in grid cells, a station is to a corner when it is near it (see [`Near`]): G1 and +/// G2 stations that near a corner whose neighbours disagree are not held to the angle tolerance. +const CLASH_CELLS: f64 = 2.0; + +/// What the face beside a station asks the patch to repeat across it. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Across { + /// The unit direction the neighbour leaves in across the edge (its tangent plane, at right + /// angles to the edge, away from the neighbour): where the patch goes. + pub leave: Vec3, + /// The neighbour's unit normal (either sign). + pub normal: Vec3, + /// G2: the neighbour's curvature vector along `leave` (unit speed). + pub bend: Option, +} + +/// One sample of the loop. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Station { + pub p: Vec3, + /// The patch's expected unit normal here, so that the loop runs counter-clockwise about it + /// with the patch on its left: the loop's tangent × the direction into the patch. It need not + /// be exact; it decides which way the loop turns in the patch. + pub up: Vec3, + /// G1 or G2 conditions (`None`: position only, G0). + pub across: Option, + /// The loop has a corner here: position only (two neighbours' conditions meet). + pub corner: bool, + /// A corner where the neighbours' tangent planes disagree: no surface is G1 to both, so the + /// angle tolerance is not asked of the stations close to it. + pub clash: bool, +} + +/// What [`fill_general`] is asked for. +#[derive(Clone, Copy, Debug)] +pub struct FillSpec<'a> { + /// How far the surface may be from the loop (and from the points). + pub tol: f64, + /// How far (radians) the surface's normal may turn from a G1/G2 neighbour's. + pub angle_tol: f64, + /// How full the middle is made (see [`crate::fill::fill_loop`]'s crown, between + /// [`crate::fill::MIN_CROWN`] and [`crate::fill::MAX_CROWN`]): the bump, `crown` times the + /// domain's diagonal at its middle, is laid along the surface's own normal, turned to the side + /// the G1/G2 neighbours' normals face, on control points clear of the loop, so neither the + /// loop nor the conditions across it move. A loop with no G1/G2 station has no side to crown + /// to and is refused one. + pub crown: f64, + /// Points the surface passes through. + pub points: &'a [Vec3], +} + +/// The fill and what it achieved. +#[derive(Clone, Debug, PartialEq)] +pub struct GeneralFill { + pub surface: NurbsSurface, + /// Whether the loop, in the order given, runs counter-clockwise about the surface's normal + /// `∂S/∂u × ∂S/∂v`. + pub counter_clockwise: bool, + /// The farthest station from the surface. + pub deviation: f64, + /// The largest angle (radians) between the surface's normal and a neighbour's at a G1/G2 + /// station held to the angle tolerance, and at any G1/G2 station. + pub normal_deviation: f64, + pub normal_deviation_all: f64, + /// The largest difference (1/length) between the surface's normal curvature across the + /// boundary and the neighbour's, at the G2 stations held to the angle tolerance. + pub curvature_mismatch: f64, + /// The farthest point from the surface. + pub point_deviation: f64, + /// Spans of the grid each way, and stations sampled. + pub spans: [usize; 2], + pub stations: usize, +} + +/// Fill the closed loop of `sides` sides with one surface (see the module documentation). +/// `stations(k, ts)` samples side `k` at the parameters `ts` (0 at its start, which is where side +/// `k−1` ends, towards 1 at its end); it is asked for the same parameters again only when the +/// sampling grows. +pub fn fill_general( + sides: usize, + stations: &mut dyn FnMut(usize, &[f64]) -> Result, String>, + spec: &FillSpec, +) -> SplineResult { + if sides == 0 || sides.saturating_mul(16) > MAX_LOOP_POINTS / 2 { + return Err(SplineError::TooLarge(format!("a loop of {sides} sides (at most {})", MAX_LOOP_POINTS / 32))); + } + if !(spec.tol.is_finite() && spec.tol > 0.0) || !(spec.angle_tol.is_finite() && spec.angle_tol > 0.0) { + return invalid("the fill tolerances must be positive"); + } + // A range test also throws out NaN. + if !(crate::fill::MIN_CROWN..=crate::fill::MAX_CROWN).contains(&spec.crown) { + return invalid("the crown must be between -0.9 and 3"); + } + if spec.points.len() > MAX_FILL_POINTS || spec.points.iter().any(|p| !p.is_finite()) { + return invalid(format!("a patch passes through at most {MAX_FILL_POINTS} finite points")); + } + let mut sample = |ts: &[f64]| -> SplineResult>> { + (0..sides) + .map(|k| { + let got = stations(k, ts).map_err(SplineError::Invalid)?; + if got.len() == ts.len() && got.iter().all(finite_station) { + Ok(got) + } else { + Err(SplineError::Invalid(format!("side {k} has no finite stations"))) + } + }) + .collect() + }; + let mut per = 16usize; + let mut fit = sample(&(0..per).map(|i| i as f64 / per as f64).collect::>())?; + let mut best = Shortfall::default(); + // Each round starts at the grid the last one stopped at: the coarser grids fell short with + // fewer stations, and more stations do not make a coarse grid hold the loop closer. + let mut level = MIN_SPANS; + loop { + let check = sample(&(0..per).map(|i| (i as f64 + 0.5) / per as f64).collect::>())?; + // Each side in order: fitted and checked stations alternating. + let all: Vec> = fit.iter().zip(&check).map(|(f, c)| f.iter().zip(c).flat_map(|(a, b)| [*a, *b]).collect()).collect(); + let flat: Vec = all.iter().flatten().copied().collect(); + match fit_round(&flat, spec, level)? { + Round::Done(fill) => return Ok(*fill), + Round::Denser(s, at) => { + best = best.min(s); + level = at; + } + Round::Exhausted(s) => { + best = best.min(s); + break; + } + } + if flat.len() * 2 > MAX_LOOP_POINTS { + break; + } + fit = all; + per *= 2; + } + let mut why = format!("filling this loop with one surface: the closest fit stays {:e} off the loop (tolerance {:e})", best.deviation, spec.tol); + if fit.iter().flatten().any(|s| s.across.is_some()) { + why += &format!(", turns {:.3}° from its neighbours (tolerance {:.3}°)", best.angle.to_degrees(), spec.angle_tol.to_degrees()); + } + if !spec.points.is_empty() { + why += &format!(", {:e} off its points", best.points); + } + Err(SplineError::Unsupported(why)) +} + +fn finite_station(s: &Station) -> bool { + s.p.is_finite() + && s.up.is_finite() + && s.across.is_none_or(|a| a.leave.is_finite() && a.normal.is_finite() && a.bend.is_none_or(|b| b.is_finite())) +} + +/// How close a grid came: what [`fill_general`] reports when none is close enough. +#[derive(Clone, Copy, Debug)] +struct Shortfall { + deviation: f64, + angle: f64, + points: f64, +} + +impl Default for Shortfall { + fn default() -> Shortfall { + Shortfall { deviation: f64::INFINITY, angle: f64::INFINITY, points: f64::INFINITY } + } +} + +impl Shortfall { + /// The closer of two, by how many times its tolerance each misses. + fn min(self, other: Shortfall) -> Shortfall { + let worst = |s: &Shortfall| s.deviation.max(s.points) + s.angle; + if worst(&other) < worst(&self) { other } else { self } + } +} + +enum Round { + Done(Box), + /// Denser sampling is needed for the grid of this level (the coarser ones fell short). + Denser(Shortfall, usize), + Exhausted(Shortfall), +} + +/// The loop developed into the plane (see the module documentation). +struct Developed { + /// The stations' places in the domain. + d: Vec<[f64; 2]>, + /// The domain's inward unit normal at each station. + nu: Vec<[f64; 2]>, + /// The boundary's speed at each station: loop length per domain length. + speed: Vec, + /// Whether the domain turns at each station by the loop's own angle there (to + /// [`CORNER_ANGLE_TOL`]): a corner where it does not cannot be met conformally. + true_angle: Vec, + counter_clockwise: bool, +} + +/// Develop the loop `st` into the plane. +fn develop(st: &[Station]) -> SplineResult { + let m = st.len(); + if m < 3 { + return invalid("a loop needs three stations"); + } + let at = |k: usize| st.get(k % m).map_or(Vec3::ZERO, |s| s.p); + let chord: Vec = (0..m).map(|k| at(k + 1) - at(k)).collect(); + let len: Vec = chord.iter().map(|c| c.len()).collect(); + let total: f64 = len.iter().sum(); + if !(total > 0.0 && total.is_finite()) { + return invalid("the loop has no length"); + } + if len.iter().any(|l| !(*l > 1e-12 * total)) { + return Err(SplineError::Unsupported("filling a loop with two stations in one place".into())); + } + // The turn at each station, split into its part in the patch (signed about `up`) and the rest. + // A G0 station's `up` is the loop's tangent × its neighbour's leave direction: the patch's + // normal if it carried on from the neighbour, which says nothing about the turns of a loop + // whose neighbours stand square to the patch (a lid on walls). It is tilted towards the + // loop's own mean plane (Newell's area vector, its part square to the tangent, at half + // weight): that decides those turns, and where a loop wraps round so far that the mean plane + // faces the other way, the neighbour's direction still wins. + let area = (0..m).fold(Vec3::ZERO, |a, k| a + at(k).cross(at(k + 1))); + let loop_normal = area.normalized().unwrap_or(Vec3::ZERO); + let up: Vec = (0..m) + .map(|k| { + let Some(s) = st.get(k) else { return Vec3::ZERO }; + if s.across.is_some() { + return s.up; + } + let (Some(a), Some(b)) = (chord.get((k + m - 1) % m), chord.get(k)) else { return s.up }; + let Some(t) = (*a + *b).normalized() else { return s.up }; + let square = loop_normal - t * loop_normal.dot(t); + (s.up + square * 0.5).normalized().unwrap_or(s.up) + }) + .collect(); + let mut full_turn = vec![0.0; m]; + let mut in_patch = vec![0.0; m]; + let mut bending = vec![0.0; m]; + for k in 0..m { + let (Some(a), Some(b), Some(&up)) = (chord.get((k + m - 1) % m), chord.get(k), up.get(k)) else { continue }; + let (la, lb) = (a.len(), b.len()); + let cos = (a.dot(*b) / (la * lb)).clamp(-1.0, 1.0); + let turn = cos.acos(); + if turn > FOLD_TURN { + return Err(SplineError::Unsupported("filling a loop that folds back on itself".into())); + } + if let Some(f) = full_turn.get_mut(k) { + *f = turn; + } + let axis = a.cross(*b); + let Some(axis) = axis.normalized().filter(|_| axis.len() > 1e-15 * la * lb) else { continue }; + let along = axis.dot(up).clamp(-1.0, 1.0); + if let (Some(g), Some(r)) = (in_patch.get_mut(k), bending.get_mut(k)) { + *g = turn * along; + *r = turn * (1.0 - along * along).max(0.0).sqrt(); + } + } + let geodesic: f64 = in_patch.iter().sum(); + let bend_total: f64 = bending.iter().sum(); + let sense = if geodesic >= 0.0 { 1.0 } else { -1.0 }; + let full = std::f64::consts::TAU; + // A saddle's boundary turns more than once round in the patch: the excess comes off the + // smooth stretches, so the corners keep their angles (a G1 patch is conformal at its + // boundary, and a map conformal at a corner keeps its angle); only when the corners alone + // turn too far are they scaled too. + let is_corner = |k: usize| st.get(k).is_some_and(|s| s.corner || s.clash); + let at_corners: f64 = (0..m).filter(|k| is_corner(*k)).filter_map(|k| in_patch.get(k)).sum::() * sense; + let smooth = geodesic * sense - at_corners; + let turns: Vec = if geodesic.abs() >= full && at_corners <= full && smooth > 0.0 { + let f = (full - at_corners) / smooth; + in_patch.iter().enumerate().map(|(k, g)| if is_corner(k) { *g } else { g * f }).collect() + } else if geodesic.abs() >= full { + in_patch.iter().map(|g| g * full / geodesic.abs()).collect() + } else if bend_total > 1e-9 { + let mu = (full - geodesic.abs()) / bend_total; + in_patch.iter().zip(&bending).map(|(g, r)| g + sense * mu * r).collect() + } else { + return Err(SplineError::Unsupported("filling a loop that does not turn round".into())); + }; + // Steps at their lengths along the headings, the gap spread by length. + let mut heading = 0.0; + let mut steps: Vec<[f64; 2]> = Vec::with_capacity(m); + for k in 0..m { + heading += turns.get(k).copied().unwrap_or(0.0); + let l = len.get(k).copied().unwrap_or(0.0); + steps.push([l * heading.cos(), l * heading.sin()]); + } + let gap = steps.iter().fold([0.0, 0.0], |g, s| [g[0] + s[0], g[1] + s[1]]); + for (s, l) in steps.iter_mut().zip(&len) { + s[0] -= gap[0] * l / total; + s[1] -= gap[1] * l / total; + } + let mut d = Vec::with_capacity(m); + let mut q = [0.0, 0.0]; + for s in &steps { + d.push(q); + q = [q[0] + s[0], q[1] + s[1]]; + } + // Centred and turned so the domain, and with it the grid and the fill, does not depend on + // where the loop starts: a mostly rectilinear domain (most of its length runs along two + // directions at right angles, as a polygon of straight sides often does) with those + // directions along the grid's lines, which hold data along them exactly; any other onto its + // principal axes. The larger extent is along u. + let centre = d.iter().fold([0.0, 0.0], |c, p| [c[0] + p[0] / m as f64, c[1] + p[1] / m as f64]); + let (mut sxx, mut sxy, mut syy) = (0.0, 0.0, 0.0); + for p in &d { + let (x, y) = (p[0] - centre[0], p[1] - centre[1]); + sxx += x * x; + sxy += x * y; + syy += y * y; + } + let (mut c4, mut s4) = (0.0, 0.0); + for st in &steps { + let (l, a) = (st[0].hypot(st[1]), st[1].atan2(st[0])); + c4 += l * (4.0 * a).cos(); + s4 += l * (4.0 * a).sin(); + } + let angle = if c4.hypot(s4) > RECTILINEAR * total { 0.25 * s4.atan2(c4) } else { 0.5 * (2.0 * sxy).atan2(sxx - syy) }; + let turned = |angle: f64| -> Vec<[f64; 2]> { + let (c, s) = (angle.cos(), angle.sin()); + d.iter() + .map(|p| { + let (x, y) = (p[0] - centre[0], p[1] - centre[1]); + [c * x + s * y, -s * x + c * y] + }) + .collect() + }; + let mut d = turned(angle); + let spread = + |d: &[[f64; 2]], k: usize| d.iter().map(|p| p[k]).fold(f64::NEG_INFINITY, f64::max) - d.iter().map(|p| p[k]).fold(f64::INFINITY, f64::min); + if spread(&d, 1) > spread(&d, 0) * (1.0 + 1e-9) { + d = turned(angle + std::f64::consts::FRAC_PI_2); + } + let skew: f64 = d.iter().map(|p| p[0] * p[0] * p[0] + p[1] * p[1] * p[1]).sum(); + if skew < 0.0 { + for p in &mut d { + *p = [-p[0], -p[1]]; + } + } + check_simple(&d)?; + // The domain runs the way the loop turned in the patch. + let area: f64 = (0..m).map(|k| cross2(get2(&d, k), get2(&d, k + 1))).sum::() * 0.5; + if !(area * sense > 0.0) { + return Err(SplineError::Unsupported("filling a loop that develops inside out".into())); + } + let mut nu = Vec::with_capacity(m); + let mut speed = Vec::with_capacity(m); + for k in 0..m { + let (prev, next) = (get2(&d, k + m - 1), get2(&d, k + 1)); + let here = get2(&d, k); + let (a, b) = ([here[0] - prev[0], here[1] - prev[1]], [next[0] - here[0], next[1] - here[1]]); + let (la, lb) = (a[0].hypot(a[1]), b[0].hypot(b[1])); + if !(la > 0.0 && lb > 0.0) { + return Err(SplineError::Unsupported("filling a loop with two stations in one place".into())); + } + let t = [a[0] / la + b[0] / lb, a[1] / la + b[1] / lb]; + let lt = t[0].hypot(t[1]); + let t = if lt > 1e-9 { [t[0] / lt, t[1] / lt] } else { [b[0] / lb, b[1] / lb] }; + // The left normal is inward on a counter-clockwise domain. + nu.push([-t[1] * sense, t[0] * sense]); + let space = len.get((k + m - 1) % m).copied().unwrap_or(0.0) + len.get(k).copied().unwrap_or(0.0); + speed.push(space / (la + lb)); + } + let true_angle = turns.iter().zip(&full_turn).map(|(t, f)| (t.abs() - f).abs() <= CORNER_ANGLE_TOL).collect(); + Ok(Developed { d, nu, speed, true_angle, counter_clockwise: sense > 0.0 }) +} + +fn get2(d: &[[f64; 2]], k: usize) -> [f64; 2] { + d.get(k % d.len().max(1)).copied().unwrap_or([0.0; 2]) +} + +fn cross2(a: [f64; 2], b: [f64; 2]) -> f64 { + a[0] * b[1] - a[1] * b[0] +} + +/// Refuse a domain polygon whose non-adjacent sides meet. +fn check_simple(d: &[[f64; 2]]) -> SplineResult<()> { + let m = d.len(); + let q = |k: usize| get2(d, k); + let mut segs: Vec<(f64, f64, usize)> = (0..m) + .map(|k| { + let (a, b) = (q(k), q(k + 1)); + (a[0].min(b[0]), a[0].max(b[0]), k) + }) + .collect(); + segs.sort_by(|a, b| a.0.total_cmp(&b.0)); + for (i, &(_, hi, k)) in segs.iter().enumerate() { + for &(lo2, _, l) in segs.get(i + 1..).into_iter().flatten() { + if lo2 > hi { + break; + } + let adjacent = (k + 1) % m == l || (l + 1) % m == k; + if !adjacent && segments_meet(q(k), q(k + 1), q(l), q(l + 1)) { + return Err(SplineError::Unsupported("filling a loop whose development crosses itself".into())); + } + } + } + Ok(()) +} + +/// The basis of a grid at `q`: the first indices and the values of the four non-zero functions +/// and their first two derivatives, each way. +struct Basis { + first: [usize; 2], + /// `[direction][derivative][function]`. + b: [[[f64; 4]; 3]; 2], +} + +impl Basis { + fn at(grid: &Grid, q: [f64; 2]) -> Option { + let mut first = [0usize; 2]; + let mut b = [[[0.0; 4]; 3]; 2]; + for d in 0..2 { + for r in 0..3 { + let (f, vals) = grid.basis_derivative(d, *q.get(d)?, r)?; + *first.get_mut(d)? = f; + *b.get_mut(d)?.get_mut(r)? = vals; + } + } + Some(Basis { first, b }) + } + + /// The 16 (i, j, coefficient) terms of the derivative `Σ w[ru][rv] ∂^ru_u ∂^rv_v` at the point. + fn cells(&self, w: &[(usize, usize, f64)]) -> [(usize, usize, f64); 16] { + let mut out = [(0usize, 0usize, 0.0f64); 16]; + for (t, slot) in out.iter_mut().enumerate() { + let (a, b) = (t % 4, t / 4); + let c: f64 = w.iter().map(|&(ru, rv, k)| k * self.b[0][ru.min(2)][a] * self.b[1][rv.min(2)][b]).sum(); + *slot = (self.first[0] + a, self.first[1] + b, c); + } + out + } + + /// The 16 (index, coefficient) terms of the derivative `Σ w[ru][rv] ∂^ru_u ∂^rv_v` at the + /// point, over control points numbered i fastest in a grid `cu` wide. + fn terms(&self, cu: usize, w: &[(usize, usize, f64)]) -> [(usize, f64); 16] { + let mut out = [(0usize, 0.0f64); 16]; + for (t, slot) in out.iter_mut().enumerate() { + let (a, b) = (t % 4, t / 4); + let c: f64 = w.iter().map(|&(ru, rv, k)| k * self.b[0][ru.min(2)][a] * self.b[1][rv.min(2)][b]).sum(); + *slot = ((self.first[1] + b) * cu + self.first[0] + a, c); + } + out + } +} + +/// The normal equations of the fit, shared by the three coordinates. +struct System { + band: Band, + rhs: [Vec; 3], + rows: usize, +} + +impl System { + fn new(n: usize, bw: usize) -> System { + System { band: Band::new(n, bw), rhs: [vec![0.0; n], vec![0.0; n], vec![0.0; n]], rows: 0 } + } + + /// One least-squares row: `Σ coefficient · control = target`. + fn row(&mut self, terms: &[(usize, f64); 16], target: Vec3) { + for &(r, wr) in terms { + for (rhs, t) in self.rhs.iter_mut().zip([target.x, target.y, target.z]) { + if let Some(slot) = rhs.get_mut(r) { + *slot += wr * t; + } + } + for &(s, ws) in terms { + self.band.add(r, s, wr * ws); + } + } + self.rows += 1; + } +} + +/// A fitted surface's net on its grid. +struct Net<'a> { + grid: &'a Grid, + ctrl: Vec, +} + +impl Net<'_> { + /// `∂^ru_u ∂^rv_v S` at `q`. + fn der(&self, q: [f64; 2], ru: usize, rv: usize) -> Option { + let b = Basis::at(self.grid, q)?; + let cu = *self.grid.count.first()?; + let mut s = Vec3::ZERO; + for (i, c) in b.terms(cu, &[(ru, rv, 1.0)]) { + s += *self.ctrl.get(i)? * c; + } + Some(s) + } + + fn point(&self, q: [f64; 2]) -> Option { + self.der(q, 0, 0) + } + + /// The unit normal `S_u × S_v` and the first derivatives. + fn frame(&self, q: [f64; 2]) -> Option<(Vec3, Vec3, Vec3)> { + let (su, sv) = (self.der(q, 1, 0)?, self.der(q, 0, 1)?); + Some((su.cross(sv).normalized()?, su, sv)) + } +} + +/// The weights `(ru, rv, w)` of the first and second derivatives along the domain direction `nu`. +fn along(nu: [f64; 2], order: usize) -> Vec<(usize, usize, f64)> { + let [a, b] = nu; + match order { + 1 => vec![(1, 0, a), (0, 1, b)], + _ => vec![(2, 0, a * a), (1, 1, 2.0 * a * b), (0, 2, b * b)], + } +} + +/// One round: develop the loop `st` (the fitted stations at even places), then try the grids. +fn fit_round(st: &[Station], spec: &FillSpec, start: usize) -> SplineResult { + let dev = develop(st)?; + let m = st.len(); + let (mut lo, mut hi) = ([f64::INFINITY; 2], [f64::NEG_INFINITY; 2]); + for q in &dev.d { + for k in 0..2 { + lo[k] = lo[k].min(q[k]); + hi[k] = hi[k].max(q[k]); + } + } + let pad = 0.05 * (hi[0] - lo[0]).max(hi[1] - lo[1]); + for k in 0..2 { + lo[k] -= pad; + hi[k] += pad; + } + let extent = [hi[0] - lo[0], hi[1] - lo[1]]; + if !(extent[0] > 0.0 && extent[1] > 0.0 && extent[0].is_finite() && extent[1].is_finite()) { + return invalid("the loop's development has no extent"); + } + // The largest step between fitted stations, per direction. + let mut step = [0.0f64; 2]; + for k in (0..m).step_by(2) { + let (p, q) = (get2(&dev.d, k), get2(&dev.d, k + 2)); + for i in 0..2 { + step[i] = step[i].max((q[i] - p[i]).abs()); + } + } + // Where the corners are in the domain, and those whose neighbours disagree. + let corners: Vec<[f64; 2]> = st.iter().zip(&dev.d).filter(|(s, _)| s.corner || s.clash).map(|(_, d)| *d).collect(); + let rough: Vec<[f64; 2]> = st + .iter() + .zip(&dev.d) + .zip(&dev.true_angle) + .filter(|((s, _), true_angle)| s.clash || (s.corner && !**true_angle)) + .map(|((_, d), _)| *d) + .collect(); + let clashes: Vec<[f64; 2]> = st.iter().zip(&dev.d).filter(|(s, _)| s.clash).map(|(_, d)| *d).collect(); + let fade = fades(st); + let longer = extent[0].max(extent[1]); + let mut best = Shortfall::default(); + let mut level = start.max(MIN_SPANS); + while level <= MAX_GENERAL_SPANS { + let spans = [0, 1].map(|d| ((extent[d] / longer * level as f64).ceil() as usize).clamp(1, level)); + let holds = (0..2).all(|d| extent[d] / spans[d] as f64 >= 2.0 * step[d]); + if !holds && level > start { + return Ok(Round::Denser(best, level)); + } + let grid = Grid::new(spans, lo, hi); + let cell = (grid.cell(0) * grid.cell(1)).sqrt(); + let near = Near::of(&dev.d, &corners, &rough, &clashes, CLASH_CELLS * cell); + let (net, located) = fit_grid(&grid, st, &dev, spec, &near, &fade)?; + let got = measure(&net, st, &dev, spec, &located, &near, &fade)?; + best = best.min(got.shortfall); + let within = got.shortfall.deviation <= spec.tol && got.shortfall.angle <= spec.angle_tol && got.shortfall.points <= spec.tol; + if within && holds { + check_folds(&net, &dev.d, lo, hi)?; + let ctrl: Vec = net.ctrl.iter().map(|p| HPoint::weighted(*p, 1.0)).collect(); + let surface = NurbsSurface::from_homogeneous([3, 3], grid.knots.clone(), ctrl)?; + let surface = crown(surface, &net, st, &dev, spec.crown)?; + return Ok(Round::Done(Box::new(GeneralFill { + surface, + counter_clockwise: dev.counter_clockwise, + deviation: got.shortfall.deviation, + normal_deviation: got.shortfall.angle, + normal_deviation_all: got.angle_all, + curvature_mismatch: got.curvature, + point_deviation: got.shortfall.points, + spans, + stations: m, + }))); + } + if !holds { + return Ok(Round::Denser(best, level)); + } + level *= 2; + } + Ok(Round::Exhausted(best)) +} + +/// Which stations are close to a corner, and to a corner whose neighbours disagree, on one grid. +/// +/// At a corner the two sides' conditions meet and overdetermine the surface (the compatibility +/// of Hermite boundary data at a vertex): there one side's cross-boundary derivatives are the +/// other side's boundary curve, and the development's corner angle only approximates the loop's. +/// Close to any corner the second-derivative rows are left out, and close to a corner the domain +/// does not turn by the loop's own angle (a conformal map keeps angles), or where the neighbours' +/// tangent planes disagree, the first-derivative rows too: the positions decide there. The +/// neighbourhood is a few grid cells, so it shrinks as the grid refines. (Where the neighbours' +/// tangent planes disagree no surface is tangent to both; see [`fades`].) +struct Near { + corner: Vec, + rough: Vec, + clash: Vec, +} + +impl Near { + fn of(d: &[[f64; 2]], corners: &[[f64; 2]], rough: &[[f64; 2]], clashes: &[[f64; 2]], reach: f64) -> Near { + let close = |set: &[[f64; 2]]| -> Vec { d.iter().map(|q| set.iter().any(|c| (c[0] - q[0]).hypot(c[1] - q[1]) <= reach)).collect() }; + Near { corner: close(corners), rough: close(rough), clash: close(clashes) } + } + + fn rough(&self, k: usize) -> bool { + self.rough.get(k).copied().unwrap_or(true) + } + + fn corner(&self, k: usize) -> bool { + self.corner.get(k).copied().unwrap_or(true) + } + + fn clash(&self, k: usize) -> bool { + self.clash.get(k).copied().unwrap_or(true) + } +} + +/// Whether each station's G1 rows are asked (1) or not (0): not near a corner whose neighbours' +/// tangent planes disagree. +struct Fade { + weight: Vec, +} + +impl Fade { + fn at(&self, k: usize) -> f64 { + self.weight.get(k).copied().unwrap_or(0.0) + } +} + +/// The G1 rows' weights. At a corner whose neighbours' tangent planes disagree no surface is +/// tangent to both: its tangent plane there must hold both boundary tangents, so it is singular +/// at the corner however it is shaped. Within [`BLEND_SHARE`] of the run to the next corner of +/// such a corner the G1 rows are left out and the positions and the smoothness decide the +/// surface; those stations are not held to the angle tolerance. (Rows fading out towards the +/// corner instead keep asking for a tangent plane the corner cannot have, and the fit folds +/// there.) +fn fades(st: &[Station]) -> Fade { + let m = st.len(); + let mut weight = vec![1.0f64; m]; + let clashes: Vec = (0..m).filter(|k| st.get(*k).is_some_and(|s| s.clash)).collect(); + if clashes.is_empty() || m < 3 { + return Fade { weight }; + } + let p = |k: usize| st.get(k % m).map_or(Vec3::ZERO, |s| s.p); + let step = |k: usize| p(k + 1).dist(p(k)); + let is_corner = |k: usize| st.get(k % m).is_some_and(|s| s.corner || s.clash); + for &c in &clashes { + if let Some(w) = weight.get_mut(c) { + *w = 0.0; + } + // Back along the incoming side, then on along the outgoing one. + for back in [true, false] { + let next = |k: usize| if back { (k + m - 1) % m } else { (k + 1) % m }; + // The run to the next corner. + let (mut run, mut k) = (0.0, c); + for _ in 0..m { + let n = next(k); + run += if back { step(n) } else { step(k) }; + k = n; + if is_corner(k) { + break; + } + } + let reach = BLEND_SHARE * run; + let (mut dist, mut k) = (0.0, c); + for _ in 0..m { + let n = next(k); + dist += if back { step(n) } else { step(k) }; + k = n; + if dist >= reach || is_corner(k) { + break; + } + if let Some(w) = weight.get_mut(k) { + *w = 0.0; + } + } + } + } + Fade { weight } +} + +/// Fit the grid to the even stations' conditions and the points: the net, and where the points +/// ended up in the domain. +fn fit_grid<'a>( + grid: &'a Grid, + st: &[Station], + dev: &Developed, + spec: &FillSpec, + near: &Near, + fade: &Fade, +) -> SplineResult<(Net<'a>, Vec<[f64; 2]>)> { + let cell = (grid.cell(0) * grid.cell(1)).sqrt(); + let mut rows: Vec = Vec::new(); + for (k, s) in st.iter().enumerate().step_by(2) { + let Some(across) = s.across.filter(|_| !s.corner) else { continue }; + let (Some(&q), Some(&nu), Some(&speed)) = (dev.d.get(k), dev.nu.get(k), dev.speed.get(k)) else { continue }; + let weight = fade.at(k); + if !(weight > 1e-3) { + continue; + } + let full = weight >= 1.0; + let bend = across.bend.filter(|_| full && !near.corner(k)); + rows.push(Hermite { + q, + nu, + normal: across.normal, + leave: across.leave, + hermite: full && !near.rough(k), + weight, + curvature: bend.map(|b| across.normal.dot(b)), + speed, + }); + } + // The fitted stations (every other one), then the points where they were last found. + let fitted: Vec<([f64; 2], Vec3)> = st.iter().zip(&dev.d).step_by(2).map(|(s, q)| (*q, s.p)).collect(); + let mut extra: Vec<([f64; 2], Vec3)> = Vec::new(); + let rounds = if spec.points.is_empty() { 0 } else { POINT_ROUNDS }; + let mut round = 0; + loop { + let data: Vec<([f64; 2], Vec3)> = fitted.iter().chain(&extra).copied().collect(); + let net = if rows.is_empty() { fit_positions(grid, &data)? } else { fit_coupled(grid, &data, &mut rows, cell)? }; + let located = spec.points.iter().map(|p| locate(&net, &dev.d, *p)).collect::>>()?; + if round >= rounds { + return Ok((net, located)); + } + extra = located.iter().copied().zip(spec.points.iter().copied()).collect(); + round += 1; + } +} + +/// The fit to positions alone (G0 everywhere): the three coordinates apart, sharing one matrix. +fn fit_positions<'a>(grid: &'a Grid, data: &[([f64; 2], Vec3)]) -> SplineResult> { + let [cu, cv] = grid.count; + let n = cu * cv; + let mut sys = System::new(n, 3 * cu + 3); + for (q, p) in data { + let b = Basis::at(grid, *q).ok_or_else(|| SplineError::Invalid("a station is off the grid".into()))?; + sys.row(&b.terms(cu, &[(0, 0, 1.0)]), *p); + } + let lambda = SMOOTHNESS * sys.rows as f64 / n as f64 * grid.cell(0) * grid.cell(1); + grid.add_energy(&mut sys.band, lambda); + let factor = sys.band.factor_rows()?; + let [x, y, z] = sys.rhs; + let (x, y, z) = (factor.substitute(x)?, factor.substitute(y)?, factor.substitute(z)?); + net_of(grid, x.iter().zip(&y).zip(&z).map(|((x, y), z)| Vec3::new(*x, *y, *z)).collect()) +} + +fn net_of(grid: &Grid, ctrl: Vec) -> SplineResult> { + if ctrl.iter().any(|p| !p.is_finite()) { + return invalid("the fit is not finite"); + } + Ok(Net { grid, ctrl }) +} + +/// One G1/G2 station's rows (see [`fit_coupled`]). +struct Hermite { + q: [f64; 2], + nu: [f64; 2], + /// The normal of the tangent plane asked for, and the direction the patch should leave in. + normal: Vec3, + leave: Vec3, + /// Whether the weak Hermite rows are asked (not close to a corner the development does not + /// keep the angle of, unless in a blend). + hermite: bool, + /// How much the G1 row weighs (see [`fades`]). + weight: f64, + /// G2: the neighbour's normal curvature along `leave`. + curvature: Option, + /// How fast the patch leaves: first the boundary's own speed, then what the fit gives. + speed: f64, +} + +/// How much the Hermite rows (the full first derivative `s·w`) weigh against the G1 rows: enough +/// to say which way and about how fast the patch leaves, little enough to give way where the +/// full data cannot all hold. +const HERMITE_WEIGHT: f64 = 0.05; +/// Rounds of updating the speed a G2 row's curvature is asked at. +const G2_ROUNDS: usize = 4; + +/// The fit with G1 and G2 rows: the three coordinates together, since the rows that hold the +/// continuity are geometric and couple them. +/// +/// - Positions, as for [`fit_positions`]. +/// - G1: `n·∂S/∂ν = 0` — the patch's cross-boundary derivative lies in the tangent plane asked +/// for, which with the boundary tangent in it makes the normals agree. That is all G1 asks. +/// - Weak Hermite rows `∂S/∂ν = s·w` (weight [`HERMITE_WEIGHT`]): which way the patch leaves +/// (into the patch, not folding back over the neighbour) and about how fast. A whole set of +/// full Hermite data is consistent only with a boundary map that is exactly conformal at the +/// boundary's own speed, which a curved patch's is not, so these give way. +/// - G2: `n·∂²S/∂ν² = (∂S/∂ν·w)²·κ` — the patch's normal curvature across the edge is the +/// neighbour's. The speed is the fit's own, so the row is solved again with the speed from the +/// last solution (the matrix stays, so each round is a back substitution). +/// +/// The unknowns are numbered with the grid's shorter direction fastest and the coordinates +/// innermost, which keeps the band narrow. +fn fit_coupled<'a>(grid: &'a Grid, data: &[([f64; 2], Vec3)], rows: &mut [Hermite], cell: f64) -> SplineResult> { + let [cu, cv] = grid.count; + let n = cu * cv; + let u_fast = cu <= cv; + let fast = if u_fast { cu } else { cv }; + let index = |i: usize, j: usize, c: usize| if u_fast { (j * cu + i) * 3 + c } else { (i * cv + j) * 3 + c }; + let mut band = Band::new(3 * n, 3 * (3 * fast + 3) + 2); + let mut rhs = vec![0.0; 3 * n]; + let mut count = 0usize; + // A row on the three coordinates: `Σ_c Σ_t coefficient_t · weight_c · control_{t,c} = target`. + let mut row = |band: &mut Band, rhs: &mut [f64], cells: &[(usize, usize, f64); 16], weights: [f64; 3], target: f64| { + let mut terms = [(0usize, 0.0f64); 48]; + for (c, w) in weights.iter().enumerate() { + for (t, &(i, j, k)) in cells.iter().enumerate() { + if let Some(slot) = terms.get_mut(c * 16 + t) { + *slot = (index(i, j, c), k * w); + } + } + } + for &(r, wr) in &terms { + if wr == 0.0 { + continue; + } + if let Some(slot) = rhs.get_mut(r) { + *slot += wr * target; + } + for &(s, ws) in &terms { + band.add(r, s, wr * ws); + } + } + count += 1; + }; + let axes = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]; + let coords = |p: Vec3| [p.x, p.y, p.z]; + for &(q, p) in data { + let b = Basis::at(grid, q).ok_or_else(|| SplineError::Invalid("a station is off the grid".into()))?; + let cells = b.cells(&[(0, 0, 1.0)]); + for (axis, t) in axes.iter().zip(coords(p)) { + row(&mut band, &mut rhs, &cells, *axis, t); + } + } + // The right-hand side so far; the G2 rows' targets are added per round. + let mut second: Vec<([(usize, usize, f64); 16], usize)> = Vec::new(); + for (h, r) in rows.iter().enumerate() { + let b = Basis::at(grid, r.q).ok_or_else(|| SplineError::Invalid("a station is off the grid".into()))?; + let first = b.cells(&along(r.nu, 1)).map(|(i, j, k)| (i, j, k * cell)); + let held = first.map(|(i, j, k)| (i, j, k * r.weight)); + row(&mut band, &mut rhs, &held, coords(r.normal), 0.0); + if r.hermite { + let w = HERMITE_WEIGHT * r.weight; + let weak = first.map(|(i, j, k)| (i, j, k * w)); + for (axis, t) in axes.iter().zip(coords(r.leave * (r.speed * cell * w))) { + row(&mut band, &mut rhs, &weak, *axis, t); + } + } + if r.curvature.is_some() { + let cells = b.cells(&along(r.nu, 2)).map(|(i, j, k)| (i, j, k * cell * cell)); + row(&mut band, &mut rhs, &cells, coords(r.normal), 0.0); + second.push((cells, h)); + } + } + let lambda = SMOOTHNESS * count as f64 / (3 * n) as f64 * grid.cell(0) * grid.cell(1); + grid.energy(lambda, &mut |(i, j), (k, l), e| { + for c in 0..3 { + band.add(index(i, j, c), index(k, l, c), e); + } + }); + let factor = band.factor_rows()?; + let rounds = if second.is_empty() { 1 } else { G2_ROUNDS }; + let mut net = None; + for _ in 0..rounds { + let mut b = rhs.clone(); + for (cells, h) in &second { + let Some(r) = rows.get(*h) else { continue }; + let target = r.curvature.unwrap_or(0.0) * r.speed * r.speed * cell * cell; + for (c, w) in coords(r.normal).iter().enumerate() { + for &(i, j, k) in cells { + if let Some(slot) = b.get_mut(index(i, j, c)) { + *slot += k * w * target; + } + } + } + } + let x = factor.substitute(b)?; + let ctrl: Vec = (0..cv) + .flat_map(|j| (0..cu).map(move |i| (i, j))) + .map(|(i, j)| { + Vec3::new( + x.get(index(i, j, 0)).copied().unwrap_or(f64::NAN), + x.get(index(i, j, 1)).copied().unwrap_or(f64::NAN), + x.get(index(i, j, 2)).copied().unwrap_or(f64::NAN), + ) + }) + .collect(); + let current = net_of(grid, ctrl)?; + for r in rows.iter_mut().filter(|r| r.curvature.is_some()) { + let (Some(su), Some(sv)) = (current.der(r.q, 1, 0), current.der(r.q, 0, 1)) else { + return invalid("the fit cannot be evaluated on the loop"); + }; + let ahead = (su * r.nu[0] + sv * r.nu[1]).dot(r.leave); + if ahead > 0.0 { + r.speed = ahead; + } + } + net = Some(current); + } + net.ok_or_else(|| SplineError::Invalid("the fit made no rounds".into())) +} + +/// Where in the domain the surface comes nearest `p`: the nearest of a grid of samples inside +/// the domain, then Gauss–Newton. A point nearest the surface outside the loop is refused. +fn locate(net: &Net, d: &[[f64; 2]], p: Vec3) -> SplineResult<[f64; 2]> { + let kn = &net.grid.knots; + let (Some(u0), Some(u1), Some(v0), Some(v1)) = (kn[0].first(), kn[0].last(), kn[1].first(), kn[1].last()) else { + return invalid("the grid has no knots"); + }; + let (lo, hi) = ([*u0, *v0], [*u1, *v1]); + let samples = 48; + let mut best: Option<(f64, [f64; 2])> = None; + for i in 0..=samples { + for j in 0..=samples { + let q = [lo[0] + (hi[0] - lo[0]) * i as f64 / samples as f64, lo[1] + (hi[1] - lo[1]) * j as f64 / samples as f64]; + if !inside_loop(q, d) { + continue; + } + let Some(s) = net.point(q) else { continue }; + let dist = s.dist(p); + if best.is_none_or(|(b, _)| dist < b) { + best = Some((dist, q)); + } + } + } + let Some((_, mut q)) = best else { return Err(SplineError::Unsupported("passing through points: the loop's development is too thin".into())) }; + for _ in 0..40 { + let (Some(s), Some(su), Some(sv)) = (net.point(q), net.der(q, 1, 0), net.der(q, 0, 1)) else { break }; + let r = s - p; + let (a, b, c) = (su.dot(su), su.dot(sv), sv.dot(sv)); + let det = a * c - b * b; + if !(det.abs() > 1e-300) { + break; + } + let (ga, gb) = (su.dot(r), sv.dot(r)); + let next = [(q[0] - (c * ga - b * gb) / det).clamp(lo[0], hi[0]), (q[1] - (a * gb - b * ga) / det).clamp(lo[1], hi[1])]; + let moved = (next[0] - q[0]).hypot(next[1] - q[1]); + q = next; + if moved < 1e-13 * (hi[0] - lo[0]) { + break; + } + } + if !inside_loop(q, d) { + return Err(SplineError::Unsupported("passing through a point that is not over the inside of the loop".into())); + } + Ok(q) +} + +/// What a fit achieved. +struct Measured { + shortfall: Shortfall, + angle_all: f64, + curvature: f64, +} + +/// Measure a fit against every station and the points. +fn measure(net: &Net, st: &[Station], dev: &Developed, spec: &FillSpec, located: &[[f64; 2]], near: &Near, fade: &Fade) -> SplineResult { + let off = |what: &str| SplineError::Invalid(format!("the fit cannot be evaluated {what}")); + let sense = if dev.counter_clockwise { 1.0 } else { -1.0 }; + let (mut deviation, mut angle, mut angle_all, mut curvature) = (0.0f64, 0.0f64, 0.0f64, 0.0f64); + for (k, s) in st.iter().enumerate() { + let (Some(&q), Some(&nu)) = (dev.d.get(k), dev.nu.get(k)) else { continue }; + deviation = deviation.max(net.point(q).ok_or_else(|| off("on the loop"))?.dist(s.p)); + let Some(across) = s.across.filter(|_| !s.corner) else { continue }; + let (normal, su, sv) = net.frame(q).ok_or_else(|| off("on the loop"))?; + // The neighbour's normal turned to the side the patch's normal should be on. + let want = if across.normal.dot(s.up) * sense >= 0.0 { across.normal } else { across.normal * -1.0 }; + let a = normal.dot(want).clamp(-1.0, 1.0).acos(); + angle_all = angle_all.max(a); + if near.clash(k) || fade.at(k) < 1.0 { + continue; + } + angle = angle.max(a); + if let Some(bend) = across.bend.filter(|_| !near.corner(k)) { + let first = su * nu[0] + sv * nu[1]; + let second = net.der(q, 2, 0).zip(net.der(q, 1, 1)).zip(net.der(q, 0, 2)).ok_or_else(|| off("on the loop"))?; + let ((uu, uv), vv) = second; + let second = uu * (nu[0] * nu[0]) + uv * (2.0 * nu[0] * nu[1]) + vv * (nu[1] * nu[1]); + let speed2 = first.dot(first); + if speed2 > 0.0 { + curvature = curvature.max((normal.dot(second) / speed2 - want.dot(bend)).abs()); + } + } + } + let mut points = 0.0f64; + for (q, p) in located.iter().zip(spec.points) { + points = points.max(net.point(*q).ok_or_else(|| off("at a point"))?.dist(*p)); + } + if ![deviation, angle, angle_all, curvature, points].iter().all(|x| x.is_finite()) { + return invalid("the fit is not finite"); + } + Ok(Measured { shortfall: Shortfall { deviation, angle, points }, angle_all, curvature }) +} + +/// The fill with its crown (see [`FillSpec::crown`]). +fn crown(surface: NurbsSurface, net: &Net, st: &[Station], dev: &Developed, crown: f64) -> SplineResult { + if crown == 0.0 { + return Ok(surface); + } + // The side the neighbours' normals face, by the patch's normals at the G1/G2 stations. + let mut side = 0.0; + for (s, q) in st.iter().zip(&dev.d) { + if let (Some(a), Some((n, _, _))) = (s.across, net.frame(*q)) { + side += n.dot(a.normal).signum(); + } + } + if side == 0.0 { + return Err(SplineError::Unsupported("crowning a patch with no g1 or g2 edge to say which side is out".into())); + } + let sign = side.signum(); + let along = |g: [f64; 2]| net.frame(g).map(|(n, _, _)| n * sign); + let (surface, _, _) = crate::fill::crown_fill_along(surface, &dev.d, crown, &along, MAX_GENERAL_SPANS)?; + Ok(surface) +} + +/// Refuse a fit that folds over or pinches inside the loop: on a grid of samples inside the +/// domain, the normal must be defined everywhere (the surface's parameter speed nowhere collapses +/// against the rest) and must not turn over between neighbouring samples. +fn check_folds(net: &Net, d: &[[f64; 2]], lo: [f64; 2], hi: [f64; 2]) -> SplineResult<()> { + let n = FOLD_SAMPLES; + let mut cells: Vec> = Vec::with_capacity((n + 1) * (n + 1)); + for j in 0..=n { + for i in 0..=n { + let q = [lo[0] + (hi[0] - lo[0]) * i as f64 / n as f64, lo[1] + (hi[1] - lo[1]) * j as f64 / n as f64]; + if !inside_loop(q, d) { + cells.push(None); + continue; + } + let (Some(su), Some(sv)) = (net.der(q, 1, 0), net.der(q, 0, 1)) else { return invalid("the fit cannot be evaluated inside") }; + let c = su.cross(sv); + cells.push(Some((c, c.len()))); + } + } + let mut sizes: Vec = cells.iter().flatten().map(|(_, l)| *l).collect(); + if sizes.is_empty() { + return Ok(()); + } + sizes.sort_by(f64::total_cmp); + let median = sizes.get(sizes.len() / 2).copied().unwrap_or(0.0); + let smallest = sizes.first().copied().unwrap_or(0.0); + if !(smallest > 1e-3 * median) { + return Err(SplineError::Unsupported("filling this loop: the surface pinches inside it".into())); + } + let at = |i: usize, j: usize| cells.get(j * (n + 1) + i).copied().flatten(); + for j in 0..=n { + for i in 0..=n { + let Some((a, _)) = at(i, j) else { continue }; + for (x, y) in [(i + 1, j), (i, j + 1)] { + if let Some((b, _)) = at(x, y) + && a.dot(b) <= 0.0 + { + return Err(SplineError::Unsupported("filling this loop: the surface folds over inside it".into())); + } + } + } + } + Ok(()) +} + +#[cfg(test)] +#[path = "fill_general_tests.rs"] +mod tests; diff --git a/crates/geom/src/fill_general_tests.rs b/crates/geom/src/fill_general_tests.rs new file mode 100644 index 00000000..35ff6b6c --- /dev/null +++ b/crates/geom/src/fill_general_tests.rs @@ -0,0 +1,191 @@ +use super::*; +use std::f64::consts::{PI, TAU}; + +const TOL: f64 = 1e-5; +const ANGLE: f64 = 0.05 * PI / 180.0; + +/// The surface at fractions `(x, y)` of its domain. +fn at(s: &NurbsSurface, x: f64, y: f64) -> Vec3 { + let [du, dv] = s.domain(); + s.eval((du[0] + (du[1] - du[0]) * x).min(du[1]), (dv[0] + (dv[1] - dv[0]) * y).min(dv[1])).unwrap() +} + +fn spec(points: &[Vec3]) -> FillSpec<'_> { + FillSpec { tol: TOL, angle_tol: ANGLE, crown: 0.0, points } +} + +/// The radial and tangential unit vectors of a cylinder about z at angle `a`. +fn radial(a: f64) -> Vec3 { + Vec3::new(a.cos(), a.sin(), 0.0) +} +fn around(a: f64) -> Vec3 { + Vec3::new(-a.sin(), a.cos(), 0.0) +} + +/// The open end of a tube about z (radius `r`, its wall below z = 0), one side: G1 (or G2) to +/// the wall, which leaves straight up. +fn tube_end(r: f64, g2: bool) -> impl FnMut(usize, &[f64]) -> Result, String> { + move |_, ts| { + Ok(ts + .iter() + .map(|t| { + let a = t * TAU; + let across = Across { leave: Vec3::Z, normal: radial(a), bend: g2.then_some(Vec3::ZERO) }; + Station { p: radial(a) * r, up: around(a).cross(Vec3::Z), across: Some(across), corner: false, clash: false } + }) + .collect()) + } +} + +/// A window in the wall of a cylinder about z (radius `r`): `a` either side of the x axis, `h` +/// high, its four sides G1 to the wall around it (the same cylinder). +fn window(r: f64, a: f64, h: f64) -> impl FnMut(usize, &[f64]) -> Result, String> { + move |k, ts| { + Ok(ts + .iter() + .map(|&t| { + let (angle, z, tangent, leave) = match k { + 0 => (-a + 2.0 * a * t, 0.0, around(-a + 2.0 * a * t), Vec3::Z), + 1 => (a, h * t, Vec3::Z, around(a) * -1.0), + 2 => (a - 2.0 * a * t, h, around(a - 2.0 * a * t) * -1.0, Vec3::Z * -1.0), + _ => (-a, h - h * t, Vec3::Z * -1.0, around(-a)), + }; + let p = radial(angle) * r + Vec3::Z * z; + let up = tangent.cross(leave); + Station { p, up, across: Some(Across { leave, normal: radial(angle), bend: None }), corner: t == 0.0, clash: false } + }) + .collect()) + } +} + +#[test] +fn a_tube_end_capped_g1_meets_the_wall_and_rises_over_it() { + for g2 in [false, true] { + let f = fill_general(1, &mut tube_end(10.0, g2), &spec(&[])).unwrap(); + eprintln!( + "tube end (g2 {g2}): {:?} spans, {} stations, off the rim {:e}, normals within {:.4}°, curvature mismatch {:e}", + f.spans, + f.stations, + f.deviation, + f.normal_deviation.to_degrees(), + f.curvature_mismatch + ); + assert!(f.deviation <= TOL && f.normal_deviation <= ANGLE); + // The middle rises over the rim (the cap leaves the wall upwards). + let mid = at(&f.surface, 0.5, 0.5); + eprintln!("the cap's middle: {mid:?}"); + assert!(mid.z > 2.0 && mid.x.hypot(mid.y) < 1.0, "{mid:?}"); + } +} + +#[test] +fn a_window_in_a_cylinder_wrapping_past_half_way_fills_close_to_the_cylinder() { + // 240° around: no plane sees it as a height field. + let (r, a, h) = (10.0, 2.0 * PI / 3.0, 10.0); + let f = fill_general(4, &mut window(r, a, h), &spec(&[])).unwrap(); + eprintln!("window: {:?} spans, off {:e}, normals within {:.4}°", f.spans, f.deviation, f.normal_deviation.to_degrees()); + assert!(f.deviation <= TOL && f.normal_deviation <= ANGLE); + // Inside, the fill stays near the cylinder it was cut from (it is the least-bending surface + // with the cylinder's boundary data, not the cylinder itself). + let mut worst = 0.0f64; + for k in 1..20 { + for j in 1..20 { + let (ang, z) = (-a + 2.0 * a * k as f64 / 20.0, h * j as f64 / 20.0); + // The fill's point nearest this cylinder point, by a coarse search. + let target = radial(ang) * r + Vec3::Z * z; + let mut best = f64::INFINITY; + for u in 0..=120 { + for v in 0..=120 { + best = best.min(at(&f.surface, u as f64 / 120.0, v as f64 / 120.0).dist(target)); + } + } + worst = worst.max(best); + } + } + eprintln!("window: the cylinder is within {worst} of the fill"); + assert!(worst < 0.25, "{worst}"); +} + +#[test] +fn a_planar_polygon_is_its_plane_and_a_point_lifts_it() { + // An L-shaped hexagon (one reflex corner) in z = 0, its neighbours standing up around it. + let c = [[0.0, 0.0], [20.0, 0.0], [20.0, 10.0], [10.0, 10.0], [10.0, 20.0], [0.0, 20.0]]; + let mut stations = |k: usize, ts: &[f64]| -> Result, String> { + let (Some(a), Some(b)) = (c.get(k), c.get((k + 1) % 6)) else { return Err("no side".into()) }; + let tangent = Vec3::new(b[0] - a[0], b[1] - a[1], 0.0).normalized().ok_or("no tangent")?; + let inward = Vec3::Z.cross(tangent); + Ok(ts + .iter() + .map(|t| { + let p = Vec3::new(a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t, 0.0); + Station { p, up: Vec3::Z, across: Some(Across { leave: inward, normal: Vec3::Z, bend: None }), corner: *t == 0.0, clash: false } + }) + .collect()) + }; + let f = fill_general(6, &mut stations, &spec(&[])).unwrap(); + let mut off = 0.0f64; + for i in 0..=20 { + for j in 0..=20 { + off = off.max(at(&f.surface, i as f64 / 20.0, j as f64 / 20.0).z.abs()); + } + } + eprintln!("L: {:?} spans, the plane within {off:e}", f.spans); + assert!(off < 1e-8, "{off}"); + // G0 with a point 3 above the inside of the L: the surface passes through it. + let mut g0 = |k: usize, ts: &[f64]| stations(k, ts).map(|v| v.into_iter().map(|s| Station { across: None, ..s }).collect()); + let point = [Vec3::new(5.0, 5.0, 3.0)]; + let f = fill_general(6, &mut g0, &spec(&point)).unwrap(); + eprintln!("L through a point: {:?} spans, off the loop {:e}, off the point {:e}", f.spans, f.deviation, f.point_deviation); + assert!(f.deviation <= TOL && f.point_deviation <= TOL); +} + +#[test] +fn hostile_loops_and_specs_are_refused() { + let ok = |s: &FillSpec| fill_general(1, &mut tube_end(10.0, false), s); + for (tol, angle, crown) in + [(0.0, ANGLE, 0.0), (f64::NAN, ANGLE, 0.0), (TOL, -1.0, 0.0), (TOL, ANGLE, f64::NAN), (TOL, ANGLE, -0.95), (TOL, ANGLE, 1e300)] + { + assert!(ok(&FillSpec { tol, angle_tol: angle, crown, points: &[] }).is_err(), "{tol} {angle} {crown}"); + } + assert!(ok(&spec(&[Vec3::new(f64::NAN, 0.0, 0.0)])).is_err()); + assert!(ok(&spec(&vec![Vec3::ZERO; MAX_FILL_POINTS + 1])).is_err()); + // A point not over the inside of the loop. + assert!(ok(&spec(&[Vec3::new(30.0, 0.0, 0.0)])).is_err()); + assert!(fill_general(0, &mut tube_end(10.0, false), &spec(&[])).is_err()); + assert!(fill_general(10_000, &mut tube_end(10.0, false), &spec(&[])).is_err()); + assert!(fill_general(1, &mut |_, _| Err("no".into()), &spec(&[])).is_err()); + assert!( + fill_general( + 1, + &mut |_, ts| Ok(vec![Station { p: Vec3::ZERO, up: Vec3::Z, across: None, corner: false, clash: false }; ts.len()]), + &spec(&[]) + ) + .is_err() + ); + // A loop that runs out and straight back. + let mut spike = |_: usize, ts: &[f64]| -> Result, String> { + Ok(ts + .iter() + .map(|t| Station { p: Vec3::X * (1.0 - (2.0 * t - 1.0f64).abs()) * 10.0, up: Vec3::Z, across: None, corner: false, clash: false }) + .collect()) + }; + assert!(fill_general(1, &mut spike, &spec(&[])).is_err()); + // A figure eight: its development crosses itself. + let mut eight = |_: usize, ts: &[f64]| -> Result, String> { + Ok(ts + .iter() + .map(|t| { + let a = t * TAU; + Station { p: Vec3::new(10.0 * a.sin(), 5.0 * (2.0 * a).sin(), 0.0), up: Vec3::Z, across: None, corner: false, clash: false } + }) + .collect()) + }; + assert!(fill_general(1, &mut eight, &spec(&[])).is_err()); + // Huge and tiny loops do not panic. + for scale in [1e300, 1e-300] { + let mut s = |_: usize, ts: &[f64]| -> Result, String> { + Ok(ts.iter().map(|t| Station { p: radial(t * TAU) * scale, up: Vec3::Z, across: None, corner: false, clash: false }).collect()) + }; + let _ = fill_general(1, &mut s, &spec(&[])); + } +} diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index 21e44948..850f471f 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -14,6 +14,7 @@ pub mod bspline; pub mod compare; pub mod env; pub mod fill; +pub mod fill_general; mod mesh; mod plane; mod profile; diff --git a/crates/kernel/src/bridge.rs b/crates/kernel/src/bridge.rs index 1b55c40a..23db3874 100644 --- a/crates/kernel/src/bridge.rs +++ b/crates/kernel/src/bridge.rs @@ -191,44 +191,115 @@ pub fn frames(face: &mt::Face, edges: &[mt::Edge], leave: Vec3, count: usize) -> let leave = leave.normalized().filter(|l| l.is_finite()).ok_or_else(|| KernelError::Invalid("the leave direction must be a finite vector".into()))?; let (pts, length) = chain_samples(edges, count)?; - let surf = face.surface(); - let (Some(range_u), Some(range_v)) = surf.try_range_tuple() else { - return Err(KernelError::Failed("bridge: the support surface is unbounded".into())); + let support = Support::of(face, length)?; + let mut walk = Walk::new(length); + let along = |n: Vec3| -> Result { + (leave - n * leave.dot(n)) + .normalized() + .ok_or_else(|| KernelError::Invalid("the leave direction is normal to the support surface: it has no tangent plane there".to_string())) }; - // A face's own edges may sit a fit tolerance off its surface (a cylinder's wall as a NURBS - // next to its exact ring), so the projection is only asked to be close, not exact. - let tol = (length * 1e-4).max(1e-6); - // Finite-difference step: a thousandth of the chain, so curvature over the chain is read - // locally but round-off stays far below the geometry. - let mut step = (length * 1e-3).clamp(1e-6, 0.5); let mut out = Vec::with_capacity(pts.len()); - // The station before this one, as the hint for the projection: it keeps the search on the - // same side of a periodic or split surface. - let mut hint: Option<(f64, f64)> = None; for p in &pts { - let project = |hint: Option<(f64, f64)>| surf.search_nearest_parameter(p3(*p), hint, 100); - let close = |uv: (f64, f64)| from_p3(surf.subs(uv.0, uv.1)).dist(*p) <= tol; + out.push(support.station(*p, &along, &mut walk, true)?.0); + } + Ok(out) +} + +/// The frames of `face` at stations on its own boundary, each leaving the face straight across +/// that boundary: `stations` are points on the boundary with the boundary's direction there *as +/// the face's own boundary runs* (its wires as [`mt::Face::boundaries`] gives them), and the leave +/// direction is that direction × the face's outward normal — in the face's tangent plane, at +/// right angles to the boundary, away from the face. `length` is the size of the boundary the +/// stations are on (it sets how closely a station must lie on the surface and the step the shape +/// is read with). Each frame comes with the face's unit normal there (as the face is oriented). +/// With `shape` false only the leave direction and the normal are read (`bend` and +/// `bend_change` are zero), which needs no room on the surface beside the station. +pub fn frames_across(face: &mt::Face, stations: &[(Vec3, Vec3)], length: f64, shape: bool) -> Result> { + if stations.len() > 1 << 16 || !(length.is_finite() && length > 0.0) { + return Err(KernelError::Invalid("frames across a boundary need at most 65536 stations on a boundary of finite size".into())); + } + let support = Support::of(face, length)?; + let outward = if face.orientation() { 1.0 } else { -1.0 }; + let mut walk = Walk::new(length); + let mut out = Vec::with_capacity(stations.len()); + for (p, tangent) in stations { + if !p.is_finite() || !tangent.is_finite() { + return Err(KernelError::Invalid("a boundary station must be finite".into())); + } + let away = |n: Vec3| -> Result { + tangent.cross(n * outward).normalized().ok_or_else(|| KernelError::Invalid("the boundary runs along the face's normal".to_string())) + }; + let (frame, normal) = support.station(*p, &away, &mut walk, shape)?; + out.push((frame, normal * outward)); + } + Ok(out) +} + +/// A support surface read at stations: the surface, its parameter box, and how close a station +/// must lie to it. +struct Support { + surf: mt::Surface, + range_u: (f64, f64), + range_v: (f64, f64), + tol: f64, +} + +/// What reading one station hands on to the next: the station's parameters (the hint for the +/// next projection, which keeps the search on the same side of a periodic or split surface) and +/// the finite-difference step that worked. +struct Walk { + hint: Option<(f64, f64)>, + step: f64, +} + +impl Walk { + fn new(length: f64) -> Walk { + // Finite-difference step: a thousandth of the chain, so curvature over the chain is read + // locally but round-off stays far below the geometry. + Walk { hint: None, step: (length * 1e-3).clamp(1e-6, 0.5) } + } +} + +impl Support { + fn of(face: &mt::Face, length: f64) -> Result { + let surf = face.surface(); + let (Some(range_u), Some(range_v)) = surf.try_range_tuple() else { + return Err(KernelError::Failed("bridge: the support surface is unbounded".into())); + }; + // A face's own edges may sit a fit tolerance off its surface (a cylinder's wall as a NURBS + // next to its exact ring), so the projection is only asked to be close, not exact. + let tol = (length * 1e-4).max(1e-6); + Ok(Support { surf, range_u, range_v, tol }) + } + + /// The frame at `p`, leaving in `leave(n)` (`n` the surface's unit normal there), and that + /// normal. Without `read_shape`, the frame's shape is left at zero. + fn station(&self, p: Vec3, leave: &dyn Fn(Vec3) -> Result, walk: &mut Walk, read_shape: bool) -> Result<(Frame, Vec3)> { + let (surf, range_u, range_v, tol) = (&self.surf, self.range_u, self.range_v, self.tol); + let project = |hint: Option<(f64, f64)>| surf.search_nearest_parameter(p3(p), hint, 100); + let close = |uv: (f64, f64)| from_p3(surf.subs(uv.0, uv.1)).dist(p) <= tol; // The station before this one is the hint where it lands close; a chain that wraps the // surface (or a periodic one) falls back to a free search. - let (u, v) = match project(hint).filter(|uv| close(*uv)).or_else(|| project(None).filter(|uv| close(*uv))) { + let (u, v) = match project(walk.hint).filter(|uv| close(*uv)).or_else(|| project(None).filter(|uv| close(*uv))) { Some(uv) => uv, None => return Err(KernelError::Invalid("a picked edge is not on the support surface".into())), }; - hint = Some((u, v)); + walk.hint = Some((u, v)); let at = from_p3(surf.subs(u, v)); - if at.dist(*p) > tol { + if at.dist(p) > tol { return Err(KernelError::Invalid("a picked edge is not on the support surface".into())); } let normal = surf.normal(u, v); let n = Vec3::new(normal.x, normal.y, normal.z) .normalized() .ok_or_else(|| KernelError::Failed("bridge: the support surface has no normal".into()))?; - let w = (leave - n * leave.dot(n)) - .normalized() - .ok_or_else(|| KernelError::Invalid("the leave direction is normal to the support surface: it has no tangent plane there".to_string()))?; + let w = leave(n)?; + if !read_shape { + return Ok((Frame { p, leave: w, bend: Vec3::ZERO, bend_change: Vec3::ZERO }, n)); + } // Parameter-space direction that leaves in `w` at unit 3D speed: S_u·du + S_v·dv = w. - let su = tangent(&surf, u, v, range_u, true); - let sv = tangent(&surf, u, v, range_v, false); + let su = tangent(surf, u, v, range_u, true); + let sv = tangent(surf, u, v, range_v, false); let (a, b, c) = (su.dot(su), su.dot(sv), sv.dot(sv)); let det = a * c - b * b; let scale = a.max(c).max(1e-30); @@ -260,7 +331,7 @@ pub fn frames(face: &mt::Face, edges: &[mt::Edge], leave: Vec3, count: usize) -> None } }; - let mut h = step; + let mut h = walk.step; let mut chosen = pick(h); for _ in 0..8 { if chosen.is_some() { @@ -299,10 +370,9 @@ pub fn frames(face: &mt::Face, edges: &[mt::Edge], leave: Vec3, count: usize) -> // against the leave direction: the curvature vector is the same either way, the change of // curvature turns round with the direction of travel. let bend_change = if t.dot(w) < 0.0 { -bend_change } else { bend_change }; - out.push(Frame { p: *p, leave: w, bend, bend_change }); - step = h; + walk.step = h; + Ok((Frame { p, leave: w, bend, bend_change }, n)) } - Ok(out) } /// A first derivative of the surface in one parameter, by central differences (one-sided where diff --git a/crates/kernel/src/fill.rs b/crates/kernel/src/fill.rs index 90e80d96..29ea2b2f 100644 --- a/crates/kernel/src/fill.rs +++ b/crates/kernel/src/fill.rs @@ -1,16 +1,259 @@ -//! Filling a closed loop of a body's edges that is not planar with one face: a B-spline height -//! field over the loop's best-fit plane (`solidcraft_geom::fill::fill_loop`), bounded by the -//! body's own edges, whole. Loops that do not bound a height field over a plane are refused. +//! Filling a closed loop of a body's edges that is not planar with one face, bounded by the +//! body's own edges, whole: a B-spline height field over the loop's best-fit plane +//! (`solidcraft_geom::fill::fill_loop`) when the patch only has to meet the edges (G0) and the +//! loop bounds one; otherwise the general fill (`solidcraft_geom::fill_general`), developed from +//! the loop itself, which also meets the faces beside the loop with G1 or G2 and passes through +//! given points. use solidcraft_geom::Vec3; use solidcraft_geom::bspline::SplineError; +use solidcraft_geom::fill_general::{Across, FillSpec, MAX_FILL_POINTS, Station}; use truck_modeling as mt; use crate::body::{Body, from_p3}; +use crate::bridge::{Continuity, frames_across}; use crate::plane_fit::EdgeSampler; use crate::record::{Rebuild, rebuilt_body}; use crate::{KernelError, Result}; +/// How far (degrees) a G1 or G2 patch's normal may turn from its neighbour's along the loop. +pub const PATCH_ANGLE_TOL_DEG: f64 = 0.05; + +/// A spline error as the patch reports it. +fn spline_error(e: SplineError) -> KernelError { + match e { + SplineError::Unsupported(m) => KernelError::Invalid(format!("not supported yet: {m}")), + other => KernelError::Invalid(format!("cannot fill the loop: {other}")), + } +} + +/// What a patch on a loop of edges must do besides filling it. +#[derive(Clone, Debug, PartialEq)] +pub struct PatchSpec { + /// How every edge's neighbour is met (G0, G1 or G2). + pub continuity: Continuity, + /// Edges met otherwise: a point on the edge and its continuity. + pub edges: Vec<(Vec3, Continuity)>, + /// Points the patch passes through. + pub points: Vec, + /// How full the patch's middle is, `crown` times the loop's size (−0.9…3): the bump of + /// [`patch_edges_crowned`], laid along the general fill's own normal there (to the side its + /// G1/G2 neighbours' normals face) on control points clear of the loop, so the loop and the + /// conditions across it do not move. A general fill needs a G1 or G2 edge to crown. + pub crown: f64, +} + +impl Default for PatchSpec { + fn default() -> PatchSpec { + PatchSpec { continuity: Continuity::G0, edges: Vec::new(), points: Vec::new(), crown: 0.0 } + } +} + +impl PatchSpec { + /// Whether this is the plain patch of [`patch_edges_crowned`]: G0 everywhere, no points. + pub fn is_plain(&self) -> bool { + self.continuity == Continuity::G0 && self.edges.iter().all(|(_, g)| *g == Continuity::G0) && self.points.is_empty() + } +} + +/// What a general fill achieved (all zero for a planar or height-field patch, which meet their +/// edges within [`face_tolerance`] and ask nothing else). +#[derive(Clone, Copy, Debug, Default, PartialEq)] +pub struct PatchReport { + /// Whether the general fill made the face. + pub general: bool, + /// The farthest loop station from the surface. + pub gap: f64, + /// The largest angle (degrees) between the patch's normal and a G1/G2 neighbour's, away from + /// corners where the neighbours themselves disagree, and anywhere. + pub normal_deviation_deg: f64, + pub normal_deviation_all_deg: f64, + /// The largest difference (1/mm) between the patch's normal curvature across a G2 edge and + /// the neighbour's. + pub curvature_mismatch: f64, + /// The farthest point from the surface. + pub point_gap: f64, +} + +/// [`crate::patch_edges`] meeting the faces beside the loop with the continuity `spec` asks for +/// and passing through its points; a plain spec is exactly [`patch_edges_crowned`]. Its history +/// is over `[body]`: the face is new, its edges the body's own. +pub fn patch_edges_with(body: &Body, at: &[Vec3], spec: &PatchSpec) -> Result { + patch_edges_report(body, at, spec).map(|(b, _)| b) +} + +/// [`patch_edges_with`] and what the fill achieved. +pub fn patch_edges_report(body: &Body, at: &[Vec3], spec: &PatchSpec) -> Result<(Body, PatchReport)> { + if at.is_empty() || at.len() > 10_000 { + return Err(KernelError::Invalid("pick the edges of a closed loop".into())); + } + // A range test also throws out NaN. + if !(solidcraft_geom::fill::MIN_CROWN..=solidcraft_geom::fill::MAX_CROWN).contains(&spec.crown) { + return Err(KernelError::Invalid("a patch's crown must be between -0.9 and 3".into())); + } + let g3 = |g: Continuity| g == Continuity::G3; + if g3(spec.continuity) || spec.edges.iter().any(|(_, g)| g3(*g)) { + return Err(KernelError::Invalid("not supported yet: a patch meets its neighbours with g0, g1 or g2".into())); + } + if spec.edges.len() > 10_000 || spec.edges.iter().any(|(p, _)| !p.is_finite()) { + return Err(KernelError::Invalid("pick at most 10000 edges, each by a finite point, for their continuity".into())); + } + if spec.points.len() > MAX_FILL_POINTS || spec.points.iter().any(|p| !p.is_finite()) { + return Err(KernelError::Invalid(format!("a patch passes through at most {MAX_FILL_POINTS} finite points"))); + } + if spec.is_plain() { + return patch_edges_crowned(body, at, spec.crown).map(|b| (b, PatchReport::default())); + } + body.require_brep("patch")?; + crate::guard("patch", || { + let chain = crate::surfaces::picked_loop(body, at)?; + let mut orders = vec![spec.continuity; chain.len()]; + for (p, g) in &spec.edges { + let e = crate::surfaces::edge_near(body, *p).ok_or_else(|| KernelError::Invalid(format!("no edge at {p:?}")))?; + let k = chain + .iter() + .position(|c| c.id() == e.id()) + .ok_or_else(|| KernelError::Invalid("a continuity pick is not on an edge of the loop".into()))?; + if let Some(slot) = orders.get_mut(k) { + *slot = *g; + } + } + if spec.crown != 0.0 && orders.iter().all(|g| *g == Continuity::G0) { + return Err(KernelError::Invalid("not supported yet: a crown on a patch through points with no g1 or g2 edge".into())); + } + general_fill(body, &chain, &orders, &spec.points, spec.crown) + }) +} + +/// One edge of the loop as the general fill samples it. +struct Side { + sampler: EdgeSampler, + start: Vec3, + /// The face beside the edge, and whether its boundary runs the edge the way the loop does. + face: mt::Face, + along: bool, + order: Continuity, + /// The loop has a corner at the edge's start, and the neighbours' tangent planes disagree + /// there. + corner: bool, + clash: bool, +} + +impl Side { + fn point(&mut self, t: f64) -> Option { + // The loop's corners are its vertices (an imported edge's curve may end a little off). + if t <= 0.0 { Some(self.start) } else { self.sampler.at(t) } + } + + /// The loop's unit tangent at `t`, by a central difference (one-sided at the ends). + fn tangent(&mut self, t: f64) -> Option { + const STEP: f64 = 1e-4; + let (a, b) = ((t - STEP).max(0.0), (t + STEP).min(1.0)); + let (p, q) = (self.sampler.at(a)?, self.sampler.at(b)?); + (q - p).normalized() + } + + /// The face's frames at points of this edge with the loop's tangents there. + fn frames(&self, stations: &[(Vec3, Vec3)], size: f64, shape: bool) -> Result> { + let along: Vec<(Vec3, Vec3)> = stations.iter().map(|(p, t)| (*p, if self.along { *t } else { *t * -1.0 })).collect(); + frames_across(&self.face, &along, size, shape) + } +} + +/// The general fill of the closed loop `chain` of `body`'s edges, edge `k` meeting the face beside +/// it with `orders[k]`, passing through `points`, its middle made fuller by `crown` (see +/// [`PatchSpec::crown`]). Call it under [`crate::guard`]. +pub(crate) fn general_fill(body: &Body, chain: &[mt::Edge], orders: &[Continuity], points: &[Vec3], crown: f64) -> Result<(Body, PatchReport)> { + let tol = face_tolerance(body); + let size = body.size().max(solidcraft_geom::LINEAR_TOL); + let mut sides: Vec = Vec::with_capacity(chain.len()); + for (e, order) in chain.iter().zip(orders) { + let faces: Vec<&mt::Face> = body.solid.face_iter().filter(|f| f.edge_iter().any(|x| x.id() == e.id())).collect(); + let face = match faces.as_slice() { + [f] => (*f).clone(), + [f, ..] if *order == Continuity::G0 => (*f).clone(), + [] => return Err(KernelError::Failed("patch: an edge of the loop bounds no face".into())), + _ => return Err(KernelError::Invalid("an edge of the loop lies between two faces: g1 and g2 need the one face beside it".into())), + }; + let along = face + .boundaries() + .iter() + .flat_map(|w| w.edge_iter()) + .find(|x| x.id() == e.id()) + .map(|x| x.orientation() == e.orientation()) + .ok_or_else(|| KernelError::Failed("patch: the face beside an edge does not hold it".into()))?; + let sampler = EdgeSampler::new(e).ok_or_else(|| KernelError::Invalid("an edge has no parameter range".into()))?; + sides.push(Side { sampler, start: from_p3(e.front().point()), face, along, order: *order, corner: false, clash: false }); + } + // The corners: where the loop's tangent turns, and where the faces beside it disagree. + let n = sides.len(); + for k in 0..n { + let prev = (k + n - 1) % n; + let (Some(t_in), Some(t_out)) = (sides.get_mut(prev).and_then(|s| s.tangent(1.0)), sides.get_mut(k).and_then(|s| s.tangent(0.0))) else { + return Err(KernelError::Invalid("an edge of the loop has no tangent".into())); + }; + let Some(vertex) = sides.get(k).map(|s| s.start) else { continue }; + let corner = t_in.dot(t_out) < CORNER_DEG.to_radians().cos(); + let normal = |s: Option<&Side>, t: Vec3| -> Result { + let s = s.ok_or_else(|| KernelError::Failed("patch: side".into()))?; + let f = s.frames(&[(vertex, t)], size, false)?; + f.first().map(|f| f.1).ok_or_else(|| KernelError::Failed("patch: frame".into())) + }; + let g1 = |s: Option<&Side>| s.is_some_and(|s| s.order != Continuity::G0); + let clash = if g1(sides.get(prev)) || g1(sides.get(k)) { + let (a, b) = (normal(sides.get(prev), t_in)?, normal(sides.get(k), t_out)?); + a.dot(b).abs() < (2.0 * PATCH_ANGLE_TOL_DEG).to_radians().cos() + } else { + false + }; + if let Some(s) = sides.get_mut(k) { + s.corner = corner; + s.clash = clash; + } + } + let mut stations = |k: usize, ts: &[f64]| -> std::result::Result, String> { + let side = sides.get_mut(k).ok_or("no such side")?; + let mut at = Vec::with_capacity(ts.len()); + for &t in ts { + let (Some(p), Some(tangent)) = (side.point(t), side.tangent(t)) else { return Err("an edge cannot be sampled".into()) }; + at.push((p, tangent)); + } + let frames = side.frames(&at, size, side.order == Continuity::G2).map_err(|e| e.to_string())?; + let mut out = Vec::with_capacity(ts.len()); + for ((&t, (p, tangent)), (frame, normal)) in ts.iter().zip(&at).zip(&frames) { + let up = tangent.cross(frame.leave).normalized().ok_or("the loop runs across its neighbour")?; + let across = (side.order != Continuity::G0).then(|| Across { + leave: frame.leave, + normal: *normal, + bend: (side.order == Continuity::G2).then_some(frame.bend), + }); + let first = t <= 0.0; + out.push(Station { p: *p, up, across, corner: first && side.corner, clash: first && side.clash }); + } + Ok(out) + }; + let spec = FillSpec { tol, angle_tol: PATCH_ANGLE_TOL_DEG.to_radians(), crown, points }; + let fill = solidcraft_geom::fill_general::fill_general(chain.len(), &mut stations, &spec).map_err(spline_error)?; + // The wire runs counter-clockwise about the surface's normal. + let wire: mt::Wire = + if fill.counter_clockwise { chain.to_vec().into() } else { chain.iter().rev().map(mt::Edge::inverse).collect::>().into() }; + let face = mt::Face::try_new(vec![wire], crate::spline_conv::surface(&fill.surface)?).map_err(|e| KernelError::Failed(format!("fill: {e}")))?; + let mut rec = Rebuild::default(); + rec.made(&face, Vec::new()); + let report = PatchReport { + general: true, + gap: fill.deviation, + normal_deviation_deg: fill.normal_deviation.to_degrees(), + normal_deviation_all_deg: fill.normal_deviation_all.to_degrees(), + curvature_mismatch: fill.curvature_mismatch, + point_gap: fill.point_deviation, + }; + Ok((rebuilt_body(crate::surfaces::surface_body(vec![face]), &[&body.solid], &rec), report)) +} + +/// The loop turns by more than this (degrees) at a corner: its stations hold position only. +const CORNER_DEG: f64 = 1.0; + /// A surface body of one face filling the closed loop `chain` of `body`'s edges (each edge /// oriented along the loop, each ending where the next starts), within [`face_tolerance`] of /// every edge. Its history is over `[body]`: the face is new, its edges the body's own. Call it @@ -34,10 +277,16 @@ pub(crate) fn fill_loop_crowned(body: &Body, chain: &[mt::Edge], crown: f64) -> // The loop's corners are its vertices (an imported edge's curve may end a little off). if t <= 0.0 { Some(*start) } else { sampler.at(t) } }; - let fill = solidcraft_geom::fill::fill_loop(chain.len(), &mut point, tol, crown).map_err(|e| match e { - SplineError::Unsupported(m) => KernelError::Invalid(format!("not supported yet: {m}")), - other => KernelError::Invalid(format!("cannot fill the loop: {other}")), - })?; + let fill = match solidcraft_geom::fill::fill_loop(chain.len(), &mut point, tol, crown) { + Ok(fill) => fill, + // Not a height field over its plane (or too wavy for one): the general fill, developed + // from the loop itself, when no crown asks for the height field's bump. + Err(SplineError::Unsupported(_)) if crown == 0.0 => { + let orders = vec![Continuity::G0; chain.len()]; + return general_fill(body, chain, &orders, &[], 0.0).map(|(b, _)| b); + } + Err(e) => return Err(spline_error(e)), + }; // The wire runs counter-clockwise about the surface's normal. let wire: mt::Wire = if fill.counter_clockwise { chain.to_vec().into() } else { chain.iter().rev().map(mt::Edge::inverse).collect::>().into() }; @@ -489,9 +738,11 @@ mod tests { } #[test] - fn a_loop_that_does_not_bound_a_height_field_is_refused() { + fn a_loop_that_does_not_bound_a_height_field_fills_with_one_face() { // Half a round post's side: two straight edges up and two half circles, whose plane - // stands on the post's axis; the half circles run along its normal at their ends. + // stands on the post's axis; the half circles run along its normal at their ends, so no + // height field over it bounds them (the height-field fill refuses it). The general fill + // (#759) fills it with one face on its four edges. let post = cylinder(Vec3::ZERO, Vec3::Z, 10.0, 20.0).unwrap(); let halves = crate::split_by_plane(&post, &Plane::from_normal(Vec3::ZERO, Vec3::X).unwrap()).unwrap(); let half = halves.into_iter().find(|b| measure(b).unwrap().bbox.min.x > -1e-6).unwrap(); @@ -499,9 +750,16 @@ mod tests { (-4..=4).map(|k| Vec2::from_angle(k as f64 * 0.35).extend(10.0) * 1.0).map(|p| Vec3::new(p.x * 10.0, p.y * 10.0, 10.0)).collect(); let side = copy_faces(&half, &at).unwrap(); let rim = free_edge_points(&side); - let err = patch_edges(&side, &rim).unwrap_err().to_string(); - eprintln!("half a post's side: {err}"); - assert!(err.contains("not supported yet"), "{err}"); + let fill = patch_edges(&side, &rim).unwrap(); + assert_eq!(fill.face_count(), 1, "one face"); + let face = fill.solid.face_iter().next().unwrap().clone(); + let (off, tol) = (edges_off_surface(&face), super::face_tolerance(&side)); + assert!(off <= tol, "an edge is {off} off the face (tolerance {tol})"); + // Its area is at least its shadow on the plane x = 0 (the 20 × 20 square the loop bounds + // there) and, a least-bending fill, no more than the half side it spans. + let area = measure(&fill).unwrap().area; + eprintln!("half a post's side: one face of area {area}, its edges within {off:e}"); + assert!(area > 400.0 && area < std::f64::consts::PI * 10.0 * 20.0, "{area}"); } #[test] @@ -611,3 +869,7 @@ mod tests { } } } + +#[cfg(test)] +#[path = "fill_general_tests.rs"] +mod general_tests; diff --git a/crates/kernel/src/fill_general_tests.rs b/crates/kernel/src/fill_general_tests.rs new file mode 100644 index 00000000..3ac620ca --- /dev/null +++ b/crates/kernel/src/fill_general_tests.rs @@ -0,0 +1,506 @@ +//! The general fill on bodies (#759): G1 and G2 to the faces beside the loop, loops that are not +//! height fields, N-sided loops and points the patch passes through. Every number is checked +//! against the analytic surface the loop was cut from, or against a quadrature of the fitted +//! surface itself, never against the mesher alone. + +use mt::{BoundedCurve, InnerSpace, ParametricCurve, ParametricSurface, ParametricSurface3D}; +use solidcraft_geom::{Loop2, Plane, Region2, Vec2, Vec3}; +use truck_modeling as mt; + +use super::{PATCH_ANGLE_TOL_DEG, PatchReport, PatchSpec, face_tolerance, patch_edges_report}; +use crate::body::{Body, from_p3, p3}; +use crate::bridge::Continuity; +use crate::{BoolOp, boolean, box_solid, copy_faces, cylinder, extrude, measure, patch_edges, stitch, torus}; + +/// The surface's parameters nearest `p`: the nearest of a grid of samples (or `hint`), then +/// Gauss–Newton within the domain. +fn nearest(s: &mt::Surface, p: Vec3, hint: Option<(f64, f64)>) -> (f64, f64) { + let (Some((u0, u1)), Some((v0, v1))) = s.try_range_tuple() else { panic!("a bounded surface") }; + let p = p3(p); + let (mut u, mut v) = match hint { + Some(h) => h, + None => { + let n = 100; + let mut best = (f64::MAX, u0, v0); + for i in 0..=n { + for j in 0..=n { + let (u, v) = (u0 + (u1 - u0) * i as f64 / n as f64, v0 + (v1 - v0) * j as f64 / n as f64); + let d = (s.subs(u, v) - p).magnitude(); + if d < best.0 { + best = (d, u, v); + } + } + } + (best.1, best.2) + } + }; + for _ in 0..60 { + let (q, a, b) = (s.subs(u, v), s.uder(u, v), s.vder(u, v)); + let r = q - p; + let (aa, ab, bb) = (a.dot(a), a.dot(b), b.dot(b)); + let det = aa * bb - ab * ab; + if det.abs() < 1e-300 { + break; + } + let (ga, gb) = (a.dot(r), b.dot(r)); + u = (u - (bb * ga - ab * gb) / det).clamp(u0, u1); + v = (v - (aa * gb - ab * ga) / det).clamp(v0, v1); + } + (u, v) +} + +fn v3(v: mt::Vector3) -> Vec3 { + Vec3::new(v.x, v.y, v.z) +} + +/// Points along every edge of the loop: `per` on each. +fn rim_points(face: &mt::Face, per: usize) -> Vec { + let mut out = Vec::new(); + for e in face.boundaries().iter().flat_map(|w| w.edge_iter()) { + let c = e.oriented_curve(); + let (t0, t1) = c.range_tuple(); + for k in 0..per { + out.push(from_p3(c.subs(t0 + (t1 - t0) * (k as f64 + 0.5) / per as f64))); + } + } + out +} + +/// Along the loop, densely: the farthest a rim point is from the patch, the largest angle +/// (degrees) between the patch's normal and the neighbour's exact one (`normal(p)`, one sign for +/// the whole loop: a fold would show as 180°), and the largest difference between the patch's +/// normal curvature across the loop and `curvature(p)`, the neighbour's. +fn along_the_loop(fill: &Body, normal: &dyn Fn(Vec3) -> Vec3, curvature: &dyn Fn(Vec3) -> f64) -> (f64, f64, f64) { + let face = fill.solid.face_iter().next().unwrap().clone(); + let s = face.surface(); + let pts = rim_points(&face, 200); + let mut found = Vec::with_capacity(pts.len()); + let mut hint = None; + for p in &pts { + let uv = nearest(&s, *p, hint); + hint = Some(uv); + found.push(uv); + } + // One sign for the whole loop: the one most rim points agree with. + let patch_normal = |(u, v): (f64, f64)| v3(s.normal(u, v)); + let agree: f64 = pts.iter().zip(&found).map(|(p, uv)| patch_normal(*uv).dot(normal(*p)).signum()).sum(); + let sign = agree.signum(); + let (mut gap, mut angle, mut kappa) = (0.0f64, 0.0f64, 0.0f64); + for (k, (p, &(u, v))) in pts.iter().zip(&found).enumerate() { + gap = gap.max(from_p3(s.subs(u, v)).dist(*p)); + let n = patch_normal((u, v)); + let want = normal(*p) * sign; + angle = angle.max(n.dot(want).clamp(-1.0, 1.0).acos().to_degrees()); + // Across the loop: the tangent-plane direction square to the rim. + let next = pts.get(k + 1).or(pts.get(k.wrapping_sub(1))).copied().unwrap(); + let t = (next - *p).normalized().unwrap(); + let Some(d) = n.cross(t).normalized() else { continue }; + let (su, sv) = (v3(s.uder(u, v)), v3(s.vder(u, v))); + // d = a·su + b·sv. + let (aa, ab, bb) = (su.dot(su), su.dot(sv), sv.dot(sv)); + let det = aa * bb - ab * ab; + let (a, b) = ((bb * d.dot(su) - ab * d.dot(sv)) / det, (aa * d.dot(sv) - ab * d.dot(su)) / det); + // The second derivatives by central differences of the surface's own points. + let h = 1e-3; + let e = |du: f64, dv: f64| from_p3(s.subs(u + du * h, v + dv * h)); + let suu = (e(1.0, 0.0) - e(0.0, 0.0) * 2.0 + e(-1.0, 0.0)) * (1.0 / (h * h)); + let svv = (e(0.0, 1.0) - e(0.0, 0.0) * 2.0 + e(0.0, -1.0)) * (1.0 / (h * h)); + let suv = (e(1.0, 1.0) - e(1.0, -1.0) - e(-1.0, 1.0) + e(-1.0, -1.0)) * (0.25 / (h * h)); + let second = suu * (a * a) + suv * (2.0 * a * b) + svv * (b * b); + kappa = kappa.max((second.dot(want) - curvature(*p)).abs()); + } + (gap, angle, kappa) +} + +/// `∫∫ f(S, S_u × S_v) du dv` over the face's own region of its surface (its loop taken into +/// parameter space by the nearest parameters of points along it), by the midpoint rule on a +/// `n × n` grid over the region's box: the quadrature the mesher's area is checked against. +fn quadrature(face: &mt::Face, n: usize, f: &dyn Fn(Vec3, Vec3) -> f64) -> f64 { + let s = face.surface(); + let mut poly: Vec<(f64, f64)> = Vec::new(); + let mut hint = None; + for p in rim_points(face, 400) { + let uv = nearest(&s, p, hint); + hint = Some(uv); + poly.push(uv); + } + let (mut lo, mut hi) = ((f64::MAX, f64::MAX), (f64::MIN, f64::MIN)); + for &(u, v) in &poly { + lo = (lo.0.min(u), lo.1.min(v)); + hi = (hi.0.max(u), hi.1.max(v)); + } + let inside = |u: f64, v: f64| { + let mut w = false; + for k in 0..poly.len() { + let (a, b) = (poly[k], poly[(k + 1) % poly.len()]); + if (a.1 > v) != (b.1 > v) && u < a.0 + (v - a.1) / (b.1 - a.1) * (b.0 - a.0) { + w = !w; + } + } + w + }; + let (du, dv) = ((hi.0 - lo.0) / n as f64, (hi.1 - lo.1) / n as f64); + let mut sum = 0.0; + for i in 0..n { + for j in 0..n { + let (u, v) = (lo.0 + (i as f64 + 0.5) * du, lo.1 + (j as f64 + 0.5) * dv); + if inside(u, v) { + sum += f(from_p3(s.subs(u, v)), v3(s.uder(u, v)).cross(v3(s.vder(u, v)))) * du * dv; + } + } + } + sum +} + +/// The fill is one face on the open body's own edges, within the fill's tolerance of them; its +/// mesh has no folded triangle (each faces the way its corners' normals do); and the mesher's +/// area agrees with a quadrature of the surface over its own region. Returns the area. +fn check_face(name: &str, open: &Body, fill: &Body) -> f64 { + assert!(fill.is_surface()); + assert_eq!(fill.face_count(), 1, "{name}: one face"); + let own: Vec = open.solid.edge_iter().map(|e| e.id()).collect(); + assert!(fill.solid.edge_iter().all(|e| own.contains(&e.id())), "{name}: an edge that is not the body's"); + let face = fill.solid.face_iter().next().unwrap().clone(); + let s = face.surface(); + let tol = face_tolerance(open); + let mut off = 0.0f64; + let mut hint = None; + for p in rim_points(&face, 64) { + let uv = nearest(&s, p, hint); + hint = Some(uv); + off = off.max(from_p3(s.subs(uv.0, uv.1)).dist(p)); + } + assert!(off <= tol * 1.01, "{name}: an edge is {off} off the face (tolerance {tol})"); + let m = fill.tessellate(0.01).unwrap(); + assert!(!m.triangles.is_empty(), "{name}: the face meshes"); + for t in &m.triangles { + let [a, b, c] = m.tri(t).unwrap(); + let v = (b - a).cross(c - a); + let n: Vec3 = t.iter().map(|&i| m.normals[i as usize]).fold(Vec3::ZERO, |s, x| s + x); + if v.len() > 1e-10 { + assert!(v.dot(n) > 0.0, "{name}: a folded triangle"); + } + } + let area = measure(fill).unwrap().area; + let quad = quadrature(&face, 500, &|_, n| n.len()); + eprintln!("{name}: one face, edges within {off:e} (tolerance {tol:e}), mesh area {area:.4}, quadrature {quad:.4}"); + assert!((area - quad).abs() < 3e-3 * quad, "{name}: mesh area {area} against the surface's {quad}"); + area +} + +fn spec(g: Continuity) -> PatchSpec { + PatchSpec { continuity: g, ..PatchSpec::default() } +} + +fn report(name: &str, r: &PatchReport) { + eprintln!( + "{name}: fill reports gap {:e}, normals within {:.4}° ({:.4}° near clashing corners), curvature mismatch {:e}, points within {:e}", + r.gap, r.normal_deviation_deg, r.normal_deviation_all_deg, r.curvature_mismatch, r.point_gap + ); +} + +/// A cylinder about z of radius 10 and 30 high with its top face taken off; points around the rim. +fn open_tube() -> (Body, Vec) { + let solid = cylinder(Vec3::ZERO, Vec3::Z, 10.0, 30.0).unwrap(); + let open = copy_faces(&solid, &[&around(10.0, 15.0)[..], &[Vec3::ZERO]].concat()).unwrap(); + assert!(open.is_surface()); + let rim = (0..8).map(|k| Vec2::from_angle(k as f64 * 0.785).extend(0.0) * 10.0 + Vec3::Z * 30.0).collect(); + (open, rim) +} + +/// `body` without its planar faces but those at height `z0` or `z1` (the ends of a cylinder): +/// a cut solid's walls (and ends) as a surface, the faces the cut made taken away. +fn curved_and_ends(body: &Body, ends: &[f64]) -> Body { + let drop: Vec = body + .solid + .face_iter() + .enumerate() + .filter(|(_, f)| match f.surface() { + mt::Surface::Plane(p) => !ends.iter().any(|z| (p.origin().z - z).abs() < 1e-9 && p.normal().z.abs() > 1.0 - 1e-12), + _ => false, + }) + .map(|(i, _)| i) + .collect(); + let open = crate::remove_faces(body, &drop).unwrap(); + assert!(open.is_surface()); + open +} + +/// A point on each edge that bounds only one face of `body` (its free edges): the middle of +/// each, where the edge picker samples it. +fn free_edge_points(body: &Body) -> Vec { + let mut uses: std::collections::HashMap = std::collections::HashMap::new(); + for f in body.solid.face_iter() { + for e in f.boundaries().iter().flat_map(|w| w.edge_iter()) { + uses.entry(e.id()).or_insert((0, e.clone())).0 += 1; + } + } + let mut points = uses + .values() + .filter(|(n, _)| *n == 1) + .map(|(_, e)| { + let c = e.curve(); + let (t0, t1) = c.range_tuple(); + from_p3(c.subs(0.5 * (t0 + t1))) + }) + .collect::>(); + // In a fixed order, whatever the map's. + points.sort_by(|a, b| (a.z, a.y, a.x).partial_cmp(&(b.z, b.y, b.x)).unwrap()); + points +} + +/// Points round the circle of radius `r` about z at height `z`, every 30°. +fn around(r: f64, z: f64) -> Vec { + (0..12).map(|k| Vec2::from_angle(0.1 + k as f64 * std::f64::consts::TAU / 12.0).extend(0.0) * r + Vec3::Z * z).collect() +} + +fn radial(p: Vec3) -> Vec3 { + Vec3::new(p.x, p.y, 0.0).normalized().unwrap() +} + +#[test] +fn an_opened_tube_end_capped_g1_and_g2_meets_the_wall() { + let (open, rim) = open_tube(); + let mut kappas = Vec::new(); + for g in [Continuity::G1, Continuity::G2] { + let name = format!("tube end {}", g.as_str()); + let started = std::time::Instant::now(); + let (fill, r) = patch_edges_report(&open, &rim, &spec(g)).unwrap(); + eprintln!("{name}: filled in {:?}", started.elapsed()); + report(&name, &r); + assert!(r.general); + let area = check_face(&name, &open, &fill); + // The wall's normal is radial and it does not bend along its axis. + let (gap, angle, kappa) = along_the_loop(&fill, &radial, &|_| 0.0); + eprintln!( + "{name}: along the rim (1600 points): within {gap:e} of it, normal within {angle:.5}° of the wall's, curvature across {kappa:e} /mm (the wall's 0)" + ); + assert!(angle < 0.1, "{name}: the normal turns {angle}° from the wall's"); + kappas.push(kappa); + // A dome over the rim, not a flat lid: more than the disc. + assert!(area > std::f64::consts::PI * 100.0 * 1.05, "{name}: {area}"); + // Sealed: the tube's volume plus the dome's, the dome's by the divergence theorem over the + // patch (∫∫ (z − 30)·n_z dA), independent of the mesher. + let closed = stitch(&[&open, &fill], 0.0).unwrap(); + assert!(!closed.is_surface(), "{name}: the cap closes the tube"); + let face = fill.solid.face_iter().next().unwrap().clone(); + let dome = quadrature(&face, 500, &|p, n| (p.z - 30.0) * n.z).abs(); + let (got, want) = (measure(&closed).unwrap().volume, std::f64::consts::PI * 100.0 * 30.0 + dome); + eprintln!("{name}: sealed volume {got:.3}, tube + dome {want:.3}"); + assert!((got - want).abs() < 2e-4 * want, "{name}: {got} vs {want}"); + } + eprintln!("tube end: curvature across the rim g1 {:e}, g2 {:e} /mm", kappas[0], kappas[1]); + assert!(kappas[1] < kappas[0], "g2 holds the wall's curvature closer than g1"); +} + +/// A cylinder about z (radius 10, 40 high) with a window `half` either side of the x axis +/// through its wall, from z = 15 to 25; the wall and both ends without the faces the cut made (the +/// window's edges bound only the wall), and a point on each of the window's edges. +fn windowed(half: f64) -> (Body, Body, Vec) { + let solid = cylinder(Vec3::ZERO, Vec3::Z, 10.0, 40.0).unwrap(); + let x = 10.0 * half.cos(); + let tool = box_solid(Vec3::new(x, -11.0, 15.0), Vec3::new(11.0, 11.0, 25.0)).unwrap(); + let cut = boolean(&solid, &tool, BoolOp::Cut).unwrap().unwrap(); + let open = curved_and_ends(&cut, &[0.0, 40.0]); + let rim = free_edge_points(&open); + (solid, open, rim) +} + +#[test] +fn a_window_cut_through_a_cylinder_wall_past_half_way_fills_g1_close_to_the_cylinder() { + // 240° round: seen along any direction the window folds over itself, so it bounds no height + // field; the old fill refused it. + let half = 2.0 * std::f64::consts::PI / 3.0; + let (solid, open, rim) = windowed(half); + for g in [Continuity::G0, Continuity::G1] { + let name = format!("window {}", g.as_str()); + let (fill, r) = patch_edges_report(&open, &rim, &spec(g)).unwrap(); + report(&name, &r); + let area = check_face(&name, &open, &fill); + let (gap, angle, _) = along_the_loop(&fill, &radial, &|_| 0.0); + eprintln!("{name}: along the window (1600 points): within {gap:e}, normal within {angle:.5}° of the wall's"); + // How far the fill strays from the cylinder it replaces, inside. + let face = fill.solid.face_iter().next().unwrap().clone(); + let s = face.surface(); + let (Some((u0, u1)), Some((v0, v1))) = s.try_range_tuple() else { panic!() }; + let mut stray = 0.0f64; + for i in 0..=40 { + for j in 0..=40 { + let p = from_p3(s.subs(u0 + (u1 - u0) * i as f64 / 40.0, v0 + (v1 - v0) * j as f64 / 40.0)); + let ang = p.y.atan2(p.x); + if ang.abs() < half - 0.05 && p.z > 15.5 && p.z < 24.5 { + stray = stray.max((Vec2::new(p.x, p.y).len() - 10.0).abs()); + } + } + } + let exact = 10.0 * 2.0 * half * 10.0; + let closed = stitch(&[&open, &fill], 0.0).unwrap(); + let (got, want) = (measure(&closed).unwrap().volume, measure(&solid).unwrap().volume); + eprintln!( + "{name}: area {area:.3} (the cylinder's window {exact:.3}), the fill within {stray:.4} of the cylinder inside, sealed volume {got:.2} against the whole cylinder's {want:.2}" + ); + assert!(!closed.is_surface()); + if g == Continuity::G1 { + assert!(angle < 0.1, "{name}: the normal turns {angle}° from the wall's"); + assert!((area - exact).abs() < 0.01 * exact && stray < 0.25, "{name}: {area} vs {exact}, stray {stray}"); + assert!((got - want).abs() < 0.005 * want, "{name}: {got} vs {want}"); + } + } +} + +#[test] +fn a_saddle_window_in_a_torus_fills_g0_and_g1() { + // The inner wall of a torus (major 20, minor 8) is saddle-shaped; a window 240° round it + // between z = ±3 bounds no height field. + let solid = torus(Vec3::ZERO, 20.0, 8.0).unwrap(); + let tool = box_solid(Vec3::new(-6.0, -13.0, -3.0), Vec3::new(13.0, 13.0, 3.0)).unwrap(); + let cut = boolean(&solid, &tool, BoolOp::Cut).unwrap().unwrap(); + let open = curved_and_ends(&cut, &[]); + let rim = free_edge_points(&open); + // The torus's normal: from the tube's centre circle. + let normal = |p: Vec3| { + let c = Vec3::new(p.x, p.y, 0.0).normalized().unwrap() * 20.0; + (p - c).normalized().unwrap() + }; + for g in [Continuity::G0, Continuity::G1] { + let name = format!("torus window {}", g.as_str()); + let (fill, r) = patch_edges_report(&open, &rim, &spec(g)).unwrap(); + report(&name, &r); + check_face(&name, &open, &fill); + let (gap, angle, _) = along_the_loop(&fill, &normal, &|_| 0.0); + eprintln!("{name}: along the window (1600 points): within {gap:e}, normal within {angle:.5}° of the torus's"); + if g == Continuity::G1 { + assert!(angle < 0.1, "{name}: {angle}°"); + } + } +} + +/// A cylinder about z (radius 10, 40 high) with a pentagonal pocket cut into its wall from the +/// +x side (the pentagon in the yz plane about z = 20, circumradius 6): the wall and both ends +/// copied out, and a point on each of the pocket rim's five edges. +fn pentagon_window() -> (Body, Vec) { + let solid = cylinder(Vec3::ZERO, Vec3::Z, 10.0, 40.0).unwrap(); + let corners: Vec = (0..5).map(|k| Vec2::from_angle(0.3 + k as f64 * std::f64::consts::TAU / 5.0) * 6.0 + Vec2::new(0.0, 20.0)).collect(); + let prism = extrude(&Plane::YZ, &[Region2 { outer: Loop2::polygon(&corners), holes: vec![] }], 0.0, 15.0).unwrap().pop().unwrap(); + let cut = boolean(&solid, &prism, BoolOp::Cut).unwrap().unwrap(); + let open = curved_and_ends(&cut, &[0.0, 40.0]); + let rim = free_edge_points(&open); + (open, rim) +} + +#[test] +fn a_five_sided_window_fills_as_one_face_g1_and_through_a_point() { + let (open, rim) = pentagon_window(); + let (fill, r) = patch_edges_report(&open, &rim, &spec(Continuity::G1)).unwrap(); + report("pentagon g1", &r); + assert!(fill.edge_count() >= 5, "the rim's own edges (five, or more where a seam splits one)"); + check_face("pentagon g1", &open, &fill); + let (gap, angle, _) = along_the_loop(&fill, &radial, &|_| 0.0); + eprintln!("pentagon g1: along the rim: within {gap:e}, normal within {angle:.5}° of the wall's"); + assert!(angle < 0.1, "{angle}°"); + // Through a point half a millimetre out from the wall in the middle (G0 and G1). + let point = Vec3::new(10.5, 0.0, 20.0); + for g in [Continuity::G0, Continuity::G1] { + let name = format!("pentagon {} through a point", g.as_str()); + let spec = PatchSpec { continuity: g, points: vec![point], ..PatchSpec::default() }; + let (fill, r) = patch_edges_report(&open, &rim, &spec).unwrap(); + report(&name, &r); + check_face(&name, &open, &fill); + let s = fill.solid.face_iter().next().unwrap().surface(); + let (u, v) = nearest(&s, point, None); + let off = from_p3(s.subs(u, v)).dist(point); + eprintln!("{name}: the point is {off:e} off the surface"); + assert!(off <= face_tolerance(&open), "{name}: {off}"); + } +} + +#[test] +fn a_hexagonal_rim_fills_through_a_point_and_g1_to_walls_meeting_at_angles_holds_away_from_the_corners() { + // A hexagonal prism's top cut by a sideways cylinder: six walls meeting at 120°. + let corners: Vec = (0..6).map(|k| Vec2::from_angle(0.3 + k as f64 * std::f64::consts::TAU / 6.0) * 20.0).collect(); + let prism = extrude(&Plane::XY, &[Region2 { outer: Loop2::polygon(&corners), holes: vec![] }], 0.0, 20.0).unwrap().pop().unwrap(); + let tool = cylinder(Vec3::new(0.0, -60.0, 50.0), Vec3::Y, 40.0, 120.0).unwrap(); + let solid = boolean(&prism, &tool, BoolOp::Cut).unwrap().unwrap(); + let mut keep = vec![Vec3::ZERO]; + let mut rim = Vec::new(); + for k in 0..6 { + let m = corners[k].lerp(corners[(k + 1) % 6], 0.5); + keep.push(m.extend(5.0)); + rim.push(m.extend(50.0 - (1600.0 - m.x * m.x).sqrt())); + } + let open = copy_faces(&solid, &keep).unwrap(); + let plain = patch_edges(&open, &rim).unwrap(); + let sealed = |f: &Body| measure(&stitch(&[&open, f], 0.0).unwrap()).unwrap().volume; + // Six sides, G0, through a point 3 mm above the middle of the plain fill: one face. + let s = plain.solid.face_iter().next().unwrap().surface(); + let (u, v) = nearest(&s, Vec3::new(0.0, 0.0, 12.0), None); + let point = from_p3(s.subs(u, v)) + Vec3::Z * 3.0; + let (lifted, r) = patch_edges_report(&open, &rim, &PatchSpec { points: vec![point], ..PatchSpec::default() }).unwrap(); + report("hexagon through a point", &r); + assert_eq!(lifted.edge_count(), 6, "the rim's own six edges"); + check_face("hexagon through a point", &open, &lifted); + let s = lifted.solid.face_iter().next().unwrap().surface(); + let (u, v) = nearest(&s, point, None); + let off = from_p3(s.subs(u, v)).dist(point); + eprintln!("hexagon through a point: the point is {off:e} off the surface; sealed {:.2} against {:.2} plain", sealed(&lifted), sealed(&plain)); + assert!(off <= face_tolerance(&open) && sealed(&lifted) > sealed(&plain)); + // G1 to the walls: at a corner no surface is tangent to both walls (its tangent plane there + // holds both rim edges), so G1 is asked only away from the corners. On walls standing square + // to the patch the corner is pinched too sharply for the finest grid to hold the rim within + // the fill's tolerance, and the fill says so. + match patch_edges_report(&open, &rim, &spec(Continuity::G1)) { + Ok((fill, r)) => { + report("hexagon g1", &r); + check_face("hexagon g1", &open, &fill); + assert!(r.normal_deviation_deg <= PATCH_ANGLE_TOL_DEG); + assert!(sealed(&fill) > sealed(&plain)); + } + Err(e) => { + eprintln!("hexagon g1: {e}"); + assert!(e.to_string().contains("not supported yet"), "{e}"); + } + } +} + +#[test] +fn per_edge_continuity_and_hostile_specs() { + let (open, rim) = pentagon_window(); + // G1 on two edges only. + let spec = PatchSpec { edges: vec![(rim[0], Continuity::G1), (rim[2], Continuity::G2)], ..PatchSpec::default() }; + let (fill, r) = patch_edges_report(&open, &rim, &spec).unwrap(); + report("pentagon, two edges g1/g2", &r); + check_face("pentagon, two edges g1/g2", &open, &fill); + let bad = |spec: PatchSpec| patch_edges_report(&open, &rim, &spec).is_err(); + let g1 = PatchSpec { continuity: Continuity::G1, ..PatchSpec::default() }; + assert!(bad(PatchSpec { continuity: Continuity::G3, ..g1.clone() })); + assert!(bad(PatchSpec { edges: vec![(rim[0], Continuity::G3)], ..g1.clone() })); + assert!(bad(PatchSpec { edges: vec![(Vec3::new(0.0, 0.0, 40.0), Continuity::G1)], ..g1.clone() }), "an edge off the loop"); + assert!(bad(PatchSpec { edges: vec![(Vec3::new(f64::NAN, 0.0, 0.0), Continuity::G1)], ..g1.clone() })); + assert!(bad(PatchSpec { points: vec![Vec3::new(f64::INFINITY, 0.0, 0.0)], ..g1.clone() })); + assert!(bad(PatchSpec { points: vec![Vec3::new(10.0, 0.0, 20.0); 65], ..g1.clone() })); + assert!(bad(PatchSpec { points: vec![Vec3::new(-10.0, 0.0, 20.0)], ..g1.clone() }), "a point off the far side"); + for crown in [f64::NAN, 3.5, -1.0] { + assert!(bad(PatchSpec { crown, ..g1.clone() }), "crown {crown}"); + } + assert!(bad(PatchSpec { crown: 0.1, points: vec![Vec3::new(10.2, 0.0, 20.0)], ..PatchSpec::default() }), "a crown with no g1 edge"); + // A crown fills the middle out from the wall (the side the wall's normals face) and leaves + // the loop, and G1 across it, as they were. + let (plain, _) = patch_edges_report(&open, &rim, &g1).unwrap(); + let (full, r) = patch_edges_report(&open, &rim, &PatchSpec { crown: 0.1, ..g1.clone() }).unwrap(); + report("pentagon g1 crown 0.1", &r); + check_face("pentagon g1 crown 0.1", &open, &full); + let (_, angle, _) = along_the_loop(&full, &radial, &|_| 0.0); + let middle = |b: &Body| { + let s = b.solid.face_iter().next().unwrap().surface(); + let (u, v) = nearest(&s, Vec3::new(10.0, 0.0, 20.0), None); + let p = from_p3(s.subs(u, v)); + Vec2::new(p.x, p.y).len() + }; + let (a, b) = (middle(&plain), middle(&full)); + eprintln!("pentagon g1: the middle {a:.4} from the axis plain, {b:.4} with crown 0.1; normal within {angle:.5}° of the wall's"); + assert!(b > a + 0.5 && angle < 0.1, "{a} {b} {angle}"); + // The edges of a closed solid lie between two faces: g1 has no one face to meet. + let b = box_solid(Vec3::ZERO, Vec3::new(10.0, 10.0, 10.0)).unwrap(); + let top = [Vec3::new(5.0, 0.0, 10.0), Vec3::new(10.0, 5.0, 10.0), Vec3::new(5.0, 10.0, 10.0), Vec3::new(0.0, 5.0, 10.0)]; + assert!(patch_edges_report(&b, &top, &PatchSpec { continuity: Continuity::G1, ..PatchSpec::default() }).is_err()); +} diff --git a/crates/kernel/src/lib.rs b/crates/kernel/src/lib.rs index 77ae5f6e..21248372 100644 --- a/crates/kernel/src/lib.rs +++ b/crates/kernel/src/lib.rs @@ -122,14 +122,14 @@ pub use blend_corner::fillet_radii; pub use blend_variable::{RadiusSpline, Taper, fillet_chord, fillet_tapered, fillet_variable, fillet_variable_points, fillet_variable_spline, taper}; pub use body::{Body, CapSignature, DISPLAY_ANGLE, EdgeInfo, FaceInfo, FacePaint, FaceSource, Paint}; pub use boolean_diagnostics::{BooleanDiagnostics, BooleanFallback, BooleanStrategy, capture_boolean_diagnostics}; -pub use bridge::{Continuity, Frame, blend_between, bridge, chain, chain_closes, frames, sample_grid}; +pub use bridge::{Continuity, Frame, blend_between, bridge, chain, chain_closes, frames, frames_across, sample_grid}; pub use build::{ PathSeg, box_solid, canonical_revolve_axis, cylinder, extrude, extrude_tapered, extrude_tapered_between, loft, revolve, sphere, sweep, sweep_with, torus, }; pub use delete_face::delete_faces; pub use emboss::{cylinder_emboss_anchor, emboss_cylinder, emboss_cylinder_at, has_cylinder_emboss_wall, infer_cylinder_emboss}; -pub use fill::patch_edges_crowned; +pub use fill::{PATCH_ANGLE_TOL_DEG, PatchReport, PatchSpec, patch_edges_crowned, patch_edges_report, patch_edges_with}; pub use helix::{loft_to_point, sweep_helix, sweep_helix_mesh}; pub use iges_in::{MAX_IGES_BYTES, iges_import, iges_to_step}; pub use iges_out::iges_export_bodies; diff --git a/crates/ui-egui/src/dialogs_surface.rs b/crates/ui-egui/src/dialogs_surface.rs index d235be7d..7621066b 100644 --- a/crates/ui-egui/src/dialogs_surface.rs +++ b/crates/ui-egui/src/dialogs_surface.rs @@ -9,7 +9,7 @@ use solidcraft_engine::geom::Vec3; use crate::SolveApp; use crate::dialogs::{Kind, OP_LABELS, OPS, combo, field, op_index, row_label}; -use crate::selection::{AXES, BODIES, CURVES, EDGES, FACES, PLANES, PROFILES, SelInput}; +use crate::selection::{AXES, BODIES, CURVES, EDGES, FACES, PLANES, PROFILES, SelInput, VERTICES}; /// A surface tool dialog's values. #[derive(Clone, Debug, PartialEq)] @@ -27,7 +27,11 @@ pub enum Surf { }, Sweep(crate::dialogs_sweep::SweepOpts), Loft, - Patch, + /// `continuity`: an index into [`PATCH_CONTINUITY`] (how a patch on edges meets the faces + /// beside its loop). + Patch { + continuity: usize, + }, Unstitch, ReverseNormal, Offset { @@ -70,7 +74,7 @@ impl Surf { Surf::Revolve { .. } => "SURFACE REVOLVE", Surf::Sweep(_) => "SURFACE SWEEP", Surf::Loft => "SURFACE LOFT", - Surf::Patch => "PATCH", + Surf::Patch { .. } => "PATCH", Surf::Unstitch => "UNSTITCH", Surf::ReverseNormal => "REVERSE NORMAL", Surf::Offset { .. } => "OFFSET SURFACE", @@ -97,6 +101,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")]; + /// Which holes Untrim takes away: every one, the picked ones, or none (the outer boundary only). const UNTRIM_HOLES: [&str; 3] = ["All", "Picked", "None"]; @@ -116,7 +123,7 @@ fn inputs(k: &Surf) -> Vec { crate::dialogs_sweep::direction_input(), ], Surf::Loft => vec![SelInput::new("Sections (one sketch each)", PROFILES | CURVES, true)], - Surf::Patch => vec![SelInput::new("Boundary", PROFILES | 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)], Surf::Offset { .. } => vec![SelInput::new("Faces", FACES, true)], @@ -143,7 +150,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.patch" => Surf::Patch, + "surface.patch" => Surf::Patch { continuity: 0 }, "surface.unstitch" => Surf::Unstitch, "surface.reverse_normal" => Surf::ReverseNormal, "surface.offset" => Surf::Offset { distance: "2 mm".into(), chain: false }, @@ -247,8 +254,13 @@ pub fn rows(ui: &mut egui::Ui, k: &mut Surf) -> bool { ui.checkbox(in_place, ""); ui.end_row(); } + Surf::Patch { continuity } => { + row_label(ui, "Continuity"); + combo(ui, "pa_cont", &PATCH_CONTINUITY.map(|(label, _)| label), continuity); + ui.end_row(); + } Surf::Sweep(opts) => enter |= crate::dialogs_sweep::rows(ui, opts, "ss"), - Surf::Patch | Surf::Trim | Surf::SplitFace | Surf::Unstitch | Surf::ReverseNormal | Surf::Loft => {} + Surf::Trim | Surf::SplitFace | Surf::Unstitch | Surf::ReverseNormal | Surf::Loft => {} } enter } @@ -390,7 +402,7 @@ pub fn commands(app: &SolveApp, k: &Surf, inputs: &[SelInput], extra: &Map>()?; ("surface.loft", json!({ "sections": sections })) } - Surf::Patch => { + Surf::Patch { continuity } => { need(0, "a boundary (sketch profiles or a loop of edges)")?; let profiles: Vec<(u64, usize)> = items(inputs, 0).iter().filter_map(|x| if let Sel::Profile { sketch, index } = x { Some((*sketch, *index)) } else { None }).collect(); @@ -401,7 +413,21 @@ pub fn commands(app: &SolveApp, k: &Surf, inputs: &[SelInput], extra: &Map = items(inputs, 0).iter().filter_map(|x| if let Sel::Edge { point, .. } = x { Some(pt(*point)) } else { None }).collect(); let body = items(inputs, 0).iter().find_map(|x| if let Sel::Edge { body, .. } = x { Some(body.clone()) } else { None }); - ("surface.patch", json!({"edges": edges, "body": body})) + let mut p = json!({"edges": edges, "body": body}); + if let Some((_, word)) = PATCH_CONTINUITY.get(*continuity).filter(|_| *continuity > 0) { + p["continuity"] = json!(word); + } + let points: Vec = items(inputs, 1) + .iter() + .filter_map(|x| match x { + Sel::Vertex { point, .. } | Sel::Face { point, .. } => Some(pt(*point)), + _ => None, + }) + .collect(); + if !points.is_empty() { + p["points"] = json!(points); + } + ("surface.patch", p) } } Surf::Offset { distance, chain } => { @@ -614,7 +640,7 @@ pub fn for_feature(app: &SolveApp, kind: &FeatureKind, frame: &Mat) -> Option<(K extra.insert("convention".into(), json!(extent.convention)); return Some((Kind::Surface(k), ins, extra)); } - FeatureKind::Patch { sketch, edges, .. } => { + FeatureKind::Patch { sketch, edges, crown, continuity, points, .. } => { let mut sels: Vec = edges.iter().filter_map(edge).collect(); if let Some(sk) = sketch { let st = s.world_state(); @@ -622,7 +648,26 @@ pub fn for_feature(app: &SolveApp, kind: &FeatureKind, frame: &Mat) -> Option<(K sels.extend((0..ss.profiles.len()).map(|index| Sel::Profile { sketch: *sk, index })); } } - (Surf::Patch, vec![sels]) + let index = PATCH_CONTINUITY.iter().position(|(_, w)| *w == continuity.all.as_str()).unwrap_or(0); + // What the dialog does not show is kept as it was: the crown, edges with a + // continuity of their own, and points that are not on a face any more. + let through: Vec = points.iter().filter_map(|p| crate::dialogs::face_sel(s, *p)).collect(); + let mut extra = Map::new(); + if *crown != 0.0 { + extra.insert("crown".into(), json!(crown)); + } + if !continuity.edges.is_empty() { + extra.insert("continuity".into(), serde_json::to_value(continuity).unwrap_or(Value::Null)); + } + if through.len() != points.len() { + extra.insert("points".into(), json!(points.iter().map(|p| pt(*p)).collect::>())); + } + let k = Surf::Patch { continuity: index }; + let mut ins = inputs(&k); + for (input, picks) in ins.iter_mut().zip([sels, through]) { + input.items = picks; + } + return Some((Kind::Surface(k), ins, extra)); } FeatureKind::SurfaceOffset { faces, distance, chain, .. } => { (Surf::Offset { distance: distance.clone(), chain: *chain }, vec![faces.iter().filter_map(|p| crate::dialogs::face_sel(s, *p)).collect()]) diff --git a/crates/ui-egui/src/preselect_tests.rs b/crates/ui-egui/src/preselect_tests.rs index 63c1b00b..adebe32a 100644 --- a/crates/ui-egui/src/preselect_tests.rs +++ b/crates/ui-egui/src/preselect_tests.rs @@ -1216,3 +1216,54 @@ fn a_primitive_without_a_pick_or_with_its_position_cleared_stays_at_the_origin() let (p, _) = primitive_on(&mut app, vec![v], "solid.sphere", |_| {}); assert_eq!(p["center"], json!([40.0, 30.0, 20.0]), "{p}"); } + +#[test] +fn patch_dialog_takes_continuity_and_points_and_brings_them_back() { + // A tube with its top face taken off: the rim's edges, G1, through a point above it. + let mut app = app_with(json!([ + {"command": "solid.cylinder", "params": {"radius": 10, "height": 30, "body_name": "Tube"}}, + {"command": "surface.unstitch", "params": {"body": "Tube", "faces": [[0, 0, 30]]}} + ])); + let rim: Vec = { + let st = app.session.model.state(); + let m = st.body("Tube").unwrap().mesh(); + m.edges + .iter() + .enumerate() + .filter(|(_, e)| e.iter().all(|p| (p.z - 30.0).abs() < 1e-6)) + .map(|(i, e)| Sel::Edge { body: "Tube".into(), index: i, point: e[e.len() / 2] }) + .collect() + }; + assert!(!rim.is_empty()); + let mut d = start(&mut app, vec![], "surface.patch"); + d.inputs[0].items = rim; + d.inputs[1].items = vec![Sel::Vertex { body: "Tube".into(), point: Vec3::new(0.0, 0.0, 36.0) }]; + if let Kind::Surface(crate::dialogs_surface::Surf::Patch { continuity }) = &mut d.kind { + *continuity = 1; + } + let (id, p) = apply_commands(&app, &d).unwrap().remove(0); + assert_eq!((id.as_str(), &p["continuity"], &p["points"]), ("surface.patch", &json!("g1"), &json!([[0.0, 0.0, 36.0]])), "{p}"); + app.session.execute(&id, &p).unwrap(); + let top = app.session.model.state().bodies.last().unwrap().mesh().bounds().max.z; + assert!(top > 35.9, "the cap reaches the point: {top}"); + // Edit Feature: G1 again, and the point (on no face) kept as it was. + let fid = app.session.doc.features.last().unwrap().id; + let back = crate::dialogs::for_feature(&app, fid, None).unwrap(); + assert!(matches!(back.kind, Kind::Surface(crate::dialogs_surface::Surf::Patch { continuity: 1 }))); + let (_, p) = apply_commands(&app, &back).unwrap().remove(0); + assert_eq!((&p["params"]["continuity"], &p["params"]["points"]), (&json!("g1"), &json!([[0.0, 0.0, 36.0]])), "{p}"); +} + +#[test] +fn patch_dialog_rows_are_translated() { + for (lang, continuity, through) in [("en", "Continuity", "Pass through"), ("de", "Stetigkeit", "Durch Punkte"), ("ja", "連続性", "通過点")] + { + let mut h = crate::scenario::Harness::new(); + h.app.preferences.language = lang.into(); + h.app.home.open = false; + h.call("ui.start", json!({"command": "surface.patch"})); + h.frames(3); + assert!(h.widget(continuity).is_some(), "{lang}: the Continuity row must show {continuity:?}"); + assert!(h.widget(through).is_some(), "{lang}: the points input must show {through:?}"); + } +}