// Copyright (c) 2024-2026 Tarik Moussa.
// SPDX-License-Identifier: MIT

// test_euclidean_functional.cpp
//
// Phase 3d — EuclideanCyclicFunctional ported to ConformalMesh.
//
// Corresponds to de.varylab.discreteconformal.functional.EuclideanCyclicFunctionalTest.
//
// Test map (Java → C++)
// ──────────────────────
//   testHessian (Ignored)       → GradientCheck_Hessian           (ported)
//   testGradient…Triangle       → GradientCheck_TriangleVertex     (ported)
//   testGradient…QuadStrip      → GradientCheck_QuadStripVertex    (ported)
//   testGradient…Tetrahedron    → GradientCheck_TetrahedronVertex  (ported)
//   testGradient…AllDofs        → GradientCheck_TetrahedronAllDofs (ported)
//   testFunctionalAtNaNValue    → AnglesFiniteAtKnownPoint          (ported)
//
// Energy model
// ────────────
//   Uses the Schläfli path integral E(x) = ∫₀¹⟨G(tx),x⟩dt (10-point GL).
//   The gradient check verifies G is curl-free.

#include "conformal_mesh.hpp"
#include "mesh_builder.hpp"
#include "euclidean_geometry.hpp"
#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;

// ════════════════════════════════════════════════════════════════════════════
// Cross-module Hessian check: euclidean_gradient() ↔ euclidean_hessian()
//
// Java @Ignore reason: "no Hessian implemented yet" — the Java functional
// test was written before the Hessian existed.  In C++ the analytic
// cotangent-Laplace Hessian (euclidean_hessian.hpp, Phase 3f) is complete.
//
// This test verifies cross-module consistency:
//   H[i,j] ≈ (G_i(x+ε·eⱼ) − G_i(x−ε·eⱼ)) / (2ε)
// using the gradient from euclidean_functional.hpp and the Hessian from
// euclidean_hessian.hpp.  A bug in DOF-index mapping or sign convention
// that affects both modules independently would only be caught here.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, GradientCheck_Hessian)
{
    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);

    // hessian_check_euclidean: H[i,j] ≈ FD(G)[i,j] using euclidean_gradient()
    EXPECT_TRUE(hessian_check_euclidean(mesh, x, maps))
        << "Cross-module: euclidean_gradient() and euclidean_hessian() are inconsistent";
}

// ════════════════════════════════════════════════════════════════════════════
// Angle formula: equilateral triangle → all angles = π/3
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, EquilateralTriangleAnglesArePiOver3)
{
    // All sides equal: l = 1.0, log-length = 0.
    auto fa = euclidean_angles(0.0, 0.0, 0.0);

    ASSERT_TRUE(fa.valid) << "Equilateral triangle must be valid";
    constexpr double PI_3 = 3.14159265358979323846 / 3.0;
    EXPECT_NEAR(fa.alpha1, PI_3, 1e-12);
    EXPECT_NEAR(fa.alpha2, PI_3, 1e-12);
    EXPECT_NEAR(fa.alpha3, PI_3, 1e-12);
}

