Accurately solving coupled systems of partial differential equations is often limited by high computational costs and significant implementation complexity, particularly due to the need to linearize coupled nonlinear equations, which represents a major bottleneck in material design. We propose a finite element–based physics-informed operator learning framework for coupled partial differential equations on arbitrary domains. The framework learns a PDE-parameterized solution operator directly from discretized weak forms, enabling discretization-independent predictions beyond the training resolution while strictly enforcing the governing physics. The training procedure does not require labeled simulation data or precomputed solutions. As the operator-learning backbone, various architectures from the literature can be seamlessly incorporated; in this work, we examine Fourier neural operators, Deep operator networks, and a newly proposed implicit finite operator learning approach [1]. The proposed approach is validated on coupled thermo-mechanical problems across multiple scales. At the microscale, nonlinear two- and three-dimensional representative volume elements with heterogeneous microstructures are examined, while at the macroscale, a close-to-reality industrial casting example is studied under varying boundary conditions. We also investigate the impact of different training schemes, including monolithic and staggered approaches, as well as training sample quality, on prediction performance. The results demonstrate the potential of physics-informed operator learning with a finite element–based loss as a unified and scalable framework for coupled multiphysics simulations. [1] R. N. Asl, Y. Yamazaki, K. Taghikhani, M. Muramatsu, M. Apel, A Physics-Informed Meta-Learning Framework for the Continuous Solution of Parametric PDEs on Arbitrary Geometries, arXiv preprint arXiv:2504.02459, 2025.