3-D boundary-flux recovery supports only P1/P2 triangular traces, so dynamic topography cannot be cross-validated

[lmoresi]

In 3-D, _desmear — the primitive behind boundary_flux, boundary_normal_traction
and dynamic_topography — accepts only P1 and P2 triangular traces. Anything else
raises:

• src/underworld3/utilities/boundary_flux.py:327 — non-triangular facets (so any hex
  or prism mesh, whose boundary facets are quadrilaterals)
• src/underworld3/utilities/boundary_flux.py:331 — a facet carrying interior DOFs
  (so P3 and above)

The 2-D path has no such restriction: it handles arbitrary trace degree, building the
element mass from the actual node parameters along the edge chord.

Two consequences.

Dynamic topography has exactly one supported discretisation in 3-D. Velocity must
be P2 on tets. There is no second discretisation to check the first against, so a
recovery error that is a property of the P2 triangle — the vertex-integral checkerboard
the mass="consistent" docstring already warns about — cannot be distinguished from a
property of the physics by varying the element.

It blocked a cross-element control in #633. The vertex-versus-midpoint split under
the grad-div penalty was to be tested at P3/P2 and on hexes, to see whether it tracked
the node or the polynomial order. Both arms raise, so that discrimination was done on a
2-D annulus instead (~/+Simulations/topography_penalty_633/), where the trace degree
is free.

Worth having, in rough order of value:

1. Quadrilateral facets (hexes), P1/P2 — the tensor-product mass matrices are standard
   and it opens a genuinely different element family.
2. P3 and above on triangles — needs the facet-interior node coordinates, which the 2-D
   path already builds via _trace_interior_coords.

Neither is needed for correctness of what ships today; both are needed to validate it.