From 8c353bb884179f6c7dac2b4db9796fc394846957 Mon Sep 17 00:00:00 2001 From: Tarik Moussa Date: Tue, 12 May 2026 07:39:14 +0200 Subject: [PATCH] feat(phase3d+3e): port EuclideanCyclicFunctional + add SphericalFunctional gauge-fix MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Phase 3d — EuclideanCyclicFunctional: • euclidean_geometry.hpp: t-value / atan2 corner-angle formula with centering trick (μ = (Λ̃₁₂+Λ̃₂₃+Λ̃₃₁)/6) for numerical stability • euclidean_functional.hpp: EuclideanMaps bundle, gradient (G_v = Θ_v − Σα_v, G_e = α_opp⁺ + α_opp⁻ − φ_e), 10-point GL path-integral energy, gradient_check_euclidean — identical halfedge convention to SphericalFunctional • test_euclidean_functional.cpp: 11 tests (1 skip) covering angle formula, right-isosceles triangle, angle sum = π, degenerate detection, gradient checks on triangle/quad-strip/tetrahedron/fan-5/mixed-pinned, NaN check Phase 3e — Spherical gauge-fix: • spherical_gauge_shift(): Newton + backtracking line search to find t* where ΣG_v(x + t·1) = 0 (maximises E along the global scale direction); bisection used when sign change is detectable, Newton+backtrack otherwise • apply_spherical_gauge(): in-place wrapper • 3 new tests: GaugeFix_SpherTetVertexZerosSumGv, GaugeFix_ApplyInPlace, GaugeFix_AlreadyAtGaugeReturnsTNearZero Total: 45 cgal tests pass, 3 skipped (@Ignore Hessian stubs, one per functional) Co-Authored-By: Claude Sonnet 4.6 --- README.md | 41 ++- code/include/euclidean_functional.hpp | 325 ++++++++++++++++++ code/include/euclidean_geometry.hpp | 91 +++++ code/include/spherical_functional.hpp | 118 +++++++ code/tests/cgal/CMakeLists.txt | 5 +- code/tests/cgal/test_euclidean_functional.cpp | 269 +++++++++++++++ code/tests/cgal/test_spherical_functional.cpp | 96 +++++- 7 files changed, 924 insertions(+), 21 deletions(-) create mode 100644 code/include/euclidean_functional.hpp create mode 100644 code/include/euclidean_geometry.hpp create mode 100644 code/tests/cgal/test_euclidean_functional.cpp diff --git a/README.md b/README.md index dc38685..5f4c94f 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ conformallab++ is a modern C++ reimplementation of the [ConformalLab](https://gi The long-term goal is a **CGAL package** that brings discrete conformal maps (hyper-ideal, spherical, Euclidean) to the CGAL ecosystem using `CGAL::Surface_mesh` as the underlying half-edge data structure. -> **Status:** Phase 3 abgeschlossen. Kern-Geometrie (Clausen, Hyper-Ideal, Sphärisch) und CGAL-Mesh-Infrastruktur sind vorhanden. Zum Abschluss von Phase 3 und als Vorbereitung für Phase 4 stehen noch sechs konkrete Schritte aus — siehe [Roadmap](#roadmap). +> **Status:** Phasen 3d + 3e abgeschlossen. Alle drei Kern-Geometrien (Hyper-Ideal, Sphärisch, Euklidisch) plus Sphärischer Gauge-Fix sind auf `ConformalMesh` portiert. **45 Tests, 3 skipped** (Hessian-Stubs). Nächster Schritt: Phase 3f (Hessian) → Phase 4 (Newton-Solver) — siehe [Roadmap](#roadmap). --- @@ -17,7 +17,8 @@ The long-term goal is a **CGAL package** that brings discrete conformal maps (hy | Hyper-ideal functional (energy + gradient, CGAL mesh) | ✅ Phase 3b | | Spherical functional (energy + gradient, CGAL mesh) | ✅ Phase 3c | | CGAL `Surface_mesh` infrastructure + mesh builders | ✅ Phase 3a | -| Euclidean functional | 🔜 Phase 4 | +| Euclidean functional (energy + gradient, CGAL mesh) | ✅ Phase 3d | +| Spherical gauge-fix (scale gauge for closed surfaces) | ✅ Phase 3e | | Solvers (Newton, gradient flow) | 🔜 Phase 4 | | XML serialisation / mesh I/O | 🔜 Phase 5 | | Full CLI app | 🔜 Phase 5 | @@ -72,7 +73,7 @@ cmake --build build --target conformallab_cgal_tests -j$(nproc) ctest --test-dir build -R "^cgal\." --output-on-failure ``` -Expected output: **30 tests pass, 2 skipped** (the two `@Ignore` Hessian stubs). +Expected output: **45 tests pass, 3 skipped** (the three `@Ignore` Hessian stubs, one per functional). ### Full CLI app (CGAL + viewer) @@ -93,7 +94,9 @@ cmake --build build -j$(nproc) | `hyper_ideal_utility.hpp` | Tetrahedron volume (Meyerhoff / Kolpakov–Mednykh) | | `hyper_ideal_functional.hpp` | Hyper-ideal energy + gradient on `ConformalMesh` | | `spherical_geometry.hpp` | Spherical arc length, half-angle angle formula | -| `spherical_functional.hpp` | Spherical energy + gradient on `ConformalMesh` | +| `spherical_functional.hpp` | Spherical energy + gradient + gauge-fix on `ConformalMesh` | +| `euclidean_geometry.hpp` | Euclidean corner-angle formula (t-value / atan2) | +| `euclidean_functional.hpp` | Euclidean energy + gradient on `ConformalMesh` | | `conformal_mesh.hpp` | `ConformalMesh` = `CGAL::Surface_mesh`, index types, property-map helpers | | `mesh_builder.hpp` | Factory meshes: triangle, tetrahedron, quad strip, fan, spherical tetrahedron, octahedron face | | `discrete_elliptic_utility.hpp` | Discrete elliptic integrals | @@ -112,7 +115,9 @@ code/ │ ├── hyper_ideal_geometry.hpp # ζ, lᵢⱼ, αᵢⱼ — pure math │ ├── hyper_ideal_functional.hpp # HyperIdealFunctional on ConformalMesh │ ├── spherical_geometry.hpp # spherical arc length + angles -│ ├── spherical_functional.hpp # SphericalFunctional on ConformalMesh +│ ├── spherical_functional.hpp # SphericalFunctional + gauge-fix on ConformalMesh +│ ├── euclidean_geometry.hpp # Euclidean corner-angle formula (t-value) +│ ├── euclidean_functional.hpp # EuclideanCyclicFunctional on ConformalMesh │ ├── clausen.hpp # Clausen / Lobachevsky / ImLi₂ │ └── hyper_ideal_utility.hpp # Tetrahedron volumes ├── src/ @@ -125,7 +130,8 @@ code/ │ ├── CMakeLists.txt │ ├── test_conformal_mesh.cpp # 14 mesh infrastructure tests │ ├── test_hyper_ideal_functional.cpp # 6 hyper-ideal gradient checks -│ └── test_spherical_functional.cpp # 8 spherical gradient checks +│ ├── test_spherical_functional.cpp # 11 spherical gradient + gauge-fix checks +│ └── test_euclidean_functional.cpp # 11 Euclidean gradient checks └── deps/ ├── tarballs/ # bundled dependency archives ├── eigen-3.4.0/ # header-only linear algebra (always extracted) @@ -159,7 +165,8 @@ CGAL `Surface_mesh` tests (test prefix `cgal.`): | `ConformalMeshProperties` | 5 | Property maps: λ, θ, idx, α, geometry type | | `ConformalMeshValidity` | 1 | CGAL validity check for all factory meshes | | `HyperIdealFunctional` | 6 | Finite-difference gradient checks on triangle, tetrahedron, quad strip, fan | -| `SphericalFunctional` | 8 | Angle formula correctness + gradient checks on spherical meshes | +| `SphericalFunctional` | 11 | Angle formula + gradient checks on spherical meshes + 3 gauge-fix tests | +| `EuclideanFunctional` | 11 | Angle formula + gradient checks on Euclidean meshes (triangle, quad strip, fan, tetrahedron) | --- @@ -200,19 +207,15 @@ bereiten Phase 4 (Solver) vor. Die Zeitangaben sind grobe Schätzungen. --- ``` -Phase 3d EuclideanCyclicFunctional portieren (2–3 Tage) - → Wichtigster fehlender Baustein: Euklidische Uniformisierung - (genus-0 Flachparameterisierung, Kotangenten-Laplace) - → Energie + Gradient auf ConformalMesh - → Entsperrt die verbleibenden HDS-blockierten Tests - → Voraussetzung für jeden end-to-end Optimiererlauf +Phase 3d EuclideanCyclicFunctional portieren ✅ abgeschlossen + → euclidean_geometry.hpp: t-Wert / atan2 Winkelformel + → euclidean_functional.hpp: Energie + Gradient auf ConformalMesh + → 11 Tests: Winkelformel, Gradient-Checks, NaN-Test, Fan, Mixed-Pinned -Phase 3e Gauge-Fix für SphericalFunctional (1 Tag) - → Java-Äquivalent: maximizeInNegativeDirection() vor jeder - Auswertung — 1D-Liniensuche (Brent) entlang der - globalen Additionskonstante in u_v - → Ohne Gauge-Fix ist die Energie auf geschlossenen Flächen - nach unten unbeschränkt; kein Optimierer kann konvergieren +Phase 3e Gauge-Fix für SphericalFunctional ✅ abgeschlossen + → spherical_gauge_shift() + apply_spherical_gauge() in + spherical_functional.hpp — Newton + Backtracking + → 3 neue Tests: ZerosSumGv, ApplyInPlace, AlreadyAtGauge Phase 3f Hessian (Sphärisch → Euklidisch) (2–3 Tage) → Analytischer Hessian = Kotangenten-Laplace-Operator diff --git a/code/include/euclidean_functional.hpp b/code/include/euclidean_functional.hpp new file mode 100644 index 0000000..8ff6203 --- /dev/null +++ b/code/include/euclidean_functional.hpp @@ -0,0 +1,325 @@ +#pragma once +// euclidean_functional.hpp +// +// Energy and gradient of the Euclidean discrete conformal functional +// (EuclideanCyclicFunctional) evaluated on a ConformalMesh. +// +// Ported from de.varylab.discreteconformal.functional.EuclideanCyclicFunctional. +// +// ┌──────────────────────────────────────────────────────────────────────────┐ +// │ DOFs │ +// │ x[v_idx[v]] = u_v – conformal factor at vertex v │ +// │ x[e_idx[e]] = λ_e – edge log-length variable (optional) │ +// │ -1 means "pinned" (u_v = 0 / λ_e = 0, only λ° contributes) │ +// │ │ +// │ Effective log-length (always additive, unlike SphericalFunctional): │ +// │ Λ̃_ij = λ°_ij + u_i + u_j + (x[e_idx[e]] if variable, else 0) │ +// │ │ +// │ Side length: l_ij = exp(Λ̃_ij / 2) │ +// │ │ +// │ Gradient: │ +// │ ∂E/∂u_v = Θ_v − Σ_{faces adj. v} α_v(face) │ +// │ ∂E/∂λ_e = α_opp(face⁺) + α_opp(face⁻) − φ_e │ +// │ │ +// │ Energy: │ +// │ Computed as the path integral E(x) = ∫₀¹ ⟨G(tx), x⟩ dt │ +// │ using 10-point Gauss-Legendre quadrature (same as SphericalFunctional)│ +// │ This is E(0)=0 by construction and exact for conservative G. │ +// └──────────────────────────────────────────────────────────────────────────┘ +// +// Halfedge convention (identical to SphericalFunctional): +// h0 = mesh.halfedge(f), h1 = next(h0), h2 = next(h1) +// v1 = source(h0), v2 = source(h1), v3 = source(h2) +// h_alpha[h0] = α3 (angle at v3, opposite edge h0 = v1v2) +// h_alpha[h1] = α1 (angle at v1, opposite edge h1 = v2v3) +// h_alpha[h2] = α2 (angle at v2, opposite edge h2 = v3v1) +// +// Property-map name prefix: "ev:" (vertex) and "ee:" (edge). + +#include "conformal_mesh.hpp" +#include "euclidean_geometry.hpp" +#include +#include +#include +#include + +namespace conformallab { + +// ── Property-map type aliases ───────────────────────────────────────────────── + +using EuclVMapD = ConformalMesh::Property_map; +using EuclVMapI = ConformalMesh::Property_map; +using EuclEMapD = ConformalMesh::Property_map; +using EuclEMapI = ConformalMesh::Property_map; + +// ── Persistent map bundle ───────────────────────────────────────────────────── + +struct EuclideanMaps { + EuclVMapI v_idx; ///< DOF index per vertex (-1 = pinned / u_v = 0) + EuclEMapI e_idx; ///< DOF index per edge (-1 = no edge DOF) + EuclVMapD theta_v; ///< target cone angle Θ_v (default 2π) + EuclEMapD phi_e; ///< target edge turn angle φ_e (default π) + EuclEMapD lambda0; ///< base log-length λ°_e (default 0.0) +}; + +// Create and attach property maps with sensible defaults. +// theta_v = 2π (flat vertex), phi_e = π (interior edge, flat surface). +inline EuclideanMaps setup_euclidean_maps(ConformalMesh& mesh) +{ + constexpr double TWO_PI = 2.0 * 3.14159265358979323846; + constexpr double PI = 3.14159265358979323846; + + EuclideanMaps m; + m.v_idx = mesh.add_property_map ("ev:idx", -1 ).first; + m.e_idx = mesh.add_property_map ("ee:idx", -1 ).first; + m.theta_v= mesh.add_property_map("ev:theta", TWO_PI ).first; + m.phi_e = mesh.add_property_map("ee:phi", PI ).first; + m.lambda0= mesh.add_property_map("ee:lam0", 0.0 ).first; + return m; +} + +// Assign DOF indices 0..n-1 for all vertices only (no edge DOFs). +inline int assign_euclidean_vertex_dof_indices(ConformalMesh& mesh, EuclideanMaps& m) +{ + int idx = 0; + for (auto v : mesh.vertices()) m.v_idx[v] = idx++; + return idx; +} + +// Assign DOF indices for all vertices AND all edges. +inline int assign_euclidean_all_dof_indices(ConformalMesh& mesh, EuclideanMaps& 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; +} + +// Count variable DOFs (vertices + edges). +inline int euclidean_dimension(const ConformalMesh& mesh, const EuclideanMaps& 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; +} + +// Set lambda0 from mesh vertex positions (Euclidean): +// λ°_e = 2·log(|p_i − p_j|) (natural log of Euclidean edge length squared) +// +// This gives exp(Λ̃_ij / 2) = l_ij at x=0. +inline void compute_euclidean_lambda0_from_mesh(ConformalMesh& mesh, EuclideanMaps& m) +{ + for (auto e : mesh.edges()) { + auto h = mesh.halfedge(e); + auto p1 = mesh.point(mesh.source(h)); + auto p2 = mesh.point(mesh.target(h)); + double dx = p1.x() - p2.x(); + double dy = p1.y() - p2.y(); + double dz = p1.z() - p2.z(); + double len = std::sqrt(dx*dx + dy*dy + dz*dz); + if (len > 1e-15) + m.lambda0[e] = 2.0 * std::log(len); + else + m.lambda0[e] = -30.0; // degenerate edge + } +} + +// ── Internal helpers ────────────────────────────────────────────────────────── + +static inline double eucl_dof_val(int idx, const std::vector& x) +{ + return idx >= 0 ? x[static_cast(idx)] : 0.0; +} + +static inline std::size_t eucl_hidx(Halfedge_index h) +{ + return static_cast(static_cast(h)); +} + +// ── Gradient ────────────────────────────────────────────────────────────────── +// +// G_v = Θ_v − Σ_{faces adj. v} α_v(face) +// G_e = α_opp(face⁺) + α_opp(face⁻) − φ_e +// +// Corner-angle storage (h_alpha): +// h_alpha[h] = corner angle OPPOSITE to the edge of halfedge h in its face. +// h_alpha[h0] = α3, h_alpha[h1] = α1, h_alpha[h2] = α2 +// (same convention as SphericalFunctional) +inline std::vector euclidean_gradient( + ConformalMesh& mesh, + const std::vector& x, + const EuclideanMaps& m) +{ + const int n = euclidean_dimension(mesh, m); + std::vector G(static_cast(n), 0.0); + + // Per-halfedge corner-angle storage. + const std::size_t nh = mesh.number_of_halfedges(); + std::vector h_alpha(nh, 0.0); + + // ── Pass 1: compute corner angles 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); + + 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); + + // Effective log-lengths (always additive: u_i + u_j regardless of DOF status) + double u1 = eucl_dof_val(m.v_idx[v1], x); + double u2 = eucl_dof_val(m.v_idx[v2], x); + double u3 = eucl_dof_val(m.v_idx[v3], x); + + double lam12 = m.lambda0[e12] + u1 + u2 + eucl_dof_val(m.e_idx[e12], x); + double lam23 = m.lambda0[e23] + u2 + u3 + eucl_dof_val(m.e_idx[e23], x); + double lam31 = m.lambda0[e31] + u3 + u1 + eucl_dof_val(m.e_idx[e31], x); + + auto fa = euclidean_angles(lam12, lam23, lam31); + if (!fa.valid) continue; // degenerate triangle: contributes 0 + + // h_alpha[h] = corner angle OPPOSITE to h's edge: + // h0 (edge v1v2) → opposite corner at v3 → α3 + // h1 (edge v2v3) → opposite corner at v1 → α1 + // h2 (edge v3v1) → opposite corner at v2 → α2 + h_alpha[eucl_hidx(h0)] = fa.alpha3; + h_alpha[eucl_hidx(h1)] = fa.alpha1; + h_alpha[eucl_hidx(h2)] = fa.alpha2; + } + + // ── Pass 2: accumulate vertex gradient ─────────────────────────────────── + // G_v = Θ_v − Σ h_alpha[prev(h)] for each incoming non-border h to v. + for (auto v : mesh.vertices()) { + int iv = m.v_idx[v]; + if (iv < 0) continue; + double sum_alpha = 0.0; + for (auto h : CGAL::halfedges_around_target(v, mesh)) { + if (mesh.is_border(h)) continue; + sum_alpha += h_alpha[eucl_hidx(mesh.prev(h))]; + } + G[static_cast(iv)] = m.theta_v[v] - sum_alpha; + } + + // ── Pass 3: accumulate edge gradient ───────────────────────────────────── + // G_e = α_opp(f⁺) + α_opp(f⁻) − φ_e + // α_opp of edge e in face f = h_alpha[halfedge h of e pointing INTO f]. + 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); + double sum = -m.phi_e[e]; + if (!mesh.is_border(h)) sum += h_alpha[eucl_hidx(h)]; + if (!mesh.is_border(ho)) sum += h_alpha[eucl_hidx(ho)]; + G[static_cast(ie)] = sum; + } + + return G; +} + +// ── Energy via Gauss-Legendre path integral ─────────────────────────────────── +// +// E(x) = ∫₀¹ ⟨G(tx), x⟩ dt (10-point GL quadrature, same as SphericalFunctional) +inline double euclidean_energy( + ConformalMesh& mesh, + const std::vector& x, + const EuclideanMaps& m) +{ + static const double gl_s[10] = { + -0.9739065285171717, -0.8650633666889845, + -0.6794095682990244, -0.4333953941292472, + -0.1488743389816312, 0.1488743389816312, + 0.4333953941292472, 0.6794095682990244, + 0.8650633666889845, 0.9739065285171717 + }; + static const double gl_w[10] = { + 0.0666713443086881, 0.1494513491505806, + 0.2190863625159820, 0.2692667193099963, + 0.2955242247147529, 0.2955242247147529, + 0.2692667193099963, 0.2190863625159820, + 0.1494513491505806, 0.0666713443086881 + }; + + const std::size_t n = x.size(); + double E = 0.0; + + for (int k = 0; k < 10; ++k) { + double t = (1.0 + gl_s[k]) * 0.5; + double wt = gl_w[k] * 0.5; + + std::vector tx(n); + for (std::size_t i = 0; i < n; ++i) tx[i] = t * x[i]; + + auto G = euclidean_gradient(mesh, tx, m); + + double dot = 0.0; + for (std::size_t i = 0; i < n; ++i) dot += G[i] * x[i]; + E += wt * dot; + } + return E; +} + +// ── Full evaluation (energy + gradient) ────────────────────────────────────── + +struct EuclideanResult { + double energy = 0.0; + std::vector gradient; +}; + +inline EuclideanResult evaluate_euclidean( + ConformalMesh& mesh, + const std::vector& x, + const EuclideanMaps& m, + bool need_energy = true, + bool need_gradient = true) +{ + EuclideanResult res; + if (need_gradient) + res.gradient = euclidean_gradient(mesh, x, m); + if (need_energy) + res.energy = euclidean_energy(mesh, x, m); + return res; +} + +// ── Finite-difference gradient check ───────────────────────────────────────── +// +// Tests |G[i] − (E(x+εeᵢ) − E(x−εeᵢ))/(2ε)| / max(1,|G[i]|) < tol +// for all variable DOFs. +inline bool gradient_check_euclidean( + ConformalMesh& mesh, + const std::vector& x0, + const EuclideanMaps& m, + double eps = 1e-5, + double tol = 1e-4) +{ + auto G = euclidean_gradient(mesh, x0, m); + const int n = static_cast(G.size()); + + std::vector xp = x0, xm = x0; + bool ok = true; + + for (int i = 0; i < n; ++i) { + std::size_t si = static_cast(i); + xp[si] = x0[si] + eps; + xm[si] = x0[si] - eps; + + double Ep = euclidean_energy(mesh, xp, m); + double Em = euclidean_energy(mesh, xm, m); + + xp[si] = xm[si] = x0[si]; // restore + + 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 diff --git a/code/include/euclidean_geometry.hpp b/code/include/euclidean_geometry.hpp new file mode 100644 index 0000000..da636da --- /dev/null +++ b/code/include/euclidean_geometry.hpp @@ -0,0 +1,91 @@ +#pragma once +// euclidean_geometry.hpp +// +// Corner-angle formula for Euclidean triangles in the discrete conformal +// (log-length) parametrisation. +// +// Ported from de.varylab.discreteconformal.functional.EuclideanCyclicFunctional. +// +// In the discrete conformal parametrisation a Euclidean triangle is described by +// its three effective log-lengths Λ̃_ij = λ°_ij + u_i + u_j (+ edge DOF). +// The corresponding side lengths are l_ij = exp(Λ̃_ij / 2). +// +// Vertex ordering convention (matches EuclideanCyclicFunctional.java): +// v1 is opposite edge l23, v2 is opposite l31, v3 is opposite l12. +// +// t-value trick (Springborn 2008 §3): +// t12 = −l12 + l23 + l31 = 2(s − l12) +// t23 = +l12 − l23 + l31 = 2(s − l23) +// t31 = +l12 + l23 − l31 = 2(s − l31) +// denom = sqrt(t12 · t23 · t31 · l123) = 4 · Area +// +// α_v = 2 · atan2( product of t-values adjacent to v, denom ) +// +// The centering trick (l_ij ← exp((Λ̃_ij − 2·μ)/2), μ = (Λ̃12+Λ̃23+Λ̃31)/6) +// rescales all three sides by the same factor, leaving angles unchanged but +// keeping the arguments of exp in a safe numerical range. + +#include + +namespace conformallab { + +struct EuclideanFaceAngles { + double alpha1; ///< corner angle at v1 (opposite l23) + double alpha2; ///< corner angle at v2 (opposite l31) + double alpha3; ///< corner angle at v3 (opposite l12) + bool valid; +}; + +// ── From side lengths ───────────────────────────────────────────────────────── +// +// Given three Euclidean side lengths l12, l23, l31 > 0 satisfying the triangle +// inequality, compute the corner angles. +// +// Returns valid=false if the triangle inequality is violated (any t-value ≤ 0). +inline EuclideanFaceAngles euclidean_angles_from_lengths( + double l12, double l23, double l31) +{ + const double t12 = -l12 + l23 + l31; // 2*(s − l12) + const double t23 = +l12 - l23 + l31; // 2*(s − l23) + const double t31 = +l12 + l23 - l31; // 2*(s − l31) + + if (t12 <= 0.0 || t23 <= 0.0 || t31 <= 0.0) + return {0.0, 0.0, 0.0, false}; + + const double l123 = l12 + l23 + l31; + const double denom2 = t12 * t23 * t31 * l123; // = (4·Area)² + if (denom2 <= 0.0) + return {0.0, 0.0, 0.0, false}; + + const double denom = std::sqrt(denom2); + + // α at v1 (opposite l23): adjacent t-values are t12 and t31 + // α at v2 (opposite l31): adjacent t-values are t12 and t23 + // α at v3 (opposite l12): adjacent t-values are t23 and t31 + return { + 2.0 * std::atan2(t12 * t31, denom), + 2.0 * std::atan2(t12 * t23, denom), + 2.0 * std::atan2(t23 * t31, denom), + true + }; +} + +// ── From effective log-lengths Λ̃ ───────────────────────────────────────────── +// +// Converts to side lengths l_ij = exp(Λ̃_ij / 2), applying the centering +// trick for numerical safety, then delegates to euclidean_angles_from_lengths. +// +// The centering constant μ = (Λ̃12 + Λ̃23 + Λ̃31) / 6 ensures +// l12 · l23 · l31 = 1 (geometric mean = 1) +// which keeps all l values near 1 and prevents float overflow for large |Λ̃|. +inline EuclideanFaceAngles euclidean_angles( + double lam12, double lam23, double lam31) +{ + const double mu = (lam12 + lam23 + lam31) / 6.0; + const double l12 = std::exp((lam12 - 2.0 * mu) * 0.5); + const double l23 = std::exp((lam23 - 2.0 * mu) * 0.5); + const double l31 = std::exp((lam31 - 2.0 * mu) * 0.5); + return euclidean_angles_from_lengths(l12, l23, l31); +} + +} // namespace conformallab diff --git a/code/include/spherical_functional.hpp b/code/include/spherical_functional.hpp index d8f1015..36ee0aa 100644 --- a/code/include/spherical_functional.hpp +++ b/code/include/spherical_functional.hpp @@ -354,4 +354,122 @@ inline bool gradient_check_spherical( return ok; } +// ── Gauge-fix for closed spherical surfaces ─────────────────────────────────── +// +// On a closed (boundaryless) spherical surface the energy E(u + t·1) has a +// unique maximum w.r.t. t ∈ ℝ (the "global scale" gauge mode). Without +// fixing this gauge the functional is unbounded below and no solver converges. +// +// This function returns the scalar shift t* such that +// Σ_v G_v(x + t*·1_v) = 0 +// i.e., the derivative of E(u+t·1) w.r.t. t is zero at t = t*. +// +// Apply the shift by adding t* to every vertex DOF in x. +// +// Implementation: bisection on f(t) = Σ_v G_v(x + t·1_v). +// f is strictly monotone decreasing (second derivative < 0) for a convex +// functional, so bisection converges in O(log₂(2·bracket/tol)) iterations. +// +// Parameters: +// bracket – initial search interval [−bracket, +bracket] (default 50) +// tol – absolute tolerance on t* (default 1e-8) +// +// Returns 0.0 if the zero cannot be bracketed (already at gauge maximum, +// or open surface — no shift needed). +inline double spherical_gauge_shift( + ConformalMesh& mesh, + const std::vector& x, + const SphericalMaps& m, + double bracket = 50.0, + double tol = 1e-8) +{ + // Helper: sum of all vertex gradient components at x + t·1_v. + auto sum_Gv = [&](double t) -> double { + // Build a shifted copy of x (only vertex DOFs shifted). + std::vector xt = x; + for (auto v : mesh.vertices()) { + int iv = m.v_idx[v]; + if (iv >= 0) + xt[static_cast(iv)] += t; + } + auto G = spherical_gradient(mesh, xt, m); + double sum = 0.0; + for (auto v : mesh.vertices()) { + int iv = m.v_idx[v]; + if (iv >= 0) + sum += G[static_cast(iv)]; + } + return sum; + }; + + // ── Try bisection first (works when there is a sign change in [−bracket, +bracket]) ── + double a = -bracket, b = bracket; + double fa = sum_Gv(a), fb = sum_Gv(b); + + if (fa * fb <= 0.0) { + for (int iter = 0; iter < 120; ++iter) { + double c = 0.5 * (a + b); + double fc = sum_Gv(c); + if (std::abs(fc) < tol || (b - a) < tol) return c; + if (fa * fc < 0.0) { b = c; fb = fc; } + else { a = c; fa = fc; } + } + return 0.5 * (a + b); + } + + // ── No sign change (zero may lie at a domain boundary). ─────────────────── + // Use damped Newton's method with backtracking line search. + // f'(t) estimated by forward finite difference. + // When the Newton step overshoots the valid domain (ΣG_v jumps back up + // because faces become degenerate), backtracking halves the step until + // |f| strictly decreases. + const double fd_eps = 1e-5; + double t = 0.0; + double ft = sum_Gv(t); + + for (int iter = 0; iter < 120; ++iter) { + if (std::abs(ft) < tol) return t; + + double ftp = sum_Gv(t + fd_eps); + double dft = (ftp - ft) / fd_eps; + if (std::abs(dft) < 1e-14) return t; // gradient flat — give up + + double dt_raw = -ft / dft; + + // Backtracking line search: halve dt until |f| decreases. + double alpha = 1.0; + bool improved = false; + for (int back = 0; back < 40; ++back) { + double t_try = t + alpha * dt_raw; + t_try = std::max(-bracket, std::min(bracket, t_try)); + double ft_try = sum_Gv(t_try); + if (std::abs(ft_try) < std::abs(ft)) { + t = t_try; + ft = ft_try; + improved = true; + break; + } + alpha *= 0.5; + } + if (!improved) return t; // cannot reduce further + } + return t; +} + +// Apply the gauge shift in-place: x_v ← x_v + t* for all variable vertices. +inline void apply_spherical_gauge( + ConformalMesh& mesh, + std::vector& x, + const SphericalMaps& m, + double bracket = 50.0, + double tol = 1e-8) +{ + double t = spherical_gauge_shift(mesh, x, m, bracket, tol); + for (auto v : mesh.vertices()) { + int iv = m.v_idx[v]; + if (iv >= 0) + x[static_cast(iv)] += t; + } +} + } // namespace conformallab diff --git a/code/tests/cgal/CMakeLists.txt b/code/tests/cgal/CMakeLists.txt index 8ad0330..2cc999b 100644 --- a/code/tests/cgal/CMakeLists.txt +++ b/code/tests/cgal/CMakeLists.txt @@ -15,8 +15,11 @@ add_executable(conformallab_cgal_tests # ── Phase 3b: HyperIdealFunctional ───────────────────────────────────── test_hyper_ideal_functional.cpp - # ── Phase 3c: SphericalFunctional ───────────────────────────────────── + # ── Phase 3c + 3e: SphericalFunctional + gauge-fix ──────────────────── test_spherical_functional.cpp + + # ── Phase 3d: EuclideanCyclicFunctional ─────────────────────────────── + test_euclidean_functional.cpp ) target_include_directories(conformallab_cgal_tests SYSTEM PRIVATE diff --git a/code/tests/cgal/test_euclidean_functional.cpp b/code/tests/cgal/test_euclidean_functional.cpp new file mode 100644 index 0000000..9337f4c --- /dev/null +++ b/code/tests/cgal/test_euclidean_functional.cpp @@ -0,0 +1,269 @@ +// 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 (SKIPPED) +// 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 +#include +#include + +using namespace conformallab; + +// ════════════════════════════════════════════════════════════════════════════ +// @Ignore in Java: no Hessian implemented yet +// ════════════════════════════════════════════════════════════════════════════ + +TEST(EuclideanFunctional, GradientCheck_Hessian) +{ + GTEST_SKIP() << "@Ignore in Java – Hessian not yet implemented"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// 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 x(static_cast(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 x(static_cast(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 x(static_cast(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 x(static_cast(n), 0.0); + // Set vertex DOFs slightly negative to keep triangles non-degenerate. + for (int i = 0; i < 4; ++i) + x[static_cast(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 x(static_cast(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 x(static_cast(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 x = {-0.1, -0.3, -0.2}; + + EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps)) + << "Gradient check failed for mixed pinned/variable vertices"; +} diff --git a/code/tests/cgal/test_spherical_functional.cpp b/code/tests/cgal/test_spherical_functional.cpp index 95d632c..4e5533d 100644 --- a/code/tests/cgal/test_spherical_functional.cpp +++ b/code/tests/cgal/test_spherical_functional.cpp @@ -1,4 +1,4 @@ -// test_spherical_functional.cpp +// test_spherical_functional.cpp (Phase 3c + 3e) // // Phase 3c — SphericalFunctional ported to ConformalMesh. // @@ -244,3 +244,97 @@ TEST(SphericalFunctional, GradientCheck_MixedPinnedVertices) 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 x(static_cast(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(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 x_fixed = x; + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv >= 0) + x_fixed[static_cast(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(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 x(static_cast(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(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 x(static_cast(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; +}