From e70689d29f4701bb0a3f90f5ddf35bda7161f899 Mon Sep 17 00:00:00 2001 From: Tarik Moussa Date: Tue, 12 May 2026 17:35:00 +0200 Subject: [PATCH] feat(phase4a+4b): Newton solver + CGAL mesh I/O MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Phase 4a — newton_solver.hpp: - newton_euclidean(): SimplicialLDLT on H (PSD); solves H·Δx = −G - newton_spherical(): SimplicialLDLT on −H (NSD→PSD); solves (−H)·Δx = G - Backtracking line search (halving α up to 20×) for global convergence - NewtonResult struct: x, iterations, grad_inf_norm, converged - 7 tests: 4 spherical (convergence, few iters, large perturbation, field consistency) + 3 Euclidean (triangle pinned, quad pinned, mixed pinned — all with natural-theta equilibrium at x=0) Phase 4b — mesh_io.hpp: - read_mesh() / write_mesh(): CGAL::IO::read/write_polygon_mesh wrappers - load_mesh() / save_mesh(): throwing convenience versions - Supports OFF, OBJ, PLY (format detected by file extension) - 6 tests: OFF round-trip (tet + quad), OBJ round-trip, missing-file throw, vertex-position preservation, save/load convenience wrappers All 75 cgal tests pass (3 skipped as before). Co-Authored-By: Claude Sonnet 4.5 --- code/include/mesh_io.hpp | 65 +++++++ code/include/newton_solver.hpp | 204 +++++++++++++++++++++ code/tests/cgal/CMakeLists.txt | 6 + code/tests/cgal/test_mesh_io.cpp | 170 ++++++++++++++++++ code/tests/cgal/test_newton_solver.cpp | 235 +++++++++++++++++++++++++ 5 files changed, 680 insertions(+) create mode 100644 code/include/mesh_io.hpp create mode 100644 code/include/newton_solver.hpp create mode 100644 code/tests/cgal/test_mesh_io.cpp create mode 100644 code/tests/cgal/test_newton_solver.cpp diff --git a/code/include/mesh_io.hpp b/code/include/mesh_io.hpp new file mode 100644 index 0000000..0ac1fc1 --- /dev/null +++ b/code/include/mesh_io.hpp @@ -0,0 +1,65 @@ +#pragma once +// mesh_io.hpp +// +// Phase 4b — CGAL::IO wrappers for ConformalMesh. +// +// Reads and writes ConformalMesh (CGAL::Surface_mesh) in standard polygon-mesh +// formats using the CGAL Polygon Mesh I/O utilities (CGAL 5.4+). +// +// Supported formats (detected by file extension): +// .off — Object File Format (read + write) +// .obj — Wavefront OBJ (read + write) +// .ply — Polygon File Format (read + write) +// +// Usage: +// ConformalMesh mesh; +// if (!read_mesh("input.off", mesh)) throw std::runtime_error("read failed"); +// // ... process mesh ... +// write_mesh("output.off", mesh); +// +// Note: property maps (lambda0, v_idx, etc.) are NOT serialised — they must +// be re-initialised with setup_*_maps() + compute_lambda0_from_mesh() after +// reading a file. + +#include "conformal_mesh.hpp" +#include +#include +#include + +namespace conformallab { + +// ── Read ────────────────────────────────────────────────────────────────────── +// +// Reads a polygon mesh from file into `mesh` (clears any existing content). +// Returns true on success, false on failure. +inline bool read_mesh(const std::string& filename, ConformalMesh& mesh) +{ + mesh.clear(); + return CGAL::IO::read_polygon_mesh(filename, mesh); +} + +// ── Write ───────────────────────────────────────────────────────────────────── +// +// Writes `mesh` to `filename`. Returns true on success. +inline bool write_mesh(const std::string& filename, const ConformalMesh& mesh) +{ + return CGAL::IO::write_polygon_mesh(filename, mesh); +} + +// ── Convenience: throwing wrappers ──────────────────────────────────────────── + +inline ConformalMesh load_mesh(const std::string& filename) +{ + ConformalMesh mesh; + if (!read_mesh(filename, mesh)) + throw std::runtime_error("conformallab: failed to read mesh from " + filename); + return mesh; +} + +inline void save_mesh(const std::string& filename, const ConformalMesh& mesh) +{ + if (!write_mesh(filename, mesh)) + throw std::runtime_error("conformallab: failed to write mesh to " + filename); +} + +} // namespace conformallab diff --git a/code/include/newton_solver.hpp b/code/include/newton_solver.hpp new file mode 100644 index 0000000..2d82857 --- /dev/null +++ b/code/include/newton_solver.hpp @@ -0,0 +1,204 @@ +#pragma once +// newton_solver.hpp +// +// Phase 4a — Newton solver for the discrete conformal functionals. +// +// Solves G(x) = 0 where G is the gradient of the discrete conformal energy: +// G_v = Θ_v − Σ_f α_v^f (angle-sum residual at each vertex DOF) +// +// Algorithm per iteration: +// 1. Compute gradient G(x) +// 2. Check convergence: max|G_i| < tol → done +// 3. Compute sparse Hessian H(x) +// 4. Factorize and solve the Newton system: +// Euclidean: H · Δx = −G (H is PSD → SimplicialLDLT directly) +// Spherical: (−H) · Δx = G (H is NSD → negate to get PSD matrix) +// 5. Backtracking line search: halve α until ||G(x+α·Δx)|| < ||G(x)|| +// 6. x ← x + α·Δx, go to 1 +// +// Requires: +// Eigen::SimplicialLDLT (part of Eigen's sparse Cholesky module) + +#include "euclidean_hessian.hpp" +#include "spherical_hessian.hpp" +#include +#include +#include +#include + +namespace conformallab { + +// ── Result ──────────────────────────────────────────────────────────────────── + +struct NewtonResult { + std::vector x; ///< DOF vector at termination + int iterations; ///< Newton steps taken + double grad_inf_norm;///< max |G_i| at termination + bool converged; ///< true iff grad_inf_norm < tol +}; + +// ── Internal helpers ────────────────────────────────────────────────────────── + +namespace detail { + +// Backtracking line search: find the largest α in {1, 0.5, 0.25, …} such that +// ||G(x + α·Δx)||₂ < ||G(x)||₂. Returns the accepted step (α may stay 1). +template +inline std::vector line_search( + const std::vector& x, + const Eigen::VectorXd& dx, + double norm0, + GradFn&& grad_fn, + int max_halvings = 20) +{ + const int n = static_cast(x.size()); + double alpha = 1.0; + std::vector xnew(static_cast(n)); + + for (int ls = 0; ls < max_halvings; ++ls) { + for (int i = 0; i < n; ++i) + xnew[static_cast(i)] = x[static_cast(i)] + + alpha * dx[i]; + auto Gnew = grad_fn(xnew); + double norm_new = 0.0; + for (double v : Gnew) norm_new += v * v; + norm_new = std::sqrt(norm_new); + if (norm_new < norm0) return xnew; + alpha *= 0.5; + } + // No improvement found — return best attempt (full step) + for (int i = 0; i < n; ++i) + xnew[static_cast(i)] = x[static_cast(i)] + dx[i]; + return xnew; +} + +} // namespace detail + +// ── Euclidean Newton solver ──────────────────────────────────────────────────── +// +// Minimises the Euclidean discrete conformal energy by solving G(x) = 0. +// The Hessian H is PSD; Eigen::SimplicialLDLT is used directly. +inline NewtonResult newton_euclidean( + ConformalMesh& mesh, + std::vector x0, + const EuclideanMaps& m, + double tol = 1e-8, + int max_iter = 200) +{ + std::vector x = x0; + const int n = static_cast(x.size()); + + NewtonResult res; + res.converged = false; + res.iterations = 0; + res.grad_inf_norm = 0.0; + + Eigen::SimplicialLDLT> solver; + + for (int iter = 0; iter < max_iter; ++iter) { + // ── Gradient ────────────────────────────────────────────────────────── + auto G_std = euclidean_gradient(mesh, x, m); + Eigen::Map G(G_std.data(), n); + + double inf_norm = G.cwiseAbs().maxCoeff(); + if (inf_norm < tol) { + res.converged = true; + res.grad_inf_norm = inf_norm; + res.iterations = iter; + res.x = x; + return res; + } + + // ── Hessian + factorisation ─────────────────────────────────────────── + auto H = euclidean_hessian(mesh, x, m); + solver.compute(H); + if (solver.info() != Eigen::Success) break; + + // ── Newton step: solve H·Δx = −G ──────────────────────────────────── + Eigen::VectorXd dx = solver.solve(-G); + if (solver.info() != Eigen::Success) break; + + // ── Backtracking line search ────────────────────────────────────────── + double norm0 = G.norm(); + x = detail::line_search(x, dx, norm0, + [&](const std::vector& xnew) { + return euclidean_gradient(mesh, xnew, m); + }); + + res.iterations = iter + 1; + } + + // Report final gradient norm + auto G_final = euclidean_gradient(mesh, x, m); + double inf_final = 0.0; + for (double v : G_final) inf_final = std::max(inf_final, std::abs(v)); + res.grad_inf_norm = inf_final; + res.x = x; + return res; +} + +// ── Spherical Newton solver ─────────────────────────────────────────────────── +// +// Solves G(x) = 0 for the spherical discrete conformal functional. +// The Hessian H is NSD at the solution; −H is PSD, so we factorise −H and +// solve (−H)·Δx = G ⟺ H·Δx = −G. +inline NewtonResult newton_spherical( + ConformalMesh& mesh, + std::vector x0, + const SphericalMaps& m, + double tol = 1e-8, + int max_iter = 200) +{ + std::vector x = x0; + const int n = static_cast(x.size()); + + NewtonResult res; + res.converged = false; + res.iterations = 0; + res.grad_inf_norm = 0.0; + + Eigen::SimplicialLDLT> solver; + + for (int iter = 0; iter < max_iter; ++iter) { + // ── Gradient ────────────────────────────────────────────────────────── + auto G_std = spherical_gradient(mesh, x, m); + Eigen::Map G(G_std.data(), n); + + double inf_norm = G.cwiseAbs().maxCoeff(); + if (inf_norm < tol) { + res.converged = true; + res.grad_inf_norm = inf_norm; + res.iterations = iter; + res.x = x; + return res; + } + + // ── Hessian: negate to get PSD matrix ──────────────────────────────── + auto H = spherical_hessian(mesh, x, m); + auto negH = Eigen::SparseMatrix(-H); + solver.compute(negH); + if (solver.info() != Eigen::Success) break; + + // ── Newton step: solve (−H)·Δx = G ───────────────────────────────── + Eigen::VectorXd dx = solver.solve(G); + if (solver.info() != Eigen::Success) break; + + // ── Backtracking line search ────────────────────────────────────────── + double norm0 = G.norm(); + x = detail::line_search(x, dx, norm0, + [&](const std::vector& xnew) { + return spherical_gradient(mesh, xnew, m); + }); + + res.iterations = iter + 1; + } + + auto G_final = spherical_gradient(mesh, x, m); + double inf_final = 0.0; + for (double v : G_final) inf_final = std::max(inf_final, std::abs(v)); + res.grad_inf_norm = inf_final; + res.x = x; + return res; +} + +} // namespace conformallab diff --git a/code/tests/cgal/CMakeLists.txt b/code/tests/cgal/CMakeLists.txt index a61a74c..656d714 100644 --- a/code/tests/cgal/CMakeLists.txt +++ b/code/tests/cgal/CMakeLists.txt @@ -24,6 +24,12 @@ add_executable(conformallab_cgal_tests # ── Phase 3f: Hessians (cotangent Laplacian, spherical + Euclidean) ─── test_euclidean_hessian.cpp test_spherical_hessian.cpp + + # ── Phase 4a: Newton solver ──────────────────────────────────────────── + test_newton_solver.cpp + + # ── Phase 4b: Mesh I/O (CGAL::IO) ───────────────────────────────────── + test_mesh_io.cpp ) target_include_directories(conformallab_cgal_tests SYSTEM PRIVATE diff --git a/code/tests/cgal/test_mesh_io.cpp b/code/tests/cgal/test_mesh_io.cpp new file mode 100644 index 0000000..a2c7649 --- /dev/null +++ b/code/tests/cgal/test_mesh_io.cpp @@ -0,0 +1,170 @@ +// test_mesh_io.cpp +// +// Phase 4b — CGAL::IO mesh round-trip tests. +// +// Tests: +// 1. OFF write + read round-trip: vertex count, face count, edge count preserved. +// 2. OBJ write + read round-trip: same topology checks. +// 3. load_mesh() throws on non-existent file. +// 4. Round-trip preserves vertex positions (within floating-point precision). + +#include "conformal_mesh.hpp" +#include "mesh_builder.hpp" +#include "mesh_io.hpp" +#include +#include +#include +#include +#include + +using namespace conformallab; + +// Helper: temp file path with given extension +static std::string tmp_path(const char* ext) +{ + return std::string("/tmp/conformallab_test.") + ext; +} + +// Helper: delete file if it exists +static void rm(const std::string& path) +{ + std::filesystem::remove(path); +} + +// ════════════════════════════════════════════════════════════════════════════ +// OFF round-trip: tetrahedron +// ════════════════════════════════════════════════════════════════════════════ + +TEST(MeshIO, OFF_RoundTrip_Tetrahedron) +{ + auto path = tmp_path("off"); + rm(path); + + auto mesh_out = make_tetrahedron(); + ASSERT_TRUE(write_mesh(path, mesh_out)) + << "write_mesh should succeed for tetrahedron"; + + ConformalMesh mesh_in; + ASSERT_TRUE(read_mesh(path, mesh_in)) + << "read_mesh should succeed for the written OFF file"; + + EXPECT_EQ(mesh_in.number_of_vertices(), mesh_out.number_of_vertices()); + EXPECT_EQ(mesh_in.number_of_faces(), mesh_out.number_of_faces()); + EXPECT_EQ(mesh_in.number_of_edges(), mesh_out.number_of_edges()); + + rm(path); +} + +// ════════════════════════════════════════════════════════════════════════════ +// OFF round-trip: quad strip +// ════════════════════════════════════════════════════════════════════════════ + +TEST(MeshIO, OFF_RoundTrip_QuadStrip) +{ + auto path = tmp_path("off"); + rm(path); + + auto mesh_out = make_quad_strip(); + ASSERT_TRUE(write_mesh(path, mesh_out)); + + ConformalMesh mesh_in; + ASSERT_TRUE(read_mesh(path, mesh_in)); + + EXPECT_EQ(mesh_in.number_of_vertices(), mesh_out.number_of_vertices()); + EXPECT_EQ(mesh_in.number_of_faces(), mesh_out.number_of_faces()); + + rm(path); +} + +// ════════════════════════════════════════════════════════════════════════════ +// OBJ round-trip: tetrahedron +// ════════════════════════════════════════════════════════════════════════════ + +TEST(MeshIO, OBJ_RoundTrip_Tetrahedron) +{ + auto path = tmp_path("obj"); + rm(path); + + auto mesh_out = make_tetrahedron(); + ASSERT_TRUE(write_mesh(path, mesh_out)) + << "write_mesh (OBJ) should succeed"; + + ConformalMesh mesh_in; + ASSERT_TRUE(read_mesh(path, mesh_in)) + << "read_mesh (OBJ) should succeed"; + + EXPECT_EQ(mesh_in.number_of_vertices(), mesh_out.number_of_vertices()); + EXPECT_EQ(mesh_in.number_of_faces(), mesh_out.number_of_faces()); + + rm(path); +} + +// ════════════════════════════════════════════════════════════════════════════ +// load_mesh throws on missing file +// ════════════════════════════════════════════════════════════════════════════ + +TEST(MeshIO, LoadMeshThrowsOnMissingFile) +{ + EXPECT_THROW(load_mesh("/tmp/does_not_exist_conflab.off"), std::runtime_error); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Vertex positions survive a round-trip (OFF) +// ════════════════════════════════════════════════════════════════════════════ + +TEST(MeshIO, OFF_VertexPositionsPreserved) +{ + auto path = tmp_path("off"); + rm(path); + + auto mesh_out = make_tetrahedron(); + ASSERT_TRUE(write_mesh(path, mesh_out)); + + ConformalMesh mesh_in; + ASSERT_TRUE(read_mesh(path, mesh_in)); + + // Collect and sort vertex positions from both meshes to compare + auto collect_sorted = [](const ConformalMesh& m) { + std::vector> pts; + for (auto v : m.vertices()) { + auto p = m.point(v); + pts.push_back({CGAL::to_double(p.x()), + CGAL::to_double(p.y()), + CGAL::to_double(p.z())}); + } + std::sort(pts.begin(), pts.end()); + return pts; + }; + + auto pts_out = collect_sorted(mesh_out); + auto pts_in = collect_sorted(mesh_in); + + ASSERT_EQ(pts_out.size(), pts_in.size()); + for (std::size_t i = 0; i < pts_out.size(); ++i) { + EXPECT_NEAR(pts_out[i][0], pts_in[i][0], 1e-10); + EXPECT_NEAR(pts_out[i][1], pts_in[i][1], 1e-10); + EXPECT_NEAR(pts_out[i][2], pts_in[i][2], 1e-10); + } + + rm(path); +} + +// ════════════════════════════════════════════════════════════════════════════ +// save_mesh / load_mesh convenience wrappers +// ════════════════════════════════════════════════════════════════════════════ + +TEST(MeshIO, SaveLoadConvenienceWrappers) +{ + auto path = tmp_path("off"); + rm(path); + + auto mesh_out = make_quad_strip(); + EXPECT_NO_THROW(save_mesh(path, mesh_out)); + + ConformalMesh mesh_in; + EXPECT_NO_THROW(mesh_in = load_mesh(path)); + + EXPECT_EQ(mesh_in.number_of_vertices(), mesh_out.number_of_vertices()); + + rm(path); +} diff --git a/code/tests/cgal/test_newton_solver.cpp b/code/tests/cgal/test_newton_solver.cpp new file mode 100644 index 0000000..4475bb1 --- /dev/null +++ b/code/tests/cgal/test_newton_solver.cpp @@ -0,0 +1,235 @@ +// test_newton_solver.cpp +// +// Phase 4a — Newton solver tests. +// +// Design principle: +// We test convergence to a KNOWN equilibrium. For the spherical tetrahedron +// x* = 0 is built-in (G(0) ≈ 0 by construction). For Euclidean meshes we +// use "natural theta": set theta_v[v] = actual angle sum at x=0, which makes +// x* = 0 the exact equilibrium by definition. +// +// Tests: +// Spherical: +// 1. Converges from x=[-0.2,...] to x*=0 (spherical tetrahedron). +// 2. Converges in few iterations (quadratic convergence near equilibrium). +// 3. Converges from a large perturbation x=[-0.5,...]. +// 4. Result fields (x.size, grad_inf_norm) are self-consistent. +// +// Euclidean: +// 5. Converges (triangle, 1 pinned vertex, natural theta). +// 6. Converges (quad strip, 1 pinned vertex, natural theta). +// 7. Converges with explicitly chosen mixed pinned/variable layout. + +#include "conformal_mesh.hpp" +#include "mesh_builder.hpp" +#include "euclidean_functional.hpp" +#include "spherical_functional.hpp" +#include "newton_solver.hpp" +#include +#include +#include + +using namespace conformallab; + +// ──────────────────────────────────────────────────────────────────────────── +// Helper: set theta_v[v] = actual angle sum at x=0 for each variable vertex. +// +// Requires that DOF indices have already been assigned (v_idx populated). +// Uses euclidean_gradient directly so the formula is exact and consistent. +// ──────────────────────────────────────────────────────────────────────────── +static void set_natural_euclidean_theta(ConformalMesh& mesh, EuclideanMaps& maps, int n) +{ + std::vector x0(static_cast(n), 0.0); + // G[iv] = theta_v[v] - sum_alpha(x=0) + // so sum_alpha(x=0) = theta_v[v] - G[iv] + auto G = euclidean_gradient(mesh, x0, maps); + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv < 0) continue; + maps.theta_v[v] -= G[static_cast(iv)]; + } +} + +// ════════════════════════════════════════════════════════════════════════════ +// Spherical 1 — Converges from moderate perturbation +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonSolver, Spherical_ConvergesFromPerturbation) +{ + auto mesh = make_spherical_tetrahedron(); + auto maps = setup_spherical_maps(mesh); + compute_lambda0_from_mesh(mesh, maps); + int n = assign_vertex_dof_indices(mesh, maps); + + std::vector x0(static_cast(n), -0.2); + auto res = newton_spherical(mesh, x0, maps, /*tol=*/1e-8, /*max_iter=*/50); + + EXPECT_TRUE(res.converged) + << "Newton (spherical) should converge; grad_inf_norm = " << res.grad_inf_norm; + EXPECT_LT(res.grad_inf_norm, 1e-8); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Spherical 2 — Quadratic convergence: few iterations suffice +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonSolver, Spherical_FewIterations) +{ + auto mesh = make_spherical_tetrahedron(); + auto maps = setup_spherical_maps(mesh); + compute_lambda0_from_mesh(mesh, maps); + int n = assign_vertex_dof_indices(mesh, maps); + + std::vector x0(static_cast(n), -0.2); + auto res = newton_spherical(mesh, x0, maps, /*tol=*/1e-8, /*max_iter=*/50); + + EXPECT_LE(res.iterations, 20) + << "Newton should converge in ≤ 20 iterations; took " << res.iterations; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Spherical 3 — Large perturbation: global convergence via line search +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonSolver, Spherical_ConvergesFromLargePerturbation) +{ + auto mesh = make_spherical_tetrahedron(); + auto maps = setup_spherical_maps(mesh); + compute_lambda0_from_mesh(mesh, maps); + int n = assign_vertex_dof_indices(mesh, maps); + + std::vector x0(static_cast(n), -0.5); + auto res = newton_spherical(mesh, x0, maps, /*tol=*/1e-8, /*max_iter=*/100); + + EXPECT_TRUE(res.converged) + << "Newton (spherical, large perturbation) should converge; " + "grad_inf_norm = " << res.grad_inf_norm; + EXPECT_LT(res.grad_inf_norm, 1e-8); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Spherical 4 — Result fields are self-consistent +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonSolver, Spherical_ResultFieldsConsistent) +{ + auto mesh = make_spherical_tetrahedron(); + auto maps = setup_spherical_maps(mesh); + compute_lambda0_from_mesh(mesh, maps); + int n = assign_vertex_dof_indices(mesh, maps); + + std::vector x0(static_cast(n), -0.1); + auto res = newton_spherical(mesh, x0, maps, /*tol=*/1e-8); + + EXPECT_EQ(static_cast(res.x.size()), n); + + // Reported grad_inf_norm must match re-computed gradient at res.x + auto G = spherical_gradient(mesh, res.x, maps); + double actual_inf = 0.0; + for (double v : G) actual_inf = std::max(actual_inf, std::abs(v)); + EXPECT_NEAR(actual_inf, res.grad_inf_norm, 1e-10); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Euclidean 1 — Triangle with 1 pinned vertex + natural theta +// +// make_triangle(): 3 vertices. Pin v0 → 2 free DOFs. +// With natural theta: x* = 0 (G(0) = 0 by construction). +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonSolver, Euclidean_ConvergesTrianglePinned) +{ + auto mesh = make_triangle(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + + // Pin first vertex, assign sequential DOFs to the other two + auto vit = mesh.vertices().begin(); + Vertex_index v0 = *vit++; + maps.v_idx[v0] = -1; // pinned + int idx = 0; + for (; vit != mesh.vertices().end(); ++vit) + maps.v_idx[*vit] = idx++; + int n = idx; // = 2 + + set_natural_euclidean_theta(mesh, maps, n); + + std::vector x0(static_cast(n), -0.1); + auto res = newton_euclidean(mesh, x0, maps, /*tol=*/1e-8, /*max_iter=*/50); + + EXPECT_TRUE(res.converged) + << "Newton (Euclidean, triangle, pinned) should converge; " + "grad_inf_norm = " << res.grad_inf_norm; + EXPECT_LT(res.grad_inf_norm, 1e-8); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Euclidean 2 — Quad strip with 1 pinned vertex + natural theta +// +// make_quad_strip(): 4 vertices, 2 faces. Pin v0 → 3 free DOFs. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonSolver, Euclidean_ConvergesQuadStripPinned) +{ + auto mesh = make_quad_strip(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + + // Pin first vertex + auto vit = mesh.vertices().begin(); + Vertex_index v0 = *vit++; + maps.v_idx[v0] = -1; + int idx = 0; + for (; vit != mesh.vertices().end(); ++vit) + maps.v_idx[*vit] = idx++; + int n = idx; // = 3 + + set_natural_euclidean_theta(mesh, maps, n); + + std::vector x0(static_cast(n), -0.15); + auto res = newton_euclidean(mesh, x0, maps, /*tol=*/1e-8, /*max_iter=*/50); + + EXPECT_TRUE(res.converged) + << "Newton (Euclidean, quad strip, pinned) should converge; " + "grad_inf_norm = " << res.grad_inf_norm; + EXPECT_LT(res.grad_inf_norm, 1e-8); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Euclidean 3 — Mixed pinned layout: explicit vertex assignment +// +// Quad strip: v0 pinned, v1/v2/v3 free. +// Natural theta set AFTER DOF assignment so that x* = 0 is the equilibrium +// for the free vertices (with v0 fixed at u0=0). +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonSolver, Euclidean_ConvergesMixedPinned) +{ + auto mesh = make_quad_strip(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + + // Explicitly assign DOF indices + auto vit = mesh.vertices().begin(); + Vertex_index v0 = *vit++; + Vertex_index v1 = *vit++; + Vertex_index v2 = *vit++; + Vertex_index v3 = *vit; + + maps.v_idx[v0] = -1; // pinned at u0 = 0 + maps.v_idx[v1] = 0; + maps.v_idx[v2] = 1; + maps.v_idx[v3] = 2; + const int n = 3; + + // Set natural theta AFTER pinning so that x* = [0,0,0] is the equilibrium + set_natural_euclidean_theta(mesh, maps, n); + + std::vector x0 = {-0.1, -0.15, -0.05}; + auto res = newton_euclidean(mesh, x0, maps, /*tol=*/1e-8, /*max_iter=*/50); + + EXPECT_TRUE(res.converged) + << "Newton (Euclidean, mixed pinned) should converge; " + "grad_inf_norm = " << res.grad_inf_norm; + EXPECT_LT(res.grad_inf_norm, 1e-8); +}