In this work, we introduce an efficient estimation procedure for smoothing space-time functional data observed over multidimensional irregular domains. We discuss the limitations of existing approaches when applied to large datasets of spatio-temporal functional data, observed over complex supports, such as those encountered in neuroimaging applications, and propose a novel iterative procedure for the solution of physics-informed nonparametric regression problems. The proposed method combines computational efficiency with high accuracy. Moreover, it provides the basis to address more complex and large-scale functional data analysis problems, at the population level, concerning for instance functional principal component analysis and functional clustering.