\section{Introduction}

The laws of nature can be approximated and described by various means,
among which partial differential equations can reflect the relationship
of physical quantities with space and time, so they have become an important tool for describing physical laws \citet{loss_is_minimum_potential_energy}.
In real world applications, it is highly unlikely for partial differential equations to have an analytical solution, and one has to resort to numerical methods to get approximate
solutions.
There are several traditional numerical methods to obtain the approximate solutions, such as the established finite element method (FEM) \citet{zienkiewicz2005finite}, the finite difference method \citet{darwish2016finitevolumemethod}, the finite volume method \citet{leveque2007finitedifferentialmethod}, and mesh-free method \cite{zhang2016material}.
These methods have been developed and used for a long time, and they are reliable numerical methods, especially FEM in solid mechanics simulation.
Different trial and test functions produce various
finite element numerical methods. The essence of FEM is to find the
optimal solution in the approximate function space. However, constructing different elements requires a considerable brainpower cost \citet{wang2022cenn}.

Artificial intelligence has impacted many fields in the past ten years,
among which deep learning has attracted the most attention. Deep learning
is empowering various traditional fields \citet{alphago,star_game, ALEXNET,speech_recognition,machine_translation,alphafold},
and mechanics is no exception. The combination of deep learning and
mechanics has two main aspects. The first aspect is to use the strong
approximation of the neural network and the existing reliable
big data to build a reliable replacement model, e.g., the establishment
of the constitutive model \citet{kirchdoerfer2016data,kirchdoerfer2018data}. The prediction of the equivalent modulus of non-uniform materials
\citet{li2019predicting}, reverse design acoustic metamaterial \citet{li2020designing}.
This is a completely data-driven modeling method using the super-strong
approximation of neural networks \citet{super_approximation}. Once
the model is successfully trained, the efficiency and accuracy will
often be outstanding, but its disadvantage is often the lack of interpretability.
Using the neural network to fit the relationship between the data
often has poor scalability of the model, i.e., the neural network
often needs to be retrained for different problems, and even the network
structure has changed drastically \citet{ill_gradient}. Moreover, the physical mechanism behind it is often not well explained \citet{li2022ref}.
In addition, this method has challenges in solving high-dimensional
problems because it needs to use existing methods to obtain high-dimensional
data, so it faces the dimension explosion problem \citet{bc-pinn}.
The second aspect is to solve partial differential equations (PDEs)
through neural networks, and the most concerned method is physics-informed
neural network (PINNs) proposed in 2019 \citet{PINN_original_paper}.
The loss function of PINNs is composed of PDEs. PINNs have many applications in many fields, e.g., solid mechanics \citet{PINN_solid_mechanics, li2022equilibrium, goswami2022physics}, fluid mechanics \citet{cai2021physics, raissi2020hidden}, and biomechanics \citet{yin2021non, yin2022simulating}. Besides, the loss function can be changed to energy in solid mechanics, and the minimum
potential energy principle can be used for optimization. This
is the deep energy method (DEM) proposed in 2020 \citet{loss_is_minimum_potential_energy, PINN_hyperelasticity, goswami2020adaptive, goswami2020transfer}.
The advantage of PINNs is general, i.e., almost all partial differential
equations can be solved by PINNs. However, there are still many challenges to determine the hyperparameters of PINNs, especially in solid mechanics (The equations have high-order tensor and derivative), despite currently valuable research results for the selection of hyperparameters
\citet{NTK_PINN, NTK_to_get_hyperparameter_of_PINN,ill_gradient}.
The advantage of DEM is that there are few hyperparameters and the
calculation efficiency and accuracy are generally higher than those
of the strong form, but not all the partial differential equations
have corresponding energy forms, so the versatility is not as good
as the strong form, and it is necessary to construct possible field
functions of interest in advance \citet{the_comparision_of_strong_and_energy_form}. To sum up,
PINNs and DEM have paved a new way to solve partial
differential equations and have important research value.

In the linear theory of elasticity, the most critical variational principle is the dual extreme principle, i.e., the principle of minimum potential and complementary energy \citet{nonlinear_plate_theory_comple}. However, the current deep energy method lacks the form of a complementary energy principle. We believe that the advantages of the deep energy method based on the principle of potential energy and complementary energy are similar to those of the displacement finite element and stress finite element.
Stress finite element plays an important role in computational mechanics.
Therefore, it is very important to propose a deep energy method based on the principle of complementary energy. On one end of the spectrum, pure data-driven requires a lot of data. On the opposite end of the spectrum, if it is completely based on physical equations, although the theory of generalized approximate functions has great theoretical prospects, it is not easy to surpass the traditional finite element method in terms of efficiency and accuracy \citet{super_approximation}. There is a long way to go \citep{lu2022comprehensive}. Recently, DeepONet \citep{DeepOnet} based on universal approximation operator theory \citep{chen1995universal} and FNO \citep{li2020fourier} have received a lot of attention due to their practical application potential \citep{wang2021learning, wen2022u}. Since FNO is a particular case of DeepONet, we only consider DeepONet in this paper \citep{lu2022comprehensive}. Combining DeepONet with DEM to solve fracture mechanics, the computational efficiency and accuracy are greatly improved \citet{goswami2022physics}. Therefore, we believe that semi-supervised learning combined with operator and physical laws has a very broad prospect in the future of the data explosion.

In this work, we propose a deep complementary energy form (DCM) based on the principle of minimum complementary energy instead of minimum potential
energy. Our output function is the stress function that naturally satisfies the equilibrium equation. This is the first attempt to leverage the power of the PINNs energy form based on complementary energy.
We continue to expand the functionality of DCM and propose a DCM-Plus (DCM-P) algorithm where we add the terms that naturally satisfy the biharmonic equation in the Airy stress function.
Considering our current or future high-fidelity numerical results and experiment data, we develop DCM based on the DeepONet operator (DCM-O), including branch net, trunk net, basis net, and particular net. The motivation of DCM-O is to reuse the existing calculation results and not waste the previous calculation results for improving the computational efficiency of DCM. The advantages of the DCM are multifold. Since the principle of minimum complementary energy does not need to construct an admissible function field on the displacement boundary condition, it is theoretically guaranteed that DCM (based on complementary energy) is superior compared to DEM (based on potential energy) in the construction of admissible functions. At the same time, we verify that DCM is more suitable for dealing with problems dominated by displacement boundary conditions in our numerical result. The comparison between DCM and DEM is similar to the comparison between finite elements based on stress and displacement, so we propose an important supplement to the existing deep energy method. At the same time, DCM can be combined with existing data, combined with operator learning and physical equations (DCM-O). The more data, the faster the DCM calculation.

The outline of the paper is as follows. In \Cref{sec: prerequisite}, we introduce the prerequisite knowledge, including the feed-forward network, DeepONet algorithm, deep energy method, and the important stress functions in solid mechanics. In \Cref{sec: method}, we describe the methodology of the proposed method DCM, DCM-P and DCM-O. In \Cref{sec:Result}, we present the numerical result of the most common stress function, Prandtl and Airy stress function with the different boundary conditions under PINNs strong form, DEM, and DCM. The results show that DEM is more suitable for dealing with the problems dominated by the force boundary conditions, and DCM is suitable for dealing with displacement boundary conditions. Accuracy and efficiency are closely related to the derivative order of the interesting field and equation, which we should consider in selecting the energy principle. Finally, \Cref{sec:Discussion} and \Cref{sec:Conclusion} present some discussion of DCM, concluding remarks, some limitation of DCM, and possible future work. In the Appendices, we show some important proofs of DCM.
\section{Prerequisite knowledge} \label{sec: prerequisite}

\subsection{Introduction to feed-forward neural network}
There are mainly four types of neural networks, namely fully connected (feed-forward neural network), convolutional neural network \citep{lecunCNN}, recurrent neural network (LSTM \citep{LSTM} is the most used) and the most recent Transformer \citep{vaswani2017attention}. The combination of machine learning and computational mechanics mainly depends on the powerful fitting ability of the neural network \citep{lu2022comprehensive, hornik1989multilayer}, so the mainstream is still using the feed-forward neural network in scientific applications with machine learning at present. Therefore, this article mainly introduces the feed-forward neural network. The feed-forward neural network is a nonlinear multiple linear regression, and the activation function adds nonlinear capabilities. The mathematical formula of the feed-forward neuron network can be expressed as:

\begin{equation}
	\begin{split}\boldsymbol{z}^{(1)} & =\boldsymbol{w}^{(1)}\cdot\boldsymbol{x}+\boldsymbol{b}^{(1)}\\
		\boldsymbol{a}^{(1)} & =\boldsymbol{\sigma}(\boldsymbol{z}^{(1)})\\
		\vdots\\
		\boldsymbol{z}^{(L)} & =\boldsymbol{w}^{(L)}\cdot\boldsymbol{a}^{(L-1)}+\boldsymbol{b}^{(L)}\\
		\boldsymbol{a}^{(L)} & =\boldsymbol{\sigma}(\boldsymbol{z}^{(L)})\\
		\boldsymbol{y} & =\boldsymbol{w}^{(L+1)}\cdot\boldsymbol{a}^{(L)}+\boldsymbol{b}^{(L+1)},
	\end{split}
\end{equation}
where $\boldsymbol{x}$ and $\boldsymbol{y}$ are the input and output, respectively. $\boldsymbol{z}$ is the linear output. $\boldsymbol{\sigma}$ is the activation function applied to $\boldsymbol{z}$. $\boldsymbol{w}$ and $\boldsymbol{b}$ are the trainable parameters of the neural network. $\boldsymbol{z}^{(l)}$ is the linear transformation of the neurons $\boldsymbol{a}^{(l-1)}$ in the previous layer, and the  $1\leq l\leq L$ layers are  hidden layers. $\boldsymbol{w}^{(m)}$ is the linear transformation weight from layer $m-1$ to layer $m$, and $\boldsymbol{b}^{(m)}$ is the bias of layer $m$. $\boldsymbol{a}^{(l)}$ is the output of $\boldsymbol{z}^{(l)}$ through the activation function $\sigma$, and the activation function $\sigma$ is a nonlinear function, such as {\rm tanh} and {\rm sigmoid}. In this paper, we uniformly use the {\rm tanh} function as the activation function. The expression for the tanh function is:

\begin{equation}
	{\rm {\rm tanh}}=\frac{e^{x}-e^{-x}}{e^{x}+e^{-x}}.
\end{equation}

It is worth noting that we do not use the ReLU activation function that is commonly used in deep learning because it often involves high-order derivatives of space-time, but the second-order derivative of ReLU is zero, so PINNs often uses infinitely smooth functions as the activation function.

\subsection{Introduction to DeepONet} \label{sec: introduction to DeepONet}
Let's talk about the conclusion of DeepONet first. DeepONet is an algorithm based on neural network operator learning. The mathematical structure is similar to polynomial and periodic function fitting

\begin{equation}
	f(\boldsymbol{y})=\sum_{i=1}^{n}\alpha_{i}\phi(\boldsymbol{y}_{i}).\label{eq:basisfunction approximation}
\end{equation}

Trunk net uses the basis function $\phi(\boldsymbol{y})$, which is fitted with a neural network. The trunk net is similar to the basic idea of PINNs, which fits from the coordinate space to the objective function. However,  the basis function is fitted instead of the objective function in DeepONet. In addition, the constant weight $\alpha$ in front of the basis function is fitted by the branch net. It is worth noting that the same input function of the branch net determines a fixed weight $\alpha$, so this is similar to the function approximation idea in numerical analysis.
The traditional function approximation algorithm is to select the basis function in advance. In contrast, DeepONet uses the neural network to choose the proper basis function according to the data adaptively.
Since trunk net is similar to PINNs, we can use the automatic differential algorithm \citet{automatic_differential} in PINNs to construct a partial differential operator and obtain the loss function (strong form or energy principle) according to PDEs.
We develop a semi-supervised combination of data operators and partial differential equations, which will be discussed in detail in DCM-O \Cref{sec: DCM-o}. Next, we introduce the network structure of DeepONet.

The network structure of DeepONet is similar to the theoretical results of the article on universal operator approximation in 1995  \citep{DeepOnet, chen1995universal}
\begin{equation}
	|G(\boldsymbol{u})(\boldsymbol{y})-\sum_{k=1}^{p}\sum_{i=1}^{n}c_{i}^{k}\sigma(\sum_{j=1}^{m}w_{ij}^{k}\boldsymbol{u}(\boldsymbol{x}_{j})+b_{i}^{k})\sigma(W^{k}\cdot\boldsymbol{y}+B^{k})|<\epsilon,\label{eq: DeepONet theory}
\end{equation}
where we can use the neural network to replace the item in \Cref{eq: DeepONet theory}

\begin{equation}
	\begin{cases}
		NN_{Branch}^{k}(\boldsymbol{u}(\boldsymbol{x});\boldsymbol{\theta}_{NN}^{B}) & =\sum_{i=1}^{n}c_{i}^{k}\sigma(\sum_{j=1}^{m}w_{ij}^{k}\boldsymbol{u}(\boldsymbol{x}_{j})+b_{i}^{k})\\
		NN_{Trunk}^{k}(\boldsymbol{y};\boldsymbol{\theta}_{NN}^{T}) & =\sigma(W^{k}\cdot\boldsymbol{y}+B^{k})\\
		DeepONet(\boldsymbol{x};\boldsymbol{u}) & =\sum_{k=1}^{p}NN_{Branch}^{k}(\boldsymbol{u}(\boldsymbol{x});\boldsymbol{\theta}_{NN}^{B})\cdot NN_{Trunk}^{k}(\boldsymbol{y};\boldsymbol{\theta}_{NN}^{T}),
	\end{cases}
\end{equation}
where $\boldsymbol{u}(\boldsymbol{x})$ can be approximated from the sensor $\boldsymbol{x}$, i.e., discrete points approximate continuous function. For example, the input $\boldsymbol{u}(\boldsymbol{x})$ of $NN_{Branch}$  is the gravitational potential energy of a three-dimensional cube material ($[0,1]^{3}$). The gravitational potential energy field in $[0,1]^{3}$ takes $101^{3}$ at equal intervals of 0.01. The data structure of $\boldsymbol{u}(\boldsymbol{x})$ is a variety of different structures, such as the data structure close to the image. Therefore, we can use different network structures as the network architecture of the branch net, e.g.,  CNN with local features and time-related RNN. Noting that the network structure of branch net must consider the essence of the problem to select the appropriate network structure closest to the essence of the problem. Hence, we can find that $\boldsymbol{u}(\boldsymbol{x})$ determines the output of the branch net in DeepONet, which is related to the weight $\alpha$ in \Cref{eq:basisfunction approximation}, and has no correlation with the coordinates $\boldsymbol{y}$ we are interested in. Similarly,  $\boldsymbol{y}$ completely determines the output of trunk net, which has the same characteristics as the basis function $\phi(\boldsymbol{y})$ in \Cref{eq:basisfunction approximation}. Therefore, DeepONet has many similar ideas to traditional function approximation.

For many details and extensions about DeepONet, please refer to \citep{lu2022comprehensive}.

\subsection{Introduction to deep energy method}

The deep energy method (DEM) \citet{loss_is_minimum_potential_energy} is based on the principle of minimum
potential energy in order to obtain the true displacement solution. The principle
of minimum potential energy means that the true displacement field
minimizes the potential energy functional $J(\boldsymbol{u})$ among
all the admissible displacement fields satisfying the displacement
boundary conditions in advance. Note that the theory here is not only
for the case of linear elasticity, but also for some nonlinear hyperelastic
problems \citep{PINN_hyperelasticity}. The formulation of DEM can be written as:
\begin{equation}
\begin{aligned}
	u^{ture} & =\argmin_{u}J(u)\\
	s.t. \quad   \boldsymbol{u}(\boldsymbol{x}) & =\boldsymbol{\bar{u}}(\boldsymbol{x}),\boldsymbol{x}\subseteq\Gamma^{u},
\end{aligned}
\end{equation}
where $\boldsymbol{\bar{u}}(\boldsymbol{x})$ is the given Dirichlet
boundary $\Gamma^{u}$ condition and
\begin{equation}
J(\boldsymbol{u})=\int_{\Omega}\psi(\boldsymbol{u})dV-\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{u}dV-\int_{\Gamma^{t}}\bar{\boldsymbol{t}}\cdot\boldsymbol{u}dA,
\end{equation}
where $\psi$, $\Omega$, $f$, $\bar{t}$ and $\Gamma^{t}$ is
the strain energy density function, the domain, the body force, the
surface force and Neumann boundary respectively. To be specific, the
strain energy density in linear elasticity is
\begin{equation}
\begin{aligned}\psi & =\frac{1}{2}\boldsymbol{\varepsilon}:\boldsymbol{C}:\boldsymbol{\varepsilon}\\
	\boldsymbol{\varepsilon} & =\frac{1}{2}(\nabla\boldsymbol{u}+\boldsymbol{u}\nabla).
