The dense-Jacobian branch of _update_gradient_JTCJ_island skips an
iteration when J[ic1, idof_i] is zero, but in the dim1 != dim2 case
the swap-pair contribution hcone * J[ic2, idof_i] * J[ic1, idof_j]
still needs to write — the swap-pair-only block under 'if dim1id !=
dim2id' only ran when J[ic1, i] AND J[ic2, j] were both nonzero. Cells
where exactly one of those four J entries was zero silently dropped
the corresponding cone contribution.

Per cell H[a, b] for an off-diagonal cone pair (dim1, dim2) the
contribution should be
  hcone * (J[ic1, a] * J[ic2, b] + J[ic2, a] * J[ic1, b])
which the monolithic _update_gradient_JTCJ_dense computes correctly
(solver.py:2700-2704). This change keeps the original loop (which
already handles the first term) and adds a second loop with the same
shape but ic1/ic2 swapped, gated on dim1id != dim2id, for the second
term. Both loops keep the aggressive J==0 early-skip from the
original.

The existing scalar _cholesky_factorize_solve_island absorbs the
resulting indefinite per-island H via its mid-factorization clamp
(s <= 1e-6 -> s = 1e-6), so no existing test fails. The bug becomes
visible when the per-island solve uses any Cholesky variant without
diagonal clamping (e.g. wp.tile_cholesky on the per-island H
sub-block), producing NaN qacc on configurations such as
constraints.xml keyframe 2 with elliptic cone + Newton + dense
Jacobian.