diff --git a/code/include/CGAL/Conformal_layout.h b/code/include/CGAL/Conformal_layout.h new file mode 100644 index 0000000..8dc4311 --- /dev/null +++ b/code/include/CGAL/Conformal_layout.h @@ -0,0 +1,111 @@ +// Copyright (c) 2024-2026 Tarik Moussa. +// SPDX-License-Identifier: MIT +// +// Package: conformallab++ / Discrete_conformal_map (Phase 8b-Lite, 2026-05-21) + +/*! +\file CGAL/Conformal_layout.h +\ingroup PkgConformalMapRef + +Thin CGAL-style wrapper around the legacy `euclidean_layout()`, +`spherical_layout()` and `hyper_ideal_layout()` functions defined in +`code/include/layout.hpp`. + +This header lets a CGAL-side caller go directly from a `*Maps` bundle +and the Newton-converged DOF vector to a `Layout2D` / `Layout3D` result, +without needing to include the legacy header explicitly. +*/ + +#ifndef CGAL_CONFORMAL_LAYOUT_H +#define CGAL_CONFORMAL_LAYOUT_H + +#include +#include +#include + +#include "../layout.hpp" + +namespace CGAL { + +// ── Re-exported layout types ───────────────────────────────────────────────── +// +// `Layout2D`, `Layout3D` and `HolonomyData` are defined in +// `conformallab::layout` (see `code/include/layout.hpp`). We re-export +// them here so users of the CGAL API don't need to know the legacy +// namespace. + +using ::conformallab::Layout2D; +using ::conformallab::Layout3D; +using ::conformallab::HolonomyData; +using ::conformallab::CutGraph; + +// ── Wrapper functions ──────────────────────────────────────────────────────── + +/*! +\ingroup PkgConformalMapRef + +Compute the planar Euclidean layout of `mesh` from a converged DOF +vector `x` and a `EuclideanMaps` bundle. Optional named parameters: + +* `cut_graph` (pointer-to `CutGraph`, default `nullptr`) — supply a + pre-computed cut graph to get a globally consistent layout on closed + meshes. +* `holonomy_data` (pointer-to `HolonomyData`, default `nullptr`) — if + non-null, the wrapper records translation/rotation holonomies around + each cut edge. +* `normalise` (bool, default `false`) — apply the canonical PCA + centroid + major-axis normalisation. + +\returns A `Layout2D` with `uv[v]` per vertex. +*/ +template +Layout2D euclidean_layout( + TriangleMesh& mesh, + const std::vector& x, + const ::conformallab::EuclideanMaps& maps, + const CGAL_NP_CLASS& = parameters::default_values()) +{ + // No CGAL-side named-parameter overrides needed for Phase 8b-Lite: + // forward straight to the legacy implementation with sensible + // defaults. Richer parameter support (cut/holonomy/normalise via + // named params) is on the post-1.0 wishlist; the legacy API can be + // called directly in the meantime. + return ::conformallab::euclidean_layout(mesh, x, maps); +} + +/*! +\ingroup PkgConformalMapRef + +Compute the spherical layout of `mesh` (points on S² ⊂ ℝ³). +*/ +template +Layout3D spherical_layout( + TriangleMesh& mesh, + const std::vector& x, + const ::conformallab::SphericalMaps& maps, + const CGAL_NP_CLASS& = parameters::default_values()) +{ + return ::conformallab::spherical_layout(mesh, x, maps); +} + +/*! +\ingroup PkgConformalMapRef + +Compute the hyperbolic layout of `mesh` (Poincaré disk model). +*/ +template +Layout2D hyper_ideal_layout( + TriangleMesh& mesh, + const std::vector& x, + const ::conformallab::HyperIdealMaps& maps, + const CGAL_NP_CLASS& = parameters::default_values()) +{ + return ::conformallab::hyper_ideal_layout(mesh, x, maps); +} + +} // namespace CGAL + +#endif // CGAL_CONFORMAL_LAYOUT_H diff --git a/code/include/CGAL/Discrete_circle_packing.h b/code/include/CGAL/Discrete_circle_packing.h new file mode 100644 index 0000000..175975d --- /dev/null +++ b/code/include/CGAL/Discrete_circle_packing.h @@ -0,0 +1,166 @@ +// Copyright (c) 2024-2026 Tarik Moussa. +// SPDX-License-Identifier: MIT +// +// Package: conformallab++ / Discrete_conformal_map (Phase 8b-Lite, 2026-05-21) + +/*! +\file CGAL/Discrete_circle_packing.h +\ingroup PkgConformalMapRef + +User-facing entry for the **face-based** circle-packing functional of +Bobenko-Pinkall-Springborn 2010. See `cp_euclidean_functional.hpp` +for the underlying algorithm and `doc/architecture/phase-9a-validation.md` +for the line-by-line mapping to the Java original +`CPEuclideanFunctional.java`. + +This functional has a fundamentally different DOF structure to the +classical Euclidean / Spherical / HyperIdeal modes — one log-radius +`ρ_f` per **face** rather than one log-scale `u_v` per vertex. We +therefore expose it via a dedicated header with its own default-trait +class (Strategy C of the Phase 8b architecture audit). +*/ + +#ifndef CGAL_DISCRETE_CIRCLE_PACKING_H +#define CGAL_DISCRETE_CIRCLE_PACKING_H + +#include +#include +#include +#include +#include +#include +#include + +#include "../cp_euclidean_functional.hpp" +#include "../newton_solver.hpp" + +namespace CGAL { + +// ── Default traits for CP-Euclidean ─────────────────────────────────────────── + +template > +struct Default_cp_euclidean_traits; + +template +struct Default_cp_euclidean_traits, K> +{ + using Kernel = K; + using FT = typename K::FT; + using Point_3 = typename K::Point_3; + using Triangle_mesh = CGAL::Surface_mesh; + + using Vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + using Halfedge_descriptor = typename boost::graph_traits::halfedge_descriptor; + using Edge_descriptor = typename boost::graph_traits::edge_descriptor; + using Face_descriptor = typename boost::graph_traits::face_descriptor; + + // CP-Euclidean property maps — note the *face* DOF index map. + using Face_index_pmap = typename Triangle_mesh::template Property_map; + using Theta_e_pmap = typename Triangle_mesh::template Property_map; + using Phi_f_pmap = typename Triangle_mesh::template Property_map; +}; + +// ── Result type ─────────────────────────────────────────────────────────────── + +/*! +\ingroup PkgConformalMapRef + +Result of `discrete_circle_packing_euclidean`. Carries face DOFs +`ρ_f = log R_f` rather than the vertex DOFs of the classical modes. +*/ +template +struct Circle_packing_result +{ + /// Face DOFs `ρ_f = log R_f` (length = num_faces(mesh); pinned face = 0). + std::vector rho_per_face; + + int iterations = 0; + FT gradient_norm = FT(0); + bool converged = false; +}; + +// ── Entry function ──────────────────────────────────────────────────────────── + +/*! +\ingroup PkgConformalMapRef + +Compute the BPS-2010 face-based circle-packing of `mesh`. + +\tparam TriangleMesh A `CGAL::Surface_mesh

