Completing multidimensional tensor-structured data with missing entries is a fundamental task for many real-world applications involving incomplete or corrupted datasets. For data with spatial or temporal side information, low-rank factorization models with smoothness constraints have demonstrated strong performance. Although effective at capturing global and long-range correlations, these models often struggle to capture short-scale, high-frequency variations in the data. To address this limitation, we propose the Generalized Least Squares Kernelized Tensor Factorization (GLSKF) framework for tensor completion. GLSKF integrates smoothness-constrained low-rank factorization with a locally correlated residual process; the resulting additive structure enables effective characterization of both global dependencies and local variations. Specifically, we define the covariance norm to enforce the smoothness of factor matrices in the global low-rank factorization, and use structured covariance/kernel functions to model the local processes. For model estimation, we develop an alternating least squares (ALS) procedure with closed-form solutions for each subproblem. GLSKF utilizes zero-padding and slicing operations based on projection matrices which preserve the Kronecker structure of covariances, facilitating efficient computations through the conjugate gradient (CG) method. The proposed framework is evaluated on four real-world datasets across diverse tasks. Experimental results demonstrate that GLSKF achieves superior performance and scalability, establishing it as a novel solution for multidimensional tensor completion.