\end{aligned}
\end{equation}

In nonlinear elasticity, e.g., hyperelastic case, $\psi$ can be
expressed in terms of the deformation gradient $\boldsymbol{F}$.
For detailed conceptualization of the strain energy density form in
hyperelastic, readers may refer to \citep{wang2022cenn,gasser2006hyperelastic,chagnon2015hyperelastic,belytschko2013nonlinear}.

From a physical point of view, it is a necessary condition for using
the energy principle that the physical process is reversible because the
theorem generally cannot be directly applied to some irreversible
systems, e.g., elastic-plastic dissipative systems. Mathematically,
all PDEs have the corresponding weak forms obtained from the integration by parts. Only PDEs with good properties (differential operators
are linear self-adjoint) can use the variational principle to change the stationary problem to the
extreme value problem.

The idea in DEM is simple but effective. DEM uses a neural network to
replace the approximation function like the trail function in FEM.
Due to the huge function space of the neural network, DEM can have a better performance compared to the traditional methods theoretically, regardless of the optimization
error. DEM is the energy form of physics-informed neural networks
(PINNs). In the traditional PINNs, we construct the physical loss
by the strong form of the PDEs directly \citep{PINN_solid_mechanics}.
However, the differential equations in solid mechanics are high-order
tensors and its derivatives. If we use the energy form of the PINNs, the
accuracy and efficiency are commonly higher than the strong form. \Cref{fig:Schematic-of-PINN_DEM}
shows the difference and connection between the strong form of PINNs
and DEM.

In DEM, we must construct the admissible displacement in advance
before using the principle of minimum potential energy, which is a
challenge due to the reduced function space and the issue of not
being able to learn when close to the boundary if we construct it by distance
function $\boldsymbol{d}(\boldsymbol{x})$ as
\begin{equation}
	\boldsymbol{u}(\boldsymbol{x})=\boldsymbol{\bar{u}}(\boldsymbol{x})+\boldsymbol{d}(\boldsymbol{x})\odot\boldsymbol{u}_{g}(\boldsymbol{x}),
\end{equation}
where $\odot$ is the element-wise, i.e., multiplication of corresponding
elements does not change the shape, and $\boldsymbol{u}_{g}(\boldsymbol{x})$
is a general function, which can be neural networks in DEM. The distance
function is a special and currently common way to
satisfy boundary condition in advance \citep{wang2022cenn,PINN_hyperelasticity,boundary_conditions_distance_functions,lu2021physics}.
Besides, there is another easier way to satisfy boundary conditions
than using the distance constraint, i.e., the soft constraint way \citep{deep_ritz}.
The soft constraint method is usually implemented by a penalty loss,
i.e. $\beta$$\int_{\Gamma}(\boldsymbol{u}-\bar{\boldsymbol{u}})^{2}d\Gamma$.
However, it is a way that the boundary condition can not be satisfied
exactly, and the different penalty factor $\beta$ causes the different
numerical results.

\begin{figure}
	\begin{centering}
		\includegraphics[scale=0.85]{PINNandDEM_schematic}
		\par\end{centering}
	\caption{Schematic of the difference and connection between the strong form
		of PINNs and DEM, purple circles in the Neural network are the inputs
		and outputs. Blue circles in the Neural network are the hidden neurons.
		$\text{P}()$ is the energy density of the functional. $\text{Loss}_{i}$ is
		the essential boundary loss. MSE is the mean square error to let
		the interesting field to satisfy the essential boundary. If the boundary
		condition is satisfied in advance, $\text{Loss}_{i}$ can be dismissed. $\text{L}()$
		is the differential operator related to the PDEs. $\lambda_\text{PINNs}$
		and $\lambda_{i}$ are the weight of the loss of strong form in PINNs
		and $\text{Loss}_{i}$ (including initial condition if temporal problem)
		respectively. Note that the number of $\lambda_{\text{PINNs}}$ is determined
		by the number of PDEs. The dotted arrow means not all PDEs can
		be converted to the energy form.\label{fig:Schematic-of-PINN_DEM}}
\end{figure}

\subsection{Introduction to stress function}\label{sec: introduction to stress function}

Stress function, an extension of the famous stress solution of the elasticity problem in solid mechanics, is aimed to satisfy the equilibrium equation automatically. There are two main methods, displacement and stress, to solve the
elasticity problem. In fact, there is another solution called the
strain solution method, but since the strain and stress in elastic
mechanics problems are linearly related, the difficulty of the strain
solution is the same as the stress solution method. In addition, most
boundary conditions are represented by displacement or stress instead
of strain. As a result, the strain solution method is rarely used, and we only discuss the displacement and stress solution.
In the stress method, although the stress must satisfy the equilibrium equations, the stress function can automatically satisfy the equilibrium equations.
Thus, the stress function can not only reduce the number of equations
in the stress solution but also maintain the advantages of the stress
solution, so the stress function solution is very important in elastic
mechanics problems. To better understand the stress function, we briefly
introduce displacement and stress solutions of elastic mechanics. For detailed
conceptualization of methods to solve the elastic mechanic's problems,
readers may refer to \citep{the_foundation_of_solid_mechanics_feng}.

The displacement method is a solution in which displacement $\boldsymbol{u}$ is
the basic unknown quantity to be solved. However, stress is the basic
unknown quantity in stress solution. We substitute the strain tensor obtained by geometrical equations
\begin{equation}
	\boldsymbol{\varepsilon}=\frac{1}{2}(\nabla\boldsymbol{u}+\boldsymbol{u}\nabla)
\end{equation}
into the linear constitutive equations, i.e., Hooke's law
\begin{equation}
	\boldsymbol{\sigma}=\boldsymbol{C}:\boldsymbol{\varepsilon},
\end{equation}
where $\boldsymbol{u}$, $\nabla$, $\boldsymbol{\sigma}$,
$\boldsymbol{C}$ and $\boldsymbol{\varepsilon}$ is the displacement
field, gradient operator, stress tensor, tensor of the elastic constants,
and strain tensor. Thus, we can express the stress $\boldsymbol{\sigma}$ in terms
of displacement $\boldsymbol{u}$ in the equation of the equilibrium
\begin{equation}
	\nabla\cdot\boldsymbol{\sigma}+\boldsymbol{f}=\rho\ddot{\boldsymbol{u}},
\end{equation}
where $\boldsymbol{f}$ is the body force, such as gravity. $\rho$
and $\ddot{\boldsymbol{u}}$ are the density and acceleration. $\nabla\cdot()$
means divergence operator.

We solve the equilibrium equations represented by the displacement,
called Navier-Lamé equations (LN equations), which is the displacement
solution method. LN equation can be written in the form
\begin{equation}
	\begin{cases}
	Gu_{i,jj}+(\lambda+G)u_{j,ji}+f_{i}=0 & \boldsymbol{x}\in\Omega\\
	u_{i}=\bar{u}_{i} & \boldsymbol{x}\in\partial\Omega,
	\end{cases}\label{eq:the LN equations}
\end{equation}
where the comma notation ``,'' indicates partial derivatives, such
as $u_{x,x}=\partial u_{x}/\partial x$, and the subscripts $i,j$
conforms to the Einstein notation, such as $\sigma_{kk}=\sigma_{xx}+\sigma_{yy}+\sigma_{zz}$
and $u_{i}$ means $i$ can be \{x, y, z\}. $\lambda$ and $G$ are
Lame parameters
\begin{equation}
	\begin{cases}
		\lambda=\frac{\nu E}{(1+\upsilon)(1-2\upsilon)}\\
		G=\frac{E}{2(1+\upsilon)},
	\end{cases}
\end{equation}
where $E$ and $\upsilon$ represent the elastic modulus and Poisson
ratio respectively. We can get the solution of $\boldsymbol{u}$ by
solving \Cref{eq:the LN equations} with well-defined boundary conditions, which
we called the displacement method because our solution variable is
only the displacement $\boldsymbol{u}$. When the boundary conditions are all the given displacement, the displacement method is preferred over the stress method.

The stress method is a solution in which stress $\boldsymbol{\sigma}$ is the basic unknown quantity to be solved. We substitute the constitutive
equations into Saint-Venant compatibility equations
\begin{equation}
	\nabla\times\boldsymbol{\boldsymbol{\varepsilon}}\times\nabla=0\label{eq:compatibility equations},
\end{equation}
where $\nabla\times$ denotes the curl operator.
In the constitutive
equations, strain is expressed in terms of stress
\begin{equation}
	\varepsilon_{ij}=\frac{1+\upsilon}{E}\sigma_{ij}-\frac{\upsilon}{E}\sigma_{kk}\delta_{ij},\label{eq:strain-stress}
\end{equation}
where $\delta_{ij}$ is the Kronecker delta, which means $\delta_{ij}=1$, only when $i=j$ (not summed). The compatibility equations
guarantee single-valued and continuous displacement, which means there
is not cracking and overlapping physical phenomenon \citep{sadd2009elasticity}. We can rewrite \Cref{eq:compatibility equations}
in terms of stress $\boldsymbol{\sigma}$ as
\begin{equation}
	\nabla^{2}\sigma_{ij}+\frac{1}{1+\upsilon}\theta_{,ij}=-(f_{i,j}+f_{j,i})-\frac{\upsilon}{1-\upsilon}\delta_{ij}f_{k,k},\label{eq:BM equations}
\end{equation}
where $\theta=\sigma_{kk}$, which means 3 times the volumetric stress
(hydrostatic stress) changing only in volume in the physical meaning.
\Cref{eq:BM equations} are known as the Beltrami-Michell compatibility
equations (BM equations). Note that the BM equations are dependent,
which means the independent equations are less than 6 if the problem
is 3D. Because the number of the basic unknowns (stress) is 6, we
need to add the additional equations, i.e., the equilibrium equations
in terms of stress, to make the equation closed. Thus the closed equations
of elastostatics problem in stress solution can be summarized as:
\begin{equation}
	\begin{cases}
		\nabla^{2}\sigma_{ij}+\frac{1}{1+\upsilon}\theta_{,ij}+(f_{i,j}+f_{j,i})+\frac{\upsilon}{1-\upsilon}\delta_{ij}f_{k,k}=0 & \boldsymbol{x}\in\Omega\\
		\sigma_{ij,i}+f_{j}=0 & \boldsymbol{x}\in\Omega\\
		\sigma_{ij}n_{j}=\bar{t}_{i} & \boldsymbol{x}\in\partial\Omega,
	\end{cases},\label{eq:BM equtaions add equail equations}
\end{equation}
where we only consider the Neumann boundary conditions (specified
surface tractions) because the stress solution is quite difficult to
deal with displacement boundary conditions due to boundary integrals.

In order to solve the non-independence problem of the stress solution
in \Cref{eq:BM equtaions add equail equations}, we can use some tricks,
i.e., Bianchi Identity
\begin{equation}
	\begin{split}
		L_{ij,j} & =0\\
		L_{ij} & =e_{ink}e_{jml}\Phi_{mn,kl},
	\end{split}
\end{equation}
where $e_{ijk}$, the permutation symbol, is defined that $e_{ijk}=1$ when subscripts permute like 1,2,3; $e_{ijk}=-1$ when
subscripts permute like 3,2,1; $e_{ijk}=0$ when the value of any
subscript is the same. Because the characteristic of $L_{ij}$ satisfy
the similar form of equilibrium equation, the $\boldsymbol{\Phi}$
called the stress function can be used to be the basic unknown variable
instead of stress. The \Cref{eq:BM equtaions add equail equations}
without body force can be expressed by the stress function $\boldsymbol{\Phi}$
as:
\begin{equation}
	\begin{cases}
		\nabla^{2}(e_{ink}e_{jml}\Phi_{mn,kl})+\frac{1}{1+\upsilon}(e_{pnk}e_{pml}\Phi_{mn,kl}){}_{,ij}=0 & \boldsymbol{x}\in\Omega\\
		(e_{ink}e_{jml}\Phi_{mn,kl})n_{j}=\bar{t}_{i} & \boldsymbol{x}\in\partial\Omega
	\end{cases},\label{eq:stress equtaions 3d}
\end{equation}
where the stress function $\Phi_{ij}$ (symmetry) has 3 independent
components in three-dimensional elasticity. \Cref{eq:stress equtaions 3d}
above is the stress function method, which is the other form of the stress
solution. There are $C_{6}^{3}-3=17$ possible choices of the stress function.
The reason for subtracting 3 is that there are 3 combinations
that automatically satisfy \Cref{eq:stress equtaions 3d}. The most
common choices of the stress function are Maxwell stress function
and Morera stress function,
\begin{equation}
	\Phi_{Maxwell}=\left[\begin{array}{ccc}
		\Phi_{11} & 0 & 0\\
		& \Phi_{22} & 0\\
		Sym &  & \Phi_{33}
	\end{array}\right];\Phi_{Morera}=\left[\begin{array}{ccc}
		0 & \Phi_{12} & \Phi_{13}\\
		& 0 & \Phi_{23}\\
		Sym &  & 0
	\end{array}\right].
\end{equation}
To the best of our knowledge, this is the first attempt to leverage
the power of DEM based on the principle of the complementary energy
method, so we focus on the most common two-dimensional problems in
solid mechanics, i.e., the plane state of stress or strain and torsion
problems. \Cref{fig:the relation of different solution} shows the relation between displacement and stress solution.

\begin{figure}
	\begin{centering}
		\includegraphics[scale=1.8]{solution_relation}
		\par\end{centering}
	\caption{The relation between displacement and stress solution. The black solid line arrow is the derivation process;    the dotted arrow indicates that the equation is automatically satisfied \label{fig:the relation of different solution}}
\end{figure}

In the two-dimensional stress and strain problems, the Maxwell stress
function can be reduced to the Airy stress function, i.e., $\Phi_{11}=\Phi_{22}=0$,$\Phi_{33}\neq0$.
Besides, in the torsion problem of the cylindrical body, the Morera stress
function can be reduced to the Prandtl stress function, i.e., $\Phi_{12}=\Phi_{13}=0$,
$\partial\Phi_{23}/\partial x\neq0$, which is shown in \Cref{eq:Airy and prandtl}.
\begin{equation}
	\Phi_{Airy}=\left[\begin{array}{ccc}
		0 & 0 & 0\\
		& 0 & 0\\
		Sym &  & \Phi_{33}
	\end{array}\right];\Phi_{Prandtl}=\left[\begin{array}{ccc}
		0 & 0 & 0\\
		& 0 & \Phi_{23}\\
		Sym &  & 0
	\end{array}\right]\label{eq:Airy and prandtl}
\end{equation}

 Thus, to tackle the special issue,  the Airy and Prandtl stress
functions are the special and simple forms of the Maxwell and Morera
stress function, respectively. In this paper, we consider
the most common stress function in solid mechanics, i.e., Airy and Prandtl stress function,
to research the plain problems and the torsion problem by the principle
of minimum complementary energy, given the lack of deep energy method
based on the principle of minimum complementary energy. As a result,
we introduce the Airy and Prandtl stress function respectively in the next section.

\subsubsection{Airy stress function}\label{sec: Airy stress}

Airy stress function is a special and simple form of the Maxwell stress
function to tackle the plane problems. Only the stress function component
$\Phi_{33}$ of the Maxwell stress function is not zero, so we use
Airy function $\phi_{A}$ to replace the only Maxwell stress function
not equal to zero. Thus, the \Cref{eq:stress equtaions 3d} can be
simplified to
\begin{equation}
	\begin{cases}
		\nabla^{2}\nabla^{2}\phi_{A}=-(1-\upsilon)\nabla^{2}V & \boldsymbol{x}\in\Omega\\
		l(\frac{\partial^{2}\phi_{A}}{\partial y^{2}}+V)-m(\frac{\partial^{2}\phi_{A}}{\partial x\partial y})=\bar{t}_{x} & \boldsymbol{x}\in\partial\Omega\\
		-l(\frac{\partial^{2}\phi_{A}}{\partial x\partial y})+m(\frac{\partial^{2}\phi_{A}}{\partial x^{2}}+V)=\bar{t}_{y} & \boldsymbol{x}\in\partial\Omega
	\end{cases},\label{eq: The Airy equations}
