Modeling large datasets through Multivariate Functional Approximations (MFA) plays a critical role in scientific analysis and visualization workflows. However, this requires scalable data partitioning approaches to compute MFA representations in a reasonable amount of time. We propose a fully parallel and efficient method for computing MFA with B-spline bases without sacrificing the reconstructed solution accuracy or continuity. Our approach reduces the total work per task and uses a restricted Additive Schwarz (RAS) method to converge control point data across subdomain boundaries. We also provide a detailed analysis of the parallel approach with domain decomposition solvers to minimize subdomain error residuals and recover high-order continuity with optimal communication cost determined by the overlap regions in the RAS implementation. In contrast to previous methods that generally only recover $C^1$ continuity for arbitrary B-spline order $p$ or required post-processing to blend discontinuities in the reconstructed data, the accuracy of the MFA remains bounded as the number of subdomains is increased. We demonstrate the effectiveness of our approach using analytical and scientific datasets in 1, 2, and 3 dimensions and show that it is highly scalable (due to bounded outer iteration counts) and that the parallel performance at scale is directly proportional to the nearest-neighbor communication implementations.