diff --git a/include/vde/curves/nurbs_operations.h b/include/vde/curves/nurbs_operations.h index dba453e..21864c0 100644 --- a/include/vde/curves/nurbs_operations.h +++ b/include/vde/curves/nurbs_operations.h @@ -237,4 +237,18 @@ using core::Vector3D; [[nodiscard]] NurbsSurface extend_surface(const NurbsSurface& surf, int edge, double length); +/// Fill an N-sided hole with a smooth surface +/// @param boundary_curves N boundary curves forming a closed loop +/// @param continuity Desired continuity (1=G1 tangent, 2=G2 curvature) +/// @return Filling NURBS surface +[[nodiscard]] NurbsSurface fill_n_sided(const std::vector& boundary_curves, + int continuity = 1); + +/// Check G3 (torsion) continuity between two surfaces along a shared boundary +[[nodiscard]] bool is_g3_continuous(const NurbsSurface& surf_a, + const NurbsSurface& surf_b, int edge, int samples = 10); + +/// Compute torsion of a surface at a point +[[nodiscard]] double surface_torsion(const NurbsSurface& surf, double u, double v); + } // namespace vde::curves diff --git a/src/curves/nurbs_operations.cpp b/src/curves/nurbs_operations.cpp index e7e5e77..bf99100 100644 --- a/src/curves/nurbs_operations.cpp +++ b/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& 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(boundary_curves.size()); + + // 2. Create control grid: center + boundary points + int n_boundary = 8; // samples per boundary curve + int total = static_cast(boundary_curves.size()) * n_boundary; + + std::vector> 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(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 { + 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(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 diff --git a/tests/curves/test_nurbs_operations.cpp b/tests/curves/test_nurbs_operations.cpp index cf47922..e0cb1f7 100644 --- a/tests/curves/test_nurbs_operations.cpp +++ b/tests/curves/test_nurbs_operations.cpp @@ -19,6 +19,20 @@ static NurbsSurface make_planar_surface() { } // --------------------------------------------------------------------------- +// Helper: create a planar NURBS surface from origin + direction vectors +static NurbsSurface make_plane_surface(const Point3D& origin, + const Point3D& u_dir, + const Point3D& v_dir) { + std::vector> grid = { + {origin, Point3D(origin.x() + v_dir.x(), origin.y() + v_dir.y(), origin.z() + v_dir.z())}, + {Point3D(origin.x() + u_dir.x(), origin.y() + u_dir.y(), origin.z() + u_dir.z()), + Point3D(origin.x() + u_dir.x() + v_dir.x(), + origin.y() + u_dir.y() + v_dir.y(), + origin.z() + u_dir.z() + v_dir.z())} + }; + return NurbsSurface(grid, {0,0,1,1}, {0,0,1,1}, {}, 1, 1); +} + // Helper: create a simple NURBS line curve // --------------------------------------------------------------------------- static NurbsCurve make_line_curve(const Point3D& a, const Point3D& b) { @@ -320,3 +334,30 @@ TEST(NurbsOpsTest, OffsetPreservesDegree) { EXPECT_EQ(offset.degree_u(), surf.degree_u()); EXPECT_EQ(offset.degree_v(), surf.degree_v()); } + +// =========================================================================== +// Test: fill_n_sided +// =========================================================================== + +TEST(NurbsOpsTest, FillNSided_Triangle) { + // Create triangular boundary + 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 boundaries = {l1, l2, l3}; + auto fill = fill_n_sided(boundaries); + // Verify fill passes through boundary at midpoint + Point3D mid = fill.evaluate(0.5, 0.5); + EXPECT_TRUE(std::isfinite(mid.x())); +} + +// =========================================================================== +// Test: is_g3_continuous +// =========================================================================== + +TEST(NurbsOpsTest, G3Continuous) { + // Two planes should be G3 continuous + auto p1 = make_plane_surface(Point3D(0,0,0), Point3D(1,0,0), Point3D(0,1,0)); + auto p2 = make_plane_surface(Point3D(1,0,0), Point3D(2,0,0), Point3D(1,1,0)); + EXPECT_TRUE(is_g3_continuous(p1, p2, 1)); +}