833f9e7237
H3: enforce_gauss_bonnet now returns |deficit| (total absolute correction
applied) so callers can detect pathological input without a separate
check_gauss_bonnet call. Both overloads (raw property-map and Maps
template) now return double instead of void. API comment updated.
V5: load_result_xml now implements strict-subset XML rejection instead of
silently mis-reading reformatted-but-valid XML into zeros. The
function documents itself as accepting only the one-element-per-line
format written by save_result_xml. Three strict-subset checks added:
1. <ConformalResult geometry=...> attribute must be on its opening line.
2. <Solver> required attributes must be on the same line.
3. <DOFVector> tag '>' must be on the same line as the tag name.
Non-conforming files throw std::runtime_error immediately.
V6: new helper check_dof_vector_size(x, expected_dofs, context) added to
serialization.hpp. Throws std::runtime_error with a clear message when
the loaded DOF-vector size does not match the mesh's expected DOF count.
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
205 lines
9.9 KiB
C++
205 lines
9.9 KiB
C++
#pragma once
|
||
// Copyright (c) 2024-2026 Tarik Moussa.
|
||
// SPDX-License-Identifier: MIT
|
||
|
||
// gauss_bonnet.hpp
|
||
//
|
||
// 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.
|
||
//
|
||
// ┌─────────────────────────────────────────────────────────────────────────┐
|
||
// │ 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 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
|
||
// (π − Θ_v) term instead, so the identity Σ(2π−Θ_v) = 2π·χ does NOT hold and
|
||
// `check_gauss_bonnet` will (correctly) throw. For open meshes pin the
|
||
// boundary directly and skip the Gauss–Bonnet check (see the CLI's
|
||
// flattening path).
|
||
//
|
||
// API:
|
||
// int euler_characteristic(mesh)
|
||
// int genus(mesh)
|
||
// 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
|
||
// double enforce_gauss_bonnet(mesh, maps) — shifts θ_v by uniform Δ; returns |deficit|
|
||
// (HyperIdealMaps overloads are deleted — see box above)
|
||
|
||
#include "conformal_mesh.hpp"
|
||
#include "euclidean_functional.hpp"
|
||
#include "spherical_functional.hpp"
|
||
#include "hyper_ideal_functional.hpp"
|
||
#include "constants.hpp"
|
||
|
||
#include <stdexcept>
|
||
#include <sstream>
|
||
#include <cmath>
|
||
#include <string>
|
||
|
||
namespace conformallab {
|
||
|
||
// ── Topology helpers ──────────────────────────────────────────────────────────
|
||
|
||
/// Euler characteristic χ = V − E + F.
|
||
/// For closed orientable surfaces: χ = 2 − 2g.
|
||
inline int euler_characteristic(const ConformalMesh& mesh)
|
||
{
|
||
return static_cast<int>(mesh.number_of_vertices())
|
||
- static_cast<int>(mesh.number_of_edges())
|
||
+ static_cast<int>(mesh.number_of_faces());
|
||
}
|
||
|
||
/// Genus of a closed orientable surface: g = (2 − χ) / 2.
|
||
/// Returns 0 for open meshes (boundary present) — callers should check.
|
||
inline int genus(const ConformalMesh& mesh)
|
||
{
|
||
int chi = euler_characteristic(mesh);
|
||
return (2 - chi) / 2;
|
||
}
|
||
|
||
// ── Left-hand side Σ(2π − Θ_v) ─────────────────────────────────────────────
|
||
|
||
/// Sum `Σ_v (2π − Θ_v)` for a raw vertex → angle property map.
|
||
inline double gauss_bonnet_sum(
|
||
const ConformalMesh& mesh,
|
||
const ConformalMesh::Property_map<Vertex_index, double>& theta)
|
||
{
|
||
double s = 0.0;
|
||
for (auto v : mesh.vertices())
|
||
s += TWO_PI - theta[v];
|
||
return s;
|
||
}
|
||
|
||
/// `gauss_bonnet_sum` for the Euclidean-functional property bundle.
|
||
inline double gauss_bonnet_sum(const ConformalMesh& m, const EuclideanMaps& mp)
|
||
{ return gauss_bonnet_sum(m, mp.theta_v); }
|
||
/// `gauss_bonnet_sum` for the Spherical-functional property bundle.
|
||
inline double gauss_bonnet_sum(const ConformalMesh& m, const SphericalMaps& mp)
|
||
{ return gauss_bonnet_sum(m, mp.theta_v); }
|
||
|
||
// gauss_bonnet_sum for 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)`.
|
||
inline double gauss_bonnet_rhs(const ConformalMesh& mesh)
|
||
{
|
||
return TWO_PI * static_cast<double>(euler_characteristic(mesh));
|
||
}
|
||
|
||
// ── Deficit: lhs − rhs (0 = Gauss–Bonnet satisfied) ─────────────────────────
|
||
|
||
/// Gauss-Bonnet deficit `lhs − rhs`; zero iff the identity is satisfied.
|
||
template <typename Maps>
|
||
inline double gauss_bonnet_deficit(const ConformalMesh& mesh, const Maps& maps)
|
||
{
|
||
return gauss_bonnet_sum(mesh, maps) - gauss_bonnet_rhs(mesh);
|
||
}
|
||
|
||
/// Throws `std::runtime_error` if `|lhs − 2π·χ| > tol`.
|
||
/// Overload accepting a precomputed `lhs`.
|
||
inline void check_gauss_bonnet(const ConformalMesh& mesh,
|
||
double lhs,
|
||
double tol = 1e-8)
|
||
{
|
||
double rhs = gauss_bonnet_rhs(mesh);
|
||
double def = lhs - rhs;
|
||
if (std::abs(def) > tol) {
|
||
std::ostringstream msg;
|
||
msg << "Gauss–Bonnet violated:\n"
|
||
<< " Σ(2π−Θ_v) = " << lhs
|
||
<< " expected 2π·χ = " << rhs
|
||
<< " (χ = " << euler_characteristic(mesh)
|
||
<< ", genus = " << genus(mesh) << ")\n"
|
||
<< " deficit = " << def;
|
||
throw std::runtime_error(msg.str());
|
||
}
|
||
}
|
||
|
||
/// Throws `std::runtime_error` if Gauss-Bonnet is violated by more than `tol`.
|
||
template <typename Maps>
|
||
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);
|
||
}
|
||
|
||
// ── enforce_gauss_bonnet — adjust θ_v by uniform Δ ───────────────────────────
|
||
//
|
||
// Adds δ = (lhs − rhs) / V to every θ_v so that Gauss–Bonnet holds exactly.
|
||
// After this call, check_gauss_bonnet() will not throw (up to floating-point).
|
||
// Modifies ALL vertices' θ_v (no v_idx filtering) — the shift is a property
|
||
// of the target angles, independent of which vertices are free DOFs.
|
||
//
|
||
// H3 (test-coverage audit, 2026-06-01): both overloads now return the total
|
||
// absolute correction applied: |Σ(2π−Θ_v) − 2π·χ|. A large value signals
|
||
// that the input angles were far from satisfying Gauss–Bonnet.
|
||
|
||
/// Distribute the Gauss-Bonnet deficit uniformly across all `Θ_v`:
|
||
/// add `δ = (lhs − rhs) / V` to every entry so that the identity holds
|
||
/// exactly afterwards. Overload for a raw property map.
|
||
/// Returns `|lhs − rhs|` (total absolute correction applied).
|
||
inline double enforce_gauss_bonnet(
|
||
ConformalMesh& mesh,
|
||
ConformalMesh::Property_map<Vertex_index, double>& theta)
|
||
{
|
||
double lhs = gauss_bonnet_sum(mesh, theta);
|
||
double rhs = gauss_bonnet_rhs(mesh);
|
||
// Adding δ to every θ_v decreases the sum Σ(2π−θ_v) by V·δ.
|
||
// We need lhs − V·δ = rhs, so δ = (lhs − rhs) / V.
|
||
double delta = (lhs - rhs) / static_cast<double>(mesh.number_of_vertices());
|
||
for (auto v : mesh.vertices())
|
||
theta[v] += delta;
|
||
return std::abs(lhs - rhs);
|
||
}
|
||
|
||
/// 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.
|
||
/// Returns `|lhs − rhs|` (total absolute correction applied; see raw-map overload).
|
||
template <typename Maps>
|
||
inline double enforce_gauss_bonnet(ConformalMesh& mesh, Maps& maps)
|
||
{
|
||
return 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
|