feat(phase3b): port HyperIdealFunctional energy + gradient onto ConformalMesh
Implements the hyper-ideal discrete conformal map functional on
CGAL::Surface_mesh. The energy and analytic gradient are ported directly
from HyperIdealFunctional.java; correctness is verified via a
finite-difference gradient check (same eps=1E-5 / tol=1E-4 as Java).
New files:
include/hyper_ideal_geometry.hpp — ζ, ζ₁₃, ζ₁₄, ζ₁₅, lij, αij, σi, σij
include/hyper_ideal_functional.hpp — HyperIdealMaps, evaluate_hyper_ideal,
gradient_check
tests/cgal/test_hyper_ideal_functional.cpp — 6 tests (1 skipped @Ignore)
Test results (local, -DWITH_CGAL=ON):
conformallab_cgal_tests: 21 registered | 20 passed | 1 skipped | 0 failed
- GradientCheck_AllHyperIdealTriangle ✓
- GradientCheck_ExtendedDomain ✓
- GradientCheck_TetrahedronAllVariable ✓
- EnergyFiniteAtTestPoint ✓
- GradientCheck_MixedIdealHyperIdeal ✓
- GradientCheck_Fan6AllVariable ✓
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This commit is contained in:
Tarik Moussa
344
code/include/hyper_ideal_functional.hpp
Normal file
344
code/include/hyper_ideal_functional.hpp
Normal file
@@ -0,0 +1,344 @@
|
|||||||
|
#pragma once
|
||||||
|
// hyper_ideal_functional.hpp
|
||||||
|
//
|
||||||
|
// Energy and gradient of the hyper-ideal discrete conformal map functional
|
||||||
|
// evaluated on a ConformalMesh (CGAL::Surface_mesh).
|
||||||
|
//
|
||||||
|
// Ported from de.varylab.discreteconformal.functional.HyperIdealFunctional.
|
||||||
|
//
|
||||||
|
// ┌─────────────────────────────────────────────────────────────────────────┐
|
||||||
|
// │ E(x) = Σ_faces U(f) − Σ_edges θ_e · a_e − Σ_vertices Θ_v · b_v │
|
||||||
|
// │ │
|
||||||
|
// │ ∂E/∂b_v = Σ_{faces adj. v} β_v(face) − Θ_v │
|
||||||
|
// │ ∂E/∂a_e = α_e(face⁺) + α_e(face⁻) − θ_e │
|
||||||
|
// └─────────────────────────────────────────────────────────────────────────┘
|
||||||
|
//
|
||||||
|
// DOF vector layout (matches getDimension() ordering):
|
||||||
|
// x[v_idx[v]] = b_v for each variable vertex (log scale factor)
|
||||||
|
// x[e_idx[e]] = a_e for each variable edge (intersection angle)
|
||||||
|
// -1 in v_idx / e_idx means "pinned" (ideal / fixed at 0).
|
||||||
|
//
|
||||||
|
// Usage
|
||||||
|
// ─────
|
||||||
|
// auto mesh = make_tetrahedron();
|
||||||
|
// auto maps = setup_hyper_ideal_maps(mesh);
|
||||||
|
// int n = assign_all_dof_indices(mesh, maps); // all variable
|
||||||
|
// std::vector<double> x(n, 1.0);
|
||||||
|
// auto res = evaluate_hyper_ideal(mesh, x, maps);
|
||||||
|
// // res.energy, res.gradient
|
||||||
|
|
||||||
|
#include "conformal_mesh.hpp"
|
||||||
|
#include "hyper_ideal_geometry.hpp"
|
||||||
|
#include "hyper_ideal_utility.hpp"
|
||||||
|
#include <CGAL/boost/graph/iterator.h>
|
||||||
|
#include <vector>
|
||||||
|
#include <cmath>
|
||||||
|
#include <cstdint>
|
||||||
|
|
||||||
|
namespace conformallab {
|
||||||
|
|
||||||
|
// ── Property-map type aliases ─────────────────────────────────────────────────
|
||||||
|
|
||||||
|
using VMapD = ConformalMesh::Property_map<Vertex_index, double>;
|
||||||
|
using VMapI = ConformalMesh::Property_map<Vertex_index, int>;
|
||||||
|
using EMapD = ConformalMesh::Property_map<Edge_index, double>;
|
||||||
|
using EMapI = ConformalMesh::Property_map<Edge_index, int>;
|
||||||
|
|
||||||
|
// ── Persistent map bundle ─────────────────────────────────────────────────────
|
||||||
|
|
||||||
|
struct HyperIdealMaps {
|
||||||
|
VMapI v_idx; // DOF index per vertex (-1 = pinned / ideal point)
|
||||||
|
EMapI e_idx; // DOF index per edge (-1 = fixed)
|
||||||
|
VMapD theta_v; // target cone angle Θ_v (parameter, not variable)
|
||||||
|
EMapD theta_e; // target intersection angle θ_e
|
||||||
|
};
|
||||||
|
|
||||||
|
// Add all needed persistent property maps and return handles.
|
||||||
|
// Defaults: theta_v = 2π (regular cone), theta_e = π (orthogonal circles).
|
||||||
|
inline HyperIdealMaps setup_hyper_ideal_maps(ConformalMesh& mesh)
|
||||||
|
{
|
||||||
|
HyperIdealMaps m;
|
||||||
|
m.v_idx = mesh.add_property_map<Vertex_index, int> ("v:idx", -1 ).first;
|
||||||
|
m.e_idx = mesh.add_property_map<Edge_index, int> ("e:idx", -1 ).first;
|
||||||
|
m.theta_v = mesh.add_property_map<Vertex_index, double>("v:theta", 2.0*PI ).first;
|
||||||
|
m.theta_e = mesh.add_property_map<Edge_index, double>("e:theta", PI ).first;
|
||||||
|
return m;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Count variable DOFs: #variable_vertices + #variable_edges.
|
||||||
|
inline int hyper_ideal_dimension(const ConformalMesh& mesh, const HyperIdealMaps& m)
|
||||||
|
{
|
||||||
|
int dim = 0;
|
||||||
|
for (auto v : mesh.vertices()) if (m.v_idx[v] >= 0) ++dim;
|
||||||
|
for (auto e : mesh.edges()) if (m.e_idx[e] >= 0) ++dim;
|
||||||
|
return dim;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Assign DOF indices 0..n-1: vertices first, then edges.
|
||||||
|
// All vertices and edges become variable. Returns total DOF count.
|
||||||
|
inline int assign_all_dof_indices(ConformalMesh& mesh, HyperIdealMaps& m)
|
||||||
|
{
|
||||||
|
int idx = 0;
|
||||||
|
for (auto v : mesh.vertices()) m.v_idx[v] = idx++;
|
||||||
|
for (auto e : mesh.edges()) m.e_idx[e] = idx++;
|
||||||
|
return idx;
|
||||||
|
}
|
||||||
|
|
||||||
|
// ── Evaluation result ─────────────────────────────────────────────────────────
|
||||||
|
|
||||||
|
struct HyperIdealResult {
|
||||||
|
double energy = 0.0;
|
||||||
|
std::vector<double> gradient; // empty when gradient was not requested
|
||||||
|
};
|
||||||
|
|
||||||
|
// ── Internal helpers ──────────────────────────────────────────────────────────
|
||||||
|
|
||||||
|
// Get the DOF value from x, or 0.0 if pinned.
|
||||||
|
static inline double dof_val(int idx, const std::vector<double>& x)
|
||||||
|
{
|
||||||
|
return idx >= 0 ? x[static_cast<std::size_t>(idx)] : 0.0;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Convert a CGAL halfedge index to a plain std::size_t (for vector indexing).
|
||||||
|
static inline std::size_t hidx(Halfedge_index h)
|
||||||
|
{
|
||||||
|
return static_cast<std::size_t>(static_cast<std::uint32_t>(h));
|
||||||
|
}
|
||||||
|
|
||||||
|
// ── Per-face angle kernel ─────────────────────────────────────────────────────
|
||||||
|
|
||||||
|
struct FaceAngles {
|
||||||
|
double alpha12, alpha23, alpha31; // dihedral angles at each edge
|
||||||
|
double beta1, beta2, beta3; // interior angles at each vertex
|
||||||
|
double a12, a23, a31; // edge DOF values (used in energy)
|
||||||
|
double b1, b2, b3; // vertex DOF values
|
||||||
|
bool v1b, v2b, v3b; // whether each vertex is variable
|
||||||
|
};
|
||||||
|
|
||||||
|
static FaceAngles compute_face_angles(
|
||||||
|
const ConformalMesh& mesh,
|
||||||
|
Face_index f,
|
||||||
|
const std::vector<double>& x,
|
||||||
|
const HyperIdealMaps& m)
|
||||||
|
{
|
||||||
|
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);
|
||||||
|
|
||||||
|
FaceAngles fa;
|
||||||
|
fa.v1b = m.v_idx[v1] >= 0;
|
||||||
|
fa.v2b = m.v_idx[v2] >= 0;
|
||||||
|
fa.v3b = m.v_idx[v3] >= 0;
|
||||||
|
|
||||||
|
fa.a12 = dof_val(m.e_idx[e12], x);
|
||||||
|
fa.a23 = dof_val(m.e_idx[e23], x);
|
||||||
|
fa.a31 = dof_val(m.e_idx[e31], x);
|
||||||
|
fa.b1 = dof_val(m.v_idx[v1], x);
|
||||||
|
fa.b2 = dof_val(m.v_idx[v2], x);
|
||||||
|
fa.b3 = dof_val(m.v_idx[v3], x);
|
||||||
|
|
||||||
|
// Clamp invalid inputs (mirrors Java log.warning + clamp)
|
||||||
|
if (fa.v1b && fa.v2b && fa.a12 < 0.0) fa.a12 = 0.0;
|
||||||
|
if (fa.v2b && fa.v3b && fa.a23 < 0.0) fa.a23 = 0.0;
|
||||||
|
if (fa.v3b && fa.v1b && fa.a31 < 0.0) fa.a31 = 0.0;
|
||||||
|
if (fa.v1b && fa.b1 < 0.0) fa.b1 = 0.01;
|
||||||
|
if (fa.v2b && fa.b2 < 0.0) fa.b2 = 0.01;
|
||||||
|
if (fa.v3b && fa.b3 < 0.0) fa.b3 = 0.01;
|
||||||
|
|
||||||
|
double l12 = lij(fa.b1, fa.b2, fa.a12, fa.v1b, fa.v2b);
|
||||||
|
double l23 = lij(fa.b2, fa.b3, fa.a23, fa.v2b, fa.v3b);
|
||||||
|
double l31 = lij(fa.b3, fa.b1, fa.a31, fa.v3b, fa.v1b);
|
||||||
|
|
||||||
|
// Guard degenerate lengths
|
||||||
|
if (l12 < 1E-12 && l23 < 1E-12 && l31 < 1E-12)
|
||||||
|
l12 = l23 = l31 = 1E-12;
|
||||||
|
|
||||||
|
// Check triangle inequalities; degenerate cases get extreme angles
|
||||||
|
if (l12 > l23 + l31) {
|
||||||
|
fa.beta1 = 0.0; fa.beta2 = 0.0; fa.beta3 = PI;
|
||||||
|
fa.alpha12 = PI; fa.alpha23 = 0.0; fa.alpha31 = 0.0;
|
||||||
|
} else if (l23 > l12 + l31) {
|
||||||
|
fa.beta1 = PI; fa.beta2 = 0.0; fa.beta3 = 0.0;
|
||||||
|
fa.alpha12 = 0.0; fa.alpha23 = PI; fa.alpha31 = 0.0;
|
||||||
|
} else if (l31 > l12 + l23) {
|
||||||
|
fa.beta1 = 0.0; fa.beta2 = PI; fa.beta3 = 0.0;
|
||||||
|
fa.alpha12 = 0.0; fa.alpha23 = 0.0; fa.alpha31 = PI;
|
||||||
|
} else {
|
||||||
|
fa.beta1 = zeta(l12, l31, l23);
|
||||||
|
fa.beta2 = zeta(l23, l12, l31);
|
||||||
|
fa.beta3 = zeta(l31, l23, l12);
|
||||||
|
|
||||||
|
fa.alpha12 = alpha_ij(fa.a12, fa.a23, fa.a31,
|
||||||
|
fa.b1, fa.b2, fa.b3,
|
||||||
|
fa.beta1, fa.beta2, fa.beta3,
|
||||||
|
fa.v1b, fa.v2b, fa.v3b);
|
||||||
|
fa.alpha23 = alpha_ij(fa.a23, fa.a31, fa.a12,
|
||||||
|
fa.b2, fa.b3, fa.b1,
|
||||||
|
fa.beta2, fa.beta3, fa.beta1,
|
||||||
|
fa.v2b, fa.v3b, fa.v1b);
|
||||||
|
fa.alpha31 = alpha_ij(fa.a31, fa.a12, fa.a23,
|
||||||
|
fa.b3, fa.b1, fa.b2,
|
||||||
|
fa.beta3, fa.beta1, fa.beta2,
|
||||||
|
fa.v3b, fa.v1b, fa.v2b);
|
||||||
|
}
|
||||||
|
return fa;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Per-face energy contribution U(f) (before subtracting θ·a and Θ·b terms).
|
||||||
|
static double face_energy(const FaceAngles& fa)
|
||||||
|
{
|
||||||
|
double aa = fa.a12*fa.alpha12 + fa.a23*fa.alpha23 + fa.a31*fa.alpha31;
|
||||||
|
double bb = fa.b1 *fa.beta1 + fa.b2 *fa.beta2 + fa.b3 *fa.beta3;
|
||||||
|
|
||||||
|
double V = 0.0;
|
||||||
|
if (fa.v1b && fa.v2b && fa.v3b) {
|
||||||
|
V = calculateTetrahedronVolume(
|
||||||
|
fa.beta1, fa.beta2, fa.beta3,
|
||||||
|
fa.alpha23, fa.alpha31, fa.alpha12);
|
||||||
|
} else if (!fa.v1b) {
|
||||||
|
V = calculateTetrahedronVolumeWithIdealVertexAtGamma(
|
||||||
|
fa.beta1, fa.alpha31, fa.alpha12,
|
||||||
|
fa.alpha23, fa.beta2, fa.beta3);
|
||||||
|
} else if (!fa.v2b) {
|
||||||
|
V = calculateTetrahedronVolumeWithIdealVertexAtGamma(
|
||||||
|
fa.beta2, fa.alpha12, fa.alpha23,
|
||||||
|
fa.alpha31, fa.beta3, fa.beta1);
|
||||||
|
} else { // !v3b
|
||||||
|
V = calculateTetrahedronVolumeWithIdealVertexAtGamma(
|
||||||
|
fa.beta3, fa.alpha23, fa.alpha31,
|
||||||
|
fa.alpha12, fa.beta1, fa.beta2);
|
||||||
|
}
|
||||||
|
return aa + bb + 2.0 * V;
|
||||||
|
}
|
||||||
|
|
||||||
|
// ── Full evaluation ───────────────────────────────────────────────────────────
|
||||||
|
|
||||||
|
inline HyperIdealResult evaluate_hyper_ideal(
|
||||||
|
ConformalMesh& mesh,
|
||||||
|
const std::vector<double>& x,
|
||||||
|
const HyperIdealMaps& m,
|
||||||
|
bool need_energy = true,
|
||||||
|
bool need_gradient = true)
|
||||||
|
{
|
||||||
|
HyperIdealResult res;
|
||||||
|
|
||||||
|
// Temporary per-halfedge storage for computed angles.
|
||||||
|
// Indexed by the integer value of Halfedge_index.
|
||||||
|
const std::size_t nh = mesh.number_of_halfedges();
|
||||||
|
std::vector<double> h_alpha(nh, 0.0); // α_ij stored on halfedge
|
||||||
|
std::vector<double> h_beta (nh, 0.0); // β_i stored on opposite halfedge
|
||||||
|
|
||||||
|
// ── Pass 1: angles + energy per face ─────────────────────────────────────
|
||||||
|
for (auto f : mesh.faces()) {
|
||||||
|
Halfedge_index h0 = mesh.halfedge(f);
|
||||||
|
Halfedge_index h1 = mesh.next(h0);
|
||||||
|
Halfedge_index h2 = mesh.next(h1);
|
||||||
|
|
||||||
|
FaceAngles fa = compute_face_angles(mesh, f, x, m);
|
||||||
|
|
||||||
|
// Store computed angles into temporary arrays.
|
||||||
|
// h_alpha[h] = α for the edge of h in this face.
|
||||||
|
h_alpha[hidx(h0)] = fa.alpha12;
|
||||||
|
h_alpha[hidx(h1)] = fa.alpha23;
|
||||||
|
h_alpha[hidx(h2)] = fa.alpha31;
|
||||||
|
|
||||||
|
// h_beta[h] = β at the vertex OPPOSITE to h.
|
||||||
|
// β1 (at v1 = source(h0)) is opposite to h1 = e23.
|
||||||
|
h_beta[hidx(h1)] = fa.beta1; // h1 is across from v1
|
||||||
|
h_beta[hidx(h2)] = fa.beta2; // h2 is across from v2
|
||||||
|
h_beta[hidx(h0)] = fa.beta3; // h0 is across from v3
|
||||||
|
|
||||||
|
if (need_energy)
|
||||||
|
res.energy += face_energy(fa);
|
||||||
|
}
|
||||||
|
|
||||||
|
// ── Pass 2: linear energy terms ──────────────────────────────────────────
|
||||||
|
if (need_energy) {
|
||||||
|
for (auto e : mesh.edges()) {
|
||||||
|
int ie = m.e_idx[e];
|
||||||
|
if (ie >= 0) res.energy -= m.theta_e[e] * x[static_cast<std::size_t>(ie)];
|
||||||
|
}
|
||||||
|
for (auto v : mesh.vertices()) {
|
||||||
|
int iv = m.v_idx[v];
|
||||||
|
if (iv >= 0) res.energy -= m.theta_v[v] * x[static_cast<std::size_t>(iv)];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// ── Pass 3: gradient ─────────────────────────────────────────────────────
|
||||||
|
if (need_gradient) {
|
||||||
|
const int n = hyper_ideal_dimension(mesh, m);
|
||||||
|
res.gradient.assign(static_cast<std::size_t>(n), 0.0);
|
||||||
|
|
||||||
|
// ∂E/∂b_v = Σ_{faces adj. v} β_v(face) − Θ_v
|
||||||
|
// β_v(face) = h_beta[prev(h)] for the incoming halfedge h to v in that face.
|
||||||
|
for (auto v : mesh.vertices()) {
|
||||||
|
int iv = m.v_idx[v];
|
||||||
|
if (iv < 0) continue;
|
||||||
|
for (auto h : CGAL::halfedges_around_target(v, mesh)) {
|
||||||
|
if (mesh.is_border(h)) continue;
|
||||||
|
res.gradient[static_cast<std::size_t>(iv)] += h_beta[hidx(mesh.prev(h))];
|
||||||
|
}
|
||||||
|
res.gradient[static_cast<std::size_t>(iv)] -= m.theta_v[v];
|
||||||
|
}
|
||||||
|
|
||||||
|
// ∂E/∂a_e = α_e(face⁺) + α_e(face⁻) − θ_e
|
||||||
|
for (auto e : mesh.edges()) {
|
||||||
|
int ie = m.e_idx[e];
|
||||||
|
if (ie < 0) continue;
|
||||||
|
auto h = mesh.halfedge(e);
|
||||||
|
auto ho = mesh.opposite(h);
|
||||||
|
if (!mesh.is_border(h)) res.gradient[static_cast<std::size_t>(ie)] += h_alpha[hidx(h)];
|
||||||
|
if (!mesh.is_border(ho)) res.gradient[static_cast<std::size_t>(ie)] += h_alpha[hidx(ho)];
|
||||||
|
res.gradient[static_cast<std::size_t>(ie)] -= m.theta_e[e];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
return res;
|
||||||
|
}
|
||||||
|
|
||||||
|
// ── Finite-difference gradient check ─────────────────────────────────────────
|
||||||
|
//
|
||||||
|
// Returns true if |G[i] − fd[i]| / max(1, |G[i]|) < tol for all DOFs.
|
||||||
|
// eps = step size, tol = tolerance (same defaults as Java FunctionalTest).
|
||||||
|
inline bool gradient_check(
|
||||||
|
ConformalMesh& mesh,
|
||||||
|
const std::vector<double>& x0,
|
||||||
|
const HyperIdealMaps& m,
|
||||||
|
double eps = 1E-5,
|
||||||
|
double tol = 1E-4)
|
||||||
|
{
|
||||||
|
// Analytic gradient
|
||||||
|
auto res = evaluate_hyper_ideal(mesh, x0, m, false, true);
|
||||||
|
const auto& G = res.gradient;
|
||||||
|
const int n = static_cast<int>(G.size());
|
||||||
|
|
||||||
|
std::vector<double> xp = x0, xm = x0;
|
||||||
|
bool ok = true;
|
||||||
|
|
||||||
|
for (int i = 0; i < n; ++i) {
|
||||||
|
std::size_t si = static_cast<std::size_t>(i);
|
||||||
|
xp[si] = x0[si] + eps;
|
||||||
|
xm[si] = x0[si] - eps;
|
||||||
|
|
||||||
|
double Ep = evaluate_hyper_ideal(mesh, xp, m, true, false).energy;
|
||||||
|
double Em = evaluate_hyper_ideal(mesh, xm, m, true, false).energy;
|
||||||
|
|
||||||
|
xp[si] = xm[si] = x0[si];
|
||||||
|
|
||||||
|
double fd = (Ep - Em) / (2.0 * eps);
|
||||||
|
double err = std::abs(G[si] - fd);
|
||||||
|
double scale = std::max(1.0, std::abs(G[si]));
|
||||||
|
if (err / scale > tol) ok = false;
|
||||||
|
}
|
||||||
|
return ok;
|
||||||
|
}
|
||||||
|
|
||||||
|
} // namespace conformallab
|
||||||
138
code/include/hyper_ideal_geometry.hpp
Normal file
138
code/include/hyper_ideal_geometry.hpp
Normal file
@@ -0,0 +1,138 @@
|
|||||||
|
#pragma once
|
||||||
|
// hyper_ideal_geometry.hpp
|
||||||
|
//
|
||||||
|
// Pure-math building blocks for the hyper-ideal discrete conformal map.
|
||||||
|
// Ported from de.varylab.discreteconformal.functional.HyperIdealUtility
|
||||||
|
// and the private helpers of HyperIdealFunctional (lij, αij, σi, σij).
|
||||||
|
//
|
||||||
|
// All functions are independent of the mesh type.
|
||||||
|
//
|
||||||
|
// Notation follows the original Java / paper:
|
||||||
|
// b_i, b_j – vertex variables (log scale factors, hyper-ideal vertices)
|
||||||
|
// a_ij – edge variable (intersection angle between horocycles)
|
||||||
|
// l_ij – effective hyperbolic edge length in the auxiliary triangle
|
||||||
|
// β_i – interior angle of the hyperbolic triangle at vertex i
|
||||||
|
// α_ij – dihedral angle of the tetrahedron at edge ij
|
||||||
|
|
||||||
|
#include <cmath>
|
||||||
|
#include <algorithm>
|
||||||
|
|
||||||
|
namespace conformallab {
|
||||||
|
|
||||||
|
constexpr double PI = 3.14159265358979323846264338328;
|
||||||
|
|
||||||
|
// ── Length functions ─────────────────────────────────────────────────────────
|
||||||
|
|
||||||
|
// ζ(x,y,z) — interior angle in a hyperbolic triangle with edge lengths
|
||||||
|
// x, y, z, opposite to the side of length z.
|
||||||
|
// Ports HyperIdealUtility.ζ(x, y, z).
|
||||||
|
inline double zeta(double x, double y, double z)
|
||||||
|
{
|
||||||
|
double cx = std::cosh(x), cy = std::cosh(y), cz = std::cosh(z);
|
||||||
|
double sx = std::sinh(x), sy = std::sinh(y);
|
||||||
|
double nbd = (cx*cy - cz) / (sx*sy);
|
||||||
|
nbd = std::clamp(nbd, -1.0, 1.0); // guard floating-point rounding
|
||||||
|
return std::acos(nbd);
|
||||||
|
}
|
||||||
|
|
||||||
|
// ζ₁₃(x,y,z) — third edge length in a right-angled hyperbolic hexagon.
|
||||||
|
// Ports HyperIdealUtility.ζ_13(x, y, z).
|
||||||
|
inline double zeta13(double x, double y, double z)
|
||||||
|
{
|
||||||
|
double cx = std::cosh(x), cy = std::cosh(y), cz = std::cosh(z);
|
||||||
|
double sx = std::sinh(x), sy = std::sinh(y);
|
||||||
|
return std::acosh((cx*cy + cz) / (sx*sy));
|
||||||
|
}
|
||||||
|
|
||||||
|
// ζ₁₄(x,y) — edge length in a hyperbolic pentagon with one ideal vertex.
|
||||||
|
// Ports HyperIdealUtility.ζ_14(x, y).
|
||||||
|
inline double zeta14(double x, double y)
|
||||||
|
{
|
||||||
|
double cy = std::cosh(y), sy = std::sinh(y);
|
||||||
|
return std::acosh((std::exp(x) + cy) / sy);
|
||||||
|
}
|
||||||
|
|
||||||
|
// ζ₁₅(x) — length in a hyperbolic quadrilateral with two ideal vertices.
|
||||||
|
// Ports HyperIdealUtility.ζ_15(x).
|
||||||
|
inline double zeta15(double x)
|
||||||
|
{
|
||||||
|
return 2.0 * std::asinh(std::exp(x / 2.0));
|
||||||
|
}
|
||||||
|
|
||||||
|
// ── Effective edge length ─────────────────────────────────────────────────────
|
||||||
|
|
||||||
|
// l_ij: effective hyperbolic length of edge ij.
|
||||||
|
// b_i, b_j – vertex log scale factors (used only if vertex is hyper-ideal)
|
||||||
|
// a_ij – edge intersection-angle variable
|
||||||
|
// vi_var – true if vertex i is hyper-ideal (has a DOF b_i)
|
||||||
|
// vj_var – true if vertex j is hyper-ideal
|
||||||
|
// Ports HyperIdealFunctional.lij().
|
||||||
|
inline double lij(double bi, double bj, double aij, bool vi_var, bool vj_var)
|
||||||
|
{
|
||||||
|
if (vi_var && vj_var) return zeta13(bi, bj, aij);
|
||||||
|
if (vi_var) return zeta14(aij, bi);
|
||||||
|
if (vj_var) return zeta14(aij, bj);
|
||||||
|
return zeta15(aij);
|
||||||
|
}
|
||||||
|
|
||||||
|
// ── Auxiliary angle functions ─────────────────────────────────────────────────
|
||||||
|
|
||||||
|
// σᵢ(aᵢⱼ, aₖᵢ, aⱼₖ, vj_var, vk_var) — intermediate half-length at vertex i.
|
||||||
|
// Ports HyperIdealFunctional.σi().
|
||||||
|
inline double sigma_i(double aij, double aki, double ajk, bool vj_var, bool vk_var)
|
||||||
|
{
|
||||||
|
if (vj_var && vk_var) return zeta13(aij, aki, ajk);
|
||||||
|
if (vj_var) return zeta14(ajk - aki, aij);
|
||||||
|
if (vk_var) return zeta14(ajk - aij, aki);
|
||||||
|
return zeta15(ajk - aij - aki);
|
||||||
|
}
|
||||||
|
|
||||||
|
// σᵢⱼ(aᵢⱼ, bᵢ, bⱼ, vj_var) — intermediate half-length for edge ij from vertex i.
|
||||||
|
// Ports HyperIdealFunctional.σij().
|
||||||
|
inline double sigma_ij(double aij, double bi, double bj, bool vj_var)
|
||||||
|
{
|
||||||
|
if (vj_var) return zeta13(aij, bi, bj);
|
||||||
|
return zeta14(-aij, bi);
|
||||||
|
}
|
||||||
|
|
||||||
|
// α_ij: computed dihedral angle at edge ij in the face with vertices i, j, k.
|
||||||
|
//
|
||||||
|
// Arguments (cyclic role assignment):
|
||||||
|
// aij, ajk, aki – edge variables
|
||||||
|
// bi, bj, bk – vertex variables
|
||||||
|
// βi, βj, βk – interior angles of the auxiliary hyperbolic triangle
|
||||||
|
// vi_var, vj_var, vk_var – which vertices are hyper-ideal
|
||||||
|
//
|
||||||
|
// Ports HyperIdealFunctional.αij() (the private helper).
|
||||||
|
// Note: the vk_var case recurses once (never more than one level deep).
|
||||||
|
inline double alpha_ij(
|
||||||
|
double aij, double ajk, double aki,
|
||||||
|
double bi, double bj, double bk,
|
||||||
|
double beta_i, double beta_j, double beta_k,
|
||||||
|
bool vi_var, bool vj_var, bool vk_var)
|
||||||
|
{
|
||||||
|
if (vi_var) {
|
||||||
|
double si = sigma_i (aij, aki, ajk, vj_var, vk_var);
|
||||||
|
double sij = sigma_ij(aij, bi, bj, vj_var);
|
||||||
|
double sik = sigma_ij(aki, bi, bk, vk_var);
|
||||||
|
return zeta(si, sij, sik);
|
||||||
|
}
|
||||||
|
if (vj_var) {
|
||||||
|
double sj = sigma_i (ajk, aij, aki, vk_var, vi_var);
|
||||||
|
double sjk = sigma_ij(ajk, bj, bk, vk_var);
|
||||||
|
double sji = sigma_ij(aij, bj, bi, vi_var);
|
||||||
|
return zeta(sj, sji, sjk);
|
||||||
|
}
|
||||||
|
if (vk_var) {
|
||||||
|
// Derive α_ij from α_jk (one level of recursion).
|
||||||
|
double a_jk = alpha_ij(ajk, aki, aij,
|
||||||
|
bj, bk, bi,
|
||||||
|
beta_j, beta_k, beta_i,
|
||||||
|
vj_var, vk_var, vi_var);
|
||||||
|
return PI - a_jk - beta_j;
|
||||||
|
}
|
||||||
|
// All ideal: closed-form formula.
|
||||||
|
return 0.5 * (PI + beta_k - beta_i - beta_j);
|
||||||
|
}
|
||||||
|
|
||||||
|
} // namespace conformallab
|
||||||
@@ -12,8 +12,8 @@ add_executable(conformallab_cgal_tests
|
|||||||
# ── Phase 3a: mesh infrastructure ──────────────────────────────────────
|
# ── Phase 3a: mesh infrastructure ──────────────────────────────────────
|
||||||
test_conformal_mesh.cpp
|
test_conformal_mesh.cpp
|
||||||
|
|
||||||
# ── Phase 3b: HyperIdealFunctional (to be added) ──────────────────────
|
# ── Phase 3b: HyperIdealFunctional ─────────────────────────────────────
|
||||||
# test_hyper_ideal_functional.cpp
|
test_hyper_ideal_functional.cpp
|
||||||
|
|
||||||
# ── Phase 3c: SphericalFunctional (to be added) ──────────────────────
|
# ── Phase 3c: SphericalFunctional (to be added) ──────────────────────
|
||||||
# test_spherical_functional.cpp
|
# test_spherical_functional.cpp
|
||||||
|
|||||||
184
code/tests/cgal/test_hyper_ideal_functional.cpp
Normal file
184
code/tests/cgal/test_hyper_ideal_functional.cpp
Normal file
@@ -0,0 +1,184 @@
|
|||||||
|
// test_hyper_ideal_functional.cpp
|
||||||
|
//
|
||||||
|
// Phase 3b — HyperIdealFunctional ported to ConformalMesh.
|
||||||
|
//
|
||||||
|
// Corresponds to de.varylab.discreteconformal.functional.HyperIdealFunctionalTest.
|
||||||
|
//
|
||||||
|
// The Java tests use CoHDS + HyperIdealGenerator to build complex meshes and
|
||||||
|
// then verify correctness via a finite-difference gradient check (FunctionalTest).
|
||||||
|
// Here we apply the same gradient-check strategy on simple hand-crafted meshes
|
||||||
|
// (triangle, tetrahedron) to validate the port independently of the generator.
|
||||||
|
//
|
||||||
|
// Test map (Java → C++)
|
||||||
|
// ──────────────────────
|
||||||
|
// testHessian (Ignored) → GradientCheck_Hessian (SKIPPED, not implemented)
|
||||||
|
// testGradientWithHyperIdeal… → GradientCheck_AllHyperIdealTriangle (ported)
|
||||||
|
// testGradientInExtendedDomain → GradientCheck_ExtendedDomain (ported)
|
||||||
|
// testGradientWithHyperelliptic → GradientCheck_TetrahedronAllVariable (ported)
|
||||||
|
// testFunctionalAtNaNValue → EnergyFiniteAtTestPoint (ported)
|
||||||
|
|
||||||
|
#include "conformal_mesh.hpp"
|
||||||
|
#include "mesh_builder.hpp"
|
||||||
|
#include "hyper_ideal_functional.hpp"
|
||||||
|
#include <gtest/gtest.h>
|
||||||
|
#include <cmath>
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
using namespace conformallab;
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Helpers
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
// Build a mesh with all vertices and edges variable, and set reasonable
|
||||||
|
// DOF values: b_i = b_val, a_e = a_val.
|
||||||
|
static std::vector<double> make_x_all_variable(
|
||||||
|
ConformalMesh& mesh, HyperIdealMaps& maps,
|
||||||
|
double b_val, double a_val)
|
||||||
|
{
|
||||||
|
int n = assign_all_dof_indices(mesh, maps);
|
||||||
|
std::vector<double> x(static_cast<std::size_t>(n));
|
||||||
|
|
||||||
|
// Vertices first, then edges (matching assign_all_dof_indices order)
|
||||||
|
for (auto v : mesh.vertices())
|
||||||
|
x[static_cast<std::size_t>(maps.v_idx[v])] = b_val;
|
||||||
|
for (auto e : mesh.edges())
|
||||||
|
x[static_cast<std::size_t>(maps.e_idx[e])] = a_val;
|
||||||
|
|
||||||
|
return x;
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// @Ignore in Java: no Hessian implemented
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(HyperIdealFunctional, GradientCheck_Hessian)
|
||||||
|
{
|
||||||
|
GTEST_SKIP() << "@Ignore in Java – Hessian not implemented in the functional";
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Gradient check: single triangle, all vertices hyper-ideal
|
||||||
|
//
|
||||||
|
// Mirrors testGradientWithHyperIdealAndIdealPoints (simplified to single face).
|
||||||
|
// DOFs: 3 vertex b-values + 3 edge a-values = 6 total.
|
||||||
|
// Point: b_i = 1.0, a_e = 0.5.
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(HyperIdealFunctional, GradientCheck_AllHyperIdealTriangle)
|
||||||
|
{
|
||||||
|
auto mesh = make_triangle();
|
||||||
|
auto maps = setup_hyper_ideal_maps(mesh);
|
||||||
|
auto x = make_x_all_variable(mesh, maps, /*b=*/1.0, /*a=*/0.5);
|
||||||
|
|
||||||
|
EXPECT_TRUE(gradient_check(mesh, x, maps))
|
||||||
|
<< "Finite-difference gradient check failed on all-hyper-ideal triangle";
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Gradient check in the "extended domain"
|
||||||
|
//
|
||||||
|
// Mirrors testGradientInTheExtendedDomain: larger DOF values
|
||||||
|
// (Java: x_i = 1.2 + |rnd|, so typically > 1.2).
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(HyperIdealFunctional, GradientCheck_ExtendedDomain)
|
||||||
|
{
|
||||||
|
auto mesh = make_triangle();
|
||||||
|
auto maps = setup_hyper_ideal_maps(mesh);
|
||||||
|
auto x = make_x_all_variable(mesh, maps, /*b=*/2.0, /*a=*/1.5);
|
||||||
|
|
||||||
|
EXPECT_TRUE(gradient_check(mesh, x, maps))
|
||||||
|
<< "Finite-difference gradient check failed in extended domain";
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Gradient check: tetrahedron (4 faces), all vertices and edges variable
|
||||||
|
//
|
||||||
|
// Mirrors testGradientWithHyperellipticCurve — larger mesh, closed surface.
|
||||||
|
// DOFs: 4 vertex b-values + 6 edge a-values = 10 total.
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(HyperIdealFunctional, GradientCheck_TetrahedronAllVariable)
|
||||||
|
{
|
||||||
|
auto mesh = make_tetrahedron();
|
||||||
|
auto maps = setup_hyper_ideal_maps(mesh);
|
||||||
|
auto x = make_x_all_variable(mesh, maps, /*b=*/1.0, /*a=*/0.5);
|
||||||
|
|
||||||
|
EXPECT_TRUE(gradient_check(mesh, x, maps))
|
||||||
|
<< "Finite-difference gradient check failed on all-variable tetrahedron";
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Energy is finite (not NaN / Inf)
|
||||||
|
//
|
||||||
|
// Mirrors testFunctionalAtNaNValue: evaluates at a test point and checks
|
||||||
|
// the energy is a real number. Uses a quad-strip to exercise the interior
|
||||||
|
// edge path (adjacent faces sharing one non-border edge).
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(HyperIdealFunctional, EnergyFiniteAtTestPoint)
|
||||||
|
{
|
||||||
|
auto mesh = make_quad_strip();
|
||||||
|
auto maps = setup_hyper_ideal_maps(mesh);
|
||||||
|
int n = assign_all_dof_indices(mesh, maps);
|
||||||
|
|
||||||
|
// Values taken from the Java testFunctionalAtNaNValue spirit:
|
||||||
|
// large-ish but positive DOF values that could expose degenerate paths.
|
||||||
|
std::vector<double> x(static_cast<std::size_t>(n), 1.5);
|
||||||
|
|
||||||
|
auto res = evaluate_hyper_ideal(mesh, x, maps, true, false);
|
||||||
|
|
||||||
|
EXPECT_FALSE(std::isnan(res.energy)) << "Energy must not be NaN";
|
||||||
|
EXPECT_FALSE(std::isinf(res.energy)) << "Energy must not be Inf";
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Gradient check: mixed vertices (some ideal = pinned, some hyper-ideal)
|
||||||
|
//
|
||||||
|
// One vertex is ideal (v_idx = -1, b = 0), the other two are hyper-ideal.
|
||||||
|
// This exercises the lij / αij branches for the ideal-vertex case.
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(HyperIdealFunctional, GradientCheck_MixedIdealHyperIdeal)
|
||||||
|
{
|
||||||
|
auto mesh = make_triangle();
|
||||||
|
auto maps = setup_hyper_ideal_maps(mesh);
|
||||||
|
|
||||||
|
// Make v0 ideal (pinned), v1 and v2 hyper-ideal.
|
||||||
|
auto vit = mesh.vertices().begin();
|
||||||
|
Vertex_index v0 = *vit++;
|
||||||
|
Vertex_index v1 = *vit++;
|
||||||
|
Vertex_index v2 = *vit;
|
||||||
|
|
||||||
|
maps.v_idx[v0] = -1; // ideal: b_0 = 0 (fixed)
|
||||||
|
maps.v_idx[v1] = 0;
|
||||||
|
maps.v_idx[v2] = 1;
|
||||||
|
|
||||||
|
// All edges variable: indices 2, 3, 4
|
||||||
|
int eidx = 2;
|
||||||
|
for (auto e : mesh.edges()) maps.e_idx[e] = eidx++;
|
||||||
|
|
||||||
|
// DOF vector: [b1, b2, a_e0, a_e1, a_e2]
|
||||||
|
std::vector<double> x = {1.0, 1.0, 0.5, 0.5, 0.5};
|
||||||
|
|
||||||
|
EXPECT_TRUE(gradient_check(mesh, x, maps))
|
||||||
|
<< "Finite-difference gradient check failed for mixed ideal / hyper-ideal";
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Gradient check: fan mesh (n=6), all variable
|
||||||
|
//
|
||||||
|
// Exercises the full halfedge-around-vertex traversal for a high-valence
|
||||||
|
// interior vertex.
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(HyperIdealFunctional, GradientCheck_Fan6AllVariable)
|
||||||
|
{
|
||||||
|
auto mesh = make_fan(6);
|
||||||
|
auto maps = setup_hyper_ideal_maps(mesh);
|
||||||
|
auto x = make_x_all_variable(mesh, maps, /*b=*/1.0, /*a=*/0.5);
|
||||||
|
|
||||||
|
EXPECT_TRUE(gradient_check(mesh, x, maps))
|
||||||
|
<< "Finite-difference gradient check failed on fan-6 mesh";
|
||||||
|
}
|
||||||
Reference in New Issue
Block a user