We introduce geodesic finite elements as a new way to discretize the non-linear configuration space of a geometrically exact Cosserat rod. These geodesic finite elements naturally generalize standard one-dimensional finite elements to spaces of functions with values in a Riemannian manifold. For the special orthogonal group, our approach reproduces the interpolation formulas of Crisfield and Jelenić. Geodesic finite elements are conforming and lead to objective and path-independent problem formulations. We introduce geodesic finite elements for general Riemannian manifolds, discuss the relationship between geodesic finite elements and coefficient vectors, and estimate the interpolation error. Then we use them to find static equilibria of hyperelastic Cosserat rods. Using the Riemannian trust-region algorithm of Absil et al. we show numerically that the discretization error depends optimally on the mesh size.