Skip to content

Commit 9a704ea

Browse files
zferguszfergus
committed
Add edge angle hessian
1 parent fb4ed31 commit 9a704ea

3 files changed

Lines changed: 333 additions & 6 deletions

File tree

‎src/ipc/geometry/angle.cpp‎

Lines changed: 287 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -59,4 +59,291 @@ Eigen::Vector<double, 12> dihedral_angle_gradient(
5959
return dsin_dx * cos_theta - dcos_dx * sin_theta;
6060
}
6161

62+
namespace {
63+
inline Eigen::Vector3d cross(
64+
Eigen::ConstRef<Eigen::Vector3d> a, Eigen::ConstRef<Eigen::Vector3d> b)
65+
{
66+
return a.cross(b);
67+
}
68+
} // namespace
69+
70+
Matrix12d dihedral_angle_hessian(
71+
Eigen::ConstRef<Eigen::Vector3d> x0,
72+
Eigen::ConstRef<Eigen::Vector3d> x1,
73+
Eigen::ConstRef<Eigen::Vector3d> x2,
74+
Eigen::ConstRef<Eigen::Vector3d> x3)
75+
{
76+
const Eigen::Vector3d n0 = triangle_normal(x0, x1, x2);
77+
const Eigen::Vector3d n1 = triangle_normal(x1, x0, x3);
78+
79+
// -------------------------------------------------------------------------
80+
// Jacobian of n0 and n1 w.r.t. all 12 DOFs
81+
// dn0_dx: n0 depends on (x0, x1, x2), not x3 → zero last 3 cols
82+
// dn1_dx: n1 = triangle_normal(x1, x0, x3), permute cols (x1→0..2, x0→3..5,
83+
// x3→6..8) → x2 block zero
84+
// -------------------------------------------------------------------------
85+
86+
Eigen::Matrix<double, 3, 12> dn0_dx;
87+
dn0_dx.leftCols<9>() = triangle_normal_jacobian(x0, x1, x2);
88+
dn0_dx.rightCols<3>().setZero();
89+
90+
Eigen::Matrix<double, 3, 12> dn1_dx;
91+
const std::array<int, 9> idx = { { 3, 4, 5, 0, 1, 2, 9, 10, 11 } };
92+
dn1_dx(Eigen::all, idx) = triangle_normal_jacobian(x1, x0, x3);
93+
dn1_dx.middleCols<3>(6).setZero();
94+
95+
// -------------------------------------------------------------------------
96+
// Hessians of n0 and n1 — shape (27 × 9) each, then embedded in 12-DOF
97+
// space
98+
//
99+
// dn0: lives in the first 9 DOFs (x0, x1, x2).
100+
// d2n0_full[3*i, j] with i,j ∈ {0..8}
101+
//
102+
// dn1 = triangle_normal(x1, x0, x3): lives in DOFs {x1, x0, x3}
103+
// mapped back: local order (x1→3..5, x0→0..2, x3→9..11)
104+
// -------------------------------------------------------------------------
105+
106+
// Raw hessians in local (9-DOF) coordinate systems
107+
const Eigen::Matrix<double, 27, 9> d2n0_local =
108+
triangle_normal_hessian(x0, x1, x2);
109+
const Eigen::Matrix<double, 27, 9> d2n1_local =
110+
triangle_normal_hessian(x1, x0, x3);
111+
112+
// We will access them on-the-fly via the index maps rather than building
113+
// full (12-DOF) tensors. The local DOF → global DOF maps are:
114+
// n0: local col k → global col k (for k in 0..8), col 9..11 → zero
115+
// n1: local col k → global col idx[k] (idx = {3,4,5,0,1,2,9,10,11})
116+
117+
// Helper: retrieve d2nX_k / (d xi d xj) for normal component k (0..2),
118+
// global DOF indices p and q.
119+
// Returns 0 if the global index falls outside the local support.
120+
121+
// Global to local maps (-1 means "not in support → zero")
122+
// n0 support: global 0..8 → local 0..8; global 9..11 → -1
123+
auto g2l_n0 = [](int g) -> int { return (g < 9) ? g : -1; };
124+
125+
// n1 support: global {0,1,2}→local{3,4,5}, global{3,4,5}→local{0,1,2},
126+
// global{9,10,11}→local{6,7,8}, global{6,7,8}→-1
127+
const std::array<int, 12> g2l_n1_map = { { 3, 4, 5, 0, 1, 2, -1, -1, -1, 6,
128+
7, 8 } };
129+
auto g2l_n1 = [&g2l_n1_map](int g) -> int { return g2l_n1_map[g]; };
130+
131+
// -------------------------------------------------------------------------
132+
// Hessian of the normalized edge e — non-zero only for global DOFs 0..5.
133+
// We store two 3×3 blocks:
134+
// He_block[k](p,q) = ∂²e_k / (∂x_{p} ∂x_{q}) with x_p, x_q ∈ {x0,x1}
135+
//
136+
// e = f(v) = v/||v||, v = x1 - x0
137+
// ∂e_k/∂v_p = J_e[k,p]
138+
// ∂²e_k/(∂v_p ∂v_q) = He[k][p,q] (from
139+
// normalization_and_jacobian_and_hessian)
140+
//
141+
// Chain rule for v = x1 - x0:
142+
// ∂v/∂x0 = -I, ∂v/∂x1 = +I
143+
// ∂²e_k/(∂x0_p ∂x0_q) = +He[k][p,q]
144+
// ∂²e_k/(∂x1_p ∂x1_q) = +He[k][p,q]
145+
// ∂²e_k/(∂x0_p ∂x1_q) = -He[k][p,q]
146+
// ∂²e_k/(∂x1_p ∂x0_q) = -He[k][p,q]
147+
// -------------------------------------------------------------------------
148+
149+
const auto [e, Je, He] = normalization_and_jacobian_and_hessian(x1 - x0);
150+
151+
// He is std::array<MatrixMax3d, 3> where He[k] is the 3×3 Hessian for
152+
// component k of the normalized vector.
153+
154+
// -------------------------------------------------------------------------
155+
// Jacobian of the normalized edge direction e = (x1 - x0) / ||x1 - x0||
156+
// Only x0 (cols 0..2) and x1 (cols 3..5) contribute.
157+
// ∂e/∂x0 = -J_norm, ∂e/∂x1 = +J_norm
158+
// -------------------------------------------------------------------------
159+
160+
Eigen::Matrix<double, 3, 12> de_dx = Eigen::Matrix<double, 3, 12>::Zero();
161+
de_dx.leftCols<3>() = -Je; // ∂e/∂x0
162+
de_dx.middleCols<3>(3) = Je; // ∂e/∂x1
163+
164+
// -------------------------------------------------------------------------
165+
// Scalars s = (n0 × n1) · e and c = n0 · n1
166+
// -------------------------------------------------------------------------
167+
168+
const Eigen::Vector3d m = n0.cross(n1); // n0 × n1
169+
const double sin_theta = m.dot(e);
170+
const double cos_theta = n0.dot(n1);
171+
172+
// -------------------------------------------------------------------------
173+
// Jacobian of s and c (12-vectors) — re-derived consistently with gradient
174+
// -------------------------------------------------------------------------
175+
176+
// dm_dx = ∂(n0×n1)/∂x (3×12)
177+
const Eigen::Matrix<double, 3, 12> dm_dx =
178+
cross_product_matrix(n0) * dn1_dx - cross_product_matrix(n1) * dn0_dx;
179+
180+
const Eigen::Vector<double, 12> dcos_dx =
181+
dn0_dx.transpose() * n1 + dn1_dx.transpose() * n0;
182+
183+
const Eigen::Vector<double, 12> dsin_dx =
184+
dm_dx.transpose() * e + de_dx.transpose() * m;
185+
186+
const Eigen::Vector<double, 12> dtheta_dx =
187+
dsin_dx * cos_theta - dcos_dx * sin_theta;
188+
189+
// -------------------------------------------------------------------------
190+
// Hessian of c = n0 · n1
191+
//
192+
// ∂²c/(∂xp ∂xq) = (∂²n0/(∂xp ∂xq)) · n1
193+
// + (∂n0/∂xp) · (∂n1/∂xq)
194+
// + (∂n1/∂xp) · (∂n0/∂xq)
195+
// + n0 · (∂²n1/(∂xp ∂xq))
196+
// -------------------------------------------------------------------------
197+
198+
Matrix12d H_cos;
199+
200+
// Cross-Jacobian terms (symmetric)
201+
H_cos = dn0_dx.transpose() * dn1_dx;
202+
H_cos += dn1_dx.transpose() * dn0_dx;
203+
204+
// Second-derivative terms: sum over normal components k
205+
for (int k = 0; k < 3; k++) {
206+
// Contribution from n0's hessian contracted with n1[k]
207+
// and from n1's hessian contracted with n0[k]
208+
for (int p = 0; p < 12; p++) {
209+
const int lp0 = g2l_n0(p);
210+
const int lp1 = g2l_n1(p);
211+
for (int q = 0; q < 12; q++) {
212+
const int lq0 = g2l_n0(q);
213+
const int lq1 = g2l_n1(q);
214+
// d2n0_k / (dxp dxq)
215+
if (lp0 >= 0 && lq0 >= 0) {
216+
H_cos(p, q) += n1(k) * d2n0_local(3 * lp0 + k, lq0);
217+
}
218+
// d2n1_k / (dxp dxq)
219+
if (lp1 >= 0 && lq1 >= 0) {
220+
H_cos(p, q) += n0(k) * d2n1_local(3 * lp1 + k, lq1);
221+
}
222+
}
223+
}
224+
}
225+
226+
// -------------------------------------------------------------------------
227+
// Hessian of s = m · e = (n0 × n1) · e
228+
//
229+
// ∂²s/(∂xp ∂xq) = [∂²m/(∂xp ∂xq)] · e
230+
// + [∂m/∂xp] · [∂e/∂xq]
231+
// + [∂e/∂xp] · [∂m/∂xq]
232+
// + m · [∂²e/(∂xp ∂xq)]
233+
// -------------------------------------------------------------------------
234+
235+
Matrix12d H_sin;
236+
237+
// Cross-Jacobian terms (m–e interaction)
238+
H_sin = dm_dx.transpose() * de_dx;
239+
H_sin += de_dx.transpose() * dm_dx;
240+
241+
// Second derivative of e contracted with m — only DOFs 0..5 contribute
242+
// ∂²e_k/(∂v_p ∂v_q) with signs for x0/x1 blocks
243+
// Block (x0,x0): +He[k], Block (x1,x1): +He[k],
244+
// Block (x0,x1): -He[k], Block (x1,x0): -He[k]
245+
for (int k = 0; k < 3; k++) {
246+
const Eigen::Matrix3d& Hek = He[k]; // 3×3
247+
const double mk = m(k);
248+
// (x0, x0)
249+
H_sin.block<3, 3>(0, 0) += mk * Hek;
250+
// (x1, x1)
251+
H_sin.block<3, 3>(3, 3) += mk * Hek;
252+
// (x0, x1)
253+
H_sin.block<3, 3>(0, 3) -= mk * Hek;
254+
// (x1, x0)
255+
H_sin.block<3, 3>(3, 0) -= mk * Hek;
256+
}
257+
258+
// Second derivative of m = n0 × n1 contracted with e:
259+
//
260+
// ∂²(n0 × n1)/(∂xp ∂xq)
261+
// = -[∂n1/∂xq ×] · ∂n0/∂xp - [n1×] · ∂²n0/(∂xp ∂xq)
262+
// + [∂n0/∂xq ×] · ∂n1/∂xp + [n0×] · ∂²n1/(∂xp ∂xq)
263+
//
264+
// Contracted with e:
265+
// e · ∂²m/(∂xp ∂xq)
266+
// = -eᵀ [∂n1/∂xq ×] ∂n0/∂xp - eᵀ [n1×] ∂²n0/(∂xp ∂xq)
267+
// + eᵀ [∂n0/∂xq ×] ∂n1/∂xp + eᵀ [n0×] ∂²n1/(∂xp ∂xq)
268+
//
269+
// Note: eᵀ [v×] w = e · (v × w) = (e × v) · w → (-v × e) · w
270+
// so eᵀ [v×] = (-v × e)ᵀ = (e × v)ᵀ (as a row vector)
271+
272+
// Pre-compute (e × n0) for efficiency
273+
const Eigen::Vector3d e_cross_n0 = cross(e, n0);
274+
275+
// --- Terms involving first × first derivatives (cross of Jacobian cols)
276+
// --- -eᵀ [dn1/dxq ×] dn0/dxp = (e × dn1/dxq) · dn0/dxp
277+
// = (e_cross_dn1_q) · dn0_dxp
278+
// +eᵀ [dn0/dxq ×] dn1/dxp = -(e × dn0/dxq) · dn1/dxp
279+
// Wait: eᵀ [v×] w = e · (v × w). Let v = dn0/dxq, w = dn1/dxp
280+
// => e · (dn0/dxq × dn1/dxp)
281+
// Also: e · (v × w) = -(e × v) · w... let's just use dot directly.
282+
for (int q = 0; q < 12; q++) {
283+
const Eigen::Vector3d dn1_q = dn1_dx.col(q);
284+
const Eigen::Vector3d dn0_q = dn0_dx.col(q);
285+
// e · (dn1_q × dn0_p) for all p → negate the cross and dot with e
286+
// Term: -eᵀ [dn1_q×] dn0_dxp = e · (dn1_q × dn0_dxp) ... but
287+
// [v×]w = v×w, so eᵀ [v×] w = e·(v×w) = (e×v)·w
288+
// -eᵀ [dn1_q ×] dn0_dxp = -(e × dn1_q) · dn0_dxp
289+
const Eigen::Vector3d neg_e_cross_dn1_q = -(cross(e, dn1_q));
290+
// +eᵀ [dn0_q ×] dn1_dxp = (e × dn0_q) · dn1_dxp
291+
const Eigen::Vector3d e_cross_dn0_q = cross(e, dn0_q);
292+
for (int p = 0; p < 12; p++) {
293+
H_sin(p, q) += neg_e_cross_dn1_q.dot(dn0_dx.col(p));
294+
H_sin(p, q) += e_cross_dn0_q.dot(dn1_dx.col(p));
295+
}
296+
}
297+
298+
// --- Terms involving second derivatives of n0 and n1 contracted with e ---
299+
// -eᵀ [n1×] ∂²n0/(∂xp ∂xq) = -(e × n1) · ∂²n0/(∂xp ∂xq)
300+
// = e_cross_n1_neg · ∂²n0/(∂xp ∂xq)
301+
// +eᵀ [n0×] ∂²n1/(∂xp ∂xq) = (e × n0) · ∂²n1/(∂xp ∂xq)
302+
// where eᵀ [n1×] w = (e×n1)·w
303+
// eᵀ [n0×] w = (e×n0)·w
304+
const Eigen::Vector3d neg_e_cross_n1 = -cross(e, n1); // -(e × n1)
305+
306+
for (int p = 0; p < 12; p++) {
307+
const int lp0 = g2l_n0(p);
308+
const int lp1 = g2l_n1(p);
309+
for (int q = 0; q < 12; q++) {
310+
const int lq0 = g2l_n0(q);
311+
const int lq1 = g2l_n1(q);
312+
313+
// -eᵀ [n1×] d²n0/(dxp dxq): d²n0_k/(dxp dxq) is in d2n0_local
314+
if (lp0 >= 0 && lq0 >= 0) {
315+
for (int k = 0; k < 3; k++) {
316+
H_sin(p, q) +=
317+
neg_e_cross_n1(k) * d2n0_local(3 * lp0 + k, lq0);
318+
}
319+
}
320+
// +eᵀ [n0×] d²n1/(dxp dxq)
321+
if (lp1 >= 0 && lq1 >= 0) {
322+
for (int k = 0; k < 3; k++) {
323+
H_sin(p, q) += e_cross_n0(k) * d2n1_local(3 * lp1 + k, lq1);
324+
}
325+
}
326+
}
327+
}
328+
329+
// -------------------------------------------------------------------------
330+
// Assemble the Hessian of θ = atan2(s, c)
331+
//
332+
// Since s² + c² = 1 (n0, n1 are unit normals and e ⊥ plane spanned by
333+
// n0 and n1, so (n0×n1)·e = ||n0×n1|| = sin(φ):
334+
//
335+
// H_θ = c·H_s - s·H_c (second-derivative terms)
336+
// + ∇c·(∇s)ᵀ - ∇s·(∇c)ᵀ (product-rule cross terms, antisymmetric)
337+
// - 2·∇θ·(s·∇s + c·∇c)ᵀ (denominator derivative)
338+
// -------------------------------------------------------------------------
339+
340+
Matrix12d H_theta = cos_theta * H_sin - sin_theta * H_cos
341+
+ dcos_dx * dsin_dx.transpose() - dsin_dx * dcos_dx.transpose()
342+
- dtheta_dx
343+
* ((2 * sin_theta) * dsin_dx + (2 * cos_theta) * dcos_dx)
344+
.transpose();
345+
346+
return H_theta;
347+
}
348+
62349
} // namespace ipc

