I have finally reimplemented some of the core functionality to allow higher order elements (uniform order), see the 'ori' branch  (the name stands for 'orientation', because the key part that was missing had been the orientation of faces allowing to define face degrees of freedom).
I tested only the very basic things, so bugs are to be expected. Any testing is more than welcome, You can try, for example, to modify examples/diffusion/poisson.py as follows:
There are still things missing for it to be really useful:
higher order quadratures - try to set 'approx_order' : 4 in the example above and you get 'warp violation' error. This is because the approximation order is too high w.r.t. the quadrature (3 is currently our highest order for tetrahedrons).
proper visualization - the higher order degrees of freedom are just thrown away now - only vertex values are used (or cell values for zero-order approximations).
hierarchical finite element basis - the Lagrange polynomial basis that is used now has very bad properties when the approximation order is increasing, leading to a matrix with a large condition number.
Any comments/suggestions/help welcome! r.