Merge pull request 'Phase 9a-Newton + Phase 8b-Lite: complete the CGAL API surface for all 5 DCE models' (#11) from feature/phase-9a-newton into main

Reviewed-on: #11

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 <CGAL/Conformal_map/internal/parameters.h>
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+
+#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 <typename TriangleMesh,
+          typename CGAL_NP_TEMPLATE_PARAMETERS>
+Layout2D euclidean_layout(
+    TriangleMesh&                              mesh,
+    const std::vector<double>&                 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 <typename TriangleMesh,
+          typename CGAL_NP_TEMPLATE_PARAMETERS>
+Layout3D spherical_layout(
+    TriangleMesh&                              mesh,
+    const std::vector<double>&                 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 <typename TriangleMesh,
+          typename CGAL_NP_TEMPLATE_PARAMETERS>
+Layout2D hyper_ideal_layout(
+    TriangleMesh&                              mesh,
+    const std::vector<double>&                 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

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 <CGAL/Conformal_map/internal/parameters.h>
+#include <CGAL/Kernel_traits.h>
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Simple_cartesian.h>
+#include <boost/graph/graph_traits.hpp>
+
+#include "../cp_euclidean_functional.hpp"
+#include "../newton_solver.hpp"
+
+namespace CGAL {
+
+// ── Default traits for CP-Euclidean ───────────────────────────────────────────
+
+template <typename TriangleMesh,
+          typename Kernel_ = CGAL::Simple_cartesian<double>>
+struct Default_cp_euclidean_traits;
+
+template <typename K>
+struct Default_cp_euclidean_traits<CGAL::Surface_mesh<typename K::Point_3>, K>
+{
+    using Kernel        = K;
+    using FT            = typename K::FT;
+    using Point_3       = typename K::Point_3;
+    using Triangle_mesh = CGAL::Surface_mesh<Point_3>;
+
+    using Vertex_descriptor   = typename boost::graph_traits<Triangle_mesh>::vertex_descriptor;
+    using Halfedge_descriptor = typename boost::graph_traits<Triangle_mesh>::halfedge_descriptor;
+    using Edge_descriptor     = typename boost::graph_traits<Triangle_mesh>::edge_descriptor;
+    using Face_descriptor     = typename boost::graph_traits<Triangle_mesh>::face_descriptor;
+
+    // CP-Euclidean property maps — note the *face* DOF index map.
+    using Face_index_pmap = typename Triangle_mesh::template Property_map<Face_descriptor,   int>;
+    using Theta_e_pmap    = typename Triangle_mesh::template Property_map<Edge_descriptor,   FT>;
+    using Phi_f_pmap      = typename Triangle_mesh::template Property_map<Face_descriptor,   FT>;
+};
+
+// ── 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 <typename FT = double>
+struct Circle_packing_result
+{
+    /// Face DOFs `ρ_f = log R_f` (length = num_faces(mesh); pinned face = 0).
+    std::vector<FT> 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<P>`.
+\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<FT>` 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 <typename TriangleMesh,
+          typename CGAL_NP_TEMPLATE_PARAMETERS>
+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<Point_type>::Kernel;
+    using Default_traits = Default_cp_euclidean_traits<TriangleMesh, Default_kernel>;
+    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<FT> 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<double> x0(static_cast<std::size_t>(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<std::size_t>(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<std::size_t>(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

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<P>` 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<FT>` 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 <typename TriangleMesh,
+          typename CGAL_NP_TEMPLATE_PARAMETERS>
+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<Point_type>::Kernel;
+    using Default_traits = Default_conformal_map_traits<TriangleMesh, Default_kernel>;
+    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<FT> 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<double> x0(static_cast<std::size_t>(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<std::size_t>(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<std::size_t>(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 <typename FT = double>
+struct Hyper_ideal_map_result
+{
+    /// Vertex DOFs `b_v` (length = num_vertices(mesh); pinned vertices = 0).
+    std::vector<FT> b_per_vertex;
+    /// Edge DOFs `a_e` (length = num_edges(mesh); pinned edges = 0).
+    std::vector<FT> 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 <typename TriangleMesh,
+          typename CGAL_NP_TEMPLATE_PARAMETERS>
+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<Point_type>::Kernel;
+    using Default_traits = Default_conformal_map_traits<TriangleMesh, Default_kernel>;
+    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<FT> 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<double> x0(static_cast<std::size_t>(n), 0.0);
+    for (auto v : mesh.vertices()) {
+        int i = maps.v_idx[v];
+        if (i >= 0) x0[static_cast<std::size_t>(i)] = 1.0;
+    }
+    for (auto e : mesh.edges()) {
+        int i = maps.e_idx[e];
+        if (i >= 0) x0[static_cast<std::size_t>(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<std::size_t>(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<std::size_t>(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

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 <CGAL/Conformal_map/internal/parameters.h>
+#include <CGAL/Kernel_traits.h>
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Simple_cartesian.h>
+#include <boost/graph/graph_traits.hpp>
+
+#include <CGAL/Discrete_conformal_map.h>   // for Conformal_map_result<FT>
+
+#include "../inversive_distance_functional.hpp"
+#include "../newton_solver.hpp"
+
+namespace CGAL {
+
+// ── Default traits for Inversive-Distance ────────────────────────────────────
+
+template <typename TriangleMesh,
+          typename Kernel_ = CGAL::Simple_cartesian<double>>
+struct Default_inversive_distance_traits;
+
+template <typename K>
+struct Default_inversive_distance_traits<CGAL::Surface_mesh<typename K::Point_3>, K>
+{
+    using Kernel        = K;
+    using FT            = typename K::FT;
+    using Point_3       = typename K::Point_3;
+    using Triangle_mesh = CGAL::Surface_mesh<Point_3>;
+
+    using Vertex_descriptor   = typename boost::graph_traits<Triangle_mesh>::vertex_descriptor;
+    using Edge_descriptor     = typename boost::graph_traits<Triangle_mesh>::edge_descriptor;
+
+    // Inversive-distance specific property maps.
+    using Vertex_index_pmap = typename Triangle_mesh::template Property_map<Vertex_descriptor, int>;
+    using Theta_v_pmap      = typename Triangle_mesh::template Property_map<Vertex_descriptor, FT>;
+    using R0_pmap           = typename Triangle_mesh::template Property_map<Vertex_descriptor, FT>;
+    using I_e_pmap          = typename Triangle_mesh::template Property_map<Edge_descriptor,   FT>;
+};
+
+// ── 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<P>`.
+\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<FT>` 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 <typename TriangleMesh,
+          typename CGAL_NP_TEMPLATE_PARAMETERS>
+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<Point_type>::Kernel;
+    using Default_traits = Default_inversive_distance_traits<TriangleMesh, Default_kernel>;
+    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<FT> 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<double> x0(static_cast<std::size_t>(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<std::size_t>(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<std::size_t>(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

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 <Eigen/SparseCholesky>
 #include <Eigen/SparseQR>
 #include <Eigen/OrderingMethods>
@@ -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<double>        x0,
+    const CPEuclideanMaps&     m,
+    double                     tol      = 1e-8,
+    int                        max_iter = 200)
+{
+    std::vector<double> x = x0;
+    const int n = static_cast<int>(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<const Eigen::VectorXd> 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<double>& 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<double>           x0,
+    const InversiveDistanceMaps&  m,
+    double                        tol      = 1e-8,
+    int                           max_iter = 200,
+    double                        hess_eps = 1e-5)
+{
+    std::vector<double> x = x0;
+    const int n = static_cast<int>(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<double>& xc) -> Eigen::SparseMatrix<double> {
+        std::vector<Eigen::Triplet<double>> trips;
+        trips.reserve(static_cast<std::size_t>(n) * 16);   // sparse heuristic
+
+        std::vector<double> xp = xc, xm = xc;
+        for (int j = 0; j < n; ++j) {
+            const std::size_t sj = static_cast<std::size_t>(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<std::size_t>(i)]
+                            - Gm[static_cast<std::size_t>(i)])
+                             / (2.0 * hess_eps);
+                if (std::abs(val) > 1e-15)
+                    trips.emplace_back(i, j, val);
+            }
+        }
+        Eigen::SparseMatrix<double> H(n, n);
+        H.setFromTriplets(trips.begin(), trips.end());
+        // Symmetrise — FD rounding may introduce tiny asymmetries.
+        Eigen::SparseMatrix<double> 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<const Eigen::VectorXd> 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<double>& 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

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
+    #    <CGAL/Discrete_*.h> public API + Conformal_layout.h wrapper.
+    test_cgal_phase8b_lite.cpp
 )
 
 target_include_directories(conformallab_cgal_tests SYSTEM PRIVATE

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 <CGAL/Discrete_conformal_map.h>
+#include <CGAL/Discrete_circle_packing.h>
+#include <CGAL/Discrete_inversive_distance.h>
+#include <CGAL/Conformal_layout.h>
+
+#include "mesh_builder.hpp"
+#include "conformal_mesh.hpp"
+
+#include <gtest/gtest.h>
+#include <cmath>
+
+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<double> 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()));
+    }
+}

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 <gtest/gtest.h>
+#include <vector>
+
+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<double> x0(static_cast<std::size_t>(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<std::size_t>(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<double> x0_zero(static_cast<std::size_t>(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<std::size_t>(i)];
+    }
+
+    // Start from a perturbation.
+    std::vector<double> 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<double> x0(static_cast<std::size_t>(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<std::size_t>(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<double> x0(static_cast<std::size_t>(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<std::size_t>(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<double> x0(static_cast<std::size_t>(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<std::size_t>(i)];
+    }
+
+    // Perturb away from the equilibrium and watch it return.
+    std::vector<double> x_pert(static_cast<std::size_t>(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<double> x0(static_cast<std::size_t>(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<std::size_t>(i)];
+    }
+
+    std::vector<double> x_pert(static_cast<std::size_t>(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<double> x0(static_cast<std::size_t>(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<std::size_t>(i)];
+    }
+
+    // Strong perturbation — quadratic Newton with analytic Hessian
+    // should still converge in a handful of iterations.
+    std::vector<double> 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";
+}

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.
+```