fix(regularisers): correct diag_bending/diag_all corner cross-term

[balbasty]

Closes#48

The defect

The boundary-corrected diagonal in diag_bending (field 2d/3d, flow 2d/3d) and diag_all (flow 2d/3d) expands each corner weight as

w * (fx0*fy0 + fx1*fy0 + fx1*fy0 + fx1*fy1)

fx1*fy0 is counted twice and fx0*fy1 is dropped. The correct expansion of (fx0+fx1)*(fy0+fy1) is

w * (fx0*fy0 + fx1*fy0 + fx0*fy1 + fx1*fy1)

Why that is the right expansion

Not by assertion — by comparison with the matvec in the same file. matvec_bending accumulates the four corner neighbours explicitly, e.g. regularisers/field/2d.h:383:

+ w11 * (get(x0+y0, fx0*fy0) + get(x1+y0, fx1*fy0) +
get(x0+y1, fx0*fy1) + get(x1+y1, fx1*fy1))

The diagonal's boundary correction is the same four corners folded back onto the centre voxel, so it must carry the same four sign products. Term-for-term the diagonal matched matvec on the first, second and fourth, and disagreed on the third. field/3d.h:418-423 shows the same one-for-one correspondence for all three of w110/w101/w011.

Why it went unnoticed

The error cancels wherever both axes' one-sided boundary signs agree: with fx0==fx1 and fy0==fy1, buggy and correct forms are both 4*fx0*fy0. It only bites where the signs differ — e.g. fx0=-1, fx1=+1 gives fy0+fy1 instead of 0. That is a boundary-only, sign-flipping-bound-only effect, so the existing interior-only diagonal checks could never see it.

Impact: it corrupts the shipped field_diag / flow_diag output at boundary voxels under sign-flipping boundary conditions (Zero / DST1 / DST2) whenever bending is active, and therefore any Jacobi preconditioner built from that diagonal.

Sites

12 terms across 6 functions in 4 files:

file	function(s)	terms
regularisers/field/2d.h	diag_bending	1
regularisers/field/3d.h	diag_bending	3
regularisers/flow/2d.h	diag_bending, diag_all	1 + 1
regularisers/flow/3d.h	diag_bending, diag_all	3 + 3

After the change, a grep for the doubled-term form under regularisers/ returns nothing:

$ grep -rE "\(f[xyz]0\*f[xyz]0 \+ (f[xyz]1\*f[xyz]0) \+ \1 \+" regularisers/
(no matches)

Verification

Full CPU suite (fastfields-cpu-lib, the project's primary automated gate) built against this branch through the dev-tree symlinks, make -j4 test CXX=clang++ — all 13 suites pass, 0 failures:

test_distance 2352 test_distance_mesh 4622
test_distance_spline 704 test_posdef 4012
test_pushpull 308 test_pushpull_backward 6381
test_reg_field 18284 test_reg_flow 11415
test_reg_op 186 test_resize 630
test_restrict 65 test_splinc 4577
test_solve_field (pass)

No existing test changes value, which is the expected result: nothing currently on main covers the buggy boundary corners. A regression test that does fail without this change is going up separately against fastfields-cpu-lib — this repo has no test suite of its own — and only passes once this lands and the pin chain reaches it.

Only column-aligned expressions changed; git diff --check is clean and no whitespace or alignment was disturbed.

Provenance

The same defect was independently identified and fixed on the teeny branch in 5ed785a. fastfields-kernels#49 was an attempt to port that fix to main, but it was filed from the claude-fastfields-to-teeny workstream, whose PRs must always target the teeny branch — so it was mis-filed rather than wrong. Its analysis is sound and was useful prior art; credit to it. This PR lands the equivalent fix on main properly, re-derived and re-verified from scratch against current main (main has since moved these functions from static methods using Impl:: qualification onto instance methods holding a bound::dyn<b>, so the surrounding code differs).

🤖 Generated with Claude Code

https://claude.ai/code/session_016AjQcY78NgbagPSbPJRr6Z

---

Generated by Claude Code