Hi Alvar,
using
Code:
fes = HDiv(mesh, order=k, RT=True)
gives the Raviart-Thomas elements (for triangles and tetrahedra), which lie between polynomials of order k and k+1 and the divergence is then in L2 of order k.
For HDivSurface this option is currently not implemented in the main code, but it might come soon.
Best,
Michael