‎src/ipc/geometry/angle.hpp‎

Lines changed: 19 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -5,9 +5,9 @@
55
namespace ipc {
66

77
/// @brief Compute the bending angle between two triangles sharing an edge.
8-
/// x0---x2
8+
/// x0---x3
99
/// | \ |
10-
/// x1---x3
10+
/// x2---x1
1111
/// @param x0 The first vertex of the edge.
1212
/// @param x1 The second vertex of the edge.
1313
/// @param x2 The opposite vertex of the first triangle.
@@ -20,9 +20,9 @@ double dihedral_angle(
2020
Eigen::ConstRef<Eigen::Vector3d> x3);
2121

2222
/// @brief Compute the Jacobian of the bending angle between two triangles sharing an edge.
23-
/// x0---x2
23+
/// x0---x3
2424
/// | \ |
25-
/// x1---x3
25+
/// x2---x1
2626
/// @param x0 The first vertex of the edge.
2727
/// @param x1 The second vertex of the edge.
2828
/// @param x2 The opposite vertex of the first triangle.
@@ -34,4 +34,19 @@ Eigen::Vector<double, 12> dihedral_angle_gradient(
3434
Eigen::ConstRef<Eigen::Vector3d> x2,
3535
Eigen::ConstRef<Eigen::Vector3d> x3);
3636

37+
/// @brief Compute the Hessian of the bending angle between two triangles sharing an edge.
38+
/// x0---x3
39+
/// | \ |
40+
/// x2---x1
41+
/// @param x0 The first vertex of the edge.
42+
/// @param x1 The second vertex of the edge.
43+
/// @param x2 The opposite vertex of the first triangle.
44+
/// @param x3 The opposite vertex of the second triangle.
45+
/// @return The 12x12 Hessian matrix of the bending angle with respect to the input vertices.
46+
Eigen::Matrix<double, 12, 12> dihedral_angle_hessian(
47+
Eigen::ConstRef<Eigen::Vector3d> x0,
48+
Eigen::ConstRef<Eigen::Vector3d> x1,
49+
Eigen::ConstRef<Eigen::Vector3d> x2,
50+
Eigen::ConstRef<Eigen::Vector3d> x3);
51+
3752
} // namespace ipc

