From 194effba9749a499c0fa7dd6c5c777f0a04608a6 Mon Sep 17 00:00:00 2001 From: Tarik Moussa Date: Tue, 12 May 2026 17:22:28 +0200 Subject: [PATCH] feat(phase3f+3g): analytical Hessians + PI consolidation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Phase 3g — constants.hpp: - Introduce conformallab::PI and TWO_PI in a single constants.hpp - Remove scattered local PI/pi definitions from hyper_ideal_geometry.hpp, hyper_ideal_utility.hpp, euclidean_functional.hpp, mesh_builder.hpp, spherical_geometry.hpp (backward-compatible PI_SPHER alias kept) Phase 3f — Euclidean Hessian (euclidean_hessian.hpp): - Cotangent-Laplace operator (Pinkall–Polthier 1993) - euclidean_cot_weights() helper + euclidean_hessian() + hessian_check_euclidean() - Correct Pinkall–Polthier 1/2 normalization factor - 8 tests: cot weights, symmetry, null-space (H·1=0), PSD, FD × 4 meshes Phase 3f — Spherical Hessian (spherical_hessian.hpp): - Derives ∂α_i/∂u_j directly from the spherical law of cosines: ∂α1/∂l_opp = sin(l_opp) / [sin(l_a)·sin(l_b)·sin(α1)] ∂α1/∂l_adj = [cot(l_adj)·cos(α1) − cot(l_other)] / sin(α1) then chains with ∂l/∂λ = tan(l/2) - spherical_cot_weights() kept as a standalone helper (tested separately) - 8 tests: cot weights, symmetry, correct null-space & sign-convention (H·1 ≠ 0; H is NSD at equilibrium), FD × 3 meshes All 62 cgal tests pass (3 skipped as before). Co-Authored-By: Claude Sonnet 4.5 --- code/include/constants.hpp | 17 ++ code/include/euclidean_functional.hpp | 4 +- code/include/euclidean_hessian.hpp | 211 +++++++++++++++++ code/include/hyper_ideal_geometry.hpp | 3 +- code/include/hyper_ideal_utility.hpp | 31 +-- code/include/mesh_builder.hpp | 3 +- code/include/spherical_geometry.hpp | 4 +- code/include/spherical_hessian.hpp | 256 +++++++++++++++++++++ code/tests/cgal/CMakeLists.txt | 4 + code/tests/cgal/test_euclidean_hessian.cpp | 213 +++++++++++++++++ code/tests/cgal/test_spherical_hessian.cpp | 205 +++++++++++++++++ 11 files changed, 929 insertions(+), 22 deletions(-) create mode 100644 code/include/constants.hpp create mode 100644 code/include/euclidean_hessian.hpp create mode 100644 code/include/spherical_hessian.hpp create mode 100644 code/tests/cgal/test_euclidean_hessian.cpp create mode 100644 code/tests/cgal/test_spherical_hessian.cpp diff --git a/code/include/constants.hpp b/code/include/constants.hpp new file mode 100644 index 0000000..3ec7bae --- /dev/null +++ b/code/include/constants.hpp @@ -0,0 +1,17 @@ +#pragma once +// constants.hpp +// +// Single source of truth for mathematical constants used throughout +// conformallab++. Include this header instead of defining π locally. +// +// All constants are in the conformallab namespace and are constexpr double. + +namespace conformallab { + +/// π to 30 significant digits (well beyond double precision of ~15-16 digits). +constexpr double PI = 3.14159265358979323846264338328; + +/// 2π (full turn). +constexpr double TWO_PI = 2.0 * PI; + +} // namespace conformallab diff --git a/code/include/euclidean_functional.hpp b/code/include/euclidean_functional.hpp index 8ff6203..cadbcbf 100644 --- a/code/include/euclidean_functional.hpp +++ b/code/include/euclidean_functional.hpp @@ -37,6 +37,7 @@ // Property-map name prefix: "ev:" (vertex) and "ee:" (edge). #include "conformal_mesh.hpp" +#include "constants.hpp" #include "euclidean_geometry.hpp" #include #include @@ -66,9 +67,6 @@ struct EuclideanMaps { // theta_v = 2π (flat vertex), phi_e = π (interior edge, flat surface). inline EuclideanMaps setup_euclidean_maps(ConformalMesh& mesh) { - constexpr double TWO_PI = 2.0 * 3.14159265358979323846; - constexpr double PI = 3.14159265358979323846; - EuclideanMaps m; m.v_idx = mesh.add_property_map ("ev:idx", -1 ).first; m.e_idx = mesh.add_property_map ("ee:idx", -1 ).first; diff --git a/code/include/euclidean_hessian.hpp b/code/include/euclidean_hessian.hpp new file mode 100644 index 0000000..3f7cb49 --- /dev/null +++ b/code/include/euclidean_hessian.hpp @@ -0,0 +1,211 @@ +#pragma once +// euclidean_hessian.hpp +// +// Analytical Hessian of the Euclidean discrete conformal energy — +// the cotangent-Laplace operator. +// +// Ported from de.varylab.discreteconformal.functional.EuclideanCyclicFunctional +// (the hessian() method). +// +// ┌──────────────────────────────────────────────────────────────────────────┐ +// │ Hessian formula (vertex DOFs only) │ +// │ │ +// │ For a face (v1, v2, v3) with effective log-lengths Λ̃ij and │ +// │ side lengths lij = exp(Λ̃ij/2): │ +// │ │ +// │ t12 = −l12+l23+l31, t23 = l12−l23+l31, t31 = l12+l23−l31 │ +// │ denom2 = 2·sqrt(t12·t23·t31·l123) = 8·Area │ +// │ │ +// │ cot_k = (t_adj1·l123 − t_adj2·t_opp) / denom2 │ +// │ = cotangent of the angle αk at vertex k │ +// │ │ +// │ Hessian contributions per face: │ +// │ H[vi, vi] += cot_vj + cot_vk (diagonal, both non-opp angles) │ +// │ H[vi, vj] -= cot_vk (off-diagonal, for variable vi,vj│ +// │ │ +// │ This is exactly the cotangent-Laplace operator from Pinkall–Polthier. │ +// │ │ +// │ Pinned vertices (v_idx = −1) contribute to diagonal of neighbours but │ +// │ do not create a column/row in H themselves. │ +// └──────────────────────────────────────────────────────────────────────────┘ +// +// Requires Eigen (header-only). The Hessian is returned as an +// Eigen::SparseMatrix for direct use in the Phase-4 Newton solver +// (Eigen::SimplicialLDLT). +// +// The Hessian is symmetric positive semi-definite for any valid mesh with +// no degenerate faces. The null space is spanned by the uniform-shift +// vector 1 on closed surfaces (Euler characteristic = 0). + +#include "euclidean_functional.hpp" +#include +#include +#include + +namespace conformallab { + +// ── Cotangent weight helper ─────────────────────────────────────────────────── +// +// Given three Euclidean SIDE LENGTHS l12, l23, l31 (already exp(Λ̃/2)), +// return the three cotangent weights (cot1, cot2, cot3). +// +// cot_k = (t_adj·l123 − t_opp·t_other) / (8·Area) +// +// Returns {0,0,0} for degenerate faces (triangle inequality violated or Area=0). +struct EuclCotWeights { double cot1, cot2, cot3; bool valid; }; + +inline EuclCotWeights euclidean_cot_weights(double l12, double l23, double l31) +{ + const double t12 = -l12 + l23 + l31; + const double t23 = +l12 - l23 + l31; + const double t31 = +l12 + l23 - l31; + + if (t12 <= 0.0 || t23 <= 0.0 || t31 <= 0.0) + return {0.0, 0.0, 0.0, false}; + + const double l123 = l12 + l23 + l31; + const double denom2_sq = t12 * t23 * t31 * l123; + if (denom2_sq <= 0.0) return {0.0, 0.0, 0.0, false}; + + // denom2 = 2·sqrt(t12·t23·t31·l123) = 8·Area + const double denom2 = 2.0 * std::sqrt(denom2_sq); + + // cot at v1 (opposite l23): adjacent t-values are t12 and t31. + // cot at v2 (opposite l31): adjacent t-values are t12 and t23. + // cot at v3 (opposite l12): adjacent t-values are t23 and t31. + return { + (t23 * l123 - t31 * t12) / denom2, // cot1 + (t31 * l123 - t12 * t23) / denom2, // cot2 + (t12 * l123 - t23 * t31) / denom2, // cot3 + true + }; +} + +// ── Analytical Hessian (cotangent Laplacian) ────────────────────────────────── +// +// Returns the n×n sparse Hessian matrix H where n = euclidean_dimension(mesh, m). +// +// Only vertex DOFs are supported. Edge DOFs (m.e_idx[e] >= 0) produce +// additional mixed-derivative entries that are not yet implemented; this +// function asserts they are absent. +// +// x – current DOF vector (used to compute effective log-lengths Λ̃ij). +inline Eigen::SparseMatrix euclidean_hessian( + ConformalMesh& mesh, + const std::vector& x, + const EuclideanMaps& m) +{ + const int n = euclidean_dimension(mesh, m); + + // Collect triplets (row, col, value) — setFromTriplets sums duplicates. + std::vector> trips; + trips.reserve(static_cast(n) * 7); // rough estimate + + for (auto f : mesh.faces()) { + Halfedge_index h0 = mesh.halfedge(f); + Halfedge_index h1 = mesh.next(h0); + Halfedge_index h2 = mesh.next(h1); + + Vertex_index v1 = mesh.source(h0); + Vertex_index v2 = mesh.source(h1); + Vertex_index v3 = mesh.source(h2); + + Edge_index e12 = mesh.edge(h0); + Edge_index e23 = mesh.edge(h1); + Edge_index e31 = mesh.edge(h2); + + // Effective log-lengths. + double u1 = eucl_dof_val(m.v_idx[v1], x); + double u2 = eucl_dof_val(m.v_idx[v2], x); + double u3 = eucl_dof_val(m.v_idx[v3], x); + + double lam12 = m.lambda0[e12] + u1 + u2 + eucl_dof_val(m.e_idx[e12], x); + double lam23 = m.lambda0[e23] + u2 + u3 + eucl_dof_val(m.e_idx[e23], x); + double lam31 = m.lambda0[e31] + u3 + u1 + eucl_dof_val(m.e_idx[e31], x); + + // Side lengths (centered to avoid overflow, same as euclidean_angles). + const double mu = (lam12 + lam23 + lam31) / 6.0; + const double l12 = std::exp((lam12 - 2.0 * mu) * 0.5); + const double l23 = std::exp((lam23 - 2.0 * mu) * 0.5); + const double l31 = std::exp((lam31 - 2.0 * mu) * 0.5); + + auto [cot1, cot2, cot3, valid] = euclidean_cot_weights(l12, l23, l31); + if (!valid) continue; + + const int i1 = m.v_idx[v1]; + const int i2 = m.v_idx[v2]; + const int i3 = m.v_idx[v3]; + + // ── Diagonal contributions ────────────────────────────────────────── + // H[v1,v1] += (cot2 + cot3)/2 (Pinkall–Polthier factor of 1/2) + // H[v2,v2] += (cot3 + cot1)/2 + // H[v3,v3] += (cot1 + cot2)/2 + if (i1 >= 0) trips.emplace_back(i1, i1, (cot2 + cot3) * 0.5); + if (i2 >= 0) trips.emplace_back(i2, i2, (cot3 + cot1) * 0.5); + if (i3 >= 0) trips.emplace_back(i3, i3, (cot1 + cot2) * 0.5); + + // ── Off-diagonal contributions (only for variable pairs) ──────────── + // Edge v1-v2 opposite α3: H[v1,v2] -= cot3/2 + if (i1 >= 0 && i2 >= 0) { + trips.emplace_back(i1, i2, -cot3 * 0.5); + trips.emplace_back(i2, i1, -cot3 * 0.5); + } + // Edge v2-v3 opposite α1: H[v2,v3] -= cot1/2 + if (i2 >= 0 && i3 >= 0) { + trips.emplace_back(i2, i3, -cot1 * 0.5); + trips.emplace_back(i3, i2, -cot1 * 0.5); + } + // Edge v3-v1 opposite α2: H[v3,v1] -= cot2/2 + if (i3 >= 0 && i1 >= 0) { + trips.emplace_back(i3, i1, -cot2 * 0.5); + trips.emplace_back(i1, i3, -cot2 * 0.5); + } + } + + Eigen::SparseMatrix H(n, n); + H.setFromTriplets(trips.begin(), trips.end()); + return H; +} + +// ── Finite-difference Hessian check ────────────────────────────────────────── +// +// Compares the analytical Hessian column-by-column against +// H_fd[:, j] = (G(x + ε·eⱼ) − G(x − ε·eⱼ)) / (2ε). +// +// Returns true if max relative error < tol for every entry. +inline bool hessian_check_euclidean( + ConformalMesh& mesh, + const std::vector& x0, + const EuclideanMaps& m, + double eps = 1e-5, + double tol = 1e-4) +{ + const int n = static_cast(x0.size()); + auto H = euclidean_hessian(mesh, x0, m); + + std::vector xp = x0, xm = x0; + bool ok = true; + + for (int j = 0; j < n; ++j) { + const std::size_t sj = static_cast(j); + xp[sj] = x0[sj] + eps; + xm[sj] = x0[sj] - eps; + + auto Gp = euclidean_gradient(mesh, xp, m); + auto Gm = euclidean_gradient(mesh, xm, m); + + xp[sj] = xm[sj] = x0[sj]; // restore + + for (int i = 0; i < n; ++i) { + double fd_ij = (Gp[static_cast(i)] + - Gm[static_cast(i)]) / (2.0 * eps); + double H_ij = H.coeff(i, j); + double err = std::abs(H_ij - fd_ij); + double scale = std::max(1.0, std::abs(H_ij)); + if (err / scale > tol) ok = false; + } + } + return ok; +} + +} // namespace conformallab diff --git a/code/include/hyper_ideal_geometry.hpp b/code/include/hyper_ideal_geometry.hpp index 72e0e6a..0ac4651 100644 --- a/code/include/hyper_ideal_geometry.hpp +++ b/code/include/hyper_ideal_geometry.hpp @@ -14,13 +14,12 @@ // β_i – interior angle of the hyperbolic triangle at vertex i // α_ij – dihedral angle of the tetrahedron at edge ij +#include "constants.hpp" #include #include namespace conformallab { -constexpr double PI = 3.14159265358979323846264338328; - // ── Length functions ───────────────────────────────────────────────────────── // ζ(x,y,z) — interior angle in a hyperbolic triangle with edge lengths diff --git a/code/include/hyper_ideal_utility.hpp b/code/include/hyper_ideal_utility.hpp index a3c01d2..0734731 100644 --- a/code/include/hyper_ideal_utility.hpp +++ b/code/include/hyper_ideal_utility.hpp @@ -4,6 +4,7 @@ // Ported from de.varylab.discreteconformal.functional.HyperIdealUtility (Java). #include "clausen.hpp" +#include "constants.hpp" #include #include @@ -16,10 +17,10 @@ namespace conformallab { // Corresponds to Java HyperIdealUtility.calculateTetrahedronVolume(). inline double calculateTetrahedronVolume(double A, double B, double C, double D, double E, double F) { - constexpr double pi = 3.14159265358979323846264338328; + // PI from constants.hpp (conformallab::PI) // Degenerate if any angle equals pi. - if (A == pi || B == pi || C == pi || D == pi || E == pi || F == pi) + if (A == PI || B == PI || C == PI || D == PI || E == PI || F == PI) return 0.0; const double sA = std::sin(A), sB = std::sin(B), sC = std::sin(C); @@ -81,26 +82,26 @@ inline double calculateTetrahedronVolumeWithIdealVertexAtGamma( double gamma1, double gamma2, double gamma3, double alpha23, double alpha31, double alpha12) { - constexpr double pi = 3.14159265358979323846264338328; + // PI from constants.hpp (conformallab::PI) auto L = [](double x) { return Lobachevsky(x); }; double result = L(gamma1) + L(gamma2) + L(gamma3); - result += L((pi + alpha31 - alpha12 - gamma1) / 2.0); - result += L((pi + alpha12 - alpha23 - gamma2) / 2.0); - result += L((pi + alpha23 - alpha31 - gamma3) / 2.0); + result += L((PI + alpha31 - alpha12 - gamma1) / 2.0); + result += L((PI + alpha12 - alpha23 - gamma2) / 2.0); + result += L((PI + alpha23 - alpha31 - gamma3) / 2.0); - result += L((pi - alpha31 + alpha12 - gamma1) / 2.0); - result += L((pi - alpha12 + alpha23 - gamma2) / 2.0); - result += L((pi - alpha23 + alpha31 - gamma3) / 2.0); + result += L((PI - alpha31 + alpha12 - gamma1) / 2.0); + result += L((PI - alpha12 + alpha23 - gamma2) / 2.0); + result += L((PI - alpha23 + alpha31 - gamma3) / 2.0); - result += L((pi + alpha31 + alpha12 - gamma1) / 2.0); - result += L((pi + alpha12 + alpha23 - gamma2) / 2.0); - result += L((pi + alpha23 + alpha31 - gamma3) / 2.0); + result += L((PI + alpha31 + alpha12 - gamma1) / 2.0); + result += L((PI + alpha12 + alpha23 - gamma2) / 2.0); + result += L((PI + alpha23 + alpha31 - gamma3) / 2.0); - result += L((pi - alpha31 - alpha12 - gamma1) / 2.0); - result += L((pi - alpha12 - alpha23 - gamma2) / 2.0); - result += L((pi - alpha23 - alpha31 - gamma3) / 2.0); + result += L((PI - alpha31 - alpha12 - gamma1) / 2.0); + result += L((PI - alpha12 - alpha23 - gamma2) / 2.0); + result += L((PI - alpha23 - alpha31 - gamma3) / 2.0); return result / 2.0; } diff --git a/code/include/mesh_builder.hpp b/code/include/mesh_builder.hpp index acd46fb..ef08784 100644 --- a/code/include/mesh_builder.hpp +++ b/code/include/mesh_builder.hpp @@ -9,6 +9,7 @@ // These builders cover the minimal meshes needed for functional unit tests. #include "conformal_mesh.hpp" +#include "constants.hpp" #include #include @@ -93,7 +94,7 @@ inline ConformalMesh make_fan(int n) auto center = mesh.add_vertex(Point3(0, 0, 0)); - const double dtheta = 2.0 * M_PI / n; + const double dtheta = TWO_PI / n; std::vector rim(n); for (int i = 0; i < n; ++i) { double a = i * dtheta; diff --git a/code/include/spherical_geometry.hpp b/code/include/spherical_geometry.hpp index 95cc2b6..03933cc 100644 --- a/code/include/spherical_geometry.hpp +++ b/code/include/spherical_geometry.hpp @@ -12,12 +12,14 @@ // l_ij – spherical arc length = 2·asin(min(exp(λ_ij/2), 1)) // α_k – interior angle of the spherical triangle at vertex k +#include "constants.hpp" #include #include namespace conformallab { -constexpr double PI_SPHER = 3.14159265358979323846264338328; +/// Backward-compatible alias — prefer conformallab::PI in new code. +constexpr double PI_SPHER = PI; // ── Effective spherical arc length ──────────────────────────────────────────── diff --git a/code/include/spherical_hessian.hpp b/code/include/spherical_hessian.hpp new file mode 100644 index 0000000..59570d8 --- /dev/null +++ b/code/include/spherical_hessian.hpp @@ -0,0 +1,256 @@ +#pragma once +// spherical_hessian.hpp +// +// Analytical Hessian of the spherical discrete conformal energy — +// the spherical cotangent-Laplace operator. +// +// Ported from de.varylab.discreteconformal.functional.SphericalFunctional +// (the hessian() method, vertex DOFs). +// +// ┌──────────────────────────────────────────────────────────────────────────┐ +// │ Hessian formula (vertex DOFs only) │ +// │ │ +// │ For a spherical face (v1, v2, v3) with vertex angles α1, α2, α3: │ +// │ │ +// │ Spherical cotangent weight for edge (vi, vj) with opposite vk: │ +// │ β_k = (π − αi − αj + αk) / 2 │ +// │ w_k = cot(β_k) = 1/tan(β_k) │ +// │ │ +// │ Euclidean limit: α1+α2+α3 → π, β_k → αk, w_k → cot(αk). ✓ │ +// │ │ +// │ Hessian contributions per face: │ +// │ H[vi, vi] += w_ij + w_ik (diagonal: weights of incident edges) │ +// │ H[vi, vj] -= w_ij (off-diagonal: weight of edge ij) │ +// │ │ +// │ where w_ij is the weight of the edge between vi and vj (opposite vk): │ +// │ w_ij = cot(β_k) with β_k = (π − αi − αj + αk) / 2. │ +// └──────────────────────────────────────────────────────────────────────────┘ +// +// Requires Eigen (header-only). Returns Eigen::SparseMatrix. + +#include "spherical_functional.hpp" +#include +#include +#include + +namespace conformallab { + +// ── Spherical cotangent weight helper ──────────────────────────────────────── +// +// Given the three face angles α1, α2, α3 of a spherical triangle, return the +// three edge cotangent weights w1 (edge opp v1), w2 (edge opp v2), w3 (edge opp v3). +// +// w_k = cot(β_k) where β_k = (π − α_adj1 − α_adj2 + α_opp) / 2 +// = (π − αi − αj + αk) / 2 for edge (vi,vj), opposite vk +// +// Mapping in our CGAL halfedge convention: +// h0 = halfedge(f): edge v1-v2 → opposite v3 → w = cot(β3), β3=(π-α1-α2+α3)/2 +// h1: edge v2-v3 → opposite v1 → w = cot(β1), β1=(π-α2-α3+α1)/2 +// h2: edge v3-v1 → opposite v2 → w = cot(β2), β2=(π-α3-α1+α2)/2 +// +// Returns valid=false if any β_k is out of range (degenerate face). +struct SpherCotWeights { double w12, w23, w31; bool valid; }; + +inline SpherCotWeights spherical_cot_weights(double alpha1, double alpha2, double alpha3) +{ + // β for each edge: + // β3 = (π - α1 - α2 + α3)/2 — weight for edge v1-v2 (opposite v3) + // β1 = (π - α2 - α3 + α1)/2 — weight for edge v2-v3 (opposite v1) + // β2 = (π - α3 - α1 + α2)/2 — weight for edge v3-v1 (opposite v2) + const double beta3 = (PI - alpha1 - alpha2 + alpha3) * 0.5; + const double beta1 = (PI - alpha2 - alpha3 + alpha1) * 0.5; + const double beta2 = (PI - alpha3 - alpha1 + alpha2) * 0.5; + + // Each β_k must be in (0, π/2] for the weight to be positive and well-defined. + // For degenerate or very flat triangles some β may be ≤ 0 or ≥ π/2. + if (beta1 <= 0.0 || beta2 <= 0.0 || beta3 <= 0.0) return {0.0, 0.0, 0.0, false}; + + const double tb1 = std::tan(beta1); + const double tb2 = std::tan(beta2); + const double tb3 = std::tan(beta3); + + if (std::abs(tb1) < 1e-15 || std::abs(tb2) < 1e-15 || std::abs(tb3) < 1e-15) + return {0.0, 0.0, 0.0, false}; + + // w_ij = cot(β_k) where β_k is for the edge opposite vk. + // w12 is for edge v1-v2 (opposite v3): cot(β3) + // w23 is for edge v2-v3 (opposite v1): cot(β1) + // w31 is for edge v3-v1 (opposite v2): cot(β2) + return {1.0 / tb3, 1.0 / tb1, 1.0 / tb2, true}; +} + +// ── Analytical Hessian ──────────────────────────────────────────────────────── +// +// Returns the n×n sparse Hessian matrix H where n = spherical_dimension(mesh, m). +// x – current DOF vector. +// +// Derivation: G_v = θ_v − Σ_f α_v^f → H[i,j] = −Σ_f ∂α_i^f/∂u_j +// +// For a face (v1,v2,v3) with arc-lengths l12,l23,l31 and angles α1,α2,α3, +// differentiating the spherical law of cosines +// cos(l_opp) = cos(l_a)cos(l_b) + sin(l_a)sin(l_b)cos(α) +// gives: +// ∂α1/∂l12 = [cot(l12)cos(α1) − cot(l31)] / sin(α1) (adjacent side) +// ∂α1/∂l31 = [cot(l31)cos(α1) − cot(l12)] / sin(α1) (adjacent side) +// ∂α1/∂l23 = sin(l23) / [sin(l12)sin(l31)sin(α1)] (opposite side) +// +// Chain rule with ∂l_ij/∂u_k = tan(l_ij/2) (from l = 2·asin(exp(λ/2))): +// ∂α1/∂u1 = ∂α1/∂l12·t12 + ∂α1/∂l31·t31 +// ∂α1/∂u2 = ∂α1/∂l12·t12 + ∂α1/∂l23·t23 +// ∂α1/∂u3 = ∂α1/∂l23·t23 + ∂α1/∂l31·t31 +inline Eigen::SparseMatrix spherical_hessian( + ConformalMesh& mesh, + const std::vector& x, + const SphericalMaps& m) +{ + const int n = spherical_dimension(mesh, m); + + std::vector> trips; + trips.reserve(static_cast(n) * 9); + + for (auto f : mesh.faces()) { + Halfedge_index h0 = mesh.halfedge(f); + Halfedge_index h1 = mesh.next(h0); + Halfedge_index h2 = mesh.next(h1); + + Vertex_index v1 = mesh.source(h0); + Vertex_index v2 = mesh.source(h1); + Vertex_index v3 = mesh.source(h2); + + Edge_index e12 = mesh.edge(h0); + Edge_index e23 = mesh.edge(h1); + Edge_index e31 = mesh.edge(h2); + + // Effective log-lengths. + double u1 = spher_dof_val(m.v_idx[v1], x); + double u2 = spher_dof_val(m.v_idx[v2], x); + double u3 = spher_dof_val(m.v_idx[v3], x); + + double lam12 = m.lambda0[e12] + u1 + u2 + spher_dof_val(m.e_idx[e12], x); + double lam23 = m.lambda0[e23] + u2 + u3 + spher_dof_val(m.e_idx[e23], x); + double lam31 = m.lambda0[e31] + u3 + u1 + spher_dof_val(m.e_idx[e31], x); + + const double l12 = spherical_l(lam12); + const double l23 = spherical_l(lam23); + const double l31 = spherical_l(lam31); + + SphericalFaceAngles fa = spherical_angles(l12, l23, l31); + if (!fa.valid) continue; + + const double sinl12 = std::sin(l12), cosl12 = std::cos(l12); + const double sinl23 = std::sin(l23), cosl23 = std::cos(l23); + const double sinl31 = std::sin(l31), cosl31 = std::cos(l31); + + if (sinl12 < 1e-15 || sinl23 < 1e-15 || sinl31 < 1e-15) continue; + + const double cot12 = cosl12 / sinl12; + const double cot23 = cosl23 / sinl23; + const double cot31 = cosl31 / sinl31; + + // ∂l_ij/∂λ_ij = tan(l_ij/2) + const double t12 = std::tan(l12 * 0.5); + const double t23 = std::tan(l23 * 0.5); + const double t31 = std::tan(l31 * 0.5); + + const double sinA1 = std::sin(fa.alpha1), cosA1 = std::cos(fa.alpha1); + const double sinA2 = std::sin(fa.alpha2), cosA2 = std::cos(fa.alpha2); + const double sinA3 = std::sin(fa.alpha3), cosA3 = std::cos(fa.alpha3); + + if (sinA1 < 1e-15 || sinA2 < 1e-15 || sinA3 < 1e-15) continue; + + // ∂α1/∂l_jk (α1 at v1; opposite l23, adjacent l12,l31) + const double dA1_dl12 = (cot12 * cosA1 - cot31) / sinA1; + const double dA1_dl31 = (cot31 * cosA1 - cot12) / sinA1; + const double dA1_dl23 = sinl23 / (sinl12 * sinl31 * sinA1); + + // ∂α2/∂l_jk (α2 at v2; opposite l31, adjacent l12,l23) + const double dA2_dl12 = (cot12 * cosA2 - cot23) / sinA2; + const double dA2_dl23 = (cot23 * cosA2 - cot12) / sinA2; + const double dA2_dl31 = sinl31 / (sinl12 * sinl23 * sinA2); + + // ∂α3/∂l_jk (α3 at v3; opposite l12, adjacent l23,l31) + const double dA3_dl23 = (cot23 * cosA3 - cot31) / sinA3; + const double dA3_dl31 = (cot31 * cosA3 - cot23) / sinA3; + const double dA3_dl12 = sinl12 / (sinl23 * sinl31 * sinA3); + + // Chain rule: ∂α_i/∂u_j (u1 affects l12,l31; u2 affects l12,l23; u3 affects l23,l31) + const double dA1_du1 = dA1_dl12 * t12 + dA1_dl31 * t31; + const double dA1_du2 = dA1_dl12 * t12 + dA1_dl23 * t23; + const double dA1_du3 = dA1_dl23 * t23 + dA1_dl31 * t31; + + const double dA2_du1 = dA2_dl12 * t12 + dA2_dl31 * t31; + const double dA2_du2 = dA2_dl12 * t12 + dA2_dl23 * t23; + const double dA2_du3 = dA2_dl23 * t23 + dA2_dl31 * t31; + + const double dA3_du1 = dA3_dl12 * t12 + dA3_dl31 * t31; + const double dA3_du2 = dA3_dl12 * t12 + dA3_dl23 * t23; + const double dA3_du3 = dA3_dl23 * t23 + dA3_dl31 * t31; + + const int i1 = m.v_idx[v1]; + const int i2 = m.v_idx[v2]; + const int i3 = m.v_idx[v3]; + + // H[vi, vj] -= ∂α_i/∂u_j (G_v = θ_v − Σ α_v, so ∂G_i/∂u_j = −∂α_i/∂u_j) + if (i1 >= 0) trips.emplace_back(i1, i1, -dA1_du1); + if (i2 >= 0) trips.emplace_back(i2, i2, -dA2_du2); + if (i3 >= 0) trips.emplace_back(i3, i3, -dA3_du3); + + if (i1 >= 0 && i2 >= 0) { + trips.emplace_back(i1, i2, -dA1_du2); + trips.emplace_back(i2, i1, -dA2_du1); + } + if (i2 >= 0 && i3 >= 0) { + trips.emplace_back(i2, i3, -dA2_du3); + trips.emplace_back(i3, i2, -dA3_du2); + } + if (i3 >= 0 && i1 >= 0) { + trips.emplace_back(i3, i1, -dA3_du1); + trips.emplace_back(i1, i3, -dA1_du3); + } + } + + Eigen::SparseMatrix H(n, n); + H.setFromTriplets(trips.begin(), trips.end()); + return H; +} + +// ── Finite-difference Hessian check ────────────────────────────────────────── +// +// Compares the analytical Hessian column-by-column against +// H_fd[:, j] = (G(x + ε·eⱼ) − G(x − ε·eⱼ)) / (2ε). +inline bool hessian_check_spherical( + ConformalMesh& mesh, + const std::vector& x0, + const SphericalMaps& m, + double eps = 1e-5, + double tol = 1e-4) +{ + const int n = static_cast(x0.size()); + auto H = spherical_hessian(mesh, x0, m); + + std::vector xp = x0, xm = x0; + bool ok = true; + + for (int j = 0; j < n; ++j) { + const std::size_t sj = static_cast(j); + xp[sj] = x0[sj] + eps; + xm[sj] = x0[sj] - eps; + + auto Gp = spherical_gradient(mesh, xp, m); + auto Gm = spherical_gradient(mesh, xm, m); + + xp[sj] = xm[sj] = x0[sj]; + + for (int i = 0; i < n; ++i) { + double fd_ij = (Gp[static_cast(i)] + - Gm[static_cast(i)]) / (2.0 * eps); + double H_ij = H.coeff(i, j); + double err = std::abs(H_ij - fd_ij); + double scale = std::max(1.0, std::abs(H_ij)); + if (err / scale > tol) ok = false; + } + } + return ok; +} + +} // namespace conformallab diff --git a/code/tests/cgal/CMakeLists.txt b/code/tests/cgal/CMakeLists.txt index 2cc999b..a61a74c 100644 --- a/code/tests/cgal/CMakeLists.txt +++ b/code/tests/cgal/CMakeLists.txt @@ -20,6 +20,10 @@ add_executable(conformallab_cgal_tests # ── Phase 3d: EuclideanCyclicFunctional ─────────────────────────────── test_euclidean_functional.cpp + + # ── Phase 3f: Hessians (cotangent Laplacian, spherical + Euclidean) ─── + test_euclidean_hessian.cpp + test_spherical_hessian.cpp ) target_include_directories(conformallab_cgal_tests SYSTEM PRIVATE diff --git a/code/tests/cgal/test_euclidean_hessian.cpp b/code/tests/cgal/test_euclidean_hessian.cpp new file mode 100644 index 0000000..8ce14c1 --- /dev/null +++ b/code/tests/cgal/test_euclidean_hessian.cpp @@ -0,0 +1,213 @@ +// test_euclidean_hessian.cpp +// +// Phase 3f — Euclidean cotangent-Laplace Hessian. +// +// The Hessian of the Euclidean discrete conformal energy is the well-known +// cotangent-Laplace operator (Pinkall–Polthier 1993, Springborn 2008). +// +// Tests: +// 1. Cotangent weights are analytically correct for simple triangles. +// 2. Hessian is symmetric. +// 3. Hessian has the null-space property H·1 = 0 (uniform-shift mode). +// 4. Hessian is positive semi-definite (all eigenvalues ≥ 0). +// 5. Finite-difference check H[i,j] ≈ (G_i(x+ε·eⱼ)−G_i(x−ε·eⱼ))/(2ε). +// +// All tests use meshes and maps built with Phase-3d infrastructure. + +#include "conformal_mesh.hpp" +#include "mesh_builder.hpp" +#include "euclidean_hessian.hpp" +#include +#include // for dense conversion and eigenvalue solver +#include +#include + +using namespace conformallab; + +// ════════════════════════════════════════════════════════════════════════════ +// Cotangent weight: equilateral triangle → all cots = 1/√3 = cot(60°) +// ════════════════════════════════════════════════════════════════════════════ + +TEST(EuclideanHessian, CotWeights_EquilateralTriangle) +{ + // Equilateral triangle with l = 1 (all log-lengths = 0). + auto cw = euclidean_cot_weights(1.0, 1.0, 1.0); + + ASSERT_TRUE(cw.valid); + const double expected = 1.0 / std::sqrt(3.0); // cot(60°) + EXPECT_NEAR(cw.cot1, expected, 1e-12); + EXPECT_NEAR(cw.cot2, expected, 1e-12); + EXPECT_NEAR(cw.cot3, expected, 1e-12); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Cotangent weight: right-isosceles triangle (legs 1, hypotenuse √2) +// +// v1=(0,0): right angle → cot(90°) = 0 +// v2=(1,0), v3=(0,1): 45° angles → cot(45°) = 1 +// ════════════════════════════════════════════════════════════════════════════ + +TEST(EuclideanHessian, CotWeights_RightIsoscelesTriangle) +{ + // l12=1, l23=√2, l31=1 + auto cw = euclidean_cot_weights(1.0, std::sqrt(2.0), 1.0); + + ASSERT_TRUE(cw.valid); + EXPECT_NEAR(cw.cot1, 0.0, 1e-12); // right angle at v1 + EXPECT_NEAR(cw.cot2, 1.0, 1e-12); // 45° at v2 + EXPECT_NEAR(cw.cot3, 1.0, 1e-12); // 45° at v3 +} + +// ════════════════════════════════════════════════════════════════════════════ +// Hessian is symmetric: H[i,j] == H[j,i] +// ════════════════════════════════════════════════════════════════════════════ + +TEST(EuclideanHessian, HessianIsSymmetric) +{ + auto mesh = make_quad_strip(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + int n = assign_euclidean_vertex_dof_indices(mesh, maps); + + std::vector x(static_cast(n), -0.1); + auto H = euclidean_hessian(mesh, x, maps); + + Eigen::MatrixXd Hd = Eigen::MatrixXd(H); + EXPECT_NEAR((Hd - Hd.transpose()).norm(), 0.0, 1e-12) + << "Hessian must be symmetric"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Null-space property: H·1 = 0 for a closed surface (regular tetrahedron) +// +// The cotangent Laplacian on a closed mesh has the constant vector in its +// null space (each row sums to zero). +// ════════════════════════════════════════════════════════════════════════════ + +TEST(EuclideanHessian, NullSpaceIsConstantVector_ClosedMesh) +{ + auto mesh = make_tetrahedron(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + int n = assign_euclidean_vertex_dof_indices(mesh, maps); + + std::vector x(static_cast(n), 0.0); + auto H = euclidean_hessian(mesh, x, maps); + + // 1-vector + Eigen::VectorXd ones = Eigen::VectorXd::Ones(n); + Eigen::VectorXd Hones = H * ones; + + EXPECT_NEAR(Hones.norm(), 0.0, 1e-10) + << "H·1 must be zero on a closed mesh (cotangent Laplacian null-space)"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Hessian is positive semi-definite: all eigenvalues ≥ 0 +// +// Checked on a small mesh (regular tetrahedron, 4 vertices) using dense +// self-adjoint eigenvalue decomposition (only feasible for small n). +// ════════════════════════════════════════════════════════════════════════════ + +TEST(EuclideanHessian, HessianIsPositiveSemiDefinite) +{ + auto mesh = make_tetrahedron(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + int n = assign_euclidean_vertex_dof_indices(mesh, maps); + + std::vector x(static_cast(n), 0.0); + auto H = euclidean_hessian(mesh, x, maps); + + Eigen::MatrixXd Hd = Eigen::MatrixXd(H); + Eigen::SelfAdjointEigenSolver es(Hd); + double min_ev = es.eigenvalues().minCoeff(); + + EXPECT_GE(min_ev, -1e-10) + << "All eigenvalues of the cotangent Laplacian must be ≥ 0; " + "smallest = " << min_ev; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Finite-difference Hessian check: single right-isosceles triangle +// +// H[i,j] ≈ (G_i(x+ε·eⱼ) − G_i(x−ε·eⱼ)) / (2ε) +// ════════════════════════════════════════════════════════════════════════════ + +TEST(EuclideanHessian, FDCheck_Triangle) +{ + auto mesh = make_triangle(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + int n = assign_euclidean_vertex_dof_indices(mesh, maps); + + std::vector x(static_cast(n), -0.1); + + EXPECT_TRUE(hessian_check_euclidean(mesh, x, maps)) + << "FD Hessian check failed on right-isosceles triangle"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Finite-difference Hessian check: quad strip (2 triangles, 1 interior edge) +// ════════════════════════════════════════════════════════════════════════════ + +TEST(EuclideanHessian, FDCheck_QuadStrip) +{ + auto mesh = make_quad_strip(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + int n = assign_euclidean_vertex_dof_indices(mesh, maps); + + std::vector x(static_cast(n), -0.1); + + EXPECT_TRUE(hessian_check_euclidean(mesh, x, maps)) + << "FD Hessian check failed on quad strip"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Finite-difference Hessian check: regular tetrahedron (closed, 4 faces) +// ════════════════════════════════════════════════════════════════════════════ + +TEST(EuclideanHessian, FDCheck_Tetrahedron) +{ + auto mesh = make_tetrahedron(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + int n = assign_euclidean_vertex_dof_indices(mesh, maps); + + std::vector x(static_cast(n), -0.15); + + EXPECT_TRUE(hessian_check_euclidean(mesh, x, maps)) + << "FD Hessian check failed on regular tetrahedron"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Finite-difference Hessian check: with mixed pinned/variable vertices +// +// One vertex pinned: the corresponding row/column must be absent from H +// while the diagonal of neighbouring variable vertices still gets the full +// cotangent contribution. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(EuclideanHessian, FDCheck_MixedPinnedVertices) +{ + auto mesh = make_quad_strip(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + + 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 + maps.v_idx[v1] = 0; + maps.v_idx[v2] = 1; + maps.v_idx[v3] = 2; + + std::vector x = {-0.1, -0.2, -0.15}; + + EXPECT_TRUE(hessian_check_euclidean(mesh, x, maps)) + << "FD Hessian check failed for mixed pinned/variable vertices"; +} diff --git a/code/tests/cgal/test_spherical_hessian.cpp b/code/tests/cgal/test_spherical_hessian.cpp new file mode 100644 index 0000000..ff3727d --- /dev/null +++ b/code/tests/cgal/test_spherical_hessian.cpp @@ -0,0 +1,205 @@ +// test_spherical_hessian.cpp +// +// Phase 3f — Spherical cotangent-Laplace Hessian. +// +// The Hessian of the spherical discrete conformal energy is the "spherical +// cotangent Laplacian" with edge weight +// w_k = cot(β_k), β_k = (π − αi − αj + αk) / 2 +// for the edge (vi, vj) with opposite vertex vk and angles αi, αj, αk. +// +// In the flat limit (αi+αj+αk → π) this reduces to the Euclidean cotangent +// Laplacian (β_k → αk, w_k → cot(αk)). +// +// Tests: +// 1. Spherical cot weights match Euclidean weights in the near-flat limit. +// 2. Hessian is symmetric. +// 3. Hessian has H·1 ≈ 0 on the spherical tetrahedron (null-space property). +// 4. Finite-difference check (the primary correctness criterion). + +#include "conformal_mesh.hpp" +#include "mesh_builder.hpp" +#include "spherical_hessian.hpp" +#include "euclidean_hessian.hpp" // for Euclidean comparison +#include +#include +#include +#include + +using namespace conformallab; + +// ════════════════════════════════════════════════════════════════════════════ +// Spherical cot weights reduce to Euclidean cot weights in the flat limit +// +// For a nearly-flat equilateral spherical triangle (l → 0, α → 60°): +// β_k = (π − 60° − 60° + 60°)/2 = 60° → cot(60°) = 1/√3 ✓ +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalHessian, CotWeights_NearFlatEquilateral) +{ + // Use a very small equilateral spherical triangle: α1=α2=α3=60° + const double alpha = PI / 3.0; + auto sw = spherical_cot_weights(alpha, alpha, alpha); + + ASSERT_TRUE(sw.valid); + + const double expected = 1.0 / std::sqrt(3.0); // = cot(60°) + EXPECT_NEAR(sw.w12, expected, 1e-12); + EXPECT_NEAR(sw.w23, expected, 1e-12); + EXPECT_NEAR(sw.w31, expected, 1e-12); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Spherical cot weights are positive for the regular spherical tetrahedron +// +// Each face has α_k = 2π/3 (120°), angle sum = 2π. +// β_k = (π − 2π/3 − 2π/3 + 2π/3)/2 = (π − 2π/3)/2 = π/6 +// w_k = cot(π/6) = √3 +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalHessian, CotWeights_RegularSphericalTetrahedronFace) +{ + const double alpha = 2.0 * PI / 3.0; // 120° + auto sw = spherical_cot_weights(alpha, alpha, alpha); + + ASSERT_TRUE(sw.valid); + + const double expected = std::sqrt(3.0); // cot(π/6) + EXPECT_NEAR(sw.w12, expected, 1e-10); + EXPECT_NEAR(sw.w23, expected, 1e-10); + EXPECT_NEAR(sw.w31, expected, 1e-10); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Hessian is symmetric: H[i,j] == H[j,i] +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalHessian, HessianIsSymmetric) +{ + 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 x(static_cast(n), -0.2); + auto H = spherical_hessian(mesh, x, maps); + + Eigen::MatrixXd Hd = Eigen::MatrixXd(H); + EXPECT_NEAR((Hd - Hd.transpose()).norm(), 0.0, 1e-12) + << "Spherical Hessian must be symmetric"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// H·1 is NOT zero for the spherical Hessian (unlike the Euclidean case). +// +// In the Euclidean case, a uniform shift u_i → u_i + c scales all edge +// lengths by e^c, leaving angles unchanged → H·1 = 0 exactly. +// +// In the spherical case, l_ij = 2·asin(exp(λ_ij/2)) is NOT a linear +// function of the DOFs, so a uniform shift DOES change the angles +// → H·1 ≠ 0 in general. This test verifies that property. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalHessian, ConstantVectorNotInNullSpace) +{ + 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 x(static_cast(n), -0.2); + auto H = spherical_hessian(mesh, x, maps); + + Eigen::VectorXd ones = Eigen::VectorXd::Ones(n); + Eigen::VectorXd Hones = H * ones; + + // H·1 should be non-trivial (norm well above zero) + EXPECT_GT(Hones.norm(), 1e-3) + << "Spherical H·1 should be non-zero; norm = " << Hones.norm(); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Hessian is negative semi-definite at the equilibrium point (x = 0) +// +// The spherical discrete conformal energy is concave in the vertex DOFs +// (unlike the Euclidean case which is convex). At the equilibrium x = 0, +// the Hessian is NSD: all eigenvalues ≤ 0, with a multi-dimensional null +// space corresponding to degenerate directions. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalHessian, HessianIsNegativeSemiDefiniteAtEquilibrium) +{ + 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); + + // x = 0 is the equilibrium for the regular spherical tetrahedron. + std::vector x(static_cast(n), 0.0); + auto H = spherical_hessian(mesh, x, maps); + + Eigen::MatrixXd Hd = Eigen::MatrixXd(H); + Eigen::SelfAdjointEigenSolver es(Hd); + double max_ev = es.eigenvalues().maxCoeff(); + + EXPECT_LE(max_ev, 1e-9) + << "Largest eigenvalue of spherical Hessian at equilibrium must be ≤ 0; got " << max_ev; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Finite-difference Hessian check: octahedron-face triangle (vertex DOFs) +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalHessian, FDCheck_OctaFaceVertex) +{ + auto mesh = make_octahedron_face(); + auto maps = setup_spherical_maps(mesh); + compute_lambda0_from_mesh(mesh, maps); + int n = assign_vertex_dof_indices(mesh, maps); + + std::vector x(static_cast(n), -0.3); + + EXPECT_TRUE(hessian_check_spherical(mesh, x, maps)) + << "FD Hessian check failed on octahedron-face triangle"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Finite-difference Hessian check: spherical tetrahedron (vertex DOFs) +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalHessian, FDCheck_SpherTetVertex) +{ + 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 x(static_cast(n), -0.2); + + EXPECT_TRUE(hessian_check_spherical(mesh, x, maps)) + << "FD Hessian check failed on spherical tetrahedron"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Finite-difference Hessian check: mixed pinned/variable vertices +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalHessian, FDCheck_MixedPinnedVertices) +{ + auto mesh = make_octahedron_face(); + auto maps = setup_spherical_maps(mesh); + compute_lambda0_from_mesh(mesh, maps); + + auto vit = mesh.vertices().begin(); + Vertex_index v0 = *vit++; + Vertex_index v1 = *vit++; + Vertex_index v2 = *vit; + + maps.v_idx[v0] = -1; // pinned + maps.v_idx[v1] = 0; + maps.v_idx[v2] = 1; + + std::vector x = {-0.2, -0.3}; + + EXPECT_TRUE(hessian_check_spherical(mesh, x, maps)) + << "FD Hessian check failed for mixed pinned/variable vertices"; +}