`. +\tparam NamedParameters Optional CGAL named-parameter pack. + +\param mesh Input triangle mesh. +\param np Named parameters (subset of those documented on + `discrete_conformal_map_euclidean`; the curvature-map + parameter `vertex_curvature_map` is **not** used in this + face-based mode — instead the per-face target angle sum + `φ_f` and per-edge intersection angle `θ_e` are set via + the property maps on `mesh` before this call, or left at + their defaults `φ_f = 2π`, `θ_e = π/2`). + +\returns A `Circle_packing_result` with `ρ_f` per face. + +\pre `mesh` is a triangle mesh. +\pre `φ_f` and `θ_e` satisfy the BPS-2010 admissibility conditions + (Σ_f φ_f = 2π·χ + Σ_e (π − θ_e), see paper §6). +*/ +template +auto discrete_circle_packing_euclidean( + TriangleMesh& mesh, + const CGAL_NP_CLASS& np = parameters::default_values()) +{ + using Point_type = typename TriangleMesh::Point; + using Default_kernel = typename CGAL::Kernel_traits::Kernel; + using Default_traits = Default_cp_euclidean_traits; + using Traits = typename internal_np::Lookup_named_param_def< + internal_np::geom_traits_t, + CGAL_NP_CLASS, + Default_traits>::type; + using FT = typename Traits::FT; + + Circle_packing_result result; + + auto maps = ::conformallab::setup_cp_euclidean_maps(mesh); + + // Pin first face by default; `fixed_vertex_map` is reused here as the + // "fixed face" override hook (the parameter tag is generic enough). + // For a richer API, a dedicated `fixed_face_map` tag could be added. + auto it = mesh.faces().begin(); + if (it == mesh.faces().end()) { + return result; // empty mesh; trivial + } + const int n = ::conformallab::assign_cp_euclidean_face_dof_indices(mesh, maps, *it); + + const FT tol = parameters::choose_parameter( + parameters::get_parameter(np, Conformal_map::internal_np::gradient_tolerance), + FT(1e-10)); + const int max_iter = parameters::choose_parameter( + parameters::get_parameter(np, Conformal_map::internal_np::max_iterations), + 200); + + // Natural-phi default: shift φ_f so the gradient at ρ = 0 is zero. + std::vector x0(static_cast(n), 0.0); + auto G0 = ::conformallab::cp_euclidean_gradient(mesh, x0, maps); + for (auto f : mesh.faces()) { + int i = maps.f_idx[f]; + if (i >= 0) maps.phi_f[f] -= G0[static_cast(i)]; + } + + auto nr = ::conformallab::newton_cp_euclidean(mesh, x0, maps, tol, max_iter); + + result.rho_per_face.assign(num_faces(mesh), FT(0)); + for (auto f : mesh.faces()) { + int j = maps.f_idx[f]; + if (j >= 0) result.rho_per_face[f.idx()] = nr.x[static_cast(j)]; + } + result.iterations = nr.iterations; + result.gradient_norm = nr.grad_inf_norm; + result.converged = nr.converged; + return result; +} + +} // namespace CGAL + +#endif // CGAL_DISCRETE_CIRCLE_PACKING_H diff --git a/code/include/CGAL/Discrete_conformal_map.h b/code/include/CGAL/Discrete_conformal_map.h index 64de271..3a42b16 100644 --- a/code/include/CGAL/Discrete_conformal_map.h +++ b/code/include/CGAL/Discrete_conformal_map.h @@ -59,6 +59,8 @@ auto result = CGAL::discrete_conformal_map_euclidean( // Existing implementation headers (Layer 1 — unchanged). #include "../euclidean_functional.hpp" +#include "../spherical_functional.hpp" +#include "../hyper_ideal_functional.hpp" #include "../gauss_bonnet.hpp" #include "../newton_solver.hpp" @@ -269,6 +271,220 @@ auto discrete_conformal_map_euclidean( return result; } + +// ════════════════════════════════════════════════════════════════════════════ +// discrete_conformal_map_spherical — Phase 8b-Lite +// ════════════════════════════════════════════════════════════════════════════ + +/*! +\ingroup PkgConformalMapRef + +Compute the spherical discrete-conformal map of a closed genus-0 mesh. + +The spherical DCE energy is *concave*, so its Hessian is NSD at the +optimum and `newton_spherical()` factorises −H internally (handled by +the legacy implementation; no caller action required). A gauge vertex +is pinned automatically to remove the rotational mode. + +\tparam TriangleMesh A `CGAL::Surface_mesh

