FYI: Finally, the nodal (Lagrange) basis has been fixed [1]. It should be possible now to use nodal fields with approximation orders greater than 2 also in 3D.

Cheers, r.

[1] https://github.com/sfepy/sfepy/issues/205