feat(v3.8): G2 curvature + continuity implementation

- surface_curvature: mean (H) and Gaussian (K) curvature via
  first/second fundamental forms with analytic first derivatives
  and finite-difference second derivatives
- is_g2_continuous: G0/G1/G2 continuity check along shared boundary
  with proper parameter range extraction and reverse-parameter
  matching on opposite edges
- extend_surface: tangent-aligned surface extension along any edge
  (umin/umax/vmin/vmax) with knot vector propagation

include/vde/curves/nurbs_operations.h

@@ -2,6 +2,7 @@
 #include "vde/curves/nurbs_surface.h"
 #include "vde/curves/nurbs_curve.h"
 #include "vde/core/point.h"
+#include <utility>
 #include <vector>
 
 namespace vde::curves {

src/curves/nurbs_operations.cpp

@@ -364,4 +364,276 @@ NurbsSurface extrude_curve(const NurbsCurve& curve,
                          1);                       // dv — linear in v-direction
 }
 
+// ---------------------------------------------------------------------------
+//  surface_curvature
+// ---------------------------------------------------------------------------
+std::pair<double, double> surface_curvature(const NurbsSurface& surf,
+                                             double u, double v) {
+    const double h = 1e-4; // finite-difference step
+
+    // First derivatives (via central differences of evaluated points)
+    Vector3D Su  = surf.derivative_u(u, v);
+    Vector3D Sv  = surf.derivative_v(u, v);
+
+    // Second derivatives (via finite differences of first derivatives)
+    Vector3D Su_c = surf.derivative_u(u + h, v);
+    Vector3D Su_p = surf.derivative_u(u - h, v);
+    Vector3D Suu = (Su_c - Su_p) / (2.0 * h);
+
+    Vector3D Suu_c = surf.derivative_u(u, v + h);
+    Vector3D Suu_p = surf.derivative_u(u, v - h);
+    Vector3D Suv = (Suu_c - Suu_p) / (2.0 * h);
+
+    Vector3D Sv_c = surf.derivative_v(u, v + h);
+    Vector3D Sv_p = surf.derivative_v(u, v - h);
+    Vector3D Svv = (Sv_c - Sv_p) / (2.0 * h);
+
+    // Normal N = (Su × Sv).normalized()
+    Vector3D N = Su.cross(Sv);
+    double nlen = N.norm();
+    if (nlen < 1e-12) return {0.0, 0.0};
+    N = N / nlen;
+
+    // First fundamental form: E, F, G
+    double E = Su.dot(Su);
+    double F = Su.dot(Sv);
+    double G = Sv.dot(Sv);
+
+    // Second fundamental form: L, M, N2  (N2 to avoid name clash)
+    double L = N.dot(Suu);
+    double M = N.dot(Suv);
+    double N2 = N.dot(Svv);
+
+    double denom = E * G - F * F;
+    if (std::abs(denom) < 1e-15) return {0.0, 0.0};
+
+    double H = (E * N2 - 2.0 * F * M + G * L) / (2.0 * denom);
+    double K = (L * N2 - M * M) / denom;
+
+    return {H, K};
+}
+
+// ---------------------------------------------------------------------------
+//  is_g2_continuous
+// ---------------------------------------------------------------------------
+bool is_g2_continuous(const NurbsSurface& surf_a, const NurbsSurface& surf_b,
+                       int edge, int samples) {
+    const auto& ku_a = surf_a.knots_u();
+    const auto& kv_a = surf_a.knots_v();
+    const auto& ku_b = surf_b.knots_u();
+    const auto& kv_b = surf_b.knots_v();
+
+    // Parameter ranges
+    double u_min_a = ku_a[surf_a.degree_u()];
+    double u_max_a = ku_a[ku_a.size() - surf_a.degree_u() - 1];
+    double v_min_a = kv_a[surf_a.degree_v()];
+    double v_max_a = kv_a[kv_a.size() - surf_a.degree_v() - 1];
+    double u_min_b = ku_b[surf_b.degree_u()];
+    double u_max_b = ku_b[ku_b.size() - surf_b.degree_u() - 1];
+    double v_min_b = kv_b[surf_b.degree_v()];
+    double v_max_b = kv_b[kv_b.size() - surf_b.degree_v() - 1];
+
+    // Tolerance values
+    const double tol_pos   = 1e-6;
+    const double tol_norm  = 1e-3;
+    const double tol_curv  = 1e-3;
+
+    // Determine which parameter varies along the boundary and which is fixed
+    bool is_u_edge_a = (edge == 0 || edge == 1);
+    double u_fix_a = (edge == 0) ? u_min_a : ((edge == 1) ? u_max_a : 0.0);
+    double v_fix_a = (edge == 2) ? v_min_a : ((edge == 3) ? v_max_a : 0.0);
+
+    // For surf_b, assume the matching edge is opposite:
+    //   edge=0(umin)→1(umax), edge=1(umax)→0(umin), edge=2(vmin)→3(vmax), edge=3(vmax)→2(vmin)
+    int edge_b;
+    switch (edge) {
+        case 0: edge_b = 1; break; // umin ↔ umax
+        case 1: edge_b = 0; break;
+        case 2: edge_b = 3; break; // vmin ↔ vmax
+        case 3: edge_b = 2; break;
+        default: return false;
+    }
+    bool is_u_edge_b = (edge_b == 0 || edge_b == 1);
+    double u_fix_b = (edge_b == 0) ? u_min_b : ((edge_b == 1) ? u_max_b : 0.0);
+    double v_fix_b = (edge_b == 2) ? v_min_b : ((edge_b == 3) ? v_max_b : 0.0);
+
+    for (int s = 0; s < samples; ++s) {
+        double t = static_cast<double>(s) / (samples - 1);
+
+        // Evaluate on surf_a along the boundary
+        double u_a, v_a;
+        if (is_u_edge_a) { u_a = u_fix_a; v_a = v_min_a + t * (v_max_a - v_min_a); }
+        else             { u_a = u_min_a + t * (u_max_a - u_min_a); v_a = v_fix_a; }
+
+        // Evaluate on surf_b along the matching boundary (reverse parameter direction)
+        double u_b, v_b;
+        if (is_u_edge_b) { u_b = u_fix_b; v_b = v_min_b + (1.0 - t) * (v_max_b - v_min_b); }
+        else             { u_b = u_min_b + (1.0 - t) * (u_max_b - u_min_b); v_b = v_fix_b; }
+
+        // G0: position continuity
+        Point3D pa = surf_a.evaluate(u_a, v_a);
+        Point3D pb = surf_b.evaluate(u_b, v_b);
+        if ((pa - pb).norm() > tol_pos) return false;
+
+        // G1: tangent plane continuity (cross-product of normals should be ~zero)
+        Vector3D Na = surf_a.normal(u_a, v_a);
+        Vector3D Nb = surf_b.normal(u_b, v_b);
+        double norm_dot = std::abs(Na.dot(Nb));
+        if (std::abs(norm_dot - 1.0) > tol_norm) return false;
+
+        // G2: curvature continuity
+        auto [Ha, Ka] = surface_curvature(surf_a, u_a, v_a);
+        auto [Hb, Kb] = surface_curvature(surf_b, u_b, v_b);
+
+        double rel_H = (std::abs(Ha) + std::abs(Hb)) > 1e-12
+            ? std::abs(Ha - Hb) / (std::abs(Ha) + std::abs(Hb)) : std::abs(Ha - Hb);
+        double rel_K = (std::abs(Ka) + std::abs(Kb)) > 1e-12
+            ? std::abs(Ka - Kb) / (std::abs(Ka) + std::abs(Kb)) : std::abs(Ka - Kb);
+        if (rel_H > tol_curv || rel_K > tol_curv) return false;
+    }
+
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+//  extend_surface
+// ---------------------------------------------------------------------------
+NurbsSurface extend_surface(const NurbsSurface& surf, int edge, double length) {
+    if (std::abs(length) < 1e-15) return surf;
+    if (length < 0.0)
+        throw std::invalid_argument("extend_surface: length must be non-negative");
+
+    const auto& cp  = surf.control_points();
+    const auto& w   = surf.weights();
+    auto ku = surf.knots_u();
+    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());
+
+    if (nu == 0 || nv == 0)
+        throw std::invalid_argument("extend_surface: empty surface");
+
+    auto gu = greville_abscissae(ku, du);
+    auto gv = greville_abscissae(kv, dv);
+
+    std::vector<std::vector<Point3D>> new_cp;
+    std::vector<std::vector<double>> new_w;
+    std::vector<double> new_ku, new_kv;
+
+    switch (edge) {
+        case 0: { // u_min: prepend row
+            double u_val = gu[0];
+            // New row of CPs at boundary + tangent direction * length
+            std::vector<Point3D> new_row(nv);
+            std::vector<double> new_row_w(nv, 1.0);
+            for (int j = 0; j < nv; ++j) {
+                Vector3D tan = surf.derivative_u(u_val, gv[j]);
+                double tlen = tan.norm();
+                // Extend in the opposite u-direction (toward smaller u)
+                new_row[j] = Point3D(
+                    cp[0][j].x() - (tlen > 1e-12 ? length * tan.x() / tlen : 0.0),
+                    cp[0][j].y() - (tlen > 1e-12 ? length * tan.y() / tlen : 0.0),
+                    cp[0][j].z() - (tlen > 1e-12 ? length * tan.z() / tlen : 0.0)
+                );
+                if (!w.empty()) new_row_w[j] = w[0][j];
+            }
+            new_cp.push_back(new_row);
+            for (int i = 0; i < nu; ++i) new_cp.push_back(cp[i]);
+            new_w = w;
+            if (!new_w.empty()) new_w.insert(new_w.begin(), new_row_w);
+            // Extend u-knots: prepend a new interval
+            double dk = ku[du + nu - 1] - ku[du + nu - 2];
+            new_ku = {ku[0] - dk};
+            new_ku.insert(new_ku.end(), ku.begin(), ku.end());
+            new_kv = kv;
+            break;
+        }
+        case 1: { // u_max: append row
+            double u_val = gu[nu - 1];
+            std::vector<Point3D> new_row(nv);
+            std::vector<double> new_row_w(nv, 1.0);
+            for (int j = 0; j < nv; ++j) {
+                Vector3D tan = surf.derivative_u(u_val, gv[j]);
+                double tlen = tan.norm();
+                new_row[j] = Point3D(
+                    cp[nu - 1][j].x() + (tlen > 1e-12 ? length * tan.x() / tlen : 0.0),
+                    cp[nu - 1][j].y() + (tlen > 1e-12 ? length * tan.y() / tlen : 0.0),
+                    cp[nu - 1][j].z() + (tlen > 1e-12 ? length * tan.z() / tlen : 0.0)
+                );
+                if (!w.empty()) new_row_w[j] = w[nu - 1][j];
+            }
+            new_cp = cp;
+            new_cp.push_back(new_row);
+            new_w = w;
+            if (!new_w.empty()) new_w.push_back(new_row_w);
+            // Extend u-knots: append a new interval
+            double dk = ku[du + nu - 1] - ku[du + nu - 2];
+            new_ku = ku;
+            new_ku.push_back(ku.back() + dk);
+            new_kv = kv;
+            break;
+        }
+        case 2: { // v_min: prepend column
+            double v_val = gv[0];
+            for (int i = 0; i < nu; ++i) {
+                Vector3D tan = surf.derivative_v(gu[i], v_val);
+                double tlen = tan.norm();
+                Point3D new_pt(
+                    cp[i][0].x() - (tlen > 1e-12 ? length * tan.x() / tlen : 0.0),
+                    cp[i][0].y() - (tlen > 1e-12 ? length * tan.y() / tlen : 0.0),
+                    cp[i][0].z() - (tlen > 1e-12 ? length * tan.z() / tlen : 0.0)
+                );
+                std::vector<Point3D> row = {new_pt};
+                row.insert(row.end(), cp[i].begin(), cp[i].end());
+                new_cp.push_back(row);
+                if (!w.empty()) {
+                    std::vector<double> row_w = {w[i][0]};
+                    row_w.insert(row_w.end(), w[i].begin(), w[i].end());
+                    new_w.push_back(row_w);
+                }
+            }
+            if (w.empty()) new_w = w;
+            // Extend v-knots: prepend a new interval
+            double dk = kv[dv + nv - 1] - kv[dv + nv - 2];
+            new_ku = ku;
+            new_kv = {kv[0] - dk};
+            new_kv.insert(new_kv.end(), kv.begin(), kv.end());
+            break;
+        }
+        case 3: { // v_max: append column
+            double v_val = gv[nv - 1];
+            for (int i = 0; i < nu; ++i) {
+                Vector3D tan = surf.derivative_v(gu[i], v_val);
+                double tlen = tan.norm();
+                Point3D new_pt(
+                    cp[i][nv - 1].x() + (tlen > 1e-12 ? length * tan.x() / tlen : 0.0),
+                    cp[i][nv - 1].y() + (tlen > 1e-12 ? length * tan.y() / tlen : 0.0),
+                    cp[i][nv - 1].z() + (tlen > 1e-12 ? length * tan.z() / tlen : 0.0)
+                );
+                std::vector<Point3D> row = cp[i];
+                row.push_back(new_pt);
+                new_cp.push_back(row);
+                if (!w.empty()) {
+                    std::vector<double> row_w = w[i];
+                    row_w.push_back(w[i][nv - 1]);
+                    new_w.push_back(row_w);
+                }
+            }
+            if (w.empty()) new_w = w;
+            // Extend v-knots: append a new interval
+            double dk = kv[dv + nv - 1] - kv[dv + nv - 2];
+            new_ku = ku;
+            new_kv = kv;
+            new_kv.push_back(kv.back() + dk);
+            break;
+        }
+        default:
+            throw std::invalid_argument("extend_surface: edge must be 0-3");
+    }
+
+    return NurbsSurface(new_cp, new_ku, new_kv, new_w, du, dv);
+}
+
 } // namespace vde::curves