  We introduce a novel type of approximation spaces for functions with values
in a nonlinear manifold. The discrete functions are constructed by piecewise
polynomial interpolation in a Euclidean embedding space, and then projecting
pointwise onto the manifold. We show optimal interpolation error bounds with
respect to Lebesgue and Sobolev norms. Additionally, we show similar bounds for
the test functions, i.e., variations of discrete functions. Combining these
results with a nonlinear C\'ea lemma, we prove optimal $L^2$ and $H^1$
discretization error bounds for harmonic maps from a planar domain into a
smooth manifold. All these error bounds are also verified numerically.