// ════════════════════════════════════════════════════════════════════════════
// Angle formula: right isosceles triangle (legs 1, hypotenuse √2)
//
// For the 45-45-90 triangle: angles are π/4, π/4, π/2.
// From make_triangle default (v0=(0,0), v1=(1,0), v2=(0,1)):
//   e01: l=1, λ°=0
//   e12: l=√2, λ°=log(2)
//   e02: l=1, λ°=0
// Angle at v0 (opposite e12) = π/2.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, RightIsoscelesTriangleAnglesCorrect)
{
    const double log2 = std::log(2.0);
    // lam12 = 0 (v0-v1, length 1), lam23 = log(2) (v1-v2, length √2), lam31 = 0 (v2-v0, length 1)
    // v1 = v0 in our ordering → remap: l01=1, l12=√2, l20=1
    // Using euclidean_angles(lam_v1v2, lam_v2v3, lam_v3v1):
    //   v1=(0,0), v2=(1,0), v3=(0,1)
    //   lam12 = log(1²) = 0, lam23 = log(√2 ²) = log2, lam31 = log(1²) = 0
    auto fa = euclidean_angles(0.0, log2, 0.0);

    ASSERT_TRUE(fa.valid);
    constexpr double PI   = 3.14159265358979323846;
    // v1=(0,0) is at the right-angle corner (opposite the hypotenuse l23=√2) → α1 = 90°.
    // v2=(1,0) and v3=(0,1) are the 45° corners (each opposite a leg of length 1).
    EXPECT_NEAR(fa.alpha1, PI / 2.0, 1e-12);   // angle at v1 (opposite l23=√2): 90°
    EXPECT_NEAR(fa.alpha2, PI / 4.0, 1e-12);   // angle at v2 (opposite l31=1):  45°
    EXPECT_NEAR(fa.alpha3, PI / 4.0, 1e-12);   // angle at v3 (opposite l12=1):  45°
}

// ════════════════════════════════════════════════════════════════════════════
// Angle sum = π for any valid Euclidean triangle
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, AngleSumEqualsPi)
{
    // Scalene triangle with log-lengths (0, 0.5, -0.3).
    auto fa = euclidean_angles(0.0, 0.5, -0.3);

    ASSERT_TRUE(fa.valid);
    constexpr double PI = 3.14159265358979323846;
    EXPECT_NEAR(fa.alpha1 + fa.alpha2 + fa.alpha3, PI, 1e-12);
}

// ════════════════════════════════════════════════════════════════════════════
// Degenerate triangle → valid = false
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, DegenerateTriangleReturnsFalse)
{
    // l12 = l23 = 1, l31 = 3 → violates triangle inequality.
    auto fa = euclidean_angles_from_lengths(1.0, 1.0, 3.0);
    EXPECT_FALSE(fa.valid);
}

// ════════════════════════════════════════════════════════════════════════════
// Gradient check: default right-isosceles triangle, vertex DOFs only
//
// Mirrors Java testGradient…SingleTriangle.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, GradientCheck_TriangleVertex)
{
    auto mesh = make_triangle();   // (0,0)–(1,0)–(0,1)
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);
    int n = assign_euclidean_vertex_dof_indices(mesh, maps);

    // Small uniform conformal perturbation.
    std::vector<double> x(static_cast<std::size_t>(n), -0.1);

    EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps))
        << "Gradient check failed on right-isosceles triangle (vertex DOFs)";
}

// ════════════════════════════════════════════════════════════════════════════
// Gradient check: quad strip (2 triangles, 1 interior edge), vertex DOFs only
//
// Mirrors Java testGradient…QuadStrip / testGradientInExtendedDomain.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, GradientCheck_QuadStripVertex)
{
    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.2);

    EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps))
        << "Gradient check failed on quad strip (vertex DOFs)";
}

// ════════════════════════════════════════════════════════════════════════════
// Gradient check: regular tetrahedron, vertex DOFs only
//
// Closed surface (4 faces, 4 vertices, 6 interior edges).
// Exercises per-vertex angle-sum accumulation on multiple faces.
// Mirrors Java testGradient…Tetrahedron / testGradientWithHyperIdeal…
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, GradientCheck_TetrahedronVertex)
{
    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(gradient_check_euclidean(mesh, x, maps))
        << "Gradient check failed on regular tetrahedron (vertex DOFs)";
}

