feat(phase3f+3g): analytical Hessians + PI consolidation

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 <noreply@anthropic.com>

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

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 <CGAL/boost/graph/iterator.h>
 #include <vector>
@@ -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<Vertex_index, int>   ("ev:idx",   -1     ).first;
     m.e_idx  = mesh.add_property_map<Edge_index,   int>   ("ee:idx",   -1     ).first;

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<double> 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 <Eigen/Sparse>
+#include <vector>
+#include <cmath>
+
+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<double> euclidean_hessian(
+    ConformalMesh&             mesh,
+    const std::vector<double>& x,
+    const EuclideanMaps&       m)
+{
+    const int n = euclidean_dimension(mesh, m);
+
+    // Collect triplets (row, col, value) — setFromTriplets sums duplicates.
+    std::vector<Eigen::Triplet<double>> trips;
+    trips.reserve(static_cast<std::size_t>(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<double> 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<double>& x0,
+    const EuclideanMaps&       m,
+    double                     eps = 1e-5,
+    double                     tol = 1e-4)
+{
+    const int n = static_cast<int>(x0.size());
+    auto H = euclidean_hessian(mesh, x0, m);
+
+    std::vector<double> xp = x0, xm = x0;
+    bool ok = true;
+
+    for (int j = 0; j < n; ++j) {
+        const std::size_t sj = static_cast<std::size_t>(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<std::size_t>(i)]
+                           - Gm[static_cast<std::size_t>(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

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 <cmath>
 #include <algorithm>
 
 namespace conformallab {
 
-constexpr double PI = 3.14159265358979323846264338328;
-
 // ── Length functions ─────────────────────────────────────────────────────────
 
 // ζ(x,y,z)  — interior angle in a hyperbolic triangle with edge lengths

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 <Eigen/Dense>
 #include <cmath>
@@ -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 + alpha31 - alpha12 - gamma1) / 2.0);
-    result += L((pi + alpha12 - alpha23 - gamma2) / 2.0);
+    result += L((PI + alpha12 - alpha23 - gamma2) / 2.0);
-    result += L((pi + alpha23 - alpha31 - gamma3) / 2.0);
+    result += L((PI + alpha23 - alpha31 - gamma3) / 2.0);
 
-    result += L((pi - alpha31 + alpha12 - gamma1) / 2.0);
+    result += L((PI - alpha31 + alpha12 - gamma1) / 2.0);
-    result += L((pi - alpha12 + alpha23 - gamma2) / 2.0);
+    result += L((PI - alpha12 + alpha23 - gamma2) / 2.0);
-    result += L((pi - alpha23 + alpha31 - gamma3) / 2.0);
+    result += L((PI - alpha23 + alpha31 - gamma3) / 2.0);
 
-    result += L((pi + alpha31 + alpha12 - gamma1) / 2.0);
+    result += L((PI + alpha31 + alpha12 - gamma1) / 2.0);
-    result += L((pi + alpha12 + alpha23 - gamma2) / 2.0);
+    result += L((PI + alpha12 + alpha23 - gamma2) / 2.0);
-    result += L((pi + alpha23 + alpha31 - gamma3) / 2.0);
+    result += L((PI + alpha23 + alpha31 - gamma3) / 2.0);
 
-    result += L((pi - alpha31 - alpha12 - gamma1) / 2.0);
+    result += L((PI - alpha31 - alpha12 - gamma1) / 2.0);
-    result += L((pi - alpha12 - alpha23 - gamma2) / 2.0);
+    result += L((PI - alpha12 - alpha23 - gamma2) / 2.0);
-    result += L((pi - alpha23 - alpha31 - gamma3) / 2.0);
+    result += L((PI - alpha23 - alpha31 - gamma3) / 2.0);
 
     return result / 2.0;
 }

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 <cmath>
 #include <vector>
 
@@ -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<Vertex_index> rim(n);
     for (int i = 0; i < n; ++i) {
         double a = i * dtheta;

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 <cmath>
 #include <algorithm>
 
 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 ────────────────────────────────────────────
 

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<double>.
+
+#include "spherical_functional.hpp"
+#include <Eigen/Sparse>
+#include <vector>
+#include <cmath>
+
+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<double> spherical_hessian(
+    ConformalMesh&             mesh,
+    const std::vector<double>& x,
+    const SphericalMaps&       m)
+{
+    const int n = spherical_dimension(mesh, m);
+
+    std::vector<Eigen::Triplet<double>> trips;
+    trips.reserve(static_cast<std::size_t>(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<double> 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<double>& x0,
+    const SphericalMaps&       m,
+    double                     eps = 1e-5,
+    double                     tol = 1e-4)
+{
+    const int n = static_cast<int>(x0.size());
+    auto H = spherical_hessian(mesh, x0, m);
+
+    std::vector<double> xp = x0, xm = x0;
+    bool ok = true;
+
+    for (int j = 0; j < n; ++j) {
+        const std::size_t sj = static_cast<std::size_t>(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<std::size_t>(i)]
+                          - Gm[static_cast<std::size_t>(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

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

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 <gtest/gtest.h>
+#include <Eigen/Dense>   // for dense conversion and eigenvalue solver
+#include <cmath>
+#include <vector>
+
+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<double> x(static_cast<std::size_t>(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<double> x(static_cast<std::size_t>(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<double> x(static_cast<std::size_t>(n), 0.0);
+    auto H = euclidean_hessian(mesh, x, maps);
+
+    Eigen::MatrixXd Hd = Eigen::MatrixXd(H);
+    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> 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<double> x(static_cast<std::size_t>(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<double> x(static_cast<std::size_t>(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<double> x(static_cast<std::size_t>(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<double> x = {-0.1, -0.2, -0.15};
+
+    EXPECT_TRUE(hessian_check_euclidean(mesh, x, maps))
+        << "FD Hessian check failed for mixed pinned/variable vertices";
+}

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 <gtest/gtest.h>
+#include <Eigen/Dense>
+#include <cmath>
+#include <vector>
+
+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<double> x(static_cast<std::size_t>(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<double> x(static_cast<std::size_t>(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<double> x(static_cast<std::size_t>(n), 0.0);
+    auto H = spherical_hessian(mesh, x, maps);
+
+    Eigen::MatrixXd Hd = Eigen::MatrixXd(H);
+    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> 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<double> x(static_cast<std::size_t>(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<double> x(static_cast<std::size_t>(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<double> x = {-0.2, -0.3};
+
+    EXPECT_TRUE(hessian_check_spherical(mesh, x, maps))
+        << "FD Hessian check failed for mixed pinned/variable vertices";
+}