[lachlangrose]

Stacked on #345 (phase 0). Merge that first; this PR then retargets to master.

What

• Add centroid and gauss4 tet quadrature rules (_quadrature.py).
• evaluate_shape_derivatives returns NaN rows and element -1 for points outside the mesh, not LinAlgError. P1 and P2 now return the same (dN, elements, inside). Callers updated.
• New P2FoldInterpolator(P2Interpolator) with the three fold constraint families, sampled by quadrature. Row weight is w * sqrt(q_i), applied after normalisation, with no volume factor. Not yet registered in the factory (phase 2).
• Fix bug 7 in add_constraints_to_least_squares: per-row weights for 3D input use np.repeat, not np.tile.

Notes

• The normalisation target is not scaled by fold_normalisation, unlike DiscreteFoldInterpolator. A code comment records this.
• 7 failures in tests/unit/.../test_2d_p1_p2_support.py also fail on the parent commit (element_size cross-product error).

Smoke test (fold_sparse, 2560 tets)

	centroid	gauss4	DFI (P1)
value rms	0.0033	0.0165	0.112
gradient angle median / p90 (deg)	1.36 / 5.3	0.86 / 2.1	2.46 / 10.4

Memory at about 5e5 tets: the run hit a MemoryError. Extrapolated peak is about 9 GB (centroid) and 14 GB (gauss4).

Test plan

• loop_common and loop_interpolation unit tests
• Phase 4 re-runs the validation numbers