Merge pull request 'feat(euclidean): edge-DOF (cyclic) Hessian — block-FD + full analytic + Java cross-validation' (#29) from docs/java-ignore-crossvalidation into main

code/include/euclidean_hessian.hpp

@@ -43,6 +43,7 @@
 #include "euclidean_functional.hpp"
 #include <Eigen/Sparse>
 #include <vector>
+#include <array>
 #include <cmath>
 #include <stdexcept>
 
@@ -187,6 +188,232 @@ inline Eigen::SparseMatrix<double> euclidean_hessian(
     return H;
 }
 
+// ── Block-FD Hessian (cyclic: vertex + edge DOFs) ────────────────────────────
+//
+// The analytic `euclidean_hessian` above covers only the vertex-block cotangent
+// Laplacian.  The *cyclic* functional adds edge DOFs λ_e, giving vertex-edge and
+// edge-edge Hessian blocks.  We obtain the full Hessian by per-face block FD,
+// mirroring `hyper_ideal_hessian_block_fd`.
+//
+// Locality lemma: the gradient decomposes by face,
+//     G_v = Θ_v − Σ_{f∋v} α_v(f)            (vertex output  = −α_v)
+//     G_e = Σ_{f∋e} α_opp(f) − φ_e          (edge output    = +α_opp)
+// and α at face f depends ONLY on f's 6 local DOFs (u₁,u₂,u₃,λ₁₂,λ₂₃,λ₃₁), so
+//     ∂G_x/∂y = Σ_{f: x,y ∈ local(f)} ∂(output)/∂y at f.
+// Accumulating per-face 6×6 blocks therefore reproduces the full Hessian and is
+// consistent with `euclidean_gradient` by construction (it is its FD).
+
+/// Per-face contributions to the cyclic gradient, in slot order
+/// (v₁, v₂, v₃, e₁₂, e₂₃, e₃₁).  Vertex slots carry −α (since G_v = Θ−Σα),
+/// edge slots carry +α_opp (since G_e = Σα_opp − φ).
+inline std::array<double, 6> eucl_face_local_outputs(
+    double u1, double u2, double u3,
+    double d12, double d23, double d31,
+    double lam0_12, double lam0_23, double lam0_31)
+{
+    const double lam12 = lam0_12 + u1 + u2 + d12;
+    const double lam23 = lam0_23 + u2 + u3 + d23;
+    const double lam31 = lam0_31 + u3 + u1 + d31;
+    const auto fa = euclidean_angles(lam12, lam23, lam31);
+    // edge slot e12 ↔ opposite corner α3, e23 ↔ α1, e31 ↔ α2 (see gradient Pass 3).
+    return { -fa.alpha1, -fa.alpha2, -fa.alpha3,
+             +fa.alpha3, +fa.alpha1, +fa.alpha2 };
+}
+
+/// Block-FD Euclidean cyclic Hessian (vertex + edge DOFs), sparse.
+/// Supports both vertex-only and cyclic DOF layouts; cost `F·12` face-angle
+/// evaluations.  Consistent with `euclidean_gradient` to `O(ε²)`.
+inline Eigen::SparseMatrix<double> euclidean_hessian_block_fd(
+    ConformalMesh&             mesh,
+    const std::vector<double>& x,
+    const EuclideanMaps&       m,
+    double                     eps = 1e-5)
+{
+    const int n = euclidean_dimension(mesh, m);
+    std::vector<Eigen::Triplet<double>> trips;
+    trips.reserve(static_cast<std::size_t>(36) * mesh.number_of_faces());
+
+    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);
+
+        const int idx[6] = {
+            m.v_idx[v1], m.v_idx[v2], m.v_idx[v3],
+            m.e_idx[e12], m.e_idx[e23], m.e_idx[e31]
+        };
+        const double lam0_12 = m.lambda0[e12];
+        const double lam0_23 = m.lambda0[e23];
+        const double lam0_31 = m.lambda0[e31];
+
+        const double vals[6] = {
+            eucl_dof_val(idx[0], x), eucl_dof_val(idx[1], x), eucl_dof_val(idx[2], x),
+            eucl_dof_val(idx[3], x), eucl_dof_val(idx[4], x), eucl_dof_val(idx[5], x)
+        };
+
+        for (int j = 0; j < 6; ++j) {
+            if (idx[j] < 0) continue;   // pinned: no column
+
+            double vp[6], vm[6];
+            for (int k = 0; k < 6; ++k) { vp[k] = vm[k] = vals[k]; }
+            vp[j] += eps;
+            vm[j] -= eps;
+
+            auto Op = eucl_face_local_outputs(
+                vp[0], vp[1], vp[2], vp[3], vp[4], vp[5], lam0_12, lam0_23, lam0_31);
+            auto Om = eucl_face_local_outputs(
+                vm[0], vm[1], vm[2], vm[3], vm[4], vm[5], lam0_12, lam0_23, lam0_31);
+
+            for (int i = 0; i < 6; ++i) {
+                if (idx[i] < 0) continue;
+                const double val = (Op[static_cast<std::size_t>(i)]
+                                  - Om[static_cast<std::size_t>(i)]) / (2.0 * eps);
+                if (std::abs(val) > 1e-15)
+                    trips.emplace_back(idx[i], idx[j], val);
+            }
+        }
+    }
+
+    Eigen::SparseMatrix<double> H(n, n);
+    H.setFromTriplets(trips.begin(), trips.end());
+    return H;
+}
+
+/// Symmetrised block-FD Euclidean cyclic Hessian: `(H + Hᵀ)/2`.
+inline Eigen::SparseMatrix<double> euclidean_hessian_block_fd_sym(
+    ConformalMesh&             mesh,
+    const std::vector<double>& x,
+    const EuclideanMaps&       m,
+    double                     eps = 1e-5)
+{
+    auto H = euclidean_hessian_block_fd(mesh, x, m, eps);
+    Eigen::SparseMatrix<double> Ht = H.transpose();
+    return (H + Ht) * 0.5;
+}
+
+// ── Analytic cyclic Hessian (vertex + edge DOFs) ─────────────────────────────
+//
+// Closed-form Jacobian of the cyclic gradient.  With s_i = effective log-length
+// (Λ̃) of the edge OPPOSITE vertex i (s = 2·ln ℓ), the Euclidean corner-angle
+// derivatives are (derived from the law of cosines; Σ_j ∂α_i/∂s_j = 0):
+//
+//     ∂α_i/∂s_i =  ℓ_i² / (4A)
+//     ∂α_i/∂s_j =  ½ cot α_i  −  ℓ_j² / (4A)      (j ≠ i)
+//
+// where A is the triangle area (4A = denom2/2, denom2 = 8A from the cot helper).
+// Chaining through  s₁ = Λ̃(e₂₃), s₂ = Λ̃(e₃₁), s₃ = Λ̃(e₁₂)  and
+// Λ̃_e = λ⁰_e + u_i + u_j + λ_e  gives the per-face 6×6 block over the local
+// DOFs (u₁,u₂,u₃,λ₁₂,λ₂₃,λ₃₁); per-face outputs carry the gradient signs
+// (−α vertex, +α_opp edge), so the assembled matrix equals ∂G/∂x exactly.
+// This is the analytic counterpart of `euclidean_hessian_block_fd` (no FD).
+
+/// Analytic Euclidean cyclic Hessian (vertex + edge DOFs), sparse.
+/// Closed-form (no finite differences); matches `euclidean_hessian_block_fd`
+/// to round-off and the analytic vertex-only Laplacian on vertex-only layouts.
+inline Eigen::SparseMatrix<double> euclidean_hessian_analytic(
+    ConformalMesh&             mesh,
+    const std::vector<double>& x,
+    const EuclideanMaps&       m)
+{
+    const int n = euclidean_dimension(mesh, m);
+    std::vector<Eigen::Triplet<double>> trips;
+    trips.reserve(static_cast<std::size_t>(36) * mesh.number_of_faces());
+
+    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);
+
+        const double u1 = eucl_dof_val(m.v_idx[v1], x);
+        const double u2 = eucl_dof_val(m.v_idx[v2], x);
+        const double u3 = eucl_dof_val(m.v_idx[v3], x);
+
+        const double lam12 = m.lambda0[e12] + u1 + u2 + eucl_dof_val(m.e_idx[e12], x); // s3
+        const double lam23 = m.lambda0[e23] + u2 + u3 + eucl_dof_val(m.e_idx[e23], x); // s1
+        const double lam31 = m.lambda0[e31] + u3 + u1 + eucl_dof_val(m.e_idx[e31], x); // s2
+
+        // Centered side lengths (scale-invariant for ℓ²/4A and cot).
+        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);
+
+        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) continue;
+        const double l123 = l12 + l23 + l31;
+        const double denom2_sq = t12 * t23 * t31 * l123;
+        if (denom2_sq <= 0.0) continue;
+        const double denom2 = 2.0 * std::sqrt(denom2_sq);   // = 8A
+        const double fourA  = denom2 * 0.5;                 // = 4A
+
+        auto cw = euclidean_cot_weights(l12, l23, l31);
+        if (!cw.valid) continue;
+        const double cot[3] = { cw.cot1, cw.cot2, cw.cot3 };  // at v1, v2, v3
+
+        // ℓ_i² = squared length of edge opposite vertex i:
+        //   opp(v1)=e23=l23,  opp(v2)=e31=l31,  opp(v3)=e12=l12
+        const double Lsq[3] = { l23 * l23, l31 * l31, l12 * l12 };
+
+        // dA[i][j] = ∂α_i/∂s_j   (i,j ∈ {0,1,2}; s_j opposite v_{j+1})
+        double dA[3][3];
+        for (int i = 0; i < 3; ++i)
+            for (int j = 0; j < 3; ++j)
+                dA[i][j] = (i == j) ? (Lsq[i] / fourA)
+                                    : (0.5 * cot[i] - Lsq[j] / fourA);
+
+        // ∂α_i/∂(local dof), dof order (u1,u2,u3, λ12,λ23,λ31).
+        // u1 ↔ s2,s3 ; u2 ↔ s1,s3 ; u3 ↔ s1,s2 ; λ12 ↔ s3 ; λ23 ↔ s1 ; λ31 ↔ s2.
+        double g[3][6];
+        for (int i = 0; i < 3; ++i) {
+            g[i][0] = dA[i][1] + dA[i][2];  // u1
+            g[i][1] = dA[i][0] + dA[i][2];  // u2
+            g[i][2] = dA[i][0] + dA[i][1];  // u3
+            g[i][3] = dA[i][2];             // λ12 (s3)
+            g[i][4] = dA[i][0];             // λ23 (s1)
+            g[i][5] = dA[i][1];             // λ31 (s2)
+        }
+
+        // Output slot → (angle index, sign): (−α1,−α2,−α3,+α3,+α1,+α2).
+        const int    out_a[6]  = { 0, 1, 2, 2, 0, 1 };
+        const double out_sg[6] = { -1.0, -1.0, -1.0, +1.0, +1.0, +1.0 };
+
+        const int idx[6] = {
+            m.v_idx[v1], m.v_idx[v2], m.v_idx[v3],
+            m.e_idx[e12], m.e_idx[e23], m.e_idx[e31]
+        };
+
+        for (int oi = 0; oi < 6; ++oi) {
+            if (idx[oi] < 0) continue;
+            for (int dj = 0; dj < 6; ++dj) {
+                if (idx[dj] < 0) continue;
+                const double val = out_sg[oi] * g[out_a[oi]][dj];
+                if (std::abs(val) > 1e-15)
+                    trips.emplace_back(idx[oi], idx[dj], val);
+            }
+        }
+    }
+
+    Eigen::SparseMatrix<double> H(n, n);
+    H.setFromTriplets(trips.begin(), trips.end());
+    return H;
+}
+
 // ── Finite-difference Hessian check ──────────────────────────────────────────
 /// FD Hessian check for the Euclidean functional.  Compares analytic
 /// `H` column-by-column to `(G(x+εeⱼ) − G(x−εeⱼ)) / (2ε)`; returns

