diff --git a/.gitea/workflows/doxygen-pages.yml b/.gitea/workflows/doxygen-pages.yml index fe65d82..f338277 100644 --- a/.gitea/workflows/doxygen-pages.yml +++ b/.gitea/workflows/doxygen-pages.yml @@ -47,8 +47,12 @@ jobs: head -30 doc/doxygen/doxygen-warnings.log fi - - name: Report Doxygen coverage - run: bash scripts/doxygen-coverage.sh + - name: Enforce Doxygen coverage 100% + # Coverage is measured against every public symbol under + # code/include/ (the `detail::` namespaces are excluded). As of + # the `docs/doxygen-coverage-100` PR the baseline is 100 %, so + # the gate fires only on regressions. + run: bash scripts/doxygen-coverage.sh --threshold 100 - name: Regenerate doc/api/headers.md from XML run: | diff --git a/code/include/CGAL/Conformal_map_traits.h b/code/include/CGAL/Conformal_map_traits.h index 7036b44..51bd2ee 100644 --- a/code/include/CGAL/Conformal_map_traits.h +++ b/code/include/CGAL/Conformal_map_traits.h @@ -80,10 +80,13 @@ namespace CGAL { // Default_conformal_map_traits — primary template (undefined) // ════════════════════════════════════════════════════════════════════════════ // -// The undefined primary template forces specialisation per mesh type. -// MVP provides only the `Surface_mesh` specialisation below; further -// mesh types (Polyhedron_3, OpenMesh, pmp) are deferred to Phase 8a.2. - +/*! +\ingroup PkgConformalMapConcepts +\brief Primary `ConformalMapTraits` template — undefined, must be +specialised per mesh type. The MVP only ships the `Surface_mesh` +specialisation below; further mesh types (Polyhedron_3, OpenMesh, +pmp) are deferred to Phase 8a.2. +*/ template > struct Default_conformal_map_traits; diff --git a/code/include/clausen.hpp b/code/include/clausen.hpp index 29e7c5d..5337cc0 100644 --- a/code/include/clausen.hpp +++ b/code/include/clausen.hpp @@ -128,9 +128,9 @@ inline int ncl5pi6() noexcept { } // namespace detail -// Clausen's integral Cl2(x) = -integral_0^x log|2 sin(t/2)| dt. -// High-precision Chebyshev implementation. -// Corresponds to Java Clausen.clausen2(). +/// Clausen's integral `Cl₂(x) = −∫₀ˣ log|2 sin(t/2)| dt`, +/// computed via a high-precision Chebyshev expansion. +/// Same as Java `Clausen.clausen2()`. inline double clausen2(double x) noexcept { using namespace detail; constexpr double pi = 3.14159265358979323846264338328; @@ -154,8 +154,8 @@ inline double clausen2(double x) noexcept { return rh ? -f : f; } -// Milnor's Lobachevsky function Л(x) = Cl2(2x)/2. -// Corresponds to Java Clausen.Л(). +/// Milnor's Lobachevsky function `Л(x) = Cl₂(2x) / 2`. +/// Same as Java `Clausen.Л()`. inline double Lobachevsky(double x) noexcept { constexpr double pi = 3.14159265358979323846264338328; x = std::fmod(x, pi); @@ -163,8 +163,8 @@ inline double Lobachevsky(double x) noexcept { return clausen2(2.0 * x) / 2.0; } -// Imaginary part of the dilogarithm Im(Li2(z)). -// Corresponds to Java Clausen.ImLi2(). +/// Imaginary part of the dilogarithm `Im(Li₂(z))` for complex `z`. +/// Same as Java `Clausen.ImLi2()`. inline double ImLi2(std::complex z) noexcept { auto a = std::log(1.0 - std::conj(z)); // log(1 - conj(z)) auto b = std::log(1.0 - z); // log(1 - z) diff --git a/code/include/cp_euclidean_functional.hpp b/code/include/cp_euclidean_functional.hpp index 82efb31..49d6549 100644 --- a/code/include/cp_euclidean_functional.hpp +++ b/code/include/cp_euclidean_functional.hpp @@ -77,12 +77,17 @@ namespace conformallab { // ── Property-map type aliases ──────────────────────────────────────────────── +/// Property map face → `int` for the CP-Euclidean functional. using CPFMapI = ConformalMesh::Property_map; +/// Property map face → `double` for the CP-Euclidean functional. using CPFMapD = ConformalMesh::Property_map; +/// Property map edge → `double` for the CP-Euclidean functional. using CPEMapD = ConformalMesh::Property_map; // ── Persistent map bundle ───────────────────────────────────────────────────── +/// Bundle of the three property maps consumed by the CP-Euclidean +/// (Bobenko-Pinkall-Springborn 2010) circle-packing functional. struct CPEuclideanMaps { CPFMapI f_idx; ///< DOF index per face (−1 = pinned) CPEMapD theta_e; ///< intersection angle per edge (default π/2 = orthogonal) @@ -184,11 +189,8 @@ inline double dof_val(int idx, const std::vector& x) noexcept } // namespace cp_detail -// ── Energy ──────────────────────────────────────────────────────────────────── -// -// Mirrors evaluateEnergyAndGradient in the Java code (lines 170-240) for the -// energy accumulation only. The gradient is computed in a dedicated function -// below for clarity. +/// CP-Euclidean energy value at DOF vector `x` (ρ per face). +/// Mirrors `evaluateEnergyAndGradient()` in the Java original (lines 170-240). inline double cp_euclidean_energy(const ConformalMesh& mesh, const std::vector& x, const CPEuclideanMaps& m) @@ -233,10 +235,8 @@ inline double cp_euclidean_energy(const ConformalMesh& mesh, return E; } -// ── Gradient ────────────────────────────────────────────────────────────────── -// -// ∂E/∂ρ_f = φ_f − Σ_{h: face(h)=f, !is_border(h)} (p + θ*) -// OR (boundary): 2 θ* +/// CP-Euclidean gradient `∂E/∂ρ_f` (per face DOF). Interior term +/// `−(p + θ*)`, boundary term `−2 θ*`; see `setup_cp_euclidean_maps`. inline std::vector cp_euclidean_gradient(const ConformalMesh& mesh, const std::vector& x, const CPEuclideanMaps& m) @@ -280,12 +280,10 @@ inline std::vector cp_euclidean_gradient(const ConformalMesh& mes return G; } -// ── Hessian (analytic) ──────────────────────────────────────────────────────── -// -// Per interior undirected edge e with adjacent faces (j, k): -// h_jk = sin θ / (cosh(Δρ) − cos θ) -// Diagonal contributions on both endpoints; off-diagonal block is −h_jk. -// Pinned faces are skipped (their DOF index is −1 ⇒ excluded from the matrix). +/// Analytic CP-Euclidean Hessian, sparse. Per interior edge `(j,k)` +/// the contribution is `h_jk = sin θ / (cosh(Δρ) − cos θ)`, added to +/// diagonals `H_jj`, `H_kk` and subtracted off-diagonals `H_jk = H_kj`. +/// Pinned faces are excluded (DOF index −1). inline Eigen::SparseMatrix cp_euclidean_hessian(const ConformalMesh& mesh, const std::vector& x, const CPEuclideanMaps& m) @@ -323,11 +321,8 @@ inline Eigen::SparseMatrix cp_euclidean_hessian(const ConformalMesh& return H; } -// ── Finite-difference gradient check ───────────────────────────────────────── -// -// Mirrors the Java FunctionalTest pattern: -// For each DOF i, compare analytic G[i] to (E(x+ε·e_i) − E(x−ε·e_i)) / (2ε). -// Default tolerance 1e-6 with ε = 1e-5 leaves ~3 digits of margin for sane meshes. +/// FD gradient check for the CP-Euclidean functional. Mirrors the +/// Java `FunctionalTest`; default `eps = 1e-5`, `tol = 1e-6`. inline bool gradient_check_cp_euclidean(const ConformalMesh& mesh, const std::vector& x, const CPEuclideanMaps& m, @@ -355,10 +350,8 @@ inline bool gradient_check_cp_euclidean(const ConformalMesh& mesh, return true; } -// ── Finite-difference Hessian check ────────────────────────────────────────── -// -// Verifies analytic H against ( G(x+ε·e_i) − G(x−ε·e_i) ) / (2ε) column-wise. -// Symmetry is implicit in the analytic form; we check both off-diagonal entries. +/// FD Hessian check for the CP-Euclidean functional. Verifies analytic +/// `H` column-by-column against `(G(x+εe_j) − G(x−εe_j)) / (2ε)`. inline bool hessian_check_cp_euclidean(const ConformalMesh& mesh, const std::vector& x, const CPEuclideanMaps& m, diff --git a/code/include/cut_graph.hpp b/code/include/cut_graph.hpp index 0092ef1..1e802a5 100644 --- a/code/include/cut_graph.hpp +++ b/code/include/cut_graph.hpp @@ -36,6 +36,8 @@ namespace conformallab { // CutGraph // ───────────────────────────────────────────────────────────────────────────── +/// Cut-graph result of the tree-cotree algorithm: the set of `2g` edges +/// whose removal turns a closed genus-`g` surface into a topological disk. struct CutGraph { /// cut_edge_flags[e.idx()] = true ↔ this edge is a cut edge. /// Size = mesh.number_of_edges(). @@ -47,6 +49,7 @@ struct CutGraph { /// Genus of the surface (0 for topological spheres and open patches). int genus = 0; + /// `true` iff edge `e` is a cut edge of this graph. bool is_cut(Edge_index e) const { return static_cast(e.idx()) < cut_edge_flags.size() @@ -54,17 +57,10 @@ struct CutGraph { } }; -// ───────────────────────────────────────────────────────────────────────────── -// compute_cut_graph -// ───────────────────────────────────────────────────────────────────────────── -// -// Implements the standard tree-cotree algorithm (Erickson–Whittlesey 2005): -// -// Step 1: BFS primal spanning tree T (V−1 primal tree edges). -// Step 2: BFS dual spanning tree T* (F−1 dual/primal edges, avoiding -// edges whose primal crosses T). -// Step 3: cut edges = primal edges neither in T nor "used" by T*. - +/// Compute the cut graph of `mesh` via the standard tree-cotree +/// algorithm (Erickson–Whittlesey 2005): primal BFS spanning tree T, +/// dual BFS spanning tree T* avoiding T-primals, then the `2g` cut +/// edges are those in neither T nor T*. inline CutGraph compute_cut_graph(const ConformalMesh& mesh) { const std::size_t nv = mesh.number_of_vertices(); diff --git a/code/include/discrete_elliptic_utility.hpp b/code/include/discrete_elliptic_utility.hpp index c49dd0d..4622968 100644 --- a/code/include/discrete_elliptic_utility.hpp +++ b/code/include/discrete_elliptic_utility.hpp @@ -17,7 +17,9 @@ namespace conformallab { // 3. Re-flip: Re < 0 → Re = -Re // 4. S-invert: |tau| < 1 → tau = 1/tau // -// Corresponds to Java DiscreteEllipticUtility.normalizeModulus(Complex). +/// Normalise a complex modulus `τ` into the standard fundamental +/// domain of an elliptic curve (`|τ| ≥ 1`, `0 ≤ Re τ ≤ ½`, `Im τ ≥ 0`). +/// Same as Java `DiscreteEllipticUtility.normalizeModulus(Complex)`. inline std::complex normalizeModulus(std::complex tau) { int maxIter = 100; while (--maxIter > 0) { diff --git a/code/include/euclidean_functional.hpp b/code/include/euclidean_functional.hpp index 1d9fc86..638709c 100644 --- a/code/include/euclidean_functional.hpp +++ b/code/include/euclidean_functional.hpp @@ -48,13 +48,18 @@ namespace conformallab { // ── Property-map type aliases ───────────────────────────────────────────────── +/// Property map vertex → `double` for the Euclidean functional. using EuclVMapD = ConformalMesh::Property_map; +/// Property map vertex → `int` for the Euclidean functional. using EuclVMapI = ConformalMesh::Property_map; +/// Property map edge → `double` for the Euclidean functional. using EuclEMapD = ConformalMesh::Property_map; +/// Property map edge → `int` for the Euclidean functional. using EuclEMapI = ConformalMesh::Property_map; // ── Persistent map bundle ───────────────────────────────────────────────────── +/// Bundle of the five property maps consumed by the Euclidean functional. 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) @@ -119,10 +124,9 @@ inline int euclidean_dimension(const ConformalMesh& mesh, const EuclideanMaps& m 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. +/// Set `lambda0` from mesh vertex positions: +/// `λ°_e = 2·log(|p_i − p_j|)` (natural log of Euclidean edge length²). +/// 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()) { @@ -142,25 +146,23 @@ inline void compute_euclidean_lambda0_from_mesh(ConformalMesh& mesh, EuclideanMa // ── Internal helpers ────────────────────────────────────────────────────────── +/// Read DOF value from `x` for index `idx`; return 0 if pinned (idx < 0). static inline double eucl_dof_val(int idx, const std::vector& x) { return idx >= 0 ? x[static_cast(idx)] : 0.0; } +/// Convert a CGAL half-edge index to a plain `std::size_t` for vector indexing. 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) +/// Compute the Euclidean-functional gradient G(x): +/// * `G_v = Θ_v − Σ_faces α_v(face)` +/// * `G_e = α_opp(face⁺) + α_opp(face⁻) − φ_e` +/// +/// Same half-edge corner-angle storage convention as `spherical_gradient`. inline std::vector euclidean_gradient( ConformalMesh& mesh, const std::vector& x, @@ -238,9 +240,8 @@ inline std::vector euclidean_gradient( return G; } -// ── Energy via Gauss-Legendre path integral ─────────────────────────────────── -// -// E(x) = ∫₀¹ ⟨G(tx), x⟩ dt (10-point GL quadrature, same as SphericalFunctional) +/// Euclidean energy `E(x) = ∫₀¹ ⟨G(t·x), x⟩ dt`, evaluated with +/// 10-point Gauss-Legendre quadrature (same as the Spherical functional). inline double euclidean_energy( ConformalMesh& mesh, const std::vector& x, @@ -282,11 +283,14 @@ inline double euclidean_energy( // ── Full evaluation (energy + gradient) ────────────────────────────────────── +/// Output of `evaluate_euclidean()` — energy plus optional gradient. struct EuclideanResult { - double energy = 0.0; - std::vector gradient; + double energy = 0.0; ///< Functional value at input DOFs. + std::vector gradient; ///< Gradient ∇E (empty if not requested). }; +/// Evaluate the Euclidean functional at DOFs `x`. Returns energy and +/// gradient (toggle via `need_energy` / `need_gradient`). inline EuclideanResult evaluate_euclidean( ConformalMesh& mesh, const std::vector& x, @@ -302,10 +306,8 @@ inline EuclideanResult evaluate_euclidean( 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. +/// Finite-difference gradient check for the Euclidean functional +/// (central differences). Defaults `eps = 1e-5`, `tol = 1e-4`. inline bool gradient_check_euclidean( ConformalMesh& mesh, const std::vector& x0, diff --git a/code/include/euclidean_geometry.hpp b/code/include/euclidean_geometry.hpp index da636da..0b6e072 100644 --- a/code/include/euclidean_geometry.hpp +++ b/code/include/euclidean_geometry.hpp @@ -29,19 +29,17 @@ namespace conformallab { +/// Interior corner angles of a Euclidean triangle. 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; + double alpha1; ///< Corner angle at v₁ (opposite l₂₃). + double alpha2; ///< Corner angle at v₂ (opposite l₃₁). + double alpha3; ///< Corner angle at v₃ (opposite l₁₂). + bool valid; ///< `false` when the triangle is degenerate. }; -// ── 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). +/// Compute the corner angles of a Euclidean triangle from its three +/// side lengths. Returns `valid = false` when the triangle inequality +/// is violated. inline EuclideanFaceAngles euclidean_angles_from_lengths( double l12, double l23, double l31) { @@ -70,14 +68,9 @@ inline EuclideanFaceAngles euclidean_angles_from_lengths( }; } -// ── 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 |Λ̃|. +/// Compute the corner angles of a Euclidean triangle from its three +/// effective log-lengths `Λ̃ᵢⱼ`. Internally centres lengths so that +/// `l₁₂·l₂₃·l₃₁ = 1` to avoid float overflow for large `|Λ̃|`. inline EuclideanFaceAngles euclidean_angles( double lam12, double lam23, double lam31) { diff --git a/code/include/euclidean_hessian.hpp b/code/include/euclidean_hessian.hpp index 3f7cb49..ff76669 100644 --- a/code/include/euclidean_hessian.hpp +++ b/code/include/euclidean_hessian.hpp @@ -52,8 +52,18 @@ namespace conformallab { // cot_k = (t_adj·l123 − t_opp·t_other) / (8·Area) // // Returns {0,0,0} for degenerate faces (triangle inequality violated or Area=0). -struct EuclCotWeights { double cot1, cot2, cot3; bool valid; }; +/// Three Euclidean cotangent weights `(cot1, cot2, cot3)` for the +/// vertices opposite to edges (l₂₃, l₃₁, l₁₂) of a triangle, plus a +/// `valid` flag that is `false` when the triangle is degenerate. +struct EuclCotWeights { + double cot1; ///< Cotangent at vertex 1 (opposite to l₂₃). + double cot2; ///< Cotangent at vertex 2 (opposite to l₃₁). + double cot3; ///< Cotangent at vertex 3 (opposite to l₁₂). + bool valid;///< `false` when the triangle is degenerate (triangle inequality violated or area = 0). +}; +/// Compute the three Euclidean cotangent weights from edge lengths. +/// Returns `{0,0,0,false}` for degenerate triangles. inline EuclCotWeights euclidean_cot_weights(double l12, double l23, double l31) { const double t12 = -l12 + l23 + l31; @@ -81,15 +91,9 @@ inline EuclCotWeights euclidean_cot_weights(double l12, double l23, double l31) }; } -// ── Analytical Hessian (cotangent Laplacian) ────────────────────────────────── -// -// Returns the n×n sparse Hessian matrix H where n = euclidean_dimension(mesh, m). -// -// Only vertex DOFs are supported. Edge DOFs (m.e_idx[e] >= 0) produce -// additional mixed-derivative entries that are not yet implemented; this -// function asserts they are absent. -// -// x – current DOF vector (used to compute effective log-lengths Λ̃ij). +/// Analytical Euclidean Hessian (cotangent Laplacian), sparse. +/// Only vertex DOFs are supported — the function asserts that no edge +/// DOF is variable. `x` is used to compute effective log-lengths Λ̃ᵢⱼ. inline Eigen::SparseMatrix euclidean_hessian( ConformalMesh& mesh, const std::vector& x, @@ -168,11 +172,9 @@ inline Eigen::SparseMatrix euclidean_hessian( } // ── Finite-difference Hessian check ────────────────────────────────────────── -// -// Compares the analytical Hessian column-by-column against -// H_fd[:, j] = (G(x + ε·eⱼ) − G(x − ε·eⱼ)) / (2ε). -// -// Returns true if max relative error < tol for every entry. +/// FD Hessian check for the Euclidean functional. Compares analytic +/// `H` column-by-column to `(G(x+εeⱼ) − G(x−εeⱼ)) / (2ε)`; returns +/// `true` iff max relative error is below `tol`. inline bool hessian_check_euclidean( ConformalMesh& mesh, const std::vector& x0, diff --git a/code/include/fundamental_domain.hpp b/code/include/fundamental_domain.hpp index 0ea180a..e3c1b8f 100644 --- a/code/include/fundamental_domain.hpp +++ b/code/include/fundamental_domain.hpp @@ -43,6 +43,9 @@ namespace conformallab { // FundamentalDomain // ───────────────────────────────────────────────────────────────────────────── +/// Fundamental polygon of a closed surface obtained by cutting along a +/// `CutGraph`: corner vertices, paired-edge identifications and holonomy +/// generators. For genus-1 the polygon is a parallelogram with 4 corners. struct FundamentalDomain { /// Polygon corners in order (CCW). Size = 4 for genus-1. std::vector vertices; @@ -56,6 +59,7 @@ struct FundamentalDomain { /// For genus-1: generators[0] = ω_1, generators[1] = ω_2. std::vector generators; + /// `true` iff the polygon has at least 3 vertices. bool is_valid() const { return vertices.size() >= 3; } }; @@ -75,6 +79,8 @@ struct FundamentalDomain { // bottom (v0→v1) ≡ top (v3→v2) by ω_2 // left (v3→v0) ≡ right (v2→v1) by ω_1 (reversed convention) // ───────────────────────────────────────────────────────────────────────────── +/// Build the parallelogram fundamental domain from genus-1 Euclidean +/// holonomy data (`hol.translations[0] = ω₁`, `hol.translations[1] = ω₂`). inline FundamentalDomain compute_fundamental_domain_genus1( const HolonomyData& hol) { @@ -148,6 +154,8 @@ inline FundamentalDomain compute_fundamental_domain_genus1( // this is intentionally deferred and NOT implemented here. // See period_matrix.hpp for the genus-1 case (τ = ω_2/ω_1 ∈ ℍ). // ───────────────────────────────────────────────────────────────────────────── +/// Dispatcher: for genus 1 returns `compute_fundamental_domain_genus1`, +/// for higher genus returns an empty domain (4g-polygon not yet implemented). inline FundamentalDomain compute_fundamental_domain( const HolonomyData& hol) { @@ -165,6 +173,8 @@ inline FundamentalDomain compute_fundamental_domain( // return a translated copy of the layout shifted by m·ω_1 + n·ω_2. // Useful for visualising the tiled universal cover. // ───────────────────────────────────────────────────────────────────────────── +/// Return a translated copy of `layout` shifted by `m·ω₁ + n·ω₂`. +/// Useful for visualising the tiled universal cover. inline Layout2D tiling_copy(const Layout2D& layout, const Eigen::Vector2d& w1, const Eigen::Vector2d& w2, @@ -183,6 +193,8 @@ inline Layout2D tiling_copy(const Layout2D& layout, // Returns a vector of tiling copies for (m, n) with |m| ≤ m_max, |n| ≤ n_max. // The result includes the original (m=0, n=0) at index (m_max)(2*n_max+1)+n_max. // ───────────────────────────────────────────────────────────────────────────── +/// Build a tiling neighbourhood: all `(m, n)` with `|m| ≤ m_max`, +/// `|n| ≤ n_max`. The original tile `(0, 0)` is included. inline std::vector tiling_neighbourhood( const Layout2D& layout, const HolonomyData& hol, diff --git a/code/include/gauss_bonnet.hpp b/code/include/gauss_bonnet.hpp index 0d5bc61..eed2daf 100644 --- a/code/include/gauss_bonnet.hpp +++ b/code/include/gauss_bonnet.hpp @@ -56,6 +56,7 @@ inline int genus(const ConformalMesh& mesh) // ── Left-hand side Σ(2π − Θ_v) ───────────────────────────────────────────── +/// Sum `Σ_v (2π − Θ_v)` for a raw vertex → angle property map. inline double gauss_bonnet_sum( const ConformalMesh& mesh, const ConformalMesh::Property_map& theta) @@ -66,15 +67,19 @@ inline double gauss_bonnet_sum( return s; } +/// `gauss_bonnet_sum` for the Euclidean-functional property bundle. inline double gauss_bonnet_sum(const ConformalMesh& m, const EuclideanMaps& mp) { return gauss_bonnet_sum(m, mp.theta_v); } +/// `gauss_bonnet_sum` for the Spherical-functional property bundle. inline double gauss_bonnet_sum(const ConformalMesh& m, const SphericalMaps& mp) { return gauss_bonnet_sum(m, mp.theta_v); } +/// `gauss_bonnet_sum` for the HyperIdeal-functional property bundle. inline double gauss_bonnet_sum(const ConformalMesh& m, const HyperIdealMaps& mp) { return gauss_bonnet_sum(m, mp.theta_v); } // ── Right-hand side 2π · χ(M) ─────────────────────────────────────────────── +/// Right-hand side of Gauss-Bonnet: `2π · χ(M)`. inline double gauss_bonnet_rhs(const ConformalMesh& mesh) { return TWO_PI * static_cast(euler_characteristic(mesh)); @@ -82,14 +87,15 @@ inline double gauss_bonnet_rhs(const ConformalMesh& mesh) // ── Deficit: lhs − rhs (0 = Gauss–Bonnet satisfied) ───────────────────────── +/// Gauss-Bonnet deficit `lhs − rhs`; zero iff the identity is satisfied. template inline double gauss_bonnet_deficit(const ConformalMesh& mesh, const Maps& maps) { return gauss_bonnet_sum(mesh, maps) - gauss_bonnet_rhs(mesh); } -// ── check_gauss_bonnet — throws std::runtime_error if |deficit| > tol ───────── - +/// Throws `std::runtime_error` if `|lhs − 2π·χ| > tol`. +/// Overload accepting a precomputed `lhs`. inline void check_gauss_bonnet(const ConformalMesh& mesh, double lhs, double tol = 1e-8) @@ -108,6 +114,7 @@ inline void check_gauss_bonnet(const ConformalMesh& mesh, } } +/// Throws `std::runtime_error` if Gauss-Bonnet is violated by more than `tol`. template inline void check_gauss_bonnet(const ConformalMesh& mesh, const Maps& maps, @@ -123,6 +130,9 @@ inline void check_gauss_bonnet(const ConformalMesh& mesh, // Only modifies free vertices (v_idx[v] >= 0 for EuclideanMaps / SphericalMaps; // always all vertices for the raw property-map overload). +/// Distribute the Gauss-Bonnet deficit uniformly across all `Θ_v`: +/// add `δ = (lhs − rhs) / V` to every entry so that the identity holds +/// exactly afterwards. Overload for a raw property map. inline void enforce_gauss_bonnet( ConformalMesh& mesh, ConformalMesh::Property_map& theta) @@ -136,6 +146,7 @@ inline void enforce_gauss_bonnet( theta[v] += delta; } +/// Distribute the Gauss-Bonnet deficit uniformly across `maps.theta_v`. template inline void enforce_gauss_bonnet(ConformalMesh& mesh, Maps& maps) { diff --git a/code/include/hyper_ideal_functional.hpp b/code/include/hyper_ideal_functional.hpp index 2a50868..a449dd5 100644 --- a/code/include/hyper_ideal_functional.hpp +++ b/code/include/hyper_ideal_functional.hpp @@ -39,18 +39,23 @@ namespace conformallab { // ── Property-map type aliases ───────────────────────────────────────────────── +/// Property map vertex → `double` (HyperIdeal scalar-per-vertex data). using VMapD = ConformalMesh::Property_map; +/// Property map vertex → `int` (HyperIdeal DOF indices). using VMapI = ConformalMesh::Property_map; +/// Property map edge → `double` (HyperIdeal scalar-per-edge data). using EMapD = ConformalMesh::Property_map; +/// Property map edge → `int` (HyperIdeal DOF indices). using EMapI = ConformalMesh::Property_map; // ── Persistent map bundle ───────────────────────────────────────────────────── +/// Bundle of the four property maps consumed by the HyperIdeal functional. 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 + 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 Θᵥ (parameter, not variable). + EMapD theta_e; ///< Target intersection angle θₑ. }; /// Attach the four HyperIdeal property maps to `mesh` and return their @@ -105,20 +110,22 @@ inline int assign_all_dof_indices(ConformalMesh& mesh, HyperIdealMaps& m) // ── Evaluation result ───────────────────────────────────────────────────────── +/// Output of `evaluate_hyper_ideal()` — the energy value and (optionally) +/// its gradient evaluated at the current DOF vector. struct HyperIdealResult { - double energy = 0.0; - std::vector gradient; // empty when gradient was not requested + double energy = 0.0; ///< Functional value at the input DOFs. + std::vector gradient; ///< Gradient ∇E; empty when not requested. }; // ── Internal helpers ────────────────────────────────────────────────────────── -// Get the DOF value from x, or 0.0 if pinned. +/// Read the DOF value from `x` for index `idx`; return 0 if pinned (idx < 0). static inline double dof_val(int idx, const std::vector& x) { return idx >= 0 ? x[static_cast(idx)] : 0.0; } -// Convert a CGAL halfedge index to a plain std::size_t (for vector indexing). +/// Convert a CGAL half-edge index to a plain `std::size_t` for vector indexing. static inline std::size_t hidx(Halfedge_index h) { return static_cast(static_cast(h)); @@ -144,11 +151,21 @@ static inline std::size_t hidx(Halfedge_index h) // HyperIdealFunctional.java's defensive behaviour (lines 122-127 of the // Java original); this keeps the FD perturbation regime well-defined. +/// Six per-face angle outputs computed from local DOFs (see +/// `face_angles_from_local_dofs`). Used by the block-FD Hessian. struct FaceAngleOutputs { - double beta1, beta2, beta3; ///< interior angles at v₁,v₂,v₃ - double alpha12, alpha23, alpha31; ///< dihedral angles at e₁₂,e₂₃,e₃₁ + double beta1; ///< Interior angle at v₁. + double beta2; ///< Interior angle at v₂. + double beta3; ///< Interior angle at v₃. + double alpha12; ///< Dihedral angle at edge e₁₂. + double alpha23; ///< Dihedral angle at edge e₂₃. + double alpha31; ///< Dihedral angle at edge e₃₁. }; +/// Pure-math 6→6 kernel: given the six local DOFs (b₁,b₂,b₃,a₁₂,a₂₃,a₃₁) +/// of one face plus the per-vertex variability flags, return the six +/// HyperIdeal angle outputs. No mesh, no property maps — used by the +/// per-face block-FD Hessian in `hyper_ideal_hessian.hpp`. inline FaceAngleOutputs face_angles_from_local_dofs( double b1, double b2, double b3, double a12, double a23, double a31, @@ -196,14 +213,29 @@ inline FaceAngleOutputs face_angles_from_local_dofs( // ── Per-face angle kernel ───────────────────────────────────────────────────── +/// Per-face angle bundle returned by `compute_face_angles()`. Carries +/// the six output angles plus the six input DOFs (so the energy and +/// gradient kernels can reuse them without re-reading the mesh). 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 + double alpha12; ///< Dihedral angle at edge e₁₂. + double alpha23; ///< Dihedral angle at edge e₂₃. + double alpha31; ///< Dihedral angle at edge e₃₁. + double beta1; ///< Interior angle at vertex v₁. + double beta2; ///< Interior angle at vertex v₂. + double beta3; ///< Interior angle at vertex v₃. + double a12; ///< Edge DOF value at e₁₂. + double a23; ///< Edge DOF value at e₂₃. + double a31; ///< Edge DOF value at e₃₁. + double b1; ///< Vertex DOF value at v₁. + double b2; ///< Vertex DOF value at v₂. + double b3; ///< Vertex DOF value at v₃. + bool v1b; ///< `true` iff vertex v₁ is variable (not pinned). + bool v2b; ///< `true` iff vertex v₂ is variable. + bool v3b; ///< `true` iff vertex v₃ is variable. }; +/// Compute the six per-face angles (+ remember the input DOFs) for face +/// `f` of `mesh`, given the current DOF vector `x` and DOF-index maps. static FaceAngles compute_face_angles( const ConformalMesh& mesh, Face_index f, @@ -281,7 +313,7 @@ static FaceAngles compute_face_angles( return fa; } -// Per-face energy contribution U(f) (before subtracting θ·a and Θ·b terms). +/// Per-face energy contribution U(f) before subtracting the θ·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; @@ -310,6 +342,14 @@ static double face_energy(const FaceAngles& fa) // ── Full evaluation ─────────────────────────────────────────────────────────── +/// Evaluate the HyperIdeal functional at DOF vector `x`. Returns the +/// energy value and (optionally) the gradient in a `HyperIdealResult`. +/// +/// \param mesh Triangle mesh carrying the DOF-index property maps. +/// \param x Current DOF vector (length = `hyper_ideal_dimension(...)`). +/// \param m Property-map bundle from `setup_hyper_ideal_maps(...)`. +/// \param need_energy If `true`, fill `result.energy` (default: `true`). +/// \param need_gradient If `true`, fill `result.gradient` (default: `true`). inline HyperIdealResult evaluate_hyper_ideal( ConformalMesh& mesh, const std::vector& x, @@ -393,10 +433,10 @@ inline HyperIdealResult evaluate_hyper_ideal( 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). +/// Finite-difference gradient check (central differences). +/// +/// Returns `true` iff `|G[i] − fd[i]| / max(1, |G[i]|) < tol` for every +/// DOF. Defaults `eps = 1e-5`, `tol = 1e-4` match the Java `FunctionalTest`. inline bool gradient_check( ConformalMesh& mesh, const std::vector& x0, diff --git a/code/include/hyper_ideal_geometry.hpp b/code/include/hyper_ideal_geometry.hpp index 0ac4651..34f8b6a 100644 --- a/code/include/hyper_ideal_geometry.hpp +++ b/code/include/hyper_ideal_geometry.hpp @@ -22,9 +22,9 @@ namespace conformallab { // ── 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). +/// `ζ(x,y,z)` — interior angle (in radians) 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); @@ -34,8 +34,8 @@ inline double zeta(double x, double y, double z) return std::acos(nbd); } -// ζ₁₃(x,y,z) — third edge length in a right-angled hyperbolic hexagon. -// Ports HyperIdealUtility.ζ_13(x, y, z). +/// `ζ₁₃(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); @@ -43,16 +43,16 @@ inline double zeta13(double x, double y, double z) 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). +/// `ζ₁₄(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). +/// `ζ₁₅(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)); @@ -60,12 +60,11 @@ inline double zeta15(double x) // ── 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(). +/// `l_ij`: effective hyperbolic length of edge ij. +/// * `bi`, `bj` — vertex log scale factors (used only when vertex is hyper-ideal). +/// * `aij` — edge intersection-angle variable. +/// * `vi_var` / `vj_var` — `true` iff the corresponding vertex 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); @@ -76,8 +75,8 @@ inline double lij(double bi, double bj, double aij, bool vi_var, bool vj_var) // ── Auxiliary angle functions ───────────────────────────────────────────────── -// σᵢ(aᵢⱼ, aₖᵢ, aⱼₖ, vj_var, vk_var) — intermediate half-length at vertex i. -// Ports HyperIdealFunctional.σi(). +/// `σᵢ(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); @@ -86,24 +85,24 @@ inline double sigma_i(double aij, double aki, double ajk, bool vj_var, bool vk_v return zeta15(ajk - aij - aki); } -// σᵢⱼ(aᵢⱼ, bᵢ, bⱼ, vj_var) — intermediate half-length for edge ij from vertex i. -// Ports HyperIdealFunctional.σij(). +/// `σᵢⱼ(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). +/// `α_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. +/// * `beta_i, beta_j, beta_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, diff --git a/code/include/hyper_ideal_hessian.hpp b/code/include/hyper_ideal_hessian.hpp index 9843a41..b809b43 100644 --- a/code/include/hyper_ideal_hessian.hpp +++ b/code/include/hyper_ideal_hessian.hpp @@ -54,13 +54,9 @@ namespace conformallab { -// ── Full finite-difference Hessian (baseline, Phase 4a) ────────────────────── -// -// Returns the n×n sparse Hessian, where n = hyper_ideal_dimension(mesh, m). -// eps: finite-difference step size (default 1e-5 gives ~1e-10 relative error). -// -// Cost: n × full-gradient evaluations ≈ O(n·F). Use for small meshes or -// as a correctness reference for the block-FD variant below. +/// Full finite-difference HyperIdeal Hessian (baseline, Phase 4a). +/// Cost: `n` full-gradient evaluations ≈ `O(n·F)`. Use for small +/// meshes or as a correctness reference for the block-FD variant. inline Eigen::SparseMatrix hyper_ideal_hessian( ConformalMesh& mesh, const std::vector& x, @@ -96,10 +92,8 @@ inline Eigen::SparseMatrix hyper_ideal_hessian( return H; } -// ── Symmetrised Hessian (full-FD variant) ──────────────────────────────────── -// -// The FD Hessian is symmetric in exact arithmetic; floating-point rounding -// can introduce tiny asymmetries. This helper returns (H + Hᵀ)/2. +/// Symmetrised full-FD HyperIdeal Hessian: returns `(H + Hᵀ) / 2` to +/// scrub the tiny asymmetries introduced by floating-point rounding. inline Eigen::SparseMatrix hyper_ideal_hessian_sym( ConformalMesh& mesh, const std::vector& x, @@ -137,6 +131,10 @@ inline Eigen::SparseMatrix hyper_ideal_hessian_sym( // vs full-FD ≈ 99,200 → ~33× speed-up. // On brezel.obj (F=13824, n≈14000): 165 888 face evaluations // vs full-FD ≈ 193 M → ~1166× speed-up. +/// Per-face block-FD HyperIdeal Hessian (Phase 9b). Uses the locality +/// lemma `∂G_x/∂y = Σ_{f: x,y ∈ local(f)} ∂(β or α)/∂y` to perturb only +/// the 6 face-local DOFs at a time, giving an `F·12` face-evaluation +/// budget vs `n·F` for full-FD (~96× speed-up on brezel.obj). inline Eigen::SparseMatrix hyper_ideal_hessian_block_fd( ConformalMesh& mesh, const std::vector& x, @@ -213,11 +211,9 @@ inline Eigen::SparseMatrix hyper_ideal_hessian_block_fd( return H; } -// ── Symmetrised block-FD Hessian ───────────────────────────────────────────── -// -// The per-face block is symmetric to within FD rounding, but the -// accumulation may amplify tiny asymmetries. This helper returns -// (H + Hᵀ)/2, identical in spirit to `hyper_ideal_hessian_sym`. +/// Symmetrised block-FD HyperIdeal Hessian: returns `(H + Hᵀ) / 2` of +/// `hyper_ideal_hessian_block_fd(...)` for downstream solvers that +/// require strict symmetry. inline Eigen::SparseMatrix hyper_ideal_hessian_block_fd_sym( ConformalMesh& mesh, const std::vector& x, diff --git a/code/include/hyper_ideal_utility.hpp b/code/include/hyper_ideal_utility.hpp index 0734731..bbaed62 100644 --- a/code/include/hyper_ideal_utility.hpp +++ b/code/include/hyper_ideal_utility.hpp @@ -12,9 +12,9 @@ namespace conformallab { -// Volume of a generalized hyperbolic tetrahedron with dihedral angles A..F. -// Formula: Meyerhoff / Ushijima (Springer 2006). -// Corresponds to Java HyperIdealUtility.calculateTetrahedronVolume(). +/// Volume of a generalized hyperbolic tetrahedron with dihedral +/// angles `A,…,F` via the Meyerhoff / Ushijima 2006 formula. +/// Same as Java `HyperIdealUtility.calculateTetrahedronVolume()`. inline double calculateTetrahedronVolume(double A, double B, double C, double D, double E, double F) { // PI from constants.hpp (conformallab::PI) @@ -73,11 +73,9 @@ inline double calculateTetrahedronVolume(double A, double B, double C, return (U(z1) - U(z2)) / 2.0; } -// Volume of a hyperideal tetrahedron with one ideal vertex (at gamma). -// Dihedral angles at the ideal vertex: gamma1, gamma2, gamma3. -// Dihedral angles at opposite edges: alpha23, alpha31, alpha12. -// Formula: Kolpakov–Mednykh (arxiv math/0603097). -// Corresponds to Java HyperIdealUtility.calculateTetrahedronVolumeWithIdealVertexAtGamma(). +/// Volume of a hyperideal tetrahedron with one ideal vertex at γ via +/// the Kolpakov-Mednykh formula (arxiv math/0603097). Same as Java +/// `HyperIdealUtility.calculateTetrahedronVolumeWithIdealVertexAtGamma()`. inline double calculateTetrahedronVolumeWithIdealVertexAtGamma( double gamma1, double gamma2, double gamma3, double alpha23, double alpha31, double alpha12) diff --git a/code/include/hyper_ideal_visualization_utility.hpp b/code/include/hyper_ideal_visualization_utility.hpp index f5ec00a..f16394c 100644 --- a/code/include/hyper_ideal_visualization_utility.hpp +++ b/code/include/hyper_ideal_visualization_utility.hpp @@ -32,9 +32,7 @@ namespace conformallab { -// --------------------------------------------------------------------------- -// Circumcenter of three 2-D points -// --------------------------------------------------------------------------- +/// Circumcenter of three 2-D points (`a`, `b`, `c`) in the Euclidean plane. inline Eigen::Vector2d circumcenter2d( const Eigen::Vector2d& a, const Eigen::Vector2d& b, @@ -54,10 +52,8 @@ inline Eigen::Vector2d circumcenter2d( return {ux, uy}; } -// --------------------------------------------------------------------------- -// 4×4 Lorentz boost: maps the hyperboloid origin e₄=(0,0,0,1) to `center`. -// `center` must lie on the hyperboloid: center[3]² - ‖center.head<3>()‖² = 1. -// --------------------------------------------------------------------------- +/// 4×4 Lorentz boost: maps the hyperboloid origin `e₄ = (0,0,0,1)` to +/// `center`. Precondition: `center` lies on the hyperboloid. inline Eigen::Matrix4d hyperboloidTranslation(const Eigen::Vector4d& center) { Eigen::Vector3d p = center.head<3>(); @@ -73,10 +69,8 @@ inline Eigen::Matrix4d hyperboloidTranslation(const Eigen::Vector4d& center) return T; } -// --------------------------------------------------------------------------- -// Project a hyperboloid point to the Poincaré disk (jReality convention: -// add 1 to the w-coordinate, then dehomogenize the spatial part). -// --------------------------------------------------------------------------- +/// Project a hyperboloid point `x` onto the Poincaré disk (jReality +/// convention: add 1 to the w-coordinate, then dehomogenise spatial part). inline Eigen::Vector2d toPoincareDisk(const Eigen::Vector4d& x) { double w = x(3) + 1.0; @@ -97,6 +91,10 @@ inline Eigen::Vector2d toPoincareDisk(const Eigen::Vector4d& x) // // Port of HyperIdealVisualizationPlugin.getEuclideanCircleFromHyperbolic() // --------------------------------------------------------------------------- +/// Convert a hyperbolic circle (`center` on the hyperboloid, hyperbolic +/// `radius`) to the corresponding Euclidean circle in the Poincaré disk; +/// returns `{cx, cy, r}`. Port of `HyperIdealVisualizationPlugin +/// .getEuclideanCircleFromHyperbolic()`. inline std::array getEuclideanCircleFromHyperbolic( const Eigen::Vector4d& center, double radius) { diff --git a/code/include/inversive_distance_functional.hpp b/code/include/inversive_distance_functional.hpp index a932fa5..dd2ad5d 100644 --- a/code/include/inversive_distance_functional.hpp +++ b/code/include/inversive_distance_functional.hpp @@ -76,12 +76,17 @@ namespace conformallab { // ── Property-map type aliases ──────────────────────────────────────────────── +/// Property map vertex → `int` for the Inversive-Distance functional. using IDVMapI = ConformalMesh::Property_map; +/// Property map vertex → `double` for the Inversive-Distance functional. using IDVMapD = ConformalMesh::Property_map; +/// Property map edge → `double` for the Inversive-Distance functional. using IDEMapD = ConformalMesh::Property_map; // ── Persistent map bundle ───────────────────────────────────────────────────── +/// Bundle of the four property maps consumed by the Inversive-Distance +/// circle-packing functional (Luo 2004 / Bowers-Stephenson 2004). struct InversiveDistanceMaps { IDVMapI v_idx; ///< DOF index per vertex (−1 = pinned / u_v = 0) IDVMapD theta_v; ///< target cone angle Θ_v (default 2π) @@ -225,12 +230,8 @@ inline double edge_length_squared(double u_i, double u_j, double I_ij) noexcept } // namespace id_detail -// ── Gradient ────────────────────────────────────────────────────────────────── -// -// G_v = Θ_v − Σ_{f ∋ v} α_v(f) -// -// Halfedge convention (identical to euclidean_functional.hpp): -// h_alpha[h] = corner angle OPPOSITE to the edge of halfedge h in its face. +/// Inversive-Distance gradient `G_v = Θ_v − Σ_faces α_v(face)`. Same +/// half-edge corner-angle storage convention as `euclidean_gradient`. inline std::vector inversive_distance_gradient( const ConformalMesh& mesh, const std::vector& x, @@ -291,9 +292,8 @@ inline std::vector inversive_distance_gradient( return G; } -// ── Energy via Gauss-Legendre path integral ───────────────────────────────── -// -// E(u) = ∫₀¹ ⟨G(tu), u⟩ dt (10-point GL, identical constants to euclidean_functional) +/// Inversive-Distance energy `E(u) = ∫₀¹ ⟨G(t·u), u⟩ dt`, evaluated +/// with 10-point Gauss-Legendre (constants shared with `euclidean_energy`). inline double inversive_distance_energy( const ConformalMesh& mesh, const std::vector& x, @@ -329,7 +329,7 @@ inline double inversive_distance_energy( return E; } -// ── Finite-difference gradient check ───────────────────────────────────────── +/// FD gradient check for the Inversive-Distance functional (central diff). inline bool gradient_check_inversive_distance( const ConformalMesh& mesh, const std::vector& x, @@ -358,7 +358,8 @@ inline bool gradient_check_inversive_distance( return true; } -// ── Newton equilibrium check: gradient vanishes at converged u ────────────── +/// Newton equilibrium check: returns `true` iff the gradient at `x` +/// is below `tol` in infinity norm (Σ adj-face angles equal Θ_v). inline bool is_inversive_distance_equilibrium( const ConformalMesh& mesh, const std::vector& x, diff --git a/code/include/layout.hpp b/code/include/layout.hpp index 1ee2fab..90de1fb 100644 --- a/code/include/layout.hpp +++ b/code/include/layout.hpp @@ -75,28 +75,38 @@ namespace conformallab { // For hyperbolic holonomy the map is an orientation-preserving isometry of // the Poincaré disk (SU(1,1) element). // ───────────────────────────────────────────────────────────────────────────── +/// Möbius transformation `T(z) = (a·z + b) / (c·z + d)` of the Riemann +/// sphere; restricted to SU(1,1) for hyperbolic holonomy on the +/// Poincaré disk. struct MobiusMap { + /// Complex scalar used for all entries. using C = std::complex; - C a{1.0, 0.0}; - C b{0.0, 0.0}; - C c{0.0, 0.0}; - C d{1.0, 0.0}; + C a{1.0, 0.0}; ///< Top-left coefficient. + C b{0.0, 0.0}; ///< Top-right coefficient. + C c{0.0, 0.0}; ///< Bottom-left coefficient. + C d{1.0, 0.0}; ///< Bottom-right coefficient. + /// Apply the transformation to a complex point. C apply(C z) const { return (a * z + b) / (c * z + d); } + /// Apply the transformation to a 2-D real point (interpreted as `x + iy`). Eigen::Vector2d apply(const Eigen::Vector2d& p) const { C w = apply(C(p.x(), p.y())); return Eigen::Vector2d(w.real(), w.imag()); } + /// Identity map. static MobiusMap identity() { return {C(1), C(0), C(0), C(1)}; } + /// Inverse map. MobiusMap inverse() const { return {d, -b, -c, a}; } + /// Composition `(*this) ∘ T`, i.e. apply `T` first then `*this`. MobiusMap compose(const MobiusMap& T) const { return { a*T.a + b*T.c, a*T.b + b*T.d, c*T.a + d*T.c, c*T.b + d*T.d }; } + /// `true` iff the map is the identity up to tolerance `tol`. bool is_identity(double tol = 1e-9) const { if (std::abs(d) < 1e-14) return false; C a_ = a/d, b_ = b/d, c_ = c/d; @@ -121,6 +131,9 @@ struct MobiusMap { // ── Result types ────────────────────────────────────────────────────────────── +/// Result of a 2-D layout (`euclidean_layout`, `hyper_ideal_layout`): +/// per-vertex UV coordinates plus a per-half-edge UV atlas for seamed +/// textures. struct Layout2D { /// uv[v.idx()] — primary 2-D position (first / shallowest-BFS-depth visit). std::vector uv; @@ -136,14 +149,15 @@ struct Layout2D { /// Size = mesh.number_of_halfedges(). Border halfedges = (0,0). std::vector halfedge_uv; - bool success = false; - bool has_seam = false; ///< true when a vertex was reached via two paths + bool success = false; ///< `true` iff the BFS placed every vertex. + bool has_seam = false; ///< `true` when a vertex was reached via two paths. }; +/// Result of a 3-D layout (`spherical_layout`): per-vertex positions on S². struct Layout3D { - std::vector pos; - bool success = false; - bool has_seam = false; + std::vector pos; ///< Per-vertex spherical positions. + bool success = false; ///< `true` iff the BFS placed every vertex. + bool has_seam = false; ///< `true` when a vertex was reached via two paths. }; /// Per-cut-edge holonomy. @@ -156,9 +170,9 @@ struct Layout3D { /// trilaterated virtual position obtained by continuing the unfolding across /// the cut. struct HolonomyData { - std::vector translations; ///< Euclidean / spherical - std::vector mobius_maps; ///< hyperbolic (Phase 7) - std::vector cut_edge_indices; + std::vector translations; ///< Euclidean / spherical translation per cut edge. + std::vector mobius_maps; ///< Hyperbolic Möbius isometry per cut edge (Phase 7). + std::vector cut_edge_indices; ///< Index (in the cut-graph edge list) of each holonomy entry. }; // ── Internal helpers ────────────────────────────────────────────────────────── @@ -343,7 +357,8 @@ inline void center_poincare_disk_weighted( } // namespace detail -// ── Vertex Voronoi area weights ─────────────────────────────────────────────── +/// Compute per-vertex area weights (sum of 1/3 of each adjacent triangle area). +/// Used by area-weighted layout normalisation routines. inline std::vector compute_vertex_area_weights(const ConformalMesh& mesh) { std::vector w(mesh.number_of_vertices(), 0.0); @@ -357,6 +372,8 @@ inline std::vector compute_vertex_area_weights(const ConformalMesh& mesh // ── Layout normalisation ────────────────────────────────────────────────────── +/// Euclidean canonical normalisation: translate centroid to the origin +/// and rotate the principal axis of the UV cloud onto the x-axis. inline void normalise_euclidean(Layout2D& layout) { if (!layout.success || layout.uv.empty()) return; @@ -388,6 +405,8 @@ inline void normalise_hyperbolic(Layout2D& layout, const ConformalMesh& mesh) // halfedge_uv follows the same Möbius map detail::center_poincare_disk(layout.halfedge_uv); } +/// Hyperbolic canonical normalisation, mesh-free fallback: uniform +/// (unweighted) iterative Möbius centring of the Poincaré disk. inline void normalise_hyperbolic(Layout2D& layout) // fallback without mesh { if (!layout.success || layout.uv.empty()) return; @@ -395,6 +414,9 @@ inline void normalise_hyperbolic(Layout2D& layout) // fallback without mesh detail::center_poincare_disk(layout.halfedge_uv); } +/// Spherical canonical normalisation: rotate the layout so that the +/// per-vertex centroid (projected back to S²) coincides with the north +/// pole (Rodrigues rotation). inline void normalise_spherical(Layout3D& layout) { if (!layout.success || layout.pos.empty()) return; @@ -852,6 +874,8 @@ inline Layout2D hyper_ideal_layout( // ── Convenience: save layout as OFF ────────────────────────────────────────── +/// Write a 2-D layout to disk in OFF format with z = 0. Convenience +/// helper for quickly inspecting the UV result in any OFF viewer. inline void save_layout_off( const std::string& path, ConformalMesh& mesh, const Layout2D& layout) { @@ -865,6 +889,7 @@ inline void save_layout_off( } } +/// Write a 3-D (spherical) layout to disk in OFF format. inline void save_layout_off( const std::string& path, ConformalMesh& mesh, const Layout3D& layout) { diff --git a/code/include/matrix_utility.hpp b/code/include/matrix_utility.hpp index f18aee6..9f75045 100644 --- a/code/include/matrix_utility.hpp +++ b/code/include/matrix_utility.hpp @@ -7,12 +7,10 @@ namespace conformallab { -// Find the 4×4 matrix R that maps source points to target points. -// Each row of `from` / `to` is a homogeneous 4-vector (one point per row). -// Post-condition: R * from.row(i).T == to.row(i).T for all i. -// -// Implementation: R = to^T * (from^T)^{-1} -// Corresponds to Java MatrixUtility.makeMappingMatrix(). +/// Find the 4×4 matrix `R` that maps each row of `from` (homogeneous +/// 4-vector) to the corresponding row of `to`: `R · fromᵀ = toᵀ`. +/// Computed as `R = toᵀ · (fromᵀ)⁻¹`. Same as Java +/// `MatrixUtility.makeMappingMatrix()`. inline Eigen::Matrix4d makeMappingMatrix(const Eigen::Matrix4d& from, const Eigen::Matrix4d& to) { return to.transpose() * from.transpose().inverse(); diff --git a/code/include/mesh_builder.hpp b/code/include/mesh_builder.hpp index ef08784..d728b35 100644 --- a/code/include/mesh_builder.hpp +++ b/code/include/mesh_builder.hpp @@ -22,7 +22,7 @@ namespace conformallab { // | \ // v0 ─ v1 // -// Returns a mesh with 1 face, 3 vertices, 3 edges. +/// Build a single right-angle triangle in the xy-plane (1 face, 3 vertices, 3 edges). // The triangle lies in the xy-plane with a right angle at v0. inline ConformalMesh make_triangle( double x0=0, double y0=0, @@ -39,7 +39,7 @@ inline ConformalMesh make_triangle( // ── Regular tetrahedron ────────────────────────────────────────────────────── // -// 4 vertices, 4 faces, 6 edges. +/// Build a regular tetrahedron (4 vertices, 4 faces, 6 edges; sphere topology). // Euler characteristic: V - E + F = 4 - 6 + 4 = 2 (sphere topology). // Used to test closed-surface traversal. inline ConformalMesh make_tetrahedron() @@ -67,8 +67,8 @@ inline ConformalMesh make_tetrahedron() // | \ | // v0 ─ v1 // -// 4 vertices, 2 faces, 5 edges (1 interior edge v1–v2 shared by both faces). -// Useful for testing edge-interior vs edge-boundary distinction. +/// Build a two-triangle strip (4 vertices, 2 faces, 5 edges; 1 interior edge). +/// Useful for testing interior- vs boundary-edge distinction. inline ConformalMesh make_quad_strip() { ConformalMesh mesh; @@ -85,8 +85,8 @@ inline ConformalMesh make_quad_strip() // ── Regular flat polygon fan ───────────────────────────────────────────────── // -// n triangles sharing a central vertex; forms a disk topology (boundary). -// Used to verify valence-n vertex traversal. +/// Build a regular flat polygon fan: `n` triangles sharing a central +/// vertex, with rim vertices on the unit circle (disk topology). inline ConformalMesh make_fan(int n) { CGAL_precondition(n >= 3); @@ -109,10 +109,9 @@ inline ConformalMesh make_fan(int n) // ── Spherical tetrahedron (vertices on the unit sphere) ─────────────────────── // -// The four vertices of a regular tetrahedron projected onto the unit sphere. -// Starting from (±1,±1,±1), dividing by √3 gives unit-length positions. -// All edge lengths equal arccos(−1/3) ≈ 1.9106 radians. -// Used for SphericalFunctional tests (all four faces are valid spherical triangles). +/// Build a regular tetrahedron with vertices on the unit sphere. +/// All edge lengths equal `arccos(−1/3) ≈ 1.9106 rad`; used by the +/// SphericalFunctional tests. inline ConformalMesh make_spherical_tetrahedron() { ConformalMesh mesh; @@ -133,10 +132,9 @@ inline ConformalMesh make_spherical_tetrahedron() // ── Octahedron face triangle (vertices on the unit sphere) ──────────────────── // -// One face of a regular octahedron: the triangle (1,0,0)→(0,1,0)→(0,0,1). -// All edge lengths equal arccos(0) = π/2. -// The corner angles are all π/2 (right-angled spherical triangle). -// base log-length: λ° = 2·log(sin(π/4)) = 2·log(1/√2) = −log(2) ≈ −0.6931. +/// Build one face of a regular octahedron `(1,0,0)→(0,1,0)→(0,0,1)`: +/// a right-angled spherical triangle with edge length `π/2` and base +/// log-length `λ° = −log 2 ≈ −0.6931`. inline ConformalMesh make_octahedron_face() { ConformalMesh mesh; diff --git a/code/include/mesh_io.hpp b/code/include/mesh_io.hpp index 0ac1fc1..7e50ba8 100644 --- a/code/include/mesh_io.hpp +++ b/code/include/mesh_io.hpp @@ -28,26 +28,22 @@ namespace conformallab { -// ── Read ────────────────────────────────────────────────────────────────────── -// -// Reads a polygon mesh from file into `mesh` (clears any existing content). -// Returns true on success, false on failure. +/// Read a polygon mesh from `filename` into `mesh` (clears existing content). +/// Returns `true` on success, `false` on failure. inline bool read_mesh(const std::string& filename, ConformalMesh& mesh) { mesh.clear(); return CGAL::IO::read_polygon_mesh(filename, mesh); } -// ── Write ───────────────────────────────────────────────────────────────────── -// -// Writes `mesh` to `filename`. Returns true on success. +/// Write `mesh` to `filename`. Returns `true` on success. inline bool write_mesh(const std::string& filename, const ConformalMesh& mesh) { return CGAL::IO::write_polygon_mesh(filename, mesh); } -// ── Convenience: throwing wrappers ──────────────────────────────────────────── - +/// Throwing wrapper around `read_mesh`: returns the mesh by value +/// or throws `std::runtime_error` on read failure. inline ConformalMesh load_mesh(const std::string& filename) { ConformalMesh mesh; @@ -56,6 +52,8 @@ inline ConformalMesh load_mesh(const std::string& filename) return mesh; } +/// Throwing wrapper around `write_mesh`; throws `std::runtime_error` on +/// write failure. inline void save_mesh(const std::string& filename, const ConformalMesh& mesh) { if (!write_mesh(filename, mesh)) diff --git a/code/include/newton_solver.hpp b/code/include/newton_solver.hpp index bf290e1..7fe0c2a 100644 --- a/code/include/newton_solver.hpp +++ b/code/include/newton_solver.hpp @@ -44,11 +44,12 @@ namespace conformallab { // ── Result ──────────────────────────────────────────────────────────────────── +/// Result of `newton_solve(...)` — converged DOF vector + diagnostics. struct NewtonResult { - std::vector x; ///< DOF vector at termination - int iterations; ///< Newton steps taken - double grad_inf_norm;///< max |G_i| at termination - bool converged; ///< true iff grad_inf_norm < tol + std::vector x; ///< DOF vector at termination. + int iterations; ///< Newton steps taken. + double grad_inf_norm;///< max |Gᵢ| at termination. + bool converged; ///< `true` iff `grad_inf_norm < tol`. }; // ── Internal helpers ────────────────────────────────────────────────────────── @@ -96,6 +97,10 @@ inline Eigen::VectorXd solve_with_fallback( // // fallback_used – if non-null, set to true iff SparseQR was invoked // Returns Eigen::VectorXd::Zero(rhs.size()) if both solvers fail. +/// Solve `A·x = rhs` with the same SimplicialLDLT → SparseQR fallback +/// strategy used inside all three Newton solvers. If `fallback_used` +/// is non-null, it is set to `true` iff the SparseQR fallback ran. +/// Returns `Eigen::VectorXd::Zero(rhs.size())` if both solvers fail. inline Eigen::VectorXd solve_linear_system( const Eigen::SparseMatrix& A, const Eigen::VectorXd& rhs, diff --git a/code/include/p2_utility.hpp b/code/include/p2_utility.hpp index d4f88ca..241f8d8 100644 --- a/code/include/p2_utility.hpp +++ b/code/include/p2_utility.hpp @@ -13,9 +13,9 @@ namespace conformallab { // ── Point / line duality ────────────────────────────────────────────────────── -// Intersection of two lines l1, l2 (or line through two points p1, p2) -// via the cross product. Works for any P2 element. -// Corresponds to Java P2.pointFromLines / P2.lineFromPoints. +/// Cross-product point–line duality in P²: returns the intersection +/// of two lines (or the line through two points). Same as Java +/// `P2.pointFromLines` / `P2.lineFromPoints`. inline Eigen::Vector3d pointFromLines(const Eigen::Vector3d& l1, const Eigen::Vector3d& l2) { return l1.cross(l2); @@ -23,11 +23,9 @@ inline Eigen::Vector3d pointFromLines(const Eigen::Vector3d& l1, // ── Euclidean perpendicular bisector ───────────────────────────────────────── -// Returns the homogeneous line coordinates (a, b, c) of the perpendicular -// bisector of the segment [p, q] in the Euclidean plane. -// Coordinates: ax + by + c = 0 (after dehomogenizing p and q). -// -// Corresponds to Java P2.perpendicularBisector(p, q, Pn.EUCLIDEAN). +/// Homogeneous line coordinates `(a, b, c)` of the perpendicular +/// bisector of `[p, q]` in the Euclidean plane (`ax + by + c = 0`). +/// Same as Java `P2.perpendicularBisector(p, q, Pn.EUCLIDEAN)`. inline Eigen::Vector3d perpendicularBisectorEuclidean(const Eigen::Vector3d& p_h, const Eigen::Vector3d& q_h) { // Dehomogenize @@ -46,8 +44,7 @@ inline Eigen::Vector3d perpendicularBisectorEuclidean(const Eigen::Vector3d& p_h return {d(0), d(1), c}; } -// ── Euclidean distance between two P2 homogeneous points ───────────────────── - +/// Euclidean distance between two P² homogeneous points (dehomogenises both). inline double euclideanDistanceP2(const Eigen::Vector3d& p_h, const Eigen::Vector3d& q_h) { Eigen::Vector2d p = p_h.head<2>() / p_h(2); @@ -57,11 +54,9 @@ inline double euclideanDistanceP2(const Eigen::Vector3d& p_h, // ── Direct Euclidean isometry from two point-frames ────────────────────────── -// Build the 3×3 projective matrix that represents the coordinate frame -// anchored at p0 with p1 defining the positive x-direction. -// Euclidean case: columns are [dehom(p0), unit_dir(p0→p1), perp_dir]. -// -// Template parameter S allows float / double / long double. +/// Build the 3×3 projective frame matrix anchored at `p0` with `p1` +/// defining the positive x-direction (Euclidean case). Columns: +/// `[dehom(p0), unit_dir(p0→p1), perp_dir]`. template Eigen::Matrix makeFrameMatrix(Eigen::Matrix p0_h, Eigen::Matrix p1_h) { @@ -84,11 +79,9 @@ Eigen::Matrix makeFrameMatrix(Eigen::Matrix p0_h, return M; } -// Find the 3×3 Euclidean isometry (as a projective matrix) that maps -// the frame (s1, s2) to the frame (t1, t2). -// -// Corresponds to Java P2.makeDirectIsometryFromFrames(s1, s2, t1, t2, Pn.EUCLIDEAN) -// and P2Big.makeDirectIsometryFromFrames(...) (the BigDecimal / high-precision variant). +/// 3×3 Euclidean isometry (as a projective matrix) that maps the +/// frame `(s1, s2)` to the frame `(t1, t2)`. Same as Java +/// `P2.makeDirectIsometryFromFrames(..., Pn.EUCLIDEAN)`. template Eigen::Matrix makeDirectIsometryFromFramesEuclidean( Eigen::Matrix s1, Eigen::Matrix s2, diff --git a/code/include/period_matrix.hpp b/code/include/period_matrix.hpp index 86647e8..2489fc3 100644 --- a/code/include/period_matrix.hpp +++ b/code/include/period_matrix.hpp @@ -50,6 +50,8 @@ namespace conformallab { // PeriodData // ───────────────────────────────────────────────────────────────────────────── +/// Period-matrix data for a genus-g closed surface. For genus 1 the +/// conformal type is fully captured by `τ = ω₂ / ω₁ ∈ ℍ`. struct PeriodData { /// Lattice generators as complex numbers (one per cut edge). /// omega[i] = translations[i].x() + i·translations[i].y() @@ -63,6 +65,7 @@ struct PeriodData { /// True if τ has been reduced to the standard fundamental domain. bool in_fundamental_domain = false; + /// Genus of the surface = `|omega| / 2`. int genus() const { return static_cast(omega.size()) / 2; } }; @@ -74,6 +77,9 @@ struct PeriodData { // // Returns the reduced τ. Throws if Im(τ) ≤ 0 (not in upper half-plane). // ───────────────────────────────────────────────────────────────────────────── +/// Reduce `τ ∈ ℍ` to the standard SL(2,ℤ) fundamental domain +/// `F = { τ ∈ ℍ : |τ| ≥ 1, −½ ≤ Re τ < ½ }` via the generators +/// `S: τ↦−1/τ` and `T: τ↦τ+1`. Throws if `Im τ ≤ 0`. inline std::complex reduce_to_fundamental_domain(std::complex tau) { if (tau.imag() <= 0.0) { @@ -103,6 +109,8 @@ inline std::complex reduce_to_fundamental_domain(std::complex ta // ───────────────────────────────────────────────────────────────────────────── // is_in_fundamental_domain — check membership in F with tolerance tol. // ───────────────────────────────────────────────────────────────────────────── +/// `true` iff `τ` lies inside the standard SL(2,ℤ) fundamental domain +/// with tolerance `tol`. inline bool is_in_fundamental_domain(std::complex tau, double tol = 1e-9) { if (tau.imag() <= 0.0) return false; @@ -117,6 +125,9 @@ inline bool is_in_fundamental_domain(std::complex tau, double tol = 1e-9 // Computes the period data from the Euclidean holonomy translations. // For genus-1 surfaces, also reduces τ to the fundamental domain. // ───────────────────────────────────────────────────────────────────────────── +/// Compute the period data from the Euclidean holonomy translations. +/// For genus 1, also reduces `τ` to the SL(2,ℤ) fundamental domain +/// when `reduce` is `true` (default). inline PeriodData compute_period_matrix(const HolonomyData& hol, bool reduce = true) { PeriodData pd; diff --git a/code/include/projective_math.hpp b/code/include/projective_math.hpp index 2c8f297..abb5f67 100644 --- a/code/include/projective_math.hpp +++ b/code/include/projective_math.hpp @@ -12,17 +12,15 @@ namespace conformallab { -// Divide a homogeneous vector by its last component. -// Corresponds to Java Pn.dehomogenize(). +/// Dehomogenise: divide a homogeneous vector by its last component. +/// Same as Java `Pn.dehomogenize()`. inline Eigen::VectorXd dehomogenize(const Eigen::VectorXd& p) { return p / p(p.size() - 1); } -// Hyperbolic distance between two homogeneous vectors of the same dimension. -// The last component is the "timelike" coordinate (jReality convention). -// Inner product: = -sum_i p_i*q_i + p_last * q_last -// Distance: arcosh() where p̂ normalises to the hyperboloid. -// Corresponds to Java Pn.distanceBetween(p, q, Pn.HYPERBOLIC). +/// Hyperbolic distance `arcosh(⟨p̂, q̂⟩)` between two homogeneous +/// vectors (last component = timelike coordinate; jReality convention). +/// Same as Java `Pn.distanceBetween(p, q, Pn.HYPERBOLIC)`. inline double hyperbolicDistance(const Eigen::VectorXd& p, const Eigen::VectorXd& q) { int n = static_cast(p.size()); @@ -34,10 +32,8 @@ inline double hyperbolicDistance(const Eigen::VectorXd& p, return std::acosh(std::max(1.0, inner)); } -// Check whether a homogeneous point p lies on the segment [s[0], s[1]]. -// Works for n-dimensional homogeneous coords; cross product uses the first -// 3 spatial components after dehomogenization (matching jReality's Rn behaviour). -// Corresponds to Java SurfaceCurveUtility.isOnSegment(). +/// `true` iff the homogeneous point `p_h` lies on the segment +/// `[s0_h, s1_h]`. Same as Java `SurfaceCurveUtility.isOnSegment()`. inline bool isOnSegment(const Eigen::VectorXd& p_h, const Eigen::VectorXd& s0_h, const Eigen::VectorXd& s1_h) { @@ -63,10 +59,10 @@ inline bool isOnSegment(const Eigen::VectorXd& p_h, return true; } -// Find the point on `target` that corresponds to `p` on `source`. -// The parameter t is determined by hyperbolic distance ratios on `source`, -// then applied as a linear interpolation on the dehomogenized `target`. -// Corresponds to Java SurfaceCurveUtility.getPointOnCorrespondingSegment(). +/// Find the point on the target segment `(tgt0, tgt1)` corresponding +/// to `p` on the source segment `(src0, src1)`, parametrised by +/// hyperbolic distance ratios on the source. Same as Java +/// `SurfaceCurveUtility.getPointOnCorrespondingSegment()`. inline Eigen::VectorXd getPointOnCorrespondingSegment( const Eigen::VectorXd& p, const Eigen::VectorXd& src0, diff --git a/code/include/serialization.hpp b/code/include/serialization.hpp index dce3048..7ab6154 100644 --- a/code/include/serialization.hpp +++ b/code/include/serialization.hpp @@ -43,7 +43,7 @@ namespace conformallab { // JSON // ════════════════════════════════════════════════════════════════════════════ -// Save solver result (+ optional 2D layout) to a JSON file. +/// Save the Newton-solver result (+ optional 2-D layout) to a JSON file. inline void save_result_json( const std::string& path, const NewtonResult& res, @@ -80,8 +80,8 @@ inline void save_result_json( ofs << std::setw(2) << j << "\n"; } -// Load DOF vector from a JSON result file. -// Returns the DOF vector; fills out res fields if non-null. +/// Load a DOF vector from a JSON result file written by +/// `save_result_json`. If `res` is non-null its fields are filled too. inline std::vector load_result_json( const std::string& path, NewtonResult* res = nullptr, @@ -181,7 +181,7 @@ inline std::vector parse_doubles(const std::string& s) } // namespace detail_xml -// Save solver result (+ optional layout) to an XML file. +/// Save the Newton-solver result (+ optional layout) to an XML file. inline void save_result_xml( const std::string& path, const NewtonResult& res, @@ -229,8 +229,9 @@ inline void save_result_xml( ofs << "\n"; } -// Load solver result from an XML file written by save_result_xml. -// Returns the DOF vector; fills res/geom/layout2d if non-null. +/// Load a DOF vector from an XML result file written by +/// `save_result_xml`. If `res`, `geom`, `layout2d` are non-null they +/// are filled as well. inline std::vector load_result_xml( const std::string& path, NewtonResult* res = nullptr, diff --git a/code/include/spherical_functional.hpp b/code/include/spherical_functional.hpp index 7ae020b..4abf34e 100644 --- a/code/include/spherical_functional.hpp +++ b/code/include/spherical_functional.hpp @@ -40,19 +40,24 @@ namespace conformallab { // ── Property-map type aliases ───────────────────────────────────────────────── +/// Property map vertex → `double` for the Spherical functional. using SpherVMapD = ConformalMesh::Property_map; +/// Property map vertex → `int` for the Spherical functional. using SpherVMapI = ConformalMesh::Property_map; +/// Property map edge → `double` for the Spherical functional. using SpherEMapD = ConformalMesh::Property_map; +/// Property map edge → `int` for the Spherical functional. using SpherEMapI = ConformalMesh::Property_map; // ── Persistent map bundle ───────────────────────────────────────────────────── +/// Bundle of the five property maps consumed by the Spherical functional. struct SphericalMaps { - SpherVMapI v_idx; // DOF index per vertex (-1 = pinned / u_v = 0) - SpherEMapI e_idx; // DOF index per edge (-1 = no edge DOF) - SpherVMapD theta_v; // target cone angle Θ_v (default 2π) - SpherEMapD theta_e; // target edge angle θ_e (default π) - SpherEMapD lambda0; // base log-length λ°_e (default 0.0) + SpherVMapI v_idx; ///< DOF index per vertex (−1 = pinned / u_v = 0). + SpherEMapI e_idx; ///< DOF index per edge (−1 = no edge DOF). + SpherVMapD theta_v; ///< Target cone angle Θᵥ (default 2π). + SpherEMapD theta_e; ///< Target edge angle θₑ (default π). + SpherEMapD lambda0; ///< Base log-length λ⁰ₑ (default 0). }; // Defaults: theta_v = 2π, theta_e = π, lambda0 = 0. @@ -133,18 +138,21 @@ inline void compute_lambda0_from_mesh(ConformalMesh& mesh, SphericalMaps& m) // ── Evaluation result ───────────────────────────────────────────────────────── +/// Output of `evaluate_spherical()` — energy plus optional gradient. struct SphericalResult { - double energy = 0.0; - std::vector gradient; + double energy = 0.0; ///< Functional value at input DOFs. + std::vector gradient; ///< Gradient ∇E (empty if not requested). }; // ── Internal helpers ────────────────────────────────────────────────────────── +/// Read DOF value from `x` for index `idx`; return 0 if pinned (idx < 0). static inline double spher_dof_val(int idx, const std::vector& x) { return idx >= 0 ? x[static_cast(idx)] : 0.0; } +/// Convert a CGAL half-edge index to a plain `std::size_t` for vector indexing. static inline std::size_t spher_hidx(Halfedge_index h) { return static_cast(static_cast(h)); @@ -152,14 +160,13 @@ static inline std::size_t spher_hidx(Halfedge_index h) // ── Gradient only (no energy) ───────────────────────────────────────────────── -// Compute gradient G(x). -// G_v = Θ_v − Σ_faces α_v(face) -// G_e = α_opp(face+) + α_opp(face−) − θ_e -// -// The corner angle α_v is stored on halfedges using the convention: -// h_alpha[h] = corner angle at source(prev(h)) = corner angle at the vertex -// ACROSS FROM the edge of halfedge h in its face. -// This convention makes both the vertex and edge gradient accumulators natural. +/// Compute the Spherical-functional gradient G(x): +/// * `G_v = Θ_v − Σ_faces α_v(face)` +/// * `G_e = α_opp(face⁺) + α_opp(face⁻) − θ_e` +/// +/// The corner angle α_v is stored on half-edges via the convention +/// `h_alpha[h] = corner angle at the vertex ACROSS FROM the edge of h +/// in its face`, which makes both gradient accumulators natural. inline std::vector spherical_gradient( ConformalMesh& mesh, const std::vector& x, @@ -274,6 +281,9 @@ inline std::vector spherical_gradient( // // 10-point GL nodes and weights on [0, 1] (transformed from [-1, 1]): // t_k = (1 + s_k) / 2, w_k = w_GL_k / 2 +/// Spherical energy `E(x) = ∫₀¹ ⟨G(t·x), x⟩ dt`, evaluated with +/// 10-point Gauss-Legendre quadrature. This is the correct potential +/// for any conservative `G = ∇E`; error ≈ O(h²⁰) for smooth G. inline double spherical_energy( ConformalMesh& mesh, const std::vector& x, @@ -318,6 +328,8 @@ inline double spherical_energy( // ── Full evaluation (energy + gradient) ────────────────────────────────────── +/// Evaluate the Spherical functional at DOFs `x`. Returns energy and +/// gradient (toggle via `need_energy` / `need_gradient`). inline SphericalResult evaluate_spherical( ConformalMesh& mesh, const std::vector& x, @@ -333,10 +345,8 @@ inline SphericalResult evaluate_spherical( return res; } -// ── Finite-difference gradient check ───────────────────────────────────────── -// -// Tests |G[i] − fd[i]| / max(1, |G[i]|) < tol for all DOFs. -// Same defaults as the hyper-ideal gradient check (Java FunctionalTest). +/// Finite-difference gradient check for the Spherical functional +/// (central differences). Same defaults as the Java `FunctionalTest`. inline bool gradient_check_spherical( ConformalMesh& mesh, const std::vector& x0, @@ -390,6 +400,8 @@ inline bool gradient_check_spherical( // // Returns 0.0 if the zero cannot be bracketed (already at gauge maximum, // or open surface — no shift needed). +/// Find the global-scale gauge shift `t*` for the closed-spherical case +/// (see comment block above for the maths). Apply via `apply_spherical_gauge`. inline double spherical_gauge_shift( ConformalMesh& mesh, const std::vector& x, @@ -470,7 +482,8 @@ inline double spherical_gauge_shift( return t; } -// Apply the gauge shift in-place: x_v ← x_v + t* for all variable vertices. +/// Apply the spherical gauge shift in-place: `x_v ← x_v + t*` for every +/// variable vertex, where `t* = spherical_gauge_shift(mesh, x, m, ...)`. inline void apply_spherical_gauge( ConformalMesh& mesh, std::vector& x, diff --git a/code/include/spherical_geometry.hpp b/code/include/spherical_geometry.hpp index 03933cc..2185e71 100644 --- a/code/include/spherical_geometry.hpp +++ b/code/include/spherical_geometry.hpp @@ -23,8 +23,8 @@ constexpr double PI_SPHER = PI; // ── Effective spherical arc length ──────────────────────────────────────────── -// l(λ) = 2·asin(min(exp(λ/2), 1)). -// Clamps exp(λ/2) to [0, 1] so the arcsin stays in domain. +/// Spherical arc length `l(λ) = 2·asin(min(exp(λ/2), 1))`. +/// Clamps `exp(λ/2)` to `[0, 1]` so `asin` stays in domain. inline double spherical_l(double lambda) { double half = std::exp(lambda * 0.5); @@ -35,29 +35,18 @@ inline double spherical_l(double lambda) // ── Interior angles of a spherical triangle ────────────────────────────────── +/// Interior angles of a spherical triangle, plus a `valid` flag. struct SphericalFaceAngles { - double alpha1, alpha2, alpha3; // corner angles at v1, v2, v3 - bool valid; // false when the three lengths fail the - // spherical triangle inequality + double alpha1; ///< Corner angle at vertex v₁. + double alpha2; ///< Corner angle at vertex v₂. + double alpha3; ///< Corner angle at vertex v₃. + bool valid; ///< `false` when the three lengths violate the spherical triangle inequality. }; -// Compute corner angles from spherical arc lengths using the half-angle formula. -// -// Convention (matching the halfedge cycle h0→v1→v2, h1→v2→v3, h2→v3→v1): -// l12 – arc length of edge opposite v3 (edge e12) -// l23 – arc length of edge opposite v1 (edge e23) -// l31 – arc length of edge opposite v2 (edge e31) -// -// Half-angle formula (spherical law of cosines): -// α_k = 2·atan2(sqrt(sin(s-a)·sin(s-b)), sqrt(sin(s)·sin(s-c))) -// where a,b are the two edges ADJACENT to vertex k, c is the opposite edge. -// -// Equivalently (in terms of s-deficiencies): -// α1 = 2·atan2( sqrt(sin(s12)·sin(s31)), sqrt(sin(s)·sin(s23)) ) -// α2 = 2·atan2( sqrt(sin(s12)·sin(s23)), sqrt(sin(s)·sin(s31)) ) -// α3 = 2·atan2( sqrt(sin(s23)·sin(s31)), sqrt(sin(s)·sin(s12)) ) -// -// where s = (l12+l23+l31)/2 and s_ij = s - l_ij. +/// Compute the spherical-triangle corner angles `(α₁, α₂, α₃)` from +/// the three arc lengths `(l₁₂, l₂₃, l₃₁)` using the half-angle form +/// of the spherical law of cosines. Returns `valid = false` for +/// degenerate or out-of-range triangles. inline SphericalFaceAngles spherical_angles(double l12, double l23, double l31) { double s = (l12 + l23 + l31) * 0.5; diff --git a/code/include/spherical_hessian.hpp b/code/include/spherical_hessian.hpp index 59570d8..57027ab 100644 --- a/code/include/spherical_hessian.hpp +++ b/code/include/spherical_hessian.hpp @@ -49,8 +49,18 @@ namespace conformallab { // h2: edge v3-v1 → opposite v2 → w = cot(β2), β2=(π-α3-α1+α2)/2 // // Returns valid=false if any β_k is out of range (degenerate face). -struct SpherCotWeights { double w12, w23, w31; bool valid; }; +/// Three spherical "cotangent" weights for the three edges of a face, +/// derived from the per-vertex interior angles `α₁, α₂, α₃` via +/// `w_ij = cot(β_k)` with `β_k = (π − α_i − α_j + α_k) / 2`. +struct SpherCotWeights { + double w12; ///< Weight for edge v₁-v₂ (opposite vertex v₃). + double w23; ///< Weight for edge v₂-v₃ (opposite vertex v₁). + double w31; ///< Weight for edge v₃-v₁ (opposite vertex v₂). + bool valid; ///< `false` when any β_k is out of `(0, π/2]` (degenerate face). +}; +/// Compute the three spherical cot weights from the three interior +/// angles `(α₁, α₂, α₃)` of a spherical triangle. See `SpherCotWeights`. inline SpherCotWeights spherical_cot_weights(double alpha1, double alpha2, double alpha3) { // β for each edge: @@ -79,25 +89,10 @@ inline SpherCotWeights spherical_cot_weights(double alpha1, double alpha2, doubl return {1.0 / tb3, 1.0 / tb1, 1.0 / tb2, true}; } -// ── Analytical Hessian ──────────────────────────────────────────────────────── -// -// Returns the n×n sparse Hessian matrix H where n = spherical_dimension(mesh, m). -// x – current DOF vector. -// -// Derivation: G_v = θ_v − Σ_f α_v^f → H[i,j] = −Σ_f ∂α_i^f/∂u_j -// -// For a face (v1,v2,v3) with arc-lengths l12,l23,l31 and angles α1,α2,α3, -// differentiating the spherical law of cosines -// cos(l_opp) = cos(l_a)cos(l_b) + sin(l_a)sin(l_b)cos(α) -// gives: -// ∂α1/∂l12 = [cot(l12)cos(α1) − cot(l31)] / sin(α1) (adjacent side) -// ∂α1/∂l31 = [cot(l31)cos(α1) − cot(l12)] / sin(α1) (adjacent side) -// ∂α1/∂l23 = sin(l23) / [sin(l12)sin(l31)sin(α1)] (opposite side) -// -// Chain rule with ∂l_ij/∂u_k = tan(l_ij/2) (from l = 2·asin(exp(λ/2))): -// ∂α1/∂u1 = ∂α1/∂l12·t12 + ∂α1/∂l31·t31 -// ∂α1/∂u2 = ∂α1/∂l12·t12 + ∂α1/∂l23·t23 -// ∂α1/∂u3 = ∂α1/∂l23·t23 + ∂α1/∂l31·t31 +/// Analytical Spherical Hessian via `∂α/∂u` from the spherical law of +/// cosines + chain rule `∂l/∂u = tan(l/2)`; returns an n×n sparse +/// matrix with `n = spherical_dimension(mesh, m)`. See block comment +/// inside the body for the per-face derivation. inline Eigen::SparseMatrix spherical_hessian( ConformalMesh& mesh, const std::vector& x, @@ -214,10 +209,8 @@ inline Eigen::SparseMatrix spherical_hessian( return H; } -// ── Finite-difference Hessian check ────────────────────────────────────────── -// -// Compares the analytical Hessian column-by-column against -// H_fd[:, j] = (G(x + ε·eⱼ) − G(x − ε·eⱼ)) / (2ε). +/// FD Hessian check for the Spherical functional. Compares analytic +/// `H` column-by-column to `(G(x+εeⱼ) − G(x−εeⱼ)) / (2ε)`. inline bool hessian_check_spherical( ConformalMesh& mesh, const std::vector& x0, diff --git a/code/include/viewer_utils.h b/code/include/viewer_utils.h index 4140d11..1f00e11 100644 --- a/code/include/viewer_utils.h +++ b/code/include/viewer_utils.h @@ -5,7 +5,9 @@ namespace viewer_utils { -// Deklaration (Implementation in viewer.cpp) +/// Open an interactive libigl OpenGL viewer window showing the mesh +/// `(V, F)`. Built only when `WITH_VIEWER=ON`; declaration here, body +/// in `viewer.cpp`. void simple_visualize(Eigen::MatrixXd& V, Eigen::MatrixXi& F); } \ No newline at end of file