// ════════════════════════════════════════════════════════════════════════════
// Gradient check: tetrahedron, all DOFs (vertex + edge)
//
// Exercises the edge-gradient branch G_e = α_opp⁺ + α_opp⁻ − π.
// Mirrors Java testGradientWithHyperellipticCurve.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, GradientCheck_TetrahedronAllDofs)
{
    auto mesh = make_tetrahedron();
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);
    int n = assign_euclidean_all_dof_indices(mesh, maps);

    // 4 vertex DOFs + 6 edge DOFs = 10 total.
    std::vector<double> x(static_cast<std::size_t>(n), 0.0);
    // Set vertex DOFs slightly negative to keep triangles non-degenerate.
    for (int i = 0; i < 4; ++i)
        x[static_cast<std::size_t>(i)] = -0.15;

    EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps))
        << "Gradient check failed on regular tetrahedron (all DOFs)";
}

// ════════════════════════════════════════════════════════════════════════════
// Angles are finite at a known interior point
//
// Mirrors Java testFunctionalAtNaNValue: stress-test the angle formula with
// large negative conformal factors (compressed triangle) to ensure no NaN/Inf.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, AnglesFiniteAtKnownPoint)
{
    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);

    // Very compressed: u_i = -3 (all sides shrunk by exp(-3) ≈ 0.05).
    // Triangle stays well-formed (equilateral shrinks uniformly).
    std::vector<double> x(static_cast<std::size_t>(n), -3.0);
    auto G = euclidean_gradient(mesh, x, maps);

    for (std::size_t i = 0; i < G.size(); ++i) {
        EXPECT_FALSE(std::isnan(G[i])) << "Gradient component " << i << " is NaN";
        EXPECT_FALSE(std::isinf(G[i])) << "Gradient component " << i << " is Inf";
    }
}

// ════════════════════════════════════════════════════════════════════════════
// Gradient check: fan of 5 flat triangles, vertex DOFs only
//
// High-valence central vertex: exercises per-vertex angle accumulation
// across 5 incident faces.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, GradientCheck_Fan5Vertex)
{
    auto mesh = make_fan(5);
    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.05);

    EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps))
        << "Gradient check failed on flat fan-5 mesh";
}

// ════════════════════════════════════════════════════════════════════════════
// Gradient check: mixed pinned/variable vertices
//
// Pins the first vertex (u_v0 = 0 fixed), lets the rest be variable.
// Verifies that the gradient accumulator skips pinned vertices correctly.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, GradientCheck_MixedPinnedVertices)
{
    auto mesh = make_quad_strip();
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);

    // Manually pin v0; assign v1, v2, v3 as DOFs 0, 1, 2.
    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.3, -0.2};

    EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps))
        << "Gradient check failed for mixed pinned/variable vertices";
}

// ─────────────────────────────────────────────────────────────────────────────
// Golden-value oracle — pin the Euclidean angle formula and the 2·Л(α) energy
// term bit-for-bit against the upstream Java reference (EuclideanCyclicFunctional
// .triangleEnergyAndAlphas, lines 341-361), captured by running the compiled Java
// library (openjdk 17) with the real de.varylab…Clausen.Л on these exact edge
// lengths.  Companion to HyperIdealGoldenJava and SphericalGoldenJava: locks the
// absolute angle/energy values against an independent implementation, catching
// silent index/sign drift the curl-free path-integral gradient check cannot see.
//
// To regenerate: /tmp/oracle/EucOracle.java.  Values are Java printf %.17g.
// ─────────────────────────────────────────────────────────────────────────────
TEST(EuclideanGoldenJava, AngleAndLobachevskyEnergyFromLengths)
{
    auto check = [](double l12, double l23, double l31,
                    double a1_g, double a2_g, double a3_g, double L_g) {
        auto fa = euclidean_angles_from_lengths(l12, l23, l31);
        EXPECT_TRUE(fa.valid);
        EXPECT_NEAR(fa.alpha1, a1_g, 1e-12);
        EXPECT_NEAR(fa.alpha2, a2_g, 1e-12);
        EXPECT_NEAR(fa.alpha3, a3_g, 1e-12);
        const double Lterm = 2.0 * Lobachevsky(fa.alpha1)
                           + 2.0 * Lobachevsky(fa.alpha2)
                           + 2.0 * Lobachevsky(fa.alpha3);
        EXPECT_NEAR(Lterm, L_g, 1e-12);
    };
    check(1.0, 1.2, 0.9,
          1.3637649752769678, 0.82416964552030680, 0.95365803279251860,
          1.9456273836230942);
    check(1.0, 1.0, 1.0,
          1.0471975511965979, 1.0471975511965979, 1.0471975511965979,
          2.0298832128193070);
}