code/include/newton_solver.hpp

@@ -258,7 +258,14 @@ inline NewtonResult newton_euclidean(
         }
 
         // ── Hessian + solve H·Δx = −G  (SparseQR fallback for singular H) ──
-        auto H = euclidean_hessian(mesh, x, m);
+        // Cyclic layout (edge DOFs present) → analytic cyclic Hessian, covering
+        // the vertex-edge / edge-edge blocks the vertex-only cotangent Laplacian
+        // omits.  Vertex-only layout → analytic cotangent Laplacian.
+        bool has_edge_dof = false;
+        for (auto e : mesh.edges())
+            if (m.e_idx[e] >= 0) { has_edge_dof = true; break; }
+        auto H = has_edge_dof ? euclidean_hessian_analytic(mesh, x, m)
+                              : euclidean_hessian(mesh, x, m);
         bool ok = false;
         Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok);
         if (!ok) break;

code/tests/cgal/test_euclidean_functional.cpp

@@ -27,9 +27,12 @@
 #include "euclidean_functional.hpp"
 #include "euclidean_hessian.hpp"
 #include "clausen.hpp"
+#include "mesh_io.hpp"
+#include "newton_solver.hpp"
 #include <gtest/gtest.h>
 #include <cmath>
 #include <vector>
+#include <string>
 
 using namespace conformallab;
 
