ConformalLabpp/code/tests/cgal/test_euclidean_functional.cpp

// 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 <gtest/gtest.h>
#include <cmath>
#include <vector>

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";
}