// ─────────────────────────────────────────────────────────────────────────────
// FULL-MESH golden oracle — the strongest cross-check: drives the REAL upstream
// EuclideanCyclicFunctional (openjdk 17) on a tetrahedron loaded from a shared
// OBJ (identical topology + geometry to make_tetrahedron()), and pins BOTH the
// per-vertex gradient G_v = Θ_v − Σα AND the energy difference ΔE = E(x) − E(0)
// bit-for-bit against it.
//
// This closes audit missing-test item 5 (full-mesh energy + gradient at a known
// x).  It is genuinely independent of the C++ implementation in two ways:
//   • the gradient is the upstream library's own analytic gradient (not an FD
//     check, which only proves curl-freeness of the C++ self-consistent energy);
//   • the energy is Java's CLOSED-FORM functional value, whereas C++ computes it
//     as a Gauss-Legendre PATH INTEGRAL of its gradient — two different methods
//     that must agree on ΔE (the initialEnergy constant and the φ·λ⁰ term cancel
//     in the difference, and there are no edge DOFs here).
//
// Setup parity (verified against UnwrapUtility.prepareInvariantDataEuclidean):
//   closed mesh, ALL 4 vertices variable (no pin), Θ_v = 2π, no edge DOFs,
//   λ°_e = 2·log(|p_i − p_j|), per-vertex u(P) = 0.10·X − 0.07·Y + 0.13·Z.
//
// To regenerate: /tmp/oracle/{tet.obj,EucMeshOracle.java}.  Values are Java %.17g.
// ─────────────────────────────────────────────────────────────────────────────
TEST(EuclideanGoldenJava, FullMeshGradientAndEnergy_Tetrahedron)
{
    constexpr double TWO_PI = 2.0 * 3.14159265358979323846264338328;

    auto mesh = make_tetrahedron();   // same 4 vertices as /tmp/oracle/tet.obj
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);

    // All four vertices are DOFs (Java prepareInvariantData makes every interior
    // vertex variable on a closed mesh); Θ_v = 2π; no edge DOFs.
    int idx = 0;
    for (auto v : mesh.vertices()) {
        maps.v_idx[v]   = idx++;
        maps.theta_v[v] = TWO_PI;
    }

    auto u_of = [](const Point3& p) {
        return 0.10 * p.x() - 0.07 * p.y() + 0.13 * p.z();
    };

    std::vector<double> x(static_cast<std::size_t>(idx), 0.0);
    for (auto v : mesh.vertices())
        x[static_cast<std::size_t>(maps.v_idx[v])] = u_of(mesh.point(v));

    auto G = euclidean_gradient(mesh, x, maps);

    // Java golden gradients keyed by vertex position (order-independent lookup).
    struct GoldRow { double X, Y, Z, G; };
    const GoldRow gold[4] = {
        { 1,  1,  1, 3.5277511803984396},
        { 1, -1, -1, 3.2837611905358440},
        {-1,  1, -1, 2.3437485291237100},
        {-1, -1,  1, 3.4111097143011780},
    };
    for (auto v : mesh.vertices()) {
        const auto& p = mesh.point(v);
        const double g = G[static_cast<std::size_t>(maps.v_idx[v])];
        bool matched = false;
        for (const auto& row : gold) {
            if (std::abs(p.x() - row.X) < 1e-9 &&
                std::abs(p.y() - row.Y) < 1e-9 &&
                std::abs(p.z() - row.Z) < 1e-9) {
                EXPECT_NEAR(g, row.G, 1e-12)
                    << "gradient mismatch at (" << p.x() << "," << p.y()
                    << "," << p.z() << ")";
                matched = true;
                break;
            }
        }
        EXPECT_TRUE(matched) << "unexpected vertex position";
    }

    // ΔE = E(x) − E(0).  C++ path integral vs Java closed-form functional value.
    std::vector<double> x0(static_cast<std::size_t>(idx), 0.0);
    const double dE = euclidean_energy(mesh, x, maps)
                    - 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";
}