@@ -405,3 +408,240 @@ TEST(EuclideanGoldenJava, FullMeshGradientAndEnergy_Tetrahedron)
                     - euclidean_energy(mesh, x0, maps);
     EXPECT_NEAR(dE, 0.15962619236187336, 1e-12);
 }
+
+// ════════════════════════════════════════════════════════════════════════════
+// Java cross-validation (Tier 1, GREEN) — circular-edge φ wiring (no solver)
+//
+// The "circular hole edge" of EuclideanCyclicConvergenceTest works by setting a
+// non-default edge turn angle φ_e. Since the cyclic edge gradient is exactly
+//     G_e = α_opp(f⁺) + α_opp(f⁻) − φ_e,
+// lowering φ_e by 0.1 must raise G_e by exactly 0.1 — independent of geometry —
+// and must leave every other gradient component untouched. This pins the φ
+// wiring without needing the (not-yet-implemented) edge-DOF Hessian, so it runs
+// today and is the evaluation-level prerequisite of the DISABLED convergence
+// test below.
+// ════════════════════════════════════════════════════════════════════════════
+TEST(EuclideanFunctional, CyclicCircularEdge_PhiEntersGradient_CatHead)
+{
+    const std::string path = std::string(CONFORMALLAB_DATA_DIR) + "/obj/cathead.obj";
+    ConformalMesh mesh;
+    ASSERT_NO_THROW(mesh = load_mesh(path)) << "cathead.obj not found: " << path;
+
+    auto maps = setup_euclidean_maps(mesh);
+    compute_euclidean_lambda0_from_mesh(mesh, maps);
+
+    int idx = 0;
+    for (auto v : mesh.vertices())
+        maps.v_idx[v] = mesh.is_border(v) ? -1 : idx++;
+    for (auto e : mesh.edges())
+        maps.e_idx[e] = idx++;
+    const int n = idx;
+    ASSERT_GT(n, 0);
+
+    Edge_index e_circ{};
+    bool found = false;
+    for (auto e : mesh.edges()) {
+        auto h  = mesh.halfedge(e);
+        auto ho = mesh.opposite(h);
+        if (mesh.is_border(h) || mesh.is_border(ho)) continue;
+        e_circ = e; found = true; break;
+    }
+    ASSERT_TRUE(found) << "no interior edge found on cathead";
+    const std::size_t ie = static_cast<std::size_t>(maps.e_idx[e_circ]);
+
+    std::vector<double> x(static_cast<std::size_t>(n), 0.0);
+
+    maps.phi_e[e_circ] = PI;
+    auto G1 = euclidean_gradient(mesh, x, maps);
+    maps.phi_e[e_circ] = PI - 0.1;
+    auto G2 = euclidean_gradient(mesh, x, maps);
+
+    // Lowering φ_e by 0.1 raises exactly this edge's gradient component by 0.1.
+    EXPECT_NEAR(0.1, G2[ie] - G1[ie], 1e-12)
+        << "circular-edge φ target not wired into the cyclic gradient";
+
+    // No other gradient component changes.
+    double max_other = 0.0;
+    for (std::size_t k = 0; k < G1.size(); ++k)
+        if (k != ie) max_other = std::max(max_other, std::abs(G2[k] - G1[k]));
+    EXPECT_LT(max_other, 1e-12)
+        << "changing one φ_e perturbed unrelated gradient components";
+}
+
+// ════════════════════════════════════════════════════════════════════════════
+// Java cross-validation (Tier 1) — EuclideanCyclicConvergenceTest
+//
+// Ports de.varylab.discreteconformal.functional.EuclideanCyclicConvergenceTest:
+// prescribe a non-default edge turn angle φ = π − 0.1 on one interior edge of
+// cathead.obj ("circular hole edge"), solve the cyclic Euclidean functional
+// (vertex + edge DOFs), then assert the realised opposite-corner-angle sum
+// across that edge equals π − 0.1.
+//
+// The C++ edge gradient is  G_e = α_opp(f⁺) + α_opp(f⁻) − φ_e, so at the
+// solution (G_e = 0) the geometric angle sum equals φ_e — exactly the Java
+// assertion `circularEdge.getAlpha() + opposite.getAlpha() == π − 0.1`.
+//
+// "Natural targets" first make x = 0 the equilibrium (so the *only* deviation
+// is the prescribed φ); the test therefore FAILS if the solver ignores a
+// non-default φ (the sum would stay at its natural value, not π − 0.1).
+//
+// Enabled 2026-05-30: `newton_euclidean` now uses the block-FD edge-DOF Hessian
+// (`euclidean_hessian_block_fd_sym`) for cyclic layouts, so the full solve runs.
+// ════════════════════════════════════════════════════════════════════════════
+TEST(EuclideanFunctional, CyclicCircularEdge_CatHead_JavaXVal)
+{
+    const std::string path = std::string(CONFORMALLAB_DATA_DIR) + "/obj/cathead.obj";
+    ConformalMesh mesh;
+    ASSERT_NO_THROW(mesh = load_mesh(path)) << "cathead.obj not found: " << path;
+
+    auto maps = setup_euclidean_maps(mesh);
+    compute_euclidean_lambda0_from_mesh(mesh, maps);
+
+    // DOFs: interior vertices (border pinned) + exactly ONE edge DOF on the
+    // "circular" edge (Java marks a single circularHoleEdge). Giving every edge
+    // a DOF would make λ_e redundant with u_i+u_j (a V-dim null space) and stall
+    // Newton; the single extra edge variable keeps the system well-posed.
+    int idx = 0;
+    for (auto v : mesh.vertices())
+        maps.v_idx[v] = mesh.is_border(v) ? -1 : idx++;
+
+    // Pick one interior edge: both incident faces present, both endpoints interior.
+    Edge_index circular{};
+    bool found = false;
+    for (auto e : mesh.edges()) {
+        auto h  = mesh.halfedge(e);
+        auto ho = mesh.opposite(h);
+        if (mesh.is_border(h) || mesh.is_border(ho)) continue;
+        if (mesh.is_border(mesh.source(h)) || mesh.is_border(mesh.target(h))) continue;
+        circular = e; found = true; break;
+    }
+    ASSERT_TRUE(found) << "no interior edge found on cathead";
+    maps.e_idx[circular] = idx++;          // the single circular-edge DOF
+    const int n = idx;
+    ASSERT_GT(n, 0);
+    const std::size_t ie = static_cast<std::size_t>(maps.e_idx[circular]);
+
+    // Natural targets: set Θ_v / φ_e so that x = 0 is the equilibrium (G(0)=0),
+    // so the only deviation is the prescribed circular φ below.
+    std::vector<double> x0(static_cast<std::size_t>(n), 0.0);
+    auto G0 = euclidean_gradient(mesh, x0, maps);
+    for (auto v : mesh.vertices()) {
+        int iv = maps.v_idx[v];
+        if (iv >= 0) maps.theta_v[v] -= G0[static_cast<std::size_t>(iv)];
+    }
+    maps.phi_e[circular] += G0[ie];
+
+    // Prescribe the circular edge turn angle φ = π − 0.1 (Java CustomEdgeInfo.phi).
+    const double phi_target = PI - 0.1;
+    maps.phi_e[circular] = phi_target;
+
+    auto res = newton_euclidean(mesh, x0, maps, /*tol=*/1e-11, /*max_iter=*/200);
+    ASSERT_TRUE(res.converged)
+        << "Newton did not converge; ||G||=" << res.grad_inf_norm;
+
+    // Realised geometric opposite-corner-angle sum = φ_e + G_e(x*)  (= α_opp+α_opp).
+    auto Gf = euclidean_gradient(mesh, res.x, maps);
+    const double realised = maps.phi_e[circular] + Gf[ie];
+    EXPECT_NEAR(phi_target, realised, 1e-9)
+        << "prescribed circular edge turn angle π−0.1 not realised at the solution";
+}
+
+// ════════════════════════════════════════════════════════════════════════════
+// Correctness of the block-FD edge-DOF (cyclic) Hessian
+//
+// Validates `euclidean_hessian_block_fd` directly: on a tetrahedron with the
+// full cyclic DOF layout (4 vertex + 6 edge), every entry must match the
+// column-wise finite difference of `euclidean_gradient` (the true Jacobian of
+// G). This pins the per-face output sign mapping (−α vertex / +α_opp edge) and
+// the scatter independently of the convergence test above.
+// ════════════════════════════════════════════════════════════════════════════
+TEST(EuclideanFunctional, CyclicHessian_BlockFD_MatchesGradientFD_Tetrahedron)
+{
+    auto mesh = make_tetrahedron();
+    auto maps = setup_euclidean_maps(mesh);
+    compute_euclidean_lambda0_from_mesh(mesh, maps);
+    const int n = assign_euclidean_all_dof_indices(mesh, maps);  // 4 + 6 = 10
+
+    std::vector<double> x(static_cast<std::size_t>(n), 0.0);
+    for (int i = 0; i < 4; ++i) x[static_cast<std::size_t>(i)] = -0.15; // vertices
+    for (int i = 4; i < n; ++i) x[static_cast<std::size_t>(i)] =  0.05; // edges
+
+    auto H = euclidean_hessian_block_fd(mesh, x, maps);
+
+    const double eps = 1e-6;
+    std::vector<double> xp = x, xm = x;
+    double max_err = 0.0;
+    for (int j = 0; j < n; ++j) {
+        const std::size_t sj = static_cast<std::size_t>(j);
+        xp[sj] = x[sj] + eps;
+        xm[sj] = x[sj] - eps;
+        auto Gp = euclidean_gradient(mesh, xp, maps);
+        auto Gm = euclidean_gradient(mesh, xm, maps);
+        xp[sj] = xm[sj] = x[sj];
+        for (int i = 0; i < n; ++i) {
+            const double fd = (Gp[static_cast<std::size_t>(i)]
+                             - Gm[static_cast<std::size_t>(i)]) / (2.0 * eps);
+            max_err = std::max(max_err, std::abs(H.coeff(i, j) - fd));
+        }
+    }
+    EXPECT_LT(max_err, 1e-5)
+        << "block-FD cyclic Hessian disagrees with the gradient finite difference";
+}
+
+// ════════════════════════════════════════════════════════════════════════════
+// Analytic cyclic Hessian == block-FD (and == gradient FD)
+//
+// `euclidean_hessian_analytic` is the closed-form counterpart of
+// `euclidean_hessian_block_fd`.  On the full cyclic layout they must agree to
+// round-off, and both must match the gradient finite difference.  This validates
+// the analytic derivation (∂α_i/∂s_i = ℓ_i²/4A, ∂α_i/∂s_j = ½cot α_i − ℓ_j²/4A).
+// ════════════════════════════════════════════════════════════════════════════
+TEST(EuclideanFunctional, CyclicHessian_Analytic_MatchesBlockFD_Tetrahedron)
+{
+    auto mesh = make_tetrahedron();
+    auto maps = setup_euclidean_maps(mesh);
+    compute_euclidean_lambda0_from_mesh(mesh, maps);
+    const int n = assign_euclidean_all_dof_indices(mesh, maps);  // 4 + 6 = 10
+
+    std::vector<double> x(static_cast<std::size_t>(n), 0.0);
+    for (int i = 0; i < 4; ++i) x[static_cast<std::size_t>(i)] = -0.15;
+    for (int i = 4; i < n; ++i) x[static_cast<std::size_t>(i)] =  0.05;
+
+    auto Ha = euclidean_hessian_analytic(mesh, x, maps);
+    auto Hb = euclidean_hessian_block_fd(mesh, x, maps);
+
+    // Analytic vs block-FD.
+    double max_ab = 0.0;
+    for (int i = 0; i < n; ++i)
+        for (int j = 0; j < n; ++j)
+            max_ab = std::max(max_ab, std::abs(Ha.coeff(i, j) - Hb.coeff(i, j)));
+    EXPECT_LT(max_ab, 1e-6)
+        << "analytic cyclic Hessian disagrees with block-FD";
+
+    // Analytic vs gradient FD (independent ground truth).
+    const double eps = 1e-6;
+    std::vector<double> xp = x, xm = x;
+    double max_af = 0.0;
+    for (int j = 0; j < n; ++j) {
+        const std::size_t sj = static_cast<std::size_t>(j);
+        xp[sj] = x[sj] + eps;
+        xm[sj] = x[sj] - eps;
+        auto Gp = euclidean_gradient(mesh, xp, maps);
+        auto Gm = euclidean_gradient(mesh, xm, maps);
+        xp[sj] = xm[sj] = x[sj];
+        for (int i = 0; i < n; ++i) {
+            const double fd = (Gp[static_cast<std::size_t>(i)]
+                             - Gm[static_cast<std::size_t>(i)]) / (2.0 * eps);
+            max_af = std::max(max_af, std::abs(Ha.coeff(i, j) - fd));
+        }
+    }
+    EXPECT_LT(max_af, 1e-5)
+        << "analytic cyclic Hessian disagrees with the gradient finite difference";
+
+    // Symmetry (Hessian of a scalar energy).
+    double max_sym = 0.0;
+    for (int i = 0; i < n; ++i)
+        for (int j = 0; j < n; ++j)
+            max_sym = std::max(max_sym, std::abs(Ha.coeff(i, j) - Ha.coeff(j, i)));
+    EXPECT_LT(max_sym, 1e-9) << "analytic cyclic Hessian is not symmetric";
+}