\end{equation}
where $l$ and $m$ are the direction cosines of the outer normal
to the boundary curve; $\bar{t}_{x}$ and $\bar{t}_{y}$ are
the surface tractions acting on the boundary surface; $\nabla^{2}\nabla^{2}$
is a biharmonic operator, which is expressed in the Cartesian coordinate
system as :
\begin{equation}
	\nabla^{2}\nabla^{2}=\frac{\partial}{\partial x^{4}}+2\frac{\partial}{\partial x^{2}\partial y^{2}}+\frac{\partial}{\partial y^{4}}.
\end{equation}
Note $V$ is the potential of the body force and the limitation of
the Airy stress function is that the Airy stress method only can solve
the problem with the potential of the body force, and can not directly
deal with the displacement boundary conditions. If the problem has
the displacement boundary conditions, we can use Saint Venant's
theory to transform the displacement (Dirichlet) boundary into the force boundary
condition. The stress is expressed by the Airy stress function as
\begin{equation}
	\begin{cases}
		\sigma_{xx}=\frac{\partial^{2}\phi_{A}}{\partial y^{2}}+V\\
		\sigma_{yy}=\frac{\partial^{2}\phi_{A}}{\partial x^{2}}+V\\
		\sigma_{xy}=-\frac{\partial^{2}\phi_{A}}{\partial x\partial y}
	\end{cases}.
\end{equation}

The Airy stress function has the following important properties :
\begin{enumerate}
	\item A difference of a linear function in Airy stress function does not
	affect the solution of the Airy stress function. The reason is that
	the lowest derivative order of the \Cref{eq: The Airy equations} is
	two without body force, so the linear function will not affect the
	solution result.
	\item In the physical perspective, the value of the Airy stress function
	is the moment of the external force from the reference point to the
	current position at the boundary, and the directional derivative of
	the Airy stress function is the integral of the external force vector.
	This property is significant to guess the formation of the Airy stress
	function, as we explain in \Cref{fig:The-value-of stress function}a.
	\item The value of the Airy stress function is independent of the reference
	frame because the stress is objective.
\end{enumerate}
\begin{figure}
	\begin{centering}
		\includegraphics{Airy_stress_illustration}
		\par\end{centering}
	\caption{The illustration of Airy and Prandtl stress function: (a) The value
	of the Airy stress function: the reference point means $\phi_{A}=\partial\phi/\partial x=\partial\phi/\partial y=0$.
	$\bar{\boldsymbol{t}}$ is the surface force. $\boldsymbol{n}$ and
	$\boldsymbol{s}$ is the normal and tangential direction, respectively.
	$R_{x}$ ($R_{y}$) means the integral of the surface force in the
	x (y) direction at the boundary from the reference point to the current
	position B. $M$ is the moment of the external force about the current
	point at the boundary from the reference point. Counterclockwise is
	the positive direction of the moment. When counterclockwise at the boundary,
	$\phi$ is equal to $M$. However, $\phi$ is equal to $-M$ as clockwise.
	(b) The value of the Prandtl stress function: the red lines represent
	the contours of the Prandtl stress function. The value of shear stress
	is equal to the negative gradient of the Prandtl stress function.
	The direction of the shear stress points is the tangent direction
	of the contour of the Prandle stress function. The Prandtl stress
	is equal to the constant (we usually choose zero) on the boundary. \label{fig:The-value-of stress function}}
\end{figure}

\subsubsection{Prandtl stress function}

Prandtl stress function is a special form of the Morera stress function
to deal with the torsion problem of a cylindrical body, as shown in
\Cref{fig:The-value-of stress function}b. Only the stress function
component $\Phi_{23}$ of the Morea stress function is not zero, so
we use Prandtl function $\phi_{P}$ to replace the only Morea stress
function not equal to zero. Thus, the \Cref{eq:stress equtaions 3d}
can be simplified to
\begin{equation}
	\begin{cases}
		\nabla^{2}\phi_{p}=-2G\alpha & \boldsymbol{x}\in\Omega\\
		M_{t}=2\iint\phi_{P}dxdy & \boldsymbol{x}\in\Omega\\
		\phi_{P}=C & \boldsymbol{x}\in\partial\Omega
	\end{cases},\label{eq: The Airy equations-1}
\end{equation}
where $\alpha$ is the rate of twist, i.e. $\alpha=d\theta/dz$. $G$
is one of the Lame parameters, i.e., the shear modulus, and $M_{t}$
is the torque acting on the boundary condition at the ends $z=H$.
$M_{t}=2\iint\phi_{P}dxdy$ is the boundary condition at the ends of
the cylindrical. $\phi_{P}=C$ is the lateral surface boundary condition
without surface force, and no loss of generality is involved in setting
$\phi_{p}=0$ in a connected region ($\phi_{p}=0$ in a connected
region by default unless otherwise stated). Note that $\alpha$ is
unknown, so we usually solve
\begin{equation}
	\begin{cases}
		\nabla^{2}\phi_{p}=-2G\alpha & \boldsymbol{x}\in\Omega\\
		\phi_{P}=0 & \boldsymbol{x}\in\partial\Omega
	\end{cases}\label{eq: The Airy equations_at_first}
\end{equation}
first with an unknown constant $\alpha$. We can readily get $\phi_{p}$ with a unknown constant $\alpha$,
because the equation is the standard Poisson equation (there are many
well-establish methods in potential theory \citep{kellogg1953foundations}
to help us to solve it). Then we substitute $\phi_{p}$ with a unknown
constant $\alpha$ into $M_{t}=2\iint\phi_{P}dxdy$, and we can solve
for $\alpha$ eventually.
The relation between Prantl stress function and the stress in the Cartesian coodinate system can be expressed as:
\begin{equation}
	\begin{cases}
		\tau_{zx} = \frac{\partial \phi_{p}}{\partial y}   \\
		\tau_{zy}  =- \frac{\partial \phi_{p}}{\partial x}.   \\
	\end{cases}
\end{equation}

The Prantl stress function has important properties as shown
in \Cref{fig:The-value-of stress function}b. We can use the minimum complementary energy variation principle to convert it into the corresponding energy form. In the next section, we introduce the principle of minimum complementary energy.

\section{Method} \label{sec: method}
In this section, we introduce the principle of minimum complementary energy.
The proposed DCM mainly uses neural networks to fit the admissible stress function field. The construction of the admissible Airy stress function field is  more complicated than the admissible displacement field based on the principle of minimum potential energy, so we give a detailed procedure for constructing the admissible stress function field.

\subsection{Introduction to the principle of minimum complementary energy} \label{sec: minimum complementary energy}

The above displacement and stress solutions from the infinitesimal parallelepiped perspective are strong forms (differential formulation).   From the energy form, we can get the variational formulation. The traditional deep energy form is from the principle of minimum potential energy. There is another important energy principle in solid mechanics called the principle of minimum complementary energy. However, there is currently no research on deep energy algorithms for this energy principle (complementary energy).

The minimum complementary energy principle is a variational principle with stress as the basic unknown variable. On one hand, the principle of minimum potential energy is proposed by Willard Gibbs \citep{the_foundation_of_solid_mechanics_feng}. There are many practical applications and algorithms based on this principle, such as the celebrated type of displacement element in FEM \citep{finite_element_book}, the approximation theory about plates and shells with the assumption of Kirchhoff, and the vanilla\footnote{Vanilla means the unmodified version in the original paper} deep energy method \citep{loss_is_minimum_potential_energy}.
On the other hand, the concept of complementary energy was first introduced by F.Z.Engesser in 1889 \citep{complementary_energy} and was elaborated by H.M.Westergaard in 1942 \citep{principleofminimumcomplementary_energy}. Also, some influential works involve the complementary principle, such as Reissner's principle (the part of Reissner's functional is expressed using the concept of complementary strain energy) and the stress element in FEM.
The complementary energy $\Pi_{c}$ is composed of the complementary strain energy $U_{c}$ and the complementary potential $V_{c}$, and its expression is:

\begin{equation}
	\Pi_{c}=U_{c}+V_{c}.\label{eq:=006700=005C0F=004F59=0080FD=00539F=007406}
\end{equation}
The strain complementary energy $U_{c}$ is the integrated overall energy of the  complementary strain energy density $W_{c}$ on the solution domain, and its expression is:

\begin{equation}
	U_{c}=\int_{\Omega}W_{c}d\Omega
\end{equation}

The density of complementary strain energy $W_{c}$ is the energy density in the integral form of stress as the basic unknown variable, and its expression is:
\begin{equation}
	W_{c}=\int_{\sigma_{ij}}\varepsilon_{ij}d\sigma_{ij}.
\end{equation}
The complementary potential is expressed as
\begin{equation}
	V_{c}=-\int_{\Gamma^{u}}\bar{u}_{i}p_{i}d\Gamma,
\end{equation}
where its calculation method is the given essential displacement boundary condition $\boldsymbol{\bar{u}}$ multiplied by the constraint reaction force $\boldsymbol{p}=\boldsymbol{n}\cdot\boldsymbol{\sigma}$ with the stress field as the basic unknown quantity, and its physical meaning is the excess energy absorbed by the support system or transmitted to other objects by the support system.

The admissible stress field must satisfy the given force boundary conditions and equilibrium equations
\begin{equation}
	\nabla\cdot\boldsymbol{\sigma}+\boldsymbol{f}=0
\end{equation}
beforehand. The stress field can be constructed in a way similar to the Ritz method:
\begin{equation}
	\sigma_{ij}=\sigma_{ij}^{0}+\sum_{m=1}^{n}a_{m}\sigma_{ij}^{m},
\end{equation}
where $\sigma_{ij}^{0}$ is a non-homogeneous boundary condition satisfying a given force boundary condition and equilibrium equations, including body force items,
\begin{equation}
	\begin{split}
		\nabla\cdot\boldsymbol{\sigma^{0}}+\boldsymbol{f}=0\\
		\boldsymbol{n}\cdot\boldsymbol{\sigma^{0}}=\boldsymbol{\bar{t}}.
	\end{split}
\end{equation}
$\sigma_{ij}^{m}$ is a stress field that satisfies a given homogeneous force boundary at the given force boundary condition and satisfies the equilibrium equation without body force term,
\begin{equation}
	\begin{split}
		\nabla\cdot\boldsymbol{\sigma^{m}}=0\\
		\boldsymbol{n}\cdot\boldsymbol{\sigma^{m}}=0.
	\end{split}
\end{equation}
$a_{m}$ are undetermined coefficients. However, it is challenging to construct the admissible stress field in advance because the equilibrium equations are not easy to satisfy. Thus, the common method is to use the stress function method to replace the stress field because the stress function naturally satisfies the equilibrium equations where the body force $\boldsymbol{f}$ is constant. Only  the force boundary conditions instead of equilibrium equations should be considered. The expression of the admissible stress function field is

\begin{equation}
	\phi_{i}(\boldsymbol{x})=\phi_{i}^{0}(\boldsymbol{x})+\sum_{m=1}^{n}a_{i}^{m}\phi_{im}(\boldsymbol{x})\label{eq:basis function of admissible stress function},
\end{equation}
where $\phi_{i}^{0}$ is the stress function field satisfying the inhomogeneous boundary condition at the given force boundary condition, and $\phi_{im}$ is the homogeneous boundary in the given force boundary condition. $a_{i}^{m}$ is the undetermined coefficient, noting that the i index here is a free index without summation. $a_{i}^{m}$ represents the undetermined coefficient of different stress function components. We put  \Cref{eq:basis function of admissible stress function} into  \Cref{eq:=006700=005C0F=004F59=0080FD=00539F=007406}, and optimize the complementary energy functional to minimize it:
\begin{equation}
	\phi_{i}=\argmin_{a_{i}^{m}}\Pi_{c}(\phi_{i}).
\end{equation}

