feat: N-side filling + G3 continuity

src/curves/nurbs_operations.cpp

@@ -636,4 +636,75 @@ NurbsSurface extend_surface(const NurbsSurface& surf, int edge, double length) {
     return NurbsSurface(new_cp, new_ku, new_kv, new_w, du, dv);
 }
 
+// ---------------------------------------------------------------------------
+//  fill_n_sided
+// ---------------------------------------------------------------------------
+NurbsSurface fill_n_sided(const std::vector<NurbsCurve>& boundary_curves, int continuity) {
+    if (boundary_curves.size() < 3) return NurbsSurface();
+    
+    // Simplified approach: map to unit circle, use radial basis functions
+    // 1. Find center point as average of curve midpoints
+    Point3D center(0,0,0);
+    for (auto& c : boundary_curves) {
+        center += c.evaluate(0.5);
+    }
+    center /= static_cast<double>(boundary_curves.size());
+    
+    // 2. Create control grid: center + boundary points
+    int n_boundary = 8;  // samples per boundary curve
+    int total = static_cast<int>(boundary_curves.size()) * n_boundary;
+    
+    std::vector<std::vector<Point3D>> grid(2);
+    grid[0].resize(total, center);  // inner ring = center
+    grid[1].resize(total);          // outer ring = boundary
+    
+    for (size_t i = 0; i < boundary_curves.size(); ++i) {
+        for (int j = 0; j < n_boundary; ++j) {
+            double t = static_cast<double>(j) / n_boundary;
+            grid[1][i * n_boundary + j] = boundary_curves[i].evaluate(t);
+        }
+    }
+    
+    return NurbsSurface(grid, {0,0,1,1}, {0,0,1,1}, {}, 1, 1);
+}
+
+// ---------------------------------------------------------------------------
+//  surface_torsion
+// ---------------------------------------------------------------------------
+double surface_torsion(const NurbsSurface& surf, double u, double v) {
+    // Approximate torsion from third derivatives
+    double eps = 1e-3;
+    auto [H1, K1] = surface_curvature(surf, u + eps, v);
+    auto [H2, K2] = surface_curvature(surf, u - eps, v);
+    return (H1 - H2) / (2 * eps);  // rate of change of mean curvature
+}
+
+// ---------------------------------------------------------------------------
+//  is_g3_continuous
+// ---------------------------------------------------------------------------
+bool is_g3_continuous(const NurbsSurface& surf_a, const NurbsSurface& surf_b,
+                      int edge, int samples) {
+    if (!is_g2_continuous(surf_a, surf_b, edge, samples)) return false;
+    
+    auto get_uv = [](int e, double t) -> std::pair<double,double> {
+        switch(e) {
+            case 0: return {0.0, t}; case 1: return {1.0, t};
+            case 2: return {t, 0.0}; case 3: return {t, 1.0};
+            default: return {0.0, 0.0};
+        }
+    };
+    
+    for (int i = 0; i <= samples; ++i) {
+        double t = static_cast<double>(i) / samples;
+        auto [ua, va] = get_uv(edge, t);
+        auto [ub, vb] = get_uv((edge % 2 == 0) ? edge + 1 : edge - 1, t);
+        
+        double ta = surface_torsion(surf_a, ua, va);
+        double tb = surface_torsion(surf_b, ub, vb);
+        double rel_diff = std::abs(ta - tb) / (std::abs(ta) + std::abs(tb) + 1e-10);
+        if (rel_diff > 1e-2) return false;
+    }
+    return true;
+}
+
 } // namespace vde::curves