doc/reviewer/java-ignore-crossvalidation.md

@@ -0,0 +1,257 @@
+# Cross-validation potential of the Java `@Ignore`d tests
+
+> **Purpose.** The upstream Java library (`de.varylab.discreteconformal`) has
+> 22 test files with `@Ignore`d tests.  This document classifies them by *why*
+> they are ignored and by their *cross-validation value for the C++ port*, so
+> that the highest-value golden oracles get ported first.
+>
+> **Created:** 2026-05-29 (external-reviewer pass).  See also
+> `doc/reviewer/java-port-audit.md` and PR #27 (`tests/java-golden-oracles`).
+
+---
+
+## 1 — Two root causes for `@Ignore`
+
+| Cause | Count | Meaning for the C++ port |
+|---|---|---|
+| **`requires ARM64 PETSc native lib`** (class-level) | 15 files | **Not a logic problem.** Java solves via PETSc/TAO (native lib) which does not run on ARM64. The C++ port uses **Eigen** (no PETSc) and *can* run these scenarios — **iff** the feature is ported **and** the Java test carries a hard-coded expected value. |
+| **Logic / known-issue** (method-level) | ~8 methods | Mostly mark **Java gaps** (no Hessian, layout unfinished) or **ARM64 FP flakiness**. Usually **not** reliable oracles. |
+
+A critical distinction: PR #27 cross-validates **functional evaluations**
+(energy/gradient — no solver, so no PETSc needed). The **full-solve** golden
+values exist **only** as hard-coded vectors recorded from a historical x86
+PETSc run — and those live exclusively in `HyperIdealConvergenceTest`.
+
+---
+
+## 2 — Already locked by PR #27 (do not duplicate)
+
+`tests/java-golden-oracles` pins, bit-for-bit (1e-12), against compiled Java:
+- `HyperIdealGoldenJava` — Clausen/Л/ImLi₂, ζ₁₃/₁₄/₁₅/ζ, both tetrahedron-volume formulas.
+- `EuclideanGoldenJava` / `SphericalGoldenJava` — angle formulas, β relations, Л energy terms, **small full-mesh oracles** on a tetrahedron (gradient `Θ−Σα`, `ΔE = E(x)−E(0)`).
+- `SphericalGoldenJava.FullMeshEdgeDofGradient` — edge-DOF gradient.
+- `PeriodMatrix.NormalizeModulus_GoldenJava` — τ-reduction fold convention.
+
+These cover the **unit math cores** and **functional evaluation**. They do
+**not** cover **solver/convergence** or **full uniformization**.
+
+---
+
+## 3 — Priority ranking for new cross-validation
+
+### 🟢 HIGH — real golden oracle, feature already in C++
+
+| Item | Java source | Oracle | Blocker / effort |
+|---|---|---|---|
+| **HyperIdeal Lawson convergence** ✅ | `HyperIdealConvergenceTest.testHyperIdealConvergence` | **hard-coded golden u\*** (perfectly symmetric — robust to DOF order) | **DONE — PR #30** via low-level half-edge `make_lawson_square_tiled` (CGAL `add_face` can't: multi-edges, §5). |
+
+Golden solution (Lawson square-tiled, BLMVM, tol 1e-10), 22 DOFs:
+```
+vertices (×4,  θ=2π)   : 1.1462158341786262
+edges    (×12, θ=π/2)  : 1.7627471737467797 / …7866   (symmetric pair)
+edges    (×6,  θ=π aux): 2.633915794495759
+```
+Assertion strategy: because the values are symmetric per DOF-class, assert
+*per class* (every vertex DOF ≈ 1.14621…, every π/2-edge ≈ 1.76274…, every
+π-edge ≈ 2.63391…) to 1e-6 — **no DOF-order fidelity required**.
+
+### 🟡 MEDIUM — golden/strong oracle, needs asset or near-term phase
+
+| Item | Java source | Oracle | Blocker |
+|---|---|---|---|
+| **Genus-1 uniformization** | `UniformizerTest.testGenus1` (Wente, 1240 V) | angle-sum = 2π + vertex count (self-consistency) | needs `wente.obj` asset + a "load OBJ → uniformize" path; C++ Euclidean solver exists. |
+| **HyperIdeal Lawson + branch points** ✅ | `…WithBranchPoints` | hard-coded golden u\* | **DONE — PR #30** (stellar subdivision + ideal centres; per-class symmetric @1e-4). |
+| **HyperIdeal layout** | `HyperIdealLayoutTest` | layout positions | needs the generator + `layout.hpp` hyper-ideal path exercised. |
+| **Euclidean big-model scale** | `EuclideanUnwrapperPetscTest.testStaticUnwrapBigModel` | mostly self-check | scale/robustness cross-check; low oracle strength. |
+
+### 🔵 LATER — feature not yet in C++ (Phase 9d.2 / 10c / 13 / 11)
+
+| Item | Java source | Maps to |
+|---|---|---|
+| Genus-2 uniformization | `UniformizerTest.testGenus2` (Lawson 2498) | Phase 10c / 13 |
+| HyperIdeal hyperelliptic convergence | `…HyperEllipticCurveLawson` (has golden u\*) | Phase 13 |
+| Hyperbolic disk uniformization | `HyperbolicDiskUniformizationTest` | Phase 10c (gated on holonomy bug) |
+| Möbius centring (grad/Hessian/convergence golden) | `MobiusCenteringFunctionalTest` | Phase 9d.4 |
+| Quasi-isothermic suite (SinCondition, SNES, Utility, Layout, ConformalStructure) | `…quasiisothermic.*` | Phase 10e |
+| Schottky uniformisation | `SchottkyUtilityTest` | Phase 11a |
+| Circle-domain (multiply-connected) | `CircleDomainUnwrapperTest` | Phase 11c |
+| Electrostatic-sphere baseline | `ElectrostaticSphereFunctionalTest` | Phase 10c′ (optional) |
+
+### ⚪ NON-ORACLE — do **not** port as cross-validation
+
+| Item | Why it is not an oracle |
+|---|---|
+| `HyperIdealFunctionalTest.testHessian` (`@Ignore@Override`, empty) | Java has **no** HyperIdeal Hessian (`hasHessian()==false`). C++ *adds* it (block-FD, Phase 9b) → validate against C++ FD, not Java. Documents the C++ delta. |
+| `ElectrostaticSphereFunctionalTest.testHessian` | same — Hessian not implemented in Java. |
+| `HyperbolicLayoutTest.testDoLayout` | throws `"Not implemented correctly yet!"` — **Java itself unfinished**; no reference to validate against. |
+| `SurfaceCurveUtilityTest.testCreateSurfaceCurves` | known issue: non-deterministic edge count on ARM64 — **flaky**, not a stable oracle. |
+| `HyperbolicCyclicConvergenceTest.testEuclideanConvergence` | known issue: ARM64 FP accumulation, angle sum outside 1e-8 — **flaky**. |
+| `HyperIdealUtilityTest.testZeta13` | asserts ζ₁₃(0.1,0.1,0.1)=0.1, which looks **wrong**; the other volume/ζ tests already run and are ported. Verify the correct value before any use. |
+| `HyperIdealPluginTest` | plugin/TAO-application wiring, not pure math. |
+| `TwoHoleExample.testLetterB` | demo/example, no assertion of record. |
+| `SphericalUnwrapperTest` PETSc variants | the **MTJ** variants (not ignored) already cover the same cases on `cathead`. |
+
+---
+
+## 3b — Effort-adjusted tiering (2026-05-29 inventory)
+
+> A broader scan (not just `@Ignore`d tests) revealed that the §3 ranking was
+> **effort-blind**: cheaper, equally-valuable oracles exist than the Lawson
+> generator. Re-tiered by **value / cost**:
+
+| Tier | Oracle | Mesh / asset | Solver | Golden | C++ status | Cost |
+|---|---|---|---|---|---|---|
+| **1** | `EuclideanCyclicConvergenceTest` | `cathead.obj` *(already in C++)* | MTJ BiCGstab (no PETSc) | `α_opp+α_opp = π−0.1` @1e-12 (per-edge φ target) | `euclidean_functional` ✅, `phi_e` ✅, cathead ✅ | **near-zero** |
+| **1** | `SphericalConvergenceTest` | octahedron | MTJ Newton (no PETSc) | self-consistency (const. curvature) | spherical ✅, `make_octahedron_face` ✅ | low |
+| **2** | Wente genus-1 uniformization | `wente_torus02.obj` + `wente_uniformization.xml` | — (stored result) | exact `IsometryPSL2R` uniformizing-group generators | holonomy/`MobiusMap` ✅; **needs XML reader** | medium |
+| ~~2~~ → **4** | `regular_uniformization.xml`, `LetterB01_canonical_group.xml` | (stored) | — | full `IsometryPSL2R` Fuchsian group + 12-edge `FundamentalPolygon` | **genus ≥ 2 hyperbolic** (not genus-1) → Phase 13 | future |
+| **3** ✅ | HyperIdeal Lawson convergence — plain + branch-points | low-level half-edge `make_lawson_square_tiled` | `newton_hyper_ideal` | hard-coded golden u\* | **DONE — PR #30** (`test_lawson_hyperideal.cpp`) | done |
+| **3** ⏸️ | HyperIdeal Lawson **hyperelliptic** | `lawson_curve_source.xml` + `HyperIdealHyperellipticUtility` | — | non-symmetric golden u\* | **deferred** — needs XML-format reader + jReality-geometry θ port | future |
+| **4** | genus-2 (`genus2.xml`/`lawson2498.obj`), holomorphic forms (`HolomorphicEuclideanUniformization_*.xml`), Schottky (`genus2.xml`) | assets present in Java | — | various | Phases 10a/10c/11/13 | future |
+
+**Tier-2 status (2026-05-30): exhausted.** The only genus-1 case is Wente
+(deferred to Phase 9f, quad/cyclic). `regular_uniformization.xml` and
+`LetterB01_canonical_group.xml` are **genus ≥ 2 hyperbolic** uniformizations
+(full PSL(2,ℝ) Fuchsian groups + a 12-edge fundamental polygon) → Phase 13, not
+Tier 2.
+
+**Tier-3 status (2026-05-30): DONE for plain + branch-points (PR #30).** The
+Lawson base is a genus-2 surface with **4 vertices / 12 edges → multi-edges**, so
+CGAL `add_face`/OFF/polygon-soup cannot build it; `make_lawson_square_tiled`
+(`test_lawson_hyperideal.cpp`) uses the **low-level Surface_mesh half-edge API**
+(`add_edge` + `set_target`/`set_next`/`set_face`/`set_halfedge`), then
+`triangulate_faces` (plain) or `Euler::add_center_vertex` (branch-points,
+stellar). Two key findings made this clean:
+1. The golden vectors are **symmetric per geometric class**, so the external jtem
+   `Triangulator`/`StellarLinear` orderings need **not** be replicated — any
+   consistent triangulation/stellar reproduces the values.
+2. Java's index-based θ split is actually **geometric** (original vs diagonal;
+   base vs spoke); branch-point centres are **ideal vertices** (b = 0).
+
+`newton_hyper_ideal` converges (from x0 = 1.0, **unconstrained** — Java needed
+bounded BLMVM) to both golden vectors:
+- plain: vertices 1.1462158341786262 / original 1.7627471737467797 / aux 2.633915794495759
+- branch-points: vertices 1.3169579 / base 2.2924317 / spokes 0
+
+**Hyperelliptic (the 3rd Lawson variant) is deferred (option c).** It needs a
+reader for `lawson_curve_source.xml` (Java conformal-data `HalfedgeEmbedding`,
+26 vtx / 144 half-edges, a *cut* surface) **and** a port of
+`HyperIdealHyperellipticUtility.calculateCircleIntersections`, which derives θ via
+jReality `Pn` (elliptic/hyperbolic), `MatrixBuilder` (rotate / hyperbolic
+translate) and stereographic projection. Its golden vector is **not**
+class-symmetric (mixed ±values), so no symmetry shortcut — mesh and θ must match
+exactly. Viable approaches: (A) port the jReality geometry, or (B) run a Java
+dumper for θ + connectivity (no PETSc needed) and hard-code them. Both are
+multi-step with numerical-replication risk → its own future task.
+
+**Key correction:** the Lawson generator (originally flagged HIGH) is actually
+**Tier 3** — higher effort than Tier 1/2. Cheapest first wins:
+
+1. **Tier 1 — `EuclideanCyclicConvergenceTest` on cathead** (asset already
+   present; the C++ edge gradient `G_e = α_opp(f⁺)+α_opp(f⁻) − φ_e` makes the
+   Java assertion `α_opp+α_opp = π−0.1` exactly the converged `G_e = 0`
+   condition with `phi_e = π−0.1`). Fails if the solver ignores a non-default
+   `φ`, so it is a genuine check, not a tautology.
+2. **Tier 1 — `SphericalConvergenceTest`** on the octahedron.
+3. **Tier 2 — Wente genus-1 uniformization** — deferred to Phase 9f (quad/cyclic
+   pipeline; `wente_torus02.obj` is a quad mesh, not triangle).
+4. **Tier 3 — Lawson generator** — ✅ DONE (PR #30): plain + branch-points golden
+   oracles GREEN; hyperelliptic deferred (option c).
+
+### Status summary (2026-05-30)
+- **Tier 1** ✅ done (PR #29): cyclic φ oracle + full analytic edge-DOF Hessian.
+- **Tier 2** — exhausted; Wente → Phase 9f, the rest → Phase 13.
+- **Tier 3** ✅ done (PR #30): plain + branch-points Lawson golden oracles;
+  hyperelliptic deferred.
+
+### Build-verification finding (2026-05-29) → resolved (2026-05-30)
+
+Attempting the Tier-1 `EuclideanCyclicConvergenceTest` port **build-verified** a
+blocker first: `newton_euclidean` → `euclidean_hessian` threw *"edge DOFs are
+not supported"*. The cyclic functional needs **vertex + edge** DOFs (the
+*gradient* supported them; the *analytic cotangent Hessian* did not). That is
+why PR #27 cross-validated only the cyclic *evaluation*, never a solve.
+
+**Resolved (2026-05-30):** implemented a **block-FD edge-DOF Euclidean Hessian**
+(`euclidean_hessian_block_fd` / `_sym` in `euclidean_hessian.hpp`, mirroring
+`hyper_ideal_hessian_block_fd`); `newton_euclidean` now routes cyclic layouts
+(edge DOFs present) through it, vertex-only layouts still use the analytic
+cotangent Laplacian.
+
+**Landed (`test_euclidean_functional.cpp`, all GREEN, 243/243 cgal tests pass):**
+- ✅ `CyclicCircularEdge_PhiEntersGradient_CatHead` — solver-free check that the
+  circular-edge φ enters the gradient exactly (`ΔG_e = −Δφ_e` @1e-12).
+- ✅ `CyclicCircularEdge_CatHead_JavaXVal` — **the full Java convergence
+  oracle**: prescribe φ = π−0.1 on one interior ("circular") edge of cathead,
+  solve the cyclic Newton, assert the realised `α_opp+α_opp = π−0.1` @1e-9.
+- ✅ `CyclicHessian_BlockFD_MatchesGradientFD_Tetrahedron` — the block-FD edge-DOF
+  Hessian matches the column-wise gradient FD on the full cyclic layout.
+
+> DOF note: only the **single** circular edge gets an edge DOF (as in Java's
+> one `circularHoleEdge`). Giving *every* edge a DOF makes λ_e redundant with
+> u_i+u_j (a V-dim null space) and stalls Newton.
+
+Spherical convergence (Tier-1 #2) was found **already covered** by
+`test_newton_solver` (`Spherical_ConvergesFromPerturbation` et al.), and
+`make_octahedron_face()` is a single face, not a closed octahedron — so it adds
+little.
+
+### Tier-2 (Wente) finding (2026-05-30) — needs a quad/cyclic pipeline
+
+`wente_torus02.obj` is a **quad mesh** (1240 quads, genus 1), and the Java
+`wente_uniformization.xml` golden (uniformizing group → τ = (√2+i√6)/(2√2) =
+½ + i·√3/2, the hexagonal torus) comes from Java's **cyclic (quad) net**
+uniformization. The C++ period-matrix pipeline (`run_torus_pipeline`) is
+**triangle-based**, so a faithful bit-vs-Java comparison would require triangulating
+(which changes the discrete conformal structure → no exact match) or a
+**quad/cyclic period-matrix pipeline** — that belongs to **Phase 9f** (general
+polygon meshes), not the triangle path. → Tier 2 deferred to Phase 9f; the
+hexagonal golden τ is recorded here for when that lands.
+
+### Done instead (2026-05-30): full analytic edge-DOF Hessian
+With Tier-2 gated on a quad pipeline, the higher-value clean win was upgrading
+the cyclic Hessian from block-FD to **fully analytic** (closed-form
+`∂α_i/∂s_i = ℓ_i²/4A`, `∂α_i/∂s_j = ½cot α_i − ℓ_j²/4A`), satisfying the
+`novelty-statement.md §3.2` "analytic Hessians, not finite difference" claim for
+the Euclidean path. Verified against block-FD and gradient-FD.
+
+### Other generators in the Java tree
+`HyperellipticCurveGenerator` (→ Phase 13), `SchottkyGenerator` (→ Phase 11a),
+`FunctionalTest.createOctahedron/createTetrahedron` (trivial — C++ has
+`make_octahedron_face` / `make_tetrahedron`).
+
+---
+
+## 4 — Recommended order of work
+
+1. ✅ **DONE (PR #29)** — Tier 1: cyclic φ oracle + full analytic edge-DOF Hessian.
+2. ✅ **DONE (PR #30)** — Tier 3: `make_lawson_square_tiled()` (low-level
+   half-edge; §5) + plain & branch-points golden oracles (per-class symmetric,
+   no DOF-order matching needed).
+3. ⏸️ **Deferred** — Lawson **hyperelliptic** variant (XML-format reader +
+   `HyperIdealHyperellipticUtility` / jReality θ; non-symmetric golden).
+4. ⏸️ **Phase 9f** — Wente genus-1 uniformization (quad/cyclic pipeline).
+5. **Phase 10a/10c/11/13** — LATER items (genus-2, holomorphic forms, Schottky);
+   golden u\* / data recorded in the Java tests for reuse then.
+
+---
+
+## 5 — Why the Lawson mesh needs a low-level generator
+
+`HyperIdealGenerator.createLawsonSquareTiledBase()` builds a **genus-2**
+surface with **4 vertices and 12 edges** (24 half-edges), then triangulates
+6 quad faces (adding 6 diagonal "aux" edges, θ=π; the 12 originals get θ=π/2).
+
+With only 4 vertices and 12 edges, **most vertex pairs carry more than one
+edge** (multi-edges). CGAL `Surface_mesh::add_face` and polygon-soup builders
+key edges by their endpoint vertex pair and therefore **cannot represent
+multi-edges** — they would either fail or glue the surface into the wrong
+topology. Java sidesteps this with explicit half-edge IDs and
+`linkOppositeEdge` / `linkNextEdge`.
+
+The C++ port must replicate this with the **low-level half-edge API**
+(`add_vertex`, `add_edge`, and explicit `set_next` / `set_target` /
+`set_face` connectivity), mirroring the exact opposite/next links from
+`createLawsonSquareTiledBase()` (edges 0–23). The symmetric golden vector
+then makes the resulting test robust even though CGAL's DOF ordering differs
+from Java's.