docs: v3.8 engineering drawing + G2 continuity plan

src/CMakeLists.txt

@@ -151,6 +151,7 @@ add_library(vde_brep STATIC
     brep/iges_export.cpp
     brep/brep_boolean.cpp
     brep/brep_validate.cpp
+    brep/brep_heal.cpp
     brep/brep_face_split.cpp
     brep/feature_tree.cpp
     brep/assembly_constraints.cpp

src/brep/brep_heal.cpp

@@ -0,0 +1,348 @@
+#include "vde/brep/brep_heal.h"
+#include <algorithm>
+#include <set>
+#include <map>
+#include <cmath>
+
+namespace vde::brep {
+
+using core::Point3D;
+using core::Vector3D;
+
+// ═════════════════════════════════════════════════════
+//  heal_merge_vertices
+// ═════════════════════════════════════════════════════
+
+int heal_merge_vertices(BrepModel& body, double tolerance) {
+    if (body.vertices_.size() < 2) return 0;
+
+    int merged = 0;
+    size_t n = body.vertices_.size();
+
+    // Build a mapping: vertex array index → which array index to replace with
+    // -1 means keep original, otherwise it's the merge-target index
+    std::vector<int> replacement(n, -1);
+
+    for (size_t i = 0; i < n; ++i) {
+        if (replacement[i] != -1) continue;
+        for (size_t j = i + 1; j < n; ++j) {
+            if (replacement[j] != -1) continue;
+            if ((body.vertices_[i].point - body.vertices_[j].point).norm() < tolerance) {
+                replacement[j] = static_cast<int>(i);
+                merged++;
+            }
+        }
+    }
+
+    if (merged == 0) return 0;
+
+    // Build ID mapping: old vertex ID → kept vertex ID
+    // Edge references use vertex IDs, not array indices
+    // Find root for each merged vertex chain
+    auto find_root = [&](int idx) -> int {
+        int r = replacement[idx];
+        while (replacement[r] != -1) r = replacement[r];
+        return r;
+    };
+
+    std::map<int, int> id_map;
+    for (size_t i = 0; i < n; ++i) {
+        if (replacement[i] == -1) {
+            id_map[body.vertices_[i].id] = body.vertices_[i].id;
+        } else {
+            int root = find_root(static_cast<int>(i));
+            id_map[body.vertices_[i].id] = body.vertices_[root].id;
+        }
+    }
+
+    // Rewrite edge vertex references (e.v_start / e.v_end are vertex IDs)
+    for (auto& e : body.edges_) {
+        auto it_s = id_map.find(e.v_start);
+        if (it_s != id_map.end()) e.v_start = it_s->second;
+        auto it_e = id_map.find(e.v_end);
+        if (it_e != id_map.end()) e.v_end = it_e->second;
+    }
+
+    // Compact vertices vector: keep only those where replacement == -1
+    size_t write = 0;
+    for (size_t read = 0; read < n; ++read) {
+        if (replacement[read] == -1) {
+            if (write != read) body.vertices_[write] = std::move(body.vertices_[read]);
+            write++;
+        }
+    }
+    body.vertices_.resize(write);
+
+    return merged;
+}
+
+// ═════════════════════════════════════════════════════
+//  heal_merge_edges
+// ═════════════════════════════════════════════════════
+
+namespace {
+
+/// Check whether two NURBS curves represent (nearly) the same geometry
+bool curves_equivalent(const curves::NurbsCurve& a, const curves::NurbsCurve& b, double tol) {
+    // Quick check: same control point count
+    if (a.control_points().size() != b.control_points().size()) return false;
+
+    // Check control points are within tolerance
+    const auto& cpa = a.control_points();
+    const auto& cpb = b.control_points();
+    for (size_t i = 0; i < cpa.size(); ++i) {
+        if ((cpa[i] - cpb[i]).norm() > tol) return false;
+    }
+    return true;
+}
+
+} // namespace
+
+int heal_merge_edges(BrepModel& body, double tolerance) {
+    if (body.edges_.size() < 2) return 0;
+
+    int merged = 0;
+    size_t n = body.edges_.size();
+    std::vector<int> replacement(n, -1);
+
+    for (size_t i = 0; i < n; ++i) {
+        if (replacement[i] != -1) continue;
+        const auto& ei = body.edges_[i];
+        for (size_t j = i + 1; j < n; ++j) {
+            if (replacement[j] != -1) continue;
+            const auto& ej = body.edges_[j];
+
+            // Two edges are coincident if they share the same endpoints
+            // (accounting for possible reversal) and have equivalent curve geometry
+            bool same_dir = (ei.v_start == ej.v_start && ei.v_end == ej.v_end);
+            bool opp_dir  = (ei.v_start == ej.v_end  && ei.v_end == ej.v_start);
+
+            if (!same_dir && !opp_dir) continue;
+
+            // Check curve equivalence
+            bool curves_match = false;
+            if (ei.curve && ej.curve) {
+                curves_match = curves_equivalent(*ei.curve, *ej.curve, tolerance);
+            } else if (!ei.curve && !ej.curve) {
+                // Both are straight-line edges — check if their endpoint geometry matches
+                // (endpoints already match by vertex ID, so curves are same line segment)
+                curves_match = true;
+            }
+            // If one has a curve and the other doesn't, they're not equivalent
+
+            if (curves_match) {
+                replacement[j] = static_cast<int>(i);
+                merged++;
+            }
+        }
+    }
+
+    if (merged == 0) return 0;
+
+    // Build array-index mapping: old array index → final array index after compaction
+    // Loop references store array indices (from add_edge return values)
+    auto find_root = [&](int idx) -> int {
+        int r = replacement[idx];
+        while (replacement[r] != -1) r = replacement[r];
+        return r;
+    };
+
+    // Compute where each original index ends up after compaction
+    std::vector<int> final_pos(n, -1);
+    int pos = 0;
+    for (size_t i = 0; i < n; ++i) {
+        if (replacement[i] == -1) {
+            final_pos[i] = pos++;
+        }
+    }
+    for (size_t i = 0; i < n; ++i) {
+        if (replacement[i] != -1) {
+            final_pos[i] = final_pos[find_root(static_cast<int>(i))];
+        }
+    }
+
+    // Rewrite loop edge references (store array indices)
+    for (auto& loop : body.loops_) {
+        for (size_t k = 0; k < loop.edges.size(); ++k) {
+            int old_idx = loop.edges[k];
+            if (old_idx >= 0 && old_idx < static_cast<int>(n)) {
+                loop.edges[k] = final_pos[old_idx];
+            }
+        }
+    }
+
+    // Compact edges vector
+    size_t write = 0;
+    for (size_t read = 0; read < n; ++read) {
+        if (replacement[read] == -1) {
+            if (write != read) {
+                body.edges_[write] = std::move(body.edges_[read]);
+            }
+            write++;
+        }
+    }
+    body.edges_.resize(write);
+
+    return merged;
+}
+
+// ═════════════════════════════════════════════════════
+//  heal_close_gaps
+// ═════════════════════════════════════════════════════
+
+int heal_close_gaps(BrepModel& body, double max_gap) {
+    if (body.edges_.size() < 2) return 0;
+
+    int closed = 0;
+
+    // Collect all edge endpoints as (vertex_index, vertex_id, point)
+    // vertex_index is the position in vertices_ array
+    struct Endpoint {
+        size_t vert_idx;
+        int vert_id;
+        Point3D point;
+    };
+
+    // Build a map from vertex_id → vertices_ index
+    std::map<int, size_t> vert_id_to_idx;
+    for (size_t i = 0; i < body.vertices_.size(); ++i) {
+        vert_id_to_idx[body.vertices_[i].id] = i;
+    }
+
+    // For each vertex, find all other vertices within max_gap
+    // and snap them to the first one (greedy)
+    for (size_t i = 0; i < body.vertices_.size(); ++i) {
+        const auto& vi = body.vertices_[i];
+        for (size_t j = i + 1; j < body.vertices_.size(); ++j) {
+            auto& vj = body.vertices_[j];
+            double dist = (vi.point - vj.point).norm();
+            if (dist > 0 && dist < max_gap) {
+                // Snap vj to vi
+                vj.point = vi.point;
+                // Update edges that reference vj.id to reference vi.id instead
+                int old_id = vj.id;
+                int keep_id = vi.id;
+                for (auto& e : body.edges_) {
+                    if (e.v_start == old_id) e.v_start = keep_id;
+                    if (e.v_end   == old_id) e.v_end   = keep_id;
+                }
+                closed++;
+            }
+        }
+    }
+
+    return closed;
+}
+
+// ═════════════════════════════════════════════════════
+//  heal_orientation
+// ═════════════════════════════════════════════════════
+
+namespace {
+
+/// Compute a face's centroid by averaging its edge endpoints
+Point3D face_centroid(const BrepModel& body, const TopoFace& face) {
+    Point3D sum(0, 0, 0);
+    int count = 0;
+    for (int li : face.loops) {
+        for (const auto& loop : body.all_loops()) {
+            if (loop.id != li) continue;
+            for (int ei : loop.edges) {
+                if (ei < 0 || ei >= static_cast<int>(body.num_edges())) continue;
+                const auto& e = body.edge(ei);
+                for (size_t vi = 0; vi < body.num_vertices(); ++vi) {
+                    const auto& v = body.vertex(static_cast<int>(vi));
+                    if (v.id == e.v_start) { sum = sum + v.point; count++; }
+                    if (v.id == e.v_end)   { sum = sum + v.point; count++; }
+                }
+            }
+        }
+    }
+    if (count == 0) return Point3D(0, 0, 0);
+    return Point3D(sum.x() / count, sum.y() / count, sum.z() / count);
+}
+
+/// Estimate the face normal using the surface normal at face center
+Vector3D face_normal(const BrepModel& body, const TopoFace& face) {
+    if (face.surface_id < 0 || face.surface_id >= static_cast<int>(body.num_surfaces())) {
+        return Vector3D(0, 0, 1);
+    }
+    const auto& surf = body.surface(face.surface_id);
+    // Evaluate surface normal at center of parametric domain
+    double u = 0.5, v = 0.5;
+    Vector3D n = surf.normal(u, v);
+    if (face.reversed) n = -n;
+    return n;
+}
+
+/// Compute approximate centroid of entire body
+Point3D body_centroid(const BrepModel& body) {
+    Point3D sum(0, 0, 0);
+    if (body.num_vertices() == 0) return sum;
+    for (size_t i = 0; i < body.num_vertices(); ++i) {
+        sum = sum + body.vertex(static_cast<int>(i)).point;
+    }
+    double n = static_cast<double>(body.num_vertices());
+    return Point3D(sum.x() / n, sum.y() / n, sum.z() / n);
+}
+
+} // namespace
+
+int heal_orientation(BrepModel& body) {
+    if (body.faces_.empty()) return 0;
+
+    int flipped = 0;
+    Point3D centroid = body_centroid(body);
+
+    for (auto& face : body.faces_) {
+        // Compute face centroid and normal direction (outward-pointing)
+        Point3D fc = face_centroid(body, face);
+        Vector3D n = face_normal(body, face);
+
+        // Vector from body centroid to face centroid
+        Vector3D to_face(fc.x() - centroid.x(), fc.y() - centroid.y(), fc.z() - centroid.z());
+
+        // For outward-pointing normals: normal should point in same direction
+        // as the vector from body center to face center
+        double dot = n.x() * to_face.x() + n.y() * to_face.y() + n.z() * to_face.z();
+
+        if (dot < 0) {
+            // Normal points inward — flip the face
+            face.reversed = !face.reversed;
+            flipped++;
+        }
+    }
+
+    return flipped;
+}
+
+// ═════════════════════════════════════════════════════
+//  heal_topology
+// ═════════════════════════════════════════════════════
+
+bool heal_topology(BrepModel& body, double tolerance) {
+    // 1. Merge coincident vertices first — this simplifies subsequent steps
+    int vert_merged = heal_merge_vertices(body, tolerance);
+    (void)vert_merged;
+
+    // 2. Merge coincident edges
+    int edge_merged = heal_merge_edges(body, tolerance);
+    (void)edge_merged;
+
+    // 3. Close small gaps between faces
+    int gaps_closed = heal_close_gaps(body, tolerance * 100); // larger gap tolerance
+    (void)gaps_closed;
+
+    // 4. Re-merge vertices after gap closing (may create new coincidences)
+    vert_merged += heal_merge_vertices(body, tolerance);
+
+    // 5. Unify face orientations
+    int faces_flipped = heal_orientation(body);
+    (void)faces_flipped;
+
+    // Validate the healed model
+    auto result = validate(body);
+    return result.valid;
+}
+
+} // namespace vde::brep

src/brep/modeling.cpp

@@ -556,6 +556,219 @@ BrepModel fillet(const BrepModel& body, int edge_id, double radius) {
     return result;
 }
 
+// ── Variable-radius fillet ──
+BrepModel fillet_variable(const BrepModel& body, int edge_id,
+                          double r_start, double r_end, int samples) {
+    int num_edges = static_cast<int>(body.num_edges());
+    if (edge_id < 0 || edge_id >= num_edges) return body;
+    if (r_start < 0.0 || r_end < 0.0) return body;
+    if (r_start == 0.0 && r_end == 0.0) return body;
+
+    auto efs = body.edge_faces(edge_id);
+    if (efs.empty()) return body;
+    int face_a = efs[0];
+
+    const auto& edge = body.edge(edge_id);
+    Point3D ep0 = body.vertex_by_id(edge.v_start).point;
+    Point3D ep1 = body.vertex_by_id(edge.v_end).point;
+    int face_b = -1;
+    for (size_t fi = 0; fi < body.num_faces(); ++fi) {
+        if (static_cast<int>(fi) == face_a) continue;
+        auto f_edges = body.face_edges(static_cast<int>(fi));
+        for (int ei : f_edges) {
+            const auto& e2 = body.edge(ei);
+            Point3D q0 = body.vertex_by_id(e2.v_start).point;
+            Point3D q1 = body.vertex_by_id(e2.v_end).point;
+            bool same_edge = ((ep0-q0).norm() < 1e-7 && (ep1-q1).norm() < 1e-7) ||
+                             ((ep0-q1).norm() < 1e-7 && (ep1-q0).norm() < 1e-7);
+            if (same_edge) { face_b = static_cast<int>(fi); break; }
+        }
+        if (face_b >= 0) break;
+    }
+    if (face_b < 0) return body;
+
+    const auto& curve = *edge.curve;
+    Vector3D n_a = compute_face_normal(body, face_a);
+    Vector3D n_b = compute_face_normal(body, face_b);
+
+    double cos_angle = n_a.dot(n_b);
+    if (std::abs(cos_angle) > 0.9999) return body;
+
+    Vector3D bisector = (n_a + n_b).normalized();
+    double cos_half = bisector.dot(n_a);
+    if (std::abs(cos_half) < 1e-8) return body;
+
+    auto [t_min, t_max] = curve.domain();
+
+    int N = samples;
+    std::vector<double> radii(N + 1);
+    std::vector<Point3D> edge_pts(N + 1);
+    std::vector<Point3D> t1_pts(N + 1);
+    std::vector<Point3D> t2_pts(N + 1);
+    std::vector<Point3D> mid_pts(N + 1);
+
+    for (int i = 0; i <= N; ++i) {
+        double t_val = static_cast<double>(i) / N;
+        double r = r_start + (r_end - r_start) * t_val;
+        double t = t_min + (t_max - t_min) * t_val;
+        radii[i] = r;
+        edge_pts[i] = curve.evaluate(t);
+
+        double offset = r / cos_half;
+        Point3D C = edge_pts[i] + offset * bisector;
+
+        t1_pts[i] = C - r * n_a;
+        t2_pts[i] = C - r * n_b;
+
+        Vector3D to_mid = (t1_pts[i] + t2_pts[i]) * 0.5 - C;
+        double mid_len = to_mid.norm();
+        if (mid_len > 1e-10) {
+            mid_pts[i] = C + to_mid.normalized() * r;
+        } else {
+            mid_pts[i] = (t1_pts[i] + t2_pts[i]) * 0.5;
+        }
+    }
+
+    BrepModel result;
+
+    // ── 1. Copy non-affected faces verbatim ──
+    for (size_t fi = 0; fi < body.num_faces(); ++fi) {
+        if (static_cast<int>(fi) == face_a || static_cast<int>(fi) == face_b)
+            continue;
+        const auto& f_src = body.face(static_cast<int>(fi));
+        std::vector<int> new_loop_ids;
+        for (int li : f_src.loops) {
+            const auto& lop = body.loop_by_id(li);
+            std::vector<int> new_es;
+            for (int ei : lop.edges) {
+                const auto& e_src = body.edge(ei);
+                Point3D p0 = body.vertex_by_id(e_src.v_start).point;
+                Point3D p1 = body.vertex_by_id(e_src.v_end).point;
+                int nv0 = result.add_vertex(p0);
+                int nv1 = result.add_vertex(p1);
+                new_es.push_back(result.add_edge(nv0, nv1));
+            }
+            new_loop_ids.push_back(result.add_loop(new_es, lop.is_outer));
+        }
+        int new_surf = result.add_surface(body.surface(f_src.surface_id));
+        result.add_face(new_surf, new_loop_ids);
+    }
+
+    // ── 2. Rebuild adjacent faces with variable tangent offset ──
+    auto rebuild_face = [&](int face_id, const Vector3D& n,
+                            const std::vector<Point3D>& tangent_pts) {
+        const auto& f_src = body.face(face_id);
+        const auto& curve_copy = body.surface(f_src.surface_id);
+        const auto& edge_src = body.edge(edge_id);
+
+        auto f_edges = body.face_edges(face_id);
+        int fillet_edge_idx = -1;
+        for (size_t ei = 0; ei < f_edges.size(); ++ei) {
+            if (f_edges[ei] == edge_id) {
+                fillet_edge_idx = static_cast<int>(ei);
+                break;
+            }
+        }
+        if (fillet_edge_idx < 0) return;
+
+        std::vector<int> new_es;
+        for (const auto& lop : body.all_loops()) {
+            bool found = false;
+            for (int li : f_src.loops) {
+                if (lop.id == li) { found = true; break; }
+            }
+            if (!found) continue;
+
+            for (size_t ei = 0; ei < lop.edges.size(); ++ei) {
+                int eidx = lop.edges[ei];
+                if (eidx != edge_id) continue;
+
+                bool edge_reversed = (edge_src.v_start != body.edge(eidx).v_start);
+
+                std::vector<int> sorted_edges(lop.edges.size());
+                int start_offset = static_cast<int>((ei + 1) % lop.edges.size());
+                for (size_t j = 0; j < lop.edges.size(); ++j) {
+                    sorted_edges[j] = lop.edges[(start_offset + static_cast<int>(j)) % lop.edges.size()];
+                }
+                sorted_edges.insert(sorted_edges.begin(), edge_id);
+
+                Point3D t0 = edge_reversed ? tangent_pts.back() : tangent_pts.front();
+                Point3D t1 = edge_reversed ? tangent_pts.front() : tangent_pts.back();
+                new_es.push_back(result.add_edge(result.add_vertex(t0),
+                                                  result.add_vertex(t1)));
+
+                for (size_t j = 1; j < sorted_edges.size(); ++j) {
+                    int old_ei = sorted_edges[j];
+                    const auto& e_src = body.edge(old_ei);
+                    Point3D p0 = body.vertex_by_id(e_src.v_start).point;
+                    Point3D p1 = body.vertex_by_id(e_src.v_end).point;
+
+                    bool v0_on = (e_src.v_start == edge_src.v_start ||
+                                  e_src.v_start == edge_src.v_end);
+                    bool v1_on = (e_src.v_end == edge_src.v_start ||
+                                  e_src.v_end == edge_src.v_end);
+
+                    Point3D q0 = p0, q1 = p1;
+                    if (v0_on) q0 = (e_src.v_start == edge_src.v_start) ?
+                                   tangent_pts.front() : tangent_pts.back();
+                    if (v1_on) q1 = (e_src.v_end == edge_src.v_start) ?
+                                   tangent_pts.front() : tangent_pts.back();
+
+                    new_es.push_back(result.add_edge(
+                        result.add_vertex(q0), result.add_vertex(q1)));
+                }
+                break;
+            }
+            break;
+        }
+
+        if (!new_es.empty()) {
+            int new_loop = result.add_loop(new_es, true);
+            int new_surf = result.add_surface(curve_copy);
+            result.add_face(new_surf, {new_loop});
+        }
+    };
+
+    rebuild_face(face_a, n_a, t1_pts);
+    rebuild_face(face_b, n_b, t2_pts);
+
+    // ── 3. Build variable-radius fillet surface faces ──
+    for (int i = 0; i < N; ++i) {
+        Point3D q00 = t1_pts[i],     q01 = t2_pts[i];
+        Point3D q10 = t1_pts[i + 1], q11 = t2_pts[i + 1];
+
+        std::vector<std::vector<Point3D>> grid = {
+            {q00, mid_pts[i],     q01},
+            {q10, mid_pts[i + 1], q11}
+        };
+        curves::NurbsSurface fs(grid,
+            {0, 0, 1, 1},
+            {0, 0, 0, 1, 1, 1},
+            {}, 1, 2);
+
+        int v0 = result.add_vertex(q00), v1 = result.add_vertex(q10);
+        int v2 = result.add_vertex(q11), v3 = result.add_vertex(q01);
+        int e1 = result.add_edge(v0, v1);
+        int e2 = result.add_edge(v1, v2);
+        int e3 = result.add_edge(v2, v3);
+        int e4 = result.add_edge(v3, v0);
+        int surf_id = result.add_surface(fs);
+        result.add_face(surf_id,
+                        {result.add_loop({e1, e2, e3, e4}, true)});
+    }
+
+    // ── 4. Build shell and body ──
+    std::vector<int> all_face_ids;
+    for (size_t fi = 0; fi < result.num_faces(); ++fi) {
+        all_face_ids.push_back(result.face(static_cast<int>(fi)).id);
+    }
+    if (!all_face_ids.empty()) {
+        int sh = result.add_shell(all_face_ids, true);
+        result.add_body({sh}, "fillet_variable");
+    }
+    return result;
+}
+
 // ── Chamfer ──
 BrepModel chamfer(const BrepModel& body, int edge_id, double dist) {
     // Validation