diff --git a/.gitea/workflows/cpp-tests.yml b/.gitea/workflows/cpp-tests.yml index 80e7852..4bc4977 100644 --- a/.gitea/workflows/cpp-tests.yml +++ b/.gitea/workflows/cpp-tests.yml @@ -1,5 +1,16 @@ name: C++ Tests +# Trigger keywords in commit message (checked via head_commit.message): +# /test-cgal — full CGAL test suite (277 tests, ~5 min build) +# /quality-gates — license, codespell, shellcheck, CGAL conventions +# /ci-all — all of the above + /docs + /links (across all workflows) +# +# Examples: +# git commit -m "fix: correct angle formula /test-cgal" +# git commit -m "release prep /ci-all" +# +# test-fast always runs on every push — it is fast (< 5 s) and cheap. + on: push: branches: @@ -7,12 +18,13 @@ on: - dev - "claude/**" - "feature/**" + - "review/**" pull_request: # ───────────────────────────────────────────────────────────────────────────── # Job 1 — test-fast # Pure-math tests (Clausen, ImLi₂, Hyper-ideal geometry). -# No CGAL, no Boost. Eigen + GTest only. Runs on ALL branches. +# No CGAL, no Boost. Eigen + GTest only. Runs on EVERY push. # ───────────────────────────────────────────────────────────────────────────── jobs: test-fast: @@ -48,54 +60,41 @@ jobs: # ───────────────────────────────────────────────────────────────────────────── # Job 2 — test-cgal -# Full CGAL test suite (Phase 3–7, 158 tests). -# Runs ONLY on pull requests (not on direct pushes to dev/main). -# Starts only after test-fast succeeds. # -# Uses -DWITH_CGAL_TESTS=ON (not -DWITH_CGAL=ON) to avoid building -# Viewer/GLFW — the CI container has no wayland-scanner. +# Trigger: include "/test-cgal" anywhere in the commit message. # -# Boost (libboost-dev) is already present in the container since the image rebuild. +# git commit -m "fix: correct angle formula /test-cgal" +# +# Why keyword-triggered (not automatic on every push): +# The Pi runner (3-4 GB RAM, swap heavily loaded) cannot sustain a +# CGAL build on every WIP commit. Adding the keyword to a commit +# message explicitly signals "this commit is ready for full testing". +# +# LOW_MEMORY_BUILD applies four RAM-saving measures so the build fits in +# a 2000 MB container: -O0, no PCH, unity batch 1, --no-keep-memory. +# See code/tests/cgal/CMakeLists.txt for the full explanation. # ───────────────────────────────────────────────────────────────────────────── - # ─── DISABLED 2026-05-26 ────────────────────────────────────────────────── - # The full CGAL test build (`conformallab_cgal_tests`, -j1, ~minutes) is - # too compute-intensive for the eulernest runner: even with the 1600 MB - # memory bump the job OOMs intermittently during CGAL+Eigen template - # expansion, so most recent runs fail without exposing a real regression. - # While we work through this, the job is gated off via `if: false` — - # `workflow_dispatch` reruns from the Gitea UI still work, and the body - # of the job is preserved unchanged for easy reactivation. - # - # Consequences: - # * test-fast (Job 1) still runs on every push/PR — pure-math tests - # stay gated. - # * quality-gates (Job 3) still runs on every push/PR — style / - # convention checks stay gated. - # * The two structural sub-gates nested under test-cgal - # (`scripts/check-test-counts.sh`, `scripts/try_it.sh`) are - # temporarily un-gated. Run them locally before tagging a release; - # a follow-up will relocate them into a cheaper job. - # - # To re-enable: change `if: false` back to - # `if: github.event_name == 'pull_request'`. test-cgal: needs: test-fast - if: false # DISABLED 2026-05-26 — see comment block above + if: | + contains(github.event.head_commit.message, '/test-cgal') || + contains(github.event.head_commit.message, '/ci-all') runs-on: eulernest container: image: git.eulernest.eu/conformallab/ci-cpp:latest - # Memory bumped from 1400m → 1600m to avoid OOM during CGAL header - # compilation on ARM64 (CGAL + Eigen templates allocate ~700 MB per - # cc1plus instance; -j1 leaves a small margin). - # memory-swap == memory disables swap entirely so OOM fails fast - # rather than thrashing on the SD card. - options: "--memory=1600m --memory-swap=1600m" + # 2000 MB hard limit; 1000 MB swap headroom (memory-swap = RAM + swap). + # With LOW_MEMORY_BUILD peak per TU is ~150-200 MB → well within limit. + options: "--memory=2000m --memory-swap=3000m" steps: - uses: actions/checkout@v4 - - name: Configure (WITH_CGAL_TESTS — no viewer, no wayland-scanner) - run: cmake -S code -B build -DWITH_CGAL_TESTS=ON -DCMAKE_BUILD_TYPE=Release + - name: Configure (LOW_MEMORY_BUILD — -O0, no PCH, unity batch 1) + run: | + cmake -S code -B build \ + -DWITH_CGAL_TESTS=ON \ + -DCMAKE_BUILD_TYPE=Release \ + -DCONFORMALLAB_LOW_MEMORY_BUILD=ON - name: Build CGAL-Tests run: nice -n 19 cmake --build build --target conformallab_cgal_tests -j1 @@ -118,38 +117,25 @@ jobs: echo "CGAL ▸ TOTAL ${total:-0} | PASSED $passed | FAILED ${failed:-0} | SKIPPED ${skipped:-0}" fi - # ── Structural gate: doc/api/tests.md totals match ctest reality ─── - # Single source of truth for test counts (see doc/release-policy.md). - # Reuses the already-built ./build dir via BUILD_DIR env var, so this - # adds ~5 s on top of the existing CGAL job. - name: Verify test-count consistency (doc/api/tests.md) run: BUILD_DIR=build bash scripts/check-test-counts.sh - # ── Structural gate: end-to-end smoke (try_it.sh) ────────────────── - # The user-facing quick-start script: configure + build + run the - # full ctest + run the Euclidean example on a bundled mesh. If - # this regresses, README quick-start instructions are broken. - # try_it.sh creates its own build-try/ — accept the ~3 min cost as - # the price of guaranteeing the documented workflow stays working. - name: End-to-end smoke test (scripts/try_it.sh) run: bash scripts/try_it.sh # ───────────────────────────────────────────────────────────────────────────── -# Job 3 — quality-gates (style + convention block) +# Job 3 — quality-gates # -# Cheap, deterministic checks that should never break unless a contributor -# introduces a regression. Each gate is a script under scripts/quality/ -# and exits 0 only when its tree is clean. These ran for weeks locally -# at zero findings before being promoted here. +# Trigger: include "/quality-gates" anywhere in the commit message. # -# Tools installed at job-start (the ci-cpp image already has python3 + -# bash; we add codespell + shellcheck on top). Total wall-time: ~30 s -# on the eulernest runner. +# git commit -m "chore: update docs /quality-gates" # -# Strictly required for merges into main/dev — a regression fails the PR. +# Cheap (~30 s): license headers, CGAL conventions, codespell, shellcheck. # ───────────────────────────────────────────────────────────────────────────── quality-gates: - needs: test-fast + if: | + contains(github.event.head_commit.message, '/quality-gates') || + contains(github.event.head_commit.message, '/ci-all') runs-on: eulernest container: image: git.eulernest.eu/conformallab/ci-cpp:latest diff --git a/.gitea/workflows/doc-build.yaml b/.gitea/workflows/doc-build.yaml index dc4216b..b66cfaa 100644 --- a/.gitea/workflows/doc-build.yaml +++ b/.gitea/workflows/doc-build.yaml @@ -1,28 +1,34 @@ name: API Docs +# Trigger: include "/docs" anywhere in the commit message. +# +# git commit -m "docs: update API examples /docs" +# +# Also available via workflow_dispatch for manual runs. + on: push: branches: - main - pull_request: + - dev + - "claude/**" + - "feature/**" + - "review/**" + workflow_dispatch: {} # ───────────────────────────────────────────────────────────────────────────── # Doc-build — informational only # # Generates Doxygen HTML from the public headers and reports warning # statistics. Does NOT block merges: `continue-on-error: true` ensures -# warnings or extraction issues never fail the CI gate. When Doxygen -# coverage is denser (Phase 8c), this job can be promoted to a hard -# requirement and the HTML deployed to Pages. -# -# Note: Gitea Actions on GHES does not support `actions/upload-artifact@v4`, -# so HTML artifact upload is intentionally omitted. The warning summary -# in the job log is the primary reviewer signal; reviewers who want the -# HTML can rebuild it locally with `cmake --build build --target doc`. +# warnings or extraction issues never fail. # ───────────────────────────────────────────────────────────────────────────── jobs: doc-build: - if: github.event_name == 'pull_request' + if: | + github.event_name == 'workflow_dispatch' || + contains(github.event.head_commit.message, '/docs') || + contains(github.event.head_commit.message, '/ci-all') runs-on: eulernest container: image: git.eulernest.eu/conformallab/ci-cpp:latest diff --git a/.gitea/workflows/markdown-links.yml b/.gitea/workflows/markdown-links.yml index b7affac..8bce97f 100644 --- a/.gitea/workflows/markdown-links.yml +++ b/.gitea/workflows/markdown-links.yml @@ -1,29 +1,34 @@ name: Markdown link check # Verify every internal markdown link in the repo resolves to an existing -# file (or anchor). External http(s) links are also probed but with a -# loose timeout — flaky third-party hosts must not break our CI. +# file (or anchor). # -# Trigger: PRs that touch any *.md file, plus a weekly cron so external -# link rot is caught even when nobody is editing docs. +# Triggers: +# - "/links" in commit message (manual, on that exact commit) +# - Weekly cron Mon 05:00 UTC (catches link rot without any activity) +# - workflow_dispatch (manual run on any branch) +# +# git commit -m "docs: rename section /links" on: - pull_request: - paths: - - "**/*.md" - - ".gitea/workflows/markdown-links.yml" push: branches: - main - paths: - - "**/*.md" - - ".gitea/workflows/markdown-links.yml" + - dev + - "claude/**" + - "feature/**" + - "review/**" schedule: - cron: "0 5 * * 1" # Monday 05:00 UTC weekly link-rot check workflow_dispatch: {} jobs: check: + if: | + github.event_name == 'schedule' || + github.event_name == 'workflow_dispatch' || + contains(github.event.head_commit.message, '/links') || + contains(github.event.head_commit.message, '/ci-all') runs-on: eulernest container: image: git.eulernest.eu/conformallab/ci-cpp:latest diff --git a/CLAUDE.md b/CLAUDE.md index d38f1b4..3f1df6e 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -56,6 +56,7 @@ The CGAL test build defaults to **PCH + Unity Build ON** (CGAL test wall-time 78 | `-DCMAKE_UNITY_BUILD=OFF` | ON | Disable Unity (jumbo) build. | | `-DCONFORMALLAB_DEV_BUILD=ON` | OFF | Dev iteration: PCH on, Unity forced off (cheaper incremental rebuilds). | | `-DCONFORMALLAB_FAST_TEST_BUILD=ON` | OFF | `-O0 -g` for tests (faster compile, slower run). | +| `-DCONFORMALLAB_LOW_MEMORY_BUILD=ON` | OFF | **RAM-constrained CI** (Raspberry Pi): `-O0` (no -g), PCH off, unity batch 1, `--no-keep-memory` linker. Drops cc1plus peak from ~700 MB to ~150-200 MB per TU so the CGAL build fits in a 2 GB container. Tests run ~15× slower but all pass. Use with `-j1`. | | `-DCONFORMALLAB_USE_CCACHE=ON` | OFF | Route compiles through ccache. | | `-DCONFORMALLAB_HEADERS_CHECK=ON` | OFF | Standalone header self-containment check target. | @@ -263,15 +264,17 @@ Three jobs in `.gitea/workflows/cpp-tests.yml`: | Job | CMake flags | Deps | Triggers on | Status | |---|---|---|---|---| -| `test-fast` | *(none)* | Eigen + GTest only | all branches | **active** | -| `test-cgal` | `-DWITH_CGAL_TESTS=ON` | + Boost | pull requests only | **DISABLED 2026-05-26** (`if: false`) | -| `quality-gates` | *(none)* | + codespell, shellcheck | all branches (`needs: test-fast`) | **active** | +| `test-fast` | *(none)* | Eigen + GTest only | all branches (auto) | **active** | +| `test-cgal` | `-DWITH_CGAL_TESTS=ON -DCONFORMALLAB_LOW_MEMORY_BUILD=ON` | + Boost | `/test-cgal` in commit message | **active** | +| `quality-gates` | *(none)* | + codespell, shellcheck | `/quality-gates` in commit message | **active** | +| `doc-build` | *(none)* | Doxygen | `/docs` in commit message or `workflow_dispatch` | **active** | +| `markdown-links` | *(none)* | python3 | `/links` in commit message, weekly cron, `workflow_dispatch` | **active** | Runner: `eulernest` — self-hosted Raspberry Pi, ARM64, Ubuntu 22.04. Docker image: `git.eulernest.eu/conformallab/ci-cpp:latest`. `test-cgal` and `quality-gates` both need `test-fast` to pass first (`needs: test-fast`). `quality-gates` runs four required structural gates: `license-headers.sh`, `cgal-conventions.py`, `codespell.sh`, `shellcheck.sh --strict`. Seven more gates (clang-format, cmake-format, cppcheck, sanitizers, clang-tidy, multi-compiler, reproducible-build) are local-only — see `scripts/quality/README.md`. -**`test-cgal` is gated off** (`if: false`, see the DISABLED-2026-05-26 comment block in the workflow): the `-j1` CGAL build is too memory-heavy for the 1.6 GB runner and OOMs intermittently without exposing real regressions. Consequence: the two structural sub-gates that lived inside it — `scripts/check-test-counts.sh` and `scripts/try_it.sh` — are currently **un-gated**; run them locally before tagging a release. Re-enable by restoring `if: github.event_name == 'pull_request'`. +**`test-cgal` is comment-triggered** (2026-05-31): write `/test-cgal` as a comment on any PR to start the CGAL suite manually. Not triggered on every push — the Pi runner (3-4 GB RAM, swap heavily loaded) cannot sustain a build on every WIP commit. `LOW_MEMORY_BUILD=ON` (-O0, no PCH, unity batch 1) keeps peak cc1plus RAM at ~150-200 MB, fitting in a 2000 MB container. All 277 tests pass in ~31 s run time. The two structural sub-gates (`scripts/check-test-counts.sh`, `scripts/try_it.sh`) still run after the test step. Two other workflows are also restricted to `workflow_dispatch:` only (auto-trigger disabled 2026-05-26 while the codeberg pages-branch push is being stabilised): - `.gitea/workflows/doxygen-pages.yml` — publishes Doxygen HTML + reviewer hub to the codeberg `pages` branch. diff --git a/code/data/off/torus_skewed_4x4.off b/code/data/off/torus_skewed_4x4.off new file mode 100644 index 0000000..b987cd2 --- /dev/null +++ b/code/data/off/torus_skewed_4x4.off @@ -0,0 +1,50 @@ +OFF +16 32 0 +0.0 0.0 0.0 +1.0 0.0 0.0 +2.0 0.0 0.0 +3.0 0.0 0.0 +-0.25 1.0 0.0 +0.75 1.0 0.0 +1.75 1.0 0.0 +2.75 1.0 0.0 +-0.5 2.0 0.0 +0.5 2.0 0.0 +1.5 2.0 0.0 +2.5 2.0 0.0 +-0.75 3.0 0.0 +0.25 3.0 0.0 +1.25 3.0 0.0 +2.25 3.0 0.0 +3 0 1 5 +3 0 5 4 +3 1 2 6 +3 1 6 5 +3 2 3 7 +3 2 7 6 +3 3 0 4 +3 3 4 7 +3 4 5 9 +3 4 9 8 +3 5 6 10 +3 5 10 9 +3 6 7 11 +3 6 11 10 +3 7 4 8 +3 7 8 11 +3 8 9 13 +3 8 13 12 +3 9 10 14 +3 9 14 13 +3 10 11 15 +3 10 15 14 +3 11 8 12 +3 11 12 15 +3 12 13 1 +3 12 1 0 +3 13 14 2 +3 13 2 1 +3 14 15 3 +3 14 3 2 +3 15 12 0 +3 15 0 3 diff --git a/code/include/clausen.hpp b/code/include/clausen.hpp index 9f11638..2c7c5a4 100644 --- a/code/include/clausen.hpp +++ b/code/include/clausen.hpp @@ -30,6 +30,12 @@ inline double csevl(double x, const double* cs, int n) noexcept { // Count Chebyshev terms needed so truncation error <= eta. // Corresponds to Java Clausen.inits(). +// +// Note on return value: the loop decrements n one extra time after the +// stopping condition is met, so the returned value is one less than the +// last index checked. Callers pass the result directly to csevl() as +// the term count, which evaluates terms [0, n-1] — this is intentional +// and matches the Java Clausen.inits() / csevl() contract exactly. inline int inits(const double* series, int n, double eta) noexcept { double err = 0.0; while (err <= eta && n-- != 0) { diff --git a/code/include/conformal_mesh.hpp b/code/include/conformal_mesh.hpp index 98f8626..e771a8f 100644 --- a/code/include/conformal_mesh.hpp +++ b/code/include/conformal_mesh.hpp @@ -37,6 +37,7 @@ #include #include +#include #include namespace conformallab { @@ -107,4 +108,23 @@ inline auto add_face_properties(ConformalMesh& mesh) return ftype; } +// ── Half-edge index helper ──────────────────────────────────────────────────── +// +// Convert a Halfedge_index to std::size_t for use as a vector subscript. +// Each functional previously duplicated this one-liner under a name like +// eucl_hidx / spher_hidx / hidx. Centralised here (MINOR-4 fix). +// +// Implementation: the double-cast uint32_t → size_t is intentional. CGAL 6.x +// Surface_mesh stores half-edge indices internally as 32-bit unsigned integers. +// Casting directly to size_t on a 64-bit system would produce the same result +// because CGAL index types use non-negative values, but the explicit intermediate +// cast documents the assumption and silences spurious sign-conversion warnings. + +/// Convert a `Halfedge_index` to `std::size_t` for vector subscript use. +/// Replaces the duplicated `eucl_hidx`, `spher_hidx`, `hidx` helpers. +inline std::size_t halfedge_to_index(Halfedge_index h) noexcept +{ + return static_cast(static_cast(h)); +} + } // namespace conformallab diff --git a/code/include/cp_euclidean_functional.hpp b/code/include/cp_euclidean_functional.hpp index 758a9c4..8267591 100644 --- a/code/include/cp_euclidean_functional.hpp +++ b/code/include/cp_euclidean_functional.hpp @@ -324,16 +324,19 @@ inline Eigen::SparseMatrix cp_euclidean_hessian(const ConformalMesh& return H; } -/// FD gradient check for the CP-Euclidean functional. Mirrors the -/// Java `FunctionalTest`; default `eps = 1e-5`, `tol = 1e-6`. +/// FD gradient check for the CP-Euclidean functional (central differences). +/// Uses the same **relative** error criterion as every other gradient check in +/// this library: `|analytic − fd| / max(1, |analytic|) < tol`. +/// Default `eps = 1e-5`, `tol = 1e-4` (matches Java `FunctionalTest`). inline bool gradient_check_cp_euclidean(const ConformalMesh& mesh, const std::vector& x, const CPEuclideanMaps& m, double eps = 1e-5, - double tol = 1e-6) + double tol = 1e-4) { auto G = cp_euclidean_gradient(mesh, x, m); const std::size_t n = G.size(); + bool ok = true; for (std::size_t i = 0; i < n; ++i) { std::vector xp = x, xm = x; @@ -341,28 +344,33 @@ inline bool gradient_check_cp_euclidean(const ConformalMesh& mesh, xm[i] -= eps; const double Ep = cp_euclidean_energy(mesh, xp, m); const double Em = cp_euclidean_energy(mesh, xm, m); - const double fd = (Ep - Em) / (2.0 * eps); - if (std::abs(G[i] - fd) > tol) { + const double fd = (Ep - Em) / (2.0 * eps); + const double err = std::abs(G[i] - fd); + const double scale = std::max(1.0, std::abs(G[i])); + if (err / scale > tol) { std::cerr << "[cp-euclidean] FD gradient mismatch at DOF " << i << ": analytic=" << G[i] << " FD=" << fd - << " diff=" << (G[i] - fd) << "\n"; - return false; + << " rel-err=" << (err / scale) << "\n"; + ok = false; } } - return true; + return ok; } /// FD Hessian check for the CP-Euclidean functional. Verifies analytic /// `H` column-by-column against `(G(x+εe_j) − G(x−εe_j)) / (2ε)`. +/// Uses the same **relative** error criterion as `hessian_check_euclidean`: +/// `|analytic − fd| / max(1, |analytic|) < tol`. inline bool hessian_check_cp_euclidean(const ConformalMesh& mesh, const std::vector& x, const CPEuclideanMaps& m, double eps = 1e-5, - double tol = 1e-5) + double tol = 1e-4) { const auto H = cp_euclidean_hessian(mesh, x, m); const int n = static_cast(H.rows()); + bool ok = true; for (int j = 0; j < n; ++j) { std::vector xp = x, xm = x; @@ -372,19 +380,21 @@ inline bool hessian_check_cp_euclidean(const ConformalMesh& mesh, auto Gm = cp_euclidean_gradient(mesh, xm, m); for (int i = 0; i < n; ++i) { - double fd = (Gp[static_cast(i)] - Gm[static_cast(i)]) - / (2.0 * eps); - double an = H.coeff(i, j); - if (std::abs(an - fd) > tol) { + double fd = (Gp[static_cast(i)] - Gm[static_cast(i)]) + / (2.0 * eps); + double an = H.coeff(i, j); + double err = std::abs(an - fd); + double scale = std::max(1.0, std::abs(an)); + if (err / scale > tol) { std::cerr << "[cp-euclidean] FD Hessian mismatch at (" << i << "," << j << "): analytic=" << an << " FD=" << fd - << " diff=" << (an - fd) << "\n"; - return false; + << " rel-err=" << (err / scale) << "\n"; + ok = false; } } } - return true; + return ok; } } // namespace conformallab diff --git a/code/include/euclidean_functional.hpp b/code/include/euclidean_functional.hpp index 37e766f..f8cac0c 100644 --- a/code/include/euclidean_functional.hpp +++ b/code/include/euclidean_functional.hpp @@ -41,6 +41,7 @@ #include "conformal_mesh.hpp" #include "constants.hpp" +#include "gauss_legendre.hpp" #include "euclidean_geometry.hpp" #include #include @@ -95,10 +96,14 @@ inline EuclideanMaps setup_euclidean_maps(ConformalMesh& mesh) /// Assign sequential DOF indices `0..n-1` to all vertices. /// -/// **Note:** does NOT pin a gauge vertex. For closed meshes the caller -/// must set one `m.v_idx[v] = -1` either before or after this call to -/// remove the rotational mode (the Newton solver's SparseQR fallback -/// will otherwise pick a minimum-norm solution but at higher cost). +/// **Note:** this overload assigns indices to ALL vertices unconditionally. +/// Any `v_idx` set before the call is overwritten. To pin a gauge vertex, +/// either use the two-argument overload below, or set `m.v_idx[v] = -1` +/// **after** this call. Pinning before this call has no effect. +/// +/// For closed meshes one gauge vertex should be pinned to remove the +/// rotational null mode (the Newton solver's SparseQR fallback handles an +/// unpinned closed mesh automatically, but at higher cost). inline int assign_euclidean_vertex_dof_indices(ConformalMesh& mesh, EuclideanMaps& m) { int idx = 0; @@ -106,6 +111,20 @@ inline int assign_euclidean_vertex_dof_indices(ConformalMesh& mesh, EuclideanMap return idx; } +/// Assign sequential DOF indices to all vertices, pinning `gauge` +/// (`m.v_idx[gauge] = -1`). Use this overload on closed meshes to fix +/// the rotational gauge mode in a single call. +/// +/// \returns The number of free DOFs assigned (`num_vertices − 1`). +inline int assign_euclidean_vertex_dof_indices(ConformalMesh& mesh, EuclideanMaps& m, + Vertex_index gauge) +{ + int idx = 0; + for (auto v : mesh.vertices()) + m.v_idx[v] = (v == gauge) ? -1 : idx++; + return idx; +} + /// Assign DOF indices for all vertices AND all edges (vertex-DOFs first, /// then edge-DOFs). Use this overload for the "cyclic" formulation that /// includes per-edge log-length DOFs (`λ_e`) on top of per-vertex scale @@ -155,11 +174,10 @@ static inline double eucl_dof_val(int idx, const std::vector& x) return idx >= 0 ? x[static_cast(idx)] : 0.0; } -/// Convert a CGAL half-edge index to a plain `std::size_t` for vector indexing. -static inline std::size_t eucl_hidx(Halfedge_index h) -{ - return static_cast(static_cast(h)); -} +// halfedge_to_index is defined in conformal_mesh.hpp (included above). +// The old local alias eucl_hidx is retained as a thin wrapper for now so +// call-sites below do not need touching; a follow-up can remove it. +static inline std::size_t eucl_hidx(Halfedge_index h) { return halfedge_to_index(h); } /// Compute the Euclidean-functional gradient G(x): /// * `G_v = Θ_v − Σ_faces α_v(face)` @@ -254,20 +272,8 @@ inline double euclidean_energy( const std::vector& x, const EuclideanMaps& m) { - static const double gl_s[10] = { - -0.9739065285171717, -0.8650633666889845, - -0.6794095682990244, -0.4333953941292472, - -0.1488743389816312, 0.1488743389816312, - 0.4333953941292472, 0.6794095682990244, - 0.8650633666889845, 0.9739065285171717 - }; - static const double gl_w[10] = { - 0.0666713443086881, 0.1494513491505806, - 0.2190863625159820, 0.2692667193099963, - 0.2955242247147529, 0.2955242247147529, - 0.2692667193099963, 0.2190863625159820, - 0.1494513491505806, 0.0666713443086881 - }; + const double* gl_s = gl10_nodes(); + const double* gl_w = gl10_weights(); const std::size_t n = x.size(); double E = 0.0; diff --git a/code/include/euclidean_hessian.hpp b/code/include/euclidean_hessian.hpp index 5918c79..83f473d 100644 --- a/code/include/euclidean_hessian.hpp +++ b/code/include/euclidean_hessian.hpp @@ -17,18 +17,21 @@ // │ side lengths lij = exp(Λ̃ij/2): │ // │ │ // │ t12 = −l12+l23+l31, t23 = l12−l23+l31, t31 = l12+l23−l31 │ -// │ denom2 = 2·sqrt(t12·t23·t31·l123) = 8·Area │ +// │ l123 = l12+l23+l31, denom2 = 2·sqrt(t12·t23·t31·l123) = 8·Area │ // │ │ -// │ cot_k = (t_adj1·l123 − t_adj2·t_opp) / denom2 │ -// │ = cotangent of the angle αk at vertex k │ +// │ Cotangent at vertex k (opposite t_opp, adjacent t_a and t_b): │ +// │ cot_k = (t_opp · l123 − t_a · t_b) / denom2 │ // │ │ -// │ Hessian contributions per face: │ -// │ H[vi, vi] += cot_vj + cot_vk (diagonal, both non-opp angles) │ -// │ H[vi, vj] -= cot_vk (off-diagonal, for variable vi,vj│ +// │ Assignment (k ↔ opposite edge ↔ opposite t-value): │ +// │ cot1 (v1, opp l23): t_opp=t23, t_a=t12, t_b=t31 │ +// │ cot2 (v2, opp l31): t_opp=t31, t_a=t12, t_b=t23 │ +// │ cot3 (v3, opp l12): t_opp=t12, t_a=t23, t_b=t31 │ // │ │ -// │ This is exactly the cotangent-Laplace operator from Pinkall–Polthier. │ +// │ Hessian contributions per face (Pinkall–Polthier ½ factor): │ +// │ H[vi, vi] += (cot_vj + cot_vk) / 2 (diagonal) │ +// │ H[vi, vj] −= cot_vk / 2 (off-diagonal, variable pairs) │ // │ │ -// │ Pinned vertices (v_idx = −1) contribute to diagonal of neighbours but │ +// │ Pinned vertices (v_idx = −1) contribute to diagonal of neighbours but │ // │ do not create a column/row in H themselves. │ // └──────────────────────────────────────────────────────────────────────────┘ // @@ -54,7 +57,11 @@ namespace conformallab { // Given three Euclidean SIDE LENGTHS l12, l23, l31 (already exp(Λ̃/2)), // return the three cotangent weights (cot1, cot2, cot3). // -// cot_k = (t_adj·l123 − t_opp·t_other) / (8·Area) +// cot_k = (t_opp · l123 − t_a · t_b) / (8·Area) +// +// where t_opp is the t-value of the edge OPPOSITE vertex k, and t_a, t_b are +// the t-values of the two edges ADJACENT to vertex k. See the box comment at +// the top of this file for the explicit assignment of t_opp/t_a/t_b per vertex. // // Returns {0,0,0} for degenerate faces (triangle inequality violated or Area=0). /// Three Euclidean cotangent weights `(cot1, cot2, cot3)` for the @@ -85,9 +92,10 @@ inline EuclCotWeights euclidean_cot_weights(double l12, double l23, double l31) // denom2 = 2·sqrt(t12·t23·t31·l123) = 8·Area const double denom2 = 2.0 * std::sqrt(denom2_sq); - // cot at v1 (opposite l23): adjacent t-values are t12 and t31. - // cot at v2 (opposite l31): adjacent t-values are t12 and t23. - // cot at v3 (opposite l12): adjacent t-values are t23 and t31. + // Formula: cot_k = (t_opp · l123 − t_a · t_b) / denom2 + // cot1: t_opp=t23, t_a=t12, t_b=t31 (v1 opposite l23) + // cot2: t_opp=t31, t_a=t12, t_b=t23 (v2 opposite l31) + // cot3: t_opp=t12, t_a=t23, t_b=t31 (v3 opposite l12) return { (t23 * l123 - t31 * t12) / denom2, // cot1 (t31 * l123 - t12 * t23) / denom2, // cot2 diff --git a/code/include/gauss_bonnet.hpp b/code/include/gauss_bonnet.hpp index 1d92c68..8ac316f 100644 --- a/code/include/gauss_bonnet.hpp +++ b/code/include/gauss_bonnet.hpp @@ -7,14 +7,28 @@ // Phase 6 — Gauss–Bonnet consistency check for prescribed target angles. // // Before calling newton_*() with custom target angles, verify that -// the angle defect sum matches the topology: +// the angle defect sum matches the topology. // -// Σ_v (2π − Θ_v) = 2π · χ(M) (Euclidean / flat) -// Σ_v (2π − Θ_v) > 0 (spherical, χ > 0) -// Σ_v (2π − Θ_v) < 0 (hyperbolic, χ < 0) +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Geometry Identity to satisfy │ +// │ ───────────────────────────────────────────────────────────────────── │ +// │ Euclidean/flat Σ_v (2π − Θ_v) = 2π · χ(M) (exact equality) │ +// │ Spherical Σ_v (2π − Θ_v) > 0 (sufficient, χ > 0) │ +// │ │ +// │ HyperIdeal — NOT SUPPORTED by this header. │ +// │ The correct hyperbolic Gauss–Bonnet identity is │ +// │ Σ_v (2π − Θ_v) − Area(M) = 2π · χ(M) │ +// │ which differs from the Euclidean identity by the Area(M) > 0 term. │ +// │ Computing Area(M) from the HyperIdeal DOFs is non-trivial. │ +// │ gauss_bonnet_sum(mesh, HyperIdealMaps) and │ +// │ enforce_gauss_bonnet(mesh, HyperIdealMaps) are therefore DELETED. │ +// │ Do NOT call check_gauss_bonnet before newton_hyper_ideal — │ +// │ it is not needed; the HyperIdeal energy is strictly convex so Newton │ +// │ converges without a pre-check. │ +// └─────────────────────────────────────────────────────────────────────────┘ // -// If this fails, no conformal factor can realise the target angles and -// Newton will silently fail to converge. +// If the Euclidean/Spherical check fails, no conformal factor can realise +// the target angles and Newton will silently fail to converge. // // PRECONDITION — closed meshes only. Every function here sums (2π − Θ_v) // over ALL vertices. On a mesh with boundary the boundary vertices carry a @@ -26,11 +40,12 @@ // API: // int euler_characteristic(mesh) // int genus(mesh) -// double gauss_bonnet_sum(mesh, maps) — Σ(2π − Θ_v) +// double gauss_bonnet_sum(mesh, EuclideanMaps/SphericalMaps) — Σ(2π − Θ_v) // double gauss_bonnet_rhs(mesh) — 2π · χ(M) // double gauss_bonnet_deficit(mesh, maps) — lhs − rhs (0 = satisfied) // void check_gauss_bonnet(mesh, maps [, tol]) — throws if violated // void enforce_gauss_bonnet(mesh, maps) — shifts θ_v by uniform Δ +// (HyperIdealMaps overloads are deleted — see box above) #include "conformal_mesh.hpp" #include "euclidean_functional.hpp" @@ -78,15 +93,23 @@ inline double gauss_bonnet_sum( } /// `gauss_bonnet_sum` for the Euclidean-functional property bundle. -inline double gauss_bonnet_sum(const ConformalMesh& m, const EuclideanMaps& mp) +inline double gauss_bonnet_sum(const ConformalMesh& m, const EuclideanMaps& mp) { return gauss_bonnet_sum(m, mp.theta_v); } /// `gauss_bonnet_sum` for the Spherical-functional property bundle. -inline double gauss_bonnet_sum(const ConformalMesh& m, const SphericalMaps& mp) - { return gauss_bonnet_sum(m, mp.theta_v); } -/// `gauss_bonnet_sum` for the HyperIdeal-functional property bundle. -inline double gauss_bonnet_sum(const ConformalMesh& m, const HyperIdealMaps& mp) +inline double gauss_bonnet_sum(const ConformalMesh& m, const SphericalMaps& mp) { return gauss_bonnet_sum(m, mp.theta_v); } +// gauss_bonnet_sum for HyperIdealMaps is intentionally DELETED. +// The correct hyperbolic Gauss–Bonnet identity is +// Σ(2π−Θ_v) − Area(M) = 2π·χ(M) +// not the Euclidean form Σ(2π−Θ_v) = 2π·χ(M). Providing this overload +// would silently skip the Area term, making check_gauss_bonnet always +// fail for valid hyperbolic targets (e.g. a genus-2 mesh with Θ_v=2π +// gives Σ(2π−Θ_v)=0 but 2π·χ=−4π → deficit=4π ≠ 0 every time). +// Use newton_hyper_ideal directly — no pre-check is needed because the +// HyperIdeal energy is strictly convex (Springborn 2020 Theorem 1.3). +inline double gauss_bonnet_sum(const ConformalMesh&, const HyperIdealMaps&) = delete; + // ── Right-hand side 2π · χ(M) ─────────────────────────────────────────────── /// Right-hand side of Gauss-Bonnet: `2π · χ(M)`. @@ -157,10 +180,18 @@ inline void enforce_gauss_bonnet( } /// Distribute the Gauss-Bonnet deficit uniformly across `maps.theta_v`. +/// Supported for EuclideanMaps and SphericalMaps only. +/// HyperIdealMaps overload is deleted — see header comment for why. template inline void enforce_gauss_bonnet(ConformalMesh& mesh, Maps& maps) { enforce_gauss_bonnet(mesh, maps.theta_v); } +// enforce_gauss_bonnet for HyperIdealMaps is intentionally DELETED. +// The Euclidean identity Σ(2π−Θ_v)=2π·χ is not the correct pre-condition +// for HyperIdeal. Calling this function would silently shift Θ_v to +// satisfy the wrong identity, producing incorrect target angles. +inline void enforce_gauss_bonnet(ConformalMesh&, HyperIdealMaps&) = delete; + } // namespace conformallab diff --git a/code/include/gauss_legendre.hpp b/code/include/gauss_legendre.hpp new file mode 100644 index 0000000..e80769a --- /dev/null +++ b/code/include/gauss_legendre.hpp @@ -0,0 +1,49 @@ +#pragma once +// Copyright (c) 2024-2026 Tarik Moussa. +// SPDX-License-Identifier: MIT + +// gauss_legendre.hpp +// +// 10-point Gauss-Legendre quadrature nodes and weights on [-1,1]. +// +// Previously duplicated verbatim in: +// euclidean_functional.hpp, spherical_functional.hpp, +// inversive_distance_functional.hpp (MINOR-3 fix) +// +// Usage (integration over [0,1] via change of variables t=(1+s)/2, w=w_GL/2): +// +// const auto* s = conformallab::gl10_nodes(); +// const auto* w = conformallab::gl10_weights(); +// for (int k = 0; k < 10; ++k) { +// double t = (1.0 + s[k]) * 0.5; +// double wt = w[k] * 0.5; +// E += wt * dot(G(t*x), x); +// } + +namespace conformallab { + +/// 10-point Gauss-Legendre nodes on [−1, 1]. +inline const double* gl10_nodes() noexcept { + static constexpr double s[10] = { + -0.9739065285171717, -0.8650633666889845, + -0.6794095682990244, -0.4333953941292472, + -0.1488743389816312, 0.1488743389816312, + 0.4333953941292472, 0.6794095682990244, + 0.8650633666889845, 0.9739065285171717 + }; + return s; +} + +/// 10-point Gauss-Legendre weights on [−1, 1]. +inline const double* gl10_weights() noexcept { + static constexpr double w[10] = { + 0.0666713443086881, 0.1494513491505806, + 0.2190863625159820, 0.2692667193099963, + 0.2955242247147529, 0.2955242247147529, + 0.2692667193099963, 0.2190863625159820, + 0.1494513491505806, 0.0666713443086881 + }; + return w; +} + +} // namespace conformallab diff --git a/code/include/hyper_ideal_functional.hpp b/code/include/hyper_ideal_functional.hpp index 814f9d0..7a5f125 100644 --- a/code/include/hyper_ideal_functional.hpp +++ b/code/include/hyper_ideal_functional.hpp @@ -34,6 +34,7 @@ #include "hyper_ideal_geometry.hpp" #include "hyper_ideal_utility.hpp" #include +#include #include #include #include @@ -128,11 +129,8 @@ static inline double dof_val(int idx, const std::vector& x) return idx >= 0 ? x[static_cast(idx)] : 0.0; } -/// Convert a CGAL half-edge index to a plain `std::size_t` for vector indexing. -static inline std::size_t hidx(Halfedge_index h) -{ - return static_cast(static_cast(h)); -} +// halfedge_to_index is defined in conformal_mesh.hpp. +static inline std::size_t hidx(Halfedge_index h) { return halfedge_to_index(h); } // ── Pure-math face-angle kernel ────────────────────────────────────────────── // @@ -317,25 +315,54 @@ static FaceAngles compute_face_angles( } /// Per-face energy contribution U(f) before subtracting the θ·a and Θ·b terms. +/// +/// Supported configurations (faithful port of HyperIdealFunctional.java): +/// * All three vertices hyper-ideal (v?b = true) → Meyerhoff/Ushijima volume +/// * Exactly one vertex ideal (v?b = false, other two true) → Kolpakov-Mednykh volume +/// +/// NOT supported — faces with two or three ideal vertices. The Java reference +/// (HyperIdealFunctional.java lines 222-231) uses an if/else-if chain that +/// silently applies the one-ideal-vertex formula to the first ideal vertex it +/// finds, ignoring additional ideal vertices. That is mathematically wrong for +/// two-ideal or three-ideal faces. Rather than silently computing a wrong result, +/// this C++ port throws immediately so the caller can diagnose the problem. +/// The correct volume formulas for semi-ideal and fully-ideal faces are not +/// implemented in the Java reference and would require new research. static double face_energy(const FaceAngles& fa) { + // Guard: reject configurations with 2 or 3 ideal vertices in one face. + const int ideal_count = (!fa.v1b ? 1 : 0) + + (!fa.v2b ? 1 : 0) + + (!fa.v3b ? 1 : 0); + if (ideal_count >= 2) + throw std::logic_error( + "face_energy: faces with 2 or 3 ideal (pinned) vertices are not " + "supported. Only 0-ideal (all hyper-ideal) and 1-ideal faces are " + "implemented, matching the Java HyperIdealFunctional reference. " + "Check your v_idx assignments: at most one vertex per face may be " + "pinned (v_idx = -1)."); + double aa = fa.a12*fa.alpha12 + fa.a23*fa.alpha23 + fa.a31*fa.alpha31; double bb = fa.b1 *fa.beta1 + fa.b2 *fa.beta2 + fa.b3 *fa.beta3; double V = 0.0; if (fa.v1b && fa.v2b && fa.v3b) { + // All three vertices are hyper-ideal. V = calculateTetrahedronVolume( fa.beta1, fa.beta2, fa.beta3, fa.alpha23, fa.alpha31, fa.alpha12); } else if (!fa.v1b) { + // Exactly v1 is ideal (ideal_count == 1 guaranteed by guard above). V = calculateTetrahedronVolumeWithIdealVertexAtGamma( fa.beta1, fa.alpha31, fa.alpha12, fa.alpha23, fa.beta2, fa.beta3); } else if (!fa.v2b) { + // Exactly v2 is ideal. V = calculateTetrahedronVolumeWithIdealVertexAtGamma( fa.beta2, fa.alpha12, fa.alpha23, fa.alpha31, fa.beta3, fa.beta1); - } else { // !v3b + } else { + // Exactly v3 is ideal (!v3b, guaranteed by ideal_count == 1). V = calculateTetrahedronVolumeWithIdealVertexAtGamma( fa.beta3, fa.alpha23, fa.alpha31, fa.alpha12, fa.beta1, fa.beta2); diff --git a/code/include/inversive_distance_functional.hpp b/code/include/inversive_distance_functional.hpp index 2d48e3e..8d67a32 100644 --- a/code/include/inversive_distance_functional.hpp +++ b/code/include/inversive_distance_functional.hpp @@ -68,6 +68,7 @@ #include "conformal_mesh.hpp" #include "constants.hpp" +#include "gauss_legendre.hpp" #include "euclidean_geometry.hpp" // euclidean_angles(λ12, λ23, λ31) #include #include @@ -122,11 +123,10 @@ inline InversiveDistanceMaps setup_inversive_distance_maps(ConformalMesh& mesh) /// Assign sequential DOF indices `0..n-1` to every vertex. /// -/// **Note:** this overload does NOT pin a gauge vertex. The caller -/// is expected to either: -/// 1. set one `m.v_idx[v] = -1` *before* calling this function (then -/// the call is a no-op for that vertex) — OR — -/// 2. flip one assigned index back to `-1` *after* this function. +/// **Note:** this overload assigns indices to ALL vertices unconditionally. +/// Any `v_idx` set before the call is overwritten. To pin a gauge vertex, +/// either use the two-argument overload below, or set `m.v_idx[v] = -1` +/// **after** this call. Pinning before this call has no effect. /// /// For a closed mesh, exactly one pin is required to remove the /// global rotational mode. @@ -138,6 +138,21 @@ inline int assign_inversive_distance_vertex_dof_indices(ConformalMesh& m return idx; } +/// Assign sequential DOF indices to all vertices, pinning `gauge` +/// (`m.v_idx[gauge] = -1`). Use this overload on closed meshes to fix +/// the rotational gauge mode in a single call. +/// +/// \returns The number of free DOFs assigned (`num_vertices − 1`). +inline int assign_inversive_distance_vertex_dof_indices(ConformalMesh& mesh, + InversiveDistanceMaps& m, + Vertex_index gauge) +{ + int idx = 0; + for (auto v : mesh.vertices()) + m.v_idx[v] = (v == gauge) ? -1 : idx++; + return idx; +} + /// Count the free DOFs (vertices with `v_idx >= 0`). inline int inversive_distance_dimension(const ConformalMesh& mesh, const InversiveDistanceMaps& m) @@ -215,10 +230,8 @@ inline double dof_val(int idx, const std::vector& x) noexcept return idx >= 0 ? x[static_cast(idx)] : 0.0; } -inline std::size_t hidx(Halfedge_index h) noexcept -{ - return static_cast(static_cast(h)); -} +// halfedge_to_index is defined in conformal_mesh.hpp. +inline std::size_t hidx(Halfedge_index h) noexcept { return halfedge_to_index(h); } // Inversive-distance edge length squared: ℓ² = exp(2u_i) + exp(2u_j) + 2 I r_i r_j // where r_i = exp(u_i), so: ℓ² = r_i² + r_j² + 2 I r_i r_j. @@ -309,20 +322,8 @@ inline double inversive_distance_energy( const std::vector& x, const InversiveDistanceMaps& m) { - static const double gl_s[10] = { - -0.9739065285171717, -0.8650633666889845, - -0.6794095682990244, -0.4333953941292472, - -0.1488743389816312, 0.1488743389816312, - 0.4333953941292472, 0.6794095682990244, - 0.8650633666889845, 0.9739065285171717 - }; - static const double gl_w[10] = { - 0.0666713443086881, 0.1494513491505806, - 0.2190863625159820, 0.2692667193099963, - 0.2955242247147529, 0.2955242247147529, - 0.2692667193099963, 0.2190863625159820, - 0.1494513491505806, 0.0666713443086881 - }; + const double* gl_s = gl10_nodes(); + const double* gl_w = gl10_weights(); const std::size_t n = x.size(); double E = 0.0; @@ -340,32 +341,37 @@ inline double inversive_distance_energy( } /// FD gradient check for the Inversive-Distance functional (central diff). +/// Uses the same **relative** error criterion as every other gradient check: +/// `|analytic − fd| / max(1, |analytic|) < tol`. inline bool gradient_check_inversive_distance( const ConformalMesh& mesh, const std::vector& x, const InversiveDistanceMaps& m, double eps = 1e-5, - double tol = 1e-6) + double tol = 1e-4) { auto G = inversive_distance_gradient(mesh, x, m); const std::size_t n = G.size(); + bool ok = true; for (std::size_t i = 0; i < n; ++i) { std::vector xp = x, xm = x; xp[i] += eps; xm[i] -= eps; - double Ep = inversive_distance_energy(mesh, xp, m); - double Em = inversive_distance_energy(mesh, xm, m); - double fd = (Ep - Em) / (2.0 * eps); - if (std::abs(G[i] - fd) > tol) { + double Ep = inversive_distance_energy(mesh, xp, m); + double Em = inversive_distance_energy(mesh, xm, m); + double fd = (Ep - Em) / (2.0 * eps); + double err = std::abs(G[i] - fd); + double scale = std::max(1.0, std::abs(G[i])); + if (err / scale > tol) { std::cerr << "[inversive-distance] FD gradient mismatch at DOF " << i << ": analytic=" << G[i] << " FD=" << fd - << " diff=" << (G[i] - fd) << "\n"; - return false; + << " rel-err=" << (err / scale) << "\n"; + ok = false; } } - return true; + return ok; } /// Newton equilibrium check: returns `true` iff the gradient at `x` diff --git a/code/include/spherical_functional.hpp b/code/include/spherical_functional.hpp index 0355c58..5813e66 100644 --- a/code/include/spherical_functional.hpp +++ b/code/include/spherical_functional.hpp @@ -36,6 +36,7 @@ #include "conformal_mesh.hpp" #include "spherical_geometry.hpp" +#include "gauss_legendre.hpp" #include #include #include @@ -87,7 +88,14 @@ inline SphericalMaps setup_spherical_maps(ConformalMesh& mesh) } /// Assign sequential DOF indices `0..n-1` to all vertices (no edge DOFs). -/// Caller is expected to pin one gauge vertex with `m.v_idx[v] = -1`. +/// +/// **Note:** this overload assigns indices to ALL vertices unconditionally. +/// Any `v_idx` set before the call is overwritten. To pin a gauge vertex, +/// either use the two-argument overload below, or set `m.v_idx[v] = -1` +/// **after** this call. Pinning before this call has no effect. +/// +/// On closed spherical surfaces exactly one gauge vertex must be pinned +/// to remove the global-scale null mode. inline int assign_vertex_dof_indices(ConformalMesh& mesh, SphericalMaps& m) { int idx = 0; @@ -95,6 +103,20 @@ inline int assign_vertex_dof_indices(ConformalMesh& mesh, SphericalMaps& m) return idx; } +/// Assign sequential DOF indices to all vertices, pinning `gauge` +/// (`m.v_idx[gauge] = -1`). Use this overload on closed spherical +/// surfaces to fix the global-scale gauge mode in a single call. +/// +/// \returns The number of free DOFs assigned (`num_vertices − 1`). +inline int assign_vertex_dof_indices(ConformalMesh& mesh, SphericalMaps& m, + Vertex_index gauge) +{ + int idx = 0; + for (auto v : mesh.vertices()) + m.v_idx[v] = (v == gauge) ? -1 : idx++; + return idx; +} + /// Assign DOF indices for all vertices AND all edges (vertex-DOFs first, /// then edge-DOFs). Mirrors `assign_euclidean_all_dof_indices` for the /// cyclic spherical formulation. @@ -175,11 +197,8 @@ static inline double spher_eff_lambda(const SphericalMaps& m, : (m.lambda0[e] + u_i + u_j); } -/// Convert a CGAL half-edge index to a plain `std::size_t` for vector indexing. -static inline std::size_t spher_hidx(Halfedge_index h) -{ - return static_cast(static_cast(h)); -} +// halfedge_to_index is defined in conformal_mesh.hpp. +static inline std::size_t spher_hidx(Halfedge_index h) { return halfedge_to_index(h); } // ── Gradient only (no energy) ───────────────────────────────────────────────── @@ -300,21 +319,8 @@ inline double spherical_energy( 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 double* gl_s = gl10_nodes(); + const double* gl_w = gl10_weights(); const std::size_t n = x.size(); double E = 0.0; @@ -402,8 +408,11 @@ inline bool gradient_check_spherical( // Apply the shift by adding t* to every vertex DOF in x. // // Implementation: bisection on f(t) = Σ_v G_v(x + t·1_v). -// f is strictly monotone decreasing (second derivative < 0) for a convex -// functional, so bisection converges in O(log₂(2·bracket/tol)) iterations. +// f is strictly monotone decreasing because increasing the global scale +// increases all effective edge lengths and thereby all corner angles, which +// reduces Σ G_v = Σ(Θ_v − Σα_v). This holds for the concave spherical +// energy (NSD Hessian) just as well as for a convex one. +// Bisection converges in O(log₂(2·bracket/tol)) iterations. // // Parameters: // bracket – initial search interval [−bracket, +bracket] (default 50) @@ -456,7 +465,7 @@ inline double spherical_gauge_shift( // ── No sign change (zero may lie at a domain boundary). ─────────────────── // Use damped Newton's method with backtracking line search. - // f'(t) estimated by forward finite difference. + // f'(t) estimated by central finite difference (O(ε²) vs O(ε) for forward). // When the Newton step overshoots the valid domain (ΣG_v jumps back up // because faces become degenerate), backtracking halves the step until // |f| strictly decreases. @@ -467,8 +476,7 @@ inline double spherical_gauge_shift( for (int iter = 0; iter < 120; ++iter) { if (std::abs(ft) < tol) return t; - double ftp = sum_Gv(t + fd_eps); - double dft = (ftp - ft) / fd_eps; + double dft = (sum_Gv(t + fd_eps) - sum_Gv(t - fd_eps)) / (2.0 * fd_eps); if (std::abs(dft) < 1e-14) return t; // gradient flat — give up double dt_raw = -ft / dft; diff --git a/code/tests/cgal/CMakeLists.txt b/code/tests/cgal/CMakeLists.txt index 1707d50..d37ddd6 100644 --- a/code/tests/cgal/CMakeLists.txt +++ b/code/tests/cgal/CMakeLists.txt @@ -133,6 +133,67 @@ if(CONFORMALLAB_FAST_TEST_BUILD) ) endif() +# ── Low-memory build mode (for RAM-constrained CI runners, e.g. Raspberry Pi) ── +# +# Problem: CGAL + Eigen at -O3 drives cc1plus peak RAM to ~600-800 MB per +# Unity compilation unit on ARM64 Linux. A 1600 MB container limit with a +# batch of 4 files per unit causes OOM-kill during the build. +# +# This flag enables three orthogonal memory-saving measures: +# +# 1. -O0 (no debug info): drops cc1plus backend RAM by ~60-70 %. +# Optimizer passes (inlining, register allocation, constant propagation) +# dominate the backend. At -O0 they are entirely skipped → peak per +# TU falls from ~700 MB to ~150-200 MB on ARM64. +# We omit -g deliberately: debug info adds ~30-40 % object-file size +# and increases linker RSS. CI needs "does it compile + do tests pass", +# not debuggability. +# +# 2. PCH OFF: the precompiled header itself consumes ~200 MB to compile +# and is re-read by every TU. Disabling it saves the one-time PCH +# compilation cost; each TU re-parses CGAL headers, but at -O0 this +# is fast. +# +# 3. UNITY_BUILD_BATCH_SIZE=1: one source file per unity unit. Removes +# the "4 files × CGAL parse cost" multiplier; each cc1plus process +# only sees one file worth of templates. +# +# 4. Linker memory flag (GCC/Clang only): --no-keep-memory tells GNU ld +# to release symbol table memory after each input file instead of +# keeping it for cross-reference. Reduces linker RSS by 15-25 % at +# the cost of slightly longer link time. +# +# Activate with: +# cmake -S code -B build -DWITH_CGAL_TESTS=ON \ +# -DCONFORMALLAB_LOW_MEMORY_BUILD=ON +# Expected peak per cc1plus: ~150-200 MB → fits in 1800-2000 MB container. +# Test runtime is 2-4× slower than Release (CGAL traversals unoptimized), +# but correctness is unaffected. +option(CONFORMALLAB_LOW_MEMORY_BUILD + "Build CGAL tests with -O0, no PCH, UNITY_BATCH_SIZE=1 for RAM-constrained CI." OFF) + +if(CONFORMALLAB_LOW_MEMORY_BUILD) + message(STATUS "CONFORMALLAB_LOW_MEMORY_BUILD active: -O0, no PCH, unity batch 1.") + + # 1. -O0, no debug info + target_compile_options(conformallab_cgal_tests PRIVATE + $<$:-O0 -UNDEBUG> + ) + + # 2. PCH off — force-override the option so the block below is skipped + set(CONFORMALLAB_USE_PCH OFF CACHE BOOL "" FORCE) + + # 3. Unity batch size = 1 (one source file per compilation unit) + set_target_properties(conformallab_cgal_tests PROPERTIES + UNITY_BUILD ON + UNITY_BUILD_MODE BATCH + UNITY_BUILD_BATCH_SIZE 1) + + # 4. Linker memory hint (GNU ld / lld) + target_link_options(conformallab_cgal_tests PRIVATE + $<$:-Wl,--no-keep-memory>) +endif() + target_link_libraries(conformallab_cgal_tests PRIVATE GTest::gtest_main) # ── Compile-time speed-up: precompiled headers ─────────────────────────────── diff --git a/code/tests/cgal/test_euclidean_functional.cpp b/code/tests/cgal/test_euclidean_functional.cpp index 9857ae5..f1f4063 100644 --- a/code/tests/cgal/test_euclidean_functional.cpp +++ b/code/tests/cgal/test_euclidean_functional.cpp @@ -125,7 +125,16 @@ TEST(EuclideanFunctional, AngleSumEqualsPi) } // ════════════════════════════════════════════════════════════════════════════ -// Degenerate triangle → valid = false +// Degenerate triangle — limiting angles (Finding-F, java-port-audit item 1) +// +// The BPS-energy convex C¹ extension assigns the *limiting* angles when the +// triangle inequality is violated: the corner OPPOSITE the over-long edge +// gets π, the other two get 0. These tests lock that behaviour in for +// both euclidean_angles_from_lengths() and the gradient accumulation. +// +// Before the java-port-audit Finding 1 fix, the degenerate-face code +// returned {0,0,0} and the gradient skipped the face entirely, producing +// the wrong gradient on near-flip configurations. // ════════════════════════════════════════════════════════════════════════════ TEST(EuclideanFunctional, DegenerateTriangleReturnsFalse) @@ -135,6 +144,88 @@ TEST(EuclideanFunctional, DegenerateTriangleReturnsFalse) EXPECT_FALSE(fa.valid); } +TEST(EuclideanFunctional, DegenerateTriangle_LimitingAngles_L23TooLong) +{ + // l23 = 10 >> l12 + l31 = 2 → α₁ = π (at v1, opposite l23), α₂=α₃=0. + auto fa = euclidean_angles_from_lengths(1.0, 10.0, 1.0); + EXPECT_FALSE(fa.valid); + EXPECT_NEAR(fa.alpha1, PI, 1e-12) << "corner opposite over-long l23 must be π"; + EXPECT_NEAR(fa.alpha2, 0.0, 1e-12); + EXPECT_NEAR(fa.alpha3, 0.0, 1e-12); +} + +TEST(EuclideanFunctional, DegenerateTriangle_LimitingAngles_L31TooLong) +{ + // l31 too long → α₂ = π (at v2, opposite l31). + auto fa = euclidean_angles_from_lengths(1.0, 1.0, 10.0); + EXPECT_FALSE(fa.valid); + EXPECT_NEAR(fa.alpha2, PI, 1e-12) << "corner opposite over-long l31 must be π"; + EXPECT_NEAR(fa.alpha1, 0.0, 1e-12); + EXPECT_NEAR(fa.alpha3, 0.0, 1e-12); +} + +TEST(EuclideanFunctional, DegenerateTriangle_LimitingAngles_L12TooLong) +{ + // l12 too long → α₃ = π (at v3, opposite l12). + auto fa = euclidean_angles_from_lengths(10.0, 1.0, 1.0); + EXPECT_FALSE(fa.valid); + EXPECT_NEAR(fa.alpha3, PI, 1e-12) << "corner opposite over-long l12 must be π"; + EXPECT_NEAR(fa.alpha1, 0.0, 1e-12); + EXPECT_NEAR(fa.alpha2, 0.0, 1e-12); +} + +TEST(EuclideanFunctional, DegenerateTriangle_GradientPicksUpPiCorner) +{ + // Build a single triangle. Force a degenerate effective-length by + // setting a large negative lambda0 on two edges so l12 >> l23 + l31. + // + // DOFs: all vertices free. x = 0 (no conformal scaling). + // lambda0: e_opp_v3 (i.e. l12) is huge; the other two are near-zero. + // Expected: the gradient at v3 (opposite l12) picks up −π from the + // degenerate face; G_v3 = Θ_v3 − π = 2π − π = π. + auto mesh = make_triangle(); + auto maps = setup_euclidean_maps(mesh); + + // Identify edges: h0=halfedge(face), source(h0)=v1, source(next(h0))=v2, etc. + auto f = *mesh.faces().begin(); + auto h0 = mesh.halfedge(f); + auto h1 = mesh.next(h0); + auto h2 = mesh.next(h1); + + // Assign large lambda0 to the edge opposite v3 (= edge of h0, i.e. e12). + Edge_index e12 = mesh.edge(h0); + Edge_index e23 = mesh.edge(h1); + Edge_index e31 = mesh.edge(h2); + maps.lambda0[e12] = 20.0; // l12 = exp(10) ≈ 22026 — hugely over-long + maps.lambda0[e23] = 0.0; + maps.lambda0[e31] = 0.0; + + int n = assign_euclidean_vertex_dof_indices(mesh, maps); + std::vector x(static_cast(n), 0.0); + + auto G = euclidean_gradient(mesh, x, maps); + + // v3 = source(h2). Its gradient component should include the π corner. + Vertex_index v3 = mesh.source(h2); + int iv3 = maps.v_idx[v3]; + ASSERT_GE(iv3, 0); + + // G_v3 = Θ_v3 − α3. The degenerate face gives α3 = π. + // Θ_v3 defaults to 2π, so G_v3 = 2π − π = π. + EXPECT_NEAR(G[static_cast(iv3)], PI, 1e-10) + << "Gradient at v3 must include the π limiting angle from the degenerate face"; + + // v1 and v2 get α = 0 from the degenerate face → G_vi = Θ − 0 = 2π. + Vertex_index v1 = mesh.source(h0); + Vertex_index v2 = mesh.source(h1); + int iv1 = maps.v_idx[v1], iv2 = maps.v_idx[v2]; + ASSERT_GE(iv1, 0); ASSERT_GE(iv2, 0); + EXPECT_NEAR(G[static_cast(iv1)], TWO_PI, 1e-10) + << "Gradient at v1 must be 2π (angle contribution = 0)"; + EXPECT_NEAR(G[static_cast(iv2)], TWO_PI, 1e-10) + << "Gradient at v2 must be 2π (angle contribution = 0)"; +} + // ════════════════════════════════════════════════════════════════════════════ // Gradient check: default right-isosceles triangle, vertex DOFs only // @@ -645,3 +736,87 @@ TEST(EuclideanFunctional, CyclicHessian_Analytic_MatchesBlockFD_Tetrahedron) max_sym = std::max(max_sym, std::abs(Ha.coeff(i, j) - Ha.coeff(j, i))); EXPECT_LT(max_sym, 1e-9) << "analytic cyclic Hessian is not symmetric"; } + +// ════════════════════════════════════════════════════════════════════════════ +// DOF assignment gauge-vertex overload (Finding-D, external-audit-2026-05-30) +// +// The single-argument assign_euclidean_vertex_dof_indices() overwrites ALL +// v_idx unconditionally — setting a pin *before* the call has no effect. +// The two-argument overload (gauge vertex) pins the requested vertex in a +// single pass. These tests verify: +// (a) single-argument: gauge must be set AFTER, not before. +// (b) two-argument: the gauge vertex gets v_idx = -1, others are sequential. +// (c) Newton converges correctly when the gauge overload is used. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(EuclideanDOFAssignment, SingleArg_PinBeforeHasNoEffect) +{ + // Pin first vertex before the call → the call overwrites it → not pinned. + auto mesh = make_tetrahedron(); + auto maps = setup_euclidean_maps(mesh); + auto first = *mesh.vertices().begin(); + + maps.v_idx[first] = -1; // set pin BEFORE — should have no effect + assign_euclidean_vertex_dof_indices(mesh, maps); + + EXPECT_GE(maps.v_idx[first], 0) + << "Pre-call pin was overwritten: v_idx[first] must be >= 0 after the call"; + EXPECT_EQ(euclidean_dimension(mesh, maps), + static_cast(mesh.number_of_vertices())) + << "All vertices should be free after single-arg assign"; +} + +TEST(EuclideanDOFAssignment, TwoArg_GaugeIsPinnedOthersAreSequential) +{ + auto mesh = make_tetrahedron(); + auto maps = setup_euclidean_maps(mesh); + auto first = *mesh.vertices().begin(); + + int n = assign_euclidean_vertex_dof_indices(mesh, maps, first); + + EXPECT_EQ(maps.v_idx[first], -1) + << "Gauge vertex must have v_idx = -1"; + EXPECT_EQ(n, static_cast(mesh.number_of_vertices()) - 1) + << "Returned DOF count must be num_vertices - 1"; + EXPECT_EQ(euclidean_dimension(mesh, maps), n) + << "euclidean_dimension must equal returned count"; + + // All non-gauge vertices must have distinct indices in [0, n). + std::vector seen; + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (v == first) continue; + EXPECT_GE(iv, 0); + EXPECT_LT(iv, n); + seen.push_back(iv); + } + std::sort(seen.begin(), seen.end()); + for (int i = 0; i < n; ++i) + EXPECT_EQ(seen[static_cast(i)], i) + << "DOF indices must be sequential 0..n-1"; +} + +TEST(EuclideanDOFAssignment, TwoArg_NewtonConvergesWithGaugeOverload) +{ + // End-to-end: use the gauge overload, then run Newton — confirms the + // pinned-vertex DOF layout is consistent with the solver. + auto mesh = make_tetrahedron(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + auto first = *mesh.vertices().begin(); + + int n = assign_euclidean_vertex_dof_indices(mesh, maps, first); + + // Natural-theta: set targets = actual angle sums at x=0 so x*=0 is the solution. + 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)]; + } + + auto res = newton_euclidean(mesh, x0, maps, 1e-10, 50); + EXPECT_TRUE(res.converged) + << "Newton did not converge with gauge-overload DOF assignment"; + EXPECT_LT(res.grad_inf_norm, 1e-9); +} diff --git a/code/tests/cgal/test_euclidean_hessian.cpp b/code/tests/cgal/test_euclidean_hessian.cpp index 94829b7..80305c6 100644 --- a/code/tests/cgal/test_euclidean_hessian.cpp +++ b/code/tests/cgal/test_euclidean_hessian.cpp @@ -214,3 +214,42 @@ TEST(EuclideanHessian, FDCheck_MixedPinnedVertices) EXPECT_TRUE(hessian_check_euclidean(mesh, x, maps)) << "FD Hessian check failed for mixed pinned/variable vertices"; } + +// ════════════════════════════════════════════════════════════════════════════ +// Edge-DOF guard (Finding-G, java-port-audit item 2) +// +// euclidean_hessian() (vertex-only cotangent Laplacian) must throw +// std::logic_error when any edge DOF is active. Without this guard the +// function would silently return a Hessian with zero rows/cols for the +// edge DOFs, causing SimplicialLDLT to fail in a hard-to-diagnose way. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(EuclideanHessian, EdgeDOFGuard_Throws) +{ + auto mesh = make_tetrahedron(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + assign_euclidean_all_dof_indices(mesh, maps); // assigns vertex + edge DOFs + + const int n = euclidean_dimension(mesh, maps); + std::vector x(static_cast(n), 0.0); + + EXPECT_THROW(euclidean_hessian(mesh, x, maps), std::logic_error) + << "euclidean_hessian must throw when edge DOFs are present"; +} + +TEST(EuclideanHessian, EdgeDOFGuard_VertexOnlyDoesNotThrow) +{ + // Vertex-only layout must NOT trigger the guard. + auto mesh = make_tetrahedron(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + auto gauge = *mesh.vertices().begin(); + assign_euclidean_vertex_dof_indices(mesh, maps, gauge); + + const int n = euclidean_dimension(mesh, maps); + std::vector x(static_cast(n), 0.0); + + EXPECT_NO_THROW(euclidean_hessian(mesh, x, maps)) + << "euclidean_hessian must not throw for vertex-only DOF layout"; +} diff --git a/code/tests/cgal/test_hyper_ideal_functional.cpp b/code/tests/cgal/test_hyper_ideal_functional.cpp index fabc536..85aa825 100644 --- a/code/tests/cgal/test_hyper_ideal_functional.cpp +++ b/code/tests/cgal/test_hyper_ideal_functional.cpp @@ -197,3 +197,94 @@ TEST(HyperIdealFunctional, GradientCheck_Fan6AllVariable) EXPECT_TRUE(gradient_check(mesh, x, maps)) << "Finite-difference gradient check failed on fan-6 mesh"; } + +// ════════════════════════════════════════════════════════════════════════════ +// Guard test: face_energy() must throw for 2+ ideal vertices in one face. +// +// This tests Finding-A from doc/reviewer/external-audit-2026-05-30.md. +// +// The Java reference (HyperIdealFunctional.java lines 222-231) silently applies +// the one-ideal-vertex volume formula to the first ideal vertex found, ignoring +// any additional ideal vertices in the same face. That is mathematically wrong +// for two-ideal / three-ideal faces. The C++ port detects this at runtime and +// throws std::logic_error instead of silently producing a wrong energy value. +// +// Tests cover: +// (a) Two ideal vertices in the same face (v1+v2 ideal, v3 hyper-ideal) +// (b) All three vertices ideal +// (c) Exactly one ideal vertex — must NOT throw (valid configuration) +// ════════════════════════════════════════════════════════════════════════════ + +TEST(HyperIdealFunctional, MultiIdealGuard_TwoIdealVertices_Throws) +{ + // Triangle mesh: 3 vertices, 3 edges, 1 face (open mesh, single face). + auto mesh = make_triangle(); + auto maps = setup_hyper_ideal_maps(mesh); + + // Assign edge DOFs to all three edges. + int eidx = 0; + for (auto e : mesh.edges()) maps.e_idx[e] = eidx++; + + // Pin v1 and v2 (ideal), make only v3 hyper-ideal. + auto vit = mesh.vertices().begin(); + Vertex_index v1 = *vit++; + Vertex_index v2 = *vit++; + // v3 remains pinned (default v_idx = -1, i.e. ideal too — see below). + + maps.v_idx[v1] = -1; // ideal + maps.v_idx[v2] = -1; // ideal + maps.v_idx[*vit] = 3; // hyper-ideal: DOF index 3 (after 3 edge DOFs) + + // DOF vector: [a_e0, a_e1, a_e2, b_v3] + std::vector x = {0.5, 0.5, 0.5, 1.0}; + + // evaluate_hyper_ideal calls face_energy() which must detect 2 ideal vertices + // and throw std::logic_error. + EXPECT_THROW( + evaluate_hyper_ideal(mesh, x, maps, /*energy=*/true, /*gradient=*/false), + std::logic_error) + << "Expected std::logic_error for face with two ideal vertices"; +} + +TEST(HyperIdealFunctional, MultiIdealGuard_AllThreeIdealVertices_Throws) +{ + auto mesh = make_triangle(); + auto maps = setup_hyper_ideal_maps(mesh); + + // Only edge DOFs — all vertices remain ideal (default v_idx = -1). + int eidx = 0; + for (auto e : mesh.edges()) maps.e_idx[e] = eidx++; + + // DOF vector: [a_e0, a_e1, a_e2] + std::vector x = {0.5, 0.5, 0.5}; + + EXPECT_THROW( + evaluate_hyper_ideal(mesh, x, maps, /*energy=*/true, /*gradient=*/false), + std::logic_error) + << "Expected std::logic_error for face with all three ideal vertices"; +} + +TEST(HyperIdealFunctional, MultiIdealGuard_ExactlyOneIdeal_DoesNotThrow) +{ + // Exactly one ideal vertex per face must NOT throw — it is the supported + // one-ideal-vertex configuration (Kolpakov-Mednykh formula). + auto mesh = make_triangle(); + auto maps = setup_hyper_ideal_maps(mesh); + + // All edges variable. + int eidx = 0; + for (auto e : mesh.edges()) maps.e_idx[e] = eidx++; + + // Pin only v1 (ideal); v2 and v3 are hyper-ideal. + auto vit = mesh.vertices().begin(); + maps.v_idx[*vit] = -1; ++vit; // v1: ideal + maps.v_idx[*vit] = 3; ++vit; // v2: hyper-ideal, DOF 3 + maps.v_idx[*vit] = 4; // v3: hyper-ideal, DOF 4 + + // DOF vector: [a_e0, a_e1, a_e2, b_v2, b_v3] + std::vector x = {0.5, 0.5, 0.5, 1.0, 1.0}; + + EXPECT_NO_THROW( + evaluate_hyper_ideal(mesh, x, maps, /*energy=*/true, /*gradient=*/false)) + << "Unexpected throw for valid one-ideal-vertex configuration"; +} diff --git a/code/tests/cgal/test_phase6.cpp b/code/tests/cgal/test_phase6.cpp index f77dfaf..f271207 100644 --- a/code/tests/cgal/test_phase6.cpp +++ b/code/tests/cgal/test_phase6.cpp @@ -137,6 +137,78 @@ TEST(GaussBonnet, ManuallySetAnalyticalTheta_PassesCheck) EXPECT_NO_THROW(check_gauss_bonnet(m, maps)); } +// ════════════════════════════════════════════════════════════════════════════ +// GaussBonnet — HyperIdeal API guard (Finding-B from external-audit-2026-05-30) +// +// gauss_bonnet_sum(mesh, HyperIdealMaps) and +// enforce_gauss_bonnet(mesh, HyperIdealMaps) are intentionally DELETED. +// +// Reason: the correct hyperbolic Gauss–Bonnet identity is +// Σ_v (2π − Θ_v) − Area(M) = 2π · χ(M) +// not the Euclidean form Σ(2π−Θ_v) = 2π·χ. For a regular (Θ_v=2π) genus-2 +// surface: Σ(2π−2π)=0 but 2π·χ=−4π, so the Euclidean check would always +// throw "deficit = 4π" for a perfectly valid HyperIdeal target. +// +// Compile-time enforcement: gauss_bonnet_sum / enforce_gauss_bonnet with +// HyperIdealMaps are = delete, so any accidental call is a compile error. +// The static_asserts below confirm this is wired correctly. +// +// The runtime test shows the discrepancy numerically: even for a regular +// tetrahedron where all Θ_v=2π (a valid all-hyper-ideal starting point), +// the "Euclidean sum" is 0 but the correct RHS for χ=2 is 4π — the +// Euclidean check would be off by 4π. +// ════════════════════════════════════════════════════════════════════════════ + +// Invocability check via SFINAE: tries to form the call expression in an +// unevaluated context; a deleted function causes substitution failure. +namespace { + template + auto try_gb_sum(int) -> decltype(gauss_bonnet_sum( + std::declval(), + std::declval()), std::true_type{}); + template + std::false_type try_gb_sum(...); +} // namespace + +static_assert(!decltype(try_gb_sum(0))::value, + "gauss_bonnet_sum must NOT be invocable with HyperIdealMaps"); +static_assert( decltype(try_gb_sum(0))::value, + "gauss_bonnet_sum must still be invocable with EuclideanMaps"); +static_assert( decltype(try_gb_sum(0))::value, + "gauss_bonnet_sum must still be invocable with SphericalMaps"); + +TEST(GaussBonnet, HyperIdeal_EuclideanSumDiscrepancy_DocumentsWhyCheckIsDeleted) +{ + // On a tetrahedron (χ=2) with all Θ_v = 2π: + // Euclidean sum Σ(2π−2π) = 0 + // but 2π·χ = 4π + // deficit (Euclidean formula) = 0 − 4π = −4π ← WRONG check for HyperIdeal + // + // The HyperIdeal identity is Σ(2π−Θ_v) − Area = 2π·χ. + // Area > 0 for any non-degenerate hyperbolic metric, so the real deficit + // would be much smaller. This test documents the mismatch numerically + // so any future re-introduction of the HyperIdeal overload is caught. + auto m = make_tetrahedron(); + auto hi_maps = setup_hyper_ideal_maps(m); + // All theta_v default to 2π (regular vertex target). + + // Access the raw property map directly (not via the deleted bundle overload) + // to compute the Euclidean-style sum — just for documentation purposes. + double euclid_sum = gauss_bonnet_sum(m, hi_maps.theta_v); // raw map: OK + double rhs = gauss_bonnet_rhs(m); // 2π·χ = 4π + + EXPECT_NEAR(euclid_sum, 0.0, 1e-12) // Σ(2π−2π) = 0 + << "Expected Euclidean sum = 0 for all-regular HyperIdeal targets"; + EXPECT_NEAR(rhs, 4.0 * M_PI, 1e-12) // 2π·χ(tetrahedron) = 4π + << "Expected RHS = 4π for tetrahedron (χ=2)"; + + // The Euclidean deficit would be 0 − 4π = −4π: completely wrong for HyperIdeal. + // If check_gauss_bonnet were called with HyperIdealMaps it would ALWAYS throw + // here, even though Θ_v=2π is a valid regular-vertex HyperIdeal target. + EXPECT_NEAR(euclid_sum - rhs, -4.0 * M_PI, 1e-10) + << "Euclidean deficit for HyperIdeal target = −4π: confirms the deleted API is correct"; +} + // ════════════════════════════════════════════════════════════════════════════ // CutGraph — tree-cotree algorithm // ════════════════════════════════════════════════════════════════════════════ diff --git a/code/tests/cgal/test_phase7.cpp b/code/tests/cgal/test_phase7.cpp index 06c9647..79d45a3 100644 --- a/code/tests/cgal/test_phase7.cpp +++ b/code/tests/cgal/test_phase7.cpp @@ -502,6 +502,109 @@ TEST(HolonomyEndToEnd, Torus8x8_TauMatchesRevolutionModulus) check_torus("torus_8x8.off", /*R=*/3.0, /*r=*/1.0, /*rel_tol=*/0.05); } +// ════════════════════════════════════════════════════════════════════════════ +// Finding-H (java-port-audit item 7, external-audit-2026-05-30): +// End-to-end torus with Re(τ) < 0 before normalizeModulus +// +// torus_skewed_4x4.off is a flat torus on a parallelogram lattice +// ω₁ = (4, 0) ω₂ = (−1, 4) +// The raw τ = ω₂/ω₁ = (−0.25 + i), Re < 0. +// After normalizeModulus the mirror fold gives τ = (0.25 + i), Re ≥ 0. +// +// This guards against a regression where compute_period_matrix uses +// reduce_to_fundamental_domain (old code, no mirror fold) instead of +// normalizeModulus (Java-faithful, finding 6 fix) — in that case the +// pipeline would silently report τ with Re < 0 instead of Re ≥ 0. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(HolonomyEndToEnd, SkewedTorus_ReTauNegativeBeforeNorm_FoldedToPositive) +{ + // ── Load the skewed flat torus ──────────────────────────────────────── + const std::string path = + std::string(CONFORMALLAB_DATA_DIR) + "/off/torus_skewed_4x4.off"; + ConformalMesh mesh = load_mesh(path); + ASSERT_GT(mesh.number_of_vertices(), 0u) << "Failed to load torus_skewed_4x4.off"; + ASSERT_EQ(conformallab::euler_characteristic(mesh), 0) + << "Mesh must be a torus (χ=0)"; + + // ── Run the full pipeline ───────────────────────────────────────────── + EuclideanMaps maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + + int idx = 0; + bool pinned = false; + for (auto v : mesh.vertices()) { + if (!pinned) { maps.v_idx[v] = -1; pinned = true; } + else maps.v_idx[v] = idx++; + } + enforce_gauss_bonnet(mesh, maps); + + std::vector x0(static_cast(idx), 0.0); + auto res = newton_euclidean(mesh, x0, maps); + ASSERT_TRUE(res.converged) << "Newton did not converge on skewed flat torus"; + + CutGraph cg = compute_cut_graph(mesh); + HolonomyData hol; + euclidean_layout(mesh, res.x, maps, &cg, &hol, /*normalise=*/false); + + ASSERT_EQ(hol.translations.size(), 2u) << "Expected exactly 2 holonomy generators"; + + // ── Raw τ (no normalization) must have Re < 0 ───────────────────────── + // This confirms the mesh geometry does produce a τ with negative real + // part, making the normalizeModulus step non-trivial. + PeriodData pd_raw = compute_period_matrix(hol, /*reduce=*/false); + EXPECT_LT(pd_raw.tau.real(), 0.0) + << "Raw τ must have Re < 0 for this skewed lattice" + << " (got Re = " << pd_raw.tau.real() << ")"; + + // ── Normalized τ must have Re ≥ 0 (normalizeModulus was applied) ───── + PeriodData pd = compute_period_matrix(hol, /*reduce=*/true); + EXPECT_GE(pd.tau.real(), -1e-10) + << "Normalized τ must have Re ≥ 0 (normalizeModulus mirror fold)" + << " (got Re = " << pd.tau.real() << ")"; + EXPECT_GT(pd.tau.imag(), 0.0) + << "τ must lie in the upper half-plane"; + EXPECT_GE(std::abs(pd.tau), 1.0 - 1e-9) + << "|τ| ≥ 1 (fundamental domain condition)"; + + // ── Additional fundamental-domain conditions ─────────────────────────── + // These are the normalizeModulus guarantees (Finding 6 / java-port-audit). + EXPECT_LE(pd.tau.real(), 0.5 + 1e-9) + << "normalizeModulus must produce Re(τ) ≤ ½"; + // The exact value depends on which generators tree-cotree finds; + // we do NOT assert a specific numeric value here (generator choice is + // an implementation detail of the tree-cotree algorithm, not of + // normalizeModulus). The assertions above are sufficient to confirm + // that the mirror fold was applied. +} + +// ════════════════════════════════════════════════════════════════════════════ +// Finding-H synthetic sanity: compute_period_matrix with explicit Re(τ)<0 +// holonomy verifies the mirror fold numerically (no mesh, no tree-cotree). +// ════════════════════════════════════════════════════════════════════════════ + +TEST(HolonomyEndToEnd, SyntheticHolonomy_NegativeReTau_NormalizedToPositive) +{ + // Lattice: ω₁=(4,0), ω₂=(-1,4) → τ_raw = (-1+4i)/4 = -0.25+i + // normalizeModulus: Re=-0.25 < 0 → mirror: τ = -conj(τ) = +0.25+i + HolonomyData hol; + hol.translations = { + Eigen::Vector2d(4.0, 0.0), + Eigen::Vector2d(-1.0, 4.0) + }; + + PeriodData pd_raw = compute_period_matrix(hol, /*reduce=*/false); + EXPECT_NEAR(pd_raw.tau.real(), -0.25, 1e-10) << "Raw Re(τ) must be -0.25"; + EXPECT_NEAR(pd_raw.tau.imag(), 1.0, 1e-10) << "Raw Im(τ) must be 1.0"; + + PeriodData pd = compute_period_matrix(hol, /*reduce=*/true); + EXPECT_GE(pd.tau.real(), 0.0 - 1e-9) << "Normalized Re(τ) ≥ 0"; + EXPECT_LE(pd.tau.real(), 0.5 + 1e-9) << "Normalized Re(τ) ≤ ½"; + EXPECT_NEAR(pd.tau.real(), 0.25, 1e-9) << "Mirror fold: Re = -0.25 → +0.25"; + EXPECT_NEAR(pd.tau.imag(), 1.0, 1e-9) << "Im(τ) preserved by mirror fold"; + EXPECT_GE(std::abs(pd.tau), 1.0 - 1e-9) << "|τ| ≥ 1"; +} + // ════════════════════════════════════════════════════════════════════════════ // FundamentalDomain — genus-1 parallelogram // ════════════════════════════════════════════════════════════════════════════ diff --git a/code/tests/cgal/test_spherical_functional.cpp b/code/tests/cgal/test_spherical_functional.cpp index a8d745f..e82a6f4 100644 --- a/code/tests/cgal/test_spherical_functional.cpp +++ b/code/tests/cgal/test_spherical_functional.cpp @@ -652,3 +652,52 @@ TEST(SphericalGoldenJava, FullMeshEdgeDofGradient_Tetrahedron) EXPECT_NEAR(G[static_cast(maps.e_idx[eAB])], -0.35189517043413690, 1e-12); EXPECT_NEAR(G[static_cast(maps.e_idx[eCD])], -0.44101986058895950, 1e-12); } + +// ════════════════════════════════════════════════════════════════════════════ +// Degenerate spherical triangle — limiting angles (Finding-F, java-port-audit item 1) +// +// spherical_angles() must return the limiting angles (π opposite the +// over-long edge, 0/0 elsewhere) when the spherical triangle inequality +// is violated, with valid = false. This mirrors the Euclidean behaviour +// and is required for the convex C¹ BPS extension. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalFunctional, DegenerateTriangle_LimitingAngles_S12TooLong) +{ + // s12 > s23 + s31: s12 = 2.5, s23 = s31 = 0.5 (all < π so valid arc lengths) + // s23 < 0 → actually use s-based check + // Easier: use s12 = π − ε (nearly degenerate hemisphere edge) + // and very short s23, s31 so s12 > s23 + s31. + const double s12 = 2.0, s23 = 0.4, s31 = 0.4; // s12 > s23+s31 = 0.8 + auto fa = spherical_angles(s12, s23, s31); + EXPECT_FALSE(fa.valid); + // s12 is the edge opposite v3 → α3 = π + EXPECT_NEAR(fa.alpha3, PI, 1e-12) + << "corner opposite over-long s12 must be π"; + EXPECT_NEAR(fa.alpha1, 0.0, 1e-12); + EXPECT_NEAR(fa.alpha2, 0.0, 1e-12); +} + +TEST(SphericalFunctional, DegenerateTriangle_LimitingAngles_S23TooLong) +{ + // s23 > s12 + s31 → α1 = π + const double s12 = 0.4, s23 = 2.0, s31 = 0.4; + auto fa = spherical_angles(s12, s23, s31); + EXPECT_FALSE(fa.valid); + EXPECT_NEAR(fa.alpha1, PI, 1e-12) + << "corner opposite over-long s23 must be π"; + EXPECT_NEAR(fa.alpha2, 0.0, 1e-12); + EXPECT_NEAR(fa.alpha3, 0.0, 1e-12); +} + +TEST(SphericalFunctional, DegenerateTriangle_LimitingAngles_S31TooLong) +{ + // s31 > s12 + s23 → α2 = π + const double s12 = 0.4, s23 = 0.4, s31 = 2.0; + auto fa = spherical_angles(s12, s23, s31); + EXPECT_FALSE(fa.valid); + EXPECT_NEAR(fa.alpha2, PI, 1e-12) + << "corner opposite over-long s31 must be π"; + EXPECT_NEAR(fa.alpha1, 0.0, 1e-12); + EXPECT_NEAR(fa.alpha3, 0.0, 1e-12); +} diff --git a/doc/reviewer/external-audit-2026-05-30.md b/doc/reviewer/external-audit-2026-05-30.md new file mode 100644 index 0000000..36dc741 --- /dev/null +++ b/doc/reviewer/external-audit-2026-05-30.md @@ -0,0 +1,746 @@ +# External Code Audit — ConformalLabpp v0.10.0 + +**Date:** 2026-05-30 +**Auditor:** External reviewer (Claude Sonnet 4.6) +**Branch:** `review/external-audit-2026-05-30` +**Base:** `docs/fix-test-count-post-merge` (HEAD at audit time) + +This document is self-contained. A new session can pick up any finding +below and act on it without prior context. Each finding includes: +- exact file path + line numbers (verified by direct file read) +- a minimal reproduction of the problematic code +- the correct fix or recommended action +- acceptance criteria for "done" + +Status legend: 🔴 Bug · 🟡 API/Doc error · 🟠 Test gap · 🔵 Architectural risk + +--- + +## How to read this document in a new session + +```bash +# 1. Check out the audit branch +git checkout review/external-audit-2026-05-30 + +# 2. Build the CGAL test suite (needed for verification) +cmake -S code -B build-cgal -DWITH_CGAL_TESTS=ON +cmake --build build-cgal --target conformallab_cgal_tests -j$(nproc) + +# 3. Run the full suite before making any change +ctest --test-dir build-cgal -R '^cgal\.' --output-on-failure +# Expected: 246 passed, 0 failed + +# 4. Pick a finding below, apply the fix, re-run ctest, then commit. +``` + +The Java reference implementation lives at: +``` +/Users/tarikmoussa/Desktop/conformallab/src/de/varylab/discreteconformal/ +``` +Consult it for any port-faithfulness question. + +--- + +## FINDING-A — 🔴 CRITICAL BUG: `face_energy()` silently wrong for mixed ideal/hyper-ideal configurations + +### Location +`code/include/hyper_ideal_functional.hpp` lines 319–344 + +### Problem + +The `face_energy()` function handles only two cases: "all three vertices hyper-ideal" +and "exactly one vertex ideal (pinned)". When two or all three vertices are ideal +(`v?b = m.v_idx[v?] < 0`), the cascade falls through to a branch that calls the +**one-ideal-vertex** volume formula — which is mathematically wrong for two or +three ideal vertices. + +```cpp +// CURRENT (broken for >= 2 ideal vertices): +static double face_energy(const FaceAngles& fa) +{ + ... + if (fa.v1b && fa.v2b && fa.v3b) { // all hyper-ideal ✓ + V = calculateTetrahedronVolume( + fa.beta1, fa.beta2, fa.beta3, + fa.alpha23, fa.alpha31, fa.alpha12); + } else if (!fa.v1b) { // BUG: enters even when !v1b && !v2b + V = calculateTetrahedronVolumeWithIdealVertexAtGamma( + fa.beta1, fa.alpha31, fa.alpha12, + fa.alpha23, fa.beta2, fa.beta3); + } else if (!fa.v2b) { // BUG: enters even when !v2b && !v3b + V = calculateTetrahedronVolumeWithIdealVertexAtGamma( + fa.beta2, fa.alpha12, fa.alpha23, + fa.alpha31, fa.beta3, fa.beta1); + } else { // !v3b only // only correct for exactly 1 ideal + V = calculateTetrahedronVolumeWithIdealVertexAtGamma( + fa.beta3, fa.alpha23, fa.alpha31, + fa.alpha12, fa.beta1, fa.beta2); + } + return aa + bb + 2.0 * V; +} +``` + +**The same structural error exists in `face_angles_from_local_dofs()`** at lines 192–215: +```cpp +if (l12 > l23 + l31) { + o.beta1 = 0.0; o.beta2 = 0.0; o.beta3 = PI; + o.alpha12 = PI; o.alpha23 = 0.0; o.alpha31 = 0.0; +} else if (l23 > l12 + l31) { ... } +else if (l31 > l12 + l23) { ... } +else { /* normal */ } +``` +(this part is actually correct — no bug here; the degenerate case structure is +standard. Keeping the note for completeness.) + +### Trigger condition + +Only triggered when the `HyperIdealMaps` has **some** vertices pinned (`v_idx[v] = -1`) +and **some** variable. The default workflow `assign_all_dof_indices(mesh, maps)` makes +ALL vertices variable (v?b = true everywhere), which always takes the first branch — +safe. The bug only fires for mixed ideal/hyper-ideal configurations. + +### Fix + +The missing cases are the fully-ideal (all three ideal) and two-ideal-vertex cases. +The Java reference `HyperIdealFunctional.java` must be consulted to find the correct +volume formula for two ideal vertices. The typical fix structure is: + +```cpp +// PROPOSED FIX — verify against Java reference before applying: +if (fa.v1b && fa.v2b && fa.v3b) { + // all hyper-ideal + V = calculateTetrahedronVolume( + fa.beta1, fa.beta2, fa.beta3, + fa.alpha23, fa.alpha31, fa.alpha12); +} else if (fa.v1b && fa.v2b && !fa.v3b) { + // ideal at v3 only + V = calculateTetrahedronVolumeWithIdealVertexAtGamma( + fa.beta3, fa.alpha23, fa.alpha31, + fa.alpha12, fa.beta1, fa.beta2); +} else if (fa.v1b && !fa.v2b && fa.v3b) { + // ideal at v2 only + V = calculateTetrahedronVolumeWithIdealVertexAtGamma( + fa.beta2, fa.alpha12, fa.alpha23, + fa.alpha31, fa.beta3, fa.beta1); +} else if (!fa.v1b && fa.v2b && fa.v3b) { + // ideal at v1 only + V = calculateTetrahedronVolumeWithIdealVertexAtGamma( + fa.beta1, fa.alpha31, fa.alpha12, + fa.alpha23, fa.beta2, fa.beta3); +} else if (!fa.v1b && !fa.v2b && fa.v3b) { + // two ideal: v1, v2 — FORMULA NEEDED (see Java reference) + V = 0.0; // TODO: implement two-ideal-vertex formula +} else if (!fa.v1b && fa.v2b && !fa.v3b) { + // two ideal: v1, v3 — FORMULA NEEDED + V = 0.0; // TODO +} else if (fa.v1b && !fa.v2b && !fa.v3b) { + // two ideal: v2, v3 — FORMULA NEEDED + V = 0.0; // TODO +} else { + // all three ideal: volume = sum of Lobachevsky only + V = 0.0; // TODO: verify correct formula +} +``` + +Check Java: `HyperIdealFunctional.java` → `triangleEnergyAndAlphas()` for all cases. + +### Resolution (2026-05-30) + +**Root cause clarified:** The Java reference (`HyperIdealFunctional.java` lines 219–233) +has the **identical** cascade structure — it silently applies the one-ideal-vertex formula +to the first ideal vertex found, ignoring additional ideal vertices in the same face. +The C++ was a faithful port; it did not introduce a new divergence. + +The fix therefore does NOT change the calculation (which would diverge from Java). +Instead it: + +1. **Adds `` include** to `hyper_ideal_functional.hpp`. +2. **Adds an `ideal_count` guard** at the top of `face_energy()` that counts ideal + vertices (`!v?b`) and throws `std::logic_error` for `ideal_count >= 2`, replacing + the silent wrong result with a clear diagnostic message that names the limitation + and points to the audit document. +3. **Adds three new GTest cases** to `test_hyper_ideal_functional.cpp`: + - `MultiIdealGuard_TwoIdealVertices_Throws` — two ideal vertices → throw ✅ + - `MultiIdealGuard_AllThreeIdealVertices_Throws` — all ideal → throw ✅ + - `MultiIdealGuard_ExactlyOneIdeal_DoesNotThrow` — one ideal → no throw ✅ +4. **Adds a doc comment** above `face_energy()` explaining the supported configurations + (0-ideal / 1-ideal), the Java reference limitation, and why multi-ideal faces are + not implemented (requires new research beyond the Java reference). + +**Test result:** 262/262 CGAL tests pass (was 246 before + 3 new guard tests + 13 from +prior suite growth). No regressions. + +**Further work:** Implementing the *correct* volume formulas for 2-ideal and 3-ideal +faces is a separate research item (not a port — Java does not have them either). +It is now tracked in `doc/roadmap/research-track.md` under +"Hyper-ideal volume formulas for 2- and 3-ideal-vertex faces (Phase 9b+)". +The `throw` remains the correct safe behaviour until that research item is resolved. + +### Acceptance criteria +- [x] All 8 cases (2³ combinations) are either correct or throw explicitly +- [x] Three new tests cover the three guard cases +- [x] 262 CGAL tests pass, 0 failed + +--- + +## FINDING-B — 🟡 API CONCEPTUAL ERROR: Gauss–Bonnet check silently wrong for HyperIdeal + +### Location +`code/include/gauss_bonnet.hpp` lines 87–88, 128–134 + +### Problem + +The function `gauss_bonnet_sum()` is overloaded for all five map types including +`HyperIdealMaps`, and `check_gauss_bonnet()` is templated so it accepts any Maps: + +```cpp +// gauss_bonnet.hpp:87 +inline double gauss_bonnet_sum(const ConformalMesh& m, const HyperIdealMaps& mp) + { return gauss_bonnet_sum(m, mp.theta_v); } + +// gauss_bonnet.hpp:128-134 +template +inline void check_gauss_bonnet(const ConformalMesh& mesh, const Maps& maps, + double tol = 1e-8) +{ + check_gauss_bonnet(mesh, gauss_bonnet_sum(mesh, maps), tol); // always throws for HyperIdeal! +} +``` + +The Gauss–Bonnet identity checked here is: + +``` +Σ_v (2π − Θ_v) = 2π · χ(M) ← Euclidean / flat case only +``` + +For a **hyperbolic** metric (which is what HyperIdeal computes), the correct identity is: + +``` +Σ_v (2π − Θ_v) − Area(M) = 2π · χ(M) +``` + +For a regular (no cone singularities) genus-2 surface with Θ_v = 2π for all v: +- LHS: Σ(2π − 2π) = 0 +- RHS expected by check: 2π·χ = 2π·(2−2·2) = −4π + +So `|0 − (−4π)| = 4π ≫ tol` → `check_gauss_bonnet` always throws for valid +hyperbolic configurations. The overload exists but calling it produces a wrong result. + +### Fix options + +**Option A (recommended):** Delete the `HyperIdealMaps` overload of `gauss_bonnet_sum` +and add a compile-time or doc-level guard that prevents `check_gauss_bonnet` from +being instantiated with `HyperIdealMaps`. Add a comment explaining why. + +```cpp +// DELETE this overload: +// inline double gauss_bonnet_sum(const ConformalMesh& m, const HyperIdealMaps& mp) + +// ADD to check_gauss_bonnet doc: +// NOTE: Do NOT call with HyperIdealMaps — the hyperbolic Gauss–Bonnet identity +// includes an Area term absent from this check. For HyperIdeal configurations +// there is no simple global angle constraint analogous to the Euclidean case. +``` + +**Option B:** Keep the overload but make it compute the correct hyperbolic quantity +`Σ(2π−Θ_v) − Area(M)` and add a separate `check_hyperbolic_gauss_bonnet` that +compares against `2π·χ`. This requires computing the area from the HyperIdeal +metric, which is non-trivial. + +### Resolution (2026-05-31) + +1. **`gauss_bonnet_sum(mesh, HyperIdealMaps)` deleted** — replaced with + `= delete` overload and a multi-line comment explaining the Area-term + discrepancy. Attempting to call this is now a compile error. +2. **`enforce_gauss_bonnet(mesh, HyperIdealMaps&)` deleted** — explicit + `= delete` overload prevents the generic template from being silently + instantiated with HyperIdealMaps. +3. **Header comment block rewritten** — now has a clear box explaining + which geometries are supported (Euclidean/Spherical) and which are not + (HyperIdeal), with the correct hyperbolic Gauss–Bonnet identity shown. +4. **One new GTest** in `test_phase6.cpp`: + - `HyperIdeal_EuclideanSumDiscrepancy_DocumentsWhyCheckIsDeleted` — + verifies numerically that for a regular (Θ_v=2π) tetrahedron the + Euclidean sum = 0 while 2π·χ = 4π, i.e. the Euclidean check would + produce deficit = −4π for a valid HyperIdeal target. +5. **Three compile-time `static_assert`s** in the test file (SFINAE-based) + confirm that `gauss_bonnet_sum` is NOT invocable with HyperIdealMaps + but IS invocable with EuclideanMaps and SphericalMaps. + +**Test result:** 263/263 CGAL tests pass. No regressions. + +### Acceptance criteria +- [x] Calling `check_gauss_bonnet(mesh, hyper_ideal_maps)` → compile error +- [x] Calling `enforce_gauss_bonnet(mesh, hyper_ideal_maps)` → compile error +- [x] Comments explain why Euclidean G-B does not apply to HyperIdeal +- [x] 263 CGAL tests pass, 0 failed + +--- + +## FINDING-C — 🟡 DOCUMENTATION BUG: Cotangent formula in header comment is wrong + +### Location +`code/include/euclidean_hessian.hpp` line 26–27 (box comment) and line 57–58 +(comment above `euclidean_cot_weights`) + +### Problem + +The box comment at the top of the file states: + +``` +│ cot_k = (t_adj1·l123 − t_adj2·t_opp) / denom2 │ +│ = cotangent of the angle αk at vertex k +``` + +And the comment above `euclidean_cot_weights`: + +```cpp +// cot_k = (t_adj·l123 − t_opp·t_other) / (8·Area) +``` + +Both are **mathematically wrong**. The correct formula (verified numerically) is: + +``` +cot_k = (t_opp · l123 − t_adj1 · t_adj2) / (8·Area) +``` + +**Numerical proof** (3-4-5 right triangle, l23=3, l31=4, l12=5): +- t12=2, t23=6, t31=4, l123=12, 8·Area=48 +- cot(α₁) = 4/3 (angle at v₁ opposite l23=3) +- Code result: `(t23·l123 − t31·t12)/48 = (6·12 − 4·2)/48 = 64/48 = 4/3` ✓ +- Header formula: `(t_adj1·l123 − t_adj2·t_opp)/48 = (t12·l123 − t31·t23)/48 = (2·12 − 4·6)/48 = 0` ✗ + +The **implementation** in `euclidean_cot_weights` (lines 91–96) is correct. +Only the documentation is wrong. + +### Fix + +```cpp +// REPLACE both comment instances with the correct formula: + +// cot_k = (t_opp · l123 − t_adj1 · t_adj2) / (8·Area) +// +// where for vertex k: +// t_opp = t-value of the edge OPPOSITE to k (= 2(s − l_opp)) +// t_adj1, t_adj2 = t-values of the two edges ADJACENT to k +// l123 = l12 + l23 + l31 (perimeter) +// 8·Area = denom2 +``` + +### Acceptance criteria +- [ ] Both comment blocks corrected (box comment + function-level comment) +- [ ] A one-line comment in the function body confirms which t-value is `t_opp` + for each return value, e.g.: `// cot1: t_opp=t23 (opposite v1)` +- [ ] No code changes — only documentation + +--- + +## FINDING-D — 🟡 DOCUMENTATION BUG: DOF-assignment functions claim "pin before" works, but it doesn't + +### Location +- `code/include/euclidean_functional.hpp` lines 97–107 (`assign_euclidean_vertex_dof_indices`) +- `code/include/spherical_functional.hpp` lines 90–96 (`assign_vertex_dof_indices`) +- `code/include/inversive_distance_functional.hpp` lines 123–139 + (`assign_inversive_distance_vertex_dof_indices`) + +### Problem + +All three functions iterate unconditionally over every vertex and assign a sequential +index, overwriting any previously set `-1` pin: + +```cpp +// euclidean_functional.hpp:104-107 +inline int assign_euclidean_vertex_dof_indices(ConformalMesh& mesh, EuclideanMaps& m) +{ + int idx = 0; + for (auto v : mesh.vertices()) m.v_idx[v] = idx++; // overwrites ALL + return idx; +} +``` + +The doc comment above this function says: + +``` +/// **Note:** does NOT pin a gauge vertex. For closed meshes the caller +/// must set one `m.v_idx[v] = -1` either before or after this call +``` + +"Before" is **wrong** — the loop overwrites any pre-set pin. + +The inversive-distance version is worse: + +```cpp +// inversive_distance_functional.hpp:126-130 +/// 1. set one `m.v_idx[v] = -1` *before* calling this function (then +/// the call is a no-op for that vertex) — OR — +``` + +Explicitly claims the call is "a no-op" for a pre-pinned vertex, which is false. + +### Consequence + +A caller who pins `m.v_idx[first_vertex] = -1` and then calls +`assign_euclidean_vertex_dof_indices()` will get a **fully-free system** with no +gauge fix. On a closed mesh the Hessian is singular (gauge mode). `SimplicialLDLT` +silently falls back to `SparseQR`, which finds a minimum-norm solution — +no error is reported, but the result is not the intended pinned solution. + +### Fix + +**Option A (minimal):** Fix the doc comments only. Remove "before" as an option; +only "after" works. + +```cpp +/// NOTE: this function assigns indices to ALL vertices unconditionally. +/// To pin a gauge vertex, set `m.v_idx[v] = -1` AFTER calling this function. +``` + +**Option B (API improvement):** Add an overload that accepts a gauge vertex: + +```cpp +inline int assign_euclidean_vertex_dof_indices( + ConformalMesh& mesh, EuclideanMaps& m, Vertex_index gauge) +{ + int idx = 0; + for (auto v : mesh.vertices()) + m.v_idx[v] = (v == gauge) ? -1 : idx++; + return idx; +} +``` + +### Acceptance criteria +- [ ] Doc comments corrected in all three files to say "AFTER" only +- [ ] Optionally: overload with explicit gauge vertex added +- [ ] No code changes required for correctness (existing callers set pin after, + which already works) + +--- + +## FINDING-E — 🟡 INCONSISTENCY: `gradient_check_cp_euclidean` uses absolute error, all others use relative + +### Location +`code/include/cp_euclidean_functional.hpp` lines 338–349 (`gradient_check_cp_euclidean`) +and lines 373–387 (`hessian_check_cp_euclidean`) + +### Problem + +Every gradient check in the library normalises the error by the gradient magnitude: + +```cpp +// euclidean_functional.hpp:343 — RELATIVE error +double scale = std::max(1.0, std::abs(G[si])); +if (err / scale > tol) ok = false; +``` + +But `gradient_check_cp_euclidean` uses **absolute** error: + +```cpp +// cp_euclidean_functional.hpp:345 — ABSOLUTE error (inconsistent) +if (std::abs(G[i] - fd) > tol) { + std::cerr << "[cp-euclidean] FD gradient mismatch ..."; + return false; +} +``` + +Same issue in `hessian_check_cp_euclidean` line 378: + +```cpp +if (std::abs(an - fd) > tol) { // absolute, not relative +``` + +The default `tol = 1e-6` is acceptable for unit-scale problems but will produce +false failures if the gradient values grow large (e.g., many faces, large ρ). + +### Fix + +```cpp +// gradient_check_cp_euclidean — replace the comparison: +double err = std::abs(G[i] - fd); +double scale = std::max(1.0, std::abs(G[i])); +if (err / scale > tol) { + std::cerr << "[cp-euclidean] FD gradient mismatch at DOF " << i + << ": analytic=" << G[i] << " FD=" << fd + << " rel-err=" << (err / scale) << "\n"; + return false; +} + +// hessian_check_cp_euclidean — same pattern: +double err = std::abs(an - fd); +double scale = std::max(1.0, std::abs(an)); +if (err / scale > tol) { ... } +``` + +### Acceptance criteria +- [ ] Both `gradient_check_cp_euclidean` and `hessian_check_cp_euclidean` use + relative error (normalised by `max(1.0, |analytic|)`) +- [ ] The existing CP-Euclidean gradient check tests still pass + +--- + +## FINDING-F — 🟠 TEST GAP: Degenerate-triangle gradient (Finding 1 from java-port-audit.md) + +### Location +`code/tests/cgal/test_euclidean_functional.cpp` and +`code/tests/cgal/test_spherical_functional.cpp` + +### Problem (from java-port-audit.md item 1, still open) + +The degenerate-triangle fix (Finding 1 in java-port-audit.md) made `euclidean_angles()` +and `spherical_angles()` return the limiting angles (π opposite the over-long edge, +0/0 for the others) instead of `{0,0,0}`. This change is untested: no test builds a +triangle with one edge longer than the sum of the others and asserts the π corner. + +### Required test + +```cpp +// Add to test_euclidean_functional.cpp: +TEST(EuclideanGeometry, DegenerateTriangle_LimitingAngles) { + // l12 = 10, l23 = 1, l31 = 1 → l23+l31=2 < l12=10 → degenerate + // Expected: alpha3 = π (vertex v3 opposite l12), alpha1=alpha2=0 + auto fa = euclidean_angles_from_lengths(10.0, 1.0, 1.0); + EXPECT_FALSE(fa.valid); + EXPECT_NEAR(fa.alpha3, conformallab::PI, 1e-12); + EXPECT_NEAR(fa.alpha1, 0.0, 1e-12); + EXPECT_NEAR(fa.alpha2, 0.0, 1e-12); +} + +// Similarly for spherical_angles in test_spherical_functional.cpp +``` + +### Acceptance criteria +- [ ] Test added for Euclidean degenerate triangle → `alpha3 = π` +- [ ] Test added for Spherical degenerate triangle → same +- [ ] Both tests pass + +--- + +## FINDING-G — 🟠 TEST GAP: `euclidean_hessian` edge-DOF guard must throw (Finding 2 from java-port-audit.md) + +### Location +`code/tests/cgal/test_euclidean_hessian.cpp` (or new file) + +### Problem (from java-port-audit.md item 2, still open) + +Finding 2 added a `throw std::logic_error` guard to `euclidean_hessian()` when +edge DOFs are present. No test asserts this throw fires. + +### Required test + +```cpp +TEST(EuclideanHessian, EdgeDOFGuard_Throws) { + auto mesh = make_tetrahedron(); + auto maps = setup_euclidean_maps(mesh); + assign_euclidean_all_dof_indices(mesh, maps); // assigns edge DOFs + std::vector x(euclidean_dimension(mesh, maps), 0.0); + compute_euclidean_lambda0_from_mesh(mesh, maps); + + EXPECT_THROW( + euclidean_hessian(mesh, x, maps), + std::logic_error + ); +} +``` + +### Acceptance criteria +- [ ] Test added and passes (throw confirmed) +- [ ] Test is in the cgal suite (build with `-DWITH_CGAL_TESTS=ON`) + +--- + +## FINDING-H — 🟠 TEST GAP: Missing end-to-end torus with Re(τ) < 0 before reduction (java-port-audit.md item 7) + +### Location +`code/tests/cgal/test_pipeline.cpp` or `code/tests/cgal/test_phase7.cpp` + +### Problem (from java-port-audit.md item 7, partly open) + +`compute_period_matrix` now calls `normalizeModulus` which maps τ into the +half-strip `0 ≤ Re(τ) ≤ ½`. This is tested at the function level by the Java oracle +(`NormalizeModulus_GoldenJava`). But there is **no end-to-end test** where the +initial τ has `Re(τ) < 0` and the pipeline is verified to fold it into `Re(τ) ≥ 0`. + +If someone reverts `compute_period_matrix` to call `reduce_to_fundamental_domain` +instead of `normalizeModulus`, the function-level test still passes but the pipeline +output would silently diverge from the Java oracle. + +### Required test + +A torus whose geometry produces `Re(τ) < 0` before SL(2,ℤ) reduction. One approach: +construct a torus with an asymmetric lattice (e.g., parallelogram with obtuse angle on +the left side) where the natural τ has negative real part. + +```cpp +TEST(PeriodMatrix, EndToEnd_NegativeReTau_FoldedToPositive) { + // Build or load a torus mesh whose natural τ has Re < 0. + // Run the full pipeline: newton_euclidean → layout → compute_period_matrix. + // Assert: Re(result.tau) >= 0 (normalizeModulus was applied) + // Assert: Im(result.tau) >= 0 + // Assert: |result.tau| >= 1 +} +``` + +### Acceptance criteria +- [ ] End-to-end test added that exercises a mesh with `Re(τ) < 0` pre-reduction +- [ ] Test asserts `Re(τ) ≥ 0` after `compute_period_matrix` +- [ ] 246 CGAL tests still pass with new test included + +--- + +## FINDING-I — 🔵 ARCHITECTURAL RISK: 246/272 CGAL tests are not gated in CI + +### Location +`.gitea/workflows/cpp-tests.yml` line (see `if: false` block), CLAUDE.md lines 267–274 + +### Problem + +The CI pipeline has three jobs: +1. `test-fast` — 26 pure-math tests (no CGAL), **active** +2. `test-cgal` — 246 CGAL tests, **disabled** (`if: false` since 2026-05-26) +3. `quality-gates` — structural linting, **active** + +This means the Newton solvers, layout, holonomy, period matrix, cut graph, CGAL +public API, and all five DCE models are **not regression-tested on any push**. + +Root cause: the CGAL build OOMs on the 1.6 GB ARM64 Raspberry Pi runner at `-j` +parallel compilation. + +### Recommended actions (in priority order) + +1. **Immediate (low risk):** Run `test-cgal` with `-j1` (serial build) to avoid OOM. + The wall time increases but correctness is not compromised. + ```yaml + cmake --build build-cgal --target conformallab_cgal_tests -j1 + ``` + +2. **Short term:** Split the CGAL test binary into subsets so a failing compilation + is localised. Add a minimal subset (e.g. Newton + gradient checks only) as a + new CI job that runs on every push. + +3. **Medium term:** Add a GitHub Actions job (arm64 runner, 7 GB RAM) to mirror CI. + Self-hosted Raspberry Pi is not suitable for a library targeting CGAL submission. + +4. **Document the gap explicitly** in CHANGELOG.md and doc/release-policy.md: + "v0.10.0 ships with CGAL CI disabled — run `ctest -R '^cgal\.'` locally before + tagging any release." + +### Acceptance criteria +- [ ] At least one of the above options implemented +- [ ] The CGAL test suite runs in CI on every PR (not just locally) +- [ ] CHANGELOG.md documents the current CI limitation + +--- + +## MINOR FINDINGS (quick fixes, no architectural impact) + +### MINOR-1 — `spherical_gauge_shift` misleading comment [`spherical_functional.hpp:404`] + +```cpp +// CURRENT (wrong): +// f is strictly monotone decreasing (second derivative < 0) for a convex functional + +// FIX: +// f is strictly monotone decreasing (sum of vertex angles increases with scale) +// for any conservative gradient, including the concave spherical energy. +``` + +### MINOR-2 — `spherical_gauge_shift` uses forward FD, should use central [`spherical_functional.hpp:470–471`] + +```cpp +// CURRENT (forward difference, O(ε)): +double ftp = sum_Gv(t + fd_eps); +double dft = (ftp - ft) / fd_eps; + +// FIX (central difference, O(ε²), same cost in this context): +double dft = (sum_Gv(t + fd_eps) - sum_Gv(t - fd_eps)) / (2.0 * fd_eps); +``` + +### MINOR-3 — Gauss-Legendre constants duplicated in three files + +`gl_s[10]` and `gl_w[10]` are identically copy-pasted in: +- `code/include/euclidean_functional.hpp` lines 257–270 +- `code/include/spherical_functional.hpp` lines 304–317 +- `code/include/inversive_distance_functional.hpp` lines 312–325 + +**Fix:** Extract to a new header `gauss_legendre.hpp`: +```cpp +// code/include/gauss_legendre.hpp (new file) +namespace conformallab::detail { + inline const double* gl10_nodes() { static const double s[10] = {...}; return s; } + inline const double* gl10_weights() { static const double w[10] = {...}; return w; } +} +``` + +### MINOR-4 — `hidx()` function duplicated in four files + +Same one-liner in `euclidean_functional.hpp`, `spherical_functional.hpp`, +`hyper_ideal_functional.hpp`, `inversive_distance_functional.hpp`. Move to +`conformal_mesh.hpp` as: +```cpp +// code/include/conformal_mesh.hpp (add): +inline std::size_t halfedge_idx(Halfedge_index h) noexcept { + return static_cast(static_cast(h)); +} +``` + +### MINOR-5 — `inits()` in `clausen.hpp` has unintuitive off-by-one semantics + +The function returns `n` after decrement, which is one less than the last-checked +index. Add a comment: + +```cpp +// Returns the index one below the last term needed to reach the requested +// accuracy — callers use this as the `n` argument to csevl(), which then +// evaluates terms [0..n-1]. Matches Java Clausen.inits() semantics exactly. +``` + +--- + +## Summary table + +| ID | File | Lines | Type | Severity | Status | +|----|------|-------|------|----------|--------| +| A | `hyper_ideal_functional.hpp` | 319–344 | Bug | Critical (mixed config) | ✅ Fixed 2026-05-30 | +| B | `gauss_bonnet.hpp` | 87–88, 128–134 | API error | Medium | ✅ Fixed 2026-05-31 | +| C | `euclidean_hessian.hpp` | 26–27, 57–58 | Doc error | Medium | ✅ Fixed 2026-05-31 | +| D | `euclidean_functional.hpp` + 2 others | 97–107 | Doc error | Medium | ✅ Fixed 2026-05-31 | +| E | `cp_euclidean_functional.hpp` | 338–349, 373–387 | Inconsistency | Medium | ✅ Fixed 2026-05-31 | +| F | test files | — | Test gap | Medium | ✅ Fixed 2026-05-31 | +| G | test files | — | Test gap | Medium | ✅ Fixed 2026-05-31 | +| H | test files | — | Test gap | Medium | ✅ Fixed 2026-05-31 | +| I | CI workflow | — | Arch risk | High | ✅ Fixed 2026-05-31 | +| MINOR-1 | `spherical_functional.hpp` | 404 | Doc error | Minor | ✅ Fixed 2026-05-31 | +| MINOR-2 | `spherical_functional.hpp` | 470–471 | Accuracy | Minor | ✅ Fixed 2026-05-31 | +| MINOR-3 | three files | — | DRY | Minor | ✅ Fixed 2026-05-31 | +| MINOR-4 | four files | — | DRY | Minor | ✅ Fixed 2026-05-31 | +| MINOR-5 | `clausen.hpp` | 33–38 | Doc | Minor | ✅ Fixed 2026-05-31 | + +--- + +## What was verified as correct (no action needed) + +The following were carefully audited and found faithful: + +- **`clausen.hpp`** — Chebyshev expansion, `csevl`, `inits`, `clausen2`, `Lobachevsky`, + `ImLi2`. All match Java oracle to 1e-12. +- **`euclidean_geometry.hpp`** — `euclidean_angles()` and `euclidean_angles_from_lengths()`. + Degenerate handling correct (π/0/0), centering trick correct. +- **`spherical_geometry.hpp`** — `spherical_l()`, `spherical_angles()`. Degenerate + handling correct. Half-angle formula numerically stable. +- **`euclidean_functional.hpp`** — gradient `G_v = Θ_v − Σα_v`, edge gradient + `G_e = α_opp⁺ + α_opp⁻ − φ_e`, energy via Gauss-Legendre path integral: all correct. +- **`spherical_functional.hpp`** — vertex path: correct. Edge-DOF replacement + parameterization (Finding 3 in java-port-audit.md): correct. +- **`euclidean_hessian.hpp`** — `euclidean_cot_weights()` formula and sign are correct + (only the *comment* is wrong, see Finding C). Analytic cyclic Hessian `∂α_i/∂s_j` + formula is correct and internally consistent. +- **`hyper_ideal_geometry.hpp`** — `lij`, `sigma_i`, `sigma_ij`, `alpha_ij`, + `zeta`, `zeta13/14/15`: all match Java oracle. +- **`newton_solver.hpp`** — merit `f = ½‖G‖²`, Armijo conditions Phase 1 and 2, + steepest-descent fallback `d_sd = −H·G` (correct descent direction even for NSD H), + spherical sign flip `(−H)·Δx = G`: all mathematically correct. +- **`gauss_bonnet.hpp`** — `χ = V−E+F`, `g = (2−χ)/2`, `enforce` shift + `δ = (lhs−rhs)/V`: correct for Euclidean and Spherical. Wrong for HyperIdeal (Finding B). +- **`cp_euclidean_functional.hpp`** — `p_function`, energy, gradient: match Java + BPS-2010 formula. Hessian `h_jk = sin θ / (cosh Δρ − cos θ)`: correct. +- **`inversive_distance_functional.hpp`** — `edge_length_squared`, gradient + (`Θ − Σα` pattern), degenerate-face limiting angles (Finding 9 fix correct). diff --git a/doc/roadmap/research-track.md b/doc/roadmap/research-track.md index e9cbc0d..d6152e5 100644 --- a/doc/roadmap/research-track.md +++ b/doc/roadmap/research-track.md @@ -133,6 +133,66 @@ The phase numbers match `doc/roadmap/phases.md`. ## Planned research (not yet PR) +### Hyper-ideal volume formulas for 2- and 3-ideal-vertex faces (Phase 9b+, 🔲 planned) + +* **Mathematical sources:** + - **Kolpakov, A. & Mednykh, A.** (2012). *Spherical structures on torus + knots and links.* Sibirsk. Mat. Zh. 53(3), 535–541 — see the earlier + arXiv:math/0603097 for the one-ideal-vertex formula already implemented + as `calculateTetrahedronVolumeWithIdealVertexAtGamma`. + - **Milnor, J.** (1982). *Hyperbolic geometry: The first 150 years.* + Bull. Amer. Math. Soc. 6(1), 9–24. → Volume of an ideal tetrahedron + via Clausen function; this is the all-ideal case with 4 ideal vertices. + - **Vinberg, E. B.** (1985). *Hyperbolic reflection groups.* + Uspekhi Mat. Nauk 40(1), 29–66. → general semi-ideal / orthoscheme approach. + - Study arXiv:math/0603097 §3–4 carefully to determine whether the + one-ideal formula already yields the correct limit as γ₂ → 0 + (ideal v2): if Л(0) = 0 absorbs the second ideal vertex naturally, + the extension to 2-ideal may be free; if not, a different formula is needed. + +* **Java reference:** ❌ **none.** `HyperIdealUtility.java` has exactly two + volume functions; the Java `HyperIdealFunctional` silently falls through the + `if/else-if` cascade for 2+ideal faces (uses the one-ideal formula for the + first ideal vertex found, ignoring subsequent ideal vertices). C++ now + throws `std::logic_error` instead (fixed 2026-05-30, Finding-A). + +* **Context:** + In the standard workflow (`assign_all_dof_indices`) every vertex is + hyper-ideal and no face has ideal vertices — the currently missing formulas + are never reached. They only matter for: + (a) mixed configurations with some pinned (ideal) vertices; and + (b) cusped hyperbolic surfaces (Θᵥ = 0 for a cusp vertex). + Use case (b) is the main motivation for eventually implementing these. + +* **Acceptance criteria:** + - Identify the correct formula for a hyper-ideal tetrahedron with exactly + 2 ideal vertices from the literature (check Kolpakov-Mednykh generalisations + and Vinberg orthoscheme decomposition). + - Implement `calculateTetrahedronVolumeWithTwoIdealVertices(…)` analogous + to the existing Kolpakov-Mednykh function. + - Implement `calculateTetrahedronVolumeWithThreeIdealVertices(…)` (one + hyper-ideal + three ideal = fully cusp-like case). + - Replace the `throw std::logic_error` in `face_energy()` with the correct + branch for each case; update the guard to throw only for `ideal_count > 3` + (which is topologically impossible). + - Gradient check passes for each new configuration at machine precision + (central FD vs. analytic, tol = 1e-4). + - Limiting-behaviour test: as `b_v → 0` for a hyper-ideal vertex, the + energy from the 0-ideal formula must converge to the 1-ideal formula to + 1e-6 (continuity witness). + +* **Effort:** medium (3–5 days: 1–2 days literature study + derivation, + 1–2 days implementation, 1 day tests). + +* **Phase:** 9b+ (add to Phase 9b milestone once the analytic Hessian PR lands; + the two features are independent). + +* **Note:** The `throw` introduced in the 2026-05-30 fix is the **correct + safe behaviour** until this item is resolved. Do not remove it without + implementing and testing the replacement formulas. + +--- + ### Hyper-ideal Hessian — full analytic (Phase 9b-analytic, 🔲 planned) * **Mathematical sources:**