Physics‐informed neural networks (PINNs) are a new tool for solving boundary value problems by defining loss functions of neural networks based on governing equations, boundary conditions, and initial conditions. Recent investigations have shown that when designing loss functions for many engineering problems, using first‐order derivatives and combining equations from both strong and weak forms can lead to much better accuracy, especially when there are heterogeneity and variable jumps in the domain. This new approach is called the mixed formulation for PINNs, which takes ideas from the mixed finite element method. In this method, the PDE is reformulated as a system of equations where the primary unknowns are the fluxes or gradients of the solution, and the secondary unknowns are the solution itself. In this work, we propose applying the mixed formulation to solve multi‐physical problems, specifically a stationary thermo‐mechanically coupled system of equations. Additionally, we discuss both sequential and fully coupled unsupervised training and compare their accuracy and computational cost. To improve the accuracy of the network, we incorporate hard boundary constraints to ensure valid predictions. We then investigate how different optimizers and architectures affect accuracy and efficiency. Finally, we introduce a simple approach for parametric learning that is similar to transfer learning. This approach combines data and physics to address the limitations of PINNs regarding computational cost and improves the network's ability to predict the response of the system for unseen cases. The outcomes of this work will be useful for many other engineering applications where deep learning is employed on multiple coupled systems of equations for fast and reliable computations.

We propose a novel finite element-based physics-informed operator learning framework that allows for predicting spatiotemporal dynamics governed by partial differential equations (PDEs). The Galerkin discretized weak formulation is employed to incorporate physics into the loss function, termed finite operator learning (FOL), along with the implicit Euler time integration scheme for temporal discretization. A transient thermal conduction problem is considered to benchmark the performance, where FOL takes a temperature field at the current time step as input and predicts a temperature field at the next time step. Upon training, the network successfully predicts the temperature evolution over time for any initial temperature field at high accuracy compared to the solution by the finite element method (FEM) even with a heterogeneous thermal conductivity and arbitrary geometry. The advantages of FOL can be summarized as follows: First, the training is performed in an unsupervised manner, avoiding the need for large data prepared from costly simulations or experiments. Instead, random temperature patterns generated by the Gaussian random process and the Fourier series, combined with constant temperature fields, are used as training data to cover possible temperature cases. Additionally, shape functions and backward difference approximation are exploited for the domain discretization, resulting in a purely algebraic equation. This enhances training efficiency, as one avoids time-consuming automatic differentiation in optimizing weights and biases while accepting possible discretization errors. Finally, thanks to the interpolation power of FEM, any arbitrary geometry with heterogeneous microstructure can be handled with FOL, which is crucial to addressing various engineering application scenarios.