  We numerically benchmark methods for computing harmonic maps into the unit
sphere, with particular focus on harmonic maps with singularities. For the
discretization we compare two different approaches, both based on Lagrange
finite elements. While the first method enforces the unit-length constraint
only at the Lagrange nodes, the other one adds a pointwise projection to
fulfill the constraint everywhere. For the solution of the resulting algebraic
problems we compare a nonconforming gradient flow with a Riemannian
trust-region method. Both are energy-decreasing and can be shown to converge
globally to a stationary point of the Dirichlet energy. We observe that while
the nonconforming and the conforming discretizations both show similar
behavior, the second-order trust-region method needs less iterations than the
solver based on gradient flow.
