feat(euclidean): full analytic edge-DOF (cyclic) Hessian; Tier-2 Wente finding
Upgrades the cyclic Euclidean Hessian from block-FD to closed-form, satisfying
the novelty-statement §3.2 "analytic Hessians, not finite difference" claim for
the Euclidean path.
- euclidean_hessian.hpp: `euclidean_hessian_analytic` — closed-form cyclic
Hessian from the law-of-cosines angle derivatives
∂α_i/∂s_i = ℓ_i²/4A, ∂α_i/∂s_j = ½cot α_i − ℓ_j²/4A (Σ_j = 0),
chained to (u, λ_e) and sign-mapped to the gradient outputs (−α vertex,
+α_opp edge). Reuses euclidean_cot_weights. Block-FD kept as cross-check.
- newton_solver.hpp: newton_euclidean cyclic path now uses the analytic Hessian.
- tests: CyclicHessian_Analytic_MatchesBlockFD_Tetrahedron — analytic == block-FD
(1e-6), == gradient FD (1e-5), symmetric (1e-9). Existing cyclic convergence
oracle still GREEN with the analytic Hessian routed in.
Tier-2 (Wente) finding: wente_torus02.obj is a QUAD mesh (1240 quads) and the
Java golden comes from cyclic (quad-net) uniformization; the C++ period-matrix
pipeline is triangle-based, so a faithful bit-vs-Java τ comparison needs a
quad/cyclic pipeline (Phase 9f). Deferred and documented; golden τ = ½+i√3/2.
244/244 cgal tests pass.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Tarik Moussa
committed by
Tarik Moussa
Tarik Moussa
parent
ea5f01d73d
commit
822a27da69
@@ -298,6 +298,122 @@ inline Eigen::SparseMatrix<double> euclidean_hessian_block_fd_sym(
|
||||
return (H + Ht) * 0.5;
|
||||
}
|
||||
|
||||
// ── Analytic cyclic Hessian (vertex + edge DOFs) ─────────────────────────────
|
||||
//
|
||||
// Closed-form Jacobian of the cyclic gradient. With s_i = effective log-length
|
||||
// (Λ̃) of the edge OPPOSITE vertex i (s = 2·ln ℓ), the Euclidean corner-angle
|
||||
// derivatives are (derived from the law of cosines; Σ_j ∂α_i/∂s_j = 0):
|
||||
//
|
||||
// ∂α_i/∂s_i = ℓ_i² / (4A)
|
||||
// ∂α_i/∂s_j = ½ cot α_i − ℓ_j² / (4A) (j ≠ i)
|
||||
//
|
||||
// where A is the triangle area (4A = denom2/2, denom2 = 8A from the cot helper).
|
||||
// Chaining through s₁ = Λ̃(e₂₃), s₂ = Λ̃(e₃₁), s₃ = Λ̃(e₁₂) and
|
||||
// Λ̃_e = λ⁰_e + u_i + u_j + λ_e gives the per-face 6×6 block over the local
|
||||
// DOFs (u₁,u₂,u₃,λ₁₂,λ₂₃,λ₃₁); per-face outputs carry the gradient signs
|
||||
// (−α vertex, +α_opp edge), so the assembled matrix equals ∂G/∂x exactly.
|
||||
// This is the analytic counterpart of `euclidean_hessian_block_fd` (no FD).
|
||||
|
||||
/// Analytic Euclidean cyclic Hessian (vertex + edge DOFs), sparse.
|
||||
/// Closed-form (no finite differences); matches `euclidean_hessian_block_fd`
|
||||
/// to round-off and the analytic vertex-only Laplacian on vertex-only layouts.
|
||||
inline Eigen::SparseMatrix<double> euclidean_hessian_analytic(
|
||||
ConformalMesh& mesh,
|
||||
const std::vector<double>& x,
|
||||
const EuclideanMaps& m)
|
||||
{
|
||||
const int n = euclidean_dimension(mesh, m);
|
||||
std::vector<Eigen::Triplet<double>> trips;
|
||||
trips.reserve(static_cast<std::size_t>(36) * mesh.number_of_faces());
|
||||
|
||||
for (auto f : mesh.faces()) {
|
||||
Halfedge_index h0 = mesh.halfedge(f);
|
||||
Halfedge_index h1 = mesh.next(h0);
|
||||
Halfedge_index h2 = mesh.next(h1);
|
||||
|
||||
Vertex_index v1 = mesh.source(h0);
|
||||
Vertex_index v2 = mesh.source(h1);
|
||||
Vertex_index v3 = mesh.source(h2);
|
||||
Edge_index e12 = mesh.edge(h0);
|
||||
Edge_index e23 = mesh.edge(h1);
|
||||
Edge_index e31 = mesh.edge(h2);
|
||||
|
||||
const double u1 = eucl_dof_val(m.v_idx[v1], x);
|
||||
const double u2 = eucl_dof_val(m.v_idx[v2], x);
|
||||
const double u3 = eucl_dof_val(m.v_idx[v3], x);
|
||||
|
||||
const double lam12 = m.lambda0[e12] + u1 + u2 + eucl_dof_val(m.e_idx[e12], x); // s3
|
||||
const double lam23 = m.lambda0[e23] + u2 + u3 + eucl_dof_val(m.e_idx[e23], x); // s1
|
||||
const double lam31 = m.lambda0[e31] + u3 + u1 + eucl_dof_val(m.e_idx[e31], x); // s2
|
||||
|
||||
// Centered side lengths (scale-invariant for ℓ²/4A and cot).
|
||||
const double mu = (lam12 + lam23 + lam31) / 6.0;
|
||||
const double l12 = std::exp((lam12 - 2.0 * mu) * 0.5);
|
||||
const double l23 = std::exp((lam23 - 2.0 * mu) * 0.5);
|
||||
const double l31 = std::exp((lam31 - 2.0 * mu) * 0.5);
|
||||
|
||||
const double t12 = -l12 + l23 + l31;
|
||||
const double t23 = l12 - l23 + l31;
|
||||
const double t31 = l12 + l23 - l31;
|
||||
if (t12 <= 0.0 || t23 <= 0.0 || t31 <= 0.0) continue;
|
||||
const double l123 = l12 + l23 + l31;
|
||||
const double denom2_sq = t12 * t23 * t31 * l123;
|
||||
if (denom2_sq <= 0.0) continue;
|
||||
const double denom2 = 2.0 * std::sqrt(denom2_sq); // = 8A
|
||||
const double fourA = denom2 * 0.5; // = 4A
|
||||
|
||||
auto cw = euclidean_cot_weights(l12, l23, l31);
|
||||
if (!cw.valid) continue;
|
||||
const double cot[3] = { cw.cot1, cw.cot2, cw.cot3 }; // at v1, v2, v3
|
||||
|
||||
// ℓ_i² = squared length of edge opposite vertex i:
|
||||
// opp(v1)=e23=l23, opp(v2)=e31=l31, opp(v3)=e12=l12
|
||||
const double Lsq[3] = { l23 * l23, l31 * l31, l12 * l12 };
|
||||
|
||||
// dA[i][j] = ∂α_i/∂s_j (i,j ∈ {0,1,2}; s_j opposite v_{j+1})
|
||||
double dA[3][3];
|
||||
for (int i = 0; i < 3; ++i)
|
||||
for (int j = 0; j < 3; ++j)
|
||||
dA[i][j] = (i == j) ? (Lsq[i] / fourA)
|
||||
: (0.5 * cot[i] - Lsq[j] / fourA);
|
||||
|
||||
// ∂α_i/∂(local dof), dof order (u1,u2,u3, λ12,λ23,λ31).
|
||||
// u1 ↔ s2,s3 ; u2 ↔ s1,s3 ; u3 ↔ s1,s2 ; λ12 ↔ s3 ; λ23 ↔ s1 ; λ31 ↔ s2.
|
||||
double g[3][6];
|
||||
for (int i = 0; i < 3; ++i) {
|
||||
g[i][0] = dA[i][1] + dA[i][2]; // u1
|
||||
g[i][1] = dA[i][0] + dA[i][2]; // u2
|
||||
g[i][2] = dA[i][0] + dA[i][1]; // u3
|
||||
g[i][3] = dA[i][2]; // λ12 (s3)
|
||||
g[i][4] = dA[i][0]; // λ23 (s1)
|
||||
g[i][5] = dA[i][1]; // λ31 (s2)
|
||||
}
|
||||
|
||||
// Output slot → (angle index, sign): (−α1,−α2,−α3,+α3,+α1,+α2).
|
||||
const int out_a[6] = { 0, 1, 2, 2, 0, 1 };
|
||||
const double out_sg[6] = { -1.0, -1.0, -1.0, +1.0, +1.0, +1.0 };
|
||||
|
||||
const int idx[6] = {
|
||||
m.v_idx[v1], m.v_idx[v2], m.v_idx[v3],
|
||||
m.e_idx[e12], m.e_idx[e23], m.e_idx[e31]
|
||||
};
|
||||
|
||||
for (int oi = 0; oi < 6; ++oi) {
|
||||
if (idx[oi] < 0) continue;
|
||||
for (int dj = 0; dj < 6; ++dj) {
|
||||
if (idx[dj] < 0) continue;
|
||||
const double val = out_sg[oi] * g[out_a[oi]][dj];
|
||||
if (std::abs(val) > 1e-15)
|
||||
trips.emplace_back(idx[oi], idx[dj], val);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Eigen::SparseMatrix<double> H(n, n);
|
||||
H.setFromTriplets(trips.begin(), trips.end());
|
||||
return H;
|
||||
}
|
||||
|
||||
// ── Finite-difference Hessian check ──────────────────────────────────────────
|
||||
/// FD Hessian check for the Euclidean functional. Compares analytic
|
||||
/// `H` column-by-column to `(G(x+εeⱼ) − G(x−εeⱼ)) / (2ε)`; returns
|
||||
|
||||
Reference in New Issue
Block a user