‎tests/src/tests/geometry/test_angle.cpp‎

Lines changed: 27 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -69,8 +69,28 @@ TEST_CASE("Dihedral angle and gradient", "[angle][dihedral]")
6969
REQUIRE(grad.array().isFinite().all());
7070
CHECK(fd::compare_gradient(grad, fd_grad));
7171
if (!fd::compare_gradient(grad, fd_grad)) {
72-
std::cout << " Gradient:\n" << grad.transpose() << std::endl;
73-
std::cout << "FD Gradient:\n" << fd_grad.transpose() << std::endl;
72+
std::cout << " Gradient:\n" << grad.transpose() << "\n";
73+
std::cout << "FD Gradient:\n" << fd_grad.transpose() << "\n";
74+
}
75+
76+
// --- Check the Hessian against finite differences ---
77+
78+
Eigen::MatrixXd hess = dihedral_angle_hessian(x0, x1, x2, x3);
79+
Eigen::MatrixXd fd_hess;
80+
fd::finite_jacobian(
81+
x,
82+
[](const Eigen::VectorXd& x_fd) -> Eigen::VectorXd {
83+
return dihedral_angle_gradient(
84+
x_fd.segment<3>(0), x_fd.segment<3>(3), x_fd.segment<3>(6),
85+
x_fd.segment<3>(9));
86+
},
87+
fd_hess);
88+
89+
REQUIRE(hess.array().isFinite().all());
90+
CHECK(fd::compare_hessian(hess, fd_hess));
91+
if (!fd::compare_hessian(hess, fd_hess)) {
92+
std::cout << " Hessian:\n" << grad.transpose() << "\n";
93+
std::cout << "FD Hessian:\n" << fd_grad.transpose() << "\n";
7494
}
7595
}
7696

@@ -89,4 +109,9 @@ TEST_CASE(
89109
{
90110
return dihedral_angle_gradient(x0, x1, x2, x3);
91111
};
112+
113+
BENCHMARK("Dihedral angle Hessian")
114+
{
115+
return dihedral_angle_hessian(x0, x1, x2, x3);
116+
};
92117
}

0 commit comments

Comments
 (0)