Fully intrinsic equations and boundary conditions involve only force, moment, velocity, and angular velocity variables, but no displacement or rotation variables. This paper presents variable-order finite elements for the geometrically exact, nonlinear, fully intrinsic equations for both nonrotating and rotating beams. The finite element technique allows for hp-adaptivity. Results show that these finite elements lead to very accurate solutions for the static equilibrium state as well as for modes and frequencies for infinitesimal motions about that state. For the same number of variables, the accuracy of the finite elements increases with the order of the finite element. The results based on the Galerkin approximation (which is a special case of the present approach) are the most accurate but require evaluation of complex integrals. Cubic elements are shown to provide a near optimal combination of accuracy and complexity.