// ════════════════════════════════════════════════════════════════════════════
// DOF assignment gauge-vertex overload (Finding-D, external-audit-2026-05-30)
//
// The single-argument assign_euclidean_vertex_dof_indices() overwrites ALL
// v_idx unconditionally — setting a pin *before* the call has no effect.
// The two-argument overload (gauge vertex) pins the requested vertex in a
// single pass.  These tests verify:
//   (a) single-argument: gauge must be set AFTER, not before.
//   (b) two-argument: the gauge vertex gets v_idx = -1, others are sequential.
//   (c) Newton converges correctly when the gauge overload is used.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanDOFAssignment, SingleArg_PinBeforeHasNoEffect)
{
    // Pin first vertex before the call → the call overwrites it → not pinned.
    auto mesh = make_tetrahedron();
    auto maps = setup_euclidean_maps(mesh);
    auto first = *mesh.vertices().begin();

    maps.v_idx[first] = -1;  // set pin BEFORE — should have no effect
    assign_euclidean_vertex_dof_indices(mesh, maps);

    EXPECT_GE(maps.v_idx[first], 0)
        << "Pre-call pin was overwritten: v_idx[first] must be >= 0 after the call";
    EXPECT_EQ(euclidean_dimension(mesh, maps),
              static_cast<int>(mesh.number_of_vertices()))
        << "All vertices should be free after single-arg assign";
}

TEST(EuclideanDOFAssignment, TwoArg_GaugeIsPinnedOthersAreSequential)
{
    auto mesh = make_tetrahedron();
    auto maps = setup_euclidean_maps(mesh);
    auto first = *mesh.vertices().begin();

    int n = assign_euclidean_vertex_dof_indices(mesh, maps, first);

    EXPECT_EQ(maps.v_idx[first], -1)
        << "Gauge vertex must have v_idx = -1";
    EXPECT_EQ(n, static_cast<int>(mesh.number_of_vertices()) - 1)
        << "Returned DOF count must be num_vertices - 1";
    EXPECT_EQ(euclidean_dimension(mesh, maps), n)
        << "euclidean_dimension must equal returned count";

    // All non-gauge vertices must have distinct indices in [0, n).
    std::vector<int> seen;
    for (auto v : mesh.vertices()) {
        int iv = maps.v_idx[v];
        if (v == first) continue;
        EXPECT_GE(iv, 0);
        EXPECT_LT(iv, n);
        seen.push_back(iv);
    }
    std::sort(seen.begin(), seen.end());
    for (int i = 0; i < n; ++i)
        EXPECT_EQ(seen[static_cast<std::size_t>(i)], i)
            << "DOF indices must be sequential 0..n-1";
}

TEST(EuclideanDOFAssignment, TwoArg_NewtonConvergesWithGaugeOverload)
{
    // End-to-end: use the gauge overload, then run Newton — confirms the
    // pinned-vertex DOF layout is consistent with the solver.
    auto mesh = make_tetrahedron();
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);
    auto first = *mesh.vertices().begin();

    int n = assign_euclidean_vertex_dof_indices(mesh, maps, first);

    // Natural-theta: set targets = actual angle sums at x=0 so x*=0 is the solution.
    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)];
    }

    auto res = newton_euclidean(mesh, x0, maps, 1e-10, 50);
    EXPECT_TRUE(res.converged)
        << "Newton did not converge with gauge-overload DOF assignment";
    EXPECT_LT(res.grad_inf_norm, 1e-9);
}