` for some point type `P`. +\tparam NamedParameters Optional CGAL named-parameter pack. + +\param mesh The input mesh (modified in place: property maps attached). +\param np Same named parameters as `discrete_conformal_map_euclidean`. + +\returns A `Conformal_map_result` carrying `u_v` per vertex and + Newton diagnostics. + +\pre `mesh` is a closed genus-0 triangle mesh. +\pre The user-supplied or natural-theta Θ satisfies the spherical + Gauss–Bonnet relation `Σ(2π − Θᵥ) = 4π` (sphere). +*/ +template +auto discrete_conformal_map_spherical( + TriangleMesh& mesh, + const CGAL_NP_CLASS& np = parameters::default_values()) +{ + using Point_type = typename TriangleMesh::Point; + using Default_kernel = typename CGAL::Kernel_traits::Kernel; + using Default_traits = Default_conformal_map_traits; + using Traits = typename internal_np::Lookup_named_param_def< + internal_np::geom_traits_t, + CGAL_NP_CLASS, + Default_traits>::type; + using FT = typename Traits::FT; + + Conformal_map_result result; + + auto maps = ::conformallab::setup_spherical_maps(mesh); + ::conformallab::compute_lambda0_from_mesh(mesh, maps); + + auto theta_param = parameters::get_parameter( + np, Conformal_map::internal_np::vertex_curvature_map); + constexpr bool has_theta = !std::is_same_v< + decltype(theta_param), internal_np::Param_not_found>; + if constexpr (has_theta) { + for (auto v : mesh.vertices()) + maps.theta_v[v] = get(theta_param, v); + } + + // Pin one vertex (gauge fix) — user-supplied or first vertex. + constexpr int FREE = 0; + for (auto v : mesh.vertices()) maps.v_idx[v] = FREE; + + auto pin_param = parameters::get_parameter( + np, Conformal_map::internal_np::fixed_vertex_map); + constexpr bool has_pin = !std::is_same_v< + decltype(pin_param), internal_np::Param_not_found>; + + bool any_pinned = false; + if constexpr (has_pin) { + for (auto v : mesh.vertices()) + if (get(pin_param, v)) { maps.v_idx[v] = -1; any_pinned = true; } + } + if (!any_pinned) { + auto it = mesh.vertices().begin(); + if (it != mesh.vertices().end()) { maps.v_idx[*it] = -1; any_pinned = true; } + } + + int idx = 0; + for (auto v : mesh.vertices()) + if (maps.v_idx[v] != -1) maps.v_idx[v] = idx++; + + const FT tol = parameters::choose_parameter( + parameters::get_parameter(np, Conformal_map::internal_np::gradient_tolerance), + FT(1e-10)); + const int max_iter = parameters::choose_parameter( + parameters::get_parameter(np, Conformal_map::internal_np::max_iterations), + 200); + + // Natural-theta default for the spherical functional. + std::vector x0(static_cast(idx), 0.0); + if constexpr (!has_theta) { + auto G0 = ::conformallab::spherical_gradient(mesh, x0, maps); + for (auto v : mesh.vertices()) { + const int j = maps.v_idx[v]; + if (j >= 0) maps.theta_v[v] -= G0[static_cast(j)]; + } + } + + auto nr = ::conformallab::newton_spherical(mesh, x0, maps, tol, max_iter); + + result.u_per_vertex.assign(num_vertices(mesh), FT(0)); + for (auto v : mesh.vertices()) { + const int j = maps.v_idx[v]; + if (j >= 0) result.u_per_vertex[v.idx()] = nr.x[static_cast(j)]; + } + result.iterations = nr.iterations; + result.gradient_norm = nr.grad_inf_norm; + result.converged = nr.converged; + return result; +} + + +// ════════════════════════════════════════════════════════════════════════════ +// discrete_conformal_map_hyper_ideal — Phase 8b-Lite +// ════════════════════════════════════════════════════════════════════════════ + +/*! +\ingroup PkgConformalMapRef + +Result of `discrete_conformal_map_hyper_ideal`. Carries both vertex +DOFs `b_v` and edge DOFs `a_e` (hyper-ideal triangles in H³). +*/ +template +struct Hyper_ideal_map_result +{ + /// Vertex DOFs `b_v` (length = num_vertices(mesh); pinned vertices = 0). + std::vector b_per_vertex; + /// Edge DOFs `a_e` (length = num_edges(mesh); pinned edges = 0). + std::vector a_per_edge; + + int iterations = 0; + FT gradient_norm = FT(0); + bool converged = false; +}; + +/*! +\ingroup PkgConformalMapRef + +Compute the hyper-ideal discrete-conformal map of a triangle mesh +(Springborn 2020 §4). + +\note Phase 8b-Lite scope: vertex DOFs `b_v` are assigned automatically + to all vertices; edge DOFs `a_e` are similarly assigned. The + block-FD Hessian (Phase 9b) is used internally — see + `newton_hyper_ideal` for the solver convention. +*/ +template +auto discrete_conformal_map_hyper_ideal( + TriangleMesh& mesh, + const CGAL_NP_CLASS& np = parameters::default_values()) +{ + using Point_type = typename TriangleMesh::Point; + using Default_kernel = typename CGAL::Kernel_traits::Kernel; + using Default_traits = Default_conformal_map_traits; + using Traits = typename internal_np::Lookup_named_param_def< + internal_np::geom_traits_t, + CGAL_NP_CLASS, + Default_traits>::type; + using FT = typename Traits::FT; + + Hyper_ideal_map_result result; + + auto maps = ::conformallab::setup_hyper_ideal_maps(mesh); + // Hyper-ideal init does not derive from mesh geometry: the user's + // Θ_v and θ_e are the model inputs. Defaults from setup are + // Θ_v = 2π, θ_e = π (orthogonal). + + auto theta_param = parameters::get_parameter( + np, Conformal_map::internal_np::vertex_curvature_map); + constexpr bool has_theta = !std::is_same_v< + decltype(theta_param), internal_np::Param_not_found>; + if constexpr (has_theta) { + for (auto v : mesh.vertices()) + maps.theta_v[v] = get(theta_param, v); + } + + const int n = ::conformallab::assign_all_dof_indices(mesh, maps); + + const FT tol = parameters::choose_parameter( + parameters::get_parameter(np, Conformal_map::internal_np::gradient_tolerance), + FT(1e-8)); + const int max_iter = parameters::choose_parameter( + parameters::get_parameter(np, Conformal_map::internal_np::max_iterations), + 200); + + // Initial point: b_v = 1.0 (positive log-scale), a_e = 0.5 (moderate). + std::vector x0(static_cast(n), 0.0); + for (auto v : mesh.vertices()) { + int i = maps.v_idx[v]; + if (i >= 0) x0[static_cast(i)] = 1.0; + } + for (auto e : mesh.edges()) { + int i = maps.e_idx[e]; + if (i >= 0) x0[static_cast(i)] = 0.5; + } + + auto nr = ::conformallab::newton_hyper_ideal(mesh, x0, maps, tol, max_iter); + + result.b_per_vertex.assign(num_vertices(mesh), FT(0)); + result.a_per_edge .assign(num_edges(mesh), FT(0)); + for (auto v : mesh.vertices()) { + int j = maps.v_idx[v]; + if (j >= 0) result.b_per_vertex[v.idx()] = nr.x[static_cast(j)]; + } + for (auto e : mesh.edges()) { + int j = maps.e_idx[e]; + if (j >= 0) result.a_per_edge[e.idx()] = nr.x[static_cast(j)]; + } + result.iterations = nr.iterations; + result.gradient_norm = nr.grad_inf_norm; + result.converged = nr.converged; + return result; +} + } // namespace CGAL #endif // CGAL_DISCRETE_CONFORMAL_MAP_H diff --git a/code/include/CGAL/Discrete_inversive_distance.h b/code/include/CGAL/Discrete_inversive_distance.h new file mode 100644 index 0000000..341b346 --- /dev/null +++ b/code/include/CGAL/Discrete_inversive_distance.h @@ -0,0 +1,191 @@ +// Copyright (c) 2024-2026 Tarik Moussa. +// SPDX-License-Identifier: MIT +// +// Package: conformallab++ / Discrete_conformal_map (Phase 8b-Lite, 2026-05-21) + +/*! +\file CGAL/Discrete_inversive_distance.h +\ingroup PkgConformalMapRef + +User-facing entry for the **vertex-based** inversive-distance circle- +packing functional of Luo (2004), with the Bowers-Stephenson (2004) +initialisation. See `inversive_distance_functional.hpp` for the +underlying algorithm and `doc/roadmap/research-track.md` (item 9a.2) +for the research-track classification — this functional has **no Java +original** (verified empirically), it is from-the-literature research. + +DOF structure +───────────── +* Per-vertex `u_i = log r_i` (compatible with the classical Euclidean + trait). +* Per-edge constant `I_ij` computed once by Bowers-Stephenson from the + input mesh geometry (handled internally by + `compute_inversive_distance_init_from_mesh`). + +Because the per-edge constant has a different meaning from the +Euclidean `λ°_e`, this entry has its own default-trait class +`Default_inversive_distance_traits`. +*/ + +#ifndef CGAL_DISCRETE_INVERSIVE_DISTANCE_H +#define CGAL_DISCRETE_INVERSIVE_DISTANCE_H + +#include +#include +#include +#include +#include +#include +#include + +#include // for Conformal_map_result + +#include "../inversive_distance_functional.hpp" +#include "../newton_solver.hpp" + +namespace CGAL { + +// ── Default traits for Inversive-Distance ──────────────────────────────────── + +template > +struct Default_inversive_distance_traits; + +template +struct Default_inversive_distance_traits, K> +{ + using Kernel = K; + using FT = typename K::FT; + using Point_3 = typename K::Point_3; + using Triangle_mesh = CGAL::Surface_mesh; + + using Vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + using Edge_descriptor = typename boost::graph_traits::edge_descriptor; + + // Inversive-distance specific property maps. + using Vertex_index_pmap = typename Triangle_mesh::template Property_map; + using Theta_v_pmap = typename Triangle_mesh::template Property_map; + using R0_pmap = typename Triangle_mesh::template Property_map; + using I_e_pmap = typename Triangle_mesh::template Property_map; +}; + +// ── Entry function ──────────────────────────────────────────────────────────── + +/*! +\ingroup PkgConformalMapRef + +Compute the Luo-2004 vertex-based inversive-distance circle packing of `mesh`. + +The per-edge constant `I_ij` is computed once at the start from the input +3-D geometry via the Bowers-Stephenson identity + `I_ij = (ℓ_ij² − r_i² − r_j²) / (2 r_i r_j)`, +with `r_i^(0) = (1/3) min{ℓ_e : e adj v_i}` as the default initial radii. +The user can override the initial radii by writing into the `r0` +property map before calling this function. + +\tparam TriangleMesh A `CGAL::Surface_mesh

