feat(phase3d+3e): port EuclideanCyclicFunctional + add SphericalFunctional gauge-fix
Phase 3d — EuclideanCyclicFunctional:
• euclidean_geometry.hpp: t-value / atan2 corner-angle formula with
centering trick (μ = (Λ̃₁₂+Λ̃₂₃+Λ̃₃₁)/6) for numerical stability
• euclidean_functional.hpp: EuclideanMaps bundle, gradient (G_v = Θ_v − Σα_v,
G_e = α_opp⁺ + α_opp⁻ − φ_e), 10-point GL path-integral energy,
gradient_check_euclidean — identical halfedge convention to SphericalFunctional
• test_euclidean_functional.cpp: 11 tests (1 skip) covering angle formula,
right-isosceles triangle, angle sum = π, degenerate detection, gradient
checks on triangle/quad-strip/tetrahedron/fan-5/mixed-pinned, NaN check
Phase 3e — Spherical gauge-fix:
• spherical_gauge_shift(): Newton + backtracking line search to find t*
where ΣG_v(x + t·1) = 0 (maximises E along the global scale direction);
bisection used when sign change is detectable, Newton+backtrack otherwise
• apply_spherical_gauge(): in-place wrapper
• 3 new tests: GaugeFix_SpherTetVertexZerosSumGv, GaugeFix_ApplyInPlace,
GaugeFix_AlreadyAtGaugeReturnsTNearZero
Total: 45 cgal tests pass, 3 skipped (@Ignore Hessian stubs, one per functional)
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This commit is contained in:
Tarik Moussa
41
README.md
41
README.md
@@ -4,7 +4,7 @@ conformallab++ is a modern C++ reimplementation of the [ConformalLab](https://gi
|
|||||||
|
|
||||||
The long-term goal is a **CGAL package** that brings discrete conformal maps (hyper-ideal, spherical, Euclidean) to the CGAL ecosystem using `CGAL::Surface_mesh` as the underlying half-edge data structure.
|
The long-term goal is a **CGAL package** that brings discrete conformal maps (hyper-ideal, spherical, Euclidean) to the CGAL ecosystem using `CGAL::Surface_mesh` as the underlying half-edge data structure.
|
||||||
|
|
||||||
> **Status:** Phase 3 abgeschlossen. Kern-Geometrie (Clausen, Hyper-Ideal, Sphärisch) und CGAL-Mesh-Infrastruktur sind vorhanden. Zum Abschluss von Phase 3 und als Vorbereitung für Phase 4 stehen noch sechs konkrete Schritte aus — siehe [Roadmap](#roadmap).
|
> **Status:** Phasen 3d + 3e abgeschlossen. Alle drei Kern-Geometrien (Hyper-Ideal, Sphärisch, Euklidisch) plus Sphärischer Gauge-Fix sind auf `ConformalMesh` portiert. **45 Tests, 3 skipped** (Hessian-Stubs). Nächster Schritt: Phase 3f (Hessian) → Phase 4 (Newton-Solver) — siehe [Roadmap](#roadmap).
|
||||||
|
|
||||||
---
|
---
|
||||||
|
|
||||||
@@ -17,7 +17,8 @@ The long-term goal is a **CGAL package** that brings discrete conformal maps (hy
|
|||||||
| Hyper-ideal functional (energy + gradient, CGAL mesh) | ✅ Phase 3b |
|
| Hyper-ideal functional (energy + gradient, CGAL mesh) | ✅ Phase 3b |
|
||||||
| Spherical functional (energy + gradient, CGAL mesh) | ✅ Phase 3c |
|
| Spherical functional (energy + gradient, CGAL mesh) | ✅ Phase 3c |
|
||||||
| CGAL `Surface_mesh` infrastructure + mesh builders | ✅ Phase 3a |
|
| CGAL `Surface_mesh` infrastructure + mesh builders | ✅ Phase 3a |
|
||||||
| Euclidean functional | 🔜 Phase 4 |
|
| Euclidean functional (energy + gradient, CGAL mesh) | ✅ Phase 3d |
|
||||||
|
| Spherical gauge-fix (scale gauge for closed surfaces) | ✅ Phase 3e |
|
||||||
| Solvers (Newton, gradient flow) | 🔜 Phase 4 |
|
| Solvers (Newton, gradient flow) | 🔜 Phase 4 |
|
||||||
| XML serialisation / mesh I/O | 🔜 Phase 5 |
|
| XML serialisation / mesh I/O | 🔜 Phase 5 |
|
||||||
| Full CLI app | 🔜 Phase 5 |
|
| Full CLI app | 🔜 Phase 5 |
|
||||||
@@ -72,7 +73,7 @@ cmake --build build --target conformallab_cgal_tests -j$(nproc)
|
|||||||
ctest --test-dir build -R "^cgal\." --output-on-failure
|
ctest --test-dir build -R "^cgal\." --output-on-failure
|
||||||
```
|
```
|
||||||
|
|
||||||
Expected output: **30 tests pass, 2 skipped** (the two `@Ignore` Hessian stubs).
|
Expected output: **45 tests pass, 3 skipped** (the three `@Ignore` Hessian stubs, one per functional).
|
||||||
|
|
||||||
### Full CLI app (CGAL + viewer)
|
### Full CLI app (CGAL + viewer)
|
||||||
|
|
||||||
@@ -93,7 +94,9 @@ cmake --build build -j$(nproc)
|
|||||||
| `hyper_ideal_utility.hpp` | Tetrahedron volume (Meyerhoff / Kolpakov–Mednykh) |
|
| `hyper_ideal_utility.hpp` | Tetrahedron volume (Meyerhoff / Kolpakov–Mednykh) |
|
||||||
| `hyper_ideal_functional.hpp` | Hyper-ideal energy + gradient on `ConformalMesh` |
|
| `hyper_ideal_functional.hpp` | Hyper-ideal energy + gradient on `ConformalMesh` |
|
||||||
| `spherical_geometry.hpp` | Spherical arc length, half-angle angle formula |
|
| `spherical_geometry.hpp` | Spherical arc length, half-angle angle formula |
|
||||||
| `spherical_functional.hpp` | Spherical energy + gradient on `ConformalMesh` |
|
| `spherical_functional.hpp` | Spherical energy + gradient + gauge-fix on `ConformalMesh` |
|
||||||
|
| `euclidean_geometry.hpp` | Euclidean corner-angle formula (t-value / atan2) |
|
||||||
|
| `euclidean_functional.hpp` | Euclidean energy + gradient on `ConformalMesh` |
|
||||||
| `conformal_mesh.hpp` | `ConformalMesh` = `CGAL::Surface_mesh<Point3>`, index types, property-map helpers |
|
| `conformal_mesh.hpp` | `ConformalMesh` = `CGAL::Surface_mesh<Point3>`, index types, property-map helpers |
|
||||||
| `mesh_builder.hpp` | Factory meshes: triangle, tetrahedron, quad strip, fan, spherical tetrahedron, octahedron face |
|
| `mesh_builder.hpp` | Factory meshes: triangle, tetrahedron, quad strip, fan, spherical tetrahedron, octahedron face |
|
||||||
| `discrete_elliptic_utility.hpp` | Discrete elliptic integrals |
|
| `discrete_elliptic_utility.hpp` | Discrete elliptic integrals |
|
||||||
@@ -112,7 +115,9 @@ code/
|
|||||||
│ ├── hyper_ideal_geometry.hpp # ζ, lᵢⱼ, αᵢⱼ — pure math
|
│ ├── hyper_ideal_geometry.hpp # ζ, lᵢⱼ, αᵢⱼ — pure math
|
||||||
│ ├── hyper_ideal_functional.hpp # HyperIdealFunctional on ConformalMesh
|
│ ├── hyper_ideal_functional.hpp # HyperIdealFunctional on ConformalMesh
|
||||||
│ ├── spherical_geometry.hpp # spherical arc length + angles
|
│ ├── spherical_geometry.hpp # spherical arc length + angles
|
||||||
│ ├── spherical_functional.hpp # SphericalFunctional on ConformalMesh
|
│ ├── spherical_functional.hpp # SphericalFunctional + gauge-fix on ConformalMesh
|
||||||
|
│ ├── euclidean_geometry.hpp # Euclidean corner-angle formula (t-value)
|
||||||
|
│ ├── euclidean_functional.hpp # EuclideanCyclicFunctional on ConformalMesh
|
||||||
│ ├── clausen.hpp # Clausen / Lobachevsky / ImLi₂
|
│ ├── clausen.hpp # Clausen / Lobachevsky / ImLi₂
|
||||||
│ └── hyper_ideal_utility.hpp # Tetrahedron volumes
|
│ └── hyper_ideal_utility.hpp # Tetrahedron volumes
|
||||||
├── src/
|
├── src/
|
||||||
@@ -125,7 +130,8 @@ code/
|
|||||||
│ ├── CMakeLists.txt
|
│ ├── CMakeLists.txt
|
||||||
│ ├── test_conformal_mesh.cpp # 14 mesh infrastructure tests
|
│ ├── test_conformal_mesh.cpp # 14 mesh infrastructure tests
|
||||||
│ ├── test_hyper_ideal_functional.cpp # 6 hyper-ideal gradient checks
|
│ ├── test_hyper_ideal_functional.cpp # 6 hyper-ideal gradient checks
|
||||||
│ └── test_spherical_functional.cpp # 8 spherical gradient checks
|
│ ├── test_spherical_functional.cpp # 11 spherical gradient + gauge-fix checks
|
||||||
|
│ └── test_euclidean_functional.cpp # 11 Euclidean gradient checks
|
||||||
└── deps/
|
└── deps/
|
||||||
├── tarballs/ # bundled dependency archives
|
├── tarballs/ # bundled dependency archives
|
||||||
├── eigen-3.4.0/ # header-only linear algebra (always extracted)
|
├── eigen-3.4.0/ # header-only linear algebra (always extracted)
|
||||||
@@ -159,7 +165,8 @@ CGAL `Surface_mesh` tests (test prefix `cgal.`):
|
|||||||
| `ConformalMeshProperties` | 5 | Property maps: λ, θ, idx, α, geometry type |
|
| `ConformalMeshProperties` | 5 | Property maps: λ, θ, idx, α, geometry type |
|
||||||
| `ConformalMeshValidity` | 1 | CGAL validity check for all factory meshes |
|
| `ConformalMeshValidity` | 1 | CGAL validity check for all factory meshes |
|
||||||
| `HyperIdealFunctional` | 6 | Finite-difference gradient checks on triangle, tetrahedron, quad strip, fan |
|
| `HyperIdealFunctional` | 6 | Finite-difference gradient checks on triangle, tetrahedron, quad strip, fan |
|
||||||
| `SphericalFunctional` | 8 | Angle formula correctness + gradient checks on spherical meshes |
|
| `SphericalFunctional` | 11 | Angle formula + gradient checks on spherical meshes + 3 gauge-fix tests |
|
||||||
|
| `EuclideanFunctional` | 11 | Angle formula + gradient checks on Euclidean meshes (triangle, quad strip, fan, tetrahedron) |
|
||||||
|
|
||||||
---
|
---
|
||||||
|
|
||||||
@@ -200,19 +207,15 @@ bereiten Phase 4 (Solver) vor. Die Zeitangaben sind grobe Schätzungen.
|
|||||||
---
|
---
|
||||||
|
|
||||||
```
|
```
|
||||||
Phase 3d EuclideanCyclicFunctional portieren (2–3 Tage)
|
Phase 3d EuclideanCyclicFunctional portieren ✅ abgeschlossen
|
||||||
→ Wichtigster fehlender Baustein: Euklidische Uniformisierung
|
→ euclidean_geometry.hpp: t-Wert / atan2 Winkelformel
|
||||||
(genus-0 Flachparameterisierung, Kotangenten-Laplace)
|
→ euclidean_functional.hpp: Energie + Gradient auf ConformalMesh
|
||||||
→ Energie + Gradient auf ConformalMesh
|
→ 11 Tests: Winkelformel, Gradient-Checks, NaN-Test, Fan, Mixed-Pinned
|
||||||
→ Entsperrt die verbleibenden HDS-blockierten Tests
|
|
||||||
→ Voraussetzung für jeden end-to-end Optimiererlauf
|
|
||||||
|
|
||||||
Phase 3e Gauge-Fix für SphericalFunctional (1 Tag)
|
Phase 3e Gauge-Fix für SphericalFunctional ✅ abgeschlossen
|
||||||
→ Java-Äquivalent: maximizeInNegativeDirection() vor jeder
|
→ spherical_gauge_shift() + apply_spherical_gauge() in
|
||||||
Auswertung — 1D-Liniensuche (Brent) entlang der
|
spherical_functional.hpp — Newton + Backtracking
|
||||||
globalen Additionskonstante in u_v
|
→ 3 neue Tests: ZerosSumGv, ApplyInPlace, AlreadyAtGauge
|
||||||
→ Ohne Gauge-Fix ist die Energie auf geschlossenen Flächen
|
|
||||||
nach unten unbeschränkt; kein Optimierer kann konvergieren
|
|
||||||
|
|
||||||
Phase 3f Hessian (Sphärisch → Euklidisch) (2–3 Tage)
|
Phase 3f Hessian (Sphärisch → Euklidisch) (2–3 Tage)
|
||||||
→ Analytischer Hessian = Kotangenten-Laplace-Operator
|
→ Analytischer Hessian = Kotangenten-Laplace-Operator
|
||||||
|
|||||||
325
code/include/euclidean_functional.hpp
Normal file
325
code/include/euclidean_functional.hpp
Normal file
@@ -0,0 +1,325 @@
|
|||||||
|
#pragma once
|
||||||
|
// euclidean_functional.hpp
|
||||||
|
//
|
||||||
|
// Energy and gradient of the Euclidean discrete conformal functional
|
||||||
|
// (EuclideanCyclicFunctional) evaluated on a ConformalMesh.
|
||||||
|
//
|
||||||
|
// Ported from de.varylab.discreteconformal.functional.EuclideanCyclicFunctional.
|
||||||
|
//
|
||||||
|
// ┌──────────────────────────────────────────────────────────────────────────┐
|
||||||
|
// │ DOFs │
|
||||||
|
// │ x[v_idx[v]] = u_v – conformal factor at vertex v │
|
||||||
|
// │ x[e_idx[e]] = λ_e – edge log-length variable (optional) │
|
||||||
|
// │ -1 means "pinned" (u_v = 0 / λ_e = 0, only λ° contributes) │
|
||||||
|
// │ │
|
||||||
|
// │ Effective log-length (always additive, unlike SphericalFunctional): │
|
||||||
|
// │ Λ̃_ij = λ°_ij + u_i + u_j + (x[e_idx[e]] if variable, else 0) │
|
||||||
|
// │ │
|
||||||
|
// │ Side length: l_ij = exp(Λ̃_ij / 2) │
|
||||||
|
// │ │
|
||||||
|
// │ Gradient: │
|
||||||
|
// │ ∂E/∂u_v = Θ_v − Σ_{faces adj. v} α_v(face) │
|
||||||
|
// │ ∂E/∂λ_e = α_opp(face⁺) + α_opp(face⁻) − φ_e │
|
||||||
|
// │ │
|
||||||
|
// │ Energy: │
|
||||||
|
// │ Computed as the path integral E(x) = ∫₀¹ ⟨G(tx), x⟩ dt │
|
||||||
|
// │ using 10-point Gauss-Legendre quadrature (same as SphericalFunctional)│
|
||||||
|
// │ This is E(0)=0 by construction and exact for conservative G. │
|
||||||
|
// └──────────────────────────────────────────────────────────────────────────┘
|
||||||
|
//
|
||||||
|
// Halfedge convention (identical to SphericalFunctional):
|
||||||
|
// h0 = mesh.halfedge(f), h1 = next(h0), h2 = next(h1)
|
||||||
|
// v1 = source(h0), v2 = source(h1), v3 = source(h2)
|
||||||
|
// h_alpha[h0] = α3 (angle at v3, opposite edge h0 = v1v2)
|
||||||
|
// h_alpha[h1] = α1 (angle at v1, opposite edge h1 = v2v3)
|
||||||
|
// h_alpha[h2] = α2 (angle at v2, opposite edge h2 = v3v1)
|
||||||
|
//
|
||||||
|
// Property-map name prefix: "ev:" (vertex) and "ee:" (edge).
|
||||||
|
|
||||||
|
#include "conformal_mesh.hpp"
|
||||||
|
#include "euclidean_geometry.hpp"
|
||||||
|
#include <CGAL/boost/graph/iterator.h>
|
||||||
|
#include <vector>
|
||||||
|
#include <cmath>
|
||||||
|
#include <cstdint>
|
||||||
|
|
||||||
|
namespace conformallab {
|
||||||
|
|
||||||
|
// ── Property-map type aliases ─────────────────────────────────────────────────
|
||||||
|
|
||||||
|
using EuclVMapD = ConformalMesh::Property_map<Vertex_index, double>;
|
||||||
|
using EuclVMapI = ConformalMesh::Property_map<Vertex_index, int>;
|
||||||
|
using EuclEMapD = ConformalMesh::Property_map<Edge_index, double>;
|
||||||
|
using EuclEMapI = ConformalMesh::Property_map<Edge_index, int>;
|
||||||
|
|
||||||
|
// ── Persistent map bundle ─────────────────────────────────────────────────────
|
||||||
|
|
||||||
|
struct EuclideanMaps {
|
||||||
|
EuclVMapI v_idx; ///< DOF index per vertex (-1 = pinned / u_v = 0)
|
||||||
|
EuclEMapI e_idx; ///< DOF index per edge (-1 = no edge DOF)
|
||||||
|
EuclVMapD theta_v; ///< target cone angle Θ_v (default 2π)
|
||||||
|
EuclEMapD phi_e; ///< target edge turn angle φ_e (default π)
|
||||||
|
EuclEMapD lambda0; ///< base log-length λ°_e (default 0.0)
|
||||||
|
};
|
||||||
|
|
||||||
|
// Create and attach property maps with sensible defaults.
|
||||||
|
// theta_v = 2π (flat vertex), phi_e = π (interior edge, flat surface).
|
||||||
|
inline EuclideanMaps setup_euclidean_maps(ConformalMesh& mesh)
|
||||||
|
{
|
||||||
|
constexpr double TWO_PI = 2.0 * 3.14159265358979323846;
|
||||||
|
constexpr double PI = 3.14159265358979323846;
|
||||||
|
|
||||||
|
EuclideanMaps m;
|
||||||
|
m.v_idx = mesh.add_property_map<Vertex_index, int> ("ev:idx", -1 ).first;
|
||||||
|
m.e_idx = mesh.add_property_map<Edge_index, int> ("ee:idx", -1 ).first;
|
||||||
|
m.theta_v= mesh.add_property_map<Vertex_index, double>("ev:theta", TWO_PI ).first;
|
||||||
|
m.phi_e = mesh.add_property_map<Edge_index, double>("ee:phi", PI ).first;
|
||||||
|
m.lambda0= mesh.add_property_map<Edge_index, double>("ee:lam0", 0.0 ).first;
|
||||||
|
return m;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Assign DOF indices 0..n-1 for all vertices only (no edge DOFs).
|
||||||
|
inline int assign_euclidean_vertex_dof_indices(ConformalMesh& mesh, EuclideanMaps& m)
|
||||||
|
{
|
||||||
|
int idx = 0;
|
||||||
|
for (auto v : mesh.vertices()) m.v_idx[v] = idx++;
|
||||||
|
return idx;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Assign DOF indices for all vertices AND all edges.
|
||||||
|
inline int assign_euclidean_all_dof_indices(ConformalMesh& mesh, EuclideanMaps& m)
|
||||||
|
{
|
||||||
|
int idx = 0;
|
||||||
|
for (auto v : mesh.vertices()) m.v_idx[v] = idx++;
|
||||||
|
for (auto e : mesh.edges()) m.e_idx[e] = idx++;
|
||||||
|
return idx;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Count variable DOFs (vertices + edges).
|
||||||
|
inline int euclidean_dimension(const ConformalMesh& mesh, const EuclideanMaps& m)
|
||||||
|
{
|
||||||
|
int dim = 0;
|
||||||
|
for (auto v : mesh.vertices()) if (m.v_idx[v] >= 0) ++dim;
|
||||||
|
for (auto e : mesh.edges()) if (m.e_idx[e] >= 0) ++dim;
|
||||||
|
return dim;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Set lambda0 from mesh vertex positions (Euclidean):
|
||||||
|
// λ°_e = 2·log(|p_i − p_j|) (natural log of Euclidean edge length squared)
|
||||||
|
//
|
||||||
|
// This gives exp(Λ̃_ij / 2) = l_ij at x=0.
|
||||||
|
inline void compute_euclidean_lambda0_from_mesh(ConformalMesh& mesh, EuclideanMaps& m)
|
||||||
|
{
|
||||||
|
for (auto e : mesh.edges()) {
|
||||||
|
auto h = mesh.halfedge(e);
|
||||||
|
auto p1 = mesh.point(mesh.source(h));
|
||||||
|
auto p2 = mesh.point(mesh.target(h));
|
||||||
|
double dx = p1.x() - p2.x();
|
||||||
|
double dy = p1.y() - p2.y();
|
||||||
|
double dz = p1.z() - p2.z();
|
||||||
|
double len = std::sqrt(dx*dx + dy*dy + dz*dz);
|
||||||
|
if (len > 1e-15)
|
||||||
|
m.lambda0[e] = 2.0 * std::log(len);
|
||||||
|
else
|
||||||
|
m.lambda0[e] = -30.0; // degenerate edge
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// ── Internal helpers ──────────────────────────────────────────────────────────
|
||||||
|
|
||||||
|
static inline double eucl_dof_val(int idx, const std::vector<double>& x)
|
||||||
|
{
|
||||||
|
return idx >= 0 ? x[static_cast<std::size_t>(idx)] : 0.0;
|
||||||
|
}
|
||||||
|
|
||||||
|
static inline std::size_t eucl_hidx(Halfedge_index h)
|
||||||
|
{
|
||||||
|
return static_cast<std::size_t>(static_cast<std::uint32_t>(h));
|
||||||
|
}
|
||||||
|
|
||||||
|
// ── Gradient ──────────────────────────────────────────────────────────────────
|
||||||
|
//
|
||||||
|
// G_v = Θ_v − Σ_{faces adj. v} α_v(face)
|
||||||
|
// G_e = α_opp(face⁺) + α_opp(face⁻) − φ_e
|
||||||
|
//
|
||||||
|
// Corner-angle storage (h_alpha):
|
||||||
|
// h_alpha[h] = corner angle OPPOSITE to the edge of halfedge h in its face.
|
||||||
|
// h_alpha[h0] = α3, h_alpha[h1] = α1, h_alpha[h2] = α2
|
||||||
|
// (same convention as SphericalFunctional)
|
||||||
|
inline std::vector<double> euclidean_gradient(
|
||||||
|
ConformalMesh& mesh,
|
||||||
|
const std::vector<double>& x,
|
||||||
|
const EuclideanMaps& m)
|
||||||
|
{
|
||||||
|
const int n = euclidean_dimension(mesh, m);
|
||||||
|
std::vector<double> G(static_cast<std::size_t>(n), 0.0);
|
||||||
|
|
||||||
|
// Per-halfedge corner-angle storage.
|
||||||
|
const std::size_t nh = mesh.number_of_halfedges();
|
||||||
|
std::vector<double> h_alpha(nh, 0.0);
|
||||||
|
|
||||||
|
// ── Pass 1: compute corner angles per face ────────────────────────────────
|
||||||
|
for (auto f : mesh.faces()) {
|
||||||
|
Halfedge_index h0 = mesh.halfedge(f);
|
||||||
|
Halfedge_index h1 = mesh.next(h0);
|
||||||
|
Halfedge_index h2 = mesh.next(h1);
|
||||||
|
|
||||||
|
Vertex_index v1 = mesh.source(h0);
|
||||||
|
Vertex_index v2 = mesh.source(h1);
|
||||||
|
Vertex_index v3 = mesh.source(h2);
|
||||||
|
|
||||||
|
Edge_index e12 = mesh.edge(h0);
|
||||||
|
Edge_index e23 = mesh.edge(h1);
|
||||||
|
Edge_index e31 = mesh.edge(h2);
|
||||||
|
|
||||||
|
// Effective log-lengths (always additive: u_i + u_j regardless of DOF status)
|
||||||
|
double u1 = eucl_dof_val(m.v_idx[v1], x);
|
||||||
|
double u2 = eucl_dof_val(m.v_idx[v2], x);
|
||||||
|
double u3 = eucl_dof_val(m.v_idx[v3], x);
|
||||||
|
|
||||||
|
double lam12 = m.lambda0[e12] + u1 + u2 + eucl_dof_val(m.e_idx[e12], x);
|
||||||
|
double lam23 = m.lambda0[e23] + u2 + u3 + eucl_dof_val(m.e_idx[e23], x);
|
||||||
|
double lam31 = m.lambda0[e31] + u3 + u1 + eucl_dof_val(m.e_idx[e31], x);
|
||||||
|
|
||||||
|
auto fa = euclidean_angles(lam12, lam23, lam31);
|
||||||
|
if (!fa.valid) continue; // degenerate triangle: contributes 0
|
||||||
|
|
||||||
|
// h_alpha[h] = corner angle OPPOSITE to h's edge:
|
||||||
|
// h0 (edge v1v2) → opposite corner at v3 → α3
|
||||||
|
// h1 (edge v2v3) → opposite corner at v1 → α1
|
||||||
|
// h2 (edge v3v1) → opposite corner at v2 → α2
|
||||||
|
h_alpha[eucl_hidx(h0)] = fa.alpha3;
|
||||||
|
h_alpha[eucl_hidx(h1)] = fa.alpha1;
|
||||||
|
h_alpha[eucl_hidx(h2)] = fa.alpha2;
|
||||||
|
}
|
||||||
|
|
||||||
|
// ── Pass 2: accumulate vertex gradient ───────────────────────────────────
|
||||||
|
// G_v = Θ_v − Σ h_alpha[prev(h)] for each incoming non-border h to v.
|
||||||
|
for (auto v : mesh.vertices()) {
|
||||||
|
int iv = m.v_idx[v];
|
||||||
|
if (iv < 0) continue;
|
||||||
|
double sum_alpha = 0.0;
|
||||||
|
for (auto h : CGAL::halfedges_around_target(v, mesh)) {
|
||||||
|
if (mesh.is_border(h)) continue;
|
||||||
|
sum_alpha += h_alpha[eucl_hidx(mesh.prev(h))];
|
||||||
|
}
|
||||||
|
G[static_cast<std::size_t>(iv)] = m.theta_v[v] - sum_alpha;
|
||||||
|
}
|
||||||
|
|
||||||
|
// ── Pass 3: accumulate edge gradient ─────────────────────────────────────
|
||||||
|
// G_e = α_opp(f⁺) + α_opp(f⁻) − φ_e
|
||||||
|
// α_opp of edge e in face f = h_alpha[halfedge h of e pointing INTO f].
|
||||||
|
for (auto e : mesh.edges()) {
|
||||||
|
int ie = m.e_idx[e];
|
||||||
|
if (ie < 0) continue;
|
||||||
|
auto h = mesh.halfedge(e);
|
||||||
|
auto ho = mesh.opposite(h);
|
||||||
|
double sum = -m.phi_e[e];
|
||||||
|
if (!mesh.is_border(h)) sum += h_alpha[eucl_hidx(h)];
|
||||||
|
if (!mesh.is_border(ho)) sum += h_alpha[eucl_hidx(ho)];
|
||||||
|
G[static_cast<std::size_t>(ie)] = sum;
|
||||||
|
}
|
||||||
|
|
||||||
|
return G;
|
||||||
|
}
|
||||||
|
|
||||||
|
// ── Energy via Gauss-Legendre path integral ───────────────────────────────────
|
||||||
|
//
|
||||||
|
// E(x) = ∫₀¹ ⟨G(tx), x⟩ dt (10-point GL quadrature, same as SphericalFunctional)
|
||||||
|
inline double euclidean_energy(
|
||||||
|
ConformalMesh& mesh,
|
||||||
|
const std::vector<double>& x,
|
||||||
|
const EuclideanMaps& m)
|
||||||
|
{
|
||||||
|
static const double gl_s[10] = {
|
||||||
|
-0.9739065285171717, -0.8650633666889845,
|
||||||
|
-0.6794095682990244, -0.4333953941292472,
|
||||||
|
-0.1488743389816312, 0.1488743389816312,
|
||||||
|
0.4333953941292472, 0.6794095682990244,
|
||||||
|
0.8650633666889845, 0.9739065285171717
|
||||||
|
};
|
||||||
|
static const double gl_w[10] = {
|
||||||
|
0.0666713443086881, 0.1494513491505806,
|
||||||
|
0.2190863625159820, 0.2692667193099963,
|
||||||
|
0.2955242247147529, 0.2955242247147529,
|
||||||
|
0.2692667193099963, 0.2190863625159820,
|
||||||
|
0.1494513491505806, 0.0666713443086881
|
||||||
|
};
|
||||||
|
|
||||||
|
const std::size_t n = x.size();
|
||||||
|
double E = 0.0;
|
||||||
|
|
||||||
|
for (int k = 0; k < 10; ++k) {
|
||||||
|
double t = (1.0 + gl_s[k]) * 0.5;
|
||||||
|
double wt = gl_w[k] * 0.5;
|
||||||
|
|
||||||
|
std::vector<double> tx(n);
|
||||||
|
for (std::size_t i = 0; i < n; ++i) tx[i] = t * x[i];
|
||||||
|
|
||||||
|
auto G = euclidean_gradient(mesh, tx, m);
|
||||||
|
|
||||||
|
double dot = 0.0;
|
||||||
|
for (std::size_t i = 0; i < n; ++i) dot += G[i] * x[i];
|
||||||
|
E += wt * dot;
|
||||||
|
}
|
||||||
|
return E;
|
||||||
|
}
|
||||||
|
|
||||||
|
// ── Full evaluation (energy + gradient) ──────────────────────────────────────
|
||||||
|
|
||||||
|
struct EuclideanResult {
|
||||||
|
double energy = 0.0;
|
||||||
|
std::vector<double> gradient;
|
||||||
|
};
|
||||||
|
|
||||||
|
inline EuclideanResult evaluate_euclidean(
|
||||||
|
ConformalMesh& mesh,
|
||||||
|
const std::vector<double>& x,
|
||||||
|
const EuclideanMaps& m,
|
||||||
|
bool need_energy = true,
|
||||||
|
bool need_gradient = true)
|
||||||
|
{
|
||||||
|
EuclideanResult res;
|
||||||
|
if (need_gradient)
|
||||||
|
res.gradient = euclidean_gradient(mesh, x, m);
|
||||||
|
if (need_energy)
|
||||||
|
res.energy = euclidean_energy(mesh, x, m);
|
||||||
|
return res;
|
||||||
|
}
|
||||||
|
|
||||||
|
// ── Finite-difference gradient check ─────────────────────────────────────────
|
||||||
|
//
|
||||||
|
// Tests |G[i] − (E(x+εeᵢ) − E(x−εeᵢ))/(2ε)| / max(1,|G[i]|) < tol
|
||||||
|
// for all variable DOFs.
|
||||||
|
inline bool gradient_check_euclidean(
|
||||||
|
ConformalMesh& mesh,
|
||||||
|
const std::vector<double>& x0,
|
||||||
|
const EuclideanMaps& m,
|
||||||
|
double eps = 1e-5,
|
||||||
|
double tol = 1e-4)
|
||||||
|
{
|
||||||
|
auto G = euclidean_gradient(mesh, x0, m);
|
||||||
|
const int n = static_cast<int>(G.size());
|
||||||
|
|
||||||
|
std::vector<double> xp = x0, xm = x0;
|
||||||
|
bool ok = true;
|
||||||
|
|
||||||
|
for (int i = 0; i < n; ++i) {
|
||||||
|
std::size_t si = static_cast<std::size_t>(i);
|
||||||
|
xp[si] = x0[si] + eps;
|
||||||
|
xm[si] = x0[si] - eps;
|
||||||
|
|
||||||
|
double Ep = euclidean_energy(mesh, xp, m);
|
||||||
|
double Em = euclidean_energy(mesh, xm, m);
|
||||||
|
|
||||||
|
xp[si] = xm[si] = x0[si]; // restore
|
||||||
|
|
||||||
|
double fd = (Ep - Em) / (2.0 * eps);
|
||||||
|
double err = std::abs(G[si] - fd);
|
||||||
|
double scale = std::max(1.0, std::abs(G[si]));
|
||||||
|
if (err / scale > tol) ok = false;
|
||||||
|
}
|
||||||
|
return ok;
|
||||||
|
}
|
||||||
|
|
||||||
|
} // namespace conformallab
|
||||||
91
code/include/euclidean_geometry.hpp
Normal file
91
code/include/euclidean_geometry.hpp
Normal file
@@ -0,0 +1,91 @@
|
|||||||
|
#pragma once
|
||||||
|
// euclidean_geometry.hpp
|
||||||
|
//
|
||||||
|
// Corner-angle formula for Euclidean triangles in the discrete conformal
|
||||||
|
// (log-length) parametrisation.
|
||||||
|
//
|
||||||
|
// Ported from de.varylab.discreteconformal.functional.EuclideanCyclicFunctional.
|
||||||
|
//
|
||||||
|
// In the discrete conformal parametrisation a Euclidean triangle is described by
|
||||||
|
// its three effective log-lengths Λ̃_ij = λ°_ij + u_i + u_j (+ edge DOF).
|
||||||
|
// The corresponding side lengths are l_ij = exp(Λ̃_ij / 2).
|
||||||
|
//
|
||||||
|
// Vertex ordering convention (matches EuclideanCyclicFunctional.java):
|
||||||
|
// v1 is opposite edge l23, v2 is opposite l31, v3 is opposite l12.
|
||||||
|
//
|
||||||
|
// t-value trick (Springborn 2008 §3):
|
||||||
|
// t12 = −l12 + l23 + l31 = 2(s − l12)
|
||||||
|
// t23 = +l12 − l23 + l31 = 2(s − l23)
|
||||||
|
// t31 = +l12 + l23 − l31 = 2(s − l31)
|
||||||
|
// denom = sqrt(t12 · t23 · t31 · l123) = 4 · Area
|
||||||
|
//
|
||||||
|
// α_v = 2 · atan2( product of t-values adjacent to v, denom )
|
||||||
|
//
|
||||||
|
// The centering trick (l_ij ← exp((Λ̃_ij − 2·μ)/2), μ = (Λ̃12+Λ̃23+Λ̃31)/6)
|
||||||
|
// rescales all three sides by the same factor, leaving angles unchanged but
|
||||||
|
// keeping the arguments of exp in a safe numerical range.
|
||||||
|
|
||||||
|
#include <cmath>
|
||||||
|
|
||||||
|
namespace conformallab {
|
||||||
|
|
||||||
|
struct EuclideanFaceAngles {
|
||||||
|
double alpha1; ///< corner angle at v1 (opposite l23)
|
||||||
|
double alpha2; ///< corner angle at v2 (opposite l31)
|
||||||
|
double alpha3; ///< corner angle at v3 (opposite l12)
|
||||||
|
bool valid;
|
||||||
|
};
|
||||||
|
|
||||||
|
// ── From side lengths ─────────────────────────────────────────────────────────
|
||||||
|
//
|
||||||
|
// Given three Euclidean side lengths l12, l23, l31 > 0 satisfying the triangle
|
||||||
|
// inequality, compute the corner angles.
|
||||||
|
//
|
||||||
|
// Returns valid=false if the triangle inequality is violated (any t-value ≤ 0).
|
||||||
|
inline EuclideanFaceAngles euclidean_angles_from_lengths(
|
||||||
|
double l12, double l23, double l31)
|
||||||
|
{
|
||||||
|
const double t12 = -l12 + l23 + l31; // 2*(s − l12)
|
||||||
|
const double t23 = +l12 - l23 + l31; // 2*(s − l23)
|
||||||
|
const double t31 = +l12 + l23 - l31; // 2*(s − l31)
|
||||||
|
|
||||||
|
if (t12 <= 0.0 || t23 <= 0.0 || t31 <= 0.0)
|
||||||
|
return {0.0, 0.0, 0.0, false};
|
||||||
|
|
||||||
|
const double l123 = l12 + l23 + l31;
|
||||||
|
const double denom2 = t12 * t23 * t31 * l123; // = (4·Area)²
|
||||||
|
if (denom2 <= 0.0)
|
||||||
|
return {0.0, 0.0, 0.0, false};
|
||||||
|
|
||||||
|
const double denom = std::sqrt(denom2);
|
||||||
|
|
||||||
|
// α at v1 (opposite l23): adjacent t-values are t12 and t31
|
||||||
|
// α at v2 (opposite l31): adjacent t-values are t12 and t23
|
||||||
|
// α at v3 (opposite l12): adjacent t-values are t23 and t31
|
||||||
|
return {
|
||||||
|
2.0 * std::atan2(t12 * t31, denom),
|
||||||
|
2.0 * std::atan2(t12 * t23, denom),
|
||||||
|
2.0 * std::atan2(t23 * t31, denom),
|
||||||
|
true
|
||||||
|
};
|
||||||
|
}
|
||||||
|
|
||||||
|
// ── From effective log-lengths Λ̃ ─────────────────────────────────────────────
|
||||||
|
//
|
||||||
|
// Converts to side lengths l_ij = exp(Λ̃_ij / 2), applying the centering
|
||||||
|
// trick for numerical safety, then delegates to euclidean_angles_from_lengths.
|
||||||
|
//
|
||||||
|
// The centering constant μ = (Λ̃12 + Λ̃23 + Λ̃31) / 6 ensures
|
||||||
|
// l12 · l23 · l31 = 1 (geometric mean = 1)
|
||||||
|
// which keeps all l values near 1 and prevents float overflow for large |Λ̃|.
|
||||||
|
inline EuclideanFaceAngles euclidean_angles(
|
||||||
|
double lam12, double lam23, double lam31)
|
||||||
|
{
|
||||||
|
const double mu = (lam12 + lam23 + lam31) / 6.0;
|
||||||
|
const double l12 = std::exp((lam12 - 2.0 * mu) * 0.5);
|
||||||
|
const double l23 = std::exp((lam23 - 2.0 * mu) * 0.5);
|
||||||
|
const double l31 = std::exp((lam31 - 2.0 * mu) * 0.5);
|
||||||
|
return euclidean_angles_from_lengths(l12, l23, l31);
|
||||||
|
}
|
||||||
|
|
||||||
|
} // namespace conformallab
|
||||||
@@ -354,4 +354,122 @@ inline bool gradient_check_spherical(
|
|||||||
return ok;
|
return ok;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
// ── Gauge-fix for closed spherical surfaces ───────────────────────────────────
|
||||||
|
//
|
||||||
|
// On a closed (boundaryless) spherical surface the energy E(u + t·1) has a
|
||||||
|
// unique maximum w.r.t. t ∈ ℝ (the "global scale" gauge mode). Without
|
||||||
|
// fixing this gauge the functional is unbounded below and no solver converges.
|
||||||
|
//
|
||||||
|
// This function returns the scalar shift t* such that
|
||||||
|
// Σ_v G_v(x + t*·1_v) = 0
|
||||||
|
// i.e., the derivative of E(u+t·1) w.r.t. t is zero at t = t*.
|
||||||
|
//
|
||||||
|
// Apply the shift by adding t* to every vertex DOF in x.
|
||||||
|
//
|
||||||
|
// Implementation: bisection on f(t) = Σ_v G_v(x + t·1_v).
|
||||||
|
// f is strictly monotone decreasing (second derivative < 0) for a convex
|
||||||
|
// functional, so bisection converges in O(log₂(2·bracket/tol)) iterations.
|
||||||
|
//
|
||||||
|
// Parameters:
|
||||||
|
// bracket – initial search interval [−bracket, +bracket] (default 50)
|
||||||
|
// tol – absolute tolerance on t* (default 1e-8)
|
||||||
|
//
|
||||||
|
// Returns 0.0 if the zero cannot be bracketed (already at gauge maximum,
|
||||||
|
// or open surface — no shift needed).
|
||||||
|
inline double spherical_gauge_shift(
|
||||||
|
ConformalMesh& mesh,
|
||||||
|
const std::vector<double>& x,
|
||||||
|
const SphericalMaps& m,
|
||||||
|
double bracket = 50.0,
|
||||||
|
double tol = 1e-8)
|
||||||
|
{
|
||||||
|
// Helper: sum of all vertex gradient components at x + t·1_v.
|
||||||
|
auto sum_Gv = [&](double t) -> double {
|
||||||
|
// Build a shifted copy of x (only vertex DOFs shifted).
|
||||||
|
std::vector<double> xt = x;
|
||||||
|
for (auto v : mesh.vertices()) {
|
||||||
|
int iv = m.v_idx[v];
|
||||||
|
if (iv >= 0)
|
||||||
|
xt[static_cast<std::size_t>(iv)] += t;
|
||||||
|
}
|
||||||
|
auto G = spherical_gradient(mesh, xt, m);
|
||||||
|
double sum = 0.0;
|
||||||
|
for (auto v : mesh.vertices()) {
|
||||||
|
int iv = m.v_idx[v];
|
||||||
|
if (iv >= 0)
|
||||||
|
sum += G[static_cast<std::size_t>(iv)];
|
||||||
|
}
|
||||||
|
return sum;
|
||||||
|
};
|
||||||
|
|
||||||
|
// ── Try bisection first (works when there is a sign change in [−bracket, +bracket]) ──
|
||||||
|
double a = -bracket, b = bracket;
|
||||||
|
double fa = sum_Gv(a), fb = sum_Gv(b);
|
||||||
|
|
||||||
|
if (fa * fb <= 0.0) {
|
||||||
|
for (int iter = 0; iter < 120; ++iter) {
|
||||||
|
double c = 0.5 * (a + b);
|
||||||
|
double fc = sum_Gv(c);
|
||||||
|
if (std::abs(fc) < tol || (b - a) < tol) return c;
|
||||||
|
if (fa * fc < 0.0) { b = c; fb = fc; }
|
||||||
|
else { a = c; fa = fc; }
|
||||||
|
}
|
||||||
|
return 0.5 * (a + b);
|
||||||
|
}
|
||||||
|
|
||||||
|
// ── No sign change (zero may lie at a domain boundary). ───────────────────
|
||||||
|
// Use damped Newton's method with backtracking line search.
|
||||||
|
// f'(t) estimated by forward finite difference.
|
||||||
|
// When the Newton step overshoots the valid domain (ΣG_v jumps back up
|
||||||
|
// because faces become degenerate), backtracking halves the step until
|
||||||
|
// |f| strictly decreases.
|
||||||
|
const double fd_eps = 1e-5;
|
||||||
|
double t = 0.0;
|
||||||
|
double ft = sum_Gv(t);
|
||||||
|
|
||||||
|
for (int iter = 0; iter < 120; ++iter) {
|
||||||
|
if (std::abs(ft) < tol) return t;
|
||||||
|
|
||||||
|
double ftp = sum_Gv(t + fd_eps);
|
||||||
|
double dft = (ftp - ft) / fd_eps;
|
||||||
|
if (std::abs(dft) < 1e-14) return t; // gradient flat — give up
|
||||||
|
|
||||||
|
double dt_raw = -ft / dft;
|
||||||
|
|
||||||
|
// Backtracking line search: halve dt until |f| decreases.
|
||||||
|
double alpha = 1.0;
|
||||||
|
bool improved = false;
|
||||||
|
for (int back = 0; back < 40; ++back) {
|
||||||
|
double t_try = t + alpha * dt_raw;
|
||||||
|
t_try = std::max(-bracket, std::min(bracket, t_try));
|
||||||
|
double ft_try = sum_Gv(t_try);
|
||||||
|
if (std::abs(ft_try) < std::abs(ft)) {
|
||||||
|
t = t_try;
|
||||||
|
ft = ft_try;
|
||||||
|
improved = true;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
alpha *= 0.5;
|
||||||
|
}
|
||||||
|
if (!improved) return t; // cannot reduce further
|
||||||
|
}
|
||||||
|
return t;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Apply the gauge shift in-place: x_v ← x_v + t* for all variable vertices.
|
||||||
|
inline void apply_spherical_gauge(
|
||||||
|
ConformalMesh& mesh,
|
||||||
|
std::vector<double>& x,
|
||||||
|
const SphericalMaps& m,
|
||||||
|
double bracket = 50.0,
|
||||||
|
double tol = 1e-8)
|
||||||
|
{
|
||||||
|
double t = spherical_gauge_shift(mesh, x, m, bracket, tol);
|
||||||
|
for (auto v : mesh.vertices()) {
|
||||||
|
int iv = m.v_idx[v];
|
||||||
|
if (iv >= 0)
|
||||||
|
x[static_cast<std::size_t>(iv)] += t;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
} // namespace conformallab
|
} // namespace conformallab
|
||||||
|
|||||||
@@ -15,8 +15,11 @@ add_executable(conformallab_cgal_tests
|
|||||||
# ── Phase 3b: HyperIdealFunctional ─────────────────────────────────────
|
# ── Phase 3b: HyperIdealFunctional ─────────────────────────────────────
|
||||||
test_hyper_ideal_functional.cpp
|
test_hyper_ideal_functional.cpp
|
||||||
|
|
||||||
# ── Phase 3c: SphericalFunctional ─────────────────────────────────────
|
# ── Phase 3c + 3e: SphericalFunctional + gauge-fix ────────────────────
|
||||||
test_spherical_functional.cpp
|
test_spherical_functional.cpp
|
||||||
|
|
||||||
|
# ── Phase 3d: EuclideanCyclicFunctional ───────────────────────────────
|
||||||
|
test_euclidean_functional.cpp
|
||||||
)
|
)
|
||||||
|
|
||||||
target_include_directories(conformallab_cgal_tests SYSTEM PRIVATE
|
target_include_directories(conformallab_cgal_tests SYSTEM PRIVATE
|
||||||
|
|||||||
269
code/tests/cgal/test_euclidean_functional.cpp
Normal file
269
code/tests/cgal/test_euclidean_functional.cpp
Normal file
@@ -0,0 +1,269 @@
|
|||||||
|
// test_euclidean_functional.cpp
|
||||||
|
//
|
||||||
|
// Phase 3d — EuclideanCyclicFunctional ported to ConformalMesh.
|
||||||
|
//
|
||||||
|
// Corresponds to de.varylab.discreteconformal.functional.EuclideanCyclicFunctionalTest.
|
||||||
|
//
|
||||||
|
// Test map (Java → C++)
|
||||||
|
// ──────────────────────
|
||||||
|
// testHessian (Ignored) → GradientCheck_Hessian (SKIPPED)
|
||||||
|
// testGradient…Triangle → GradientCheck_TriangleVertex (ported)
|
||||||
|
// testGradient…QuadStrip → GradientCheck_QuadStripVertex (ported)
|
||||||
|
// testGradient…Tetrahedron → GradientCheck_TetrahedronVertex (ported)
|
||||||
|
// testGradient…AllDofs → GradientCheck_TetrahedronAllDofs (ported)
|
||||||
|
// testFunctionalAtNaNValue → AnglesFiniteAtKnownPoint (ported)
|
||||||
|
//
|
||||||
|
// Energy model
|
||||||
|
// ────────────
|
||||||
|
// Uses the Schläfli path integral E(x) = ∫₀¹⟨G(tx),x⟩dt (10-point GL).
|
||||||
|
// The gradient check verifies G is curl-free.
|
||||||
|
|
||||||
|
#include "conformal_mesh.hpp"
|
||||||
|
#include "mesh_builder.hpp"
|
||||||
|
#include "euclidean_geometry.hpp"
|
||||||
|
#include "euclidean_functional.hpp"
|
||||||
|
#include <gtest/gtest.h>
|
||||||
|
#include <cmath>
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
using namespace conformallab;
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// @Ignore in Java: no Hessian implemented yet
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(EuclideanFunctional, GradientCheck_Hessian)
|
||||||
|
{
|
||||||
|
GTEST_SKIP() << "@Ignore in Java – Hessian not yet implemented";
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Angle formula: equilateral triangle → all angles = π/3
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(EuclideanFunctional, EquilateralTriangleAnglesArePiOver3)
|
||||||
|
{
|
||||||
|
// All sides equal: l = 1.0, log-length = 0.
|
||||||
|
auto fa = euclidean_angles(0.0, 0.0, 0.0);
|
||||||
|
|
||||||
|
ASSERT_TRUE(fa.valid) << "Equilateral triangle must be valid";
|
||||||
|
constexpr double PI_3 = 3.14159265358979323846 / 3.0;
|
||||||
|
EXPECT_NEAR(fa.alpha1, PI_3, 1e-12);
|
||||||
|
EXPECT_NEAR(fa.alpha2, PI_3, 1e-12);
|
||||||
|
EXPECT_NEAR(fa.alpha3, PI_3, 1e-12);
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Angle formula: right isosceles triangle (legs 1, hypotenuse √2)
|
||||||
|
//
|
||||||
|
// For the 45-45-90 triangle: angles are π/4, π/4, π/2.
|
||||||
|
// From make_triangle default (v0=(0,0), v1=(1,0), v2=(0,1)):
|
||||||
|
// e01: l=1, λ°=0
|
||||||
|
// e12: l=√2, λ°=log(2)
|
||||||
|
// e02: l=1, λ°=0
|
||||||
|
// Angle at v0 (opposite e12) = π/2.
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(EuclideanFunctional, RightIsoscelesTriangleAnglesCorrect)
|
||||||
|
{
|
||||||
|
const double log2 = std::log(2.0);
|
||||||
|
// lam12 = 0 (v0-v1, length 1), lam23 = log(2) (v1-v2, length √2), lam31 = 0 (v2-v0, length 1)
|
||||||
|
// v1 = v0 in our ordering → remap: l01=1, l12=√2, l20=1
|
||||||
|
// Using euclidean_angles(lam_v1v2, lam_v2v3, lam_v3v1):
|
||||||
|
// v1=(0,0), v2=(1,0), v3=(0,1)
|
||||||
|
// lam12 = log(1²) = 0, lam23 = log(√2 ²) = log2, lam31 = log(1²) = 0
|
||||||
|
auto fa = euclidean_angles(0.0, log2, 0.0);
|
||||||
|
|
||||||
|
ASSERT_TRUE(fa.valid);
|
||||||
|
constexpr double PI = 3.14159265358979323846;
|
||||||
|
// v1=(0,0) is at the right-angle corner (opposite the hypotenuse l23=√2) → α1 = 90°.
|
||||||
|
// v2=(1,0) and v3=(0,1) are the 45° corners (each opposite a leg of length 1).
|
||||||
|
EXPECT_NEAR(fa.alpha1, PI / 2.0, 1e-12); // angle at v1 (opposite l23=√2): 90°
|
||||||
|
EXPECT_NEAR(fa.alpha2, PI / 4.0, 1e-12); // angle at v2 (opposite l31=1): 45°
|
||||||
|
EXPECT_NEAR(fa.alpha3, PI / 4.0, 1e-12); // angle at v3 (opposite l12=1): 45°
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Angle sum = π for any valid Euclidean triangle
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(EuclideanFunctional, AngleSumEqualsPi)
|
||||||
|
{
|
||||||
|
// Scalene triangle with log-lengths (0, 0.5, -0.3).
|
||||||
|
auto fa = euclidean_angles(0.0, 0.5, -0.3);
|
||||||
|
|
||||||
|
ASSERT_TRUE(fa.valid);
|
||||||
|
constexpr double PI = 3.14159265358979323846;
|
||||||
|
EXPECT_NEAR(fa.alpha1 + fa.alpha2 + fa.alpha3, PI, 1e-12);
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Degenerate triangle → valid = false
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(EuclideanFunctional, DegenerateTriangleReturnsFalse)
|
||||||
|
{
|
||||||
|
// l12 = l23 = 1, l31 = 3 → violates triangle inequality.
|
||||||
|
auto fa = euclidean_angles_from_lengths(1.0, 1.0, 3.0);
|
||||||
|
EXPECT_FALSE(fa.valid);
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Gradient check: default right-isosceles triangle, vertex DOFs only
|
||||||
|
//
|
||||||
|
// Mirrors Java testGradient…SingleTriangle.
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(EuclideanFunctional, GradientCheck_TriangleVertex)
|
||||||
|
{
|
||||||
|
auto mesh = make_triangle(); // (0,0)–(1,0)–(0,1)
|
||||||
|
auto maps = setup_euclidean_maps(mesh);
|
||||||
|
compute_euclidean_lambda0_from_mesh(mesh, maps);
|
||||||
|
int n = assign_euclidean_vertex_dof_indices(mesh, maps);
|
||||||
|
|
||||||
|
// Small uniform conformal perturbation.
|
||||||
|
std::vector<double> x(static_cast<std::size_t>(n), -0.1);
|
||||||
|
|
||||||
|
EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps))
|
||||||
|
<< "Gradient check failed on right-isosceles triangle (vertex DOFs)";
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Gradient check: quad strip (2 triangles, 1 interior edge), vertex DOFs only
|
||||||
|
//
|
||||||
|
// Mirrors Java testGradient…QuadStrip / testGradientInExtendedDomain.
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(EuclideanFunctional, GradientCheck_QuadStripVertex)
|
||||||
|
{
|
||||||
|
auto mesh = make_quad_strip();
|
||||||
|
auto maps = setup_euclidean_maps(mesh);
|
||||||
|
compute_euclidean_lambda0_from_mesh(mesh, maps);
|
||||||
|
int n = assign_euclidean_vertex_dof_indices(mesh, maps);
|
||||||
|
|
||||||
|
std::vector<double> x(static_cast<std::size_t>(n), -0.2);
|
||||||
|
|
||||||
|
EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps))
|
||||||
|
<< "Gradient check failed on quad strip (vertex DOFs)";
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Gradient check: regular tetrahedron, vertex DOFs only
|
||||||
|
//
|
||||||
|
// Closed surface (4 faces, 4 vertices, 6 interior edges).
|
||||||
|
// Exercises per-vertex angle-sum accumulation on multiple faces.
|
||||||
|
// Mirrors Java testGradient…Tetrahedron / testGradientWithHyperIdeal…
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(EuclideanFunctional, GradientCheck_TetrahedronVertex)
|
||||||
|
{
|
||||||
|
auto mesh = make_tetrahedron();
|
||||||
|
auto maps = setup_euclidean_maps(mesh);
|
||||||
|
compute_euclidean_lambda0_from_mesh(mesh, maps);
|
||||||
|
int n = assign_euclidean_vertex_dof_indices(mesh, maps);
|
||||||
|
|
||||||
|
std::vector<double> x(static_cast<std::size_t>(n), -0.15);
|
||||||
|
|
||||||
|
EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps))
|
||||||
|
<< "Gradient check failed on regular tetrahedron (vertex DOFs)";
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Gradient check: tetrahedron, all DOFs (vertex + edge)
|
||||||
|
//
|
||||||
|
// Exercises the edge-gradient branch G_e = α_opp⁺ + α_opp⁻ − π.
|
||||||
|
// Mirrors Java testGradientWithHyperellipticCurve.
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(EuclideanFunctional, GradientCheck_TetrahedronAllDofs)
|
||||||
|
{
|
||||||
|
auto mesh = make_tetrahedron();
|
||||||
|
auto maps = setup_euclidean_maps(mesh);
|
||||||
|
compute_euclidean_lambda0_from_mesh(mesh, maps);
|
||||||
|
int n = assign_euclidean_all_dof_indices(mesh, maps);
|
||||||
|
|
||||||
|
// 4 vertex DOFs + 6 edge DOFs = 10 total.
|
||||||
|
std::vector<double> x(static_cast<std::size_t>(n), 0.0);
|
||||||
|
// Set vertex DOFs slightly negative to keep triangles non-degenerate.
|
||||||
|
for (int i = 0; i < 4; ++i)
|
||||||
|
x[static_cast<std::size_t>(i)] = -0.15;
|
||||||
|
|
||||||
|
EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps))
|
||||||
|
<< "Gradient check failed on regular tetrahedron (all DOFs)";
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Angles are finite at a known interior point
|
||||||
|
//
|
||||||
|
// Mirrors Java testFunctionalAtNaNValue: stress-test the angle formula with
|
||||||
|
// large negative conformal factors (compressed triangle) to ensure no NaN/Inf.
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(EuclideanFunctional, AnglesFiniteAtKnownPoint)
|
||||||
|
{
|
||||||
|
auto mesh = make_tetrahedron();
|
||||||
|
auto maps = setup_euclidean_maps(mesh);
|
||||||
|
compute_euclidean_lambda0_from_mesh(mesh, maps);
|
||||||
|
int n = assign_euclidean_vertex_dof_indices(mesh, maps);
|
||||||
|
|
||||||
|
// Very compressed: u_i = -3 (all sides shrunk by exp(-3) ≈ 0.05).
|
||||||
|
// Triangle stays well-formed (equilateral shrinks uniformly).
|
||||||
|
std::vector<double> x(static_cast<std::size_t>(n), -3.0);
|
||||||
|
auto G = euclidean_gradient(mesh, x, maps);
|
||||||
|
|
||||||
|
for (std::size_t i = 0; i < G.size(); ++i) {
|
||||||
|
EXPECT_FALSE(std::isnan(G[i])) << "Gradient component " << i << " is NaN";
|
||||||
|
EXPECT_FALSE(std::isinf(G[i])) << "Gradient component " << i << " is Inf";
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Gradient check: fan of 5 flat triangles, vertex DOFs only
|
||||||
|
//
|
||||||
|
// High-valence central vertex: exercises per-vertex angle accumulation
|
||||||
|
// across 5 incident faces.
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(EuclideanFunctional, GradientCheck_Fan5Vertex)
|
||||||
|
{
|
||||||
|
auto mesh = make_fan(5);
|
||||||
|
auto maps = setup_euclidean_maps(mesh);
|
||||||
|
compute_euclidean_lambda0_from_mesh(mesh, maps);
|
||||||
|
int n = assign_euclidean_vertex_dof_indices(mesh, maps);
|
||||||
|
|
||||||
|
std::vector<double> x(static_cast<std::size_t>(n), -0.05);
|
||||||
|
|
||||||
|
EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps))
|
||||||
|
<< "Gradient check failed on flat fan-5 mesh";
|
||||||
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Gradient check: mixed pinned/variable vertices
|
||||||
|
//
|
||||||
|
// Pins the first vertex (u_v0 = 0 fixed), lets the rest be variable.
|
||||||
|
// Verifies that the gradient accumulator skips pinned vertices correctly.
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(EuclideanFunctional, GradientCheck_MixedPinnedVertices)
|
||||||
|
{
|
||||||
|
auto mesh = make_quad_strip();
|
||||||
|
auto maps = setup_euclidean_maps(mesh);
|
||||||
|
compute_euclidean_lambda0_from_mesh(mesh, maps);
|
||||||
|
|
||||||
|
// Manually pin v0; assign v1, v2, v3 as DOFs 0, 1, 2.
|
||||||
|
auto vit = mesh.vertices().begin();
|
||||||
|
Vertex_index v0 = *vit++;
|
||||||
|
Vertex_index v1 = *vit++;
|
||||||
|
Vertex_index v2 = *vit++;
|
||||||
|
Vertex_index v3 = *vit;
|
||||||
|
|
||||||
|
maps.v_idx[v0] = -1; // pinned
|
||||||
|
maps.v_idx[v1] = 0;
|
||||||
|
maps.v_idx[v2] = 1;
|
||||||
|
maps.v_idx[v3] = 2;
|
||||||
|
|
||||||
|
std::vector<double> x = {-0.1, -0.3, -0.2};
|
||||||
|
|
||||||
|
EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps))
|
||||||
|
<< "Gradient check failed for mixed pinned/variable vertices";
|
||||||
|
}
|
||||||
@@ -1,4 +1,4 @@
|
|||||||
// test_spherical_functional.cpp
|
// test_spherical_functional.cpp (Phase 3c + 3e)
|
||||||
//
|
//
|
||||||
// Phase 3c — SphericalFunctional ported to ConformalMesh.
|
// Phase 3c — SphericalFunctional ported to ConformalMesh.
|
||||||
//
|
//
|
||||||
@@ -244,3 +244,97 @@ TEST(SphericalFunctional, GradientCheck_MixedPinnedVertices)
|
|||||||
EXPECT_TRUE(gradient_check_spherical(mesh, x, maps))
|
EXPECT_TRUE(gradient_check_spherical(mesh, x, maps))
|
||||||
<< "Gradient check failed for mixed pinned/variable vertices";
|
<< "Gradient check failed for mixed pinned/variable vertices";
|
||||||
}
|
}
|
||||||
|
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
// Phase 3e — Gauge-fix for closed spherical surfaces
|
||||||
|
//
|
||||||
|
// On a closed spherical surface, the functional has a gauge mode:
|
||||||
|
// E(u + t·1) is maximised at some t*.
|
||||||
|
// At t*, the sum of all vertex gradients equals zero: Σ G_v = 0.
|
||||||
|
//
|
||||||
|
// Test: start from a point with non-zero ΣG_v, apply the gauge shift,
|
||||||
|
// and verify ΣG_v(x + t*·1) ≈ 0.
|
||||||
|
// ════════════════════════════════════════════════════════════════════════════
|
||||||
|
|
||||||
|
TEST(SphericalFunctional, GaugeFix_SpherTetVertexZerosSumGv)
|
||||||
|
{
|
||||||
|
auto mesh = make_spherical_tetrahedron();
|
||||||
|
auto maps = setup_spherical_maps(mesh);
|
||||||
|
compute_lambda0_from_mesh(mesh, maps);
|
||||||
|
int n = assign_vertex_dof_indices(mesh, maps);
|
||||||
|
|
||||||
|
// Off-gauge starting point: all u_i = -0.5
|
||||||
|
std::vector<double> x(static_cast<std::size_t>(n), -0.5);
|
||||||
|
|
||||||
|
// Compute ΣG_v before gauge shift.
|
||||||
|
{
|
||||||
|
auto G = spherical_gradient(mesh, x, maps);
|
||||||
|
double sum = 0.0;
|
||||||
|
for (auto v : mesh.vertices()) {
|
||||||
|
int iv = maps.v_idx[v];
|
||||||
|
if (iv >= 0) sum += G[static_cast<std::size_t>(iv)];
|
||||||
|
}
|
||||||
|
// At -0.5 the surface is compressed; ΣG_v should be non-zero.
|
||||||
|
EXPECT_NE(sum, 0.0) << "Pre-gauge ΣG_v should be non-zero";
|
||||||
|
}
|
||||||
|
|
||||||
|
// Compute gauge shift and apply.
|
||||||
|
double t = spherical_gauge_shift(mesh, x, maps);
|
||||||
|
std::vector<double> x_fixed = x;
|
||||||
|
for (auto v : mesh.vertices()) {
|
||||||
|
int iv = maps.v_idx[v];
|
||||||
|
if (iv >= 0)
|
||||||
|
x_fixed[static_cast<std::size_t>(iv)] += t;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Verify ΣG_v ≈ 0 at the gauge-fixed point.
|
||||||
|
{
|
||||||
|
auto G = spherical_gradient(mesh, x_fixed, maps);
|
||||||
|
double sum = 0.0;
|
||||||
|
for (auto v : mesh.vertices()) {
|
||||||
|
int iv = maps.v_idx[v];
|
||||||
|
if (iv >= 0) sum += G[static_cast<std::size_t>(iv)];
|
||||||
|
}
|
||||||
|
EXPECT_NEAR(sum, 0.0, 1e-6)
|
||||||
|
<< "After gauge fix, Σ G_v should vanish; t* = " << t;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
TEST(SphericalFunctional, GaugeFix_ApplyInPlace)
|
||||||
|
{
|
||||||
|
auto mesh = make_spherical_tetrahedron();
|
||||||
|
auto maps = setup_spherical_maps(mesh);
|
||||||
|
compute_lambda0_from_mesh(mesh, maps);
|
||||||
|
int n = assign_vertex_dof_indices(mesh, maps);
|
||||||
|
|
||||||
|
// x = -0.3: compressed but inside the valid spherical domain.
|
||||||
|
std::vector<double> x(static_cast<std::size_t>(n), -0.3);
|
||||||
|
|
||||||
|
apply_spherical_gauge(mesh, x, maps);
|
||||||
|
|
||||||
|
// After in-place gauge fix, ΣG_v must be near 0.
|
||||||
|
auto G = spherical_gradient(mesh, x, maps);
|
||||||
|
double sum = 0.0;
|
||||||
|
for (auto v : mesh.vertices()) {
|
||||||
|
int iv = maps.v_idx[v];
|
||||||
|
if (iv >= 0) sum += G[static_cast<std::size_t>(iv)];
|
||||||
|
}
|
||||||
|
EXPECT_NEAR(sum, 0.0, 1e-6)
|
||||||
|
<< "apply_spherical_gauge must drive Σ G_v to zero";
|
||||||
|
}
|
||||||
|
|
||||||
|
TEST(SphericalFunctional, GaugeFix_AlreadyAtGaugeReturnsTNearZero)
|
||||||
|
{
|
||||||
|
// A symmetric, equilateral spherical tetrahedron at x=0 is already
|
||||||
|
// at the gauge maximum (by symmetry, ΣG_v = 0).
|
||||||
|
auto mesh = make_spherical_tetrahedron();
|
||||||
|
auto maps = setup_spherical_maps(mesh);
|
||||||
|
compute_lambda0_from_mesh(mesh, maps);
|
||||||
|
int n = assign_vertex_dof_indices(mesh, maps);
|
||||||
|
|
||||||
|
std::vector<double> x(static_cast<std::size_t>(n), 0.0);
|
||||||
|
double t = spherical_gauge_shift(mesh, x, maps);
|
||||||
|
// Symmetric starting point → t* should be very close to 0.
|
||||||
|
EXPECT_NEAR(t, 0.0, 1e-5)
|
||||||
|
<< "Gauge shift from the symmetric point should be ~0; got " << t;
|
||||||
|
}
|
||||||
|
|||||||
Reference in New Issue
Block a user