An innovative nonparametric method for density estimation over general two-dimensional Riemannian manifolds is proposed. The method follows a functional data analysis approach, combining maximum likelihood estimation with a roughness penalty that involves a differential operator appropriately defined over the manifold domain, thus controlling the smoothness of the estimate. The proposed method can accurately handle point pattern data over complicated curved domains. Moreover, it is able to capture complex multimodal signals, with strongly localized and highly skewed modes, with varying directions and intensity of anisotropy. The estimation procedure exploits a discretization in finite element bases, enabling great flexibility on the spatial...