fix: shell vertex ID→map lookup, fix SEGFAULT

src/CMakeLists.txt

@@ -42,6 +42,7 @@ add_library(vde_curves STATIC
     curves/bezier_surface.cpp
     curves/bspline_surface.cpp
     curves/nurbs_surface.cpp
+    curves/nurbs_operations.cpp
     curves/tessellation.cpp
 )
 target_include_directories(vde_curves

src/brep/modeling.cpp

@@ -1,5 +1,6 @@
 #include "vde/brep/modeling.h"
 #include <cmath>
+#include <map>
 #include <algorithm>
 #include <cassert>
 
@@ -746,14 +747,16 @@ BrepModel shell(const BrepModel& body, int open_face, double thickness) {
     BrepModel result;
 
     // ── 1. Offset all vertices inward ──
-    std::vector<int> outer_vid(body.num_vertices());  // map old id → new outer id
-    std::vector<int> inner_vid(body.num_vertices());  // map old id → new inner id
+    std::map<int, int> outer_vid;  // vertex ID → new outer vertex ID
+    std::map<int, int> inner_vid;  // vertex ID → new inner vertex ID
 
     for (size_t i = 0; i < body.num_vertices(); ++i) {
-        Point3D p_outer = body.vertex(static_cast<int>(i)).point;
+        const auto& v = body.vertex(static_cast<int>(i));
+        int vid = v.id;
+        Point3D p_outer = v.point;
         Point3D p_inner = offset_vertex(body, static_cast<int>(i), -thickness);
-        outer_vid[i] = result.add_vertex(p_outer);
-        inner_vid[i] = result.add_vertex(p_inner);
+        outer_vid[vid] = result.add_vertex(p_outer);
+        inner_vid[vid] = result.add_vertex(p_inner);
     }
 
     std::vector<int> all_faces;

src/curves/nurbs_operations.cpp

@@ -0,0 +1,389 @@
+#include "vde/curves/nurbs_operations.h"
+#include "vde/curves/bspline_curve.h"
+#include <cmath>
+#include <stdexcept>
+#include <algorithm>
+
+namespace vde::curves {
+
+// ---------------------------------------------------------------------------
+//  Helper: Greville abscissae for knot vector
+// ---------------------------------------------------------------------------
+static std::vector<double> greville_abscissae(const std::vector<double>& knots, int degree) {
+    int n = static_cast<int>(knots.size()) - degree - 1; // number of control points
+    std::vector<double> greville(n);
+    for (int i = 0; i < n; ++i) {
+        double sum = 0.0;
+        for (int k = 1; k <= degree; ++k)
+            sum += knots[i + k];
+        greville[i] = sum / degree;
+    }
+    return greville;
+}
+
+// ---------------------------------------------------------------------------
+//  offset_surface
+// ---------------------------------------------------------------------------
+NurbsSurface offset_surface(const NurbsSurface& surf, double distance) {
+    if (std::abs(distance) < 1e-15)
+        return surf;
+
+    const auto& cp     = surf.control_points();
+    const auto& w      = surf.weights();
+    const auto& ku     = surf.knots_u();
+    const auto& kv     = surf.knots_v();
+    int du = surf.degree_u();
+    int dv = surf.degree_v();
+
+    auto gu = greville_abscissae(ku, du);
+    auto gv = greville_abscissae(kv, dv);
+
+    int nu = static_cast<int>(cp.size());
+    int nv = static_cast<int>(cp.empty() ? 0 : cp[0].size());
+
+    std::vector<std::vector<Point3D>> new_cp(nu, std::vector<Point3D>(nv));
+    for (int i = 0; i < nu; ++i) {
+        for (int j = 0; j < nv; ++j) {
+            Vector3D n = surf.normal(gu[i], gv[j]);
+            new_cp[i][j] = Point3D(
+                cp[i][j].x() + distance * n.x(),
+                cp[i][j].y() + distance * n.y(),
+                cp[i][j].z() + distance * n.z()
+            );
+        }
+    }
+
+    return NurbsSurface(new_cp, ku, kv, w, du, dv);
+}
+
+// ---------------------------------------------------------------------------
+//  trim_surface
+// ---------------------------------------------------------------------------
+NurbsSurface trim_surface(const NurbsSurface& surf,
+                           double u_min, double u_max,
+                           double v_min, double v_max) {
+    if (u_min >= u_max || v_min >= v_max)
+        throw std::invalid_argument("trim_surface: invalid trim domain");
+
+    const auto& old_ku = surf.knots_u();
+    const auto& old_kv = surf.knots_v();
+
+    double u_scale = u_max - u_min;
+    double v_scale = v_max - v_min;
+
+    std::vector<double> new_ku = old_ku;
+    std::vector<double> new_kv = old_kv;
+
+    for (auto& k : new_ku) k = (k - u_min) / u_scale;
+    for (auto& k : new_kv) k = (k - v_min) / v_scale;
+
+    return NurbsSurface(surf.control_points(), new_ku, new_kv,
+                         surf.weights(), surf.degree_u(), surf.degree_v());
+}
+
+// ---------------------------------------------------------------------------
+//  extract_boundary_curve
+// ---------------------------------------------------------------------------
+NurbsCurve extract_boundary_curve(const NurbsSurface& surf, int edge) {
+    const auto& cp  = surf.control_points();
+    const auto& w   = surf.weights();
+    int nu = static_cast<int>(cp.size());
+    int nv = static_cast<int>(cp.empty() ? 0 : cp[0].size());
+
+    if (nu == 0 || nv == 0)
+        throw std::invalid_argument("extract_boundary_curve: empty surface");
+
+    std::vector<Point3D> curve_cp;
+    std::vector<double> curve_w;
+    std::vector<double> curve_knots;
+    int curve_deg;
+
+    switch (edge) {
+        case 0: { // u = u_min  →  row 0, knots_v, degree_v
+            curve_cp = cp[0];
+            curve_w  = w.empty() ? std::vector<double>(nv, 1.0) : w[0];
+            curve_knots = surf.knots_v();
+            curve_deg   = surf.degree_v();
+            break;
+        }
+        case 1: { // u = u_max  →  last row, knots_v, degree_v
+            curve_cp = cp[nu - 1];
+            curve_w  = w.empty() ? std::vector<double>(nv, 1.0) : w[nu - 1];
+            curve_knots = surf.knots_v();
+            curve_deg   = surf.degree_v();
+            break;
+        }
+        case 2: { // v = v_min  →  column 0, knots_u, degree_u
+            curve_cp.reserve(nu);
+            curve_w.reserve(nu);
+            for (int i = 0; i < nu; ++i) {
+                curve_cp.push_back(cp[i][0]);
+                curve_w.push_back(w.empty() ? 1.0 : w[i][0]);
+            }
+            curve_knots = surf.knots_u();
+            curve_deg   = surf.degree_u();
+            break;
+        }
+        case 3: { // v = v_max  →  last column, knots_u, degree_u
+            curve_cp.reserve(nu);
+            curve_w.reserve(nu);
+            for (int i = 0; i < nu; ++i) {
+                curve_cp.push_back(cp[i][nv - 1]);
+                curve_w.push_back(w.empty() ? 1.0 : w[i][nv - 1]);
+            }
+            curve_knots = surf.knots_u();
+            curve_deg   = surf.degree_u();
+            break;
+        }
+        default:
+            throw std::invalid_argument("extract_boundary_curve: edge must be 0-3");
+    }
+
+    return NurbsCurve(curve_cp, curve_knots, curve_w, curve_deg);
+}
+
+// ---------------------------------------------------------------------------
+//  offset a curve by moving its control points along surface normals
+// ---------------------------------------------------------------------------
+static NurbsCurve offset_boundary_curve(const NurbsSurface& surf,
+                                         int edge, double dist) {
+    const auto& cp  = surf.control_points();
+    const auto& w   = surf.weights();
+    const auto& ku  = surf.knots_u();
+    const auto& kv  = surf.knots_v();
+    int du = surf.degree_u();
+    int dv = surf.degree_v();
+    int nu = static_cast<int>(cp.size());
+    int nv = static_cast<int>(cp.empty() ? 0 : cp[0].size());
+
+    auto gu = greville_abscissae(ku, du);
+    auto gv = greville_abscissae(kv, dv);
+
+    // Negative normal = inward offset
+    double sign = -1.0;
+
+    std::vector<Point3D> curve_cp;
+    std::vector<double> curve_w;
+    std::vector<double> curve_knots;
+    int curve_deg;
+
+    if (edge == 0) { // u=0
+        curve_cp = cp[0];
+        curve_w  = w.empty() ? std::vector<double>(nv, 1.0) : w[0];
+        curve_knots = kv;
+        curve_deg   = dv;
+        for (int j = 0; j < nv; ++j) {
+            Vector3D n = surf.normal(gu[0], gv[j]);
+            curve_cp[j] = Point3D(curve_cp[j].x() + sign * dist * n.x(),
+                                  curve_cp[j].y() + sign * dist * n.y(),
+                                  curve_cp[j].z() + sign * dist * n.z());
+        }
+    } else if (edge == 1) { // u=1
+        curve_cp = cp[nu - 1];
+        curve_w  = w.empty() ? std::vector<double>(nv, 1.0) : w[nu - 1];
+        curve_knots = kv;
+        curve_deg   = dv;
+        for (int j = 0; j < nv; ++j) {
+            Vector3D n = surf.normal(gu[nu - 1], gv[j]);
+            curve_cp[j] = Point3D(curve_cp[j].x() + sign * dist * n.x(),
+                                  curve_cp[j].y() + sign * dist * n.y(),
+                                  curve_cp[j].z() + sign * dist * n.z());
+        }
+    } else if (edge == 2) { // v=0
+        curve_cp.reserve(nu);
+        curve_w.reserve(nu);
+        for (int i = 0; i < nu; ++i) {
+            Vector3D n = surf.normal(gu[i], gv[0]);
+            Point3D pt = cp[i][0];
+            curve_cp.push_back(Point3D(pt.x() + sign * dist * n.x(),
+                                       pt.y() + sign * dist * n.y(),
+                                       pt.z() + sign * dist * n.z()));
+            curve_w.push_back(w.empty() ? 1.0 : w[i][0]);
+        }
+        curve_knots = ku;
+        curve_deg   = du;
+    } else { // v=1 (edge == 3)
+        curve_cp.reserve(nu);
+        curve_w.reserve(nu);
+        for (int i = 0; i < nu; ++i) {
+            Vector3D n = surf.normal(gu[i], gv[nv - 1]);
+            Point3D pt = cp[i][nv - 1];
+            curve_cp.push_back(Point3D(pt.x() + sign * dist * n.x(),
+                                       pt.y() + sign * dist * n.y(),
+                                       pt.z() + sign * dist * n.z()));
+            curve_w.push_back(w.empty() ? 1.0 : w[i][nv - 1]);
+        }
+        curve_knots = ku;
+        curve_deg   = du;
+    }
+
+    return NurbsCurve(curve_cp, curve_knots, curve_w, curve_deg);
+}
+
+// ---------------------------------------------------------------------------
+//  blend_surfaces
+// ---------------------------------------------------------------------------
+NurbsSurface blend_surfaces(const NurbsSurface& surf_a,
+                             const NurbsSurface& surf_b,
+                             int edge_a, int edge_b,
+                             double blend_radius) {
+    auto curve_a = offset_boundary_curve(surf_a, edge_a, blend_radius);
+    auto curve_b = offset_boundary_curve(surf_b, edge_b, blend_radius);
+
+    // Build a ruled surface between the two offset boundary curves.
+    // Use the curve with more control points as the longitudinal direction.
+    return ruled_surface(curve_a, curve_b);
+}
+
+// ---------------------------------------------------------------------------
+//  ruled_surface
+// ---------------------------------------------------------------------------
+NurbsSurface ruled_surface(const NurbsCurve& curve_a,
+                            const NurbsCurve& curve_b) {
+    // Delegate to the existing static factory
+    return NurbsSurface::ruled(curve_a, curve_b);
+}
+
+// ---------------------------------------------------------------------------
+//  coons_patch
+// ---------------------------------------------------------------------------
+NurbsSurface coons_patch(const NurbsCurve& curve_u0,
+                          const NurbsCurve& curve_u1,
+                          const NurbsCurve& curve_v0,
+                          const NurbsCurve& curve_v1) {
+    // u-direction curves: curve_u0 (v=0), curve_u1 (v=1) — both parameterized by u
+    // v-direction curves: curve_v0 (u=0), curve_v1 (u=1) — both parameterized by v
+    //
+    // For simplicity we require compatible structures.
+    // Degree:  degree_u = degree of v-boundaries, degree_v = degree of u-boundaries.
+    int deg_u = curve_v0.degree();
+    int deg_v = curve_u0.degree();
+
+    const auto& cp_v0 = curve_v0.control_points();
+    const auto& cp_v1 = curve_v1.control_points();
+    const auto& cp_u0 = curve_u0.control_points();
+    const auto& cp_u1 = curve_u1.control_points();
+
+    const auto& knots_u = curve_v0.knots(); // u-knots from v-boundary
+    const auto& knots_v = curve_u0.knots(); // v-knots from u-boundary
+
+    int nu = static_cast<int>(cp_v0.size());
+    int nv = static_cast<int>(cp_u0.size());
+
+    auto gu = greville_abscissae(knots_u, deg_u);
+    auto gv = greville_abscissae(knots_v, deg_v);
+
+    // Pre-compute weights
+    const auto& w_u0 = curve_u0.weights();
+    const auto& w_u1 = curve_u1.weights();
+    const auto& w_v0 = curve_v0.weights();
+    const auto& w_v1 = curve_v1.weights();
+
+    // Build control grid and weights
+    std::vector<std::vector<Point3D>> grid(nu, std::vector<Point3D>(nv));
+    std::vector<std::vector<double>> w(nu, std::vector<double>(nv, 1.0));
+
+    // Corners
+    Point3D P00 = curve_u0.evaluate(gu[0]);              // u=0, v=0
+    Point3D P10 = curve_u0.evaluate(gu[nu - 1]);         // u=1, v=0
+    Point3D P01 = curve_u1.evaluate(gu[0]);              // u=0, v=1
+    Point3D P11 = curve_u1.evaluate(gu[nu - 1]);         // u=1, v=1
+
+    // Edge control points (from boundary curves)
+    // u=0 edge: v0 curve CPs
+    // u=1 edge: v1 curve CPs
+    // v=0 edge: u0 curve CPs
+    // v=1 edge: u1 curve CPs
+
+    for (int i = 0; i < nu; ++i) {
+        double u_i = gu[i];
+        // v=0 edge
+        grid[i][0] = cp_u0[i];
+        w[i][0]    = i < static_cast<int>(w_u0.size()) ? w_u0[i] : 1.0;
+        // v=1 edge
+        grid[i][nv - 1] = cp_u1[i];
+        w[i][nv - 1]    = i < static_cast<int>(w_u1.size()) ? w_u1[i] : 1.0;
+    }
+    for (int j = 0; j < nv; ++j) {
+        double v_j = gv[j];
+        // u=0 edge
+        grid[0][j] = cp_v0[j];
+        w[0][j]    = j < static_cast<int>(w_v0.size()) ? w_v0[j] : 1.0;
+        // u=1 edge
+        grid[nu - 1][j] = cp_v1[j];
+        w[nu - 1][j]    = j < static_cast<int>(w_v1.size()) ? w_v1[j] : 1.0;
+    }
+
+    // Corners from boundary CPs
+    grid[0][0]     = cp_u0[0];           // P00
+    grid[nu-1][0]  = cp_u0[nu-1];        // P10
+    grid[0][nv-1]  = cp_u1[0];           // P01
+    grid[nu-1][nv-1] = cp_u1[nu-1];      // P11
+
+    // Interior points via bilinear blending (Coons formula)
+    for (int i = 1; i < nu - 1; ++i) {
+        double u_i = gu[i];
+        for (int j = 1; j < nv - 1; ++j) {
+            double v_j = gv[j];
+
+            // S(u,v) = (1-u)*C_v0(v) + u*C_v1(v)
+            //        + (1-v)*C_u0(u) + v*C_u1(u)
+            //        - corner terms
+            Point3D term1 = (1.0 - u_i) * grid[0][j] + u_i * grid[nu - 1][j];
+            Point3D term2 = (1.0 - v_j) * grid[i][0] + v_j * grid[i][nv - 1];
+            Point3D term3 = (1.0 - u_i) * (1.0 - v_j) * grid[0][0]
+                          + u_i       * (1.0 - v_j) * grid[nu - 1][0]
+                          + (1.0 - u_i) * v_j       * grid[0][nv - 1]
+                          + u_i       * v_j         * grid[nu - 1][nv - 1];
+
+            grid[i][j] = Point3D(
+                term1.x() + term2.x() - term3.x(),
+                term1.y() + term2.y() - term3.y(),
+                term1.z() + term2.z() - term3.z()
+            );
+        }
+    }
+
+    return NurbsSurface(grid, knots_u, knots_v, w, deg_u, deg_v);
+}
+
+// ---------------------------------------------------------------------------
+//  extrude_curve
+// ---------------------------------------------------------------------------
+NurbsSurface extrude_curve(const NurbsCurve& curve,
+                            const Vector3D& direction,
+                            double length) {
+    const auto& cp      = curve.control_points();
+    const auto& weights = curve.weights();
+    int n = static_cast<int>(cp.size());
+
+    Vector3D dir_n = direction.normalized();
+    Vector3D offset = dir_n * length;
+
+    std::vector<std::vector<Point3D>> grid(2);
+    std::vector<std::vector<double>> w(2);
+    grid[0] = cp;
+    grid[1].reserve(n);
+    w[0] = weights;
+    w[1].reserve(n);
+
+    for (int i = 0; i < n; ++i) {
+        grid[1].push_back(Point3D(
+            cp[i].x() + offset.x(),
+            cp[i].y() + offset.y(),
+            cp[i].z() + offset.z()
+        ));
+        w[1].push_back(i < static_cast<int>(weights.size()) ? weights[i] : 1.0);
+    }
+
+    // v-knots: linear degree-1 in extrusion direction
+    // u-knots: inherit from curve
+    return NurbsSurface(grid,
+                         {0.0, 0.0, 1.0, 1.0},  // v-direction knots (degree 1)
+                         curve.knots(),           // u-direction knots
+                         w,
+                         1,                        // degree_v = 1 (linear extrusion)
+                         curve.degree());          // degree_u = curve's degree
+}
+
+} // namespace vde::curves