diff --git a/README.md b/README.md index 3b13a3f..82e9637 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 vollständig abgeschlossen. Alle drei Kern-Geometrien (Hyper-Ideal, Sphärisch, Euklidisch) plus analytische Hessians (Kotangenten-Laplace) sind auf `ConformalMesh` portiert. **62 Tests, 3 skipped**. Nächster Schritt: Phase 4 (Newton-Solver) — siehe [Roadmap](#roadmap). +> **Status:** Phase 4 vollständig abgeschlossen. Alle drei Geometrien lösbar via Newton-Solver (SimplicialLDLT + SparseQR-Fallback). Interaktiver Viewer-Beispiel inklusive. **87 Tests, 2 skipped**. Nächster Schritt: Phase 5 (CLI-App + Serialisierung). --- @@ -14,38 +14,130 @@ The long-term goal is a **CGAL package** that brings discrete conformal maps (hy |------|--------| | Clausen / Lobachevsky / ImLi₂ functions | ✅ Phase 1 | | Hyper-ideal geometry (ζ, lᵢⱼ, αᵢⱼ, σᵢ, σᵢⱼ) | ✅ Phase 2 | -| 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 (energy + gradient, CGAL mesh) | ✅ Phase 3d | -| Spherical gauge-fix (scale gauge for closed surfaces) | ✅ Phase 3e | -| Euclidean Hessian (cotangent-Laplace operator) | ✅ Phase 3f | +| Hyper-ideal functional (energy + gradient) | ✅ Phase 3b | +| Spherical functional (energy + gradient + gauge-fix) | ✅ Phase 3c/3e | +| Euclidean functional (energy + gradient) | ✅ Phase 3d | +| Euclidean Hessian (cotangent-Laplace, Pinkall–Polthier) | ✅ Phase 3f | | Spherical Hessian (∂α/∂u from law of cosines) | ✅ Phase 3f | -| Solvers (Newton, gradient flow) | 🔜 Phase 4 | -| XML serialisation / mesh I/O | 🔜 Phase 5 | -| Full CLI app | 🔜 Phase 5 | +| Hyper-ideal Hessian (numerical FD, symmetrised) | ✅ Phase 4a | +| Newton solver — all three geometries | ✅ Phase 4a | +| **SparseQR fallback** for rank-deficient H (gauge modes) | ✅ Phase 4a | +| Mesh I/O (CGAL::IO — OFF / OBJ / PLY) | ✅ Phase 4b | +| End-to-end pipeline tests | ✅ Phase 4c | +| **Example programs** (headless + interactive viewer) | ✅ Phase 4d | +| CLI app + XML serialisation | 🔜 Phase 5 | + +--- + +## Quick start — running the examples + +```bash +cmake -S code -B build -DWITH_CGAL=ON +cmake --build build --target example_euclidean example_hyper_ideal + +# Euclidean conformal map (headless) +./build/examples/example_euclidean [input.off] [output.off] + +# HyperIdeal conformal map (headless) +./build/examples/example_hyper_ideal [input.off] [output.off] + +# Interactive viewer (requires -DWITH_VIEWER=ON) +cmake -S code -B build -DWITH_CGAL=ON -DWITH_VIEWER=ON +cmake --build build --target example_viewer +./build/examples/example_viewer [input.off] +``` + +If no input file is given each example uses a built-in test mesh. + +--- + +## Library usage — minimal Euclidean pipeline + +```cpp +#include "conformal_mesh.hpp" +#include "mesh_builder.hpp" +#include "mesh_io.hpp" +#include "euclidean_functional.hpp" +#include "newton_solver.hpp" + +using namespace conformallab; + +int main() { + // 1. Load mesh + ConformalMesh mesh = load_mesh("input.off"); + + // 2. Set up functional maps + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + + // 3. Assign DOFs — pin first vertex (gauge fix) + auto vit = mesh.vertices().begin(); + maps.v_idx[*vit++] = -1; // pinned: u[v0] = 0 + int idx = 0; + for (; vit != mesh.vertices().end(); ++vit) + maps.v_idx[*vit] = idx++; + const int n = idx; + + // 4. Set target angles (natural equilibrium: x* = 0) + std::vector x0(n, 0.0); + auto G0 = euclidean_gradient(mesh, x0, maps); + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv >= 0) maps.theta_v[v] -= G0[iv]; + } + + // 5. Solve + auto result = newton_euclidean(mesh, std::vector(n, -0.1), maps); + + // 6. Save result + if (result.converged) + save_mesh("output.off", mesh); + + return result.converged ? 0 : 1; +} +``` + +For **HyperIdeal** geometry: + +```cpp +auto maps = setup_hyper_ideal_maps(mesh); +int n = assign_all_dof_indices(mesh, maps); +// … set theta_v / theta_e targets … +auto result = newton_hyper_ideal(mesh, x0, maps); +``` + +Using **`solve_linear_system`** directly (with fallback detection): + +```cpp +#include "newton_solver.hpp" + +bool used_fallback = false; +auto dx = conformallab::solve_linear_system(H, rhs, &used_fallback); +if (used_fallback) + std::cout << "SparseQR was used (H is rank-deficient)\n"; +``` --- ## Build modes -| Mode | CMake flag | What gets built | CI | -|------|-----------|-----------------|-----| -| **Tests only** (default) | *(none)* | `conformallab_tests` · Eigen + GTest | ✅ runs automatically | -| **CGAL tests** | `-DWITH_CGAL=ON` | `conformallab_cgal_tests` · above + CGAL + system Boost | local only | -| **Viewer** | `-DWITH_VIEWER=ON` | `viewer` library · libigl / GLFW / GLAD | local only | -| **Full app** | `-DWITH_CGAL=ON` | `conformallab_core` CLI · all of the above | local only | +| Mode | CMake flags | What gets built | CI | +|------|------------|-----------------|-----| +| **Tests only** (default) | *(none)* | `conformallab_tests` — Eigen + GTest only | ✅ automatic | +| **CGAL tests + examples** | `-DWITH_CGAL=ON` | `conformallab_cgal_tests`, `example_euclidean`, `example_hyper_ideal` | local only | +| **Interactive viewer** | `-DWITH_CGAL=ON -DWITH_VIEWER=ON` | above + `example_viewer`, `conformallab_core` | local only | -External dependencies are bundled as tarballs in `code/deps/tarballs/` and extracted lazily at CMake configure time (GTest is fetched from GitHub via `FetchContent`). +External dependencies are bundled as tarballs in `code/deps/tarballs/` and extracted lazily at CMake configure time (GTest is fetched via `FetchContent`). -**Note on Boost:** `-DWITH_CGAL=ON` requires a system-installed Boost (header-only use by CGAL). The default `Tests only` mode needs no Boost. +**Boost** is required only with `-DWITH_CGAL=ON` (header-only use by CGAL 6.x). --- ## Prerequisites -| Tool | Minimum version | -|------|----------------| +| Tool | Minimum | +|------|---------| | C++ compiler (GCC or Clang) | C++17 | | CMake | 3.20 | | Boost headers | 1.70 *(only with `-DWITH_CGAL=ON`)* | @@ -59,7 +151,7 @@ git clone https://codeberg.org/TMoussa/ConformalLabpp cd ConformalLabpp ``` -### Tests only (CI default — no system deps needed) +### Tests only (CI default — no system deps) ```bash cmake -S code -B build @@ -67,22 +159,24 @@ cmake --build build --target conformallab_tests -j$(nproc) ctest --test-dir build --output-on-failure ``` -### CGAL mesh + functional tests (requires system Boost) +### CGAL tests + headless examples (needs system Boost) ```bash cmake -S code -B build -DWITH_CGAL=ON -cmake --build build --target conformallab_cgal_tests -j$(nproc) +cmake --build build --target conformallab_cgal_tests example_euclidean example_hyper_ideal -j$(nproc) ctest --test-dir build -R "^cgal\." --output-on-failure +./build/examples/example_euclidean +./build/examples/example_hyper_ideal ``` -Expected output: **45 tests pass, 3 skipped** (the three `@Ignore` Hessian stubs, one per functional). +Expected: **87 tests pass, 2 skipped** (the two `@Ignore` Hessian stubs). -### Full CLI app (CGAL + viewer) +### Interactive viewer ```bash -cmake -S code -B build -DWITH_CGAL=ON -cmake --build build -j$(nproc) -./code/bin/conformallab_core --input data/off/example.off --show +cmake -S code -B build -DWITH_CGAL=ON -DWITH_VIEWER=ON +cmake --build build --target example_viewer -j$(nproc) +./build/examples/example_viewer data/off/example.off ``` --- @@ -91,19 +185,24 @@ cmake --build build -j$(nproc) | Header | Description | |--------|-------------| -| `clausen.hpp` | Clausen integral Cl₂, Lobachevsky Л, ImLi₂ | -| `hyper_ideal_geometry.hpp` | Pure-math ζ functions, lᵢⱼ, αᵢⱼ, σᵢ, σᵢⱼ | +| `clausen.hpp` | Clausen Cl₂, Lobachevsky Л, ImLi₂ | +| `hyper_ideal_geometry.hpp` | ζ functions, lᵢⱼ, αᵢⱼ, σᵢ, σᵢⱼ | | `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 + 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 | -| `matrix_utility.hpp` | Small linear-algebra helpers | -| `projective_math.hpp` | Projective geometry utilities | +| `hyper_ideal_functional.hpp` | HyperIdeal energy + gradient on `ConformalMesh` | +| `hyper_ideal_hessian.hpp` | HyperIdeal Hessian (numerical FD, symmetrised) | +| `hyper_ideal_visualization_utility.hpp` | Poincaré disk projection, circumcircle helpers | +| `spherical_geometry.hpp` | Spherical arc length, half-angle formula | +| `spherical_functional.hpp` | Spherical energy + gradient + gauge-fix | +| `spherical_hessian.hpp` | Spherical Hessian (∂α/∂u, law of cosines) | +| `euclidean_geometry.hpp` | Euclidean corner-angle (t-value / atan2) | +| `euclidean_functional.hpp` | Euclidean energy + gradient | +| `euclidean_hessian.hpp` | Cotangent-Laplace Hessian (Pinkall–Polthier) | +| `newton_solver.hpp` | `newton_euclidean` / `newton_spherical` / `newton_hyper_ideal` + public **`solve_linear_system`** | +| `mesh_io.hpp` | `read_mesh` / `write_mesh` / `load_mesh` / `save_mesh` | +| `conformal_mesh.hpp` | `ConformalMesh` = `CGAL::Surface_mesh` + property-map helpers | +| `mesh_builder.hpp` | `make_triangle` / `make_tetrahedron` / `make_quad_strip` / `make_fan` / `make_spherical_tetrahedron` | +| `mesh_utils.hpp` | CGAL → Eigen conversion (`cgal_to_eigen`) | +| `constants.hpp` | `conformallab::PI`, `TWO_PI` | --- @@ -111,86 +210,383 @@ cmake --build build -j$(nproc) ``` code/ -├── include/ -│ ├── conformal_mesh.hpp # CGAL mesh type + property-map helpers -│ ├── mesh_builder.hpp # make_triangle / make_tetrahedron / … -│ ├── 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 + 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 +├── include/ # All public headers (header-only library) +│ ├── conformal_mesh.hpp +│ ├── mesh_builder.hpp +│ ├── mesh_io.hpp +│ ├── mesh_utils.hpp +│ ├── newton_solver.hpp # ← public solve_linear_system + 3 Newton solvers +│ ├── hyper_ideal_{functional,hessian,geometry,utility,visualization_utility}.hpp +│ ├── spherical_{functional,hessian,geometry}.hpp +│ ├── euclidean_{functional,hessian,geometry}.hpp +│ ├── clausen.hpp +│ └── constants.hpp +├── examples/ # Standalone example programs +│ ├── CMakeLists.txt +│ ├── example_euclidean.cpp # Headless Euclidean pipeline +│ ├── example_hyper_ideal.cpp # Headless HyperIdeal pipeline +│ └── example_viewer.cpp # Interactive libigl viewer (WITH_VIEWER) ├── src/ -│ ├── apps/v0/ # conformallab_core CLI (requires WITH_CGAL) -│ └── viewer/ # simple_viewer (requires WITH_VIEWER) +│ ├── apps/v0/conformallab_cli.cpp # CLI skeleton (Phase 5) +│ └── viewer/simple_viewer.cpp ├── tests/ │ ├── CMakeLists.txt -│ ├── *.cpp # conformallab_tests (no CGAL) +│ ├── *.cpp # conformallab_tests (no CGAL) │ └── cgal/ │ ├── CMakeLists.txt -│ ├── test_conformal_mesh.cpp # 14 mesh infrastructure tests -│ ├── test_hyper_ideal_functional.cpp # 6 hyper-ideal gradient checks -│ ├── test_spherical_functional.cpp # 11 spherical gradient + gauge-fix checks -│ └── test_euclidean_functional.cpp # 11 Euclidean gradient checks +│ ├── test_conformal_mesh.cpp # 14 tests +│ ├── test_hyper_ideal_functional.cpp # 7 tests (1 skipped) +│ ├── test_spherical_functional.cpp # 11 tests (1 skipped) +│ ├── test_euclidean_functional.cpp # 11 tests +│ ├── test_euclidean_hessian.cpp # 8 tests +│ ├── test_spherical_hessian.cpp # 8 tests +│ ├── test_newton_solver.cpp # 14 tests (incl. 3 SparseQR tests) +│ ├── test_mesh_io.cpp # 6 tests +│ └── test_pipeline.cpp # 5 tests └── deps/ - ├── tarballs/ # bundled dependency archives - ├── eigen-3.4.0/ # header-only linear algebra (always extracted) - ├── CGAL-6.1.1/ # header-only geometry (extracted with WITH_CGAL) - ├── libigl-2.6.0/ # header-only viewer toolkit (extracted with WITH_VIEWER) - ├── glfw-3.4/ # windowing (extracted with WITH_VIEWER) - ├── libigl-glad/ # OpenGL loader (extracted with WITH_VIEWER) - └── single_includes/ # CLI11, json.hpp + ├── eigen-3.4.0/ # always extracted + ├── CGAL-6.1.1/ # extracted with WITH_CGAL + ├── libigl-2.6.0/ # extracted with WITH_VIEWER + ├── glfw-3.4/ # extracted with WITH_VIEWER + ├── libigl-glad/ + └── single_includes/ # CLI11, json.hpp ``` --- ## Test suites -### `conformallab_tests` (always built, runs in CI) +### `conformallab_tests` (CI — always built) -Pure-math tests requiring only Eigen: +Pure-math tests requiring only Eigen: Clausen / Lobachevsky / ImLi₂, hyper-ideal geometry, tetrahedron volumes. -- Clausen function, Lobachevsky function, ImLi₂ -- Hyper-ideal geometry (ζ₁₃, ζ₁₄, ζ₁₅, lᵢⱼ, αᵢⱼ) -- Tetrahedron volume formulas (Meyerhoff, Kolpakov–Mednykh) +### `conformallab_cgal_tests` (local — `-DWITH_CGAL=ON`) -### `conformallab_cgal_tests` (built with `-DWITH_CGAL=ON`) - -CGAL `Surface_mesh` tests (test prefix `cgal.`): - -| Suite | Tests | Description | -|-------|-------|-------------| +| Suite | Tests | What it checks | +|-------|------:|----------------| | `ConformalMeshTopology` | 4 | Euler characteristic, vertex/edge/face counts | -| `ConformalMeshTraversal` | 4 | Halfedge iteration, valence, opposite-halfedge | -| `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` | 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) | +| `ConformalMeshTraversal` | 4 | Halfedge iteration, valence, opposite | +| `ConformalMeshProperties` | 5 | Property maps (λ, θ, idx, α, geometry type) | +| `ConformalMeshValidity` | 1 | CGAL validity for all factory meshes | +| `HyperIdealFunctional` | 7 | FD gradient checks + Hessian symmetry | +| `SphericalFunctional` | 11 | Angle formula + gradient + gauge-fix | +| `EuclideanFunctional` | 11 | Angle formula + gradient | +| `EuclideanHessian` | 8 | Cotangent-Laplace structure, FD agreement, PSD, null space | +| `SphericalHessian` | 8 | Derivative correctness, NSD at equilibrium | +| `NewtonSolver` | 11 | Convergence (Euclidean ×3, Spherical ×4, HyperIdeal ×4) | +| `SparseQRFallback` | 3 | Full-rank LDLT path · singular matrix triggers QR · closed-mesh gauge-mode | +| `MeshIO` | 6 | OFF/OBJ round-trips, error handling | +| `Pipeline` | 5 | End-to-end: build → setup → solve → export → reload, all three geometries | +| **Total** | **87** | 2 skipped (Hessian stubs) | + +--- + +## Newton solver & SparseQR fallback + +`newton_solver.hpp` exposes three solvers with a unified interface: + +``` +NewtonResult newton_euclidean (mesh, x0, maps [, tol, max_iter]) +NewtonResult newton_spherical (mesh, x0, maps [, tol, max_iter]) +NewtonResult newton_hyper_ideal(mesh, x0, maps [, tol, max_iter, hess_eps]) +``` + +Each iteration: +1. Evaluate gradient **G** +2. Evaluate Hessian **H** (analytical for Euclidean / Spherical; numerical FD for HyperIdeal) +3. Solve **H·Δx = −G** — try `Eigen::SimplicialLDLT`, fall back to `Eigen::SparseQR` on failure +4. Backtracking line search (up to 20 halvings) + +**Gradient sign conventions:** + +| Geometry | **G** | **H** sign | +|----------|-------|-----------| +| Euclidean | Θ_v − Σα_v | PSD → LDLT on H | +| Spherical | Θ_v − Σα_v | NSD → LDLT on **−H** | +| HyperIdeal | Σβ_v − Θ_v | PSD → LDLT on H | + +**SparseQR fallback** (`solve_linear_system`): +When `SimplicialLDLT` fails (singular/rank-deficient **H**, e.g. gauge modes on closed meshes without a pinned vertex), `SparseQR` finds the minimum-norm Newton step orthogonal to the null space. Because the gradient always lies in the row space of **H**, the solver converges correctly without requiring the caller to pin a vertex. + +The fallback is a **public API**: + +```cpp +bool used_fallback = false; +auto dx = conformallab::solve_linear_system(H, rhs, &used_fallback); +``` + +--- + +## Mathematical scope — C++ vs. Java original + +The three core functionals are fully equivalent to the Java original at the level of energy, gradient, and Hessian formulas. The table below shows where parity holds, where there is a numerical difference, and what is not yet ported. + +| Mathematical layer | Java ConformalLab | conformallab++ | +|---|---|---| +| **Euclidean functional** — energy, gradient | ✅ | ✅ | +| **Spherical functional** — energy, gradient, gauge-fix | ✅ | ✅ | +| **HyperIdeal functional** — energy, gradient | ✅ | ✅ | +| **Inversive-distance functional** (Luo 2004, Bowers–Stephenson) | ✅ | ❌ not ported | +| **Euclidean Hessian** — cotangent-Laplace (Pinkall–Polthier 1993) | ✅ analytical | ✅ analytical | +| **Spherical Hessian** — ∂α/∂u from law of cosines | ✅ analytical | ✅ analytical | +| **HyperIdeal Hessian** — through ζ → lᵢⱼ → β/α chain | ✅ analytical | ⚠️ symmetric FD | +| **Newton solver** | ✅ | ✅ | +| SparseQR fallback for gauge-mode null spaces | ? | ✅ | +| **Cone metrics** — prescribed Θ_v ≠ 2π | ✅ full pipeline | ⚠️ data structure only | +| **Boundary conditions** — Dirichlet u=f, Neumann, free boundary | ✅ | ⚠️ pin-only | +| **Layout / embedding** — DOF vector → vertex coordinates in ℝ² / H² / S² | ✅ | ❌ not implemented | +| **Gauss–Bonnet consistency check** on target angles | ✅ | ❌ | +| **Global uniformization** for genus g ≥ 1 | ✅ | ❌ | +| **Period matrices** — Teichmüller parameters for g ≥ 2 | ✅ | ❌ | +| **Holonomy / monodromy** | ✅ | ❌ | +| **HyperIdeal generator** — constructing geometrically valid meshes | ✅ | ❌ only test meshes | +| Clausen / Lobachevsky / ImLi₂ special functions | ✅ | ✅ | +| Discrete elliptic utility (modular normalisation of τ) | ✅ | ✅ (not yet wired up) | +| Poincaré disk / Lorentz boost visualisation helpers | ✅ | ✅ | +| Mesh I/O | ✅ XML/CoHDS | ✅ OFF/OBJ/PLY | +| Interactive viewer | ✅ jReality | ✅ libigl/GLFW | + +### Key numerical difference — HyperIdeal Hessian + +The analytical Hessian of the HyperIdeal functional requires differentiating through the chain + +``` +(b_i, a_e) → l_ij → ζ₁₃/ζ₁₄/ζ₁₅ → α_ij / β_i +``` + +which is feasible but involves many nested cases (four vertex-type combinations per edge). Until Phase 5 delivers the analytical version, conformallab++ uses a symmetric finite-difference Hessian: + +``` +H[i,j] = ( G(x + ε·eⱼ)[i] − G(x − ε·eⱼ)[i] ) / (2ε) +``` + +This is O(ε²) accurate (≈ 10⁻¹⁰ relative error at ε = 10⁻⁵), positive semi-definite by strict convexity of the HyperIdeal energy (Springborn 2020), and costs n extra gradient evaluations per Newton step instead of O(n). For meshes with fewer than ~500 DOFs the difference in wall time is negligible. + +### What "cone metrics" and "layout" would require + +**Cone metrics** — the property map `theta_v` is already subtracted in the gradient (`G_v = Σα_v − Θ_v`), so prescribing a cone angle is a one-liner: `maps.theta_v[v] = desired_angle`. What is missing is the *application layer*: checking Gauss–Bonnet consistency (Σ (2π − Θ_v) = 2π·χ), distributing angle defects sensibly, and special handling at boundary vertices. + +**Layout** — after solving you have the conformal scale factors u_i but the vertex positions in the mesh are unchanged. Recovering the actual flat / hyperbolic / spherical coordinates requires integrating the discrete holomorphic differential (discrete Schwarz–Christoffel for the Euclidean case, or geodesic development for the hyperbolic case). This is the single biggest missing step for a complete uniformization pipeline. + +--- + +## For mathematicians — extending the library + +This section explains how to add new functionals, test conjectures numerically, and hook into the existing solver infrastructure, with no assumed prior knowledge of the codebase. + +### Mental model + +The library is built around one central idea: a **discrete conformal functional** E(x) whose critical points are the conformally equivalent metrics. Everything else is infrastructure for evaluating E, its gradient G = ∂E/∂x, and its Hessian H = ∂²E/∂x². + +``` +ConformalMesh — half-edge mesh (CGAL::Surface_mesh) + + property maps — per-vertex / per-edge data (λ, θ, α, DOF index, …) + +Maps struct — collects all property maps for one functional + theta_v[v] — target angle at vertex v (your input) + v_idx[v] — DOF index, or −1 if pinned + e_idx[e] — DOF index for edge DOFs (HyperIdeal only) + +x ∈ ℝⁿ — the DOF vector the solver optimises + +evaluate_*(mesh, x, maps) → { energy, gradient, … } +newton_*(mesh, x0, maps) → { x*, iterations, converged, … } +``` + +The mesh geometry (vertex positions) is only used to initialise the log edge-lengths λ°. From then on the solver works entirely in the `x`-space. + +### Adding a new functional — step-by-step + +Copy `euclidean_functional.hpp` as a template (it is the simplest of the three). You need to provide: + +**1. A `Maps` struct** that holds the property maps your functional needs: + +```cpp +// my_functional.hpp +#pragma once +#include "conformal_mesh.hpp" + +namespace conformallab { + +struct MyMaps { + // property maps attached to the mesh + ConformalMesh::Property_map lambda; // log edge-lengths + ConformalMesh::Property_map theta_v; // target angles + ConformalMesh::Property_map v_idx; // DOF indices + + // any extra parameters your functional needs + double my_parameter = 1.0; +}; + +inline MyMaps setup_my_maps(ConformalMesh& mesh) { … } +``` + +**2. An energy + gradient function:** + +```cpp +struct MyResult { + double energy; + std::vector gradient; +}; + +inline MyResult evaluate_my_functional( + ConformalMesh& mesh, + const std::vector& x, + const MyMaps& m, + bool compute_energy = true) +{ + MyResult res; + res.gradient.assign(x.size(), 0.0); + + for (auto f : mesh.faces()) { + // iterate halfedges around face + // compute your per-face contribution to E and G + // accumulate: res.gradient[m.v_idx[v]] += … + } + + // subtract target-angle term + for (auto v : mesh.vertices()) { + int iv = m.v_idx[v]; + if (iv < 0) continue; + res.gradient[iv] -= m.theta_v[v]; // G_v = actual - target + } + + return res; +} +``` + +**3. A gradient check** — before trusting your formula, verify it numerically. There is a ready-made helper in `hyper_ideal_functional.hpp` you can call directly, or write your own: + +```cpp +// Finite-difference gradient check for any functional +bool my_gradient_check(ConformalMesh& mesh, + const std::vector& x, + const MyMaps& m, + double eps = 1e-6, double tol = 1e-5) +{ + auto r0 = evaluate_my_functional(mesh, x, m, false); + const int n = static_cast(x.size()); + for (int i = 0; i < n; ++i) { + auto xp = x; xp[i] += eps; + auto xm = x; xm[i] -= eps; + double fd = (evaluate_my_functional(mesh, xp, m).energy + - evaluate_my_functional(mesh, xm, m).energy) / (2*eps); + if (std::abs(fd - r0.gradient[i]) > tol * (1 + std::abs(fd))) + return false; + } + return true; +} +``` + +Add a `TEST(MyFunctional, GradientCheck_Triangle)` in `tests/cgal/` and it will be picked up automatically by CTest. + +**4. Hook into the Newton solver.** Once your gradient and Hessian are correct, plug in `solve_linear_system` or write a thin wrapper in the style of `newton_euclidean`: + +```cpp +// Use a numerical Hessian first (safe starting point) +#include "newton_solver.hpp" +#include + +// Build H by FD of your gradient, then: +bool ok = false; +auto dx = detail::solve_with_fallback(H, -G, ok); +``` + +Or supply an analytical Hessian as a sparse matrix and pass it directly. + +### Where the key mathematical objects live + +| Object | File | What to look for | +|--------|------|-----------------| +| Corner angle formula (Euclidean) | `euclidean_geometry.hpp` | `euclidean_corner_angle()` — inputs are log half-edge lengths | +| Spherical angle formula | `spherical_geometry.hpp` | `spherical_corner_angle()` — uses spherical law of cosines | +| HyperIdeal angle (ζ₁₃/ζ₁₄/ζ₁₅) | `hyper_ideal_geometry.hpp` | `zeta13/14/15()`, `alpha_ij()` — the four vertex-type cases | +| Per-face energy term | `*_functional.hpp` | the inner loop over `mesh.faces()` | +| Gradient accumulation | `*_functional.hpp` | `grad[v_idx[v]] += …` after the face loop | +| Cotangent-Laplace structure | `euclidean_hessian.hpp` | `euclidean_hessian()` — shows the sparse-triplet pattern | +| Spherical Hessian derivation | `spherical_hessian.hpp` | comments give the ∂α/∂u formula step by step | +| Special functions | `clausen.hpp` | `Cl2()`, `lobachevsky()`, `imLi2()` — all take a `double` angle | + +### How to navigate the half-edge mesh + +```cpp +for (auto f : mesh.faces()) { + // The three halfedges of face f: + auto h0 = mesh.halfedge(f); + auto h1 = mesh.next(h0); + auto h2 = mesh.next(h1); + + // Vertices opposite to each halfedge (the vertex NOT on h): + Vertex_index v0 = mesh.target(h2); // opposite to edge h0-h1 + Vertex_index v1 = mesh.target(h0); // opposite to edge h1-h2 + Vertex_index v2 = mesh.target(h1); // opposite to edge h0-h2 (= h2 target) + + // Access DOF index (−1 = pinned): + int i0 = maps.v_idx[v0]; + + // The opposite halfedge (for the adjacent face, if not on boundary): + auto h_opp = mesh.opposite(h0); + bool is_boundary = mesh.is_border(h_opp); +} +``` + +### Attaching new data to a mesh + +```cpp +// Add a per-vertex curvature field (survives mesh copy): +auto [curv, created] = mesh.add_property_map("v:my_curv", 0.0); + +// Write and read: +curv[v] = 1.234; +double k = curv[v]; + +// Pass it through your Maps struct so functions can access it. +``` + +Property maps are reference-counted and cheap to copy. Give them unique string names to avoid collision. + +### Quick-start experiment checklist + +1. **Read** `examples/example_euclidean.cpp` — it shows the full pipeline in ~80 lines with comments at every step. +2. **Build** without a viewer first: `cmake -S code -B build -DWITH_CGAL=ON && cmake --build build --target example_euclidean`. +3. **Add a gradient check test** in `tests/cgal/` — copy any `TEST(…, GradientCheck_…)` block and swap out the functional. Run with `ctest -R your_test_name`. +4. **Try different target angles** — set `maps.theta_v[v] = M_PI / 3` for all interior vertices and see how the solver responds. The constraint `Σ(2π − Θ_v) = 2π·χ(M)` (Gauss–Bonnet) must hold for a solution to exist. +5. **Inspect convergence** — `NewtonResult` carries `iterations`, `grad_inf_norm`, and the full `x` at termination. Plot `||G(xₖ)||` per iteration to verify quadratic convergence near the solution. + +### Recommended reading + +| Paper | Relevance to this codebase | +|-------|---------------------------| +| Springborn, Schröder, Pinkall — *Conformal Equivalence of Triangle Meshes* (2008) | Euclidean & spherical functionals; the Schläfli formula at the core of `spherical_functional.hpp` | +| Springborn — *Ideal Hyperbolic Polyhedra and Discrete Uniformization* (2020) | HyperIdeal functional; the ζ₁₃/ζ₁₄/ζ₁₅ functions in `hyper_ideal_geometry.hpp` | +| Pinkall, Polthier — *Computing Discrete Minimal Surfaces* (1993) | Cotangent-Laplace Hessian in `euclidean_hessian.hpp` | +| Luo — *Combinatorial Yamabe Flow on Surfaces* (2004) | Inversive-distance functional (not yet ported — good first contribution) | +| Bobenko, Springborn — *Variational Principles for Circle Patterns* (2004) | Background for the angle-sum variational framework used throughout | --- ## Key design decisions -**CGAL as CoHDS replacement.** `CGAL::Surface_mesh` replaces the Java `CoHDS` half-edge data structure. Vertex/edge/face/halfedge descriptors are `Vertex_index`, `Edge_index`, `Face_index`, `Halfedge_index` (typed integers, not raw handles). +**CGAL as CoHDS replacement.** `CGAL::Surface_mesh` replaces the Java `CoHDS` half-edge data structure. Vertex/edge/face/halfedge descriptors are typed integers — no raw handles, no RTTI. -**Property maps.** `mesh.add_property_map("v:lambda", 0.0)` replaces the Java adapter/decorator pattern. Multiple maps can be attached to one mesh without subclassing. +**Property maps.** `mesh.add_property_map("v:lambda", 0.0)` replaces the Java adapter/decorator pattern. Multiple maps attach to one mesh without subclassing. -**Energy parameterisation.** Both functionals use a **DOF vector** `x` indexed by `v_idx[v]` / `e_idx[e]` (−1 = pinned). This matches the Java `FunctionalTest` gradient-check convention. +**DOF vector convention.** All functionals use `x` indexed by `v_idx[v]` / `e_idx[e]` (−1 = pinned). This matches the Java `FunctionalTest` gradient-check convention and is uniform across all three geometries. -**Spherical edge gradient.** For the spherical parameterisation where `Λᵢⱼ = λ°ᵢⱼ + uᵢ + uⱼ + λₑ` (additive edge DOF), the Schläfli identity gives `∂E/∂λₑ = (2α_opp − S_f)/2` per face (Euclidean limit: `S_f = π`, recovering the familiar `α_opp⁺ + α_opp⁻ − π`). +**HyperIdeal Hessian via FD.** The analytical Hessian through `ζ13/14/15 → lij → β/α` is left for Phase 5. A symmetric FD Hessian `H[i,j] = (G(x+ε·eⱼ)[i] − G(x−ε·eⱼ)[i]) / (2ε)` is O(ε²) accurate, PSD by strict convexity, and sufficient for Newton on < 500 DOFs. + +**Spherical Hessian sign.** The spherical energy is **concave** (not convex) — the Hessian **H** is NSD at equilibrium. Newton solves `(−H)·Δx = G`, so the sign flip is handled transparently inside `newton_spherical`. + +**Natural theta trick.** Tests set `theta_v = Σα_v(x=x_base)` to make `x_base` the known equilibrium, avoiding any need to manufacture reference solutions. For HyperIdeal `x_base = (b=1.0, a=0.5)` is used (x=0 is degenerate in log-space). --- ## CI -Tests run automatically on push to `main`, `dev`, and `claude/**` branches via a self-hosted Gitea Actions runner (`eulernest`, ARM64 Raspberry Pi). The pipeline uses a minimal Docker image (`git.eulernest.eu/conformallab/ci-cpp:latest`) with cmake, g++, git, and Node.js 20 pre-installed. **Only `conformallab_tests` runs in CI** (no Boost/CGAL dependency in the CI image). - -The Dockerfile for the CI image lives in `.gitea/docker/Dockerfile.ci-cpp`. Build and push it once whenever the image needs updating: +Tests run automatically on push to `main`, `dev`, and `claude/**` branches via a self-hosted Gitea Actions runner (`eulernest`, ARM64). The CI image contains cmake, g++, git, and Node.js. **Only `conformallab_tests` runs in CI** (no Boost/CGAL dependency there). ```bash +# Rebuild and push the CI image when the Dockerfile changes docker buildx build \ --platform linux/arm64 \ -f .gitea/docker/Dockerfile.ci-cpp \ @@ -203,56 +599,44 @@ docker buildx build \ ## Roadmap -Phase 3 ist vollständig abgeschlossen. Phase 4 (Newton-Solver) ist der nächste Schritt. -Die Zeitangaben sind grobe Schätzungen. - ---- - ``` +Phase 1 Clausen / Lobachevsky / ImLi₂ ✅ abgeschlossen + +Phase 2 Hyper-ideal Geometrie (ζ, lᵢⱼ, αᵢⱼ, σᵢ) ✅ abgeschlossen + Phase 3a CGAL Surface_mesh Infrastruktur ✅ abgeschlossen - → conformal_mesh.hpp, mesh_builder.hpp - → 14 Tests: Topologie, Traversal, Properties, Validity - -Phase 3b HyperIdealFunctional portieren ✅ abgeschlossen - → hyper_ideal_functional.hpp: Energie + Gradient - → 6 Tests (1 skipped): Gradient-Checks, FD-Tests - -Phase 3c SphericalFunctional portieren ✅ abgeschlossen - → spherical_functional.hpp: Energie + Gradient - → Tests: Winkelformel, Gradient-Checks - -Phase 3d EuclideanCyclicFunctional portieren ✅ abgeschlossen - → euclidean_functional.hpp: Energie + Gradient auf ConformalMesh - → Tests: Winkelformel, Gradient-Checks, NaN-Test, Fan, Mixed-Pinned - +Phase 3b HyperIdealFunctional ✅ abgeschlossen +Phase 3c SphericalFunctional ✅ abgeschlossen +Phase 3d EuclideanCyclicFunctional ✅ abgeschlossen Phase 3e Gauge-Fix für SphericalFunctional ✅ abgeschlossen - → spherical_gauge_shift() + apply_spherical_gauge() in - spherical_functional.hpp — Newton + Backtracking - → Tests: ZerosSumGv, ApplyInPlace, AlreadyAtGauge +Phase 3f Analytische Hessians (Eucl. + Sphär.) ✅ abgeschlossen +Phase 3g PI-Konstante konsolidieren ✅ abgeschlossen -Phase 3f Analytische Hessians ✅ abgeschlossen - → euclidean_hessian.hpp: Kotangenten-Laplace (Pinkall–Polthier 1993) - mit korrektem 1/2-Normierungsfaktor; 8 Tests - → spherical_hessian.hpp: ∂α_i/∂u_j direkt aus sphärischem - Cosinussatz abgeleitet + Chain-Rule mit ∂l/∂λ = tan(l/2); 8 Tests - → Erkenntnis: Sphärische Energie ist konkav (H ist NSD), nicht konvex +Phase 4a Newton-Solver (alle drei Geometrien) ✅ abgeschlossen + → newton_euclidean / newton_spherical / newton_hyper_ideal + → detail::solve_with_fallback → public solve_linear_system + → Backtracking-Line-Search + → hyper_ideal_hessian.hpp (numerischer FD-Hessian) -Phase 3g PI-Konstante konsolidieren ✅ abgeschlossen - → constants.hpp mit conformallab::PI und TWO_PI - → 5 Dateien bereinigt; PI_SPHER-Alias rückwärtskompatibel +Phase 4b CGAL::IO Mesh-Import/Export ✅ abgeschlossen + → mesh_io.hpp: read/write/load/save + → Format-Erkennung aus Dateiendung (OFF, OBJ, PLY) -Phase 4a Minimaler Newton-Solver (1–2 Tage) - → Eigen SimplicialLDLT (bereits Projektabhängigkeit) - → Newton-Schritt: Δx = −H⁻¹·g mit Hessian aus 3f - → Kein externer PETSc-Solver nötig für erste Tests - → Konvergenzkriterium: ||g||∞ < tol +Phase 4c End-to-End-Pipeline Tests ✅ abgeschlossen + → test_pipeline.cpp: 5 Tests (alle 3 Geometrien, I/O, full loop) -Phase 4b CGAL::IO für Mesh-Import/Export (< 1 Tag) - → CGAL::IO::read_OBJ / write_OBJ auf Surface_mesh: - null Eigenentwicklung, sofort nutzbar - → Ermöglicht Round-Trip-Tests gegen Java-Referenz-Meshes - → Mittelfristig: CGAL::IO::read_OFF für .off-Dateien - (ersetzt Java-XML-Serialisierung für Testzwecke) +Phase 4d SparseQR-Fallback + Beispiel-Programme ✅ abgeschlossen + → solve_linear_system als öffentliche API mit fallback_used-Flag + → 3 dedizierte SparseQR-Tests (full-rank, singular, closed mesh) + → examples/example_euclidean.cpp (headless) + → examples/example_hyper_ideal.cpp (headless) + → examples/example_viewer.cpp (interaktiver Viewer, WITH_VIEWER) + +Phase 5 CLI-App + Serialisierung (2–3 Tage) + → conformallab_core CLI: --input, --geometry, --output + → JSON/XML-Export für DOF-Vektoren und Solver-Ergebnisse + → Analytischer HyperIdeal-Hessian (direkte Ableitung) + → Integration aller drei Geometrien hinter einheitlicher CLI ``` --- diff --git a/code/CMakeLists.txt b/code/CMakeLists.txt index 64f1a5f..6f6dc19 100644 --- a/code/CMakeLists.txt +++ b/code/CMakeLists.txt @@ -106,5 +106,10 @@ if(WITH_CGAL) RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/bin) endif() +# ── Example programs (require WITH_CGAL; viewer example also needs WITH_VIEWER) ─ +if(WITH_CGAL) + add_subdirectory(examples) +endif() + # ── Tests (always) ──────────────────────────────────────────────────────────── add_subdirectory(tests) diff --git a/code/examples/CMakeLists.txt b/code/examples/CMakeLists.txt new file mode 100644 index 0000000..a338815 --- /dev/null +++ b/code/examples/CMakeLists.txt @@ -0,0 +1,49 @@ +# examples/CMakeLists.txt +# +# Example programs for conformallab++. +# All examples require -DWITH_CGAL=ON. +# example_viewer additionally requires -DWITH_VIEWER=ON. +# +# Run after building: +# ./build/examples/example_euclidean [input.off] [output.off] +# ./build/examples/example_hyper_ideal [input.off] [output.off] +# ./build/examples/example_viewer [input.off] (requires WITH_VIEWER) + +# ── Shared include paths for all examples ───────────────────────────────────── +set(EXAMPLE_INCLUDES + ${CMAKE_SOURCE_DIR}/deps/eigen-3.4.0 + ${CMAKE_SOURCE_DIR}/deps/CGAL-6.1.1/include + ${Boost_INCLUDE_DIRS} +) +set(EXAMPLE_PRIVATE_INCLUDES + ${CMAKE_SOURCE_DIR}/include +) +set(EXAMPLE_DEFS + CGAL_DISABLE_GMP + CGAL_DISABLE_MPFR +) + +# ── example_euclidean ───────────────────────────────────────────────────────── +add_executable(example_euclidean example_euclidean.cpp) +target_include_directories(example_euclidean SYSTEM PRIVATE ${EXAMPLE_INCLUDES}) +target_include_directories(example_euclidean PRIVATE ${EXAMPLE_PRIVATE_INCLUDES}) +target_compile_definitions(example_euclidean PRIVATE ${EXAMPLE_DEFS}) + +# ── example_hyper_ideal ─────────────────────────────────────────────────────── +add_executable(example_hyper_ideal example_hyper_ideal.cpp) +target_include_directories(example_hyper_ideal SYSTEM PRIVATE ${EXAMPLE_INCLUDES}) +target_include_directories(example_hyper_ideal PRIVATE ${EXAMPLE_PRIVATE_INCLUDES}) +target_compile_definitions(example_hyper_ideal PRIVATE ${EXAMPLE_DEFS}) + +# ── example_viewer (requires WITH_VIEWER) ───────────────────────────────────── +if(WITH_VIEWER) + add_executable(example_viewer example_viewer.cpp) + target_include_directories(example_viewer SYSTEM PRIVATE + ${EXAMPLE_INCLUDES} + ${CMAKE_SOURCE_DIR}/deps/libigl-2.6.0/include + ${CMAKE_SOURCE_DIR}/deps/libigl-glad/include + ) + target_include_directories(example_viewer PRIVATE ${EXAMPLE_PRIVATE_INCLUDES}) + target_compile_definitions(example_viewer PRIVATE ${EXAMPLE_DEFS}) + target_link_libraries(example_viewer PRIVATE viewer) +endif() diff --git a/code/examples/example_euclidean.cpp b/code/examples/example_euclidean.cpp new file mode 100644 index 0000000..5278b94 --- /dev/null +++ b/code/examples/example_euclidean.cpp @@ -0,0 +1,117 @@ +// example_euclidean.cpp +// +// conformallab++ — Euclidean discrete conformal map (headless example) +// +// This program demonstrates the full library pipeline for the EUCLIDEAN +// discrete conformal functional: +// +// 1. Load a triangle mesh from an OFF file +// 2. Set up the Euclidean functional maps +// 3. Pin one vertex (gauge fix for open surfaces) +// 4. Set target angles via "natural equilibrium" (x* = x_input) +// 5. Solve with Newton + backtracking line search +// 6. Print per-vertex conformal factors u_i = x[v_idx[v]] +// 7. Save the result mesh (same geometry, solver state printed) +// +// Build (requires -DWITH_CGAL=ON): +// cmake -S code -B build -DWITH_CGAL=ON +// cmake --build build --target example_euclidean +// ./build/examples/example_euclidean [input.off] [output.off] +// +// If no input file is given the built-in make_quad_strip() mesh is used. + +#include "conformal_mesh.hpp" +#include "mesh_builder.hpp" +#include "mesh_io.hpp" +#include "euclidean_functional.hpp" +#include "newton_solver.hpp" +#include +#include +#include + +using namespace conformallab; + +int main(int argc, char* argv[]) +{ + // ── Step 1: obtain mesh ─────────────────────────────────────────────── + ConformalMesh mesh; + std::string input_path = (argc > 1) ? argv[1] : ""; + std::string output_path = (argc > 2) ? argv[2] : "/tmp/conformallab_euclidean_out.off"; + + if (input_path.empty()) { + std::cout << "[example_euclidean] No input file given — using make_quad_strip().\n"; + mesh = make_quad_strip(); + } else { + std::cout << "[example_euclidean] Loading mesh from: " << input_path << "\n"; + try { mesh = load_mesh(input_path); } + catch (const std::exception& e) { + std::cerr << "Error loading mesh: " << e.what() << "\n"; + return 1; + } + } + + std::cout << "[example_euclidean] Mesh: " + << mesh.number_of_vertices() << " vertices, " + << mesh.number_of_faces() << " faces.\n"; + + // ── Step 2: set up functional maps ──────────────────────────────────── + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + + // ── Step 3: pin the first vertex (gauge fix) ────────────────────────── + auto vit = mesh.vertices().begin(); + Vertex_index v_pinned = *vit++; + maps.v_idx[v_pinned] = -1; // pinned: u[v_pinned] = 0 (fixed) + int idx = 0; + for (; vit != mesh.vertices().end(); ++vit) + maps.v_idx[*vit] = idx++; + const int n = idx; + + std::cout << "[example_euclidean] DOFs: " << n << " (1 vertex pinned).\n"; + + // ── Step 4: natural equilibrium — set theta_v = actual angle sum at x=0 ─ + // After this step x* = 0 is the equilibrium (no deformation). + // In a real application you would set theta_v = desired angle (e.g. 2π + // for flat disks, or the cone angles for a cone metric). + { + std::vector x0(static_cast(n), 0.0); + auto G0 = euclidean_gradient(mesh, x0, maps); + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv >= 0) maps.theta_v[v] -= G0[static_cast(iv)]; + } + } + + // ── Step 5: solve from a small perturbation to demonstrate Newton ───── + std::vector x0(static_cast(n), -0.05); + std::cout << "[example_euclidean] Solving Newton system…\n"; + auto result = newton_euclidean(mesh, x0, maps, /*tol=*/1e-9, /*max_iter=*/100); + + // ── Step 6: report ──────────────────────────────────────────────────── + if (result.converged) { + std::cout << "[example_euclidean] Converged in " << result.iterations + << " iterations. ||G||_inf = " << result.grad_inf_norm << "\n"; + } else { + std::cout << "[example_euclidean] Did NOT converge after " << result.iterations + << " iterations. ||G||_inf = " << result.grad_inf_norm << "\n"; + } + + std::cout << "[example_euclidean] Per-vertex conformal factors u_i:\n"; + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + double u = (iv >= 0) ? result.x[static_cast(iv)] : 0.0; + std::cout << " v" << v << " u = " << u; + if (iv < 0) std::cout << " (pinned)"; + std::cout << "\n"; + } + + // ── Step 7: write output mesh ───────────────────────────────────────── + try { + save_mesh(output_path, mesh); + std::cout << "[example_euclidean] Mesh saved to: " << output_path << "\n"; + } catch (const std::exception& e) { + std::cerr << "Warning: could not write output: " << e.what() << "\n"; + } + + return result.converged ? 0 : 1; +} diff --git a/code/examples/example_hyper_ideal.cpp b/code/examples/example_hyper_ideal.cpp new file mode 100644 index 0000000..fcaea37 --- /dev/null +++ b/code/examples/example_hyper_ideal.cpp @@ -0,0 +1,147 @@ +// example_hyper_ideal.cpp +// +// conformallab++ — Hyper-ideal discrete conformal map (headless example) +// +// Demonstrates the full library pipeline for the HYPER-IDEAL discrete conformal +// functional (Springborn 2020). The hyper-ideal functional operates in +// hyperbolic geometry: vertices have "horoball radii" (DOF b_i) and edges have +// "intersection lengths" (DOF a_e). The energy is strictly convex, so Newton +// converges globally from any valid starting point. +// +// Pipeline: +// 1. Load (or synthesise) a triangle mesh +// 2. Set up HyperIdeal maps + assign all vertex and edge DOFs +// 3. Choose equilibrium base point (b=1.0, a=0.5) and set natural targets +// 4. Perturb and solve with Newton +// 5. Print DOF values at equilibrium +// 6. Save result mesh +// +// Build (requires -DWITH_CGAL=ON): +// cmake -S code -B build -DWITH_CGAL=ON +// cmake --build build --target example_hyper_ideal +// ./build/examples/example_hyper_ideal [input.off] [output.off] + +#include "conformal_mesh.hpp" +#include "mesh_builder.hpp" +#include "mesh_io.hpp" +#include "hyper_ideal_functional.hpp" +#include "newton_solver.hpp" +#include +#include +#include +#include + +using namespace conformallab; + +int main(int argc, char* argv[]) +{ + // ── Step 1: obtain mesh ─────────────────────────────────────────────── + ConformalMesh mesh; + std::string input_path = (argc > 1) ? argv[1] : ""; + std::string output_path = (argc > 2) ? argv[2] : "/tmp/conformallab_hyper_ideal_out.off"; + + if (input_path.empty()) { + std::cout << "[example_hyper_ideal] No input file — using make_triangle().\n"; + mesh = make_triangle(); + } else { + std::cout << "[example_hyper_ideal] Loading mesh from: " << input_path << "\n"; + try { mesh = load_mesh(input_path); } + catch (const std::exception& e) { + std::cerr << "Error loading mesh: " << e.what() << "\n"; + return 1; + } + } + + std::cout << "[example_hyper_ideal] Mesh: " + << mesh.number_of_vertices() << " vertices, " + << mesh.number_of_faces() << " faces.\n"; + + // ── Step 2: set up functional maps ──────────────────────────────────── + auto maps = setup_hyper_ideal_maps(mesh); + int n = assign_all_dof_indices(mesh, maps); + + std::cout << "[example_hyper_ideal] DOFs: " << n + << " (" << mesh.number_of_vertices() << " vertex + " + << mesh.number_of_edges() << " edge).\n"; + + // ── Step 3: choose equilibrium base point and set natural targets ───── + // + // x = 0 is degenerate for the HyperIdeal functional (log-space). + // We pick a valid base point (b_i = b_base, a_e = a_base), evaluate + // the gradient there, and absorb it into the target angles so that + // G(xbase) = 0. This makes xbase the equilibrium x*. + // + // In a real application you would set theta_v / theta_e to the desired + // hyperbolic angle targets (e.g. from a reference mesh). + const double b_base = 1.0; // horoball radii at equilibrium + const double a_base = 0.5; // edge-length DOFs at equilibrium + + const auto sz = static_cast(n); + std::vector xbase(sz, 0.0); + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv >= 0) xbase[static_cast(iv)] = b_base; + } + for (auto e : mesh.edges()) { + int ie = maps.e_idx[e]; + if (ie >= 0) xbase[static_cast(ie)] = a_base; + } + + // G = Σβ − theta_target; absorb G(xbase) into targets so G(xbase) = 0 + auto G0 = evaluate_hyper_ideal(mesh, xbase, maps, false).gradient; + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv >= 0) maps.theta_v[v] += G0[static_cast(iv)]; + } + for (auto e : mesh.edges()) { + int ie = maps.e_idx[e]; + if (ie >= 0) maps.theta_e[e] += G0[static_cast(ie)]; + } + + // ── Step 4: perturb and solve ───────────────────────────────────────── + const double perturb = 0.25; + std::vector x0 = xbase; + for (auto& v : x0) v += perturb; + + double g_start = 0.0; + for (double v : G0) g_start = std::max(g_start, std::abs(v)); + std::cout << "[example_hyper_ideal] Starting Newton from perturbation +" << perturb + << " (G at xbase = " << g_start << ").\n"; + + auto result = newton_hyper_ideal(mesh, x0, maps, /*tol=*/1e-9, /*max_iter=*/200); + + // ── Step 5: report ──────────────────────────────────────────────────── + if (result.converged) { + std::cout << "[example_hyper_ideal] Converged in " << result.iterations + << " iterations. ||G||_inf = " << result.grad_inf_norm << "\n"; + } else { + std::cout << "[example_hyper_ideal] Did NOT converge after " << result.iterations + << " iterations. ||G||_inf = " << result.grad_inf_norm << "\n"; + } + + std::cout << "[example_hyper_ideal] DOF values at equilibrium:\n"; + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv < 0) continue; + std::cout << " v" << v + << " b = " << result.x[static_cast(iv)] + << " (expected " << b_base << ")\n"; + } + for (auto e : mesh.edges()) { + int ie = maps.e_idx[e]; + if (ie < 0) continue; + std::cout << " e" << e + << " a = " << result.x[static_cast(ie)] + << " (expected " << a_base << ")\n"; + } + + // ── Step 6: write output mesh ───────────────────────────────────────── + try { + save_mesh(output_path, mesh); + std::cout << "[example_hyper_ideal] Mesh saved to: " << output_path << "\n"; + } catch (const std::exception& e) { + std::cerr << "Warning: could not write output: " << e.what() << "\n"; + } + + return result.converged ? 0 : 1; +} diff --git a/code/examples/example_viewer.cpp b/code/examples/example_viewer.cpp new file mode 100644 index 0000000..71ab7c7 --- /dev/null +++ b/code/examples/example_viewer.cpp @@ -0,0 +1,150 @@ +// example_viewer.cpp +// +// conformallab++ — Interactive viewer example (requires -DWITH_VIEWER=ON) +// +// This example demonstrates the full end-to-end pipeline WITH visual output: +// +// 1. Load (or synthesise) a mesh +// 2. Compute the Euclidean discrete conformal map (Newton solver) +// 3. Display the result in an interactive libigl / GLFW window +// • Left pane: input mesh, coloured by per-vertex conformal factor u_i +// • Right pane: a flat parameterisation (future, placeholder) +// +// The viewer is split into two data sets using libigl's multi-mesh API so +// users can inspect geometry and solution simultaneously. +// +// Build (requires -DWITH_CGAL=ON -DWITH_VIEWER=ON): +// cmake -S code -B build -DWITH_CGAL=ON -DWITH_VIEWER=ON +// cmake --build build --target example_viewer +// ./build/examples/example_viewer [input.off] +// +// Navigation (libigl default): +// Mouse drag — rotate +// Scroll — zoom +// C — toggle camera mode (trackball / 2D) +// Z / X / Y — snap to axis-aligned view +// Q / Esc — quit + +#include "conformal_mesh.hpp" +#include "mesh_builder.hpp" +#include "mesh_io.hpp" +#include "euclidean_functional.hpp" +#include "mesh_utils.hpp" +#include "newton_solver.hpp" +#include +#include +#include +#include +#include +#include +#include + +using namespace conformallab; + +// ── Jet colour map: scalar → RGB ───────────────────────────────────────────── +static Eigen::RowVector3d jet(double t) +{ + t = std::max(0.0, std::min(1.0, t)); + double r = std::clamp(1.5 - std::abs(4.0 * t - 3.0), 0.0, 1.0); + double g = std::clamp(1.5 - std::abs(4.0 * t - 2.0), 0.0, 1.0); + double b = std::clamp(1.5 - std::abs(4.0 * t - 1.0), 0.0, 1.0); + return {r, g, b}; +} + +int main(int argc, char* argv[]) +{ + // ── Step 1: load or synthesise mesh ─────────────────────────────────── + ConformalMesh mesh; + std::string input_path = (argc > 1) ? argv[1] : ""; + + if (input_path.empty()) { + std::cout << "[example_viewer] No input file — using make_quad_strip().\n"; + mesh = make_quad_strip(); + } else { + std::cout << "[example_viewer] Loading: " << input_path << "\n"; + try { mesh = load_mesh(input_path); } + catch (const std::exception& e) { + std::cerr << "Error: " << e.what() << "\n"; + return 1; + } + } + + std::cout << "[example_viewer] Mesh: " + << mesh.number_of_vertices() << " vertices, " + << mesh.number_of_faces() << " faces.\n"; + + // ── Step 2: solve the Euclidean discrete conformal map ──────────────── + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + + // Pin first vertex + auto vit = mesh.vertices().begin(); + maps.v_idx[*vit++] = -1; + int idx = 0; + for (; vit != mesh.vertices().end(); ++vit) + maps.v_idx[*vit] = idx++; + const int n = idx; + + // Natural equilibrium + { + std::vector x0(static_cast(n), 0.0); + auto G0 = euclidean_gradient(mesh, x0, maps); + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv >= 0) maps.theta_v[v] -= G0[static_cast(iv)]; + } + } + + // Perturb and solve + std::vector x0(static_cast(n), -0.08); + auto result = newton_euclidean(mesh, x0, maps, 1e-9, 200); + + if (result.converged) + std::cout << "[example_viewer] Converged in " << result.iterations << " iterations.\n"; + else + std::cout << "[example_viewer] Warning: did not converge fully " + "(||G||_inf = " << result.grad_inf_norm << ").\n"; + + // ── Step 3: build Eigen V / F for libigl ───────────────────────────── + using Kernel = CGAL::Simple_cartesian; + Eigen::MatrixXd V; + Eigen::MatrixXi F; + mesh_utils::cgal_to_eigen(mesh, V, F); + + // Per-vertex colour: conformal factor u_i, mapped via jet palette + const int nv = static_cast(V.rows()); + Eigen::MatrixXd C(nv, 3); + + // Collect all u values to normalise + std::vector u_all(static_cast(nv), 0.0); + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv >= 0) + u_all[static_cast(v.idx())] = result.x[static_cast(iv)]; + } + double u_min = *std::min_element(u_all.begin(), u_all.end()); + double u_max = *std::max_element(u_all.begin(), u_all.end()); + double u_range = (u_max > u_min) ? (u_max - u_min) : 1.0; + + for (int vi = 0; vi < nv; ++vi) { + double t = (u_all[static_cast(vi)] - u_min) / u_range; + C.row(vi) = jet(t); + } + + // ── Step 4: launch interactive viewer ───────────────────────────────── + igl::opengl::glfw::Viewer viewer; + viewer.data().set_mesh(V, F); + viewer.data().set_colors(C); + viewer.data().show_lines = true; + viewer.data().show_overlay = true; + + // Status text overlay + viewer.data().add_label( + Eigen::Vector3d(V.col(0).mean(), V.col(1).mean(), V.col(2).maxCoeff()), + "Euclidean conformal factor u_i (jet: blue=min, red=max)"); + + std::cout << "[example_viewer] Launching viewer. Press Q or Esc to quit.\n"; + viewer.launch(); + + return 0; +} diff --git a/code/include/hyper_ideal_hessian.hpp b/code/include/hyper_ideal_hessian.hpp new file mode 100644 index 0000000..803b9ab --- /dev/null +++ b/code/include/hyper_ideal_hessian.hpp @@ -0,0 +1,87 @@ +#pragma once +// hyper_ideal_hessian.hpp +// +// Phase 4a — Hessian of the hyper-ideal discrete conformal functional. +// +// ┌──────────────────────────────────────────────────────────────────────────┐ +// │ Implementation strategy │ +// │ │ +// │ The hyper-ideal functional involves angle functions (ζ, σ, α, β) │ +// │ composed through several nested layers (lij → ζ13/14/15 → β/α). │ +// │ Deriving closed-form Hessian entries analytically through all these │ +// │ layers is feasible but lengthy; an analytical Hessian is left for a │ +// │ future phase. │ +// │ │ +// │ Here we compute the Hessian by symmetric finite differences of the │ +// │ gradient, which is exact to O(ε²) and sufficient for Newton's method │ +// │ at the meshes typical in Phase 4 (< 500 DOFs): │ +// │ │ +// │ H[i,j] = (G(x + ε·eⱼ)[i] − G(x − ε·eⱼ)[i]) / (2ε) │ +// │ │ +// │ The hyper-ideal energy is strictly convex (Springborn 2020), so H is │ +// │ positive semi-definite everywhere and Eigen::SimplicialLDLT applies │ +// │ directly. │ +// └──────────────────────────────────────────────────────────────────────────┘ + +#include "hyper_ideal_functional.hpp" +#include +#include +#include + +namespace conformallab { + +// ── Numerical Hessian via symmetric finite differences ──────────────────────── +// +// 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). +inline Eigen::SparseMatrix hyper_ideal_hessian( + ConformalMesh& mesh, + const std::vector& x, + const HyperIdealMaps& m, + double eps = 1e-5) +{ + const int n = hyper_ideal_dimension(mesh, m); + std::vector> trips; + trips.reserve(static_cast(n * n)); // dense upper bound + + std::vector xp = x, xm = x; + + for (int j = 0; j < n; ++j) { + const std::size_t sj = static_cast(j); + xp[sj] = x[sj] + eps; + xm[sj] = x[sj] - eps; + + auto Gp = evaluate_hyper_ideal(mesh, xp, m, /*energy=*/false).gradient; + auto Gm = evaluate_hyper_ideal(mesh, xm, m, /*energy=*/false).gradient; + + xp[sj] = xm[sj] = x[sj]; // restore + + for (int i = 0; i < n; ++i) { + double val = (Gp[static_cast(i)] + - Gm[static_cast(i)]) / (2.0 * eps); + if (std::abs(val) > 1e-15) + trips.emplace_back(i, j, val); + } + } + + Eigen::SparseMatrix H(n, n); + H.setFromTriplets(trips.begin(), trips.end()); + return H; +} + +// ── Symmetrised Hessian ─────────────────────────────────────────────────────── +// +// The FD Hessian is symmetric in exact arithmetic; floating-point rounding +// can introduce tiny asymmetries. This helper returns (H + Hᵀ)/2. +inline Eigen::SparseMatrix hyper_ideal_hessian_sym( + ConformalMesh& mesh, + const std::vector& x, + const HyperIdealMaps& m, + double eps = 1e-5) +{ + auto H = hyper_ideal_hessian(mesh, x, m, eps); + Eigen::SparseMatrix Ht = H.transpose(); + return (H + Ht) * 0.5; +} + +} // namespace conformallab diff --git a/code/include/newton_solver.hpp b/code/include/newton_solver.hpp index 2d82857..dcb222c 100644 --- a/code/include/newton_solver.hpp +++ b/code/include/newton_solver.hpp @@ -1,27 +1,39 @@ #pragma once // newton_solver.hpp // -// Phase 4a — Newton solver for the discrete conformal functionals. +// Phase 4a — Newton solver for all three discrete conformal functionals. // -// Solves G(x) = 0 where G is the gradient of the discrete conformal energy: -// G_v = Θ_v − Σ_f α_v^f (angle-sum residual at each vertex DOF) +// Solves G(x) = 0 where G is the gradient of the discrete conformal energy. // -// Algorithm per iteration: -// 1. Compute gradient G(x) -// 2. Check convergence: max|G_i| < tol → done -// 3. Compute sparse Hessian H(x) -// 4. Factorize and solve the Newton system: -// Euclidean: H · Δx = −G (H is PSD → SimplicialLDLT directly) -// Spherical: (−H) · Δx = G (H is NSD → negate to get PSD matrix) -// 5. Backtracking line search: halve α until ||G(x+α·Δx)|| < ||G(x)|| -// 6. x ← x + α·Δx, go to 1 +// ┌──────────────────────────────────────────────────────────────────────────┐ +// │ Gradient sign conventions │ +// │ Euclidean / Spherical: G_v = Θ_v − Σ α_v (target − actual) │ +// │ HyperIdeal: G_v = Σ β_v − Θ_v (actual − target) │ +// │ │ +// │ All solvers use the same Newton step Δx = −H⁻¹·G │ +// │ │ +// │ Hessian sign at equilibrium │ +// │ Euclidean: H PSD → SimplicialLDLT on H │ +// │ Spherical: H NSD → SimplicialLDLT on −H (solve (−H)Δx = G) │ +// │ HyperIdeal: H PSD → SimplicialLDLT on H (analytical H: future) │ +// └──────────────────────────────────────────────────────────────────────────┘ +// +// SparseQR fallback: +// When SimplicialLDLT reports a failure (e.g. singular H on a closed mesh +// without a pinned vertex), the solver automatically retries with +// Eigen::SparseQR, which finds the minimum-norm Newton step orthogonal to +// the null space. This handles the gauge mode on closed surfaces without +// requiring the caller to pin a vertex explicitly. // // Requires: -// Eigen::SimplicialLDLT (part of Eigen's sparse Cholesky module) +// Eigen::SimplicialLDLT, Eigen::SparseQR (Eigen sparse module) #include "euclidean_hessian.hpp" #include "spherical_hessian.hpp" +#include "hyper_ideal_hessian.hpp" #include +#include +#include #include #include #include @@ -41,6 +53,58 @@ struct NewtonResult { namespace detail { +// Solve A·Δx = rhs with SimplicialLDLT; on failure fall back to SparseQR. +// Returns Δx. ok is set to false only if both solvers fail. +// If fallback_used is non-null, it is set to true iff SparseQR was needed. +inline Eigen::VectorXd solve_with_fallback( + const Eigen::SparseMatrix& A, + const Eigen::VectorXd& rhs, + bool& ok, + bool* fallback_used = nullptr) +{ + if (fallback_used) *fallback_used = false; + + Eigen::SimplicialLDLT> ldlt(A); + if (ldlt.info() == Eigen::Success) { + Eigen::VectorXd dx = ldlt.solve(rhs); + if (ldlt.info() == Eigen::Success) { ok = true; return dx; } + } + // Fallback: SparseQR — handles singular/rank-deficient H (gauge modes). + if (fallback_used) *fallback_used = true; + Eigen::SparseQR, Eigen::COLAMDOrdering> qr(A); + if (qr.info() == Eigen::Success) { + Eigen::VectorXd dx = qr.solve(rhs); + if (qr.info() == Eigen::Success) { ok = true; return dx; } + } + ok = false; + return Eigen::VectorXd::Zero(rhs.size()); +} + +} // namespace detail + +// ── Public linear-system solver (SparseQR fallback) ────────────────────────── +// +// Solve A·x = rhs with Eigen::SimplicialLDLT; if that fails (singular or +// rank-deficient A), retry with Eigen::SparseQR which finds the minimum-norm +// solution orthogonal to the null space. +// +// This is the same primitive used internally by all three Newton solvers. +// Exposing it publicly lets callers (tests, downstream code) reuse the logic +// and — via the optional fallback_used pointer — verify which code path ran. +// +// fallback_used – if non-null, set to true iff SparseQR was invoked +// Returns Eigen::VectorXd::Zero(rhs.size()) if both solvers fail. +inline Eigen::VectorXd solve_linear_system( + const Eigen::SparseMatrix& A, + const Eigen::VectorXd& rhs, + bool* fallback_used = nullptr) +{ + bool ok = false; + return detail::solve_with_fallback(A, rhs, ok, fallback_used); +} + +namespace detail { // re-open for the remaining helpers + // Backtracking line search: find the largest α in {1, 0.5, 0.25, …} such that // ||G(x + α·Δx)||₂ < ||G(x)||₂. Returns the accepted step (α may stay 1). template @@ -109,14 +173,11 @@ inline NewtonResult newton_euclidean( return res; } - // ── Hessian + factorisation ─────────────────────────────────────────── + // ── Hessian + solve H·Δx = −G (SparseQR fallback for singular H) ── auto H = euclidean_hessian(mesh, x, m); - solver.compute(H); - if (solver.info() != Eigen::Success) break; - - // ── Newton step: solve H·Δx = −G ──────────────────────────────────── - Eigen::VectorXd dx = solver.solve(-G); - if (solver.info() != Eigen::Success) break; + bool ok = false; + Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok); + if (!ok) break; // ── Backtracking line search ────────────────────────────────────────── double norm0 = G.norm(); @@ -157,8 +218,6 @@ inline NewtonResult newton_spherical( res.iterations = 0; res.grad_inf_norm = 0.0; - Eigen::SimplicialLDLT> solver; - for (int iter = 0; iter < max_iter; ++iter) { // ── Gradient ────────────────────────────────────────────────────────── auto G_std = spherical_gradient(mesh, x, m); @@ -173,15 +232,12 @@ inline NewtonResult newton_spherical( return res; } - // ── Hessian: negate to get PSD matrix ──────────────────────────────── - auto H = spherical_hessian(mesh, x, m); - auto negH = Eigen::SparseMatrix(-H); - solver.compute(negH); - if (solver.info() != Eigen::Success) break; - - // ── Newton step: solve (−H)·Δx = G ───────────────────────────────── - Eigen::VectorXd dx = solver.solve(G); - if (solver.info() != Eigen::Success) break; + // ── Hessian: negate to get PSD; solve (−H)·Δx = G ────────────────── + auto H = spherical_hessian(mesh, x, m); + auto negH = Eigen::SparseMatrix(-H); + bool ok = false; + Eigen::VectorXd dx = detail::solve_with_fallback(negH, G, ok); + if (!ok) break; // ── Backtracking line search ────────────────────────────────────────── double norm0 = G.norm(); @@ -201,4 +257,70 @@ inline NewtonResult newton_spherical( return res; } +// ── HyperIdeal Newton solver ────────────────────────────────────────────────── +// +// Solves G(x) = 0 for the hyper-ideal discrete conformal functional. +// +// Gradient sign convention (opposite to Euclidean/Spherical): +// G_v = Σ β_v − Θ_v, G_e = Σ α_e − θ_e (actual − target) +// +// The hyper-ideal energy is strictly convex (Springborn 2020), so H is PSD +// and SimplicialLDLT (with SparseQR fallback) applies directly. +// +// The Hessian is computed by symmetric finite differences of G (see +// hyper_ideal_hessian.hpp); replace with an analytical Hessian in Phase 5. +inline NewtonResult newton_hyper_ideal( + ConformalMesh& mesh, + std::vector x0, + const HyperIdealMaps& m, + double tol = 1e-8, + int max_iter = 200, + double hess_eps = 1e-5) +{ + std::vector x = x0; + const int n = static_cast(x.size()); + + NewtonResult res; + res.converged = false; + res.iterations = 0; + res.grad_inf_norm = 0.0; + + for (int iter = 0; iter < max_iter; ++iter) { + // ── Gradient ────────────────────────────────────────────────────────── + auto G_std = evaluate_hyper_ideal(mesh, x, m, /*energy=*/false).gradient; + Eigen::Map G(G_std.data(), n); + + double inf_norm = G.cwiseAbs().maxCoeff(); + if (inf_norm < tol) { + res.converged = true; + res.grad_inf_norm = inf_norm; + res.iterations = iter; + res.x = x; + return res; + } + + // ── Hessian (numerical FD) + solve H·Δx = −G ───────────────────────── + auto H = hyper_ideal_hessian_sym(mesh, x, m, hess_eps); + bool ok = false; + Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok); + if (!ok) break; + + // ── Backtracking line search ────────────────────────────────────────── + double norm0 = G.norm(); + x = detail::line_search(x, dx, norm0, + [&](const std::vector& xnew) { + return evaluate_hyper_ideal(mesh, xnew, m, false).gradient; + }); + + res.iterations = iter + 1; + } + + auto G_final = evaluate_hyper_ideal(mesh, x, m, false).gradient; + double inf_final = 0.0; + for (double v : G_final) inf_final = std::max(inf_final, std::abs(v)); + res.grad_inf_norm = inf_final; + res.x = x; + return res; +} + } // namespace conformallab diff --git a/code/tests/cgal/CMakeLists.txt b/code/tests/cgal/CMakeLists.txt index 656d714..fa13123 100644 --- a/code/tests/cgal/CMakeLists.txt +++ b/code/tests/cgal/CMakeLists.txt @@ -30,6 +30,9 @@ add_executable(conformallab_cgal_tests # ── Phase 4b: Mesh I/O (CGAL::IO) ───────────────────────────────────── test_mesh_io.cpp + + # ── Phase 4c: End-to-end pipeline + user examples ───────────────────── + test_pipeline.cpp ) target_include_directories(conformallab_cgal_tests SYSTEM PRIVATE diff --git a/code/tests/cgal/test_hyper_ideal_functional.cpp b/code/tests/cgal/test_hyper_ideal_functional.cpp index 22da339..5e6e3f8 100644 --- a/code/tests/cgal/test_hyper_ideal_functional.cpp +++ b/code/tests/cgal/test_hyper_ideal_functional.cpp @@ -20,7 +20,9 @@ #include "conformal_mesh.hpp" #include "mesh_builder.hpp" #include "hyper_ideal_functional.hpp" +#include "hyper_ideal_hessian.hpp" #include +#include #include #include @@ -52,9 +54,19 @@ static std::vector make_x_all_variable( // @Ignore in Java: no Hessian implemented // ════════════════════════════════════════════════════════════════════════════ -TEST(HyperIdealFunctional, GradientCheck_Hessian) +TEST(HyperIdealFunctional, HessianSymmetryCheck) { - GTEST_SKIP() << "@Ignore in Java – Hessian not implemented in the functional"; + // Hessian is now implemented (numerical FD). Verify it is symmetric. + auto mesh = make_triangle(); + auto maps = setup_hyper_ideal_maps(mesh); + int n = assign_all_dof_indices(mesh, maps); + + std::vector x(static_cast(n), 0.5); + auto H = hyper_ideal_hessian_sym(mesh, x, maps); + + Eigen::MatrixXd Hd(H); + EXPECT_NEAR((Hd - Hd.transpose()).norm(), 0.0, 1e-8) + << "HyperIdeal Hessian must be symmetric"; } // ════════════════════════════════════════════════════════════════════════════ diff --git a/code/tests/cgal/test_newton_solver.cpp b/code/tests/cgal/test_newton_solver.cpp index 4475bb1..dc80e54 100644 --- a/code/tests/cgal/test_newton_solver.cpp +++ b/code/tests/cgal/test_newton_solver.cpp @@ -1,12 +1,14 @@ // test_newton_solver.cpp // -// Phase 4a — Newton solver tests. +// Phase 4 — Newton solver tests. // // Design principle: // We test convergence to a KNOWN equilibrium. For the spherical tetrahedron // x* = 0 is built-in (G(0) ≈ 0 by construction). For Euclidean meshes we // use "natural theta": set theta_v[v] = actual angle sum at x=0, which makes // x* = 0 the exact equilibrium by definition. +// For HyperIdeal we use the same "natural target" trick at a valid base point +// (b=1.0, a=0.5), since x=0 is degenerate for the HyperIdeal functional. // // Tests: // Spherical: @@ -19,12 +21,26 @@ // 5. Converges (triangle, 1 pinned vertex, natural theta). // 6. Converges (quad strip, 1 pinned vertex, natural theta). // 7. Converges with explicitly chosen mixed pinned/variable layout. +// +// HyperIdeal: +// 8. Converges on triangle (all variable, natural targets). +// 9. Result fields self-consistent. +// 10. Converges on tetrahedron (10 DOFs, larger mesh). +// 11. SparseQR: result consistent (valid starting region). +// +// SparseQR fallback (direct unit tests): +// 12. solve_linear_system recovers correct solution on rank-deficient matrix. +// 13. solve_linear_system sets fallback_used=true on a singular matrix. +// 14. Euclidean Newton on closed tetrahedron (no pinned vertex) converges +// via SparseQR gauge-mode handling. #include "conformal_mesh.hpp" #include "mesh_builder.hpp" #include "euclidean_functional.hpp" #include "spherical_functional.hpp" +#include "hyper_ideal_functional.hpp" #include "newton_solver.hpp" +#include #include #include #include @@ -233,3 +249,252 @@ TEST(NewtonSolver, Euclidean_ConvergesMixedPinned) "grad_inf_norm = " << res.grad_inf_norm; EXPECT_LT(res.grad_inf_norm, 1e-8); } + +// ════════════════════════════════════════════════════════════════════════════ +// Helper: set HyperIdeal target angles to actual sums at a non-degenerate +// base point (b_base, a_base), making that point the equilibrium x*. +// +// Note: x = 0 is degenerate for the HyperIdeal functional (log-space; the +// functional requires b_i > 0 / a_e > 0). We therefore choose a valid base +// point, evaluate G there, and absorb G into the targets so that G(xbase) = 0. +// Newton tests then start from a perturbation of xbase. +// +// Returns xbase so callers can construct a perturbed starting point. +// ════════════════════════════════════════════════════════════════════════════ +static std::vector set_natural_hyper_ideal_targets( + ConformalMesh& mesh, HyperIdealMaps& maps, int n, + double b_base = 1.0, double a_base = 0.5) +{ + const auto sz = static_cast(n); + std::vector xbase(sz, 0.0); + + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv >= 0) xbase[static_cast(iv)] = b_base; + } + for (auto e : mesh.edges()) { + int ie = maps.e_idx[e]; + if (ie >= 0) xbase[static_cast(ie)] = a_base; + } + + // G = Σβ - theta_target (initial target = 0 → G = Σβ = "actual" angles) + auto G = evaluate_hyper_ideal(mesh, xbase, maps, /*energy=*/false).gradient; + + // Set target := actual so that G(xbase) = actual - target = 0 + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv < 0) continue; + maps.theta_v[v] += G[static_cast(iv)]; + } + for (auto e : mesh.edges()) { + int ie = maps.e_idx[e]; + if (ie < 0) continue; + maps.theta_e[e] += G[static_cast(ie)]; + } + + return xbase; +} + +// ════════════════════════════════════════════════════════════════════════════ +// HyperIdeal 1 — Triangle, all DOFs variable, converges from perturbation +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonSolver, HyperIdeal_ConvergesTriangleAllVariable) +{ + auto mesh = make_triangle(); + auto maps = setup_hyper_ideal_maps(mesh); + int n = assign_all_dof_indices(mesh, maps); + + // xbase = (b=1.0, a=0.5) is the equilibrium after natural-target setup. + auto xbase = set_natural_hyper_ideal_targets(mesh, maps, n); + + // Perturb by +0.2 uniformly + std::vector x0 = xbase; + for (auto& v : x0) v += 0.2; + + auto res = newton_hyper_ideal(mesh, x0, maps, /*tol=*/1e-7, /*max_iter=*/100); + + EXPECT_TRUE(res.converged) + << "Newton (HyperIdeal, triangle) should converge; " + "grad_inf_norm = " << res.grad_inf_norm; + EXPECT_LT(res.grad_inf_norm, 1e-7); +} + +// ════════════════════════════════════════════════════════════════════════════ +// HyperIdeal 2 — Result fields self-consistent +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonSolver, HyperIdeal_ResultFieldsConsistent) +{ + auto mesh = make_triangle(); + auto maps = setup_hyper_ideal_maps(mesh); + int n = assign_all_dof_indices(mesh, maps); + + auto xbase = set_natural_hyper_ideal_targets(mesh, maps, n); + + std::vector x0 = xbase; + for (auto& v : x0) v += 0.1; + + auto res = newton_hyper_ideal(mesh, x0, maps, /*tol=*/1e-7, /*max_iter=*/100); + + EXPECT_EQ(static_cast(res.x.size()), n); + + // Reported grad_inf_norm must match re-computed gradient at res.x + auto G = evaluate_hyper_ideal(mesh, res.x, maps, false).gradient; + double actual_inf = 0.0; + for (double v : G) actual_inf = std::max(actual_inf, std::abs(v)); + EXPECT_NEAR(actual_inf, res.grad_inf_norm, 1e-9); +} + +// ════════════════════════════════════════════════════════════════════════════ +// HyperIdeal 3 — Tetrahedron (10 DOFs): 4 vertex b-vals + 6 edge a-vals +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonSolver, HyperIdeal_ConvergesTetrahedron) +{ + auto mesh = make_tetrahedron(); + auto maps = setup_hyper_ideal_maps(mesh); + int n = assign_all_dof_indices(mesh, maps); + + auto xbase = set_natural_hyper_ideal_targets(mesh, maps, n); + + // Perturb by +0.15 + std::vector x0 = xbase; + for (auto& v : x0) v += 0.15; + + auto res = newton_hyper_ideal(mesh, x0, maps, /*tol=*/1e-7, /*max_iter=*/200); + + EXPECT_TRUE(res.converged) + << "Newton (HyperIdeal, tetrahedron) should converge; " + "grad_inf_norm = " << res.grad_inf_norm; + EXPECT_LT(res.grad_inf_norm, 1e-7); +} + +// ════════════════════════════════════════════════════════════════════════════ +// HyperIdeal 4 — SparseQR fallback: solver returns a result (no crash) +// +// With all targets = 0 the equilibrium is not at x=0 but the solver should +// at minimum not crash and return a consistent result struct. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonSolver, HyperIdeal_SparseQRFallbackNoCrash) +{ + auto mesh = make_triangle(); + auto maps = setup_hyper_ideal_maps(mesh); + int n = assign_all_dof_indices(mesh, maps); + + // Leave targets at their default (0): solver tries to solve but the + // "equilibrium" is at some unknown x*. With valid starting point the + // Hessian is positive-definite and the solver should not crash. + // We don't assert convergence — just that the result struct is consistent. + std::vector x0(static_cast(n), 1.0); + // Mix vertex / edge DOFs: b=1.0, a=0.5 (valid region of the functional) + for (auto e : mesh.edges()) { + int ie = maps.e_idx[e]; + if (ie >= 0) x0[static_cast(ie)] = 0.5; + } + auto res = newton_hyper_ideal(mesh, x0, maps, /*tol=*/1e-7, /*max_iter=*/50); + + // Struct fields must always be populated + EXPECT_EQ(static_cast(res.x.size()), n); + EXPECT_GE(res.iterations, 0); + EXPECT_FALSE(std::isnan(res.grad_inf_norm)); + EXPECT_FALSE(std::isinf(res.grad_inf_norm)); +} + +// ════════════════════════════════════════════════════════════════════════════ +// SparseQR fallback — Test 12: solve_linear_system recovers correct solution +// +// The public API solve_linear_system(A, rhs) must return the correct answer +// for a well-conditioned full-rank system (LDLT path taken). +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SparseQRFallback, FullRankSystem_CorrectSolution) +{ + // Build a simple 3×3 diagonal PD matrix: A = diag(1, 2, 3) + Eigen::SparseMatrix A(3, 3); + A.insert(0, 0) = 1.0; + A.insert(1, 1) = 2.0; + A.insert(2, 2) = 3.0; + A.makeCompressed(); + + Eigen::VectorXd rhs(3); + rhs << 1.0, 4.0, 9.0; // solution = [1, 2, 3] + + bool fallback = true; // expect it to be set to false (LDLT succeeds) + Eigen::VectorXd x = conformallab::solve_linear_system(A, rhs, &fallback); + + EXPECT_FALSE(fallback) << "Full-rank system: LDLT should succeed (no SparseQR needed)"; + EXPECT_NEAR(x[0], 1.0, 1e-12); + EXPECT_NEAR(x[1], 2.0, 1e-12); + EXPECT_NEAR(x[2], 3.0, 1e-12); +} + +// ════════════════════════════════════════════════════════════════════════════ +// SparseQR fallback — Test 13: fallback_used=true on a singular matrix +// +// Construct a symmetric 3×3 matrix of rank 1 where LDLT fails (the (2,2) +// pivot is zero). SparseQR finds the minimum-norm least-squares solution. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SparseQRFallback, SingularMatrix_FallbackActivated) +{ + // A = [[2, 0, 0], + // [0, 0, 0], ← zero pivot → LDLT failure + // [0, 0, 3]] + // rhs compatible with the row space: [2, 0, 3] → solution [1, 0, 1] + Eigen::SparseMatrix A(3, 3); + A.insert(0, 0) = 2.0; + // row/col 1 deliberately all-zero + A.insert(2, 2) = 3.0; + A.makeCompressed(); + + Eigen::VectorXd rhs(3); + rhs << 2.0, 0.0, 3.0; + + bool fallback = false; + Eigen::VectorXd x = conformallab::solve_linear_system(A, rhs, &fallback); + + EXPECT_TRUE(fallback) << "Singular matrix: SparseQR fallback must be triggered"; + // SparseQR min-norm solution: x[0]=1, x[1]=0, x[2]=1 + EXPECT_NEAR(x[0], 1.0, 1e-10); + EXPECT_NEAR(x[1], 0.0, 1e-10); + EXPECT_NEAR(x[2], 1.0, 1e-10); +} + +// ════════════════════════════════════════════════════════════════════════════ +// SparseQR fallback — Test 14: Euclidean Newton on a closed mesh, no pinning +// +// make_tetrahedron() is a closed surface (4 vertices, 4 faces). Without a +// pinned vertex the Euclidean Hessian has a 1-D null space (uniform scale +// gauge mode): H·1 = 0. SimplicialLDLT fails on this rank-deficient H; +// SparseQR finds the min-norm Newton step orthogonal to the null space. +// +// The gradient always lives in the row space of H (Σ G_v = 0 by angle-sum +// invariance), so the SparseQR step is also the Newton step and the solver +// converges to the natural equilibrium. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SparseQRFallback, Euclidean_ClosedMeshNoPinConverges) +{ + auto mesh = make_tetrahedron(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + + // Assign all 4 vertices as free DOFs (no pinning). + int idx = 0; + for (auto v : mesh.vertices()) + maps.v_idx[v] = idx++; + const int n = idx; // = 4 + + // Natural theta: equilibrium at x* = 0. + set_natural_euclidean_theta(mesh, maps, n); + + std::vector x0(static_cast(n), -0.1); + auto res = newton_euclidean(mesh, x0, maps, /*tol=*/1e-8, /*max_iter=*/100); + + EXPECT_TRUE(res.converged) + << "Euclidean Newton on closed tetrahedron (no pin) must converge via SparseQR; " + "grad_inf_norm = " << res.grad_inf_norm; + EXPECT_LT(res.grad_inf_norm, 1e-8); +} diff --git a/code/tests/cgal/test_pipeline.cpp b/code/tests/cgal/test_pipeline.cpp new file mode 100644 index 0000000..17c3109 --- /dev/null +++ b/code/tests/cgal/test_pipeline.cpp @@ -0,0 +1,292 @@ +// test_pipeline.cpp +// +// Phase 4c — End-to-end pipeline tests and library-user examples. +// +// These tests exercise the full conformallab++ pipeline as a user would: +// +// 1. Build (or load) a mesh +// 2. Set up maps and assign DOFs +// 3. Configure target angles / targets +// 4. Solve with Newton +// 5. Inspect / export the result +// +// Each test mirrors a realistic usage scenario documented in the README. +// +// Tests: +// 1. Pipeline_Euclidean_TriangleToEquilibrium +// Read mesh → setup Euclidean maps → solve → verify convergence +// 2. Pipeline_Spherical_TetrahedronToEquilibrium +// Setup spherical tetrahedron → solve → verify angles sum to 4π +// 3. Pipeline_HyperIdeal_TriangleRoundTrip +// Build triangle → setup HyperIdeal → solve → verify G ≈ 0 +// 4. Pipeline_MeshIO_SolveAndExport +// Build mesh → solve → write OFF → reload → verify vertex count intact +// 5. Pipeline_AllThreeGeometries_SameTopology +// Same quad-strip mesh solved under all three geometries: all converge + +#include "conformal_mesh.hpp" +#include "mesh_builder.hpp" +#include "mesh_io.hpp" +#include "euclidean_functional.hpp" +#include "spherical_functional.hpp" +#include "hyper_ideal_functional.hpp" +#include "newton_solver.hpp" +#include +#include +#include +#include + +using namespace conformallab; + +// ──────────────────────────────────────────────────────────────────────────── +// Shared helpers +// ──────────────────────────────────────────────────────────────────────────── + +static void pin_first_vertex_euclidean(ConformalMesh& mesh, EuclideanMaps& maps, int& n) +{ + auto vit = mesh.vertices().begin(); + Vertex_index v0 = *vit++; + maps.v_idx[v0] = -1; + int idx = 0; + for (; vit != mesh.vertices().end(); ++vit) + maps.v_idx[*vit] = idx++; + n = idx; +} + +static void set_natural_euclidean_theta(ConformalMesh& mesh, EuclideanMaps& maps, int n) +{ + std::vector x0(static_cast(n), 0.0); + auto G = euclidean_gradient(mesh, x0, maps); + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv < 0) continue; + maps.theta_v[v] -= G[static_cast(iv)]; + } +} + +static std::vector set_natural_hyper_ideal_targets( + ConformalMesh& mesh, HyperIdealMaps& maps, int n, + double b_base = 1.0, double a_base = 0.5) +{ + const auto sz = static_cast(n); + std::vector xbase(sz, 0.0); + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv >= 0) xbase[static_cast(iv)] = b_base; + } + for (auto e : mesh.edges()) { + int ie = maps.e_idx[e]; + if (ie >= 0) xbase[static_cast(ie)] = a_base; + } + auto G = evaluate_hyper_ideal(mesh, xbase, maps, false).gradient; + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv < 0) continue; + maps.theta_v[v] += G[static_cast(iv)]; + } + for (auto e : mesh.edges()) { + int ie = maps.e_idx[e]; + if (ie < 0) continue; + maps.theta_e[e] += G[static_cast(ie)]; + } + return xbase; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Test 1 — Euclidean full pipeline: triangle → equilibrium +// +// Simulates a user doing: +// auto mesh = make_triangle(); +// auto maps = setup_euclidean_maps(mesh); +// compute_euclidean_lambda0_from_mesh(mesh, maps); +// // … set target angles and DOF indices … +// auto result = newton_euclidean(mesh, x0, maps); +// assert(result.converged); +// ════════════════════════════════════════════════════════════════════════════ + +TEST(Pipeline, Euclidean_TriangleToEquilibrium) +{ + // ── Step 1: build mesh ──────────────────────────────────────────────── + auto mesh = make_triangle(); + + // ── Step 2: set up maps ─────────────────────────────────────────────── + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + + // ── Step 3: assign DOFs (pin v0) ────────────────────────────────────── + int n = 0; + pin_first_vertex_euclidean(mesh, maps, n); + ASSERT_EQ(n, 2); + + // ── Step 4: choose natural target angles → x* = 0 ──────────────────── + set_natural_euclidean_theta(mesh, maps, n); + + // ── Step 5: solve ───────────────────────────────────────────────────── + std::vector x0(static_cast(n), -0.1); + auto result = newton_euclidean(mesh, x0, maps); + + // ── Step 6: verify ──────────────────────────────────────────────────── + EXPECT_TRUE(result.converged) + << "Euclidean pipeline: triangle should converge; " + "grad_inf_norm = " << result.grad_inf_norm; + EXPECT_LT(result.grad_inf_norm, 1e-8); + EXPECT_EQ(static_cast(result.x.size()), n); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Test 2 — Spherical full pipeline: tetrahedron → equilibrium +// +// The spherical tetrahedron equilibrium x* = 0 is built into the maps. +// After solving, the total angle defect Σ(Θ_v − Σα_v) should be ≈ 0. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(Pipeline, Spherical_TetrahedronToEquilibrium) +{ + // ── Steps 1–3 ───────────────────────────────────────────────────────── + 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); + + // ── Step 4: solve ───────────────────────────────────────────────────── + std::vector x0(static_cast(n), -0.2); + auto result = newton_spherical(mesh, x0, maps); + + // ── Step 5: verify convergence ──────────────────────────────────────── + EXPECT_TRUE(result.converged) + << "Spherical pipeline: tetrahedron should converge; " + "grad_inf_norm = " << result.grad_inf_norm; + EXPECT_LT(result.grad_inf_norm, 1e-8); + + // ── Step 6: verify geometric invariant — total angle defect ≈ 0 ────── + auto G_final = spherical_gradient(mesh, result.x, maps); + double total_defect = 0.0; + for (double gv : G_final) total_defect += gv; + EXPECT_NEAR(total_defect, 0.0, 1e-7) + << "Spherical: total angle defect should vanish at equilibrium"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Test 3 — HyperIdeal full pipeline: triangle → equilibrium +// ════════════════════════════════════════════════════════════════════════════ + +TEST(Pipeline, HyperIdeal_TriangleRoundTrip) +{ + // ── Steps 1–3 ───────────────────────────────────────────────────────── + auto mesh = make_triangle(); + auto maps = setup_hyper_ideal_maps(mesh); + int n = assign_all_dof_indices(mesh, maps); + + // ── Step 4: natural targets (equilibrium at b=1.0, a=0.5) ──────────── + auto xbase = set_natural_hyper_ideal_targets(mesh, maps, n); + + // Verify: gradient at xbase must be ≈ 0 before solving + auto G_at_base = evaluate_hyper_ideal(mesh, xbase, maps, false).gradient; + double max_g = 0.0; + for (double v : G_at_base) max_g = std::max(max_g, std::abs(v)); + ASSERT_LT(max_g, 1e-10) << "Natural target setup: gradient at base should be ~0"; + + // ── Step 5: perturb and solve ───────────────────────────────────────── + std::vector x0 = xbase; + for (auto& v : x0) v += 0.3; + + auto result = newton_hyper_ideal(mesh, x0, maps); + + // ── Step 6: verify ──────────────────────────────────────────────────── + EXPECT_TRUE(result.converged) + << "HyperIdeal pipeline: triangle should converge; " + "grad_inf_norm = " << result.grad_inf_norm; + EXPECT_LT(result.grad_inf_norm, 1e-8); + + // Solution should be close to xbase (same equilibrium) + for (int i = 0; i < n; ++i) { + EXPECT_NEAR(result.x[static_cast(i)], + xbase[static_cast(i)], 1e-6) + << "DOF " << i << " should recover the equilibrium value"; + } +} + +// ════════════════════════════════════════════════════════════════════════════ +// Test 4 — Mesh I/O in the pipeline: solve → write → reload → check +// +// Demonstrates: compute a conformal factor on a mesh, write it to OFF, reload +// and check that the mesh topology is preserved. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(Pipeline, MeshIO_SolveAndExport) +{ + // ── Build and solve ─────────────────────────────────────────────────── + auto mesh = make_quad_strip(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + + int n = 0; + pin_first_vertex_euclidean(mesh, maps, n); + set_natural_euclidean_theta(mesh, maps, n); + + std::vector x0(static_cast(n), -0.1); + auto result = newton_euclidean(mesh, x0, maps); + ASSERT_TRUE(result.converged) << "Solver must converge before export test"; + + // ── Write mesh ──────────────────────────────────────────────────────── + const std::string tmp_path = "/tmp/conformallab_pipeline_test.off"; + ASSERT_NO_THROW(save_mesh(tmp_path, mesh)); + ASSERT_TRUE(std::filesystem::exists(tmp_path)); + + // ── Reload and verify topology ──────────────────────────────────────── + ConformalMesh mesh2; + ASSERT_NO_THROW(mesh2 = load_mesh(tmp_path)); + + EXPECT_EQ(mesh2.number_of_vertices(), mesh.number_of_vertices()) + << "Vertex count must survive OFF round-trip"; + EXPECT_EQ(mesh2.number_of_faces(), mesh.number_of_faces()) + << "Face count must survive OFF round-trip"; + + std::filesystem::remove(tmp_path); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Test 5 — All three geometries, same quad-strip topology +// +// Validates that the solver infrastructure works uniformly: the same mesh +// topology is solvable under Euclidean, Spherical, and HyperIdeal geometries. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(Pipeline, AllThreeGeometries_QuadStrip) +{ + // ── Euclidean ────────────────────────────────────────────────────────── + { + auto mesh = make_quad_strip(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + int n = 0; + pin_first_vertex_euclidean(mesh, maps, n); + set_natural_euclidean_theta(mesh, maps, n); + std::vector x0(static_cast(n), -0.1); + auto res = newton_euclidean(mesh, x0, maps); + EXPECT_TRUE(res.converged) << "Euclidean: quad strip should converge"; + } + + // ── HyperIdeal ──────────────────────────────────────────────────────── + { + auto mesh = make_quad_strip(); + auto maps = setup_hyper_ideal_maps(mesh); + int n = assign_all_dof_indices(mesh, maps); + auto xbase = set_natural_hyper_ideal_targets(mesh, maps, n); + std::vector x0 = xbase; + for (auto& v : x0) v += 0.1; + auto res = newton_hyper_ideal(mesh, x0, maps); + EXPECT_TRUE(res.converged) << "HyperIdeal: quad strip should converge"; + } + + // ── Spherical (tetrahedron: smallest closed mesh with all vertices free) ─ + { + 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 x0(static_cast(n), -0.1); + auto res = newton_spherical(mesh, x0, maps); + EXPECT_TRUE(res.converged) << "Spherical: tetrahedron should converge"; + } +}