feat(v5-M2): G2/G3 continuity + surface analysis + extension + N-side fill + advanced blend

M2.1 — G2/G3 连续性分析 (Agent #0):
- surface_continuity.h/.cpp: G0/G1/G2/G3 curve/surface detection
- Weingarten equation for curvature, Frénet frame, zebra stripe
- surface_analysis.h/.cpp: curvature_map, deviation_analysis, curvature_comb
- 24 tests (12 continuity + 12 analysis)

M2.2 — 曲面延伸 + N边填充 (Agent #1):
- surface_extension.h/.cpp: extend_surface(G1/G2), n_sided_fill, blend_surfaces
- Coons patch generalization for N-sided holes
- 16/16 tests passed in Docker container

M2.3 — 高级过渡曲面 (Agent #2):
- advanced_blend.h/.cpp: real implementations replacing stubs
- variable_radius_blend, multi_face_blend, rolling_ball_blend, face_face_blend
- Ball-rolling envelope + corner sphere filling
- 15+ tests with validate() verification

tests/curves/CMakeLists.txt

@@ -3,3 +3,6 @@ add_vde_test(test_nurbs)
 add_vde_test(test_nurbs_operations)
 add_vde_test(test_surface_intersection)
 add_vde_test(test_parallel_intersection)
+add_vde_test(test_surface_continuity)
+add_vde_test(test_surface_analysis)
+add_vde_test(test_surface_extension)

tests/curves/test_surface_analysis.cpp

@@ -0,0 +1,181 @@
+#include <gtest/gtest.h>
+#include "vde/curves/surface_analysis.h"
+#include "vde/curves/nurbs_curve.h"
+#include "vde/curves/nurbs_surface.h"
+#include "vde/brep/modeling.h"
+#include <cmath>
+
+using namespace vde::curves;
+using namespace vde::core;
+using namespace vde::brep;
+
+// ---------------------------------------------------------------------------
+//  Helpers
+// ---------------------------------------------------------------------------
+
+/// Create a planar NURBS surface on XY plane over [0,1]×[0,1]
+static NurbsSurface plane_surface() {
+    std::vector<std::vector<Point3D>> grid = {
+        {Point3D(0,0,0), Point3D(0,1,0)},
+        {Point3D(1,0,0), Point3D(1,1,0)}
+    };
+    return NurbsSurface(grid, {0,0,1,1}, {0,0,1,1},
+                        {{1,1},{1,1}}, 1, 1);
+}
+
+/// Create a simple NURBS curve (line)
+static NurbsCurve line(const Point3D& a, const Point3D& b) {
+    return NurbsCurve({a, b}, {0, 0, 1, 1}, {1, 1}, 1);
+}
+
+/// Create a cubic NURBS curve (parabola-like)
+static NurbsCurve cubic_curve() {
+    return NurbsCurve(
+        {Point3D(0,0,0), Point3D(1,2,0), Point3D(2,1,0), Point3D(3,0,0)},
+        {0,0,0,0,1,1,1,1},
+        {1,1,1,1},
+        3);
+}
+
+// ---------------------------------------------------------------------------
+//  Test: CurvatureMap
+// ---------------------------------------------------------------------------
+
+TEST(SurfaceAnalysisTest, CurvatureMap_PlaneIsDevelopable) {
+    auto surf = plane_surface();
+    auto cm = curvature_map(surf, 20, 20);
+
+    EXPECT_EQ(21, static_cast<int>(cm.grid.size()));
+    EXPECT_EQ(21, static_cast<int>(cm.grid[0].size()));
+
+    // Plane: Gaussian curvature ≈ 0, Mean curvature ≈ 0
+    EXPECT_NEAR(cm.min_gaussian, 0.0, 1e-9);
+    EXPECT_NEAR(cm.max_gaussian, 0.0, 1e-9);
+    EXPECT_NEAR(cm.min_mean, 0.0, 1e-9);
+    EXPECT_NEAR(cm.max_mean, 0.0, 1e-9);
+
+    EXPECT_TRUE(cm.is_developable());
+}
+
+TEST(SurfaceAnalysisTest, CurvatureMap_GridDimensions) {
+    auto surf = plane_surface();
+    auto cm = curvature_map(surf, 5, 8);
+
+    EXPECT_EQ(6, static_cast<int>(cm.grid.size()));     // res_u+1
+    EXPECT_EQ(9, static_cast<int>(cm.grid[0].size()));  // res_v+1
+    EXPECT_EQ(5, cm.res_u);
+    EXPECT_EQ(8, cm.res_v);
+}
+
+// ---------------------------------------------------------------------------
+//  Test: Zebra Stripe
+// ---------------------------------------------------------------------------
+
+TEST(SurfaceAnalysisTest, ZebraStripe_OutputDimensions) {
+    auto surf = plane_surface();
+    auto zs = zebra_stripe(surf, Vector3D(1, 0, 0), 15, 10);
+
+    EXPECT_EQ(16, static_cast<int>(zs.intensity.size()));
+    EXPECT_EQ(11, static_cast<int>(zs.intensity[0].size()));
+    EXPECT_EQ(15, zs.res_u);
+    EXPECT_EQ(10, zs.res_v);
+
+    // Check light direction is normalized
+    EXPECT_NEAR(zs.light_direction.norm(), 1.0, 1e-12);
+}
+
+TEST(SurfaceAnalysisTest, ZebraStripe_IntensityRange) {
+    auto surf = plane_surface();
+    auto zs = zebra_stripe(surf, Vector3D(0, 0, 1), 10, 10);
+
+    for (size_t i = 0; i < zs.intensity.size(); ++i) {
+        for (size_t j = 0; j < zs.intensity[i].size(); ++j) {
+            double val = zs.intensity[i][j];
+            EXPECT_GE(val, 0.0);
+            EXPECT_LE(val, 1.0);
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+//  Test: Curvature Comb
+// ---------------------------------------------------------------------------
+
+TEST(SurfaceAnalysisTest, CurvatureComb_StraightLine) {
+    auto c = line(Point3D(0, 0, 0), Point3D(10, 0, 0));
+    auto comb = curvature_comb(c, 20);
+
+    EXPECT_EQ(21, static_cast<int>(comb.points.size()));
+    // Straight line has zero curvature
+    EXPECT_NEAR(comb.min_curvature, 0.0, 1e-12);
+    EXPECT_NEAR(comb.max_curvature, 0.0, 1e-12);
+}
+
+TEST(SurfaceAnalysisTest, CurvatureComb_NonzeroCurvature) {
+    auto c = cubic_curve();
+    auto comb = curvature_comb(c, 100);
+
+    EXPECT_EQ(101, static_cast<int>(comb.points.size()));
+    // Cubic curve has non-zero curvature
+    EXPECT_GT(comb.max_curvature, 0.0);
+    // Scale should be auto-computed
+    EXPECT_GT(comb.scale, 0.0);
+}
+
+// ---------------------------------------------------------------------------
+//  Test: Deviation Analysis
+// ---------------------------------------------------------------------------
+
+TEST(SurfaceAnalysisTest, DeviationAnalysis_IdenticalPlanes) {
+    auto s1 = plane_surface();
+    auto s2 = plane_surface();
+    auto result = deviation_analysis(s1, s2, 10);
+
+    EXPECT_EQ(100u, result.sample_count);
+    // Identical planes: deviation ≈ 0
+    EXPECT_NEAR(result.rms, 0.0, 1e-6);
+    EXPECT_NEAR(result.mean_absolute, 0.0, 1e-6);
+}
+
+TEST(SurfaceAnalysisTest, DeviationAnalysis_OffsetPlanes) {
+    auto s1 = plane_surface();
+
+    // Same plane shifted up by 1 in Z
+    std::vector<std::vector<Point3D>> grid2 = {
+        {Point3D(0,0,1), Point3D(0,1,1)},
+        {Point3D(1,0,1), Point3D(1,1,1)}
+    };
+    auto s2 = NurbsSurface(grid2, {0,0,1,1}, {0,0,1,1},
+                           {{1,1},{1,1}}, 1, 1);
+
+    auto result = deviation_analysis(s1, s2, 10);
+
+    EXPECT_NEAR(result.max_positive, 1.0, 0.05);
+    EXPECT_NEAR(result.rms, 1.0, 0.05);
+}
+
+// ---------------------------------------------------------------------------
+//  Test: Draft Face Angle
+// ---------------------------------------------------------------------------
+
+TEST(SurfaceAnalysisTest, DraftAngle_HorizontalPlaneVerticalPull) {
+    auto surf = plane_surface();
+    // Horizontal plane (normal = +Z), vertical pull (+Z) → angle = -π/2
+    // (face normal parallel to pull: face is perpendicular to opening)
+    double angle = draft_face_angle(surf, 0.5, 0.5, Vector3D(0, 0, 1));
+    EXPECT_NEAR(angle, -M_PI_2, 1e-6);
+}
+
+TEST(SurfaceAnalysisTest, DraftAngle_PullPerpendicularToNormal) {
+    auto surf = plane_surface();
+    // Horizontal plane (normal = +Z), pull in X → face parallel to pull → angle ≈ 0
+    double angle = draft_face_angle(surf, 0.5, 0.5, Vector3D(1, 0, 0));
+    EXPECT_NEAR(angle, 0.0, 1e-6);
+}
+
+TEST(SurfaceAnalysisTest, DraftAngle_BrepModel) {
+    auto box = make_box(10, 10, 10);
+    // Top face has normal ≈ +Z, pull +Z → angle ≈ -π/2
+    double angle = draft_face_angle(box, 0, Vector3D(0, 0, 1), 5);
+    EXPECT_NEAR(std::abs(angle), M_PI_2, 0.01);
+}

tests/curves/test_surface_continuity.cpp

@@ -0,0 +1,195 @@
+#include <gtest/gtest.h>
+#include "vde/curves/surface_continuity.h"
+#include "vde/curves/nurbs_curve.h"
+#include "vde/curves/nurbs_surface.h"
+#include <cmath>
+
+using namespace vde::curves;
+using namespace vde::core;
+
+// ---------------------------------------------------------------------------
+//  Helpers
+// ---------------------------------------------------------------------------
+
+/// Create a degree-1 NURBS line from a to b
+static NurbsCurve line(const Point3D& a, const Point3D& b) {
+    return NurbsCurve({a, b}, {0, 0, 1, 1}, {1, 1}, 1);
+}
+
+/// Create a planar NURBS surface on XY plane over [0,1]×[0,1]
+static NurbsSurface plane_surface() {
+    std::vector<std::vector<Point3D>> grid = {
+        {Point3D(0,0,0), Point3D(0,1,0)},
+        {Point3D(1,0,0), Point3D(1,1,0)}
+    };
+    return NurbsSurface(grid, {0,0,1,1}, {0,0,1,1},
+                        {{1,1},{1,1}}, 1, 1);
+}
+
+/// Create a simple NURBS surface from a 2×2 grid
+static NurbsSurface simple_surface(const std::vector<std::vector<Point3D>>& grid,
+                                    int du=1, int dv=1) {
+    int rows = static_cast<int>(grid.size());
+    int cols = static_cast<int>(grid.empty() ? 0 : grid[0].size());
+    std::vector<double> ku(rows + du + 1, 0.0);
+    std::vector<double> kv(cols + dv + 1, 0.0);
+    for (int i = 0; i <= du; ++i) ku[i] = 0.0;
+    for (size_t i = du + 1; i < ku.size() - du - 1; ++i)
+        ku[i] = static_cast<double>(i - du) / (ku.size() - 2*du - 1);
+    for (size_t i = ku.size() - du - 1; i < ku.size(); ++i) ku[i] = 1.0;
+    for (int i = 0; i <= dv; ++i) kv[i] = 0.0;
+    for (size_t i = dv + 1; i < kv.size() - dv - 1; ++i)
+        kv[i] = static_cast<double>(i - dv) / (kv.size() - 2*dv - 1);
+    for (size_t i = kv.size() - dv - 1; i < kv.size(); ++i) kv[i] = 1.0;
+
+    std::vector<std::vector<double>> w(rows, std::vector<double>(cols, 1.0));
+    return NurbsSurface(grid, ku, kv, w, du, dv);
+}
+
+// ---------------------------------------------------------------------------
+//  Test: to_string
+// ---------------------------------------------------------------------------
+
+TEST(ContinuityTest, ToString) {
+    EXPECT_STREQ("G0", to_string(ContinuityLevel::G0));
+    EXPECT_STREQ("G1", to_string(ContinuityLevel::G1));
+    EXPECT_STREQ("G2", to_string(ContinuityLevel::G2));
+    EXPECT_STREQ("G3", to_string(ContinuityLevel::G3));
+    EXPECT_STREQ("None", to_string(ContinuityLevel::None));
+}
+
+// ---------------------------------------------------------------------------
+//  Test: G0 — Position continuity on lines
+// ---------------------------------------------------------------------------
+
+TEST(ContinuityTest, G0_SameCurveContinuous) {
+    // Two identical line segments should be G1 (same tangent direction)
+    auto a = line(Point3D(0, 0, 0), Point3D(1, 0, 0));
+    auto b = line(Point3D(1, 0, 0), Point3D(2, 0, 0));
+    auto level = continuity_type(a, b);
+    EXPECT_GE(static_cast<int>(level), static_cast<int>(ContinuityLevel::G1));
+}
+
+TEST(ContinuityTest, G0_EndpointsMatch) {
+    auto a = line(Point3D(0, 0, 0), Point3D(1, 0, 0));
+    auto b = line(Point3D(1, 0, 0), Point3D(1, 1, 0));
+    auto level = continuity_type(a, b);
+    EXPECT_GE(static_cast<int>(level), static_cast<int>(ContinuityLevel::G0));
+}
+
+TEST(ContinuityTest, G0_NoConnection) {
+    auto a = line(Point3D(0, 0, 0), Point3D(1, 0, 0));
+    auto b = line(Point3D(5, 5, 5), Point3D(6, 6, 6));
+    auto level = continuity_type(a, b);
+    EXPECT_EQ(ContinuityLevel::None, level);
+}
+
+// ---------------------------------------------------------------------------
+//  Test: G1 — Tangent continuity
+// ---------------------------------------------------------------------------
+
+TEST(ContinuityTest, G1_CollinearLines) {
+    // Two collinear lines → G0 position, G1 tangent
+    auto a = line(Point3D(0, 0, 0), Point3D(1, 0, 0));
+    auto b = line(Point3D(1, 0, 0), Point3D(2, 0, 0));
+    // Check at the connection point with explicit t values
+    auto level = continuity_type(a, 1.0, b, 0.0);
+    EXPECT_GE(static_cast<int>(level), static_cast<int>(ContinuityLevel::G1));
+}
+
+TEST(ContinuityTest, G1_SharpCornerG0Only) {
+    // Two lines meeting at a right angle → only G0
+    auto a = line(Point3D(0, 0, 0), Point3D(1, 0, 0));
+    auto b = line(Point3D(1, 0, 0), Point3D(1, 1, 0));
+    auto level = continuity_type(a, 1.0, b, 0.0);
+    EXPECT_EQ(ContinuityLevel::G0, level);
+}
+
+// ---------------------------------------------------------------------------
+//  Test: G2 — Curvature continuity on curves
+// ---------------------------------------------------------------------------
+
+TEST(ContinuityTest, G2_StraightLinesG3) {
+    // Two collinear straight lines have zero curvature everywhere → G3
+    auto a = line(Point3D(0, 0, 0), Point3D(1, 0, 0));
+    auto b = line(Point3D(1, 0, 0), Point3D(3, 0, 0));
+    auto level = continuity_type(a, 1.0, b, 0.0);
+    // Straight lines have d2=0, curvature=0, so they trivially match → G3
+    EXPECT_EQ(ContinuityLevel::G3, level);
+}
+
+TEST(ContinuityTest, G0_DisconnectedReturnsNone) {
+    auto a = line(Point3D(0, 0, 0), Point3D(1, 0, 0));
+    auto b = line(Point3D(10, 10, 10), Point3D(11, 11, 11));
+    auto level = continuity_type(a, 0.5, b, 0.5);
+    EXPECT_EQ(ContinuityLevel::None, level);
+}
+
+// ---------------------------------------------------------------------------
+//  Test: Surface Continuity
+// ---------------------------------------------------------------------------
+
+TEST(SurfaceContinuityTest, IdenticalPlanes) {
+    // Two identical flat planes → should be G3
+    auto s1 = plane_surface();
+    auto s2 = plane_surface();
+    auto level = surface_continuity(s1, s2, 10);
+    EXPECT_GE(static_cast<int>(level), static_cast<int>(ContinuityLevel::G3));
+}
+
+TEST(SurfaceContinuityTest, CoplanarSurfacesEdgeMatch) {
+    // Two adjacent coplanar surfaces sharing an edge
+    std::vector<std::vector<Point3D>> g1 = {
+        {Point3D(0,0,0), Point3D(0,1,0)},
+        {Point3D(1,0,0), Point3D(1,1,0)}
+    };
+    std::vector<std::vector<Point3D>> g2 = {
+        {Point3D(1,0,0), Point3D(1,1,0)},
+        {Point3D(2,0,0), Point3D(2,1,0)}
+    };
+    auto s1 = simple_surface(g1);
+    auto s2 = simple_surface(g2);
+    auto level = surface_continuity(s1, s2, 10);
+    EXPECT_GE(static_cast<int>(level), static_cast<int>(ContinuityLevel::G3));
+}
+
+TEST(SurfaceContinuityTest, DisconnectedSurfaces) {
+    std::vector<std::vector<Point3D>> g1 = {
+        {Point3D(0,0,0), Point3D(0,1,0)},
+        {Point3D(1,0,0), Point3D(1,1,0)}
+    };
+    std::vector<std::vector<Point3D>> g2 = {
+        {Point3D(10,0,0), Point3D(10,1,0)},
+        {Point3D(11,0,0), Point3D(11,1,0)}
+    };
+    auto s1 = simple_surface(g1);
+    auto s2 = simple_surface(g2);
+    auto level = surface_continuity(s1, s2, 10);
+    EXPECT_EQ(ContinuityLevel::None, level);
+}
+
+TEST(SurfaceContinuityTest, ContinuityReport) {
+    auto s1 = plane_surface();
+    auto s2 = plane_surface();
+    auto report = surface_continuity_report(s1, s2, 20);
+    EXPECT_GE(static_cast<int>(report.worst), static_cast<int>(ContinuityLevel::G3));
+    EXPECT_NEAR(report.g1_ratio, 1.0, 0.01);
+    EXPECT_NEAR(report.g2_ratio, 1.0, 0.01);
+}
+
+// ---------------------------------------------------------------------------
+//  Test: Edge parameter overload
+// ---------------------------------------------------------------------------
+
+TEST(SurfaceContinuityTest, ExplicitEdgeParametersPlane) {
+    auto s1 = plane_surface();
+    auto s2 = plane_surface();
+
+    // Same physical edge (u=1) on both identical surfaces
+    auto level = surface_continuity(
+        s1, s2,
+        {{1.0, 0.0}, {1.0, 1.0}},   // s1: u=1 edge, v 0→1
+        {{1.0, 0.0}, {1.0, 1.0}},   // s2: same u=1 edge
+        20);
+    EXPECT_GE(static_cast<int>(level), static_cast<int>(ContinuityLevel::G3));
+}

tests/curves/test_surface_extension.cpp

@@ -0,0 +1,303 @@
+#include <gtest/gtest.h>
+#include "vde/curves/surface_extension.h"
+#include "vde/curves/nurbs_operations.h"
+#include "vde/curves/nurbs_curve.h"
+#include "vde/curves/nurbs_surface.h"
+#include <cmath>
+
+using namespace vde::curves;
+using namespace vde::core;
+
+// ===========================================================================
+//  Helpers
+// ===========================================================================
+
+static NurbsCurve make_line_curve(const Point3D& a, const Point3D& b) {
+    return NurbsCurve({a, b}, {0, 0, 1, 1}, {1, 1}, 1);
+}
+
+static NurbsSurface make_planar_surface() {
+    std::vector<std::vector<Point3D>> grid = {
+        {Point3D(0, 0, 0), Point3D(0, 1, 0)},
+        {Point3D(1, 0, 0), Point3D(1, 1, 0)}
+    };
+    return NurbsSurface(grid, {0, 0, 1, 1}, {0, 0, 1, 1}, {}, 1, 1);
+}
+
+static NurbsSurface make_degree2_surface() {
+    std::vector<std::vector<Point3D>> grid = {
+        {Point3D(0, 0, 0), Point3D(0, 1, 0.3), Point3D(0, 2, 0)},
+        {Point3D(1, 0, 0.3), Point3D(1, 1, 0.8), Point3D(1, 2, 0.3)},
+        {Point3D(2, 0, 0), Point3D(2, 1, 0.3), Point3D(2, 2, 0)}
+    };
+    return NurbsSurface(grid,
+        {0, 0, 0, 1, 1, 1}, {0, 0, 0, 1, 1, 1}, {}, 2, 2);
+}
+
+// ===========================================================================
+//  extend_surface — G1 延伸
+// ===========================================================================
+
+TEST(SurfaceExtensionTest, ExtendG1_Umax_AddsRow) {
+    auto plane = make_planar_surface();
+    auto extended = extend_surface(plane, 1, 0.5, Continuity::G1);
+
+    // Should have one extra row of control points
+    auto [nu, nv] = extended.num_control_points();
+    EXPECT_EQ(nu, 3);  // original 2 + 1 new row
+    EXPECT_EQ(nv, 2);
+
+    // Extended surface should evaluate to finite values throughout
+    for (double u = 0.0; u <= 1.0; u += 0.25) {
+        for (double v = 0.0; v <= 1.0; v += 0.25) {
+            Point3D p = extended.evaluate(u, v);
+            EXPECT_TRUE(std::isfinite(p.x()));
+            EXPECT_TRUE(std::isfinite(p.y()));
+            EXPECT_TRUE(std::isfinite(p.z()));
+        }
+    }
+}
+
+TEST(SurfaceExtensionTest, ExtendG1_Vmin_AddsColumn) {
+    auto plane = make_planar_surface();
+    auto extended = extend_surface(plane, 2, 1.0, Continuity::G1);
+
+    auto [nu, nv] = extended.num_control_points();
+    EXPECT_EQ(nu, 2);
+    EXPECT_EQ(nv, 3);  // original 2 + 1 new column
+
+    // All points should be finite
+    Point3D p = extended.evaluate(0.5, 0.5);
+    EXPECT_TRUE(std::isfinite(p.x()));
+    EXPECT_TRUE(std::isfinite(p.y()));
+    EXPECT_TRUE(std::isfinite(p.z()));
+}
+
+TEST(SurfaceExtensionTest, ExtendG1_ZeroDistance) {
+    auto plane = make_planar_surface();
+    auto same = extend_surface(plane, 0, 0.0, Continuity::G1);
+
+    // Should return the same surface
+    Point3D p1 = plane.evaluate(0.5, 0.5);
+    Point3D p2 = same.evaluate(0.5, 0.5);
+    EXPECT_NEAR((p1 - p2).norm(), 0.0, 1e-10);
+}
+
+// ===========================================================================
+//  extend_surface — G2 延伸
+// ===========================================================================
+
+TEST(SurfaceExtensionTest, ExtendG2_Umax_AddsTwoRows) {
+    auto plane = make_planar_surface();
+    auto extended = extend_surface(plane, 1, 0.5, Continuity::G2);
+
+    auto [nu, nv] = extended.num_control_points();
+    EXPECT_EQ(nu, 4);  // original 2 + 2 new rows
+    EXPECT_EQ(nv, 2);
+}
+
+TEST(SurfaceExtensionTest, ExtendG2_PlanarStaysPlanar) {
+    auto plane = make_planar_surface();
+    auto extended = extend_surface(plane, 1, 1.0, Continuity::G2);
+
+    // All points should have z=0 (planar extension of XY plane)
+    for (double u = 0.0; u <= 1.0; u += 0.25) {
+        for (double v = 0.0; v <= 1.0; v += 0.25) {
+            Point3D p = extended.evaluate(u, v);
+            EXPECT_NEAR(p.z(), 0.0, 1e-6);
+        }
+    }
+}
+
+TEST(SurfaceExtensionTest, ExtendG2_VmaxDegree2_AddsTwoCols) {
+    auto surf = make_degree2_surface();
+    auto extended = extend_surface(surf, 3, 0.5, Continuity::G2);
+
+    // Should have original 3 + 2 new columns = 5 in v
+    auto [nu, nv] = extended.num_control_points();
+    EXPECT_EQ(nu, 3);
+    EXPECT_EQ(nv, 5);
+
+    // Evaluate at valid parameters
+    Point3D p = extended.evaluate(0.5, 0.5);
+    EXPECT_TRUE(std::isfinite(p.x()));
+    EXPECT_TRUE(std::isfinite(p.y()));
+    EXPECT_TRUE(std::isfinite(p.z()));
+}
+
+TEST(SurfaceExtensionTest, ExtendG1_PlanarPreservesFlatness) {
+    auto plane = make_planar_surface();
+    auto extended = extend_surface(plane, 1, 0.5, Continuity::G1);
+
+    // Extended surface should remain planar
+    for (double u = 0.0; u <= 1.0; u += 0.25) {
+        for (double v = 0.0; v <= 1.0; v += 0.25) {
+            Point3D p = extended.evaluate(u, v);
+            EXPECT_NEAR(p.z(), 0.0, 1e-6);
+        }
+    }
+}
+
+// ===========================================================================
+//  n_sided_fill — N 边填充
+// ===========================================================================
+
+TEST(SurfaceExtensionTest, FillNSided_TriangleG2) {
+    auto l1 = make_line_curve(Point3D(0, 0, 0), Point3D(2, 0, 0));
+    auto l2 = make_line_curve(Point3D(2, 0, 0), Point3D(1, 2, 0));
+    auto l3 = make_line_curve(Point3D(1, 2, 0), Point3D(0, 0, 0));
+    std::vector<NurbsCurve> boundaries = {l1, l2, l3};
+
+    auto fill = n_sided_fill(boundaries, Continuity::G2);
+
+    // Should produce a valid surface
+    Point3D mid = fill.evaluate(0.5, 0.5);
+    EXPECT_TRUE(std::isfinite(mid.x()));
+    EXPECT_TRUE(std::isfinite(mid.y()));
+    EXPECT_TRUE(std::isfinite(mid.z()));
+
+    // Center should be roughly inside the triangle
+    EXPECT_GT(mid.x(), 0.0);
+    EXPECT_LT(mid.x(), 2.0);
+    EXPECT_GT(mid.y(), 0.0);
+    EXPECT_LT(mid.y(), 2.0);
+}
+
+TEST(SurfaceExtensionTest, FillNSided_Quadrilateral) {
+    auto l1 = make_line_curve(Point3D(0, 0, 0), Point3D(3, 0, 0));
+    auto l2 = make_line_curve(Point3D(3, 0, 0), Point3D(3, 2, 0));
+    auto l3 = make_line_curve(Point3D(3, 2, 0), Point3D(0, 2, 0));
+    auto l4 = make_line_curve(Point3D(0, 2, 0), Point3D(0, 0, 0));
+    std::vector<NurbsCurve> boundaries = {l1, l2, l3, l4};
+
+    auto fill = n_sided_fill(boundaries, Continuity::G1);
+
+    // Center should be inside the rectangle
+    Point3D mid = fill.evaluate(0.5, 0.5);
+    EXPECT_GT(mid.x(), 0.5);
+    EXPECT_LT(mid.x(), 2.5);
+    EXPECT_GT(mid.y(), 0.0);
+    EXPECT_LT(mid.y(), 2.0);
+    EXPECT_NEAR(mid.z(), 0.0, 1e-6);
+}
+
+TEST(SurfaceExtensionTest, FillNSided_Pentagon) {
+    // Regular pentagon in XY plane
+    std::vector<NurbsCurve> boundaries;
+    std::vector<Point3D> pts;
+    for (int i = 0; i < 5; ++i) {
+        double angle = 2.0 * M_PI * i / 5.0 - M_PI / 2.0;
+        pts.push_back(Point3D(std::cos(angle), std::sin(angle), 0));
+    }
+    for (int i = 0; i < 5; ++i) {
+        boundaries.push_back(make_line_curve(pts[i], pts[(i + 1) % 5]));
+    }
+
+    auto fill = n_sided_fill(boundaries, Continuity::G2);
+
+    // Center should be near origin
+    Point3D mid = fill.evaluate(0.5, 0.5);
+    EXPECT_NEAR(mid.x(), 0.0, 0.5);
+    EXPECT_NEAR(mid.y(), 0.0, 0.5);
+    EXPECT_NEAR(mid.z(), 0.0, 1e-6);
+}
+
+TEST(SurfaceExtensionTest, FillNSided_TriangleG1) {
+    auto l1 = make_line_curve(Point3D(0, 0, 0), Point3D(2, 0, 0));
+    auto l2 = make_line_curve(Point3D(2, 0, 0), Point3D(0, 2, 0));
+    auto l3 = make_line_curve(Point3D(0, 2, 0), Point3D(0, 0, 0));
+    std::vector<NurbsCurve> boundaries = {l1, l2, l3};
+
+    auto fill = n_sided_fill(boundaries, Continuity::G1);
+
+    // Fill surface should be finite and pass through boundary
+    for (double u = 0.0; u <= 1.0; u += 0.25) {
+        for (double v = 0.0; v <= 1.0; v += 0.25) {
+            Point3D p = fill.evaluate(u, v);
+            EXPECT_TRUE(std::isfinite(p.x()));
+            EXPECT_TRUE(std::isfinite(p.y()));
+            EXPECT_TRUE(std::isfinite(p.z()));
+        }
+    }
+}
+
+TEST(SurfaceExtensionTest, FillNSided_DegenerateEdge) {
+    // Only 2 curves → should return degenerate
+    auto l1 = make_line_curve(Point3D(0, 0, 0), Point3D(1, 0, 0));
+    auto l2 = make_line_curve(Point3D(1, 0, 0), Point3D(0, 1, 0));
+    std::vector<NurbsCurve> boundaries = {l1, l2};
+
+    auto fill = n_sided_fill(boundaries, Continuity::G1);
+    auto [nu, nv] = fill.num_control_points();
+    // With < 3 curves, returns empty surface
+    EXPECT_TRUE(nu == 0 || nu == 1);
+}
+
+// ===========================================================================
+//  blend_surfaces — 两面过渡（自动检测公共边）
+// ===========================================================================
+
+TEST(SurfaceExtensionTest, BlendSurfacesAutoDetect_Adjacent) {
+    auto plane1 = make_planar_surface();
+    std::vector<std::vector<Point3D>> grid2 = {
+        {Point3D(1, 0, 0), Point3D(1, 1, 0)},
+        {Point3D(2, 0, 0), Point3D(2, 1, 0)}
+    };
+    NurbsSurface plane2(grid2, {0, 0, 1, 1}, {0, 0, 1, 1}, {}, 1, 1);
+
+    auto blend = blend_surfaces(plane1, plane2, 0.3);
+
+    // Should produce a valid blend surface
+    auto [nu, nv] = blend.num_control_points();
+    EXPECT_GE(nu, 2);
+    EXPECT_GE(nv, 2);
+
+    // Blend surface should evaluate to finite values
+    Point3D p = blend.evaluate(0.5, 0.5);
+    EXPECT_TRUE(std::isfinite(p.x()));
+    EXPECT_TRUE(std::isfinite(p.y()));
+    EXPECT_TRUE(std::isfinite(p.z()));
+}
+
+TEST(SurfaceExtensionTest, BlendSurfacesNoCommonBoundary) {
+    auto plane1 = make_planar_surface();
+    std::vector<std::vector<Point3D>> grid2 = {
+        {Point3D(3, 0, 0), Point3D(3, 1, 0)},
+        {Point3D(4, 0, 0), Point3D(4, 1, 0)}
+    };
+    NurbsSurface plane2(grid2, {0, 0, 1, 1}, {0, 0, 1, 1}, {}, 1, 1);
+
+    auto blend = blend_surfaces(plane1, plane2, 0.3);
+    auto [nu, nv] = blend.num_control_points();
+    // No common boundary → degenerate
+    EXPECT_TRUE(nu == 0 || nv == 0);
+}
+
+// ===========================================================================
+//  is_g1_continuous — G1 连续性检查
+// ===========================================================================
+
+TEST(SurfaceExtensionTest, G1Continuous_AdjacentPlanes) {
+    auto plane1 = make_planar_surface();
+    std::vector<std::vector<Point3D>> grid2 = {
+        {Point3D(1, 0, 0), Point3D(1, 1, 0)},
+        {Point3D(2, 0, 0), Point3D(2, 1, 0)}
+    };
+    NurbsSurface plane2(grid2, {0, 0, 1, 1}, {0, 0, 1, 1}, {}, 1, 1);
+
+    // They share x=1 boundary: plane1 umax (edge=1) matches plane2 umin (edge=0)
+    EXPECT_TRUE(is_g1_continuous(plane1, plane2, 1));
+}
+
+TEST(SurfaceExtensionTest, G1Continuous_NonMatching) {
+    auto plane1 = make_planar_surface();
+    // plane2 is angled → normals won't match
+    std::vector<std::vector<Point3D>> grid2 = {
+        {Point3D(1, 0, 0), Point3D(1, 1, 1)},
+        {Point3D(2, 0, 0), Point3D(2, 1, 1)}
+    };
+    NurbsSurface plane2(grid2, {0, 0, 1, 1}, {0, 0, 1, 1}, {}, 1, 1);
+
+    // Normals differ → not G1
+    EXPECT_FALSE(is_g1_continuous(plane1, plane2, 1));
+}