Return the whole boundary traction from the constrained solver (#607)

[lmoresi]

Closes #607.

Stokes_Constrained assembles the momentum row's boundary term as
(h + r(n·u − g))·n, so the traction holding the boundary is that sum.
multiplier() returns h, and topography() was built on h alone — short by
the augmented-Lagrangian share, which is r times the discrete constraint
residual.

With the viscosity-weighted default r = 1e4·μ(x) that share is a few per cent
of the surface traction on a uniform-viscosity annulus. Across SolCx's 1e6
viscosity step, where r reaches 1e10 on the stiff half, it is most of it:

route	max abs(h)	corr with exact	relative l2
multiplier() as returned	0.042	−0.53	1.04
traction() (this PR)	0.382	+0.999	0.047
exact SolCx.topography_top	0.383	—	—
projected n·σ·n, same solve	0.383	+0.997	0.075

The corrected read is the most accurate surface stress available on that problem
— better than the projection — and the solve itself was never at fault
(velocity error 8.8e-06, max abs(n·u) 1.7e-10).

What changed

• traction(boundary) — new, returns h + r(n·u − g) symbolically.
• topography() — built on it, so the public dynamic-topography path is
  correct with no change at the call site.
• multiplier() — unchanged behaviour (it returns h) and now documents what
  it is not.
• Docstrings that claimed accuracy is independent of r now say the CONSTRAINT
  and the traction() read are independent of it, and h alone is not.

Why it survived validation

The existing check scored corr(λ, −n·σ·n) = 0.9999. A correlation is scale-free
and cannot see a systematic amplitude deficit, which is precisely what a missing
share of the load is. tests/test_1063_constrained_traction.py scores a relative
l2 against the exact SolCx topography and carries the bare multiplier as its
negative control (it must read > 0.5), so a regression that quietly reverts to
h cannot pass. test_1061's API test asserted topography == h/scale, which
pinned the defect; it now asserts the outcome.

The identity behind it

At convergence M_Γ(h + r(n·u − g)) balances the volume residual restricted to
the boundary, which is the CBF nodal load (Zhong, Gurnis & Hulbert 1993). So the
multiplier route and boundary_normal_traction on a rotated constraint are the
same computation; measured across two solves they agree to 3–5%, inside each
route's own error. This is also why the free-surface work was right to reject the
multiplier as returned and keep the rotated lid + CBF.

Documentation

• CONSTRAINED_FREESLIP_MULTIPLIER.md: dated correction, corrected CBF identity.
• docs/advanced/curved-boundary-conditions.md: the penalty free-slip recipe now
  uses mesh.boundary_normal rather than mesh.Gamma, with the measurement that
  a facet-normal penalty does not converge — velocity error 0.60 flat and stress
  error 0.21 → 0.26 under refinement, while the leak reads 1e-5 throughout. The
  projected-normal recipe is marked superseded.
• CHANGELOG.md entry.

Tests

tests/test_1061, 1062, 1063, 1064, 1024: 26 passed, 2 xfailed.

Related, filed separately and NOT addressed here: #608 (the reaction reported at
a corner shared with an essential BC), #614 (boundary_flux returns ~1e12 on
Stokes_Constrained), #616 (rotated free slip does not reproduce the equivalent
component condition on a flat wall, 2e-3, with an exact linear solve on both
sides).

Underworld development team with AI support from Claude Code

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[lmoresi]

Extended to cover the other two paths, so the fix lands as one consistent set
rather than leaving them to follow.

The corner was never constrained (#616, and the mechanism behind #608)

build_rotation skipped any boundary node with a velocity DOF constrained out of
the global vector — every node where a rotated wall meets a wall held by an
essential condition. The node was dropped from the rotation, so the wall-normal
component was never constrained there. The wall leaked at its own end points
while every interior node was exact.

Unit box, rotated lid, component free slip on the other three walls, uniform
viscosity:

	before	after
max abs(u_y) on the lid	4.0e-03 (16% of abs(u)max, both corners)	0.0
vs a monolithic LU solve of the component-Dirichlet problem	2.0e-03	8.9e-09
same, with _rotated_use_lu on both sides	2.1e-03	8.9e-09

The second row with an exact linear solve on both sides is the point: it was
never a convergence artefact, it was a different discrete problem.

The fix: a partially constrained node keeps the constraint on what is left. With
the pinned components at zero, n·v = 0 reduces to n_F·v_F = 0 on the free
subspace, so the frame is built there and its normal rows constrained. A node
whose normal lies entirely in the pinned subspace is already implied and is still
skipped; a node carrying a prescribed v_n datum keeps the old behaviour and says
so once (reducing the affine constraint against the pinned values is not
implemented).

Guard: tests/test_1066_rotated_meets_essential.py. Existing suites unaffected —
42 rotated / Nitsche / constrained tests pass.

boundary_flux on a constrained boundary (#614) — and a correction

I reported that it returns ~1e12. It returns ~0 (rms 3.9e-13 against a
traction of 0.37); the 1e12 was a ratio whose denominator was that near-zero. And
it is expected: a multiplier-constrained boundary has no reaction left to read,
because the constraint enters the same row it constrains. That is the #607
identity from the other side.

It now raises and points at traction() rather than handing back a quiet zero.

Issue #614 is corrected on the issue itself.