30a132de60
C2: Fix coverage.sh — remove || true from lcov capture/extract steps so real
lcov failures are visible; empty coverage.info now exits with code 3.
C3: Add coverage gate to quality-gates CI job (SKIP_COVERAGE_GATE=1 ramp-up
mode until I5 is resolved — fast suite covers ~9.6% not 80%). Thresholds:
80% line / 70% branch / 90% function (agreed 2026-05-31).
V1: Wrap JSON parse + field extraction in try/catch — nlohmann parse_error and
type_error now surface as std::runtime_error with the file path.
V2: Wrap stoi/stod in XML Solver parser — missing/non-numeric attributes throw
std::runtime_error instead of leaking std::invalid_argument.
V4: Validate required JSON keys (dof_vector, solver, solver.*) before access —
missing field produces a clear named-field error message.
I2: 6 serialization negative tests (missing file, malformed JSON, missing
dof_vector, missing solver block, missing XML file, non-numeric XML attr).
I3: load_mesh throws on non-triangulated (quad) mesh — covers the
is_triangle_mesh guard that was previously untested.
I4: spherical_hessian throws on edge DOFs — covers the logic_error guard.
290/290 tests pass (+8 new).
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
237 lines
12 KiB
C++
237 lines
12 KiB
C++
// Copyright (c) 2024-2026 Tarik Moussa.
|
||
// SPDX-License-Identifier: MIT
|
||
|
||
// test_spherical_hessian.cpp
|
||
//
|
||
// Phase 3f — Spherical cotangent-Laplace Hessian.
|
||
//
|
||
// The Hessian of the spherical discrete conformal energy is the "spherical
|
||
// cotangent Laplacian" with edge weight
|
||
// w_k = cot(β_k), β_k = (π − αi − αj + αk) / 2
|
||
// for the edge (vi, vj) with opposite vertex vk and angles αi, αj, αk.
|
||
//
|
||
// In the flat limit (αi+αj+αk → π) this reduces to the Euclidean cotangent
|
||
// Laplacian (β_k → αk, w_k → cot(αk)).
|
||
//
|
||
// Tests:
|
||
// 1. Spherical cot weights match Euclidean weights in the near-flat limit.
|
||
// 2. Hessian is symmetric.
|
||
// 3. Hessian has H·1 ≈ 0 on the spherical tetrahedron (null-space property).
|
||
// 4. Finite-difference check (the primary correctness criterion).
|
||
|
||
#include "conformal_mesh.hpp"
|
||
#include "mesh_builder.hpp"
|
||
#include "spherical_hessian.hpp"
|
||
#include "euclidean_hessian.hpp" // for Euclidean comparison
|
||
#include <gtest/gtest.h>
|
||
#include <Eigen/Dense>
|
||
#include <cmath>
|
||
#include <vector>
|
||
|
||
using namespace conformallab;
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Spherical cot weights reduce to Euclidean cot weights in the flat limit
|
||
//
|
||
// For a nearly-flat equilateral spherical triangle (l → 0, α → 60°):
|
||
// β_k = (π − 60° − 60° + 60°)/2 = 60° → cot(60°) = 1/√3 ✓
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalHessian, CotWeights_NearFlatEquilateral)
|
||
{
|
||
// Use a very small equilateral spherical triangle: α1=α2=α3=60°
|
||
const double alpha = PI / 3.0;
|
||
auto sw = spherical_cot_weights(alpha, alpha, alpha);
|
||
|
||
ASSERT_TRUE(sw.valid);
|
||
|
||
const double expected = 1.0 / std::sqrt(3.0); // = cot(60°)
|
||
EXPECT_NEAR(sw.w12, expected, 1e-12);
|
||
EXPECT_NEAR(sw.w23, expected, 1e-12);
|
||
EXPECT_NEAR(sw.w31, expected, 1e-12);
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Spherical cot weights are positive for the regular spherical tetrahedron
|
||
//
|
||
// Each face has α_k = 2π/3 (120°), angle sum = 2π.
|
||
// β_k = (π − 2π/3 − 2π/3 + 2π/3)/2 = (π − 2π/3)/2 = π/6
|
||
// w_k = cot(π/6) = √3
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalHessian, CotWeights_RegularSphericalTetrahedronFace)
|
||
{
|
||
const double alpha = 2.0 * PI / 3.0; // 120°
|
||
auto sw = spherical_cot_weights(alpha, alpha, alpha);
|
||
|
||
ASSERT_TRUE(sw.valid);
|
||
|
||
const double expected = std::sqrt(3.0); // cot(π/6)
|
||
EXPECT_NEAR(sw.w12, expected, 1e-10);
|
||
EXPECT_NEAR(sw.w23, expected, 1e-10);
|
||
EXPECT_NEAR(sw.w31, expected, 1e-10);
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Hessian is symmetric: H[i,j] == H[j,i]
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalHessian, HessianIsSymmetric)
|
||
{
|
||
auto mesh = make_spherical_tetrahedron();
|
||
auto maps = setup_spherical_maps(mesh);
|
||
compute_lambda0_from_mesh(mesh, maps);
|
||
int n = assign_vertex_dof_indices(mesh, maps);
|
||
|
||
std::vector<double> x(static_cast<std::size_t>(n), -0.2);
|
||
auto H = spherical_hessian(mesh, x, maps);
|
||
|
||
Eigen::MatrixXd Hd = Eigen::MatrixXd(H);
|
||
EXPECT_NEAR((Hd - Hd.transpose()).norm(), 0.0, 1e-12)
|
||
<< "Spherical Hessian must be symmetric";
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// H·1 is NOT zero for the spherical Hessian (unlike the Euclidean case).
|
||
//
|
||
// In the Euclidean case, a uniform shift u_i → u_i + c scales all edge
|
||
// lengths by e^c, leaving angles unchanged → H·1 = 0 exactly.
|
||
//
|
||
// In the spherical case, l_ij = 2·asin(exp(λ_ij/2)) is NOT a linear
|
||
// function of the DOFs, so a uniform shift DOES change the angles
|
||
// → H·1 ≠ 0 in general. This test verifies that property.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalHessian, ConstantVectorNotInNullSpace)
|
||
{
|
||
auto mesh = make_spherical_tetrahedron();
|
||
auto maps = setup_spherical_maps(mesh);
|
||
compute_lambda0_from_mesh(mesh, maps);
|
||
int n = assign_vertex_dof_indices(mesh, maps);
|
||
|
||
std::vector<double> x(static_cast<std::size_t>(n), -0.2);
|
||
auto H = spherical_hessian(mesh, x, maps);
|
||
|
||
Eigen::VectorXd ones = Eigen::VectorXd::Ones(n);
|
||
Eigen::VectorXd Hones = H * ones;
|
||
|
||
// H·1 should be non-trivial (norm well above zero)
|
||
EXPECT_GT(Hones.norm(), 1e-3)
|
||
<< "Spherical H·1 should be non-zero; norm = " << Hones.norm();
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Hessian is negative semi-definite at the equilibrium point (x = 0)
|
||
//
|
||
// The spherical discrete conformal energy is concave in the vertex DOFs
|
||
// (unlike the Euclidean case which is convex). At the equilibrium x = 0,
|
||
// the Hessian is NSD: all eigenvalues ≤ 0, with a multi-dimensional null
|
||
// space corresponding to degenerate directions.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalHessian, HessianIsNegativeSemiDefiniteAtEquilibrium)
|
||
{
|
||
auto mesh = make_spherical_tetrahedron();
|
||
auto maps = setup_spherical_maps(mesh);
|
||
compute_lambda0_from_mesh(mesh, maps);
|
||
int n = assign_vertex_dof_indices(mesh, maps);
|
||
|
||
// x = 0 is the equilibrium for the regular spherical tetrahedron.
|
||
std::vector<double> x(static_cast<std::size_t>(n), 0.0);
|
||
auto H = spherical_hessian(mesh, x, maps);
|
||
|
||
Eigen::MatrixXd Hd = Eigen::MatrixXd(H);
|
||
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(Hd);
|
||
double max_ev = es.eigenvalues().maxCoeff();
|
||
|
||
EXPECT_LE(max_ev, 1e-9)
|
||
<< "Largest eigenvalue of spherical Hessian at equilibrium must be ≤ 0; got " << max_ev;
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Finite-difference Hessian check: octahedron-face triangle (vertex DOFs)
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalHessian, FDCheck_OctaFaceVertex)
|
||
{
|
||
auto mesh = make_octahedron_face();
|
||
auto maps = setup_spherical_maps(mesh);
|
||
compute_lambda0_from_mesh(mesh, maps);
|
||
int n = assign_vertex_dof_indices(mesh, maps);
|
||
|
||
std::vector<double> x(static_cast<std::size_t>(n), -0.3);
|
||
|
||
EXPECT_TRUE(hessian_check_spherical(mesh, x, maps))
|
||
<< "FD Hessian check failed on octahedron-face triangle";
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Finite-difference Hessian check: spherical tetrahedron (vertex DOFs)
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalHessian, FDCheck_SpherTetVertex)
|
||
{
|
||
auto mesh = make_spherical_tetrahedron();
|
||
auto maps = setup_spherical_maps(mesh);
|
||
compute_lambda0_from_mesh(mesh, maps);
|
||
int n = assign_vertex_dof_indices(mesh, maps);
|
||
|
||
std::vector<double> x(static_cast<std::size_t>(n), -0.2);
|
||
|
||
EXPECT_TRUE(hessian_check_spherical(mesh, x, maps))
|
||
<< "FD Hessian check failed on spherical tetrahedron";
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Finite-difference Hessian check: mixed pinned/variable vertices
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalHessian, FDCheck_MixedPinnedVertices)
|
||
{
|
||
auto mesh = make_octahedron_face();
|
||
auto maps = setup_spherical_maps(mesh);
|
||
compute_lambda0_from_mesh(mesh, maps);
|
||
|
||
auto vit = mesh.vertices().begin();
|
||
Vertex_index v0 = *vit++;
|
||
Vertex_index v1 = *vit++;
|
||
Vertex_index v2 = *vit;
|
||
|
||
maps.v_idx[v0] = -1; // pinned
|
||
maps.v_idx[v1] = 0;
|
||
maps.v_idx[v2] = 1;
|
||
|
||
std::vector<double> x = {-0.2, -0.3};
|
||
|
||
EXPECT_TRUE(hessian_check_spherical(mesh, x, maps))
|
||
<< "FD Hessian check failed for mixed pinned/variable vertices";
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// I4: spherical_hessian throws on edge DOFs
|
||
//
|
||
// The spherical Hessian only implements the vertex-block cotangent Laplacian;
|
||
// if any edge has a free DOF index (e_idx >= 0) it must fail loudly rather
|
||
// than returning a silently rank-deficient matrix.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SphericalHessian, ThrowsOnEdgeDOF)
|
||
{
|
||
// Build a tetrahedron and assign one edge as a free DOF.
|
||
auto mesh = make_tetrahedron();
|
||
auto maps = setup_spherical_maps(mesh);
|
||
compute_lambda0_from_mesh(mesh, maps); // spherical variant (A2 rename pending)
|
||
|
||
int idx = 0;
|
||
for (auto v : mesh.vertices())
|
||
maps.v_idx[v] = idx++;
|
||
const int n_v = idx;
|
||
|
||
// Give one edge a free DOF index to trigger the guard.
|
||
auto e0 = *mesh.edges().begin();
|
||
maps.e_idx[e0] = n_v; // first free edge DOF
|
||
|
||
std::vector<double> x(static_cast<std::size_t>(n_v + 1), 0.0);
|
||
EXPECT_THROW(spherical_hessian(mesh, x, maps), std::logic_error);
|
||
}
|