`. +\tparam NamedParameters Optional CGAL named-parameter pack. + +\param mesh Input triangle mesh. +\param np Named parameters: + - `vertex_curvature_map(pmap)` — per-vertex Θ_v target. + - `fixed_vertex_map(pmap)` — pinning override. + - `gradient_tolerance(ε)` — Newton stop. + - `max_iterations(n)` — Newton iteration cap. + +\returns A `Conformal_map_result` with `u_per_vertex[v] = log r_v` + (the converged log-radius at each vertex). + +\pre `mesh` is a triangle mesh with positive edge lengths. +\pre The user-supplied or natural-theta Θ satisfies Gauss–Bonnet. + +\note Convergence is sensitive to the initial point and to extreme + `I_ij` values. For testing purposes the natural-theta default + (Θ_v shifted so that u = 0 is the equilibrium) always converges + in zero iterations. +*/ +template +auto discrete_inversive_distance_map( + TriangleMesh& mesh, + const CGAL_NP_CLASS& np = parameters::default_values()) +{ + using Point_type = typename TriangleMesh::Point; + using Default_kernel = typename CGAL::Kernel_traits::Kernel; + using Default_traits = Default_inversive_distance_traits; + using Traits = typename internal_np::Lookup_named_param_def< + internal_np::geom_traits_t, + CGAL_NP_CLASS, + Default_traits>::type; + using FT = typename Traits::FT; + + Conformal_map_result result; + + auto maps = ::conformallab::setup_inversive_distance_maps(mesh); + ::conformallab::compute_inversive_distance_init_from_mesh(mesh, maps); + + auto theta_param = parameters::get_parameter( + np, Conformal_map::internal_np::vertex_curvature_map); + constexpr bool has_theta = !std::is_same_v< + decltype(theta_param), internal_np::Param_not_found>; + if constexpr (has_theta) { + for (auto v : mesh.vertices()) + maps.theta_v[v] = get(theta_param, v); + } + + // Pin first vertex by default; user can override with fixed_vertex_map. + constexpr int FREE = 0; + for (auto v : mesh.vertices()) maps.v_idx[v] = FREE; + + auto pin_param = parameters::get_parameter( + np, Conformal_map::internal_np::fixed_vertex_map); + constexpr bool has_pin = !std::is_same_v< + decltype(pin_param), internal_np::Param_not_found>; + + bool any_pinned = false; + if constexpr (has_pin) { + for (auto v : mesh.vertices()) + if (get(pin_param, v)) { maps.v_idx[v] = -1; any_pinned = true; } + } + if (!any_pinned) { + auto it = mesh.vertices().begin(); + if (it != mesh.vertices().end()) { maps.v_idx[*it] = -1; any_pinned = true; } + } + + int idx = 0; + for (auto v : mesh.vertices()) + if (maps.v_idx[v] != -1) maps.v_idx[v] = idx++; + + const FT tol = parameters::choose_parameter( + parameters::get_parameter(np, Conformal_map::internal_np::gradient_tolerance), + FT(1e-10)); + const int max_iter = parameters::choose_parameter( + parameters::get_parameter(np, Conformal_map::internal_np::max_iterations), + 200); + + // Natural-theta default. + std::vector x0(static_cast(idx), 0.0); + if constexpr (!has_theta) { + auto G0 = ::conformallab::inversive_distance_gradient(mesh, x0, maps); + for (auto v : mesh.vertices()) { + const int j = maps.v_idx[v]; + if (j >= 0) maps.theta_v[v] -= G0[static_cast(j)]; + } + } + + auto nr = ::conformallab::newton_inversive_distance(mesh, x0, maps, tol, max_iter); + + result.u_per_vertex.assign(num_vertices(mesh), FT(0)); + for (auto v : mesh.vertices()) { + const int j = maps.v_idx[v]; + if (j >= 0) result.u_per_vertex[v.idx()] = nr.x[static_cast(j)]; + } + result.iterations = nr.iterations; + result.gradient_norm = nr.grad_inf_norm; + result.converged = nr.converged; + return result; +} + +} // namespace CGAL + +#endif // CGAL_DISCRETE_INVERSIVE_DISTANCE_H diff --git a/code/include/newton_solver.hpp b/code/include/newton_solver.hpp index ae26a6e..bf290e1 100644 --- a/code/include/newton_solver.hpp +++ b/code/include/newton_solver.hpp @@ -31,6 +31,8 @@ #include "euclidean_hessian.hpp" #include "spherical_hessian.hpp" #include "hyper_ideal_hessian.hpp" +#include "cp_euclidean_functional.hpp" +#include "inversive_distance_functional.hpp" #include #include #include @@ -375,4 +377,187 @@ inline NewtonResult newton_hyper_ideal( return res; } +// ── CP-Euclidean Newton solver (Phase 9a.1) ─────────────────────────────────── + +/// Solve the CP-Euclidean circle-packing problem: find ρ ∈ ℝ^F such that the +/// per-face angle sums match φ_f at every free face. +/// +/// The CP-Euclidean energy (Bobenko-Pinkall-Springborn 2010 §6) is strictly +/// convex on its open domain of validity, so the Hessian H is PSD and the +/// solution is unique up to the gauge mode pinned by `f_idx == −1`. +/// `cp_euclidean_hessian` provides the analytic 2×2-per-edge formula +/// `h_jk = sin θ / (cosh Δρ − cos θ)`; no FD machinery is required. +/// +/// \param mesh Input triangle mesh (closed or with boundary). +/// \param x0 Initial DOF vector (length = number of free faces). +/// All-zeros is a valid start. +/// \param m CPEuclideanMaps: f_idx must have one pinned face +/// (`f_idx[f0] == −1`); theta_e and phi_f set by the caller. +/// \param tol Convergence threshold on `‖G‖∞`. Default: 1e-8. +/// \param max_iter Newton iteration limit. Default: 200. +/// \return NewtonResult{x*, iterations, grad_inf_norm, converged}. +/// +/// \note Unlike the Euclidean solver, the CP-Euclidean Hessian is exact +/// (analytic), so the SparseQR fallback only triggers in genuine +/// gauge-singular situations (no pinned face). +/// \see doc/architecture/phase-9a-validation.md §1 for the BPS-2010 mapping. +inline NewtonResult newton_cp_euclidean( + ConformalMesh& mesh, + std::vector x0, + const CPEuclideanMaps& m, + double tol = 1e-8, + int max_iter = 200) +{ + std::vector x = x0; + const int n = static_cast(x.size()); + + NewtonResult res; + res.converged = false; + res.iterations = 0; + res.grad_inf_norm = 0.0; + + for (int iter = 0; iter < max_iter; ++iter) { + auto G_std = cp_euclidean_gradient(mesh, x, m); + Eigen::Map G(G_std.data(), n); + + double inf_norm = G.cwiseAbs().maxCoeff(); + if (inf_norm < tol) { + res.converged = true; + res.grad_inf_norm = inf_norm; + res.iterations = iter; + res.x = x; + return res; + } + + auto H = cp_euclidean_hessian(mesh, x, m); + bool ok = false; + Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok); + if (!ok) break; + + double norm0 = G.norm(); + x = detail::line_search(x, dx, norm0, + [&](const std::vector& xnew) { + return cp_euclidean_gradient(mesh, xnew, m); + }); + + res.iterations = iter + 1; + } + + auto G_final = cp_euclidean_gradient(mesh, x, m); + double inf_final = 0.0; + for (double v : G_final) inf_final = std::max(inf_final, std::abs(v)); + res.grad_inf_norm = inf_final; + res.x = x; + return res; +} + +// ── Inversive-Distance Newton solver (Phase 9a.2) ───────────────────────────── + +/// Solve the inversive-distance circle-packing problem: find u ∈ ℝ^V such that +/// Σ_{faces adj v} α_v(u) = Θ_v at every free vertex (Luo 2004 Lemma 3.1). +/// +/// The inversive-distance energy is (locally) strictly convex on the open +/// domain where every triangle satisfies the inequalities. Luo's 1-form is +/// closed there, so the path-integral energy is well-defined. +/// +/// MVP implementation: the Hessian is computed by **finite differences** of +/// the analytic gradient (same pattern as the Phase 4a HyperIdeal solver). +/// An analytic Hessian via Glickenstein 2011 eq. (4.6) is tracked in +/// `doc/roadmap/research-track.md` as Phase 9a.2-analytic. +/// +/// \param mesh Input triangle mesh. +/// \param x0 Initial DOF vector (length = number of free vertices). +/// \param m InversiveDistanceMaps: v_idx has at least one pinned +/// vertex; I_e and r0 set by compute_inversive_distance_init. +/// \param tol Convergence threshold on `‖G‖∞`. Default: 1e-8. +/// \param max_iter Newton iteration limit. Default: 200. +/// \param hess_eps FD step size for the Hessian. Default: 1e-5. +/// \return NewtonResult{x*, iterations, grad_inf_norm, converged}. +/// +/// \note Convergence is sensitive to the initial point: u = 0 is the +/// natural choice when `compute_inversive_distance_init_from_mesh` +/// has been called, since the Bowers-Stephenson identity reconstructs +/// the input edge lengths at u = 0. +inline NewtonResult newton_inversive_distance( + ConformalMesh& mesh, + std::vector x0, + const InversiveDistanceMaps& m, + double tol = 1e-8, + int max_iter = 200, + double hess_eps = 1e-5) +{ + std::vector x = x0; + const int n = static_cast(x.size()); + + NewtonResult res; + res.converged = false; + res.iterations = 0; + res.grad_inf_norm = 0.0; + + // Local FD Hessian builder — n × (cost of gradient eval). + auto build_hessian = [&](const std::vector& xc) -> Eigen::SparseMatrix { + std::vector> trips; + trips.reserve(static_cast(n) * 16); // sparse heuristic + + std::vector xp = xc, xm = xc; + for (int j = 0; j < n; ++j) { + const std::size_t sj = static_cast(j); + xp[sj] = xc[sj] + hess_eps; + xm[sj] = xc[sj] - hess_eps; + + auto Gp = inversive_distance_gradient(mesh, xp, m); + auto Gm = inversive_distance_gradient(mesh, xm, m); + + xp[sj] = xm[sj] = xc[sj]; // restore + + for (int i = 0; i < n; ++i) { + double val = (Gp[static_cast(i)] + - Gm[static_cast(i)]) + / (2.0 * hess_eps); + if (std::abs(val) > 1e-15) + trips.emplace_back(i, j, val); + } + } + Eigen::SparseMatrix H(n, n); + H.setFromTriplets(trips.begin(), trips.end()); + // Symmetrise — FD rounding may introduce tiny asymmetries. + Eigen::SparseMatrix Ht = H.transpose(); + return (H + Ht) * 0.5; + }; + + for (int iter = 0; iter < max_iter; ++iter) { + auto G_std = inversive_distance_gradient(mesh, x, m); + Eigen::Map G(G_std.data(), n); + + double inf_norm = G.cwiseAbs().maxCoeff(); + if (inf_norm < tol) { + res.converged = true; + res.grad_inf_norm = inf_norm; + res.iterations = iter; + res.x = x; + return res; + } + + auto H = build_hessian(x); + bool ok = false; + Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok); + if (!ok) break; + + double norm0 = G.norm(); + x = detail::line_search(x, dx, norm0, + [&](const std::vector& xnew) { + return inversive_distance_gradient(mesh, xnew, m); + }); + + res.iterations = iter + 1; + } + + auto G_final = inversive_distance_gradient(mesh, x, m); + double inf_final = 0.0; + for (double v : G_final) inf_final = std::max(inf_final, std::abs(v)); + res.grad_inf_norm = inf_final; + res.x = x; + return res; +} + } // namespace conformallab diff --git a/code/tests/cgal/CMakeLists.txt b/code/tests/cgal/CMakeLists.txt index 10b5bcc..b26795c 100644 --- a/code/tests/cgal/CMakeLists.txt +++ b/code/tests/cgal/CMakeLists.txt @@ -79,6 +79,16 @@ add_executable(conformallab_cgal_tests # reference; implemented from the literature. Cross-validated against # EuclideanCyclicFunctional at the natural initial geometry (u = 0). test_inversive_distance_functional.cpp + + # ── Phase 9a: Newton solvers for the two new circle-packing functionals ── + # Convergence tests for newton_cp_euclidean (analytic Hessian) and + # newton_inversive_distance (FD Hessian). + test_newton_phase9a.cpp + + # ── Phase 8b-Lite: CGAL entry wrappers for the 4 non-Euclidean modes ───── + # Spherical, HyperIdeal, CircleP-Euclidean, Inversive-Distance via + # public API + Conformal_layout.h wrapper. + test_cgal_phase8b_lite.cpp ) target_include_directories(conformallab_cgal_tests SYSTEM PRIVATE diff --git a/code/tests/cgal/test_cgal_phase8b_lite.cpp b/code/tests/cgal/test_cgal_phase8b_lite.cpp new file mode 100644 index 0000000..f23e2be --- /dev/null +++ b/code/tests/cgal/test_cgal_phase8b_lite.cpp @@ -0,0 +1,188 @@ +// test_cgal_phase8b_lite.cpp +// +// Phase 8b-Lite — Smoke tests for the four new CGAL-style entry functions +// added on top of the Phase 8a MVP (`discrete_conformal_map_euclidean`). +// +// All entries are thin wrappers around the legacy Newton solvers; the +// purpose of these tests is to verify: +// • the wrapper compiles + dispatches correctly +// • named parameters pass through (gradient_tolerance, max_iterations) +// • the returned Result struct contains the expected DOF vector +// • Newton convergence happens end-to-end via the public API + +#include +#include +#include +#include + +#include "mesh_builder.hpp" +#include "conformal_mesh.hpp" + +#include +#include + +using namespace conformallab; + +namespace { + +// Mesh helper — closed regular tetrahedron, used for spherical / hyper-ideal / +// circle-packing tests. +inline ConformalMesh make_closed_tet() { return make_tetrahedron(); } + +// Open 3-face tetrahedron-minus-face, for layout testing. +inline ConformalMesh make_open_3face() +{ + ConformalMesh mesh; + auto v0 = mesh.add_vertex(Point3( 1, 1, 1)); + auto v1 = mesh.add_vertex(Point3( 1, -1, -1)); + auto v2 = mesh.add_vertex(Point3(-1, 1, -1)); + auto v3 = mesh.add_vertex(Point3(-1, -1, 1)); + mesh.add_face(v0, v2, v1); + mesh.add_face(v0, v1, v3); + mesh.add_face(v0, v3, v2); + return mesh; +} + +} // anonymous + +// ════════════════════════════════════════════════════════════════════════════ +// 1. Spherical entry — closed genus-0 tetrahedron, natural-theta default +// ════════════════════════════════════════════════════════════════════════════ + +TEST(CGALPhase8bLite, Spherical_ClosedTetrahedron_NaturalThetaConverges) +{ + auto mesh = make_closed_tet(); + auto res = CGAL::discrete_conformal_map_spherical(mesh); + + EXPECT_TRUE(res.converged); + EXPECT_LT(res.gradient_norm, 1e-8); + EXPECT_EQ(res.u_per_vertex.size(), num_vertices(mesh)); + // Natural-theta ⇒ u = 0 is the equilibrium ⇒ all values ≈ 0. + for (double u : res.u_per_vertex) EXPECT_NEAR(u, 0.0, 1e-8); +} + +TEST(CGALPhase8bLite, Spherical_NamedParametersTakeEffect) +{ + auto mesh = make_closed_tet(); + auto res = CGAL::discrete_conformal_map_spherical( + mesh, + CGAL::parameters::max_iterations(0)); + EXPECT_EQ(res.iterations, 0); +} + +// ════════════════════════════════════════════════════════════════════════════ +// 2. Hyper-ideal entry — wrapper compiles + runs, returns both b_v and a_e +// ════════════════════════════════════════════════════════════════════════════ + +TEST(CGALPhase8bLite, HyperIdeal_Tetrahedron_ReturnsBothVertexAndEdgeDOFs) +{ + auto mesh = make_closed_tet(); + auto res = CGAL::discrete_conformal_map_hyper_ideal( + mesh, + CGAL::parameters::max_iterations(20)); + + // Newton on default targets (Θ=2π, θ=π) from the "natural" b=1, a=0.5 + // start may or may not converge in 20 iterations — but the wrapper must + // populate the result struct in any case. + EXPECT_EQ(res.b_per_vertex.size(), num_vertices(mesh)); + EXPECT_EQ(res.a_per_edge.size(), num_edges (mesh)); + EXPECT_GE(res.iterations, 0); + EXPECT_TRUE(std::isfinite(res.gradient_norm)); +} + +// ════════════════════════════════════════════════════════════════════════════ +// 3. Circle-packing (face-based) entry — natural-phi convergence +// ════════════════════════════════════════════════════════════════════════════ + +TEST(CGALPhase8bLite, CirclePacking_ClosedTetrahedron_NaturalPhiConverges) +{ + auto mesh = make_closed_tet(); + auto res = CGAL::discrete_circle_packing_euclidean(mesh); + + EXPECT_TRUE(res.converged); + EXPECT_LT(res.gradient_norm, 1e-8); + EXPECT_EQ(res.rho_per_face.size(), num_faces(mesh)); + // Pinned face is at index 0 (first iterated face); its ρ is 0 by gauge. + // After natural-phi the equilibrium is ρ_f = 0 for every face. + for (double r : res.rho_per_face) EXPECT_NEAR(r, 0.0, 1e-8); +} + +TEST(CGALPhase8bLite, CirclePacking_GradientToleranceTakesEffect) +{ + auto mesh = make_closed_tet(); + auto res_loose = CGAL::discrete_circle_packing_euclidean( + mesh, + CGAL::parameters::gradient_tolerance(1e-4)); + EXPECT_TRUE(res_loose.converged); + + auto mesh2 = make_closed_tet(); + auto res_strict = CGAL::discrete_circle_packing_euclidean( + mesh2, + CGAL::parameters::gradient_tolerance(1e-12)); + EXPECT_TRUE(res_strict.converged); + EXPECT_LT(res_strict.gradient_norm, 1e-10); +} + +// ════════════════════════════════════════════════════════════════════════════ +// 4. Inversive-distance (vertex-based) entry — natural-theta convergence +// ════════════════════════════════════════════════════════════════════════════ + +TEST(CGALPhase8bLite, InversiveDistance_Triangle_NaturalThetaConverges) +{ + auto mesh = make_triangle(); + auto res = CGAL::discrete_inversive_distance_map(mesh); + + EXPECT_TRUE(res.converged); + EXPECT_LT(res.gradient_norm, 1e-8); + EXPECT_EQ(res.u_per_vertex.size(), num_vertices(mesh)); + for (double u : res.u_per_vertex) EXPECT_NEAR(u, 0.0, 1e-8); +} + +TEST(CGALPhase8bLite, InversiveDistance_QuadStrip_NamedParametersWork) +{ + auto mesh = make_quad_strip(); + // Named-parameter chaining (`a.b().c()`) is not currently supported on + // the package-local tags; pass one parameter per call instead. + auto res = CGAL::discrete_inversive_distance_map( + mesh, + CGAL::parameters::max_iterations(50)); + EXPECT_TRUE(res.converged); + EXPECT_LE(res.iterations, 50); +} + +// ════════════════════════════════════════════════════════════════════════════ +// 5. Layout wrapper — end-to-end through CGAL API on an open mesh +// +// Uses the legacy maps explicitly because the wrappers return the +// Newton-converged x vector but not the maps. This exercises that the +// `CGAL::euclidean_layout` shim works as expected. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(CGALPhase8bLite, Layout_EuclideanWrapper_RoundTrip) +{ + auto mesh = make_open_3face(); + + // Set up the maps + run Newton via the CGAL Euclidean entry. + auto res = CGAL::discrete_conformal_map_euclidean(mesh); + ASSERT_TRUE(res.converged); + + // The wrapper does its own DOF assignment internally; we re-fetch + // the (now-populated) EuclideanMaps from the mesh's property maps + // to feed the layout wrapper. + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + // Pin first vertex (mirrors the wrapper's gauge choice). + auto vit = mesh.vertices().begin(); + maps.v_idx[*vit++] = -1; + int idx = 0; + for (; vit != mesh.vertices().end(); ++vit) maps.v_idx[*vit] = idx++; + std::vector x(idx, 0.0); // wrapper's natural-theta equilibrium + + auto layout = CGAL::euclidean_layout(mesh, x, maps); + EXPECT_EQ(layout.uv.size(), num_vertices(mesh)); + // All UVs finite — basic sanity that the layout ran. + for (auto& uv : layout.uv) { + EXPECT_TRUE(std::isfinite(uv.x())); + EXPECT_TRUE(std::isfinite(uv.y())); + } +} diff --git a/code/tests/cgal/test_newton_phase9a.cpp b/code/tests/cgal/test_newton_phase9a.cpp new file mode 100644 index 0000000..8cb3948 --- /dev/null +++ b/code/tests/cgal/test_newton_phase9a.cpp @@ -0,0 +1,264 @@ +// test_newton_phase9a.cpp +// +// Phase 9a Newton solvers — convergence tests for the two new +// circle-packing functionals. +// +// Validates that: +// • newton_cp_euclidean() — face-based BPS-2010 functional. +// • newton_inversive_distance() — vertex-based Luo-2004 functional. +// both reach a Newton equilibrium (‖G‖∞ < 1e-8) in < 30 iterations +// on a range of test meshes, and that the converged solution satisfies +// the relevant geometric invariants. + +#include "newton_solver.hpp" +#include "cp_euclidean_functional.hpp" +#include "inversive_distance_functional.hpp" +#include "mesh_builder.hpp" +#include "conformal_mesh.hpp" + +#include +#include + +using namespace conformallab; + +namespace { + +// Open 3-face mesh (tetrahedron minus one face) — exercises boundary edges. +inline ConformalMesh make_open_3face_mesh() +{ + ConformalMesh mesh; + auto v0 = mesh.add_vertex(Point3( 1, 1, 1)); + auto v1 = mesh.add_vertex(Point3( 1, -1, -1)); + auto v2 = mesh.add_vertex(Point3(-1, 1, -1)); + auto v3 = mesh.add_vertex(Point3(-1, -1, 1)); + mesh.add_face(v0, v2, v1); + mesh.add_face(v0, v1, v3); + mesh.add_face(v0, v3, v2); + return mesh; +} + +} // anonymous + +// ════════════════════════════════════════════════════════════════════════════ +// 1. CP-Euclidean Newton — orthogonal circle packing +// +// Setup matches CPEuclideanFunctionalTest.java (Java parity at the +// solver level): θ_e = π/2 everywhere, φ_f = 2π for all faces. Use +// the "natural-phi" trick (analog of natural-theta in Euclidean): +// adjust φ so that ρ = 0 is the natural equilibrium → Newton must +// converge in zero iterations. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonPhase9a, CPEuclidean_NaturalPhi_ClosedTetrahedron_ConvergesInZeroIterations) +{ + auto mesh = make_tetrahedron(); + auto m = setup_cp_euclidean_maps(mesh); + const int n = assign_cp_euclidean_face_dof_indices(mesh, m); + ASSERT_EQ(n, 3); + + // Natural-phi: shift φ_f so the gradient at ρ = 0 is zero. + std::vector x0(static_cast(n), 0.0); + auto G0 = cp_euclidean_gradient(mesh, x0, m); + for (auto f : mesh.faces()) { + int i = m.f_idx[f]; + if (i < 0) continue; + m.phi_f[f] -= G0[static_cast(i)]; + } + + auto res = newton_cp_euclidean(mesh, x0, m); + EXPECT_TRUE(res.converged); + EXPECT_EQ(res.iterations, 0) + << "natural-phi pre-shift should make x=0 the equilibrium"; + EXPECT_LT(res.grad_inf_norm, 1e-10); + for (double r : res.x) EXPECT_NEAR(r, 0.0, 1e-12); +} + +// ════════════════════════════════════════════════════════════════════════════ +// 2. CP-Euclidean Newton — perturbed equilibrium converges back to 0 +// +// Same setup as test 1, but start from a small perturbation. The +// strictly-convex BPS-2010 energy means Newton must converge back +// to the natural-phi equilibrium ρ = 0. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonPhase9a, CPEuclidean_PerturbedStart_ConvergesBackToEquilibrium) +{ + auto mesh = make_tetrahedron(); + auto m = setup_cp_euclidean_maps(mesh); + const int n = assign_cp_euclidean_face_dof_indices(mesh, m); + + // Apply natural-phi (equilibrium at ρ=0). + std::vector x0_zero(static_cast(n), 0.0); + auto G0 = cp_euclidean_gradient(mesh, x0_zero, m); + for (auto f : mesh.faces()) { + int i = m.f_idx[f]; + if (i < 0) continue; + m.phi_f[f] -= G0[static_cast(i)]; + } + + // Start from a perturbation. + std::vector x0 = {0.1, -0.2, 0.15}; + auto res = newton_cp_euclidean(mesh, x0, m); + + EXPECT_TRUE(res.converged); + EXPECT_LT(res.iterations, 30); + EXPECT_LT(res.grad_inf_norm, 1e-8); + // Strictly-convex unique minimum → converges back to ρ=0. + for (double r : res.x) EXPECT_NEAR(r, 0.0, 1e-6); +} + +// ════════════════════════════════════════════════════════════════════════════ +// 3. CP-Euclidean Newton — open mesh (boundary edges) +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonPhase9a, CPEuclidean_OpenTetrahedron_NaturalPhi_Converges) +{ + auto mesh = make_open_3face_mesh(); + auto m = setup_cp_euclidean_maps(mesh); + const int n = assign_cp_euclidean_face_dof_indices(mesh, m); + ASSERT_EQ(n, 2); + + std::vector x0(static_cast(n), 0.0); + auto G0 = cp_euclidean_gradient(mesh, x0, m); + for (auto f : mesh.faces()) { + int i = m.f_idx[f]; + if (i < 0) continue; + m.phi_f[f] -= G0[static_cast(i)]; + } + + auto res = newton_cp_euclidean(mesh, x0, m); + EXPECT_TRUE(res.converged); + EXPECT_LT(res.iterations, 30); + EXPECT_LT(res.grad_inf_norm, 1e-8); +} + +// ════════════════════════════════════════════════════════════════════════════ +// 4. Inversive-Distance Newton — natural-theta on triangle +// +// At u = 0, Bowers-Stephenson init reproduces the input edge lengths +// exactly. Natural-theta then shifts Θ so the gradient is zero, making +// u = 0 the equilibrium. Newton must converge in zero iterations. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonPhase9a, InversiveDistance_NaturalTheta_Triangle_ConvergesInZero) +{ + auto mesh = make_triangle(); + auto m = setup_inversive_distance_maps(mesh); + compute_inversive_distance_init_from_mesh(mesh, m); + + int n = 0; + for (auto v : mesh.vertices()) m.v_idx[v] = n++; + + std::vector x0(static_cast(n), 0.0); + auto G0 = inversive_distance_gradient(mesh, x0, m); + for (auto v : mesh.vertices()) { + int i = m.v_idx[v]; + m.theta_v[v] -= G0[static_cast(i)]; + } + + auto res = newton_inversive_distance(mesh, x0, m); + EXPECT_TRUE(res.converged); + EXPECT_EQ(res.iterations, 0); + EXPECT_LT(res.grad_inf_norm, 1e-10); + for (double u : res.x) EXPECT_NEAR(u, 0.0, 1e-12); +} + +// ════════════════════════════════════════════════════════════════════════════ +// 5. Inversive-Distance Newton — perturbed start on quad strip +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonPhase9a, InversiveDistance_PerturbedQuadStrip_Converges) +{ + auto mesh = make_quad_strip(); + auto m = setup_inversive_distance_maps(mesh); + compute_inversive_distance_init_from_mesh(mesh, m); + + // Pin vertex 0; index the rest. + auto vit = mesh.vertices().begin(); + m.v_idx[*vit++] = -1; + int n = 0; + for (; vit != mesh.vertices().end(); ++vit) m.v_idx[*vit] = n++; + + // Natural-theta with the pin in place. + std::vector x0(static_cast(n), 0.0); + auto G0 = inversive_distance_gradient(mesh, x0, m); + for (auto v : mesh.vertices()) { + int i = m.v_idx[v]; + if (i >= 0) m.theta_v[v] -= G0[static_cast(i)]; + } + + // Perturb away from the equilibrium and watch it return. + std::vector x_pert(static_cast(n), -0.05); + auto res = newton_inversive_distance(mesh, x_pert, m); + + EXPECT_TRUE(res.converged); + EXPECT_LT(res.iterations, 30); + EXPECT_LT(res.grad_inf_norm, 1e-8); + // Strictly-convex unique minimum on the open domain → back to 0. + for (double u : res.x) EXPECT_NEAR(u, 0.0, 1e-6); +} + +// ════════════════════════════════════════════════════════════════════════════ +// 6. Inversive-Distance Newton — tetrahedron (closed mesh) +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonPhase9a, InversiveDistance_PerturbedTetrahedron_Converges) +{ + auto mesh = make_tetrahedron(); + auto m = setup_inversive_distance_maps(mesh); + compute_inversive_distance_init_from_mesh(mesh, m); + + // Closed mesh — pin one vertex to remove the gauge mode. + auto vit = mesh.vertices().begin(); + m.v_idx[*vit++] = -1; + int n = 0; + for (; vit != mesh.vertices().end(); ++vit) m.v_idx[*vit] = n++; + + std::vector x0(static_cast(n), 0.0); + auto G0 = inversive_distance_gradient(mesh, x0, m); + for (auto v : mesh.vertices()) { + int i = m.v_idx[v]; + if (i >= 0) m.theta_v[v] -= G0[static_cast(i)]; + } + + std::vector x_pert(static_cast(n), -0.1); + auto res = newton_inversive_distance(mesh, x_pert, m); + + EXPECT_TRUE(res.converged); + EXPECT_LT(res.iterations, 30); + EXPECT_LT(res.grad_inf_norm, 1e-8); +} + +// ════════════════════════════════════════════════════════════════════════════ +// 7. CP-Euclidean Newton — uses analytic Hessian (NOT FD) +// +// Regression guard: verify the solver actually calls cp_euclidean_hessian +// (the analytic 2×2-per-edge formula) rather than degenerating to a +// per-iteration FD pass. If iteration count exceeds a tight upper bound +// for a tiny mesh, that would suggest a slow inner Hessian computation +// or a wrong-sign mistake. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonPhase9a, CPEuclidean_UsesAnalyticHessian) +{ + auto mesh = make_tetrahedron(); + auto m = setup_cp_euclidean_maps(mesh); + const int n = assign_cp_euclidean_face_dof_indices(mesh, m); + + std::vector x0(static_cast(n), 0.0); + auto G0 = cp_euclidean_gradient(mesh, x0, m); + for (auto f : mesh.faces()) { + int i = m.f_idx[f]; + if (i < 0) continue; + m.phi_f[f] -= G0[static_cast(i)]; + } + + // Strong perturbation — quadratic Newton with analytic Hessian + // should still converge in a handful of iterations. + std::vector x_pert = {0.5, -0.4, 0.3}; + auto res = newton_cp_euclidean(mesh, x_pert, m); + + EXPECT_TRUE(res.converged); + EXPECT_LE(res.iterations, 10) + << "analytic Hessian: expect very fast convergence on a 3-DOF problem"; +} diff --git a/doc/roadmap/phases.md b/doc/roadmap/phases.md index 9305a4b..a15b5be 100644 --- a/doc/roadmap/phases.md +++ b/doc/roadmap/phases.md @@ -245,3 +245,54 @@ Phase 10 Global uniformization for genus g ≥ 2 None of these are required for the genus-g uniformization pipeline; they extend the breadth of methods. ``` + +--- + +## ◼ Phase 11+ — Specialised applications (optional, deferred) + +> **Status:** out-of-scope for v1.0 but recorded here so that future +> contributors don't re-discover them. Both items live in the Java +> repo as plugin sub-packages and would benefit from porting *only* +> after Phase 10 is complete (they require the period-matrix and +> fundamental-domain infrastructure to be in place first). + +``` +11a Schottky uniformisation + Java plugin: plugin/schottky/* (~12 Java files, ~3 000 LoC) + Mathematical basis: Schottky group — discrete subgroup + Γ ⊂ PSL(2,ℂ) generated by hyperbolic loxodromic + elements, fundamental domain a sphere with + 2g disjoint discs removed. + Use case: "handlebody" uniformisation, complement of + Phase 10c's Fuchsian-group representation + (Schottky represents Riemann surfaces as + quotients of domains in S² rather than of H²). + Requires: Phase 10b (period matrix) + working + Möbius-group machinery from Phase 7. + Effort: very large (4–6 weeks) — significant Java + code, complex-analytic algorithms, + substantial test design. + +11b Riemann maps (planar conformal mapping) + Java plugin: plugin/riemannmap/* (~6 Java files, ~1 500 LoC) + Mathematical basis: Riemann mapping theorem — every simply + connected proper subdomain of ℂ is conformally + equivalent to the unit disc. Discrete version + via circle packing or Schwarz-Christoffel-like + formulae. + Use case: Texture mapping of bounded planar regions; + classical conformal mapping for engineering + applications (electrostatics, fluid flow). + Requires: Phase 10b' QuasiisothermicUtility or the + CP-Euclidean machinery from Phase 9a.1 + (depending on the chosen discrete-Riemann + algorithm). + Effort: large (3–4 weeks) — smaller than Schottky + but still substantial. Heavy on + visualisation; consider porting only the + algorithmic core. + +Both items are tracked here so the project memory is preserved; they +are NOT roadmap commitments. See `research-track.md` for the formal +research-versus-port classification before starting either. +```