We consider the numerical investigation of surface bound orientational order using unit tangential vector fields by means of a gradient flow equation of a weak surface Frank–Oseen energy. The energy is composed of intrinsic and extrinsic contributions, as well as a penalization term to enforce the unity of the vector field. Four different numerical discretizations, namely a discrete exterior calculus approach, a method based on vector spherical harmonics, a surface finite element method, and an approach utilizing an implicit surface description, the diffuse interface method, are described and compared with each other for surfaces with Euler characteristic 2. We demonstrate the influence of geometric properties on realizations of the Poincaré–Hopf theorem and show examples where the energy is decreased by introducing additional orientational defects.