We introduce geodesic finite elements as a conforming way to discretize partial differential equations for functions v : Ω → M, where Ω is an open subset of a Euclidean space, and M is a Riemannian manifold. These geodesic finite elements naturally generalize standard first-order finite elements for Euclidean spaces. They also generalize the geodesic finite elements proposed for d = 1 in a previous publication of the author. Our formulation is equivariant under isometries of M and, hence, preserves objectivity of continuous problem formulations. We concentrate on partial differential equations that can be formulated as minimization problems. Discretization leads to algebraic minimization problems on product manifolds Mn. These can be solved efficiently using a Riemannian trust-region method. We propose a monotone multigrid method to solve the constrained inner problems with linear multigrid speed. As an example, we numerically compute harmonic maps from a domain in R³ to S².