adbf682f0f
Finding-F and Finding-G from doc/reviewer/external-audit-2026-05-30.md
(java-port-audit missing-test items 1 and 2).
Finding-F — Degenerate triangle: limiting angles (item 1)
test_euclidean_functional.cpp:
DegenerateTriangle_LimitingAngles_L{12,23,31}TooLong
— verifies α_opposite = π, other two = 0, valid = false
for all three edge-over-long cases
DegenerateTriangle_GradientPicksUpPiCorner
— end-to-end mesh test: forces effective l12 >> l23+l31 via
lambda0, evaluates gradient, asserts G_v3 = π (not 2π from
a skipped degenerate face) and G_v1=G_v2 = 2π
test_spherical_functional.cpp:
DegenerateTriangle_LimitingAngles_S{12,23,31}TooLong
— same coverage for spherical_angles()
Finding-G — euclidean_hessian edge-DOF guard (item 2)
test_euclidean_hessian.cpp:
EdgeDOFGuard_Throws
— assign_euclidean_all_dof_indices + euclidean_hessian → throw
EdgeDOFGuard_VertexOnlyDoesNotThrow
— vertex-only layout → no throw (regression guard)
275/275 CGAL tests pass, 0 failed.
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
704 lines
34 KiB
C++
704 lines
34 KiB
C++
// Copyright (c) 2024-2026 Tarik Moussa.
|
||
// SPDX-License-Identifier: MIT
|
||
|
||
// test_spherical_functional.cpp (Phase 3c + 3e)
|
||
//
|
||
// Phase 3c — SphericalFunctional ported to ConformalMesh.
|
||
//
|
||
// Corresponds to de.varylab.discreteconformal.functional.SphericalFunctionalTest.
|
||
//
|
||
// Test map (Java → C++)
|
||
// ──────────────────────
|
||
// testHessian (Ignored) → GradientCheck_Hessian (ported)
|
||
// testGradientWithHyperIdeal… → GradientCheck_OctaFaceVertex (ported)
|
||
// testGradientInExtendedDomain → GradientCheck_SpherTetVertex (ported)
|
||
// testGradientWithHyperelliptic → GradientCheck_SpherTetAllDofs (ported)
|
||
// testFunctionalAtNaNValue → AnglesFiniteAtKnownPoint (ported)
|
||
//
|
||
// Energy model
|
||
// ────────────
|
||
// The energy is computed as the Schläfli path integral E(x) = ∫₀¹⟨G(tx),x⟩dt
|
||
// using 10-point Gauss-Legendre quadrature. The gradient check therefore
|
||
// verifies that G is curl-free (the integrability / exactness condition of
|
||
// the spherical discrete conformal functional). This is equivalent to the
|
||
// Java FunctionalTest gradient check.
|
||
|
||
#include "conformal_mesh.hpp"
|
||
#include "mesh_builder.hpp"
|
||
#include "spherical_functional.hpp"
|
||
#include "spherical_hessian.hpp"
|
||
#include "spherical_geometry.hpp"
|
||
#include "clausen.hpp"
|
||
#include <gtest/gtest.h>
|
||
#include <cmath>
|
||
#include <vector>
|
||
|
||
using namespace conformallab;
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Cross-module Hessian check: spherical_gradient() ↔ spherical_hessian()
|
||
//
|
||
// Java @Ignore reason: "no Hessian implemented" — the Java functional test
|
||
// was written before the Hessian existed. In C++ the analytic spherical
|
||
// Hessian (spherical_hessian.hpp, Phase 3f) is complete.
|
||
//
|
||
// This test verifies cross-module consistency between the functional and
|
||
// the Hessian module. The spherical Hessian is NSD (negative semi-definite)
|
||
// because the spherical energy is concave — hessian_check_spherical() uses
|
||
// the sign-corrected FD check appropriate for the spherical case.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalFunctional, GradientCheck_Hessian)
|
||
{
|
||
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))
|
||
<< "Cross-module: spherical_gradient() and spherical_hessian() are inconsistent";
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Angle formula: octahedron-face triangle has all angles = π/2
|
||
//
|
||
// The triangle (1,0,0)–(0,1,0)–(0,0,1) has l_ij = π/2 for all edges.
|
||
// Half-angle formula: s = 3π/4, s_ij = π/4 for all three.
|
||
// All angles = π/2 (right-angled spherical triangle).
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalFunctional, OctaFaceAnglesAreRightAngles)
|
||
{
|
||
// l_ij = π/2 for all edges (octahedron face on unit sphere)
|
||
const double l = PI_SPHER / 2.0;
|
||
auto fa = spherical_angles(l, l, l);
|
||
|
||
ASSERT_TRUE(fa.valid) << "Equilateral spherical triangle must be valid";
|
||
EXPECT_NEAR(PI_SPHER / 2.0, fa.alpha1, 1e-12);
|
||
EXPECT_NEAR(PI_SPHER / 2.0, fa.alpha2, 1e-12);
|
||
EXPECT_NEAR(PI_SPHER / 2.0, fa.alpha3, 1e-12);
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Angle sum of a spherical triangle exceeds π (positive curvature)
|
||
//
|
||
// For the spherical tetrahedron face (arccos(−1/3) ≈ 1.9106 per edge):
|
||
// The dihedral angle = arccos(1/3) ≈ 70.53°; by symmetry the face angles
|
||
// (vertex angles of the spherical triangle) are all equal.
|
||
// Angle sum must be > π and equal 3·arccos(1/3) ≈ 3·1.2310 ≈ 3.693 rad.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalFunctional, SpherTetAngleSumExceedsPi)
|
||
{
|
||
// Edge length of spherical tetrahedron face: arccos(−1/3)
|
||
const double l = std::acos(-1.0 / 3.0);
|
||
auto fa = spherical_angles(l, l, l);
|
||
|
||
ASSERT_TRUE(fa.valid);
|
||
EXPECT_GT(fa.alpha1 + fa.alpha2 + fa.alpha3, PI_SPHER)
|
||
<< "Angle sum of spherical triangle must exceed π";
|
||
|
||
// By symmetry all three angles must be equal
|
||
EXPECT_NEAR(fa.alpha1, fa.alpha2, 1e-12);
|
||
EXPECT_NEAR(fa.alpha2, fa.alpha3, 1e-12);
|
||
|
||
// For a regular spherical tetrahedron with edge arccos(−1/3):
|
||
// half-angle: tan(α/2) = √(sin(l/2)/sin(3l/2)) = √3 → α/2 = π/3 → α = 2π/3.
|
||
// (arccos(1/3) ≈ 1.231 is the 3D dihedral angle of a Euclidean tetrahedron, not this.)
|
||
double expected = 2.0 * PI_SPHER / 3.0; // 120°
|
||
EXPECT_NEAR(fa.alpha1, expected, 1e-10);
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Gradient check: octahedron-face triangle, vertex DOFs only
|
||
//
|
||
// Sets λ° from mesh geometry (unit sphere), all u_i = −0.3 (slightly smaller).
|
||
// Mirrors Java testGradientWithHyperIdeal… on a single-triangle mesh.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalFunctional, GradientCheck_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);
|
||
|
||
// Small uniform conformal factor: shrink the triangle slightly.
|
||
std::vector<double> x(static_cast<std::size_t>(n), -0.3);
|
||
|
||
EXPECT_TRUE(gradient_check_spherical(mesh, x, maps))
|
||
<< "Gradient check failed on octahedron-face triangle (vertex DOFs)";
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Gradient check: spherical tetrahedron (4 faces), vertex DOFs only
|
||
//
|
||
// Closed surface; exercises accumulation over multiple faces per vertex.
|
||
// Mirrors Java testGradientInTheExtendedDomain.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalFunctional, GradientCheck_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(gradient_check_spherical(mesh, x, maps))
|
||
<< "Gradient check failed on spherical tetrahedron (vertex DOFs)";
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Gradient check: spherical tetrahedron, all DOFs (vertex + edge)
|
||
//
|
||
// Exercises the edge-gradient branch: G_e = α_opp⁺ + α_opp⁻ − π.
|
||
// Mirrors Java testGradientWithHyperellipticCurve.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalFunctional, GradientCheck_SpherTetAllDofs)
|
||
{
|
||
auto mesh = make_spherical_tetrahedron();
|
||
auto maps = setup_spherical_maps(mesh);
|
||
compute_lambda0_from_mesh(mesh, maps);
|
||
int n = assign_all_spherical_dof_indices(mesh, maps);
|
||
|
||
// Replacement parameterization (Finding 3): when an edge carries a DOF its
|
||
// value *replaces* λ°_ij + u_i + u_j entirely, so Λ_ij = λ_e. Here the edge
|
||
// DOFs stay at 0 and only the vertex DOFs are perturbed; this checks that the
|
||
// gradient is curl-free (energy = Schläfli path integral), not Java-faithfulness
|
||
// of the edge formula — that is locked separately by
|
||
// EdgeGradient_RegularTetClosedForm below.
|
||
std::vector<double> x(static_cast<std::size_t>(n), 0.0);
|
||
// Set vertex DOFs (indices 0..3) to -0.2 to keep triangle well-formed.
|
||
for (int i = 0; i < 4; ++i) x[static_cast<std::size_t>(i)] = -0.2;
|
||
|
||
EXPECT_TRUE(gradient_check_spherical(mesh, x, maps))
|
||
<< "Gradient check failed on spherical tetrahedron (all DOFs)";
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Closed-form oracle for the edge-DOF gradient (Finding 3, missing-test item 4)
|
||
//
|
||
// The FD gradient check above can only confirm that G is conservative — the
|
||
// spherical energy is *defined* as the path integral of G, so the energy↔gradient
|
||
// FD agreement is automatic and CANNOT detect a wrong-but-conservative edge
|
||
// formula. This test instead pins the edge gradient against an independent,
|
||
// closed-form geometric value, so it would fail if the Finding-3 formula
|
||
// (G_e = α_opp⁺ + α_opp⁻ − θ_e, dropping the additive −(S⁺+S⁻)/2 term) ever
|
||
// regressed.
|
||
//
|
||
// Geometry: the regular spherical tetrahedron has all edges a = arccos(−1/3),
|
||
// so by the spherical law of cosines every interior corner angle is
|
||
// cos α = (cos a − cos²a)/sin²a = cos a/(1+cos a) = (−1/3)/(2/3) = −1/2
|
||
// ⇒ α = 2π/3.
|
||
// Each edge is shared by two faces, so both opposite angles equal 2π/3 and
|
||
// G_e = 2π/3 + 2π/3 − θ_e with θ_e = π (default) = π/3.
|
||
//
|
||
// Setup: all edges carry DOFs, set to their λ⁰ (the replacement convention then
|
||
// reproduces the original tetrahedron metric exactly), vertex DOFs left at 0.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalFunctional, EdgeGradient_RegularTetClosedForm)
|
||
{
|
||
const double PI_ = std::acos(-1.0);
|
||
|
||
auto mesh = make_spherical_tetrahedron();
|
||
auto maps = setup_spherical_maps(mesh);
|
||
compute_lambda0_from_mesh(mesh, maps);
|
||
int n = assign_all_spherical_dof_indices(mesh, maps);
|
||
|
||
// Edge DOF = λ⁰ → Λ_ij = λ⁰ → reproduces the arccos(−1/3) tetrahedron.
|
||
// Vertex DOFs stay at 0 (ignored by the replacement convention for DOF edges).
|
||
std::vector<double> x(static_cast<std::size_t>(n), 0.0);
|
||
int n_edge_dofs = 0;
|
||
for (auto e : mesh.edges()) {
|
||
int ie = maps.e_idx[e];
|
||
if (ie >= 0) { x[static_cast<std::size_t>(ie)] = maps.lambda0[e]; ++n_edge_dofs; }
|
||
}
|
||
ASSERT_EQ(n_edge_dofs, 6) << "regular tetrahedron must have 6 edge DOFs";
|
||
|
||
auto G = spherical_gradient(mesh, x, maps);
|
||
|
||
const double expected = PI_ / 3.0; // 2·(2π/3) − π
|
||
for (auto e : mesh.edges()) {
|
||
int ie = maps.e_idx[e];
|
||
if (ie < 0) continue;
|
||
EXPECT_NEAR(G[static_cast<std::size_t>(ie)], expected, 1e-9)
|
||
<< "edge gradient at DOF " << ie
|
||
<< " must equal the closed-form value π/3 (Finding 3)";
|
||
}
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Angles are finite at a known interior point
|
||
//
|
||
// Mirrors Java testFunctionalAtNaNValue: choose DOFs that could hit
|
||
// a degenerate branch (l_ij → 0 or triangle inequality fails) and check
|
||
// the gradient vector is free of NaN/Inf.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalFunctional, AnglesFiniteAtKnownPoint)
|
||
{
|
||
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);
|
||
|
||
// u_i = -1.5: contracts the triangle heavily but stays non-degenerate.
|
||
std::vector<double> x(static_cast<std::size_t>(n), -1.5);
|
||
auto G = spherical_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-4 mesh on unit sphere, vertex DOFs only
|
||
//
|
||
// Make a fan of 4 triangles around the north pole (0,0,1);
|
||
// rim vertices projected onto the equator.
|
||
// Exercises high-valence vertex gradient accumulation.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalFunctional, GradientCheck_SpherFan4Vertex)
|
||
{
|
||
// Build a fan with 4 spherical triangles manually (can't use make_fan
|
||
// directly because those vertices are not on the unit sphere).
|
||
ConformalMesh mesh;
|
||
auto center = mesh.add_vertex(Point3(0, 0, 1)); // north pole
|
||
const int n_rim = 4;
|
||
std::vector<Vertex_index> rim(n_rim);
|
||
const double dtheta = 2.0 * PI_SPHER / n_rim;
|
||
const double phi = PI_SPHER / 4.0; // 45° colatitude
|
||
for (int i = 0; i < n_rim; ++i) {
|
||
double theta = i * dtheta;
|
||
rim[i] = mesh.add_vertex(Point3(
|
||
std::sin(phi) * std::cos(theta),
|
||
std::sin(phi) * std::sin(theta),
|
||
std::cos(phi)));
|
||
}
|
||
for (int i = 0; i < n_rim; ++i)
|
||
mesh.add_face(center, rim[i], rim[(i + 1) % n_rim]);
|
||
|
||
auto maps = setup_spherical_maps(mesh);
|
||
compute_lambda0_from_mesh(mesh, maps);
|
||
int ndof = assign_vertex_dof_indices(mesh, maps);
|
||
|
||
std::vector<double> x(static_cast<std::size_t>(ndof), -0.3);
|
||
|
||
EXPECT_TRUE(gradient_check_spherical(mesh, x, maps))
|
||
<< "Gradient check failed on spherical fan-4 mesh";
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Gradient check: mixed pinned/variable vertices
|
||
//
|
||
// One vertex pinned (u_v = 0 fixed), others variable.
|
||
// Verifies that the gradient accumulation skips pinned vertices correctly.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalFunctional, GradientCheck_MixedPinnedVertices)
|
||
{
|
||
auto mesh = make_octahedron_face();
|
||
auto maps = setup_spherical_maps(mesh);
|
||
compute_lambda0_from_mesh(mesh, maps);
|
||
|
||
// Pin v0, make v1 and v2 variable.
|
||
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.4};
|
||
|
||
EXPECT_TRUE(gradient_check_spherical(mesh, x, maps))
|
||
<< "Gradient check failed for mixed pinned/variable vertices";
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Phase 3e — Gauge-fix for closed spherical surfaces
|
||
//
|
||
// On a closed spherical surface, the functional has a gauge mode:
|
||
// E(u + t·1) is maximised at some t*.
|
||
// At t*, the sum of all vertex gradients equals zero: Σ G_v = 0.
|
||
//
|
||
// Test: start from a point with non-zero ΣG_v, apply the gauge shift,
|
||
// and verify ΣG_v(x + t*·1) ≈ 0.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalFunctional, GaugeFix_SpherTetVertexZerosSumGv)
|
||
{
|
||
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);
|
||
|
||
// Off-gauge starting point: all u_i = -0.5
|
||
std::vector<double> x(static_cast<std::size_t>(n), -0.5);
|
||
|
||
// Compute ΣG_v before gauge shift.
|
||
{
|
||
auto G = spherical_gradient(mesh, x, maps);
|
||
double sum = 0.0;
|
||
for (auto v : mesh.vertices()) {
|
||
int iv = maps.v_idx[v];
|
||
if (iv >= 0) sum += G[static_cast<std::size_t>(iv)];
|
||
}
|
||
// At -0.5 the surface is compressed; ΣG_v should be non-zero.
|
||
EXPECT_NE(sum, 0.0) << "Pre-gauge ΣG_v should be non-zero";
|
||
}
|
||
|
||
// Compute gauge shift and apply.
|
||
double t = spherical_gauge_shift(mesh, x, maps);
|
||
std::vector<double> x_fixed = x;
|
||
for (auto v : mesh.vertices()) {
|
||
int iv = maps.v_idx[v];
|
||
if (iv >= 0)
|
||
x_fixed[static_cast<std::size_t>(iv)] += t;
|
||
}
|
||
|
||
// Verify ΣG_v ≈ 0 at the gauge-fixed point.
|
||
{
|
||
auto G = spherical_gradient(mesh, x_fixed, maps);
|
||
double sum = 0.0;
|
||
for (auto v : mesh.vertices()) {
|
||
int iv = maps.v_idx[v];
|
||
if (iv >= 0) sum += G[static_cast<std::size_t>(iv)];
|
||
}
|
||
EXPECT_NEAR(sum, 0.0, 1e-6)
|
||
<< "After gauge fix, Σ G_v should vanish; t* = " << t;
|
||
}
|
||
}
|
||
|
||
TEST(SphericalFunctional, GaugeFix_ApplyInPlace)
|
||
{
|
||
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.3: compressed but inside the valid spherical domain.
|
||
std::vector<double> x(static_cast<std::size_t>(n), -0.3);
|
||
|
||
apply_spherical_gauge(mesh, x, maps);
|
||
|
||
// After in-place gauge fix, ΣG_v must be near 0.
|
||
auto G = spherical_gradient(mesh, x, maps);
|
||
double sum = 0.0;
|
||
for (auto v : mesh.vertices()) {
|
||
int iv = maps.v_idx[v];
|
||
if (iv >= 0) sum += G[static_cast<std::size_t>(iv)];
|
||
}
|
||
EXPECT_NEAR(sum, 0.0, 1e-6)
|
||
<< "apply_spherical_gauge must drive Σ G_v to zero";
|
||
}
|
||
|
||
TEST(SphericalFunctional, GaugeFix_AlreadyAtGaugeReturnsTNearZero)
|
||
{
|
||
// A symmetric, equilateral spherical tetrahedron at x=0 is already
|
||
// at the gauge maximum (by symmetry, ΣG_v = 0).
|
||
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.0);
|
||
double t = spherical_gauge_shift(mesh, x, maps);
|
||
// Symmetric starting point → t* should be very close to 0.
|
||
EXPECT_NEAR(t, 0.0, 1e-5)
|
||
<< "Gauge shift from the symmetric point should be ~0; got " << t;
|
||
}
|
||
|
||
// ─────────────────────────────────────────────────────────────────────────────
|
||
// Golden-value oracle — pin the spherical law-of-cosines angle formula, the β
|
||
// half-angle relations, and the Lobachevsky energy term bit-for-bit against the
|
||
// upstream Java reference (SphericalFunctional.triangleEnergyAndAlphas, lines
|
||
// 401-453), captured by running the compiled Java library (openjdk 17) with the
|
||
// real de.varylab…Clausen.Л on these exact arc lengths. Unlike the Schläfli
|
||
// path-integral gradient check (which only verifies curl-freeness), this locks
|
||
// the absolute angle/β/energy values against an independent implementation,
|
||
// catching any silent index/sign/convention drift in the spherical port.
|
||
//
|
||
// To regenerate: /tmp/oracle/SphereOracle.java (recipe in doc/reviewer/
|
||
// java-port-audit.md). Values are Java printf %.17g. Tolerance 1e-12.
|
||
//
|
||
// Index map (C++ spherical_angles ↔ Java): with l12=lij, l23=ljk, l31=lki the
|
||
// returned alpha1/alpha2/alpha3 are Java αi/αj/αk (angle opposite ljk/lki/lij).
|
||
// ─────────────────────────────────────────────────────────────────────────────
|
||
TEST(SphericalGoldenJava, AngleBetaEnergyFromLengths)
|
||
{
|
||
auto check = [](double lij, double ljk, double lki,
|
||
double ai_g, double aj_g, double ak_g,
|
||
double bi_g, double bj_g, double bk_g, double L_g) {
|
||
auto fa = spherical_angles(lij, ljk, lki);
|
||
EXPECT_TRUE(fa.valid);
|
||
EXPECT_NEAR(fa.alpha1, ai_g, 1e-12);
|
||
EXPECT_NEAR(fa.alpha2, aj_g, 1e-12);
|
||
EXPECT_NEAR(fa.alpha3, ak_g, 1e-12);
|
||
|
||
const double ai = fa.alpha1, aj = fa.alpha2, ak = fa.alpha3;
|
||
const double bi = 0.5 * (PI + ai - aj - ak);
|
||
const double bj = 0.5 * (PI - ai + aj - ak);
|
||
const double bk = 0.5 * (PI - ai - aj + ak);
|
||
EXPECT_NEAR(bi, bi_g, 1e-12);
|
||
EXPECT_NEAR(bj, bj_g, 1e-12);
|
||
EXPECT_NEAR(bk, bk_g, 1e-12);
|
||
|
||
const double Lterm =
|
||
Lobachevsky(ai) + Lobachevsky(aj) + Lobachevsky(ak) +
|
||
Lobachevsky(bi) + Lobachevsky(bj) + Lobachevsky(bk) +
|
||
Lobachevsky(0.5 * (PI - ai - aj - ak));
|
||
EXPECT_NEAR(Lterm, L_g, 1e-12);
|
||
};
|
||
|
||
// Scalene spherical triangle (all Δ > 0, Δ_ijk < 2π).
|
||
check(1.0, 1.2, 0.9,
|
||
1.5305813141072122, 0.99684981272794600, 1.1246069519101605,
|
||
1.2753586015294494, 0.74162710015018330, 0.86938423933239760,
|
||
1.3576100185550408);
|
||
// Equilateral spherical triangle.
|
||
check(0.7, 0.7, 0.7,
|
||
1.1225596283199812, 1.1225596283199812, 1.1225596283199812,
|
||
1.0095165126349062, 1.0095165126349060, 1.0095165126349060,
|
||
1.6806828584976297);
|
||
}
|
||
|
||
// ─────────────────────────────────────────────────────────────────────────────
|
||
// FULL-MESH golden oracle (spherical) — drives the REAL upstream SphericalFunctional
|
||
// on the shared tetrahedron (scaled so every arc length stays < π) and pins BOTH
|
||
// the per-vertex gradient G_v = Θ_v − Σα AND ΔE = E(x) − E(0) bit-for-bit.
|
||
//
|
||
// IMPORTANT — the oracle calls Java's `conformalEnergyAndGradient` (the RAW
|
||
// θ−Σα gradient), NOT `evaluate()`: `evaluate()` first runs a 1-D Brent
|
||
// maximization over the global-scale gauge direction (`maximizeInNegativeDirection`),
|
||
// which the C++ `spherical_gradient` deliberately does NOT — C++ factors that
|
||
// gauge into the Newton solver's `spherical_gauge_shift` instead. Comparing
|
||
// against `evaluate()` would (wrongly) drive every component to ~0. The raw
|
||
// gradient is the piece that corresponds 1:1 to the C++ functional.
|
||
//
|
||
// Setup parity (UnwrapUtility.prepareInvariantDataHyperbolicAndSpherical, scale):
|
||
// closed mesh, ALL 4 vertices variable, Θ_v = 2π, no edge DOFs,
|
||
// λ°_e = 2·log(SCALE·|p_i − p_j|) (Java uses chord length × scale, whereas the
|
||
// C++ compute_lambda0_from_mesh helper assumes unit-sphere vertices and uses
|
||
// ARC length — so we set λ° directly here to match Java exactly),
|
||
// per-vertex u(P) = 0.10·X − 0.07·Y + 0.13·Z, SCALE = 0.2.
|
||
//
|
||
// To regenerate: /tmp/oracle/{tet.obj,SphereMeshOracle.java}. Values are Java %.17g.
|
||
// ─────────────────────────────────────────────────────────────────────────────
|
||
TEST(SphericalGoldenJava, FullMeshGradientAndEnergy_Tetrahedron)
|
||
{
|
||
constexpr double TWO_PI = 2.0 * 3.14159265358979323846264338328;
|
||
constexpr double SCALE = 0.2;
|
||
|
||
auto mesh = make_tetrahedron(); // same 4 vertices as /tmp/oracle/tet.obj
|
||
auto maps = setup_spherical_maps(mesh);
|
||
|
||
// λ°_e = 2·log(SCALE · chord), matching Java's prepareInvariantData(scale).
|
||
for (auto e : mesh.edges()) {
|
||
auto h = mesh.halfedge(e);
|
||
const auto& p1 = mesh.point(mesh.source(h));
|
||
const auto& p2 = mesh.point(mesh.target(h));
|
||
const double dx = p1.x() - p2.x(), dy = p1.y() - p2.y(), dz = p1.z() - p2.z();
|
||
const double chord = std::sqrt(dx*dx + dy*dy + dz*dz);
|
||
maps.lambda0[e] = 2.0 * std::log(SCALE * chord);
|
||
}
|
||
|
||
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 = spherical_gradient(mesh, x, maps);
|
||
|
||
struct GoldRow { double X, Y, Z, G; };
|
||
const GoldRow gold[4] = {
|
||
{ 1, 1, 1, 2.7671034786104927},
|
||
{ 1, -1, -1, 2.5216834054857546},
|
||
{-1, 1, -1, 1.5401761310866633},
|
||
{-1, -1, 1, 2.6518569531861100},
|
||
};
|
||
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";
|
||
}
|
||
|
||
std::vector<double> x0(static_cast<std::size_t>(idx), 0.0);
|
||
const double dE = spherical_energy(mesh, x, maps)
|
||
- spherical_energy(mesh, x0, maps);
|
||
EXPECT_NEAR(dE, 0.16409141487397116, 1e-12);
|
||
}
|
||
|
||
// ─────────────────────────────────────────────────────────────────────────────
|
||
// FULL-MESH EDGE-DOF golden oracle (spherical, Finding 3) — the missing solution-
|
||
// level Java oracle for the edge-DOF gradient path. Makes two opposite edges of
|
||
// the shared tetrahedron variable (replacement parameterization: λ_e = x[e_idx],
|
||
// no u-shift) and pins BOTH the vertex gradient (Θ−Σα) AND the edge gradient
|
||
// (α_opp⁺ + α_opp⁻ − θ_e, θ_e = π) bit-for-bit against the REAL upstream
|
||
// SphericalFunctional.conformalEnergyAndGradient (raw gradient, not evaluate()).
|
||
//
|
||
// This is a pure GRADIENT oracle (no ΔE): with an edge DOF, x = 0 means λ_e = 0
|
||
// ⇒ spherical arc length = π (degenerate), so the path-integral-from-origin
|
||
// energy reference is ill-defined — the gradient at a fixed non-degenerate x is
|
||
// the unambiguous Finding-3 quantity. It does NOT touch the spherical Hessian,
|
||
// so it is unaffected by the Finding-4 edge-DOF Hessian guard.
|
||
//
|
||
// To regenerate: /tmp/oracle/SphereEdgeOracle.java. Values are Java %.17g.
|
||
// ─────────────────────────────────────────────────────────────────────────────
|
||
TEST(SphericalGoldenJava, FullMeshEdgeDofGradient_Tetrahedron)
|
||
{
|
||
constexpr double TWO_PI = 2.0 * 3.14159265358979323846264338328;
|
||
constexpr double PI_C = 3.14159265358979323846264338328;
|
||
constexpr double SCALE = 0.2;
|
||
|
||
auto mesh = make_tetrahedron();
|
||
auto maps = setup_spherical_maps(mesh);
|
||
|
||
for (auto e : mesh.edges()) {
|
||
auto h = mesh.halfedge(e);
|
||
const auto& p1 = mesh.point(mesh.source(h));
|
||
const auto& p2 = mesh.point(mesh.target(h));
|
||
const double dx = p1.x()-p2.x(), dy = p1.y()-p2.y(), dz = p1.z()-p2.z();
|
||
maps.lambda0[e] = 2.0 * std::log(SCALE * std::sqrt(dx*dx + dy*dy + dz*dz));
|
||
maps.theta_e[e] = PI_C;
|
||
maps.e_idx[e] = -1;
|
||
}
|
||
|
||
int idx = 0;
|
||
for (auto v : mesh.vertices()) { maps.v_idx[v] = idx++; maps.theta_v[v] = TWO_PI; }
|
||
|
||
auto pos = [&](Vertex_index v){ return mesh.point(v); };
|
||
auto same = [](const Point3& p, double X, double Y, double Z) {
|
||
return std::abs(p.x()-X)<1e-9 && std::abs(p.y()-Y)<1e-9 && std::abs(p.z()-Z)<1e-9;
|
||
};
|
||
// Edge connects positions (aX,aY,aZ)–(bX,bY,bZ) in either order?
|
||
auto connects = [&](Edge_index e, const double a[3], const double b[3]) {
|
||
auto h = mesh.halfedge(e);
|
||
const Point3& s = pos(mesh.source(h));
|
||
const Point3& t = pos(mesh.target(h));
|
||
return (same(s,a[0],a[1],a[2]) && same(t,b[0],b[1],b[2])) ||
|
||
(same(s,b[0],b[1],b[2]) && same(t,a[0],a[1],a[2]));
|
||
};
|
||
const double A[3]={1,1,1}, B[3]={1,-1,-1}, C[3]={-1,1,-1}, D[3]={-1,-1,1};
|
||
|
||
// Make edges A–B and C–D variable (replacement parameterization).
|
||
int edof = idx; // edge DOFs start after the 4 vertex DOFs
|
||
Edge_index eAB, eCD;
|
||
for (auto e : mesh.edges()) {
|
||
if (connects(e, A, B)) { maps.e_idx[e] = edof++; eAB = e; }
|
||
else if (connects(e, C, D)) { maps.e_idx[e] = edof++; eCD = e; }
|
||
}
|
||
ASSERT_EQ(edof, 6); // 4 vertex + 2 edge DOFs
|
||
|
||
auto u_of = [](const Point3& p){ return 0.10*p.x() - 0.07*p.y() + 0.13*p.z(); };
|
||
|
||
std::vector<double> x(6, 0.0);
|
||
for (auto v : mesh.vertices())
|
||
x[static_cast<std::size_t>(maps.v_idx[v])] = u_of(mesh.point(v));
|
||
// Edge DOF value = λ⁰_e + 0.1 (same perturbation as the Java oracle).
|
||
x[static_cast<std::size_t>(maps.e_idx[eAB])] = maps.lambda0[eAB] + 0.1;
|
||
x[static_cast<std::size_t>(maps.e_idx[eCD])] = maps.lambda0[eCD] + 0.1;
|
||
|
||
auto G = spherical_gradient(mesh, x, maps);
|
||
|
||
// Vertex gradients keyed by position.
|
||
struct VRow { double X, Y, Z, G; };
|
||
const VRow vg[4] = {
|
||
{ 1, 1, 1, 2.5036751008175546},
|
||
{ 1, -1, -1, 2.2303362601050205},
|
||
{-1, 1, -1, 1.8463710314817632},
|
||
{-1, -1, 1, 2.7499719487341570},
|
||
};
|
||
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 m = false;
|
||
for (auto& r : vg)
|
||
if (same(p, r.X, r.Y, r.Z)) { EXPECT_NEAR(g, r.G, 1e-12); m = true; break; }
|
||
EXPECT_TRUE(m);
|
||
}
|
||
// Edge gradients: A–B and C–D (Java golden values).
|
||
EXPECT_NEAR(G[static_cast<std::size_t>(maps.e_idx[eAB])], -0.35189517043413690, 1e-12);
|
||
EXPECT_NEAR(G[static_cast<std::size_t>(maps.e_idx[eCD])], -0.44101986058895950, 1e-12);
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Degenerate spherical triangle — limiting angles (Finding-F, java-port-audit item 1)
|
||
//
|
||
// spherical_angles() must return the limiting angles (π opposite the
|
||
// over-long edge, 0/0 elsewhere) when the spherical triangle inequality
|
||
// is violated, with valid = false. This mirrors the Euclidean behaviour
|
||
// and is required for the convex C¹ BPS extension.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalFunctional, DegenerateTriangle_LimitingAngles_S12TooLong)
|
||
{
|
||
// s12 > s23 + s31: s12 = 2.5, s23 = s31 = 0.5 (all < π so valid arc lengths)
|
||
// s23 < 0 → actually use s-based check
|
||
// Easier: use s12 = π − ε (nearly degenerate hemisphere edge)
|
||
// and very short s23, s31 so s12 > s23 + s31.
|
||
const double s12 = 2.0, s23 = 0.4, s31 = 0.4; // s12 > s23+s31 = 0.8
|
||
auto fa = spherical_angles(s12, s23, s31);
|
||
EXPECT_FALSE(fa.valid);
|
||
// s12 is the edge opposite v3 → α3 = π
|
||
EXPECT_NEAR(fa.alpha3, PI, 1e-12)
|
||
<< "corner opposite over-long s12 must be π";
|
||
EXPECT_NEAR(fa.alpha1, 0.0, 1e-12);
|
||
EXPECT_NEAR(fa.alpha2, 0.0, 1e-12);
|
||
}
|
||
|
||
TEST(SphericalFunctional, DegenerateTriangle_LimitingAngles_S23TooLong)
|
||
{
|
||
// s23 > s12 + s31 → α1 = π
|
||
const double s12 = 0.4, s23 = 2.0, s31 = 0.4;
|
||
auto fa = spherical_angles(s12, s23, s31);
|
||
EXPECT_FALSE(fa.valid);
|
||
EXPECT_NEAR(fa.alpha1, PI, 1e-12)
|
||
<< "corner opposite over-long s23 must be π";
|
||
EXPECT_NEAR(fa.alpha2, 0.0, 1e-12);
|
||
EXPECT_NEAR(fa.alpha3, 0.0, 1e-12);
|
||
}
|
||
|
||
TEST(SphericalFunctional, DegenerateTriangle_LimitingAngles_S31TooLong)
|
||
{
|
||
// s31 > s12 + s23 → α2 = π
|
||
const double s12 = 0.4, s23 = 0.4, s31 = 2.0;
|
||
auto fa = spherical_angles(s12, s23, s31);
|
||
EXPECT_FALSE(fa.valid);
|
||
EXPECT_NEAR(fa.alpha2, PI, 1e-12)
|
||
<< "corner opposite over-long s31 must be π";
|
||
EXPECT_NEAR(fa.alpha1, 0.0, 1e-12);
|
||
EXPECT_NEAR(fa.alpha3, 0.0, 1e-12);
|
||
}
|