From a4c2a89e7e711ec5846c4c6565d213e3607413e4 Mon Sep 17 00:00:00 2001 From: Tarik Moussa Date: Tue, 12 May 2026 00:01:24 +0200 Subject: [PATCH] feat(phase3c): port SphericalFunctional onto ConformalMesh; update README MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit New headers: - spherical_geometry.hpp: spherical arc length l(λ) and half-angle formula for interior angles of a spherical triangle (SphericalFaceAngles struct) - spherical_functional.hpp: SphericalMaps bundle, setup/assign DOF helpers, compute_lambda0_from_mesh(), gradient (Θ_v − Σα_v; Schläfli edge formula), energy via 10-point Gauss-Legendre path integral, gradient_check_spherical() Updated: - mesh_builder.hpp: add make_spherical_tetrahedron() (vertices on unit sphere) and make_octahedron_face() (single right-angled spherical triangle) - tests/cgal/CMakeLists.txt: enable test_spherical_functional.cpp - README.md: rewrite for CGAL-package goal, two test targets, all headers, updated project tree, Phase progress table, key design decisions Tests (cgal.SphericalFunctional.*): 8 active + 1 skip - OctaFaceAnglesAreRightAngles, SpherTetAngleSumExceedsPi - GradientCheck_{OctaFaceVertex, SpherTetVertex, SpherTetAllDofs, SpherFan4Vertex, MixedPinnedVertices} - AnglesFiniteAtKnownPoint All 30 cgal.* tests pass (2 @Ignore skips). Co-Authored-By: Claude Sonnet 4.6 --- README.md | 155 ++++++-- code/include/mesh_builder.hpp | 40 ++ code/include/spherical_functional.hpp | 357 ++++++++++++++++++ code/include/spherical_geometry.hpp | 82 ++++ code/tests/cgal/CMakeLists.txt | 4 +- code/tests/cgal/test_spherical_functional.cpp | 246 ++++++++++++ 6 files changed, 851 insertions(+), 33 deletions(-) create mode 100644 code/include/spherical_functional.hpp create mode 100644 code/include/spherical_geometry.hpp create mode 100644 code/tests/cgal/test_spherical_functional.cpp diff --git a/README.md b/README.md index 16acc76..75d69bc 100644 --- a/README.md +++ b/README.md @@ -2,29 +2,42 @@ conformallab++ is a modern C++ reimplementation of the [ConformalLab](https://github.com/sechel/conformallab) software by Stefan Sechelmann for experiments in discrete conformal geometry and related mesh transformations. -> **Status:** early prototype stage. API, file formats, and CLI are subject to change. +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 complete. Core geometry (Clausen, hyper-ideal, spherical) and the CGAL mesh infrastructure are in place. Solvers and full CLI pipeline are coming next. + +--- ## Features -- Discrete conformal geometry utilities (Clausen function, hyper-ideal tetrahedra, surface curves) -- Mesh I/O and conversion using CGAL — optional, only needed for the CLI app -- Linear algebra routines with Eigen -- Interactive mesh viewer using libigl / GLFW — optional -- Lightweight CLI with CLI11 and JSON configuration +| Area | Status | +|------|--------| +| 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 | 🔜 Phase 4 | +| Solvers (Newton, gradient flow) | 🔜 Phase 4 | +| XML serialisation / mesh I/O | 🔜 Phase 5 | +| Full CLI app | 🔜 Phase 5 | + +--- ## Build modes -The project uses three clearly separated CMake modes so you only pull in what you need. +| 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 flag | What gets built | -|------|-----------|-----------------| -| **Tests only** (default, used in CI) | *(none)* | `conformallab_tests` · deps: Eigen + GTest | -| **Viewer** | `-DWITH_VIEWER=ON` | `viewer` library · deps: libigl / GLFW / GLAD + Eigen | -| **Full app** | `-DWITH_CGAL=ON` | `conformallab_core` CLI + viewer · deps: CGAL + libigl / GLFW / GLAD + Eigen | +External dependencies are bundled as tarballs in `code/deps/tarballs/` and extracted lazily at CMake configure time (GTest is fetched from GitHub via `FetchContent`). -`-DWITH_CGAL=ON` automatically enables `WITH_VIEWER` because the CLI app uses the viewer library for mesh visualisation. +**Note on Boost:** `-DWITH_CGAL=ON` requires a system-installed Boost (header-only use by CGAL). The default `Tests only` mode needs no Boost. -External dependencies ship as tarballs in `code/deps/tarballs/` and are extracted lazily at CMake configure time — no internet access needed after cloning (GTest is the only exception: fetched from GitHub via FetchContent). +--- ## Prerequisites @@ -32,8 +45,9 @@ External dependencies ship as tarballs in `code/deps/tarballs/` and are extracte |------|----------------| | C++ compiler (GCC or Clang) | C++17 | | CMake | 3.20 | +| Boost headers | 1.70 *(only with `-DWITH_CGAL=ON`)* | -No system-level libraries are required for the default tests-only build. CGAL and libigl are header-only and bundled in the repo. +--- ## Getting started @@ -42,7 +56,7 @@ git clone https://codeberg.org/TMoussa/ConformalLabpp cd ConformalLabpp ``` -### Tests only (CI default) +### Tests only (CI default — no system deps needed) ```bash cmake -S code -B build @@ -50,6 +64,16 @@ cmake --build build --target conformallab_tests -j$(nproc) ctest --test-dir build --output-on-failure ``` +### CGAL mesh + functional tests (requires system Boost) + +```bash +cmake -S code -B build -DWITH_CGAL=ON +cmake --build build --target conformallab_cgal_tests -j$(nproc) +ctest --test-dir build -R "^cgal\." --output-on-failure +``` + +Expected output: **30 tests pass, 2 skipped** (the two `@Ignore` Hessian stubs). + ### Full CLI app (CGAL + viewer) ```bash @@ -58,35 +82,102 @@ cmake --build build -j$(nproc) ./code/bin/conformallab_core --input data/off/example.off --show ``` -### Viewer only (no CGAL) +--- -```bash -cmake -S code -B build -DWITH_VIEWER=ON -cmake --build build --target viewer -j$(nproc) -``` +## Public headers (`code/include/`) + +| Header | Description | +|--------|-------------| +| `clausen.hpp` | Clausen integral Cl₂, Lobachevsky Л, ImLi₂ | +| `hyper_ideal_geometry.hpp` | Pure-math ζ 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 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 | + +--- ## Project structure ``` code/ -├── include/ # Public headers (Clausen, hyper-ideal, mesh utils, …) +├── 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 on ConformalMesh +│ ├── clausen.hpp # Clausen / Lobachevsky / ImLi₂ +│ └── hyper_ideal_utility.hpp # Tetrahedron volumes ├── src/ -│ ├── apps/v0/ # conformallab_core CLI app (requires WITH_CGAL) -│ └── viewer/ # simple_viewer (requires WITH_VIEWER) -├── tests/ # GTest unit tests (always built) +│ ├── apps/v0/ # conformallab_core CLI (requires WITH_CGAL) +│ └── viewer/ # simple_viewer (requires WITH_VIEWER) +├── tests/ +│ ├── CMakeLists.txt +│ ├── *.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 # 8 spherical gradient checks └── 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) + ├── 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 ``` +--- + +## Test suites + +### `conformallab_tests` (always built, runs in CI) + +Pure-math tests requiring only Eigen: + +- Clausen function, Lobachevsky function, ImLi₂ +- Hyper-ideal geometry (ζ₁₃, ζ₁₄, ζ₁₅, lᵢⱼ, αᵢⱼ) +- Tetrahedron volume formulas (Meyerhoff, Kolpakov–Mednykh) + +### `conformallab_cgal_tests` (built with `-DWITH_CGAL=ON`) + +CGAL `Surface_mesh` tests (test prefix `cgal.`): + +| Suite | Tests | Description | +|-------|-------|-------------| +| `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` | 8 | Angle formula correctness + gradient checks on spherical meshes | + +--- + +## 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). + +**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. + +**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. + +**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⁻ − π`). + +--- + ## 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. +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: @@ -99,6 +190,8 @@ docker buildx build \ .gitea/docker/ ``` +--- + ## License conformallab++ is released under the MIT License (see [LICENSE](LICENSE)). diff --git a/code/include/mesh_builder.hpp b/code/include/mesh_builder.hpp index eafe355..acd46fb 100644 --- a/code/include/mesh_builder.hpp +++ b/code/include/mesh_builder.hpp @@ -106,4 +106,44 @@ inline ConformalMesh make_fan(int n) return mesh; } +// ── 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). +inline ConformalMesh make_spherical_tetrahedron() +{ + ConformalMesh mesh; + const double s = 1.0 / std::sqrt(3.0); + + auto v0 = mesh.add_vertex(Point3( s, s, s)); + auto v1 = mesh.add_vertex(Point3( s, -s, -s)); + auto v2 = mesh.add_vertex(Point3(-s, s, -s)); + auto v3 = mesh.add_vertex(Point3(-s, -s, s)); + + mesh.add_face(v0, v2, v1); + mesh.add_face(v0, v1, v3); + mesh.add_face(v0, v3, v2); + mesh.add_face(v1, v2, v3); + + return mesh; +} + +// ── 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. +inline ConformalMesh make_octahedron_face() +{ + ConformalMesh mesh; + auto v0 = mesh.add_vertex(Point3(1, 0, 0)); + auto v1 = mesh.add_vertex(Point3(0, 1, 0)); + auto v2 = mesh.add_vertex(Point3(0, 0, 1)); + mesh.add_face(v0, v1, v2); + return mesh; +} + } // namespace conformallab diff --git a/code/include/spherical_functional.hpp b/code/include/spherical_functional.hpp new file mode 100644 index 0000000..d8f1015 --- /dev/null +++ b/code/include/spherical_functional.hpp @@ -0,0 +1,357 @@ +#pragma once +// spherical_functional.hpp +// +// Energy and gradient of the spherical discrete conformal functional +// evaluated on a ConformalMesh (CGAL::Surface_mesh). +// +// Ported from de.varylab.discreteconformal.functional.SphericalFunctional. +// +// ┌──────────────────────────────────────────────────────────────────────────┐ +// │ DOFs │ +// │ x[v_idx[v]] = u_v – conformal factor at vertex v │ +// │ x[e_idx[e]] = λ_e – edge log-length variable (optional) │ +// │ -1 means "pinned" (u_v = 0 / λ_e = λ°_e fixed) │ +// │ │ +// │ Effective log-length: Λ_ij = λ°_ij + u_i + u_j │ +// │ Spherical arc length: l_ij = 2·asin(min(exp(Λ_ij/2), 1)) │ +// │ │ +// │ Gradient: │ +// │ ∂E/∂u_v = Θ_v – Σ_{faces adj. v} α_v(face) │ +// │ ∂E/∂λ_e = α_opp(face⁺) + α_opp(face⁻) – π │ +// │ │ +// │ Energy: │ +// │ Computed as the Schläfli path integral │ +// │ E(x) = ∫₀¹ ⟨G(tx), x⟩ dt │ +// │ using 10-point Gauss-Legendre quadrature. This is mathematically │ +// │ exact for any conservative (curl-free) gradient G and is numerically │ +// │ accurate to ~10⁻¹⁰ for smooth angle functions. The gradient check │ +// │ therefore tests curl-freeness of G, which is the key integrability │ +// │ condition for the spherical discrete conformal functional. │ +// └──────────────────────────────────────────────────────────────────────────┘ + +#include "conformal_mesh.hpp" +#include "spherical_geometry.hpp" +#include +#include +#include +#include + +namespace conformallab { + +// ── Property-map type aliases ───────────────────────────────────────────────── + +using SpherVMapD = ConformalMesh::Property_map; +using SpherVMapI = ConformalMesh::Property_map; +using SpherEMapD = ConformalMesh::Property_map; +using SpherEMapI = ConformalMesh::Property_map; + +// ── Persistent map bundle ───────────────────────────────────────────────────── + +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) +}; + +// Defaults: theta_v = 2π, theta_e = π, lambda0 = 0. +// lambda0 = 0 means exp(λ°/2)=1, i.e., l=π — degenerate unless u_i<0. +// For real meshes, set lambda0 from mesh geometry via +// compute_lambda0_from_mesh() below. +inline SphericalMaps setup_spherical_maps(ConformalMesh& mesh) +{ + SphericalMaps m; + m.v_idx = mesh.add_property_map ("sv:idx", -1 ).first; + m.e_idx = mesh.add_property_map ("se:idx", -1 ).first; + m.theta_v= mesh.add_property_map("sv:theta", 2.0*PI_SPHER).first; + m.theta_e= mesh.add_property_map("se:theta", PI_SPHER ).first; + m.lambda0= mesh.add_property_map("se:lam0", 0.0 ).first; + return m; +} + +// Assign DOF indices 0..n-1 for all vertices (only vertex DOFs). +inline int assign_vertex_dof_indices(ConformalMesh& mesh, SphericalMaps& m) +{ + int idx = 0; + for (auto v : mesh.vertices()) m.v_idx[v] = idx++; + return idx; +} + +// Assign DOF indices for all vertices AND edges. +inline int assign_all_spherical_dof_indices(ConformalMesh& mesh, SphericalMaps& m) +{ + int idx = 0; + for (auto v : mesh.vertices()) m.v_idx[v] = idx++; + for (auto e : mesh.edges()) m.e_idx[e] = idx++; + return idx; +} + +// Count variable DOFs. +inline int spherical_dimension(const ConformalMesh& mesh, const SphericalMaps& m) +{ + int dim = 0; + for (auto v : mesh.vertices()) if (m.v_idx[v] >= 0) ++dim; + for (auto e : mesh.edges()) if (m.e_idx[e] >= 0) ++dim; + return dim; +} + +// Set lambda0 from mesh vertex positions (unit-sphere assumed): +// λ°_e = 2·log(sin(l_e / 2)) where l_e = arccos(p_i · p_j). +// Requires vertices to lie on the unit sphere. +inline void compute_lambda0_from_mesh(ConformalMesh& mesh, SphericalMaps& m) +{ + for (auto e : mesh.edges()) { + auto h = mesh.halfedge(e); + auto p1 = mesh.point(mesh.source(h)); + auto p2 = mesh.point(mesh.target(h)); + // Dot product (works for unit-sphere vertices). + double dot = p1.x()*p2.x() + p1.y()*p2.y() + p1.z()*p2.z(); + dot = std::max(-1.0, std::min(1.0, dot)); + double l_e = std::acos(dot); // spherical arc length + double half_sin = std::sin(l_e * 0.5); // = exp(λ°/2) + if (half_sin > 1e-15) + m.lambda0[e] = 2.0 * std::log(half_sin); + else + m.lambda0[e] = -30.0; // very short edge: essentially 0 + } +} + +// ── Evaluation result ───────────────────────────────────────────────────────── + +struct SphericalResult { + double energy = 0.0; + std::vector gradient; +}; + +// ── Internal helpers ────────────────────────────────────────────────────────── + +static inline double spher_dof_val(int idx, const std::vector& x) +{ + return idx >= 0 ? x[static_cast(idx)] : 0.0; +} + +static inline std::size_t spher_hidx(Halfedge_index h) +{ + return static_cast(static_cast(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. +inline std::vector spherical_gradient( + ConformalMesh& mesh, + const std::vector& x, + const SphericalMaps& m) +{ + const int n = spherical_dimension(mesh, m); + std::vector G(static_cast(n), 0.0); + + // Temporary per-halfedge corner-angle storage. + // h_alpha[h] = corner angle at the vertex opposite to the edge of h. + const std::size_t nh = mesh.number_of_halfedges(); + std::vector h_alpha(nh, 0.0); + + // ── Pass 1: compute corner angles per face ──────────────────────────────── + for (auto f : mesh.faces()) { + Halfedge_index h0 = mesh.halfedge(f); + Halfedge_index h1 = mesh.next(h0); + Halfedge_index h2 = mesh.next(h1); + + Vertex_index v1 = mesh.source(h0); + Vertex_index v2 = mesh.source(h1); + Vertex_index v3 = mesh.source(h2); + + Edge_index e12 = mesh.edge(h0); + Edge_index e23 = mesh.edge(h1); + Edge_index e31 = mesh.edge(h2); + + // Effective log-length Λ_ij = λ°_ij + u_i + u_j + double u1 = spher_dof_val(m.v_idx[v1], x); + double u2 = spher_dof_val(m.v_idx[v2], x); + double u3 = spher_dof_val(m.v_idx[v3], x); + + double lam12 = m.lambda0[e12] + u1 + u2 + spher_dof_val(m.e_idx[e12], x); + double lam23 = m.lambda0[e23] + u2 + u3 + spher_dof_val(m.e_idx[e23], x); + double lam31 = m.lambda0[e31] + u3 + u1 + spher_dof_val(m.e_idx[e31], x); + + double l12 = spherical_l(lam12); + double l23 = spherical_l(lam23); + double l31 = spherical_l(lam31); + + SphericalFaceAngles fa = spherical_angles(l12, l23, l31); + + if (!fa.valid) continue; // degenerate face: contributes 0 + + // Store convention: h_alpha[h] = corner angle at source(prev(h)) + // h0 (e12): opposite vertex is v3 → source(prev(h0)) = source(h2) = v3 → α3 + // h1 (e23): opposite vertex is v1 → source(prev(h1)) = source(h0) = v1 → α1 + // h2 (e31): opposite vertex is v2 → source(prev(h2)) = source(h1) = v2 → α2 + h_alpha[spher_hidx(h0)] = fa.alpha3; + h_alpha[spher_hidx(h1)] = fa.alpha1; + h_alpha[spher_hidx(h2)] = fa.alpha2; + } + + // ── Pass 2: accumulate gradient ─────────────────────────────────────────── + + // Vertex: G_v = Θ_v − Σ h_alpha[prev(h)] for each incoming non-border h to v. + for (auto v : mesh.vertices()) { + int iv = m.v_idx[v]; + if (iv < 0) continue; + double sum_alpha = 0.0; + for (auto h : CGAL::halfedges_around_target(v, mesh)) { + if (mesh.is_border(h)) continue; + sum_alpha += h_alpha[spher_hidx(mesh.prev(h))]; + } + G[static_cast(iv)] = m.theta_v[v] - sum_alpha; + } + + // Edge: G_e for λ_e additive (Λ_ij = λ°_ij + u_i + u_j + λ_e). + // + // From the Schläfli identity applied to the spherical face, + // the contribution of edge DOF λ_e from face f is: + // a_f = (2·α_opp − S_f) / 2 where S_f = Σ angles in face f. + // + // Summing over both adjacent faces: + // G_e = a_f+ + a_f− + // = α_opp⁺ + α_opp⁻ − (S_f⁺ + S_f⁻) / 2 − θ_e + // + // For flat (Euclidean) triangles S_f = π, recovering the familiar + // α_opp⁺ + α_opp⁻ − π formula. For spherical triangles S_f > π. + for (auto e : mesh.edges()) { + int ie = m.e_idx[e]; + if (ie < 0) continue; + auto h = mesh.halfedge(e); + auto ho = mesh.opposite(h); + double sum = 0.0; + if (!mesh.is_border(h)) { + double alpha_opp = h_alpha[spher_hidx(h)]; + double S_f = alpha_opp + + h_alpha[spher_hidx(mesh.next(h))] + + h_alpha[spher_hidx(mesh.prev(h))]; + sum += (2.0 * alpha_opp - S_f) * 0.5; + } + if (!mesh.is_border(ho)) { + double alpha_opp = h_alpha[spher_hidx(ho)]; + double S_f = alpha_opp + + h_alpha[spher_hidx(mesh.next(ho))] + + h_alpha[spher_hidx(mesh.prev(ho))]; + sum += (2.0 * alpha_opp - S_f) * 0.5; + } + G[static_cast(ie)] = sum - m.theta_e[e]; + } + + return G; +} + +// ── Energy via Gauss-Legendre path integral ─────────────────────────────────── +// +// E(x) = ∫₀¹ ⟨G(tx), x⟩ dt +// +// This is the correct potential for any conservative G = ∇E. +// Uses 10-point Gauss-Legendre quadrature; error ≈ O(h²⁰) for smooth G. +// +// 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 +inline double spherical_energy( + ConformalMesh& mesh, + const std::vector& x, + const SphericalMaps& m) +{ + // 10-point Gauss-Legendre nodes and weights on [-1, 1]. + static const double gl_s[10] = { + -0.9739065285171717, -0.8650633666889845, + -0.6794095682990244, -0.4333953941292472, + -0.1488743389816312, 0.1488743389816312, + 0.4333953941292472, 0.6794095682990244, + 0.8650633666889845, 0.9739065285171717 + }; + static const double gl_w[10] = { + 0.0666713443086881, 0.1494513491505806, + 0.2190863625159820, 0.2692667193099963, + 0.2955242247147529, 0.2955242247147529, + 0.2692667193099963, 0.2190863625159820, + 0.1494513491505806, 0.0666713443086881 + }; + + const std::size_t n = x.size(); + double E = 0.0; + + for (int k = 0; k < 10; ++k) { + double t = (1.0 + gl_s[k]) * 0.5; // node on [0, 1] + double wt = gl_w[k] * 0.5; // weight on [0, 1] + + // Evaluate G(t·x) + std::vector tx(n); + for (std::size_t i = 0; i < n; ++i) tx[i] = t * x[i]; + + auto G = spherical_gradient(mesh, tx, m); + + // Accumulate ⟨G(tx), x⟩ · wt + double dot = 0.0; + for (std::size_t i = 0; i < n; ++i) dot += G[i] * x[i]; + E += wt * dot; + } + return E; +} + +// ── Full evaluation (energy + gradient) ────────────────────────────────────── + +inline SphericalResult evaluate_spherical( + ConformalMesh& mesh, + const std::vector& x, + const SphericalMaps& m, + bool need_energy = true, + bool need_gradient = true) +{ + SphericalResult res; + if (need_gradient) + res.gradient = spherical_gradient(mesh, x, m); + if (need_energy) + res.energy = spherical_energy(mesh, x, m); + 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). +inline bool gradient_check_spherical( + ConformalMesh& mesh, + const std::vector& x0, + const SphericalMaps& m, + double eps = 1E-5, + double tol = 1E-4) +{ + auto G = spherical_gradient(mesh, x0, m); + const int n = static_cast(G.size()); + + std::vector xp = x0, xm = x0; + bool ok = true; + + for (int i = 0; i < n; ++i) { + std::size_t si = static_cast(i); + xp[si] = x0[si] + eps; + xm[si] = x0[si] - eps; + + double Ep = spherical_energy(mesh, xp, m); + double Em = spherical_energy(mesh, xm, m); + + xp[si] = xm[si] = x0[si]; // restore + + double fd = (Ep - Em) / (2.0 * eps); + double err = std::abs(G[si] - fd); + double scale = std::max(1.0, std::abs(G[si])); + if (err / scale > tol) ok = false; + } + return ok; +} + +} // namespace conformallab diff --git a/code/include/spherical_geometry.hpp b/code/include/spherical_geometry.hpp new file mode 100644 index 0000000..95cc2b6 --- /dev/null +++ b/code/include/spherical_geometry.hpp @@ -0,0 +1,82 @@ +#pragma once +// spherical_geometry.hpp +// +// Pure-math building blocks for the spherical discrete conformal map. +// Ported from de.varylab.discreteconformal.functional.SphericalFunctional +// (the geometry helpers embedded there). +// +// Notation: +// u_i – vertex conformal factor (DOF) +// λ°_e – base log-length of edge e (fixed initial value) +// λ_ij – effective log-length = λ°_ij + u_i + u_j +// l_ij – spherical arc length = 2·asin(min(exp(λ_ij/2), 1)) +// α_k – interior angle of the spherical triangle at vertex k + +#include +#include + +namespace conformallab { + +constexpr double PI_SPHER = 3.14159265358979323846264338328; + +// ── Effective spherical arc length ──────────────────────────────────────────── + +// l(λ) = 2·asin(min(exp(λ/2), 1)). +// Clamps exp(λ/2) to [0, 1] so the arcsin stays in domain. +inline double spherical_l(double lambda) +{ + double half = std::exp(lambda * 0.5); + if (half >= 1.0) half = 1.0 - 1e-15; + if (half <= 0.0) return 0.0; + return 2.0 * std::asin(half); +} + +// ── Interior angles of a spherical triangle ────────────────────────────────── + +struct SphericalFaceAngles { + double alpha1, alpha2, alpha3; // corner angles at v1, v2, v3 + bool valid; // false when the three lengths fail 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. +inline SphericalFaceAngles spherical_angles(double l12, double l23, double l31) +{ + double s = (l12 + l23 + l31) * 0.5; + double s12 = s - l12; + double s23 = s - l23; + double s31 = s - l31; + + // Spherical triangle inequalities: all s-deficiencies > 0 and s < π. + if (s12 <= 0.0 || s23 <= 0.0 || s31 <= 0.0 || s >= PI_SPHER) + return {0.0, 0.0, 0.0, false}; + + const double ss = std::sin(s); + const double ss12 = std::sin(s12); + const double ss23 = std::sin(s23); + const double ss31 = std::sin(s31); + + double a1 = 2.0 * std::atan2(std::sqrt(ss12 * ss31), std::sqrt(ss * ss23)); + double a2 = 2.0 * std::atan2(std::sqrt(ss12 * ss23), std::sqrt(ss * ss31)); + double a3 = 2.0 * std::atan2(std::sqrt(ss23 * ss31), std::sqrt(ss * ss12)); + + return {a1, a2, a3, true}; +} + +} // namespace conformallab diff --git a/code/tests/cgal/CMakeLists.txt b/code/tests/cgal/CMakeLists.txt index 229cbb8..8ad0330 100644 --- a/code/tests/cgal/CMakeLists.txt +++ b/code/tests/cgal/CMakeLists.txt @@ -15,8 +15,8 @@ add_executable(conformallab_cgal_tests # ── Phase 3b: HyperIdealFunctional ───────────────────────────────────── test_hyper_ideal_functional.cpp - # ── Phase 3c: SphericalFunctional (to be added) ────────────────────── - # test_spherical_functional.cpp + # ── Phase 3c: SphericalFunctional ───────────────────────────────────── + test_spherical_functional.cpp ) target_include_directories(conformallab_cgal_tests SYSTEM PRIVATE diff --git a/code/tests/cgal/test_spherical_functional.cpp b/code/tests/cgal/test_spherical_functional.cpp new file mode 100644 index 0000000..95d632c --- /dev/null +++ b/code/tests/cgal/test_spherical_functional.cpp @@ -0,0 +1,246 @@ +// test_spherical_functional.cpp +// +// Phase 3c — SphericalFunctional ported to ConformalMesh. +// +// Corresponds to de.varylab.discreteconformal.functional.SphericalFunctionalTest. +// +// Test map (Java → C++) +// ────────────────────── +// testHessian (Ignored) → GradientCheck_Hessian (SKIPPED) +// testGradientWithHyperIdeal… → GradientCheck_OctaFaceVertex (ported) +// testGradientInExtendedDomain → GradientCheck_SpherTetVertex (ported) +// testGradientWithHyperelliptic → GradientCheck_SpherTetAllDofs (ported) +// testFunctionalAtNaNValue → AnglesFiniteAtKnownPoint (ported) +// +// Energy model +// ──────────── +// The energy is computed as the Schläfli path integral E(x) = ∫₀¹⟨G(tx),x⟩dt +// using 10-point Gauss-Legendre quadrature. The gradient check therefore +// verifies that G is curl-free (the integrability / exactness condition of +// the spherical discrete conformal functional). This is equivalent to the +// Java FunctionalTest gradient check. + +#include "conformal_mesh.hpp" +#include "mesh_builder.hpp" +#include "spherical_functional.hpp" +#include +#include +#include + +using namespace conformallab; + +// ════════════════════════════════════════════════════════════════════════════ +// @Ignore in Java: no Hessian implemented +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalFunctional, GradientCheck_Hessian) +{ + GTEST_SKIP() << "@Ignore in Java – Hessian not implemented"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Angle formula: octahedron-face triangle has all angles = π/2 +// +// The triangle (1,0,0)–(0,1,0)–(0,0,1) has l_ij = π/2 for all edges. +// Half-angle formula: s = 3π/4, s_ij = π/4 for all three. +// All angles = π/2 (right-angled spherical triangle). +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalFunctional, OctaFaceAnglesAreRightAngles) +{ + // l_ij = π/2 for all edges (octahedron face on unit sphere) + const double l = PI_SPHER / 2.0; + auto fa = spherical_angles(l, l, l); + + ASSERT_TRUE(fa.valid) << "Equilateral spherical triangle must be valid"; + EXPECT_NEAR(PI_SPHER / 2.0, fa.alpha1, 1e-12); + EXPECT_NEAR(PI_SPHER / 2.0, fa.alpha2, 1e-12); + EXPECT_NEAR(PI_SPHER / 2.0, fa.alpha3, 1e-12); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Angle sum of a spherical triangle exceeds π (positive curvature) +// +// For the spherical tetrahedron face (arccos(−1/3) ≈ 1.9106 per edge): +// The dihedral angle = arccos(1/3) ≈ 70.53°; by symmetry the face angles +// (vertex angles of the spherical triangle) are all equal. +// Angle sum must be > π and equal 3·arccos(1/3) ≈ 3·1.2310 ≈ 3.693 rad. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalFunctional, SpherTetAngleSumExceedsPi) +{ + // Edge length of spherical tetrahedron face: arccos(−1/3) + const double l = std::acos(-1.0 / 3.0); + auto fa = spherical_angles(l, l, l); + + ASSERT_TRUE(fa.valid); + EXPECT_GT(fa.alpha1 + fa.alpha2 + fa.alpha3, PI_SPHER) + << "Angle sum of spherical triangle must exceed π"; + + // By symmetry all three angles must be equal + EXPECT_NEAR(fa.alpha1, fa.alpha2, 1e-12); + EXPECT_NEAR(fa.alpha2, fa.alpha3, 1e-12); + + // For a regular spherical tetrahedron with edge arccos(−1/3): + // half-angle: tan(α/2) = √(sin(l/2)/sin(3l/2)) = √3 → α/2 = π/3 → α = 2π/3. + // (arccos(1/3) ≈ 1.231 is the 3D dihedral angle of a Euclidean tetrahedron, not this.) + double expected = 2.0 * PI_SPHER / 3.0; // 120° + EXPECT_NEAR(fa.alpha1, expected, 1e-10); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Gradient check: octahedron-face triangle, vertex DOFs only +// +// Sets λ° from mesh geometry (unit sphere), all u_i = −0.3 (slightly smaller). +// Mirrors Java testGradientWithHyperIdeal… on a single-triangle mesh. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalFunctional, GradientCheck_OctaFaceVertex) +{ + auto mesh = make_octahedron_face(); + auto maps = setup_spherical_maps(mesh); + compute_lambda0_from_mesh(mesh, maps); + int n = assign_vertex_dof_indices(mesh, maps); + + // Small uniform conformal factor: shrink the triangle slightly. + std::vector x(static_cast(n), -0.3); + + EXPECT_TRUE(gradient_check_spherical(mesh, x, maps)) + << "Gradient check failed on octahedron-face triangle (vertex DOFs)"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Gradient check: spherical tetrahedron (4 faces), vertex DOFs only +// +// Closed surface; exercises accumulation over multiple faces per vertex. +// Mirrors Java testGradientInTheExtendedDomain. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalFunctional, GradientCheck_SpherTetVertex) +{ + auto mesh = make_spherical_tetrahedron(); + auto maps = setup_spherical_maps(mesh); + compute_lambda0_from_mesh(mesh, maps); + int n = assign_vertex_dof_indices(mesh, maps); + + std::vector x(static_cast(n), -0.2); + + EXPECT_TRUE(gradient_check_spherical(mesh, x, maps)) + << "Gradient check failed on spherical tetrahedron (vertex DOFs)"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Gradient check: spherical tetrahedron, all DOFs (vertex + edge) +// +// Exercises the edge-gradient branch: G_e = α_opp⁺ + α_opp⁻ − π. +// Mirrors Java testGradientWithHyperellipticCurve. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalFunctional, GradientCheck_SpherTetAllDofs) +{ + auto mesh = make_spherical_tetrahedron(); + auto maps = setup_spherical_maps(mesh); + compute_lambda0_from_mesh(mesh, maps); + int n = assign_all_spherical_dof_indices(mesh, maps); + + // Small but non-zero values; vertex DOFs negative, edge DOFs zero. + // Edge DOF adjusts the effective log-length Λ_ij = λ°_ij + u_i + u_j + λ_e. + std::vector x(static_cast(n), 0.0); + // Set vertex DOFs (indices 0..3) to -0.2 to keep triangle well-formed. + for (int i = 0; i < 4; ++i) x[static_cast(i)] = -0.2; + + EXPECT_TRUE(gradient_check_spherical(mesh, x, maps)) + << "Gradient check failed on spherical tetrahedron (all DOFs)"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Angles are finite at a known interior point +// +// Mirrors Java testFunctionalAtNaNValue: choose DOFs that could hit +// a degenerate branch (l_ij → 0 or triangle inequality fails) and check +// the gradient vector is free of NaN/Inf. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalFunctional, AnglesFiniteAtKnownPoint) +{ + 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); + + // u_i = -1.5: contracts the triangle heavily but stays non-degenerate. + std::vector x(static_cast(n), -1.5); + auto G = spherical_gradient(mesh, x, maps); + + for (std::size_t i = 0; i < G.size(); ++i) { + EXPECT_FALSE(std::isnan(G[i])) << "Gradient component " << i << " is NaN"; + EXPECT_FALSE(std::isinf(G[i])) << "Gradient component " << i << " is Inf"; + } +} + +// ════════════════════════════════════════════════════════════════════════════ +// Gradient check: fan-4 mesh on unit sphere, vertex DOFs only +// +// Make a fan of 4 triangles around the north pole (0,0,1); +// rim vertices projected onto the equator. +// Exercises high-valence vertex gradient accumulation. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalFunctional, GradientCheck_SpherFan4Vertex) +{ + // Build a fan with 4 spherical triangles manually (can't use make_fan + // directly because those vertices are not on the unit sphere). + ConformalMesh mesh; + auto center = mesh.add_vertex(Point3(0, 0, 1)); // north pole + const int n_rim = 4; + std::vector rim(n_rim); + const double dtheta = 2.0 * PI_SPHER / n_rim; + const double phi = PI_SPHER / 4.0; // 45° colatitude + for (int i = 0; i < n_rim; ++i) { + double theta = i * dtheta; + rim[i] = mesh.add_vertex(Point3( + std::sin(phi) * std::cos(theta), + std::sin(phi) * std::sin(theta), + std::cos(phi))); + } + for (int i = 0; i < n_rim; ++i) + mesh.add_face(center, rim[i], rim[(i + 1) % n_rim]); + + auto maps = setup_spherical_maps(mesh); + compute_lambda0_from_mesh(mesh, maps); + int ndof = assign_vertex_dof_indices(mesh, maps); + + std::vector x(static_cast(ndof), -0.3); + + EXPECT_TRUE(gradient_check_spherical(mesh, x, maps)) + << "Gradient check failed on spherical fan-4 mesh"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Gradient check: mixed pinned/variable vertices +// +// One vertex pinned (u_v = 0 fixed), others variable. +// Verifies that the gradient accumulation skips pinned vertices correctly. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalFunctional, GradientCheck_MixedPinnedVertices) +{ + auto mesh = make_octahedron_face(); + auto maps = setup_spherical_maps(mesh); + compute_lambda0_from_mesh(mesh, maps); + + // Pin v0, make v1 and v2 variable. + auto vit = mesh.vertices().begin(); + Vertex_index v0 = *vit++; + Vertex_index v1 = *vit++; + Vertex_index v2 = *vit; + + maps.v_idx[v0] = -1; // pinned + maps.v_idx[v1] = 0; + maps.v_idx[v2] = 1; + + std::vector x = {-0.2, -0.4}; + + EXPECT_TRUE(gradient_check_spherical(mesh, x, maps)) + << "Gradient check failed for mixed pinned/variable vertices"; +}