For nonlinear problems, the key to use the principle of minimum complementary energy is the transformation of the first Piola stress (Lagrange's stress) and the deformation gradient $\boldsymbol{F}$ \citep{nonlinear_plate_theory_comple}. If the principle of minimum potential energy can solve the nonlinear problem, then the problem must be solved by the principle of minimum complementary energy, which is proved in \ref{sec:Appendix-A.-Proof DCM_exist}; in addition, in mathematics, strain energy and strain complementary energy are essentially Legendre transformations, which we also prove in \ref{sec:Appendix-A.-Proof DCM_exist}.

\subsection{DCM: Deep energy method based on the principle of minimum complementary energy}\label{sec: DCM_method}
Unlike the minimum potential energy using the neural network to approximate the displacement field, we use a neural network to approximate the stress function. Then the complementary energy is optimized to the minimum by the principle of minimum complementary energy.
However, the problem's difficulty is that not all stress functions satisfy the conditions that the admissible stress function field must satisfy the given force boundary conditions in advance.
Therefore, if we use the principle of minimum complementary energy,  admissible stress functions that satisfy the force boundary conditions must be constructed in advance. Therefore, the key to the DCM lies in the construction of admissible stress functions. The expressions of the admissible stress functions can be written as:
\begin{equation}
	\phi_{i}(\boldsymbol{x})=NN_{p(i)}(\boldsymbol{x})+\sum_{m=1}^{n}NN^{(m)}_{g(i)}(\boldsymbol{x})*NN^{(m)}_{b(i)}(\boldsymbol{x}), \label{eq:distance funciton of the admissible stress function}
\end{equation}
where $NN_{p}$ is the particular solution network only satisfying the non-homogeneous force boundary condition as $\boldsymbol{x}$ at the force boundary condition; $NN_{g}$ is the generalized network (no special requirement like particular solution network and basis function network $NN_{b}$), and $NN_{b}$ is the basis function network associated with different stress functions, which we will explain below. The three types are analogous to boundary networks, generalized networks, and distance networks, respectively, of admissible displacement fields in \citet{PINNstrong_form_in_elastodynamics, complex_PINN_a_method_to_construct_admissible_function, admissible_in_PINN_energy_form}. It is worth noting that distance networks are a special case of basis function networks.
Generally speaking, any function can be used as the distance function as long as the distance function is zero at the given displacement boundary in the admissible displacement field. Thus, the distance function is not unique, which will be discussed in detail in \Cref{sec:basis function}. It is possible to use the different distance function to construct the admissible field to decrease the approximation error theoretically \citet{DeepOnet}.
In the Prandtl stress function, the construction of the admissible stress function is the same as the construction of the displacement field, and the difference between the construction of the admissible displacement and stress function is the different types of boundary we should consider. The admissible displacement field should consider the Dirichlet boundary condition (the displacement boundary), while the admissible stress function should consider the Neumann boundary condition (the force boundary).
But in the Airy stress function, not only the basis function is zero at the given force boundary, but also the normal derivative of the basis function is required to be zero at the force boundary, which is shown in \ref{sec:Appendix-C.the basis function in Airy stress function}. Compared with the admissible displacement field, the basis function in the Airy stress function has the extra condition that the normal derivative is zero, mainly due to the relationship between the stress function and the stress. As a result, the basis function is related to the different stress function.
\begin{rem}
	Theoretically, there are many ways of constructing admissible stress function in theory, such as \Cref{eq:basis function of admissible stress function} (Coefficient function method) and \Cref{eq:distance funciton of the admissible stress function} (Basis function method). \Cref{eq:basis function of admissible stress function} is to construct many admissible basis functions fields (such as the orthogonal polynomials, i.e., Legendre polynomials, Chebyshev polynomials, periodic polynomials) satisfying the homogeneous boundary conditions in the given force boundary condition, then fit the constants in front of the basis functions. \Cref{eq:distance funciton of the admissible stress function} is based on the idea of distance function, but it needs to meet the condition that the normal derivative is zero additionally in Airy stress function, as explained in \ref{sec:Appendix-C.the basis function in Airy stress function}. Of course,  the "omnipotent" penalty function method can also be used to construct the admissible stress function. \Cref{tab:The-different-methods} shows the different methods for constructing admissible stress functions.
\end{rem}

\begin{table}
	\caption{The different methods for constructing admissible stress function in Airy stress function.
		$\Gamma^{t}$ means the force boundary condition. Homogeneous force
		boundary condition (HFB) means when $\boldsymbol{t}(\boldsymbol{x})=0$
		when $\boldsymbol{x}\in\Gamma^{t}$, and non-homogeneous force boundary
		condition (non-HFB) means when $\boldsymbol{t}(\boldsymbol{x})=\bar{\boldsymbol{t}}$
		when $\boldsymbol{x}\in\Gamma^{t}$ \label{tab:The-different-methods}}
	\begin{centering}
		\resizebox{\textwidth}{15mm}{
			\begin{tabular}{cccc}
				\toprule
				\multirow{2}{*}{Methods in Airy stress function} & \multicolumn{3}{c}{$\phi(\boldsymbol{x})=A(\boldsymbol{x})+B(\boldsymbol{x})*C(\boldsymbol{x})$}\tabularnewline
				\cmidrule{2-4} \cmidrule{3-4} \cmidrule{4-4}
				& $A(\boldsymbol{x})$ & $B(\boldsymbol{x})$ & $C(\boldsymbol{x})$\tabularnewline
				\midrule
				Coefficient function method & \multirow{3}{*}{Satisfy non-HFB  ($\boldsymbol{x}\in\Gamma^{t}$)} & Constant trainable parameters & Satisfy HFB when $\boldsymbol{x}\in\Gamma^{t}$\tabularnewline
				Basis function method &  & $B(\boldsymbol{x})=0;\partial B(\boldsymbol{x})/\partial\boldsymbol{n}=0$
				($\boldsymbol{x}\in\Gamma^{t}$) & Neural network\tabularnewline
				Penalty function method & \multicolumn{3}{c}{$\mathcal{L}_{p}=\beta_{1}[\hat{\phi}(\boldsymbol{x})-\phi(\boldsymbol{x})]^{2}+\beta_{2}[\partial\hat{\phi}(\boldsymbol{x})/\partial\boldsymbol{n}-\partial\phi(\boldsymbol{x})/\partial\boldsymbol{n}]^{2}$}\tabularnewline
				\bottomrule
		\end{tabular}}
		\par\end{centering}
\end{table}

The algorithm flow of DCM is:

\begin{enumerate}
	\item According to the force boundary conditions, the particular solution network is trained, and training points are $\{\boldsymbol{x}_{i}\in\Gamma^{t},\bar{\boldsymbol{t}}_{i} \}_{i=1}^{n}$, so that the particular solution network can fit the force boundary condition. For example, the force boundary condition expressed by the Airy stress function for the plane stress problem is:
	\begin{equation}
		\begin{split}n_{x}(\frac{\partial^{2}NN_{p}}{\partial y^{2}}+V)-n_{y}(\frac{\partial^{2}NN_{p}}{\partial x\partial y}) & =\bar{t}_{x}\\
			-n_{x}(\frac{\partial^{2}NN_{p}}{\partial x\partial y})+n_{y}(\frac{\partial^{2}NN_{p}}{\partial x^ {2}}+V) & =\bar{t}_{y},
		\end{split}\label{eq: particular equations}
	\end{equation}
	and the property of the Airy stress function can be used to construct the particular solution network due to the lower derivative order (increase the accuracy and efficiency of the construction of the particular net), i.e., $\phi=M$ and $\partial\phi/\partial\boldsymbol{n}=-R_{s}$, as shown in \Cref{fig:The-value-of stress function} and proof in \ref{sec:Appendix-B.-proof of the particular solution network}.
	If it is a free torsion problem, the particular solution network must ensure that the Prandtl stress function at the boundary conditions is constant. It is convenient to set the particular solution network $NN_{p}$ to be zero when dealing with the connected region. Noting that through the properties of the Airy stress function in \Cref{sec: Airy stress}, it is feasible to obtain the particular solution, which is equivalent to \Cref{eq: particular equations},  as shown in proof  \ref{sec:Appendix-B.-proof of the particular solution network}.

	\item According to the force boundary conditions, when dealing with the Airy stress function, we can choose one of the methods for constructing the admissible stress function in \Cref{tab:The-different-methods}.  In the case of Prandtl stress function, the training points are $\{\boldsymbol{x}_{i}\in\Gamma^{t},\phi|_{\boldsymbol{x}_{i} }=0\}_{i=1}^{n}$ in Prandtl stress function. The training method is the same as the particular solution network, but the labels are different because the basis function network satisfies the homogeneous boundary conditions.

	\item According to the principle of minimum complementary energy, the parameters of the particular solution network and the basis function network are fixed, and the generalized network $NN_{g}$ is trained. The training ends when the loss function (complementary energy) converges.
\end{enumerate}

It is worth noting that, according to the minimum complementary energy principle, the admissible stress function needs to satisfy the force boundary condition in advance; according to the minimum potential energy principle, the admissible displacement field in the deep energy method needs to satisfy the displacement boundary condition in advance. Therefore, in theory, if displacement boundary conditions mainly dominate the boundary conditions of this problem, it would be more appropriate to use the principle of minimum complementary energy because there are fewer constraints on the admissible stress function field than the admissible displacement field in the minimum potential energy principle, so the powerful fitting ability of the neural network can be released because the form of the admissible function field through the distance function will limit the function space of the neural network to a certain extent. In addition, if force boundary conditions dominate the problem, it will be more advantageous to use the principle of minimum potential energy because the assumptions on the displacement function are reduced.

\subsection{DCM-P: DCM with biharmonic function in Airy stress function}

In the Airy stress function, the strong form of the Airy stress function satisfies the biharmonic equation

\begin{equation}
	\nabla^{4}\phi=0. \label{eq:the biharmonic equation_first}
\end{equation}
We can add some terms $\phi_{bh}$ naturally satisfying \Cref{eq:the biharmonic equation_first}, such as $x^{2}$, $y^{2}$. The trainable parameters $a_{i}$ before these terms can be determined by the loss function (the complementary energy $\Pi_{c}$), i.e.

\begin{equation}
	\begin{split}
	\phi^{true}  =&\argmin_{a_{i}, \boldsymbol{\theta}_{NN^{(m)}_{g}}}\Pi_{c}(\phi)\\
	\phi(\boldsymbol{x})  =&NN_{p}(\boldsymbol{x}; \boldsymbol{\theta}_{NN^{(m)}_{p}})+\sum_{m=1}^{n}NN^{(m)}_{g}(\boldsymbol{x};\boldsymbol{\theta}_{NN^{(m)}_{g}})*NN^{(m)}_{b}(\boldsymbol{x};\boldsymbol{\theta}_{NN^{(m)}_{b}}) + \sum_{i=1}^{s}a_{i}*\phi_{bh}^{i}.
	\end{split}
\end{equation}
The result shows that DCM-P can converge faster and be more accurate than DCM, which is shown in \Cref{sec: circular}. Because the algorithmic flow of the DCM is similar to DCM-P, we only conclude the algorithm of DCM-P in \Cref{alg:DCM-P}.

\begin{algorithm}
	\caption{The algorithm of DCM-P}\label{alg:DCM-P}
	\begin{algorithmic}[1]
		\State \textbf{Step 1}: Training the particular net $NN_{p}$.
		\State  Allocate training points on the force boundary conditions, i.e.,   $\{\boldsymbol{x}_{i}\in\Gamma^{t},\bar{\boldsymbol{t}}_{i} \}_{i=1}^{n}$.
		\If{Prantl stress function: }
		\State $NN_{p} = 0 $ (connected region).
		\EndIf
		\If{Airy stress function:}
		\State $Loss_{p} =
		\sum_{i=1}^{n} [l(\partial^{2}NN_{p}(\boldsymbol{x}_{i})/\partial y^{2})-m(\partial^{2}NN_{p}(\boldsymbol{x}_{i}) / \partial x\partial y)- \bar{t}_{xi}]^{2}
		+\sum_{i=1}^{n} [-l(\partial^{2}NN_{p}(\boldsymbol{x}_{i}) / \partial x\partial y)+m(\partial^{2}NN_{p}(\boldsymbol{x}_{i})/\partial x^{2})- \bar{t}_{yi}]^{2}$.
		\State Minimize $Loss_{p}$ to get the optimal parameters $\boldsymbol{\theta}_{p}$ of $NN_{p}$.
		\EndIf
		\State \textbf{Step 2}: Training the basis net $NN_{b}$.
		\State  Allocate training points $\{\boldsymbol{x}_{i}\in\Gamma^{t}, d_{i}\}_{i=1}^{m}$ on the domain and boundary, where $d_{i}$ is the nearest distance from $\boldsymbol{x}_{i}$ to the force boundary conditions.
		\If{Prantl stress function: }
		\State $Loss_{b} =
		\sum_{i=1}^{m}[NN_{b}(\boldsymbol{x}_{i})- d_{i}]^{2}$.
		\EndIf
		\If{Airy stress function:}
		\State $Loss_{b} =
		\sum_{i=1}^{m}[NN_{b}(\boldsymbol{x}_{i})- d_{i}]^{2}
		+\sum_{i=1}^{r}[NN_{b}(\boldsymbol{x}_{i})]^{2}
		+\sum_{i=1}^{r}[\partial NN_{b}(\boldsymbol{x}_{i}) / \boldsymbol{n}]^{2}$, where $r$ is the number of points on the force boundary condition and $\boldsymbol{n}$ is the normal direction of the force boundary.
		\EndIf
 		\State Minimize $Loss_{b}$ to get the optimal parameters $\boldsymbol{\theta}_{b}$ of $NN_{b}$.
		\State \textbf{Step 3}: Select the terms of biharmonic function $\phi_{bh}$ in Airy stress function.
		\State \textbf{Step 4}: Training the general net $NN_{g}$.
		\State $\phi(\boldsymbol{x})  =NN_{p}(\boldsymbol{x})+NN_{g}(\boldsymbol{x};\boldsymbol{\theta}_{NN_{g}})*NN_{b}(\boldsymbol{x}) + \sum_{i=1}^{s}a_{i}*\phi_{bh}^{i}$.
		\State Minimize the complementary energy and get the optimal parameters of $\boldsymbol{\theta}_{NN_{g}}$ and $a_{i}$ (freeze the $NN_{p}$ and $NN_{b}$).
	\end{algorithmic}
\end{algorithm}

\subsection{DCM-O: Deep operator energy method based on the principle of minimum complementary energy } \label{sec: DCM-o}
This idea's starting point is to reuse the existing calculation results fully and not waste the previous calculation results to improve the computational efficiency of DCM.
The core idea of DCM-O is to use the DeepONet framework to help DCM make operator learning, as mentioned in \Cref{sec: introduction to DeepONet}. Trunk net can build partial differential operators to construct the physical loss.
The algorithmic flow of DCM-O is shown in \Cref{alg:DCM-O}

\begin{algorithm}
	\caption{The algorithm of DCM-O}\label{alg:DCM-O}
	\begin{algorithmic}[1]
		\State \textbf{Training:}
		\State \textbf{Data:} The hide-fidelity stress function
		$\{\boldsymbol{x}^{j}_{s},\phi_{s}^{j} \}_{s=1}^{n}$
		 and variable field
		 	$\{\boldsymbol{x}^{j}_{i},f^{j}_{i} \}_{i=1}^{m}$ such as geometry, material information as the input of the branch net and $j=1,2,\cdots, t$.
		\State \textbf{Step 1}: Training the branch $NN_{br}$ and trunk net $NN_{t}$.
		\State  input the variable field
		$\{\boldsymbol{x}^{j}_{i},f^{j}_{i} \}_{i=1}^{m}$ to the  branch net;
		input the $\boldsymbol{x}^{j}_{s}$ to the  trunk net and
		$\phi_{s}^{j}$ to the output stress function.
		\State Minimize the operation loss.
		\State \textbf{Test:}
		\State \textbf{Step 2}: Training the particular net $NN_{p}$ as DCM.
		\State \textbf{Step 3}: Training the basis net $NN_{b}$
	    as DCM.
		\State \textbf{Step 4}: Training the branch $NN_{br}$ and trunk net $NN_{t}$ again.
		\State Minimize the complementary energy to obtain the optimal parameters $\boldsymbol{\theta}_{br}$ and $\boldsymbol{\theta}_{t}$ (freeze the particular and basis net).
	\end{algorithmic}
\end{algorithm}

\Cref{fig:DCM-O} shows the framework of DCM-O. We first construct the particular solution and basis function network at the force boundary condition, which is the same as the steps of DCM. We only change the general network to the network structure of DeepONet. The training method is first to learn the operator in the existing data, i.e., not to construct the physical-based complementary energy functional first. It is worth noting that if the force boundary conditions remain unchanged, the particular net and basis net freeze and no longer participate in subsequent training. If the force boundary conditions are changed, removing the particular net and basis net is more convenient, and only the branch net and trunk net still exists. We use the penalty function method to satisfy the force boundary conditions instead of the particular and basis net if the boundary condition changes. The particular net and basis net can be considered a way to satisfy the boundary condition, and the penalty function method is also a special method to satisfy the boundary condition. After training based on the existing data, we consider optimization based on the principle of complementary energy based on physical laws. The input of our data in branch net here can be different geometric shapes, materials, and force boundary conditions (other field variables related to the problem can also be inputted to branch net), and the output is the existing high-fidelity stress function result, such as the stress-function element method. Finally, we optimize the parameters of the trunk and branch net based on the principle of minimum complementary energy.

\begin{figure}
	\begin{centering}
		\includegraphics[scale=0.75]{DCM_schematic}
		\par\end{centering}
	\caption{The schematic of deep operator energy method based on the principle of minimum complementary energy (DCM-O) \label{fig:DCM-O}}
\end{figure}

It is worth noting that here we can also refer to the hybrid optimization \citet{goswami2022physics}, which combines data loss and physical loss with different hyperparameters into the same loss function for optimization together. However, in this paper, we choose not to use hybrid optimization but sequential optimization instead, i.e., optimize data loss first and finally optimize physical loss.

\section{Result\label{sec:Result}}

\subsection{Prantl stress function: full natural boundary conditions}\label{sec: prantl}
A cylindrical rod is a three-dimensional geometry obtained by translating the centroid of a plane figure along the Z-axis. It is a common engineering component that can withstand different loads such as tension, compression, bending, and torsion. The free torsion problem of cylindrical rods is a common problem in engineering \citep{the_foundation_of_solid_mechanics_feng}, which is often used to transmit torque, as shown in \cref{fig:cylinder}a.

\begin{figure}
	\begin{centering}
		\includegraphics{cylindrical}
		\par\end{centering}
	\caption{The schematic of the problem of torsion of a cylindrical body (a) and square cross-section of the cylinder (b) \label{fig:cylinder}}
\end{figure}

Two common solutions are introduced in {\Cref{sec: introduction to stress function}, namely the displacement solution and the stress function solution. The strong form expression of the displacement solution is:
	\begin{equation}
		\begin{cases}
			\frac{\partial^{2}\psi}{\partial^{2}x}+\frac{\partial^{2}\psi}{\partial^{2}y}=0 & \boldsymbol{x}\in\Omega\\
			\text{\ensuremath{\frac{\partial\psi}{\partial x}n_{x}}+\ensuremath{\frac{\partial\psi}{\partial y}n_{y}}=y\ensuremath{n_{x}}-x\ensuremath{n_{y}}} & \boldsymbol{x}\in\partial\Omega \\
			\alpha G\oiint(x^{2}+y^{2}+x\frac{\partial\psi}{\partial y}-y\frac{\partial\psi}{x})dxdy=M_{t},
		\end{cases}
	\end{equation}
	where $\psi$ is the warping function and $\alpha$ is rate of twist. The displacement function of the free torsion of the cylindrical rod based on the Saint-Venant assumption can be written as:
	\begin{equation}
		\begin{split}
			u &=-\alpha zy\\
			v &=\alpha zx\\
			w  &=\alpha\psi(x,y)\label{eq:=0067F1=005F62=006746=004F4D=0079FB=0051FD=006570}.
		\end{split}
	\end{equation}
	We usually solve the Laplace equation with Neumann boundary conditions first. Then we substitute the $\psi$  into the end boundary condition $\alpha G\iint(x^{2}+y^{2}+x \partial\psi / \partial y - y \partial\psi /\partial x)dxdy=M_{t}$ to get the unknown constant $\alpha$.

	The energy form of the displacement solution based on the principle of minimum potential energy (DEM) is:

	\begin{equation}
		\Pi=L\int_{\Omega}(\frac{1}{2}\gamma_{zx}\tau_{zx}+\frac{1}{2}\gamma_{zy}\tau_{zy})dxdy-\alpha LM_{t},
	\end{equation}
	where $L$ is the length of the cylinder. We substitute the displacement expressed as \Cref{eq:=0067F1=005F62=006746=004F4D=0079FB=0051FD=006570} into the potential, and obtain
	\begin{equation}
		\Pi=L\int_{\Omega}\frac{1}{2}G\alpha^{2}[(\frac{\partial\psi}{\partial x}-y)^{2}+(\frac{\partial\psi}{\partial y}+x)^{2}]dxdy-\alpha LM_{t}.
	\end{equation}
	Remove the constant term $\alpha LM_{t}$, approximating $G\alpha^{2}/2$. The above optimization problem is transformed into:
	\begin{equation}
		\psi=\argmin_{\psi}\{\int_{\Omega}[(\frac{\partial\psi}{\partial x}-y)^{2}+(\frac{\partial\psi}{\partial y}+x)^{2}]dxdy\}.
	\end{equation}
	It is worth noting that the warping function $\psi$ will have an undetermined constant in the free torsion problem of the cylindrical rod. This undetermined constant is the displacement of the rigid body. If it involves strain and the derivative of the displacement function, this constant will not have an impact.

	The strong form of the PINNs stress function solution is (considering simply connected regions):
	\begin{equation}
		\begin{cases}
			\frac{\partial^{2}\phi}{\partial^{2}x}+\frac{\partial^{2}\phi}{\partial^{2}y}=-2G\alpha & \boldsymbol{x}\in\Omega\\
			\phi=0 & \boldsymbol{x}\in\partial\Omega\\
			M_{t}=2\int_{\Omega}\phi d\Omega.
		\end{cases}
	\end{equation}

	In the stress function solution, the $\alpha$ is unknown in advance, so it cannot be solved directly. Here, the pre-assumption method is used, assuming $C_{1}=-2G\alpha$, to solve:

	\begin{equation}
		\begin{cases}
			\frac{\partial^{2}\phi}{\partial^{2}x}+\frac{\partial^{2}\phi}{\partial^{2}y}=C_{1}  & \boldsymbol{x}\in\Omega\\
			\phi=0  & \boldsymbol{x}\in\partial\Omega.
		\end{cases}
	\end{equation}

	The solution result is recorded as $\phi_{1}$, and $\phi_{1}$ is brought into the torque formula $M_{1}=2\int_{\Omega}\phi_{1}d\Omega$ to correct
	\begin{equation}
		\phi=\frac{M_{t}}{M_{1}}\phi_{1} \label{eq: revised fai}
	\end{equation}

	According to the uniqueness of the elastic mechanics' solution, the above revised $\phi$ in \Cref{eq: revised fai} is the true solution to the problem. Bring the corrected $\phi$ into the domain equation, and correct $C_{1}$
	\begin{equation}
		\frac{\partial^{2}\phi}{\partial^{2}x}+\frac{\partial^{2}\phi}{\partial^{2}y}=\frac{M_{t}}{M_{1}}(\frac{\partial^{2}\phi_{1}}{\partial^{2}x}+\frac{\partial^{2}\phi_{1}}{\partial^{2}y})=\frac{M_{t}}{M_{1}}C_{1}.
	\end{equation}
	Thus $\alpha$ can be written as:
	\begin{equation}
		\alpha=-\frac{M_{t}C_{1}}{2M_{1}G}.
	\end{equation}

	We consider DCM, i.e., the minimum complementary energy form of the stress function solution, which can be written as (considering simply connected regions):
	\begin{equation}
		\phi=\argmin_{\phi}\{\Pi_{c}=\frac{L}{G}\int_{\Omega}\frac{1}{2}[(\frac{\partial\phi}{\partial x})^{2}+(\frac{\partial\phi}{\partial y})^{2}]dxdy-\alpha L\int_{\Omega}2\phi dxdy\}.
	\end{equation}
	The minimum complementary energy solution also encounters the same problem as the strong form of the stress function ($\alpha$ is not known in advance). Since the first-order variation of the complementary energy is equivalent to the strong form, the same approach can be adopted, setting $\alpha=\alpha_{1}$, to solve the minimum complementary energy form, and bring the obtained $\phi_{1}$ into the torque formula $M_{1}=2\int_{\Omega}\phi_{1} d\Omega$, and we correct
	\begin{equation}
		\begin{split}
			\phi & =\frac{M_{t}}{M_{1}}\phi_{1}\\
			\alpha & =\frac{M_{t}}{M_{1}}\alpha_{1}.
		\end{split}
	\end{equation}

	We take a rectangular cylinder as an example in \Cref{fig:cylinder}b to show the comparison between the prediction and analytical solution of the torsion angle $\alpha$  of the above four methods (PINNs strong form of displacement, DEM, PINNs strong form of stress function, and DCM). The four methods' neural network structure, optimization scheme and point allocation strategy are the same. The neural network structure is the x and y coordinates of 2 neurons in the input layer, 2 hidden layers, 20 neurons in each layer, the output is a neuron, the activation function is  {\rm tanh}, the output layer without activated function; the optimization scheme is Adam; the initialization method is Xavier; all the internal point allocation methods of the four methods are randomly distributed, 10,000 points and the strong displacement form needs to allocate additional boundary points on the boundary, 1000 points per side, also randomly distributed. \Cref{fig: Comparison of the numerical solution and analytical solution of the rotation angle}  show that, for different $a/b$ aspect ratios, the four methods converge to the exact solution when increasing with the number of iterations, which shows that the four methods can predict the rotation angle $\alpha$ of unit length. It is worth noting that the DEM does not require to assume an admissible displacement field because all boundaries are force boundary conditions; However, the DCM based on the minimum complementary energy principle needs to assume an admissible stress function. The stress function of the boundary is assumed to be a constant zero. Due to the regular geometry of the problem, the admissible stress function of the boundary condition is satisfied by the coordinate construction method in advance:
	\begin{equation}
		\phi=(x^{2}-\frac{a^{2}}{4})(y^{2}-\frac{b^{2}}{4})NN(x,y;\boldsymbol{\theta}).
	\end{equation}
	The loss function of the strong form of displacement in PINNs is
	\begin{equation}
		\begin{split}\mathcal{L}_{strongdis}&=  \lambda_{1}\frac{1}{N_{dom}}\sum_{i=1}^{N_{dom}}|\frac{\partial^{2}\psi(\boldsymbol{x}_{i})}{\partial^{2}x}+\frac{\partial^{2}\psi(\boldsymbol{x}_{i})}{\partial^{2}y}|^{2}\\
			&+\lambda_{2}\frac{1}{N_{b}}\sum_{i=1}^{N_{b}}|\frac{\partial\psi(\boldsymbol{x}_{i})}{\partial x}n_{x}+\frac{\partial\psi(\boldsymbol{x}_{i})}{\partial y}n_{y}-y_{i}n_{x}+x_{i}n_{y}|^{2},\end{split}
	\end{equation}
	where we set $\lambda_{1}=\lambda_{2}=1$. The strong form of the PINNs stress function does not use a penalty function to satisfy the boundary conditions but uses an admissible stress function consistent with the DCM to satisfy the force boundary conditions, so there are no additional hyperparameters in the strong form of PINNs stress function and DCM. The strong form of PINNs stress function and  DCM both assume $\alpha=0.0005$.

	\begin{figure}
		\begin{centering}
			\includegraphics[scale=2.0]{Prantl_alpha_all_methods}
			\par\end{centering}
		\caption{ Comparison of the numerical solution and analytical solution of the rotation angle $\alpha$ in the different rectangular cylinder (aspect ratio $a/b$) under the four models of PINNs displacement strong form, DEM based on the minimum potential energy, PINNs stress function strong form and DCM based on the minimum complementary energy. \label{fig: Comparison of the numerical solution and analytical solution of the rotation angle}}
	\end{figure}

	In order to accurately compare the accuracy of the four methods, we analyze more examples, different aspect ratios  $a/b$ correspond to different exact solutions $\alpha$ and maximum shear stress $\tau_{max} $, as shown in the table \ref{tab:=005728=004E0D=00540C=00957F=005BBD=006BD4=004E0B=005355=004F4D=00957F=005EA6=00626D=0089D2=007684=0089E3=006790=0089E3}. The analytical expressions of $\alpha$ and $\tau_{max}$ are:
	\begin{equation}
		\begin{split}
			\alpha =\frac{M_{t}}{\beta Gab^{3}} \\
			\tau_{max}  =\frac{M_{t}}{\beta_{1}ab^{2}}
		\end{split}
	\end{equation}

	\begin{table}
		\caption{Analytical solutions of  $\alpha$ and the maximum shear stress $\tau_{max}$ of a rectangular cylinder under different $a/b$ aspect ratios\label{tab:=005728=004E0D=00540C=00957F=005BBD=006BD4=004E0B=005355=004F4D=00957F=005EA6=00626D=0089D2=007684=0089E3=006790=0089E3}}

		\centering{
		\begin{tabular}{ccccc}
			\toprule
			Aspect ratio $a/b$ & $\beta$ & Torsion angle $\alpha$ & $\beta_{1}$ & $\tau_{max}$\tabularnewline
			\midrule
			1.0 & 0.141 & 70.922 & 0.208 & 48076.923\tabularnewline
			1.2 & 0.166 & 50.201 & 0.219 & 38051.750\tabularnewline
			1.5 & 0.196 & 34.014 & 20231 & 28860.029\tabularnewline
			2.0 & 0.229 & 21.834 & 0.246 & 20312.203\tabularnewline
			2.5 & 0.249 & 16.064 & 0.258 & 15503.876\tabularnewline
			3.0 & 0.263 & 12.672 & 0.267 & 12484.395\tabularnewline
			4.0 & 0.281 & 8.897 & 0.282 & 8865.248\tabularnewline
			5.0 & 0.291 & 6.873 & 0.291 & 6872.852\tabularnewline
			10.0 & 0.312 & 3.205 & 0.312 & 3205.128\tabularnewline
			\bottomrule
		\end{tabular}
	\end{table}

	\Cref{fig:Prantl_alpha_tau_inital}a gives the relative error $|pred-exact|/exact$ of $\alpha$ under different aspect ratio $a/b$. The number of iterations of all methods is 2000, and all configurations related to training are the same, the results show that under the same number of iterations, the accuracy of the four methods is similar. From the perspective of computational efficiency, the minimum potential energy of DEM and the minimum complementary energy of DCM can be computationally more efficient due to the lower order of the derivative than the corresponding strong form of displacement and stress function theoretically. Therefore, all four methods can obtain satisfactory results under different rectangular sizes. Since all boundaries are force boundary conditions, DEM does not require additional construction of admissible displacement fields in advance, which is a great advantage compared to the DCM. From the other extreme, if all the boundaries are displacement boundary conditions which we show in \Cref{sec: circular}, the DCM does not need to construct the admissible stress function field in advance. Herefore, different energy principles have different advantages under different types of boundary conditions.
	\begin{figure}
		\begin{centering}
			\includegraphics[scale=1.0]{Prantl_alpha_tau_inital}
			\par\end{centering}
		\caption{(a) Relative error of rotation angle $\alpha$ of a rectangular cylinder under four models under PINNs displacement strong form, DEM, PINNs stress function strong form and DCM. (b) The effect of different initial $\alpha$ on the convergence rate under DCM. (c) Relative error of the maximum shear stress $\tau$ of the rectangular cylinder under the strong form of PINNs displacement, DEM, the strong form of PINNs stress function and DCM.\label{fig:Prantl_alpha_tau_inital}}
	\end{figure}

	In order to further analyze the potential of the four methods, we cancel the limitation of the fixed number of iterations and calculate the relative error and convergence time of $\alpha$ when convergence, as shown in the table \Cref{tab:=0056DB=0079CD=0065B9=006CD5=005728=004E0D=00540C=00957F=005BBD=006BD4=004E0B=005355=004F4D=00957F=005EA6=00626D=0089D2=007CBE=005EA6=00548C=006536=00655B=0065F6=0095F4=007684=005BF9=006BD4}. The computer hardware for all the numerical experiments is RTX2060, 6G GPU, and 16G RAM. Generally the accuracy of all four methods is good. Among them, the DCM and the DEM have higher calculation efficiency because of the low order of derivatives. In addition, in terms of accuracy, the DCM is the highest.
	\begin{center}
		\begin{table}
			\caption{Comparison of accuracy and convergence time of torsion angle $\alpha$ under different aspect ratios $a/b$ under the four models of PINNs displacement strong form, DEM, PINNs stress function strong form and DCM\label{tab:=0056DB=0079CD=0065B9=006CD5=005728=004E0D=00540C=00957F=005BBD=006BD4=004E0B=005355=004F4D=00957F=005EA6=00626D=0089D2=007CBE=005EA6=00548C=006536=00655B=0065F6=0095F4=007684=005BF9=006BD4}}

			\centering{
			\begin{tabular}{ccccc}
				\toprule
				\multirow{2}{*}{Aspect ratio $a/b$} & \multicolumn{4}{c}{Torsion angle $\alpha$: relative error and convergence time}\tabularnewline
				\cmidrule{2-5} \cmidrule{3-5} \cmidrule{4-5} \cmidrule{5-5}
				& PINNs strong form: Dis & DEM & PINNs strong form: Stress & DCM\tabularnewline
				\midrule
				1.0 & 1.51\%, 7.9s & 0.73\%, 14.4s & 0.52\%, 13.4s & 0.73\%, 16.2s\tabularnewline
				1.2 & 0.34\%, 16.1s & 0.58\%, 7.1s & 0.07\%, 20.8s & 0.20\%, 13.7s\tabularnewline
				1.5 & 1.26\%, 7.8s & 0.31\%, 6.9s & 0.04\%, 22.2s & 0.23\%, 13.7s\tabularnewline
				2.0 & 1.45\%, 16.2s & 0.23\%, 11.5s & 0.55\%, 20.4s & 0.07\%, 13.7s\tabularnewline
				2.5 & 1.27\%, 16.1s & 0.65\%, 14.1s & 0.66\%, 25.7s & 0.24\%, 14.0s\tabularnewline
				3.0 & 1.14\%, 16.1s & 1.11\%, 17.4s & 0.06\%, 12.5s & 0.13\%, 13.6s\tabularnewline
				4.0 & 0.23\%, 31.6s & 1.32\%, 13.2s & 0.53\%, 38s & 0.41\%, 18.3s\tabularnewline
				5.0 & 0.84\%, 31.2s & 1.63\%, 17.1s & 1.41\%, 51.28s & 0.62\%, 18.2s\tabularnewline
				10.0 & 0.71\%, 80.1s & 2.25\%, 17.1s & 1.44\%, 102.4s & 0.32\%, 18.3s\tabularnewline
				Mean error & $0.97\pm0.44\%$ & $0.97\pm0.62\%$ & $0.59\pm0.50\%$ & $\boldsymbol{0.33\pm0.21\%}$\tabularnewline
				Mean time & $24.8\pm21.1s$ & $\boldsymbol{13.2\pm3.8s}$ & $34.1\pm26.8s$ & $15.5+2.1s$\tabularnewline
				\bottomrule
			\end{tabular}
		\end{table}
		\par
	\end{center}

	\Cref{fig:=0077E9=005F62=0067F1=005F62=006746=0081EA=007531=00626D=008F6Ctaumax=005BF9=006BD4} gives four methods to predict the maximum shear stress $\tau$,  and the maximum shear stress effect of the four methods is good. The PINNs strong form of displacement and the DEM fluctuate greatly, but the PINNs strong form of stress function and the DCM fluctuate less because the pre-setting $\alpha$ affects the result of the stress solution, including DCM and the strong form of stress function, i.e., the initial prediction value will affect the accuracy and convergence speed. All four methods can predict maximum shear stress well, as is shown in \Cref{fig:Prantl_alpha_tau_inital}c under different $a/b$.

	\begin{figure}
		\begin{centering}
			\includegraphics[scale=2.0]{Prantl_stress_all_methods}
			\par\end{centering}
		\caption{Comparison of the numerical and analytical solutions of the maximum shear stress $\tau$ of the rectangular cylinder under the strong form of PINNs displacement, DEM, the strong form of PINNs stress function and DCM\label{fig:=0077E9=005F62=0067F1=005F62=006746=0081EA=007531=00626D=008F6Ctaumax=005BF9=006BD4}}
	\end{figure}

	DEM has been discussed in detail in \citep{loss_is_minimum_potential_energy}, so here we mainly discuss the DCM based on the minimum complementary energy form, and we analyze the  DCM separately in this paper. \Cref{fig:=0077E9=005F62=0067F1=005F62=006746=0081EA=007531=00626D=008F6Ctau=004E91=0056FE}c,f show that the maximum absolute error is mainly concentrated on the boundary, because the maximum shear stress is concentrated on the boundary, so it cannot show that the learning of the boundary is worse than that of the interior simply. On the whole, the predictive shear stress effect of DCM is good.

	\begin{figure}
		\begin{centering}
			\includegraphics[scale=1.0]{DCM_prediction_prantl}
			\par\end{centering}
		\centering{}\caption{Comparison of numerical and analytical solution of the shear stress in cylinder $\tau$ in DCM \label{fig:=0077E9=005F62=0067F1=005F62=006746=0081EA=007531=00626D=008F6Ctau=004E91=0056FE}}
	\end{figure}

	Considering that the initial $\alpha$ of the DCM  may affect the convergence rate and precision, the influence of different initial values of $\alpha$ on the convergence rate and precision is analyzed, as shown in \Cref{fig:Prantl_alpha_tau_inital}b. The convergence rate is slower with increasing $\alpha$. We explain it from the perspective of the strong form, this is because of the larger constant of the Poisson equation. The neural network is relatively more difficult to fit since the Poisson equation is more apparent in ups and downs.

	The initial value of $\alpha$  will not only affect the convergence accuracy but also significantly influence the accuracy of $\alpha$ and have little effect on $\tau_{max}$. Overall, the accuracy of the four algorithms is close, while the accuracy of the DCM  is slightly higher. However, the DCM  needs to assume an admissible stress function. Since it is a full-force boundary condition problem, the DEM does not need to assume an admissible displacement field, so DEM is more convenient to deal with the full-force boundary condition problem, and the DEM based on the minimum potential energy principle can use the enormous approximate function space of the neural network to optimize and achieve good results.

	We compare the DCM-O (operator learning based on DCM) with DCM in this full-force boundary example. The input of the branch net is geometry information, $a$ and $b$. Because the problem hasn't the body field function, we deal with the $a$ and $b$ as the constant field function. The basis function is

	\begin{equation}
		\phi=(x^{2}-\frac{a^{2}}{4})(y^{2}-\frac{b^{2}}{4}),
	\end{equation}
	and the particular function is zero. The output of the DCM-O can be written as:
	\begin{equation}
		\phi=(x^{2}-\frac{a^{2}}{4})(y^{2}-\frac{b^{2}}{4})Branch(a,b;\boldsymbol{\theta})*Trunk(x,y;\boldsymbol{\theta}).
	\end{equation}
	The architecture of the trunk net consists of 6 hidden layers and 30 neurons in every hidden layer. The architecture of the branch net is 3 hidden layers and 30 neurons in every hidden layer. We use the 11 different $a/b$ ratios as the input of the branch net. The input of the trunk net is 100*100 equal spacing points. The output of DCM-O is the Prantl stress function $\phi$, and the output training dataset is from the analytical solution of the problem \citep{fung1953behavior} or the high fidelity numerical experiment (FEM based on the stress function element \citep{cen20118}). The analytical solution of the rectangle can be written as:
	\begin{equation}
		\begin{split}\phi(x,y) & =G\alpha(\frac{b^{2}}{4}-y^{2})-\frac{8G\alpha b^{2}}{\pi^{3}}\sum_{i=0}^{\infty}\frac{(-1)^{i}}{(2i+1)^{3}}\frac{cosh(\lambda_{i}x)}{cosh(\lambda_{i}a/2)}cos(\lambda_{i}y)\\
			\lambda_{i} & =(2i+1)\frac{\pi}{b}
		\end{split}
	\end{equation}
	Note that the test dataset (we only test $a=b=1$) differs from the training set. \Cref{fig: the DCM and DCM-P in Prantl stress function} shows the comparison of DCM with DCM-O, and the absolute error of DCM-O is lower than DCM in the same iteration. The absolute error of DCM is almost the same in the initial 10, 30, and 50 iterations, but DCM-O converges to the exact solution dramatically in those iterations. \Cref{fig: the DCM and DCM-P in Prantl stress function}m shows that the DCM-O converges faster than DCM apparently in the initial iteration, and DCM-O can reach the exact complementary energy in only 10 iterations with a small training dataset.

\begin{figure}

	\begin{centering}
		\includegraphics[scale=1.1]{DCM-DCM_O_prandtl}
		\par\end{centering}
	\caption{The absolute error of DCM and DCM-O in Prantl stress function of cylinder: (a, b, c, d), (e, f, g, h), and (i, j, k, l) the absolution error in the number of  iterations 10, 30, and 50 respectively; (a, e, i, c, g, k) the result of DCM ; (b, f, j, d, h, l) the result of DCM-O. (a, b, e, f, i, j) the result of $\tau_{zx}$; (c, d, g, h, k, l) the result of $\tau_{zy}$. (m) the evolution of the complementary energy by DCM and DCM-O; (n) the prediction of alpha by DCM and DCM-O \label{fig: the DCM and DCM-P in Prantl stress function}}

\end{figure}

\subsection{Airy stress function}

\subsubsection{Circular tube: full displacement boundary conditions} \label{sec: circular}
The circular tube is of interest to the question of the thick pressure vessel, which is the common container for gas and liquid in engineering \citet{vullo2014circular}.

On one end of the spectrum, all the boundary conditions are force boundary conditions, such as \Cref{sec: prantl}. If we use DEM, the admissible displacement is not required. On the other end of the spectrum, all the boundary conditions are displacement boundary conditions. If we use DCM, the admissible stress function is not required. The DEM can also be used to solve the full displacement boundary problem, but we have to require strong assumptions about admissible displacement fields (considering all boundary conditions).

In this section, we consider the full displacement boundary problem to compare the performance of DEM and DCM, demonstrating that different boundary conditions have different advantages and disadvantages with different energy principles. The DCM  does not need to assume the admissible stress function because there are no force boundary conditions, which is verified by the circular tube acting uniform internal and external pressure, as shown in the \Cref{fig:the schematic of circular tube} (the thick-walled cylinder). We use the displacement solution of the famous Lame formula to construct the displacement boundary condition. The analytical solution of the Lame stress and displacement is:

\begin{equation}
	\begin{split}
		\sigma_{r} & =\frac{a^{2}}{b^{2}-a^{2}}(1-\frac{b^{2}}{r^{2}})p_{i}-\frac{b^{2}}{b^{2}-a^{2}}(1-\frac{a^{2}}{r^{2}})p_{o}\\
		\sigma_{\theta} & =\frac{a^{2}}{b^{2}-a^{2}}(1+\frac{b^{2}}{r^{2}})p_{i}-\frac{b^{2}}{b^{2}-a^{2}}(1+\frac{a^{2}}{r^{2}})p_{o}\\
		u_{r} & =\frac{1}{E}[\frac{(1-\nu)(a^{2}p_{i}-b^{2}p_{o})}{b^{2}-a^{2}}r+\frac{(1+\nu)a^{2}b^{2}(p_{i}-p_{o})}{b^{2}-a^{2}}\frac{1}{r}]\\
		v_{\theta} & =0,
	\end{split}
\end{equation}
where $a$ and $b$ are inner and outer radius respectively, $E$ and $\nu$ are Young's  modulus and Poisson's ratio, $p_{i}$ and $p_{o}$ are  the uniform pressure at $r=a$ and $r=b$. The compressive stress is positive.

\begin{figure}
	\begin{centering}
		\includegraphics{circular_tube}
		\par\end{centering}
	\caption{Schematic diagram of the circular tube:  the value of parameters are $p_{o}=10$,  $p_{i}=5$, $a=0.5$,  $b=1.0$,  Young's modulus $E=1000$, and Poisson ratio $\nu=0.3$ \label{fig:the schematic of circular tube}}
\end{figure}

To use the DCM  to solve the problem, we need to express the complementary strain energy and complementary potential of the problem in the form of stress function. The stress expressed by the stress function of the axisymmetric problem in polar coordinates is:
\begin{equation}
	\begin{split}
		\sigma_{r}  =\frac{1}{r}\frac{\partial\phi}{\partial r}\\
		\sigma_{\theta} =\frac{\partial^{2}\phi}{\partial r^{2}}.
	\end{split}
\end{equation}

According to the constitutive law, the strain in polar coordinates is:
\begin{equation}
	\begin{aligned}\varepsilon_{r} & =\frac{1}{E}(\sigma_{r}-\nu\sigma_{\theta})=\frac{1}{E}(\frac{1}{r}\frac{\partial\phi}{\partial r}-\nu\frac{\partial^{2}\phi}{\partial r^{2}})\\
		\varepsilon_{\theta} & =\frac{1}{E}(\sigma_{\theta}-\nu\sigma_{r})=\frac{1}{E}(\frac{\partial^{2}\phi}{\partial r^{2}}-\nu\frac{1}{r}\frac{\partial\phi}{\partial r})
	\end{aligned}
\end{equation}

The complementary strain energy is expressed as a stress function:
\begin{equation}
	\begin{split}
		W_{c} & =\int_{\Omega}\frac{1}{2}\sigma_{ij}\varepsilon_{ij}d\Omega=\int_{\Omega}\frac{1}{2}(\sigma_{r}\varepsilon_{r}+\sigma_{\theta}\varepsilon_{\theta})2\pi rdr\\
		& =\int_{\Omega}\frac{1}{2}[\frac{1}{r}\frac{\partial\phi}{\partial r}\frac{1}{E}(\frac{1}{r}\frac{\partial\phi}{\partial r}-\nu\frac{\partial^{2}\phi}{\partial r^{2}})+\frac{\partial^{2}\phi}{\partial r^{2}}\frac{1}{E}(\frac{\partial^{2}\phi}{\partial r^{2}}-\nu\frac{1}{r}\frac{\partial\phi}{\partial r})]2\pi rdr
	\end{split}
\end{equation}

The complementary potential is:
\begin{equation}
	\begin{split}
		V_{c} & =u_{r}\sigma_{r}2\pi r|_{r=b}-u_{r}\sigma_{r}2\pi r|_{r=a}\\
		& =u_{r}\frac{\partial\phi}{\partial r}2\pi|_{r=b}-u_{r}\frac{\partial\phi}{\partial r}2\pi|_{r=a}
	\end{split}
\end{equation}
\begin{figure}
	\begin{centering}
		\includegraphics[scale=1.0]{DCM_prediction_tube}
		\par\end{centering}
	\centering{}\caption{The result of stress $\sigma_{r}$ and $\sigma_{\theta}$ of the circular tube by DCM: the exact solution of $\sigma_{r}$ (a) and $\sigma_{\theta}$ (d); the prediction of DCM in $\sigma_{r}$ (b) and $\sigma_{\theta}$ (e); the absolute error of DCM in $\sigma_{r}$ (c) and $\sigma_{\theta}$ (f). \label{fig:DCM result of circular tube}}
\end{figure}

In order to research the accuracy of  DCM, the analytical solution is brought into the complementary strain energy $W_{c}=79\pi/960$ and complementary potential $V_{c}=79 \pi/480$. \Cref{fig:DCM result of circular tube} shows the stress prediction of DCM compared with the analytical solution. The points are uniformly distributed. This problem is axisymmetric so the problem can be reduced to a one-dimensional problem. Thus, the number of the input and output neurons are both 1 (the input is $r$, the output is Airy stress function), 3 hidden layers, each layer 20 neurons, optimization function Adam, 2000 iterations. The absolute error of $\sigma_{r}$ is smaller than $\sigma_{\theta}$, the primary error of $\sigma_{r}$ is on the inner ring boundary, the primary error of $\sigma_{\theta}$ is on inner and outer ring boundary.
\Cref{fig:the radial directio of circular cube} shows the comparison between the stress and the analytical solution in the radial direction. The $\mathcal{L}_{2}$ errors  ($L_{2}=||pred-exact||_{2}/||exact||_{2}$) of $\sigma_{r}$ and $\sigma_{\theta}$ are 0.000266 and 0.00273, respectively. Therefore, the error of $\sigma_{r}$ is smaller than $\sigma_{\theta}$, mainly because $\sigma_{r}$ is the first derivative of the stress function, but $\sigma_{\theta}$ is the second derivative of the stress function. With increasing derivation order,  the error is larger in general. It is worth noting that since all boundaries are displacement boundaries, the admissible stress function does not need to be constructed in advance.

\begin{figure}
	\begin{centering}
		\includegraphics[scale=1.5]{DCM_tube_prediction_line}
		\par\end{centering}
	\centering{}\caption{The stress $\sigma_{r}$ (a) and $\sigma_{\theta}$ (b) of the circular tube in the radial direction under prediction solution of  the DCM compared with  analytical solution\label{fig:the radial directio of circular cube}}
\end{figure}

Considering the strong form of the stress function:
\begin{equation}
	\nabla^{4}\phi=0. \label{eq:the biharmonic equation}
\end{equation}
Adding some basis functions that naturally satisfy the biharmonic equation as shown in \Cref{eq:the biharmonic equation} will improve the accuracy and convergence speed of DCM. Thus, DCM-P is proposed by adding some basis functions to the DCM model. Take the circular tube as an example, and the biharmonic equation for the axisymmetric problem is simplified as follows:
\begin{equation}
	\frac{d^{4}\phi}{dr^{4}}+\frac{2}{r}\frac{d^{3}\phi}{dr^{3}}-\frac{1}{r^{2}}\frac{d^{2}\phi}{dr^{2}}+\frac{1}{r^{3}}\frac{d\phi}{dr}=0.
\end{equation}
The biharmonic function satisfying this equation is ${\rm ln}(r)$, $r^{2}$ and $r^{2}{\rm ln}(r)$. The stress of these biharmonic function are shown in \Cref{tab:biharmonic stress}, and we add the biharmonic function to the stress function field:
\begin{equation}
	\phi=NN(r;\theta)+a_{1}{\rm ln}(r)+a_{2}r^{2}+a_{3}r^{2}{\rm ln}(r),
\end{equation}
\begin{table}
	\caption{The stress fields (polar coordinates) are obtained by the Airy stress function satisfying the biharmonic function.\label{tab:biharmonic stress}}

	\centering{
	\begin{tabular}{cccc}
		\toprule
		Biharmonic function & ${\rm ln}(r)$ & $r^{2}$ & $r^{2}{\rm ln}(r)$\tabularnewline
		\midrule
		$\sigma_{r}$ & $1/r^{2}$ & 2 & $1+2{\rm ln}(r)$\tabularnewline
		$\sigma_{\theta}$ & $-1/r^{2}$ & 2 & $3+2{\rm ln}(r)$\tabularnewline
		$\tau_{r\theta}$ & 0 & 0 & 0\tabularnewline
		\bottomrule
	\end{tabular}
\end{table}
where $a_{1}$, $a_{2}$ and $a_{3}$ are set as optimization variables for optimization. \Cref{fig:biharmonic result of circular} compares the stress $L_{2}$ error with the biharmonic function, i.e., adding ${\rm ln}(r)$, $r^{2}$, $r^{2}{\rm ln }(r)$ separately and all biharmonic functions. After adding the biharmonic function, the stress accuracy has been improved to a certain extent, and the convergence speed has been accelerated at the same time. The improved effect of the biharmonic function ${\rm ln}(r) $ is the most obvious because the analytical solution contains the $1/r^{2}$ item. Thus it is equivalent to naturally constructing the stress function item that satisfies the analytical solution, and the accuracy and convergence speed is naturally fast. In addition, all biharmonic are added, and the accuracy and convergence speed is almost the same as just adding ${\rm ln}(r)$, which shows that DCM-P can well adjust the best biharmonic terms adaptively, making it prominent ${\rm ln}(r )$. After adding the biharmonic function, there will be a sudden change in a certain iteration step, which is caused by a dynamic game between the neural network approximation and the trainable coefficient of the biharmonic function. The change of the coefficient of the biharmonic function will bring about a larger change in accuracy than neural network parameters, and it takes a certain amount of iterations to use the strong approximation of the neural network to smooth out the sudden change in accuracy caused by the change of coefficient of the biharmonic function.

\begin{figure}
	\begin{centering}
		\includegraphics[scale=1.5]{DCMP_prediction_tube}
		\par\end{centering}

	\centering{}\caption{The overall relative $\mathcal{L}_{2}$ errors of stress $\sigma_{r}$ (a) and $\sigma_{\theta}$ (b) in different biharmonic functions of the circular tube under the DCM-P.\label{fig:biharmonic result of circular}}
\end{figure}

Next, we compare the accuracy and efficiency of the DEM, the PINNs displacement strong form ({\rm LN} equation), and the DCM  (without adding the biharmonic function) under full displacement boundary conditions, the reason why the strong form of the PINNs stress function is not compared is that the strong form of the stress function is inconvenient to deal with the displacement boundary conditions because it needs to be given in the form of boundary integrals. Thus, the strong form of the stress function is not considered. The construction method of the admissible displacement field in DEM and the strong form of PINNs displacement is:

\begin{equation}
	u=(1-\frac{r}{a})\frac{a}{a-b}u_{r}|_{r=b}+(1-\frac{r}{b})\frac{b}{b-a}u_{r}|_{r=a}+(r-a)(r-b)*NN(r;\boldsymbol{\theta}).\label{eq:=005747=00538B=005706=007B52=0053EF=0080FD=004F4D=0079FB=0089E3=005F62=005F0F}
\end{equation}

It is not difficult to verify that the above admissible displacement field satisfies the given displacement boundary condition. The configurations of these three methods are the same (10,000 points in the radial direction, optimizer Adam, 3 hidden layers of neural network, 20 neurons in each layer). When the number of iterations is 2000, the calculation time of the DEM  is 15.47s, the PINNs displacement strong form is 29.36s, and the DCM  is 34.03s. Since the strong form of PINNs displacement involves the second-order derivative, the computation time is larger than that of the DEM with only the first-order derivative. In addition, the DCM  also involves the second-order derivative, and at the same time, it must be accurately integrated like DEM. As a result, the calculation is required for the second-order derivative and precise numerical integration. Therefore, the calculation amount of the DCM is the largest compared with the other two forms.
In terms of stress fitting, \Cref{fig:stress_different method_line_circular tube}  and \Cref{fig:stress_different method_line_circular tube_error} show that DCM is the closest to the analytical solution, and the precision is the highest. The precision of the DCM is the highest because there is no assumption on the form of the stress function, and DCM completely releases the strong fitness ability of the neural network. However, the DEM and the PINNs displacement strong form are optimized under a certain displacement assumption space, and the optimization space is not as large as the DCM. Thus, the approximation error of DEM is greater than the DCM. In sum, for the problem of dealing with full-displacement boundary conditions, the DCM is theoretically more suitable than the DEM and the PINNs displacement strong form. Furthermore, it is better for unnecessary assumptions as few as possible.

\begin{figure}
	\begin{centering}
		\includegraphics[scale=1.5]{DCM_prediction_tube_all_methods_line}
		\par\end{centering}

	\centering{}\caption{Comparison of the predicted stress $\sigma_{r}$ (a) and $\sigma_{\theta}$ (b) with the analytical solution of the circular tube under the strong form of PINNs displacement, DEM,  and DCM  in the radial direction  \label{fig:stress_different method_line_circular tube}}
\end{figure}

\begin{figure}
	\begin{centering}

		\includegraphics[scale=1.5]{DCM_prediction_tube_all_methods}
		\par\end{centering}

	\centering{}\caption{Relative error $\mathcal{L}_{2}$ of the stress $\sigma_{r}$ (a) and $\sigma_{\theta}$ (b) of the circular tube under the strong form of PINNs displacement, DEM and DCM  \label{fig:stress_different method_line_circular tube_error}}
\end{figure}

We compare the DCM-O (operator learning based on DCM) with DCM in this full-displacement boundary example. The input of the branch net is pressure, $p_{i}$ and $p_{o}$. Because the problem does not have the force boundary condition, the basis and particular solution are not required. The output of the DCM-O can be written as:
\begin{equation}
	\phi=Branch(p_{i}, p_{o};\boldsymbol{\theta})*Trunk(r;\boldsymbol{\theta}).
\end{equation}
The architecture of the trunk net is 6 hidden layers and 30 neurons in every hidden layer. The architecture of the branch net is 3 hidden layers and 30 neurons in every hidden layer. We use the 121 different $p_{i}$ and $p_{o}$ as the input of the branch net. The input of trunk net is 11 equal spacing points from $p_{i} = 1.0$ to $10.0$ and from $p_{o} = 6.0$ to $15.0$, respectively. The output of DCM-O is the Airy stress function $\phi$, and the output training dataset is from the analytical solution of the problem or the high fidelity numerical experiment (FEM based on the stress function element \citep{cen20118}). The analytical solution of the circular tube can be written as:
\begin{equation}
	\phi=\frac{a^{2}}{b^{2}-a^{2}}[\frac{r^{2}}{2}-b^{2}ln(r)]p_{i}-\frac{b^{2}}{b^{2}-a^{2}}[\frac{r^{2}}{2}-a^{2}ln(r)]p_{o}.
\end{equation}
Note that the test dataset (we only test $p_{i}=5$ and $p_{o}=10$) is different from the training set. \Cref{fig: the DCM and DCM-P in circular} shows the comparison of DCM with DCM-O, and the absolute error of DCM-O is lower than DCM in the same iteration. \Cref{fig: the DCM and DCM-P in circular}m shows the DCM-O can converge faster than the DCM in the initial iteration step, and both can converge to the exact solution at the 2000 iteration.

\begin{figure}

	\begin{centering}
		\includegraphics[scale=1.1]{DCM_DCM_O_circular}
		\par\end{centering}
	\caption{Comparison of DCM and DCM-O in Airy stress function of circular tube: (a, b, c, d), (e, f, g, h), and (i, j, k, l) the absolution error in the number of iterations 10, 30, and 50 respectively; (a, e, i, c, g, k) $\sigma_{rr}$; (b, f, j, d, h, l) $\sigma_{\theta \theta}$; (a, b, e, f, i, j) the result of DCM; (c, d, g, h, k, l) the result of DCM-O; (m) the prediction of  $\sigma_{rr}$ by DCM and DCM-O of initial 50 and converged 2000 iterations; (n) the prediction of  $\sigma_{\theta \theta}$ by DCM and DCM-O of initial 50 and converged 2000 iterations \label{fig: the DCM and DCM-P in circular}}

\end{figure}

\subsubsection{Wedge: mixed boundary conditions} \label{sec: Wedge}
The wedge problem is another a common problem in engineering, such as the curling up of the sharp leading edge of a supersonic wing \citep{fung1953behavior}.
The problem is a mixed boundary condition (both displacement and force boundary conditions exist), and the normal pressure distribution is $p=qr^{m}$ acts on the wedge, as shown in \Cref{fig:wedge}a.
The Airy stress function can be analytically expressed as:
\begin{equation}
	\phi=r^{m+2}\{a*cos[(m+2)\theta]+b*sin[(m+2)\theta]+c*cos(m\theta)+d*sin(m\theta)\},
\end{equation}
where the shape and load of the boundary determine the constants a, b, c and d. We consider the specific situation $m=0$, which means the constant pressure on the upper boundary. The Airy stress function can analytically be written as:
\begin{equation}
	\phi=cr^{2}[\alpha-\theta+{\rm sin}(\theta){\rm cos}(\theta)-{\rm cos}^{2}(\theta){\rm tan}(\alpha)], \label{eq: Airy analytical stress function in wedge}
\end{equation}
where $c=q/[2({\rm tan}(\alpha)-\alpha)]$, $\alpha=\pi/6$, $q=5$, $L=\sqrt{3}$ , $E=1000$, $\mu=0.3$, as shown in \Cref{fig:wedge}b.

\begin{figure}
	\begin{centering}
		\includegraphics{wedge}
		\par\end{centering}
	\caption{The schematic of the problem of the wedge.\label{fig:wedge}}
\end{figure}

According to the stress formula in polar coordinates, the analytical solution can be obtained:
\begin{equation}
	\begin{aligned}\sigma_{r} & = \frac{1}{r^{2}} \frac{\partial^{2} \phi}{\partial\theta^{2}}+\frac{1}{r}\frac{\partial\phi}{\partial r}=2c[\alpha-\theta- \text{sin}^{2}(\theta){\rm tan}(\alpha)-{\rm sin}(\theta){\rm cos}(\theta)]\\
		\sigma_{\theta} & = \frac{\partial^{2}\phi}{\partial r^{2}} =2c[\alpha-\theta+{\rm sin}(\theta){\rm cos}(\theta)-{\rm cos}^{2}(\theta){\rm tan}(\alpha)]\\
		\tau_{r\theta} & =-\frac{\partial}{\partial r}(\frac{1}{r}\frac{\partial\phi}{\partial\theta}) =c[1-{\rm cos}(2\theta)-{\rm sin}(2\theta){\rm tan}(\alpha)].
	\end{aligned}
\end{equation}
We use the analytical solution to verify the accuracy of the DCM algorithm and compare it with DEM.

The key to DCM lies in the construction of the admissible stress function. Different from the admissible displacement field, the admissible stress function not only needs to satisfy the boundary value condition of the stress function at the force boundary but also the basis function network needs to satisfy the derivative condition. According to the boundary condition, the form of the stress function is:
\begin{equation}
	\phi=r^{2}f(\theta) \label{eq:Airy stress function in Widge assumption}
\end{equation}
We can use the neural network to approximate $f(\theta)$. Since the stress function must satisfy the force boundary condition in advance,  $f(\theta)$ satisfies the boundary condition according to the property of the Airy stress function:
\begin{equation}
	\begin{cases}
		f(0)=-\frac{1}{2}q\\
		f(\alpha)=0\\
		f^{'}(0)=0\\
		f'(\alpha)=0.
	\end{cases}\label{eq:=007F5A=0051FD=006570=0053EF=0080FD=005E94=00529B=0051FD=006570=007684=0065BD=0052A0}
\end{equation}

To analyze the convergence of DCM, we analyze the closeness of the loss function (complementary energy) to the analytical solution, and it is not difficult to obtain the accurate strain complementary energy $W_{c}=0.447$ and complementary potential $V_{c}=0$. However, different from the admissible displacement field, the admissible stress function  not only satisfy the admissible stress function value of the boundary but also satisfy the admissible stress function value in the derivative form, which will cause certain difficulties for the construction of the admissible stress function.
We can choose one of the methods from \Cref{tab:The-different-methods}.
Since there is no relevant research on the deep energy method based on the minimum complementary energy principle, for simplicity, we apply the penalty method  to satisfy the natural boundary conditions in this example:
\begin{equation}
	\mathcal{L}=\int_{\Omega}w_{c}d\Omega+\beta_{1}|\hat{f}(0)-f(0)|^{2}+\beta_{2}|\hat{f}(\alpha)-f(\alpha)|^{2}+\beta_{1}|\hat{f}'(0)-f'(0)|^{2}+\beta_{1}|\hat{f}'(\alpha)-f'(\alpha)|^{2},
\end{equation}
where $\hat{f}$ and $f$ are the predicted value and exact boundary value of the neural network, respectively, and $w_{c}$ is the minimum complementary energy density. The neural network input is a single variable $\theta$, 3 hidden layers, 20 neurons in each layer, and the output is a scalar Airy stress function; the optimizer Adam; the number of the uniformly spaced points is 100 points in $[0,\alpha]$. Here we have empirically selected $\beta_{1}=\beta_{2}=\beta_{3}=\beta_{4}=30$.

To compare the accuracy of DEM and DCM on mixed boundary conditions,  we use DEM to calculate the problem. Since the stress function of this problem can not accurately satisfy the fixed constraint of the essential boundary condition, it satisfies the integral sense. Therefore, under the framework of the principle of minimum potential energy, we need to treat the displacement boundary conditions as force boundary conditions through the total inverse method according to Saint Venant's Principle so that the wedge-shaped body subject to uniform loads has no displacement boundary conditions, and all boundary conditions are the force boundary condition. Thus there is no need to construct the admissible displacement field.
Considering that it is more convenient to express the displacement boundary condition in the Cartesian coordinate system, the stress field on the displacement boundary condition (right end) is expressed in the Cartesian coordinate system by using the coordinate conversion formula of tensor analysis:
\begin{equation}
	\begin{split}
		\boldsymbol{\sigma}_{xy} & =\left[\begin{array}{cc}
			{\rm cos}(\theta) & -{\rm sin}(\theta)\\
			{\rm sin}(\theta) & {\rm cos}(\theta)
		\end{array}\right]\left[\begin{array}{cc}
			\sigma_{r} & \tau_{r\theta}\\
			\tau_{r\theta} & \sigma_{\theta}
		\end{array}\right]\left[\begin{array}{cc}
			{\rm cos}(\theta) & -{\rm sin}(\theta)\\
			{\rm sin}(\theta) & {\rm cos}(\theta)
		\end{array}\right]^{T}\\
		& =\left[\begin{array}{cc}
			\sigma_{r}{\rm cos}^{2}(\theta)+\sigma_{\theta}{\rm sin}^{2}(\theta)-\tau_{r\theta}{\rm sin}(2\theta) & (\sigma_{r}-\sigma_{\theta}){\rm sin}(\theta){\rm cos}(\theta)+\tau_{r\theta}{\rm cos}(2\theta)\\
			(\sigma_{r}-\sigma_{\theta}){\rm sin}(\theta){\rm cos}(\theta)+\tau_{r\theta}{\rm cos}(2\theta) & \sigma_{r}{\rm sin}^{2}(\theta)+\sigma_{\theta}{\rm cos}^{2}(\theta)+\tau_{r\theta}{\rm sin}(2\theta)
		\end{array}\right]
	\end{split}
\end{equation}

The surface force on the displacement boundary at x=L is:

\begin{equation}
	\boldsymbol{\bar{t}}_{xy}=\boldsymbol{n}\cdot\boldsymbol{\sigma}_{xy}=\left[\begin{array}{c}
		\sigma_{r}{\rm cos}^{2}(\theta)+\sigma_{\theta}{\rm sin}^{2}(\theta)-\tau_{r\theta}{\rm sin}(2\theta)\\
		\sigma_{r}{\rm sin}(\theta){\rm cos}(\theta)-\sigma_{\theta}{\rm sin}(\theta){\rm cos}(\theta)+\tau_{r\theta}{\rm cos}(2\theta)
	\end{array}\right]
\end{equation}

On the other hand, the displacement in the xy direction at x=L is:
\begin{equation}
	\boldsymbol{u}_{xy}=\left[\begin{array}{cc}
		{\rm cos}(\theta) & -{\rm sin}(\theta)\\
		{\rm sin}(\theta) & {\rm cos}(\theta)
	\end{array}\right]\left[\begin{array}{c}
		u_{r}\\
		u_{\theta}
	\end{array}\right]=\left[\begin{array}{c}
		u_{r}{\rm cos}(\theta)-u_{\theta}{\rm sin}(\theta)\\
		u_{r}{\rm sin}(\theta)+u_{\theta}{\rm cos}(\theta)
	\end{array}\right]
\end{equation}

So the loss function in DEM can be expressed as:
\begin{equation}
	\mathcal{L}=W_{int}-W_{ex},
\end{equation}
where
\begin{equation}
	\begin{cases}
		W_{int}=\int_{\Omega}\frac{1}{2}(\sigma_{r}\varepsilon_{r}+\sigma_{\theta}\varepsilon_{\theta}+\tau_{r\theta}\gamma_{r\theta})d\Omega\\
		W_{ex}=\int_{\Gamma^{t}}\boldsymbol{\bar{t}}_{xy}\cdot\boldsymbol{u}_{xy}.
	\end{cases}
\end{equation}

Since the stress field is only related to $\theta$, according to the linear constitutive relation of polar coordinates, it is not difficult to obtain that the strain field is only related to $\theta$. Then according to the geometric equation, the displacement field limiting the displacement of the rigid body can be obtained, and the analytical solution of displacement is:
\begin{equation}
	\begin{cases}
		u_{r}=\frac{2c}{E}r[(1-\mu)(\alpha-\theta)-(1+\mu){\rm sin}(\theta){\rm cos}(\theta)-tan(\alpha)({\rm sin}^{2}(\theta)-\mu {\rm cos}^{2}(\theta))]\\
		u_{\theta}=-\frac{c}{E}r(1+\mu)[{\rm cos}(2\theta)+{\rm tan}(\alpha){\rm sin}(2\theta)]+\frac{4c}{E}r{\rm ln}(r).
	\end{cases}
\end{equation}
Considering the form of the displacement field, we can use a neural network to construct the displacement field as
\begin{equation}
	\begin{cases}
		u_{r}=rg_{1}(\theta)\\
		u_{\theta}=rg_{2}(\theta)+\frac{4c}{E}r{\rm ln}(r).
	\end{cases}
\end{equation}
Among them, $g_{1}(\theta)$ and $g_{2}(\theta)$ are both approximate functions of the neural network, both of which are input as  $\theta$, with 3 hidden layers, with 20 neurons in each layer. The output displacement fields are respectively $u_{r}$ and $u_{\theta}$; the optimizer Adam; the method of points distribution is uniformly spaced 100 points in $[0,\alpha]$.

To better quantify the error of the DEM, we analyze the variation trend of its $\mathcal{L}_{2}$ error with the number of iterations. \Cref{fig:DEM and DCM in wedge L2 stress}b indicates that the relative error of  $\sigma_{r}$, $\sigma_{\theta}$, and $\tau_{r\theta}$ are 0.915\%, 28.6\% and 2.86\% in DEM. Compared with DCM, DEM performs better in $\sigma_{r}$, but DCM performs better in $\sigma_{\theta}$. It is because DCM uses the Airy stress function to find $\sigma_{\theta}$ directly. However, when DEM calculates $\sigma_{\theta}$ through displacement, it calculates the derivative of the neural network (auto differential), which leads to the low accuracy and performance of DEM in $\sigma_{\theta}$; similarly, DEM does not involve the derivative of the neural network to get $\sigma_{r}$, but DCM calculates the second derivative of the neural network. Thus, DEM in $\sigma_ {r}$ performs better than DCM. Therefore, the choice of DCM and DEM not only needs to consider whether it involves the construction of admissible fields (when the displacement boundary conditions are complex, it is suitable for DCM, and when the force boundary conditions are complex, it is suitable for DEM), but also consider the order of the derivative. The above analysis provides an idea for selecting PINNs based on the energy method to solve the problem. It is worth noting that DEM has many assumptions in this example, but DCM does not. Therefore, in practical applications, DCM is more suitable for this example and only needs to make assumptions about the stress function.

\begin{figure}
	\begin{centering}
		\includegraphics[scale=1.5]{DCM_DEM_stress_wedge}
		\par\end{centering}

	\centering{}\caption{The overall relative error $\mathcal{L}_{2}$  of the wedge in $\sigma_{r}$,$\sigma_{\theta}$, and $\tau_{r\theta}$ under DCM (a) and DEM (b). \label{fig:DEM and DCM in wedge L2 stress}}

\end{figure}

For further comparison, we compare the PINNs strong form of the stress function with DCM and DEM due to the difficulty of the displacement assumption ($4cr{\rm ln}(r)/E$ are forcibly added according to the displacement analytical solution in $u_{\theta}$), so here we no longer consider the PINNs strong form of the displacement. We still assume the form of the Airy stress function in \Cref{eq:Airy stress function in Widge assumption}.
According to the biharmonic equation in strong form polar coordinates, we can get:
\begin{equation}
	\begin{cases}
		Domain: & f^{(4)}(\theta)+4f^{(2)}(\theta)=0\\
		Boundary: & f(0)=-\frac{q}{2}\\
		& f'(0)=0\\
		& f(\alpha)=0\\
		& f'(\alpha)=0,
	\end{cases}
\end{equation}
where the superscript (2) and (4) mean the second and fourth-order derivatives. The point allocation method, optimization scheme, and neural network configuration are completely consistent with DCM, and the boundary weight function is consistent with DCM ($\beta_{1}=\beta_{2}=\beta_{3}=\beta_{4}=30$). \Cref{fig:Comparison of different stress of different models in L2 error} and \Cref{fig:Comparison of different stress of different models with analytical solution} shows the stress result of the strong form of the PINNs stress function, the DCM, and DEM. The performance of DEM on $\sigma_{r}$ is the best because DEM  does not need derivatives to obtain $\sigma_{r}$; the performance of DCM in $\sigma_{\theta}$ is the best, also because the prediction of DCM in $\sigma_{\theta}$ does not require derivatives to obtain $\sigma_{theta}$; In the $\tau_{r \theta}$, the prediction accuracy of the three methods is almost the same.

\begin{figure}
	\begin{centering}
		\includegraphics[scale=1.0]{DCM_all_methods_line_wedge}
		\par\end{centering}

	\centering{}\caption{Comparison of the stresses $\sigma_{r}$ (a), $\sigma_{\theta}$ (b) and $\sigma_{r\theta}$ (c) of the DCM, DEM and PINNs stress functions of a wedge-shaped body subjected to a uniform load and the analytical solution \label{fig:Comparison of different stress of different models with analytical solution}}
\end{figure}

\begin{figure}
	\begin{centering}
		\includegraphics[scale=1.0]{DCM_all_methods_wedge}
		\par\end{centering}
	\centering{}\caption{Comparison of the overall relative error  $\mathcal{L}_{2}$ in stresses $\sigma_{r}$ (a), $\sigma_{\theta}$ (b) and $\tau_{r\theta}$ (c) of the DCM, DEM and PINNs stress functions of a wedge-shaped body subjected to a uniform load \label{fig:Comparison of different stress of different models in L2 error}}
\end{figure}

\begin{table}
	\caption{Accuracy and efficiency comparison of DCM, DEM and strong form of stress function for wedge with uniform load\label{tab:Accuracy and efficiency comparison of DCM}}
	\centering{
	\begin{tabular}{ccccc}
		\toprule
		Model & Relative error $\sigma_{r}$ & Relative error$\sigma_{\theta}$ & Relative error$\tau_{r\theta}$ & Time(S)\tabularnewline
		\midrule
		DCM & 11.7\% & \textbf{1.99\%} & 4.12\% & \textbf{494s}\tabularnewline
		DEM & \textbf{0.915\%} & 28.6\% & \textbf{2.86\%} & 571s\tabularnewline
		PINNs strong form: Airy stress & 5.32\% & 2.93\% & 5.11\% & 940s\tabularnewline
		\bottomrule
	\end{tabular}
\end{table}

\Cref{tab:Accuracy and efficiency comparison of DCM} indicates that the calculation efficiency of DCM is the highest because of the low order of derivation of DCM. The reason why DCM is more efficient than DEM is that DCM has only one neural network to fit the stress function, but DEM requires two neural networks to fit $u_{ r}$ and $u_{\theta}$, so DEM has more optimization parameters than DCM, so the calculation efficiency of DCM is better than DEM in this example.

We compare the DCM-O (operator learning based on DCM) with DCM. The input of the branch net is geometry information $\alpha$ and external force $q$. We satisfy the boundary condition by the penalty method, so the basis and particular functions are not required. The output of the DCM-O can be written as:
\begin{equation}
	\phi = r^{2} * Branch(\alpha, q;\boldsymbol{\theta})*Trunk(\theta ;\boldsymbol{\theta}).
\end{equation}
The architecture of the trunk net is 6 hidden layers and 30 neurons in every hidden layer. The architecture of the branch net is 3 hidden layers and 30 neurons in every hidden layer. We use the 110 different $\alpha$ and external force $q$  as the input of the branch net. The input of the trunk net is 100 equal spacing points. The output of DCM-O is the Airy stress function $\phi$, and the output training dataset is from the analytical solution of the problem in \Cref{eq: Airy analytical stress function in wedge}.
Note that the test dataset (we only test $\alpha = \pi/6$ and $q=5$) differs from the training set. \Cref{fig: Comparison of DCM and DCM-O in Airy stress function of wedge} shows the comparison of DCM with DCM-O. DCM-O demonstrates better accuracy when the algorithm converges, and DCM-O converges faster than DCM, comparing the stress result in the same iteration.

\begin{figure}

	\begin{centering}
		\includegraphics[scale=1.1]{DCM_DCM_O_wedge}
		\par\end{centering}
	\caption{Comparison of DCM and DCM-O in Airy stress function of the wedge: (a, b, c), (d, e, f), and (g, h, i) the prediction of $\sigma_{r}$,  $\sigma_{\theta}$, and  $\tau_{r \theta}$ in the number of iterations 100, 500, and 1000 respectively;  The evolution of the overall relative error $\mathcal{L}_{2}$ by DCM and DCM-O:  $\sigma_{r}$ (a, d, g, j);  $\sigma_{\theta}$ (b, e, h, k);  $\tau_{r \theta}$ (c, f, i, l). \label{fig: Comparison of DCM and DCM-O in Airy stress function of wedge}}

\end{figure}

\section{Discussion\label{sec:Discussion}}
\subsection{The construction of the basis function} \label{sec:basis function}

In this section, we discuss the construction of the basis function. The distance function is a special form of the basis function, as explained in \Cref{sec: DCM_method}. If the deep energy method is based on the principle of the minimum potential energy, the basis function is zero only if the position is on the Dirichlet boundary condition. Thus, we can construct different basis functions, as shown in \Cref{fig: distance function}. The idea of the contour line is adopted in the construction of the basis function.

\begin{figure}

	\begin{centering}
		\includegraphics[scale=1.1]{parrel}
		\par\end{centering}
	\caption{The schematic of the different basis function. \label{fig: distance function}}
\end{figure}

For the sake of simplicity, we use three different basis function $NN_{b}^{1}$, $NN_{b}^{2}$, and $NN_{b}^{3}$ to illustrate the way.
We assume that  the direction of the derivative of the $NN_{b}^{1}$, $NN_{b}^{2}$, and $NN_{b}^{3}$ are zero, i.e.

\begin{equation}
	\begin{aligned}\frac{\partial NN_{b}^{(1)}}{\partial \boldsymbol{l}_{1}} & =0\\
		\frac{\partial NN_{b}^{(2)}}{\partial \boldsymbol{l}_{2}} & =0\\
		\frac{\partial NN_{b}^{(3)}}{\partial \boldsymbol{l}_{3}} & =0,
	\end{aligned}
	\label{eq:=0065B9=005411=005BFC=0065700}
\end{equation}
where $\boldsymbol{l}_{1}$, $\boldsymbol{l}_{2}$, and $\boldsymbol{l}_{3}$ are the corresponding direction of the derivative. Suppose the partial differential operator is near the boundary. In that case, it is hard for the general net to learn because the $NN_{b}=0$ and $\partial NN_{b}/ \boldsymbol{l}=0$ when the direction is parallel to the boundary. The direction of the derivative can be written as follows:
\begin{equation}
	\begin{aligned}\frac{\partial\phi}{\partial \boldsymbol{l}_{1}} & =\frac{\partial NN_{p}}{\partial \boldsymbol{l}_{1}}+\sum_{i=1}^{3}\frac{\partial NN_{g}^{(i)}}{\partial \boldsymbol{l}_{1}}*NN_{b}^{(i)}+\sum_{i=1}^{3}NN_{g}^{(i)}*\frac{\partial NN_{b}^{(i)}}{\partial \boldsymbol{l}_{1}}\\
		\frac{\partial\phi}{\partial \boldsymbol{l}_{2}} & =\frac{\partial NN_{p}}{\partial \boldsymbol{l}_{2}}+\sum_{i=1}^{3}\frac{\partial NN_{g}^{(i)}}{\partial \boldsymbol{l}_{2}}*NN_{b}^{(i)}+\sum_{i=1}^{3}NN_{g}^{(i)}*\frac{\partial NN_{b}^{(i)}}{\partial \boldsymbol{l}_{2}}\\
		\frac{\partial\phi}{\partial \boldsymbol{l}_{3}} & =\frac{\partial NN_{p}}{\partial \boldsymbol{l}_{3}}+\sum_{i=1}^{3}\frac{\partial NN_{g}^{(i)}}{\partial \boldsymbol{l}_{3}}*NN_{b}^{(i)}+\sum_{i=1}^{3}NN_{g}^{(i)}*\frac{\partial NN_{b}^{(i)}}{\partial \boldsymbol{l}_{3}}.
	\end{aligned}
	\label{eq:the dirivative of the basis function}
\end{equation}
Because $\partial NN_{b}^{1}/ \boldsymbol{l}_{1}=0$,  $\partial NN_{b}^{2}/ \boldsymbol{l}_{2}=0$, and  $\partial NN_{b}^{3}/ \boldsymbol{l}_{3}=0$,
\Cref{eq:the dirivative of the basis function} can be further written as
\begin{equation}
	\begin{aligned}\frac{\partial\phi}{\partial l_{1}} & =\frac{\partial NN{p}}{\partial l_{1}}+NN{g}^{(2)}*\frac{\partial NN{b}^{(2)}}{\partial l_{1}}+NN{g}^{(3)}*\frac{\partial NN{b}^{(3)}}{\partial l_{1}}\\
		\frac{\partial\phi}{\partial l_{2}} & =\frac{\partial NN{p}}{\partial l_{2}}+NN{g}^{(1)}*\frac{\partial NN{b}^{(1)}}{\partial l_{2}}+NN{g}^{(3)}*\frac{\partial NN{b}^{(3)}}{\partial l_{2}}\\
		\frac{\partial\phi}{\partial l_{3}} & =\frac{\partial NN{p}}{\partial l_{3}}+NN{g}^{(1)}*\frac{\partial NN{b}^{(1)}}{\partial l_{3}}+NN{g}^{(2)}*\frac{\partial NN{b}^{(2)}}{\partial l_{3}}
	\end{aligned}
	\label{eq:=005316=007B80=00540E=007684=0053EF=0080FD=005E94=00529B=0051FD=006570=007684=0065B9=005411=005BFC=006570}
\end{equation}

In the energy principle training, the basis function network and the particular net are frozen and not trained, and only the general network can be trained. The formula \ref{eq:=005316=007B80=00540E=007684=0053EF=0080FD=005E94=00529B=0051FD=006570=007684=0065B9=005411=005BFC=006570} shows that although one of the basis function networks makes the corresponding general network impossible training, the other two generalized networks are still trainable. Since the contours of the other two basis function networks must have intersections (except parallel), they will be at the intersections and cannot be trained. If there are three basis function networks, the three contours may not have corresponding intersections (if accuracy requirements are not high, two basis function networks should also achieve good results).

\subsection{The comparison of DEM, PINNs strong form, and DCM}
PINNs strong form of stress function is suitable for dealing with full-force boundary conditions. However, it is not convenient to solve the problem with the displacement boundary because displacement boundary conditions are formulated in the form of boundary integral equations. On the other hand, DCM based on the minimum complementary energy principle only needs to convert the displacement boundary conditions into the complementary potential for optimization when dealing with the displacement boundary conditions, which is difficult for the PINN strong form of the stress function. From the construction of the admissible field, DEM is more suitable for dealing with the problems dominated by the force boundary conditions because it does not require the construction of the admissible displacement field on the force boundary condition. In the same way, DCM is suitable for dealing with displacement boundary conditions.

Therefore, from the perspective of admissible field construction, different problems are suitable for different deep energy methods (DEM or DCM). One needs to understand the problem's nature so we can choose the suitable deep energy method based on potential or complementary energy. Note that the DEM and DCM can solve the same problem, not to say that only the principle of minimum potential energy or the minimum complementary energy can be used.

\subsection{The input of the branch net in DCM-O}
If the input of DeepONet's branch net is a constant field, the constant value can be input to the branch net instead of the variable field function considering the computational efficiency.

Generally speaking, the input of the branch net is a variable field. We can choose the suitable architecture of the branch net to be consistent with the nature of the input of the trunk net. For example, if the input of the trunk net is the image structure, we can use the CNN in the trunk net; if the input of the trunk net is the data associated with the time, we can use the RNN to be the architect of the branch net. The key to choosing the architecture of the branch net is to consider the nature of the input of the branch net.

\section{Conclusion\label{sec:Conclusion}}
We propose the deep complementary energy method (DCM) based on the minimum complementary energy principle, which is suitable for solving problems dominated by displacement boundary conditions. We compare the DEM, DCM, and PINNs strong forms for verification on three types of boundary conditions (full displacement boundary condition, full force boundary condition, and mixed boundary condition problem). From the perspective of constructing admissible function fields, DEM is more suitable for dealing with full-force boundary conditions, and DCM is more suitable for dealing with full-displacement boundary conditions. In terms of computational efficiency, DCM and DEM are better than the strong form of PINNs due to low-order derivatives. Furthermore, the proposed DCM algorithm is extended to DCM-Plus (DCM-P), adding the terms that naturally satisfy the biharmonic equation in the Airy stress function. Deep complementary energy operator method (DCM-O) is proposed by combining operator learning with physical equations, including the branch net, trunk net, basis net, and general net. The result of the Prandtl and Airy stress function shows that DCM-P and DCM-O can improve the accuracy and efficiency of DCM.

At present, the DCM algorithm only analyzes two common stress functions (Airy and Prantl stress functions) and has not yet been tested on a three-dimensional problem, which involves Morera and Maxwell stress functions. In addition, currently, the force boundary condition is constructed by the penalty function method in the Airy stress function. Although the basis function theory was proposed and tested on the Prantl stress function, more research on the basis function method on the Airy function is needed in the future. We have proved that DCM can solve the same nonlinear energy problem theoretically, such as hyperelasticity, if DEM can solve the problem in solid mechanics, but we have yet to test it in the numerical experiment, which we will further study in the future. This work proposes an important supplementary energy form of the deep energy method. Complementary energy exists in solid mechanics and other physical fields, such as Helmholtz free energy, enthalpy, Gibbs free energy and internal energy in thermodynamics. The relationship between energy conforms to the Legendre transformation in mathematics. Therefore, DCM can also apply to the energy principles of Legendre transformation in other physics. The work in this paper reflects the broad prospect of combining deep learning with computational mechanics. It is not easy for a purely physics-based energy method to surpass traditional finite elements, and a purely data-driven model will bring about the problem of excessive data requirement. With the in-depth research on neural networks and the advancement of computing power, combined with more and more calculation data results, there are strong reasons to believe that operator learning based on the energy method can bring a balance between data and physical equation.  If more high-precision calculation results are stored in the future, the operator learning based on the energy method has the potential to exceed the traditional algorithms, which will give computational mechanics the new and broad research prospects.

\section*{Acknowledgement}
The study was supported by the Major Project of the National Natural Science Foundation of China (12090033). The authors would like to thank Jingyun Bai and Yuqing Du for helpful discussions.