Phase transition problems on curved surfaces can lead to a panopticon of fascinating patterns. In this paper we consider finite element approximations of phase field models with a spatially inhomogeneous and anisotropic surface energy density. The problems are either posed in $\mathbb{R}^3$ or on a two-dimensional hypersurface in $\mathbb{R}^3.$ In the latter case, a fundamental choice regarding the anisotropic energy density has to be made. One possibility is to use a density defined in the ambient space $\mathbb{R}^3.$ However, we propose and advocate for an alternative, where a density is defined on a fixed chosen tangent space, and is then moved along geodesics to the other tangent spaces. Our numerical method can be employed in all of the above situations, where for the problems on hypersurfaces the algorithm uses parametric finite elements. We prove an unconditional stability result for our schemes and present several numerical experiments, including for the modelling of ice crystal growth on a sphere.