# Comparison of Neural FEM and Neural Operator Methods for applications in
Solid Mechanics

Stefan Hildebrand  TU Berlin, Berlin, Germany<EMAIL_ADDRESS>Sandra Klinge  TU Berlin, Berlin, Germany

(1 April 2022; 1 May 2022; 10 May 2022)

###### Abstract

Machine Learning methods belong to the group of most up-to-date approaches for
solving partial differential equations. The current work investigates two
classes, Neural FEM and Neural Operator Methods, for the use in elastostatics
by means of numerical experiments. The Neural Operator methods require
expensive training but then allow for solving multiple boundary value problems
with the same Machine Learning model. Main differences between the two classes
are the computational effort and accuracy. Especially the accuracy requires
more research for practical applications.

###### keywords:

Neural FEM; Neural Operator Methods; Machine Learning; Partial Differential
Equation, Elastostatics

20 May 2022 Editor Name First Reviewer, Second Reviewer, Third Reviewer

##  1 Introduction

Induced by ever-rising compute power and successful applications in several
domains, Artificial Intelligence (AI) systems and especially Machine Learning
(ML) methods attract growing attention for advanced tasks in mechanical
engineering [1, 2, 3]. This is supported by well-established and flexible
Machine Learning frameworks like PyTorch [4] and Tensorflow [5].

One particular application of ML is the solution of Parameterized Partial
Differential Equations (PPDE), which are traditionally solved by numerical
discretization methods like Finite Element Method (FEM), Finite Difference
Method (FDM), Finite Volume Method (FVM) or Boundary Element Method (BEM).
Based on ML techniques, two new classes of methods arose, namely the Neural
FEM and Neural Operator methods [6]. The aim of the current work is to compare
these two classes of methods for applications in solid body mechanics.
Therefor, their common representatives are applied to case studies, where the
well-established FEM can serve as a benchmark.

The mathematical problem to solve with either method can be described as
follows:

Let an arbitrary Parameterized Partial Differential Equation be given on an
open domain $B$ with piecewise smooth boundary $\Gamma$ in the form:

$\displaystyle\mathcal{N}[\boldsymbol{u}(\boldsymbol{y});\boldsymbol{y}]=\mathbf{0}\quad\text{on}\,\,B,\quad\mathcal{B}[\boldsymbol{u}(\boldsymbol{y});\boldsymbol{y}]=\mathbf{0}\quad\text{on}\,\,\Gamma\,\,,$
(1)

where $\mathcal{N}$ is a nonlinear operator on the domain $B$, $\mathcal{B}$
an operator on $\Gamma$ that determines the boundary conditions, and
$\boldsymbol{u}(\boldsymbol{y})\in\mathbb{R}^{d}$ the solutions of the PDE.
All quantities are parameterized by $\boldsymbol{y}\in\mathbb{R}^{n}$. The
mapping

$\displaystyle G:\quad
B\cup\Gamma\,\,\times\,\,\mathbb{R}^{n}\rightarrow\mathbb{R}^{d},\quad\left(\boldsymbol{X},\boldsymbol{y}\right)\mapsto\boldsymbol{u},\quad\boldsymbol{X}\in
B\cup\Gamma,\quad n,d\in\mathbb{N}$ (2)

is called the solution operator of the PPDE.

Neural FEM resembles a conventional FEM implementation. The artificial Neural
Network (NN) approximates the solution function of a particular realization of
the PPDE. Based on the conventional form of Physics-Informed Neural Networks
(PINN, [7]), the Deep Energy Method (DEM) and competitive PINNs (cPINN) are
proposed in [8, 9] . All these approaches are independent from a spatial
discretization of the domain $B$ (grid-independent) and can realize high
accuracies, but must be retrained for each new set of parameters.

In the case of Neural Operator methods, an NN is trained to behave like the
solution operator of a PPDE. Then, the network can be applied to arbitrary
combinations of parameters and boundary conditions, to solve Boundary Value
Problems (BVP). These methods are particularly characterized by a
discretization-independent error, allowing zero-shot super resolution
(training on coarse grid, inference on fine grid). In this work, the Deep
Operator Network (DeepONet) [10] and the Fourier Neural Operator (FNO) [6] are
studied as representatives of Neural Operator methods. Typically, Neural
Operator Methods require a large amount of training data, which may need to be
computed in a numerically expensive way [11].

Physics-Informed variants of Neural Operator methods address this drawback by
incorporating knowledge on the underlying PDE as a regularizing mechanism in
the loss function [7]. This can increase accuracy, generalizability, and data
efficiency [12]. The present contribution investigates Physics-Informed
DeepONet (PIDeepONet) [12] and Physics Informed Neural Operator (PINO) [13] as
representatives of Physics-Informed Neural Operator methods.

The paper is structured as follows. First, we give an insight in details of
selected ML methods (Sections 3 and 4). Then, we apply these methods to three
example problems (Section 5). Finally, we discuss the results in comparison to
the reference solution from FEM and highlight necessary steps for a future use
of these NN methods in elastostatics (Section 6).

##  2 Neural Network Nomenclature

Artificial Neural Networks are constructed as layers of neurons [14, 1]. Each
neuron carries out a (typically nonlinear) activation function. In case of a
Fully Connected Neural Network (FCNN), each neuron receives its input as
linear transformation of all outputs of the neurons on the layer before. The
output $\boldsymbol{\mathcal{R}}^{i}$ of the $i^{th}$ layer is thus calculated
by:

$\displaystyle\boldsymbol{W}^{i}$
$\displaystyle=\boldsymbol{w}^{i}\boldsymbol{\mathcal{R}}^{i-1}+\boldsymbol{b}^{i}$
(3) $\displaystyle\boldsymbol{\mathcal{R}}^{i}$
$\displaystyle=\boldsymbol{a}^{i}(\boldsymbol{\Theta}^{i},\boldsymbol{W}^{i})$
(4)

The weights $\boldsymbol{w}^{i}$ and biases $\boldsymbol{b}^{i}$ of the
$i^{th}$ layer, together with the parameters $\boldsymbol{\Theta}^{i}$ of the
activation functions $\boldsymbol{a}$, form the set of parameters $\theta$ of
the Neural Network. Other parameters of the network architecture (like number
of layers and layer widths) and the optimizer algorithm (e.g. step width) are
called hyperparameters and have to be chosen by the NN user or an outer
optimization strategy.

All the layers between the input an output layer are called hidden layers. The
number of (hidden) layers is referred to as the network’s depth, whereas the
number of neurons within a layer is called the width of the layer.

The whole network represents an arbitrary (continuous) mapping (Universal
approximation theorem, [15]) $\mathcal{R}:\mathbb{R}^{n}\mapsto\mathbb{R}^{m}$
from the input to the output side, with $n$ and $m$ the input and output layer
width, respectively.

To make the network approximate a mapping $\tilde{\mathcal{R}}$ on a subset
$D\subset\mathbb{R}^{n}$, the mapping is given indirectly by a set $T_{tr}$ of
training tuples
$t_{tr}^{k}\left(P^{k},\tilde{\mathcal{R}}(P^{k})\right),\,\,P^{k}\in
D,\,\,k\in\mathbb{N}$, which together form the training data set. $P^{k}$ are
the input samples, $\tilde{\mathcal{R}}(P^{k})$ the corresponding target
outputs. A set $T_{te}$ of testing tuples $t_{te}^{l},\,\,l\in\mathbb{N}$ is
required to check the quality of the approximation the network has learned so
far. Usually, $T_{tr}\cap T_{te}=\emptyset$. In the application of a network,
arbitrary sets of data within $D$ can be the input, but the exact output is
usually unknown and only approximated by the net.

Conventionally, the parameters of the network are adapted by an optimization
algorithm like Adam [16]. This algorithm minimizes the empirical risk $F$
(alternatively called loss value) that is calculated as the output of a loss
function $\mathcal{L}$. The latter is commonly defined as the discrepancy
between the target outputs and the outputs of the network with the current
parameters. A frequent choice for the loss function is Mean Square Error (MSE)
for $N\in\mathbb{N}$ tuples $t^{k}$

$\displaystyle\mathrm{MSE}=\frac{1}{N}\sum_{k=1}^{N}\left(\tilde{\mathcal{R}}(P^{k})-\mathcal{R}(P^{k})\right)^{2}\quad.$
(5)

The optimization employs (partial) derivatives of the loss function w.r.t. the
parameters in the network. Machine learning networks like PyTorch therefor
record all operations acting on a variable from input to output symbolically,
so that fast, highly accurate derivation becomes possible. This feature is
referred to as autograd [17]. Its use, however, is not limited to derivations
w.r.t. network parameters.

Figure 1: Information flow in a Fully Connected Neural Network (FCNN).

By default, the NN parameters are initialized randomly before the first
optimizer step. This leads to varying results of effort and achieved accuracy
during training.

##  3 Neural FEM

In this class of methods, the output of the NN is chosen as the unknown
function of the PDE. The loss function is then computed either from the
residual of the PDE (classical Physics-informed Neural Networks (PINN) [7] and
competitive PINN (cPINN) [9]), or from the potential energy if the minimum
principle applies (Direct Energy Method (DEM) [8], mixed DEM (mDEM) [18]),
both incorporating the outputs of the NN. Hence, no training data set is
needed.

###  3.1 Physics Informed Neural Networks (PINN)

The original PINN formulation uses an FCNN where the loss function is applied
to the squared residual of the PDE on specified collocation points, [7],

$\displaystyle
F^{\mathrm{PINN},B}=\frac{1}{N_{f}}\sum_{i=1}^{N_{f}}\left(\mathcal{N}\left[\boldsymbol{u}_{\theta}\right]\left(\boldsymbol{X}_{i}^{f}\right)\right)^{2}\quad.$
(6)

The boundary conditions are accounted for by an additional term in the form

$\displaystyle
F^{\mathrm{PINN},\Gamma}=\frac{\lambda_{b}}{N_{b}}\sum_{i=1}^{N_{b}}\left(\mathcal{B}\left[\boldsymbol{u}_{\theta}\right]\left(\boldsymbol{X}_{i}^{b}\right)\right)^{2},$
(7)

where $\lambda_{b}$ is a hyperparameter that weighs the error proportions,
since the propagated gradients can be of different magnitudes, thus driving
the optimization procedure towards an incorrect solution [19]. The total
empirical risk is then calculated by

$\displaystyle
F^{\mathrm{PINN}}\left(\boldsymbol{u}_{\theta}\right)=F^{\mathrm{PINN},B}+F^{\mathrm{PINN},\Gamma}\quad.$
(8)

According to [20], the training of PINNs fails even on very simple problems,
such as the 1D convection or the reaction-diffusion equation. An analysis of
the occurrence of comparable pathologies in elastostatic or elastodynamic
contexts requires further research. In the survey at hand, they did not
manifest.

####  3.1.1 Deep Collocation Method (DCM)

The Deep Collocation Method (DCM) [14] is a representative of the classical
PINN, where the empirical risk is built up by the squared residual at random
collocation points. The constraints are typically accounted for by additional
penalty terms but not enforced by a transformation of the output data.

###  3.2 The Deep Energy Method (DEM)

The DEM was originally introduced to calculate finite deformation
hyperelasticity [8]. This method as well as methods derived from it require
only first derivatives to compute the loss function, thus reducing the
numerical complexity. In return, errors are generated by the numerical
integration of the energy function.

The solution is sought in form of the displacement field
$\boldsymbol{u}(\boldsymbol{X})$ that corresponds to the minimum total
potential energy $\Pi$. This minimization can be accomplished by choosing the
loss function $F$ to calculate the total potential energy $F:=\Pi$. The input
of the NN are points in the reference configuration $\boldsymbol{X}\in
B\cup\Gamma$, in the domain $B$ or on its boundary $\Gamma$.

A transformation is applied to integrate the geometric boundary conditions:
Let the output of the NN be given by $\boldsymbol{z}(\theta,\boldsymbol{X})$.
To retrieve the displacement field $\boldsymbol{u}_{\theta}(\boldsymbol{X})$
based on the parameter set $\theta$ of the NN, the displacements on the
boundary are introduced in a separate term
$\boldsymbol{u}_{g}(\boldsymbol{X})$. Additionally, a mapping
$\boldsymbol{A}(\boldsymbol{X})$ with
$\boldsymbol{A}(\boldsymbol{X})=\mathbf{0}\text{ for }\boldsymbol{X}\in\Gamma$
is introduced. Then, the output is constructed as:

$\displaystyle\boldsymbol{u}_{\theta}(\boldsymbol{X})=\boldsymbol{u}_{g}(\boldsymbol{X})+\boldsymbol{A}(\boldsymbol{X})\boldsymbol{z}(\theta,\boldsymbol{X})$
(9)

Now, $\boldsymbol{u}_{\theta}$ automatically fulfills the boundary constraints
and a nonlinear optimization problem without constraints is obtained. This
problem can be solved with an optimization procedure such as the Limited-
memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm.

The corresponding loss function to minimize is

$\displaystyle
F^{\mathrm{DEM}}\left(\boldsymbol{u}_{\theta}\right)=\int_{B}\left(W\left(\boldsymbol{F}_{\theta}\right)-\boldsymbol{b}\cdot\boldsymbol{u}_{\theta}\right)\,\mathrm{d}V-\int_{\Gamma_{N}}\overline{\boldsymbol{T}}\cdot\boldsymbol{u}_{\theta}\,\mathrm{d}A$
(10)

with the elastic strain energy $W$, the deformation gradient
$\boldsymbol{F}_{\theta}=\frac{\partial(\boldsymbol{X}+\boldsymbol{u}_{\theta})}{\partial\boldsymbol{X}}$,
body forces $\boldsymbol{b}$ and the traction
$\overline{\boldsymbol{T}}=\boldsymbol{P}\cdot\boldsymbol{N}$ defined as the
$1^{st}$ Piola-Kirchhoff stress tensor $\boldsymbol{P}$ projection on the
outward normal of $\Gamma$. The deformation gradient $\boldsymbol{F}_{\theta}$
contains only first derivatives and is retrieved by the autograd feature of
the NN framework. The described procedure is shown in Fig. 2.

Figure 2: Information flow in DEM.

Several numerical methods to calculate an approximation of the integrals in
the loss function have been suggested [8]. Some examples are the Monte Carlo
integration and the trapezoidal rule, that we use for the examples presented
in Section 5.

An alternative variant of the same concept is the Shallow Energy Method (SEM),
where the deep NN is replaced by a shallow NN with a single hidden layer,
activated by Radial Basis Functions (RBF) [21].

An other enhancement of the basic DEM is the mixed DEM (mDEM) [18], where both
displacements $\boldsymbol{u}_{\theta}$ and stresses $\boldsymbol{P}_{\theta}$
are calculated by the NN. The deviation from the constitutive law derived from
the strain energy $W$ must then be integrated into the loss function to
represent the correct material behavior

$\displaystyle
F^{\mathrm{mDEM}}=F^{\mathrm{DEM}}+\frac{V}{N_{f}}\sum_{j=1}^{N_{f}}\left|\left|\boldsymbol{P}_{\theta}\left(\boldsymbol{X}_{j}^{f}\right)-\left.\frac{\partial
W\left(\boldsymbol{F}_{\theta}\right)}{\partial\boldsymbol{F}}\right|_{\boldsymbol{X}_{j}^{f}}\right|\right|_{2}^{2}\quad.$
(11)

The prescribed forces on the Neumann boundary part $\Gamma_{N}$ can be
accounted for directly by a transformation similar to Eq. (9) holding for the
geometric boundary conditions on the Dirichlet boundary part $\Gamma_{D}$.
Alternatively, an additional error term can be included to penalize the
squared deviation from the prescribed forces.

###  3.3 competitive PINN (cPINN)

cPINNs extend the idea of PINNs by formulating the learning problem as a zero-
sum game in the style of a Generative Adversarial Network (GAN) [9] with a
Nash equilibrium that corresponds to the analytical solution of the PDE. This
avoids the use of the square of the residual, which aims to improve the
learning performance. Compared to a classical PINN, cPINN introduces an
additional discriminator FCNN with NN parameter set $\phi$ which is trained to
predict errors of the PINN. Let $\boldsymbol{u}_{\theta}$ be furthermore the
output of the PINN and additionally
$\boldsymbol{d}_{\phi}=\left(\boldsymbol{d}_{\phi}^{B},\boldsymbol{d}_{\phi}^{\Gamma}\right)$
the output of the discriminator network. Then the minimax formulation of the
game is given by

$\displaystyle\max_{\phi}\min_{\theta}F^{\mathrm{cPINN}}\left(\boldsymbol{u}_{\theta},\boldsymbol{d}_{\phi}\right)=\left.\frac{1}{N_{f}}\sum_{i=1}^{N_{f}}\mathcal{N}\left[\boldsymbol{u}_{\theta}\right]\cdot\boldsymbol{d}_{\phi}^{B}\right|_{\boldsymbol{X}_{i}^{f}}+\left.\frac{\lambda_{b}}{N_{b}}\sum_{i=1}^{N_{b}}\mathcal{B}\left[\boldsymbol{u}_{\theta}\right]\cdot\boldsymbol{d}_{\phi}^{\Gamma}\right|_{\boldsymbol{X}_{i}^{b}}\,\,\,\,.$
(12)

[9] [9] solve this optimization problem by using the Adam based Competitive
Gradient Descent (ACGD) method.

##  4 Neural Operator Methods

In a multi-query context, where a PDE must be evaluated for a large number of
parameters, classical methods are computationally intensive. This includes
both the conventional FEM and the Neural FEM methods explained previously.
That drawback has motivated a large body of work on model order reduction,
which deals with the tradeoff towards reduced accuracy, stability, and
generalizability. Learning solution operators between (infinite-dimensional)
function spaces using NNs is a comparatively young field.

In this class of methods, the NN approximates the solution operator of the
PPDE, i.e. for a given parameter set, the NN shall output the solution of the
PDE at the points of interest. Now, the objective is defined as a risk
functional that takes the probability distribution $\chi$ of the parameter set
$\boldsymbol{y}$ into account [6]

$\displaystyle
F=\int_{\chi}\mathcal{L}\left(G_{\theta}(\boldsymbol{y}),G(\boldsymbol{y})\right)\mathrm{d}\chi\,\,,$
(13)

where $G$ is the solution operator and $G_{\theta}$ its approximation by the
NN.

A possibility to represent operators by means of NNs is a finite-dimensional
approximation of the function spaces and interpolation of these spaces by the
NN. Let e.g. the Boundary Value Problem (BVP) be given with approximations to
the solutions calculated by traditional FEM on the node points. Now, the
solution operator can be sought that maps the volume forces to the
displacements on the node points, hence training a discrete operator [22, 23].
However, this approach introduces a grid-dependency since the results are only
obtained on the initially chosen node points.

Alternatively, the points of interest in space can be included as input
parameter into the mapping that shall be represented by the NN,

$\displaystyle\boldsymbol{u}_{\theta}:\quad
B\cup\Gamma\times\mathbb{R}^{\left(d+N_{\text{dof}}\right)}\rightarrow\mathbb{R}^{N_{\text{dof}}}\,\,,\quad[\boldsymbol{X},\mathbf{y}]\mapsto\boldsymbol{u}_{\theta}(\boldsymbol{X},\mathbf{y})\,\,,$
(14)

where $\mathbf{y}$ is taken as the parameter of the displacement field
$\boldsymbol{u}_{\theta}$.

In the following, Deep Operator Network (DeepONet) and Fourier Neural Operator
(FNO) are discussed. Both approaches are discretization independent and allow
for small generalization errors.

###  4.1 Deep Operator Network (DeepONet) and Physics Informed DeepONet
(PIDeepONet)

NN can be employed as universal approximators of continuous functions, as well
as for nonlinear continuous operators [12, 10]

$\displaystyle\quad|G(\boldsymbol{y})(\boldsymbol{X})-\underbrace{\boldsymbol{g}\left(\boldsymbol{y}\left(\boldsymbol{X}_{1}\right),\ldots,\boldsymbol{y}\left(\boldsymbol{X}_{m}\right)\right)}_{\text{branch}}\cdot\underbrace{\boldsymbol{f}(\boldsymbol{X})}_{\text{trunk}}|<\varepsilon,$
(15)

where $\boldsymbol{f}$ (trunk) and $\boldsymbol{g}$ (branch) can be
represented by various classes of neural networks that satisfy the
requirements of the classical universal approximation theorem [15]. It is
assumed that the parameter $\boldsymbol{y}$ is known on sufficiently many grid
points $m$. On this basis, the stacked and unstacked Deep Operator Network
(DeepONet) are proposed in [10]. The stacked DeepONet differs from the
unstacked one only in the definition of the branch networks. In the unstacked
DeepONet, these are combined into one net to facilitate training (Fig. 3).

Let
$\boldsymbol{y}=\left[\boldsymbol{y}\left(\boldsymbol{X}_{1}\right),\ldots,\boldsymbol{y}\left(\boldsymbol{X}_{m}\right)\right]^{\top}$
be the parameter function at discrete points in space. Moreover, let the
output of DeepONet be $G_{\theta}(\mathbf{y})(\boldsymbol{X})$. Then, the
empirical risk functional for DeepONet is given by

$\displaystyle
F^{\text{DeepONet}}=\frac{1}{P\,N}\sum_{i=1}^{N}\sum_{j=1}^{P}\left(G_{\theta}\left(\boldsymbol{y}_{i}\right)\left(\boldsymbol{X}_{j}\right)-G\left(\boldsymbol{y}_{i}\right)\left(\boldsymbol{X}_{j}\right)\right)^{2}\quad.$
(16)

Here, $N$ is the number of realizations of the parameter input
$\boldsymbol{y}$ available for training and $P$ is the number of training data
known per realization of the input function. The reference solution
$G\left(\boldsymbol{y}_{i}\right)$ is determined by FEM-simulations or
measurements. A single data point of the training data set consists of a
triple of the form
$(\boldsymbol{y},\boldsymbol{X},\boldsymbol{u}(\boldsymbol{X}))$. If $P>1$,
the discrete parameter $\boldsymbol{y}$ must be repeated an appropriate number
of times in the data set. For $N$ realizations of the parameter and $P$
evaluations per realization, the training dataset has a total of $N\times P$
entries.

Despite of its simple structure, DeepONet can represent a wide range of
mappings (i.e. it is very expressive) and allows to achieve small
generalization errors. Furthermore, it can be applied very easily to arbitrary
parametrizations.

The physics informed variant PIDeepONet [12] extends the empirical risk
functional by adding a physically motivated term ensuring the compliance with
the PDE in the weak form. For this purpose, the risk functional is extended by
adding the squared residuals of the (nonlinear) differential operators
$\mathcal{N},\mathcal{B}$

$\displaystyle F^{\text{PIDeepONet }}=F^{\text{DeepONet }}$
$\displaystyle+\frac{1}{mN_{y}}\sum_{i=1}^{N_{y}}\sum_{j=1}^{m}\mathcal{N}^{2}\left[G_{\theta}\left(\mathbf{y}_{i}\right)\right]\left(\boldsymbol{X}_{j}\right)$
(17)
$\displaystyle+\frac{1}{N_{b}N_{y}}\sum_{i=1}^{N_{y}}\sum_{j=1}^{N_{b}}\mathcal{B}^{2}\left[G_{\theta}\left(\mathbf{y}_{i}\right)\right]\left(\boldsymbol{X}_{j}\right).$

In return, no FEM reference data is necessary.

Eq. (17) is formulated for the case, where the residual can only be evaluated
at the $m$ discretization points, where the parameters $\boldsymbol{y}$ are
known. Alternatively, an extended data set can be generated which contains
$\left\\{\boldsymbol{y}(X_{1}),\ldots,\boldsymbol{y}(X_{m}),X_{1},\boldsymbol{y}(X_{1}),\ldots,X_{r},\boldsymbol{y}(X_{r})\right\\}$
with $r$ the number of (randomly determined) additional gridpoints.

Figure 3: Information flow in PIDeepONet. Branch and trunk network are linked
by the scalar product.

###  4.2 Fourier-Neural-Operator (FNO)

An FNO [6] represents the solution operator of a PPDE with the help of a
series of Fourier Blocks.

A Fourier block includes several operations as shown in Fig. 4: i) It applies
the Fourier transform $\mathcal{F}$ (in form of the Fast Fourier Transform,
FFT) on its input $\boldsymbol{v}^{t}(\boldsymbol{X})$. ii) It applies a
linear transform $R_{\theta}^{t}$ (parameterized by the NN parameter set
$\theta$) on the lower Fourier modes and filters out the higher modes. iii) It
applies the inverse Fourier transform $\mathcal{F}^{-1}$. iv) In parallel, the
Fourier block applies another linear transform $W_{\theta}^{t}$ on the input
$\boldsymbol{v}^{t}(\boldsymbol{X})$. v) Results of both branches are summed
up and forwarded to the nonlinear activation function $\sigma$.

The input
$[\boldsymbol{X},\boldsymbol{y}(\boldsymbol{X})]\in\mathbb{R}^{N_{dof}+d}$ of
the network is lifted up to a higher dimension $d_{v}$ by a shallow (e.g.
single-layer) FCNN $P$ with linear activation function. Another FCNN $Q$
projects the output of the last Fourier Block onto the output space
$\mathbb{R}^{N_{dof}}$ which results in the FNO output
$\boldsymbol{u}_{\theta}(\boldsymbol{X})$. The whole process is visualized in
Fig. 4 and can be described by the iterative architecture

$\displaystyle\boldsymbol{v}_{0}(\boldsymbol{X})$
$\displaystyle=P(\boldsymbol{y}(\boldsymbol{X}),\boldsymbol{X})$ (18)
$\displaystyle\boldsymbol{v}_{t+1}(\boldsymbol{X})$
$\displaystyle=\sigma\left(W_{\theta}^{t}\boldsymbol{v}_{t}(\boldsymbol{X})+\mathcal{F}^{-1}\left[R_{\theta}^{t}\cdot\mathcal{F}\left[\boldsymbol{v}_{t}\right]\right](\boldsymbol{X})\right)$
$\displaystyle\boldsymbol{u}_{\theta}(\boldsymbol{X})$
$\displaystyle=Q\left(\boldsymbol{v}_{K}(\boldsymbol{X})\right)\,,$

where $K$ is the number of sequential Fourier Blocks.

Figure 4: Schematic representation of the FNO (adapted from [6]). Layers with
nonlinear activation (GELU) are marked in blue.

####  4.2.1 Physically Informed Neural Operator (PINO)

An extension of the FNO to a physically informed neural operator (PINO) [13]
has been investigated as well. FNO computes displacements on an equidistant
grid, so that a DEM-like extension is straight forward. Therefor, the same
potential energy as for the DEM (see Eq. 10) was added to the loss function

$\displaystyle
F^{\text{PINO}}=F^{\mathrm{FNO}}+F^{\mathrm{DEM}}=F^{\mathrm{FNO}}+\int_{B}\left(W\left(\boldsymbol{F}_{\theta}\right)-\boldsymbol{b}\cdot\boldsymbol{u}\right)\,\mathrm{d}V-\int_{\Gamma_{N}}\overline{\boldsymbol{T}}\cdot\boldsymbol{u}\,\mathrm{d}A\quad.$
(19)

##  5 Application to examples in elastostatics

The present section compares the NN methods previously explained by means of
example boundary value problems, including several 1D examples and one 2D
example with two load cases.

###  5.1 The 1D tensile bar

The first example is based on the setup shown in Fig. 5. The bar is clamped at
the left edge $u(-1)=0$ and loaded along its entire length with the force
density $f(X)$. A Neumann boundary condition $P(1)=T$ is applied at the right
edge.

Figure 5: Illustration of the 1D BVP with BC: $u(X=-1)=0$ and $P(X=1)=T$ .

####  5.1.1 Example A

The following energy density is considered

$\displaystyle W(F)=F^{\frac{3}{2}}-\frac{3}{2}F+\frac{1}{2}\quad\text{ with
}\quad F=1+u^{\prime}(X)\,\,.$ (20)

From Eq. (20), the first Piola-Kirchhoff stress reads

$\displaystyle P=\frac{\partial W}{\partial
F}=\frac{3}{2}\left(F^{\frac{1}{2}}-1\right)\Rightarrow-\frac{\partial
P}{\partial
X}=-\frac{3}{4}\frac{1}{\sqrt{1+u^{\prime}}}u^{\prime\prime}(X)=f(X)\,\,.$
(21)

Specifically, the force density $f(X)=X$ is chosen and the load at the free
end is set to zero: $T=0$

$\displaystyle-\frac{3}{4}\frac{1}{\sqrt{1+u^{\prime}}}u^{\prime\prime}(X)=X\quad\text{
with }\quad u(-1)=0,\,\,T=0\Rightarrow u^{\prime}(1)=0\,.$ (22)

The example has the following analytical solution which will be used to
validate the results obtained by the Neural FEM methods

$\displaystyle u(X)$
$\displaystyle=\frac{1}{135}\left(3X^{5}-40X^{3}+105X+68\right)$ (23)
$\displaystyle u^{\prime}(X)$
$\displaystyle=\frac{1}{9}\left(X^{4}-8X^{2}+7\right)\,\,.$ (24)

####  5.1.2 Example B

A linear elastic material is investigated:

$\displaystyle W=\frac{1}{2}\left(u^{\prime}\right)^{2}\quad\Rightarrow\quad
P(X)=u^{\prime}(X)\quad\Rightarrow\quad-u^{\prime\prime}(X)=f(X)$ (25)

##### Example B1

A single load case is analyzed in examples related to PINN and DEM.

$\displaystyle-Eu^{\prime\prime}(X)=f(X)=Q\cdot A\quad\text{ with }\quad
u(-1)=0\,\,\text{ and }\,\,Eu^{\prime}(1)=T.$ (26)

Here, distributed forces take the values $Q=9.395\cdot
10^{4}\mathrm{\leavevmode\nobreak\ Nm}^{-1}\text{ and }T=1.015\cdot
10^{8}\mathrm{\leavevmode\nobreak\ Nm}$. Young’s modulus corresponds to steel
($E=210\cdot 10^{9}\mathrm{\leavevmode\nobreak\ N}\mathrm{\leavevmode\nobreak\
m}^{-2}$) and the cross section surface is $A=1\mathrm{\leavevmode\nobreak\
m}^{2}$.

##### Example B2

For the Neural Operator models, the PPDE is normalized and a parameterization
of the force density $f$ as well as a parameterization of the Neumann boundary
condition are studied. Since the boundary only consists of one point, the
boundary condition can be described by a scalar $\pi_{2}$ for which a uniform
distribution between $[0,1]$ is assumed. The reference solutions are computed
using FEniCS. The BVP is described by

$\displaystyle-\frac{\partial^{2}u}{\partial X^{2}}=f(X)\quad\text{ with
}\quad\left\\{\begin{array}[]{c}u(-1)=0\\\
u^{\prime}(1)=\pi_{2}\end{array}\right.\,\,.$ (29)

###  5.2 The plate – Example C

The selected two-dimensional example deals with a plate made of a Neo-Hookean
material with the energy density

$\displaystyle
W(\boldsymbol{F})=\frac{\mu}{2}\left(I_{1}-2-\ln\mathsfit{J}\right)+\frac{\lambda}{2}(\mathsfit{J}-1)^{2}\,\,.$
(30)

$I_{1}=\operatorname{tr}\left(\boldsymbol{C}\right)$ is the first invariant of
the right Cauchy-Green deformation tensor
$\boldsymbol{C}=\boldsymbol{F}^{T}\boldsymbol{F}$ and
$\mathsfit{J}=\operatorname{det}(\boldsymbol{F})$ the determinant of the
deformation gradient. The corresponding derivatives are
$\frac{\partial\mathsfit{J}}{\partial\boldsymbol{F}}=\mathsfit{J}\boldsymbol{F}^{-1}$
and
$\frac{\partial\operatorname{tr}\left(\boldsymbol{F}^{\top}\boldsymbol{F}\right)}{\partial\boldsymbol{F}}=2\boldsymbol{F}$.
The symbols $\lambda$ and $\mu$ denote the Lamé constants.

For the energy density, Eq. (30), the 1st and 2nd Piola-Kirchhoff stress
tensor are given by:

$\displaystyle\boldsymbol{P}=\frac{\partial
W}{\partial\boldsymbol{F}}=\mu\boldsymbol{F}+(\lambda\ln\mathsfit{J}-\mu)\boldsymbol{F}^{-\top}\quad\text{
and
}\quad\boldsymbol{S}=\boldsymbol{F}^{-1}\cdot\boldsymbol{P}=\mu\boldsymbol{I}+(\lambda\ln\mathsfit{J}-\mu)\boldsymbol{C}^{-1}\quad.$
(31)

Figure 6: Geometry of the 2D BVP.

The studied example is shown in Fig. 6. The plate is clamped at the left edge
and the Neumann boundary conditions are set at the right edge. The components
of the kinematic fields in Cartesian coordinates are given by

$\displaystyle[\boldsymbol{u}]=\left[\begin{array}[]{l}u_{x}\\\
u_{y}\end{array}\right]\quad\Rightarrow\quad[\boldsymbol{F}]=\left[\begin{array}[]{ll}F_{xx}&F_{xy}\\\
F_{yx}&F_{yy}\end{array}\right]=\left[\begin{array}[]{cc}1+u_{x,x}&u_{x,y}\\\
u_{y,x}&1+u_{y,y}\end{array}\right]\quad.$ (38)

From Eq. (31), the first and second Piola-Kirchhoff stress tensor are
calculated as follows

$\displaystyle[\boldsymbol{P}]$
$\displaystyle=\left[\begin{array}[]{cc}P_{xx}&P_{xy}\\\
P_{yx}&P_{yy}\end{array}\right]=\mu\left[\begin{array}[]{cc}F_{xx}&F_{xy}\\\
F_{yx}&F_{yy}\end{array}\right]+\frac{\lambda\ln\mathsfit{J}-\mu}{\mathsfit{J}}\left[\begin{array}[]{cc}F_{yy}&-F_{yx}\\\
-F_{xy}&F_{xx}\end{array}\right]\quad,$ (45) $\displaystyle´[\boldsymbol{S}]$
$\displaystyle=\left[\begin{array}[]{ll}S_{xx}&S_{xy}\\\
S_{yx}&S_{yy}\end{array}\right]=\left[\begin{array}[]{cc}\mu&0\\\
0&\mu\end{array}\right]+\frac{\lambda\ln\mathsfit{J}-\mu}{\mathsfit{J}^{2}}\left[\begin{array}[]{cc}C_{yy}&-C_{xy}\\\
-C_{yx}&C_{xx}\end{array}\right]\quad.$ (52)

Moreover, an equivalent stress is calculated as in [8] and used in contour
plots (Section 5.5.2)
$S_{E}=\sqrt{0.5\left(\left(S_{xx}-S_{yy}\right)^{2}+S_{xx}^{2}+S_{yy}^{2}\right)+3S_{xy}}$.
Two load cases are investigated for the setup described.

##### Example C1

The first load case deals with the vertical load
$\overline{\boldsymbol{T}}=-5\boldsymbol{e}_{y}$ .

##### Example C2

The second load case is uniaxial tension with
$\overline{\boldsymbol{T}}=50\boldsymbol{e}_{x}$.

###  5.3 Error measure

With the solution operator of the PPDE $G:\mathcal{Y}\rightarrow\mathcal{S}$
and its NN approximation $G_{\theta}$, the average relative $L_{2}$ error for
the $N$ test data sets is calculated for the Neural Operator methods [22, 6,
10] as:

$\displaystyle\epsilon_{\mathrm{rel}}=\frac{1}{N}\sum_{j=1}^{N}\frac{||G_{\theta}(\boldsymbol{y}_{j})-G(\boldsymbol{y}_{j})||_{L_{2}}}{||G(\boldsymbol{y}_{j})||_{L_{2}}}\quad.$
(53)

With the Neural FEM, only one concrete BVP can be analyzed at once. Then, the
relative $L_{2}$ error is calculated based on the solution for the
displacement field

$\displaystyle\epsilon_{\mathrm{rel}}=\frac{||\boldsymbol{u}_{\theta}-\boldsymbol{u}||_{L_{2}}}{||\boldsymbol{u}||_{L_{2}}}\quad.$
(54)

In both cases, the determination of the relative $L_{2}$-error requires the
computation of the $L_{2}$-norm which is approximated by the discrete
$L_{2}$-norm. On an equidistant lattice
$\left\\{\boldsymbol{X}_{i}^{\text{equi }}\right\\}_{i=1}^{N}$, the discrete
$L_{2}$-norm is calculated as

$\displaystyle||\boldsymbol{f}||_{L_{2,d}}^{2}=\Delta
V\sum_{i=1}^{N}||\boldsymbol{f}(\boldsymbol{X}_{i}^{\mathrm{equi}})||_{2}^{2}=\Delta
V\sum_{i=1}^{N}\sum_{j=1}^{d}f_{j}^{2}(\boldsymbol{X}_{i}^{\mathrm{equi}})\,\,,$
(55)

with the volume (in 2D: surface area) of each lattice unit $\Delta V$.

###  5.4 Numerical integration

The discrete $L_{2}$-norm is based on a simple Riemann sum with error order
$\mathcal{O}(\Delta V)$. For the approximation of the risk functional, e.g.,
in connection with the calculation of the potential energy in the DEM, other
integration methods have to be considered. Two classical methods are the Monte
Carlo (MC) integration and the trapezoidal rule. The trapezoidal rule requires
partitioning of the integration domain into polytopes. In the simplest case,
these would be hypercubes on an equidistant grid. In [18], an integration
method based on the Delaunay triangulation is proposed. The trapezoidal rule
still applies, where $\bar{f}_{i}$ is the average value over the $i$-th
simplex (e.g. triangles). The two polytopes for integration in 2D are shown as
examples in Fig. 7. Let $V$ be the volume of the integration domain and
$\bar{f}_{i}$ the average of $f$ over the corners of the $i$-th polytope with
$i\in[1,N]$ and $N$ the number of vertices. Then, the integral approximations
are generally given by Eq. (56a) for Monte-Carlo and Eq. (LABEL:eq:4-13b) for
the trapezoidal rule.

$\displaystyle
I_{\mathrm{MC}}(f)=\frac{V}{N}\sum_{i=1}^{N}f\left(\boldsymbol{X}_{i}\right)\quad\text{(a)}$
$\displaystyle\qquad I_{\mathrm{T}}(f)=\sum_{i=1}^{N}\bar{f}_{i}\,\,\Delta
V_{i}\quad\text{(b)}\quad.$ (56)

Figure 7: Examples for polytopes.

Three methods are investigated to select the grid points: equidistant grid
points, pseudo-random numbers and quasi-random numbers (Latin Hypercube
Sampling, LHS). Exemplarily, we compare the absolute error of the potential
energies in the nonlinear 1D setup (Example A). The trapezoidal rule with 100
000 grid points is assumed as a quasi-exact comparison value. The MC
integration with equidistant grid points reduces to a simple Riemann sum. The
results for 100 and 1000 grid points, respectively, are summarized in Fig. 8.
A characteristic distribution of grid points is shown in Fig. 9.

Figure 8: Comparison of Monte Carlo (MC) integration and the trapezoidal rule
with $n=100$ and $n=1000$ grid points for equidistant sampling (equi), pseudo-
random uniform sampling (unif) and Latin Hypercube Sampling (LHS). Figure 9:
Quasi-random (top), pseudo-random (bottom), and equidistant (middle) grid
points for numerical integration (1D case).

Due to the larger integration error, uniform pseudo-random sampling is not
considered further. In the following, "random" sampling always refers to
quasi-random LHS. Fig. 8 also shows that the trapezoidal rule is consistently
more accurate than MC integration.

####  5.4.1 Technical Implementation

In this work, we use PyTorch (version: 1.11.0), which contains the optimizer
L-BFGS that is employed in some of the investigated methods. The computations
have been run on an Intel Core i5-7200U mobile processor. A mobile NVIDIA
GeForce GTX 950m is used as graphics card.

The parameters of the NN are always randomly initialized, so that the results
underly statistical variations.

###  5.5 Neural FEM

####  5.5.1 PINN

##### Example A

In the numerical experiments related to Example A, the NN architecture
$[1,10,1]$ is always used. Two optimizers (L-BFGS, SGD), different numbers of
collocation points (100 and 1000 points) and different computational
accuracies (single precision FP32 and double precision FP64) are compared
(Fig. 10).

The SGD is repeated for 10000 epochs and the L-BFGS for 15 epochs, each with
the default parameters of the methods. The discrepancy in the required number
of epochs is reflected by the run time, which is about
$43\mathrm{\leavevmode\nobreak\ s}$ for the SGD compared to about
$0.300\mathrm{\leavevmode\nobreak\ s}$ for the L-BFGS. Obviously, a second-
order method (such as the L-BFGS) can greatly reduce the number of necessary
iterations. In each epoch, the complete data set is used (Full-Batch). Despite
the same information being provided to the NN, the L-BFGS method performs
better on average than the SGD. Fig. 10 also shows that the reduction in total
error relative to the number of collocation points quickly goes into
saturation. The difference between $N=100$ and $N=1000$ is only about
$6.500\cdot 10^{-6}$ for the L-BFGS.

Moreover, no significant increase of the total error measured in the relative
$L_{2}$ norm is shown when computing on single precision. This can greatly
reduce the computation time on commercially available graphics cards that are
optimized for single precision computing. However, this needs to be confirmed
in the further research for more complex problems.

The best results were obtained with the tangent hyperbolic (Tanh) activation
function. Other activation functions, such as the Rectified Linear Unit (ReLU)
or the Exponential Linear Unit (ELU) do not converge or converge very poorly
against the analytical solution of the problem. The calculated displacements
for different activation functions are comparatively shown in Fig. 11. ReLU
and ELU could not be optimized with L-BFGS. Therefore, only the results after
optimization with Adam are shown. The second derivative of the approximation
with ReLU activation is everywhere constantly zero, (except for the point at
the kink). This destroys the information in the residual and a training of the
network must necessarily fail. Hence, for the following studies, Tanh was
always applied as the activation function. Other activation functions are not
considered in the present contribution. However, in the literature, the
composite function $\max\left(0,x^{3}\right)$ [24], the Swish activation
$zS(\beta z)$ (where $S$ denotes the sigmoid activation function) [25] and
GELU [6] have been successfully employed.

The training of PINN with SGD took about $40\mathrm{\leavevmode\nobreak\ s}$
for 10 000 epochs.

Figure 10: Relative $L_{2}$ error and run times in s for two optimization
methods L-BFGS, SGD with 100 and 1000 collocation points and FP32/ FP64
accuracy.

(a) Displacements

(b) Absolute error

Figure 11: Calculated displacements with 100 collocation points and different
activation functions.

##### Example C

Amongst the conventional PINN representatives, the DCM is the easiest to
implement for 2D problems and thus chosen to apply to Example C (Section
3.1.1).

Within the domain, the balance of linear momentum reads
$\nabla\cdot\boldsymbol{P}=\mathbf{0}$, which already is a residual form for
approximations of $\boldsymbol{P}$. The Neumann boundary conditions are given
by $\boldsymbol{P}\cdot\boldsymbol{N}=\overline{\boldsymbol{T}}$, the
Dirichlet constraints are incorporated directly by the application of a
transformation on the output of the NN (Eq. (9)).

The architecture of the network is specified as $[2,30,30,2]$ and on the
Neumann boundary part, 900 random collocations points are chosen. 4000
collocation points are used within the body. L-BFGS with learning rate 1.0 and
Line Search with Wolfe condition is applied as optimizer.

Other than the DEM, the DCM does not converge towards the reference solution
for load case C1. A comparison with the DEM shows that the boundary conditions
are not appropriately learned by the DCM.

[*]Wang2021a [19] discuss that the training of PINNs may fail due to numerical
inaccuracies if the contributions to the loss value – one portion from the
residual and the other portion from the Neumann boundary part – or their
gradients w.r.t to the NN parameters are in vastly different orders of
magnitude. In the case of Example C, the loss value in the DCM consists of the
portion from the residual with the value $0.119$ and the portion from the
Neumann boundary part with the value $1.121$. The gradients of each
contribution w.r.t. the NN parameters are relatively uniformly distributed
(Fig. 12). Moreover, the correct boundary conditions are not learned even if
the loss portion of the residual is excluded (manually set to 0).

Figure 12: Histogram over the gradients of empirical risk. First hidden layer
(left), output layer (right).

Hence, the failure of the DCM on the example cannot be explained by this kind
of numerical inaccuracies. It is possible that the optimizer gets stuck in a
local minimum, but this behavior needs a further investigation.

####  5.5.2 cPINN

##### Example A

For many applications, the accuracy that can be achieved with a classical PINN
is not sufficient. cPINN was developed to improve the accuracy by avoiding the
squaring of the residual [9, 26]. Instead, it employs Adaptive Competitive
Gradient Descent (ACGD) as optimization procedure whose Python implementation
is publicly available [27]. Furthermore, Tanh in the PINN part and ReLU in the
discriminator network are chosen as the activation functions.

The architecture of the PINN is chosen with 10 neurons in the hidden layer
(architecture: $[1,10,1]$ ) for comparability with the conventional PINN/ DCM.
The layer width $h$ of the discriminator, on the other hand, was varied
(architecture: $\left[1,20,2\right]$ vs. $\left[1,50,2\right]$). Initial
experiments have shown that the output of the discriminator
$\boldsymbol{d}_{\phi}$ must be separated for points in the domain and on the
boundary
$\boldsymbol{d}_{\phi}=\left(\boldsymbol{d}_{\phi}^{B},\boldsymbol{d}_{\phi}^{\Gamma}\right)$.
Therefore, the output layer contains 2 neurons. The option of separating both
outputs of the discriminator into independent subnetworks is also tested.
However, this did not result in any improvement. Based on these results, only
the first variant with the smaller number of NN parameters is considered
further. All calculations are performed with double precision.

Figure 13: Relative $L_{2}$ errors and run times in s of cPINN with 100 and
1000 collocation points and 20 and 50 neurons in the hidden layer of the
discriminator.

Fig. 13 shows the results for accuracy and run times, from 100 runs with
random NN parameter initializations in form of a box plot. It can be seen that
the accuracy is only moderately affected by the number of collocation points.
However, the width of the discriminator has a significant impact on the
training result. Increasing the number of neurons in the hidden layer from 20
to 50 reduces the error by an order of magnitude. It is not entirely clear why
such a large discriminator network is necessary. Moreover, the training is
relatively slow, taking about 2 min (up to 3 min for 1000 collocation points)
for about 6000 epochs.

The improvement of up to 2 orders of magnitude reported by [9] cannot be
demonstrated here. This may be because the pathologies related to training
PINNs [20, 19], which cPINN addresses, do not arise in this simple example. On
the other hand, the regularization by the residual leads to a complex energy
landscape of the optimization procedure, making optimization more difficult
[20]. Moreover, the material law sometimes causes the optimization process to
abort if the network parameters have been initialized unfavorably. This
problem can be solved by reducing the range for random sampling of the initial
parameter values.

In addition, the different weighting of the summands of the empirical risk can
lead to different magnitudes of the gradients of the loss function w.r.t. the
NN parameters, which can impair the NN parameter optimization [19, 12]. In the
present case, the gradients for the residual within the domain
$\nabla_{\theta}F^{\mathrm{PINN},B}$ are much larger than the gradients for
the constraints $\nabla_{\theta}F^{\mathrm{PINN},\Gamma}$. The optimization
procedure is therefore driven more strongly toward a solution that reduces the
residual while allowing for deviations from the constraints. As a result, the
optimization procedure converges toward a plausible solution, but one that
does not satisfy the boundary conditions. For the complex architecture
$[1,50,50,1]$, the gradients of the residual and boundary portion of the loss
function w.r.t. the NN parameters of the first hidden layer and the output
layer are shown in Fig. 14. However, not much discrepancy is detected in the
order of magnitude of the gradient values for both portions of the loss
function.

Figure 14: Histogram of the gradients of the optimizable parameters. First
hidden layer (left), output layer (right).

####  5.5.3 DEM

With the architecture $[1,10,1]$, DEM training runs on average twice as fast
compared to classical PINN. For deeper networks, this fact is amplified, as
will be shown in the analysis of Example C. Our own DEM implementation is
based on the public source code [8] and is extended to allow training with
Monte Carlo integration on quasi-random grid points. Tanh is used as the
activation function.

##### Example A

We enforced the geometric boundary conditions by the transformation

$\displaystyle u(X)=(1+X)z_{\theta}(X)\,\,,$ (57)

where $z_{\theta}(X)$ is the output of the NN. This way, always zero
displacement is calculated at the clamped end ($X=-1$). A comparison of the
resulting relative errors in displacements and strains as well as the run
times is shown in Fig. 15.

Figure 15: Relative $L_{2}$ errors and run times of DEM for trapezoidal rule
(T) and Monte Carlo integration (MC) on 100 and 1000 grid points each.

Two integration methods (Monte Carlo integration and trapezoidal rule) are
compared, each on different sets of randomly selected support points (100 and
1000, respectively). For the trapezoidal rule, the support points are sorted
and the boundary points are explicitly considered. The optimizer is the L-BFGS
with learning rate $1.000$. This proves to be very efficient and approaches
the solution after only 15 epochs.

The results furthermore illustrate that the use of single floating point
accuracy (FP32) leads to only a slight decrease of accuracy, similar as seen
with the conventional PINN. However, the run time even increases with FP32,
what indicates slower convergence of the optimization procedure.

##### Example B

The study of Example B reveals a pathology of the DEM that has not appeared in
Example A. Its effect can be seen in Fig. 16. This pathology can be attributed
to overfitting [11], since the potential energy, unlike the squared residual,
has no regularizing effect. In fact, early stopping significantly reduced the
influence of overfitting. Alternatively, overfitting could be avoided by
increasing the number of grid points. For 1000 grid points, hardly any
overfitting occurred without a need for early stopping.

Figure 16: Displacement results of DEM compared to FEM. Displacements with
overfitting (left), Absolute error in displacements (right); rel. $L_{2}$
error: $1.358\mathrm{e}+00$ .

##### Example C

The approximation of the energy functional can be done analogously to the 1D
example by means of different integration techniques. However, the Monte Carlo
integration is the simplest option to implement. The incorporation of the
boundary conditions, network architecture and optimizer algorithm as well as
the load case are chosen similarly as for the DCM (Section 5.5.1).

10 000 collocation points within the body bulk are selected. The external
energy is evaluated by using 200 random points on the right edge of the plate.
The results for load case C1 and error are presented in Figs. 17 and 18,
respectively.

(a) Displacements

(b) Equivalent stresses

Figure 17: Calculated displacements and equivalent stresses for the vertical
load case.

The relative $L_{2}$ error in the displacements is only 0.002 and in the
equivalent stress 0.053. Fig. 18 shows that the error in the equivalent
stresses is concentrated at the restraint. The stress peaks at the critical
points are not completely resolved by the NN.

(a) Error in displacements

(b) Error in equivalent stresses

Figure 18: Absolute error in displacements and equivalent stresses for the
vertical load case.

The integration method does not cause this error, which is confirmed by a
second calculation that uses trapezoidal rule. The relative $L_{2}$ errors in
the case are 0.003 for displacements and 0.039 for equivalent stresses. Again,
an equidistant grid with 10 000 collocation points is used. No significant
improvement is obtained by the more accurate integration procedure, so the
influence must be considered small.

For load case C2, the same procedure with MC integration is carried out. The
results and errors are presented in Figs. 19 and 18, respectively. The
relative $L_{2}$ errors are 0.005 and 0.019 in this case.

In conclusion, the NN is able to approximate the character of the solution of
the BVP. However, relatively large errors are found at the restraint.
According to [18], the same problems arise for a PINN trained with the squared
residual. The resolution of fine features of the stress and displacement
fields seem still to be a challenge for future work.

(a) Displacements

(b) Equivalent stresses

Figure 19: Calculated displacements and equivalent stresses for uniaxial
tension.

(a) Error in Displacements

(b) Error in Equivalent stresses

Figure 20: Absolute error in the calculated displacements and equivalent
stresses for uniaxial tension.

####  5.5.4 Transfer learning (TF)

One possibility to enhance the performance of Neural FEM is by means of
transfer learning e.g. in case of a varying Neumann boundary condition. This
is illustrated on Example B2, where $\pi_{2}$ is changed by only a small
amount in each iteration. Then the NN trained on the previous $\pi_{2}$ value
is already a good approximation for its subsequent value. Hence, the NN
parameter values can be copied to the new NN to reduce the number of required
learning epochs.

Applied to Example B2 (linear elastic material), the average run time could be
reduced by a factor of four. Similar time savings have been documented in [14]
in the study of plastic deformations. The relative $L_{2}$ errors are shown in
Fig. 21. The average $L_{2}$ error is 3.145 $\cdot 10^{-5}$, which is about an
order of magnitude smaller than the error with DeepONet. The training duration
is reduced from $23\mathrm{\leavevmode\nobreak\ s}$ to about
$5\mathrm{\leavevmode\nobreak\ s}$.

Figure 21: Effect of transfer learning (TF). Comparison of $L_{2}$ errors over
the parameter space.

####  5.5.5 Initialization of a conventional FE solver

As indicated in the literature and the results above, relative $L_{2}$ errors
of approx. $10^{-4}$ are usually achieved. The computed solution could then be
submitted as initialization to a traditional FEM solver in order to improve
the accuracy.

Applied to Example A, a FEniCS calculation, that is conventionally initialized
with all displacements as zero, runs for four iterations. With the solution of
the DEM with trapezoidal rule and 1000 collocation points as initialization
values for the FEM solver, the calculation is accelerated by a factor of two,
only half of the iterations until convergence are needed.

Thus, Neural FEM results can be employed as a potential way to speed up an FEM
simulation in settings where the Neural FEM is not yet able to completely
replace the FE simulation.

###  5.6 Neural Operators

The operator methods are examined on example Example B2 (Section 5.1.2) – a
tensile bar with clamping restraint at the left side and a free end at the
right side. In this simple test case, analytical solutions are available to
calculate the error for the parameterization of the Neumann boundary condition
in case of fixed force density $f(X)\equiv 1$. Combinations of both varying
the force field and the boundary conditions have not been carried out in the
present work.

Training and test data for the varying force field are generated by means of
Gaussian Processes with squared exponential covariance (correlation length
$l=0.100$). The free end of the bar at the right side yields the parameter
$\pi_{2}=0$. For training, 1000 data sets and for the tests 100 data sets have
been generated with FEniCS on an equidistant grid for $X\in[-1,1]$ with 1024
grid points and quadratic shape functions. The relative $L_{2}$ errors of the
FEM simulation are several orders of magnitude smaller than expected from the
NN methods (displacements approx. $10^{-10}$; strains approx. $10^{-7}$) and
should not influence the survey. One of the resulting data sets is illustrated
in Fig. 22.

(a) Displacements

(b) Stresses $P=u^{\prime}$ and force density $f$

Figure 22: FEM solution for a random realization of the force density $f$.

####  5.6.1 DeepONet and PIDeepONet

##### Numerical setup

For the DeepONet, the data sets need to be preprocessed since the Neural
Operator methods work with $P$ random collocation points that change between
the evaluations of the loss function, whereas the reference solutions are
produced be FEM on a fix mesh with $m$ equidistant points. Hence, the
realizations of the force fields are projected onto the FE mesh. Then, the FE
results are interpolated and evaluated at the random collocation points. The
parametrization of Neumann boundary condition has only been investigated with
DeepONet, where $\pi_{2}\sim U[0,1]$ is assumed.

The source codes for DeepONet and PIDeepONet have both been published on
Github [28, 29]. The architectures for the subnets were specified as
$[20,100,100]$ for the branch net and $[1,100,100]$ for the trunk net. The
branch is set up with layer width $m=20$. The chosen activation function are
ReLU in DeepONet and Tanh in PIDeepONet. A similar architecture has also been
suggested in [10, 12]. The L-BFGS optimizer with linesearch (strong Wolfe
condition) and learning rate $1.0$ is applied as suggested in [20], for 120
epochs. However, L-BFGS is very memory consuming, so $historysize=50$
(default: 100) is set.

1000 load cases are used for training and 100 for testing. In order to reduce
the training effort, from the 1024 grid points only 8 are randomly chosen per
load case. Hence, the training data set size reduces to 8000 data points. The
error reduction of considering more grid points per load case quickly
saturates so that this reduction is admissible. Approximately 50 random points
can be estimated as the saturation limit for 1000 training data sets [10].
Training is conducted for both methods, Full-Batch and Mini-Batch (batch size:
1000) .

Since PIDeepONet includes the whole information about the PPDE (similar to the
PINN models) in the loss function, no reference solutions by means of FEM are
necessary. In exchange, the loss functions needs to be constructed anew for
each PDE.

For the DeepONet, additionally an alternative loss function is investigated.
It employs the relative $L_{2}$ error instead of the Mean Square Error – MSE
(Eq. 5).

##### Results

Representative results for the displacement $u$ and strains $u^{\prime}$ over
the bar length obtained with DeepONet and PIDeepONet are shown in Fig. 23.
Both methods match the displacements relatively well, but DeepONet has visible
deviations in the strains. In particular, the non-smooth curve of the strain,
which is the spatial derivative of the displacement, can be attributed to the
ReLU activation function, which has a discontinuous derivative.

(a) Displacements

(b) Absolute error for displacements

(c) Strains

(d) Absolute error for strains

Figure 23: Displacements and strains with corresponding errors for an
exemplary element of the test data set with
DeepONet.

The errors for the displacements and strains are shown in Table 1. The usage
of the residual in the empirical risk in PIDeepONet improves the accuracy in
the strains about one order of magnitude compared to DeepONet.

The effect of floating point accuracy on $\epsilon_{\mathrm{rel}}$ is small,
similar as seen for Neural FEM. Overall, training with Full-Batch on FP32
performs best.

The runtimes of models, each trained with Full-Batch and Mini-Batch (batch
size 1000), are compared in Table 2. The run times for the DeepONet are
significantly lower than for the FNO with comparable accuracies. Using the
relative $L_{2}$ error instead of the MSE reduces the convergence rate,
requiring more iterations and increasing the run time. The difference between
the mean and median is reduced, but no significant effect on the error is
found.

| Mean(FP32) | Median(FP32) | Mean(FP64) | Median(FP64)
---|---|---|---|---
DeepONet(MSE), displacements | $5.070\cdot 10^{-3}$ | $3.144\cdot 10^{-3}$ | $5.076\cdot 10^{-3}$ | $3.023\cdot 10^{-3}$
DeepONet(L2), displacements | $-$ | $-$ | $2.109\cdot 10^{-3}$ | $1.329\cdot 10^{-3}$
PIDeepONet, displacements | $3.268\cdot 10^{-3}$ | $2.008\cdot 10^{-3}$ | $2.450\cdot 10^{-3}$ | $1.343\cdot 10^{-3}$
DeepONet(MSE), strains | $4.022\cdot 10^{-2}$ | $2.831\cdot 10^{-2}$ | $4.023\cdot 10^{-2}$ | $2.833\cdot 10^{-2}$
DeepONet(L2), strains | $-$ | $-$ | $5.393\cdot 10^{-2}$ | $4.492\cdot 10^{-2}$
PIDeepONet, strains | $3.055\cdot 10^{-3}$ | $2.394\cdot 10^{-3}$ | $2.561\cdot 10^{-3}$ | $1.936\cdot 10^{-3}$
Table 1: Mean error values over the test data set with DeepONet and PIDeepONet. | Full (FP64) | Mini (FP64) | Full (FP32) | Mini (FP32)
---|---|---|---|---
DeepONet(MSE) | $82\mathrm{\leavevmode\nobreak\ s}$ | $188\mathrm{\leavevmode\nobreak\ s}$ | $70\mathrm{\leavevmode\nobreak\ s}$ | $120\mathrm{\leavevmode\nobreak\ s}$
DeepONet(L2) | $112\mathrm{\leavevmode\nobreak\ s}$ | $200\mathrm{\leavevmode\nobreak\ s}$ | $-$ | $-$
PIDeepONet | $933\mathrm{\leavevmode\nobreak\ s}$ | $951\mathrm{\leavevmode\nobreak\ s}$ | $607\mathrm{\leavevmode\nobreak\ s}$ | $-$
Table 2: Run times of the optimization procedure with Full-Batch and Mini-
Batch (1000 elements), respectively.

Fig. 24 shows the loss histories from the optimization with the DeepONet and
PIDeepONet, respectively. Overall, the relative $L_{2}$ error on the test data
set is smaller for Full Batch training. The difference in resulting accuracy
between the two methods can be attributed to the training error alone.
PIDeepONet converges significantly slower and yields worse accuracy than
DeepONet. With Full-Batches, PIDeepONet even converges to a local instead of
the global minimum. The poor convergence of the PIDeepONet demonstrates the
significantly more complex optimization task, where the DeepONet makes use of
the explicitly obtained FEM results as training data sets.

Figure 24: Training loss histories after optimizing PIDeepONets with Full-
Batch and Mini-Batch.

##### Potential Energy in the loss function

The results of Neural FEM (Section 5.5.2) and [25, 30] suggest, that replacing
the squared residual in the loss function of PIDeepONet by the potential
energy can make the optimization problem easier to solve. Hence, such a method
should be more robust and efficient and make training feasible even where
training of PINNs fails. However, with 100 random realizations of the force
field and 15 collocation points per realization as suggested in [25], this
method does not converge for Example B2.

##### Influence of intialization

The default initialization by PyTorch sets the weights and bias by randomly
sampling from a uniform distribution. For
$\boldsymbol{W}\in\mathbb{R}^{N_{k}\times N_{k-1}}$ and
$\boldsymbol{b}\in\mathbb{R}^{N_{k}}$, it holds:

$\displaystyle W_{ij},b_{i}\sim U[-\sqrt{k},\sqrt{k}]\quad\text{ with }\quad
k=\frac{1}{N_{k-1}},$ (58)

where $N_{k}$ denotes the width of the $k$-th layer. [*]Glorot2010 [31] and
[*]Wang2021b [12] suggest that the convergence of the NN can be accelerated by
the Glorot initialization. Let $N\left[\mu,\sigma^{2}\right]$ be a normal
distribution with mean $\mu=0$ and variance $\sigma^{2}$. This yields

$\displaystyle b_{i}=0\quad\text{ and }\quad W_{ij}\sim
N\left[0,\sigma^{2}\right]\quad\text{ with
}\quad\sigma=\sqrt{\frac{2}{N_{k}+N_{k-1}}}$ (59)

for the parameters.

The errors of the models on single precision with Full-Batch optimization and
L-BFGS are shown in Table 3. The error of DeepONet for the displacements
becomes only slightly smaller. Also, no large difference in run times was
observed. The Glorot initialization was developed mainly for deep learning
applications and does not have much impact on the shallow networks (with one
hidden layer) used here.

| displacements | strains | running time
---|---|---|---
DeepONet | $1.629\cdot 10^{-3}$ | $3.594\cdot 10^{-2}$ | $64\mathrm{\leavevmode\nobreak\ s}$
PIDeepONet | $5.338\cdot 10^{-3}$ | $4.501\cdot 10^{-3}$ | $601\mathrm{\leavevmode\nobreak\ s}$
Table 3: Mean error values over the test data set, with Glorot initialization.

##### Neumann boundary parameterization

The (PI-)DeepONet and the FNO have been specifically designed for
approximating mappings between function spaces. One advantage of the DeepONet
over the FNO is that it is very easy to apply to arbitrary parameterizations.
For example, the PIDeepONet can be used to parameterize the Neumann boundary.
Let again be given the dimensionless problem

$\displaystyle-u^{\prime\prime}(X)=1\quad\text{ with }\quad
u(-1)=0,\,\,u^{\prime}(1)=\pi_{2}\quad.$ (60)

The analytical solution is given by

$\displaystyle u(X)=\frac{3}{2}+X-\frac{1}{2}X^{2}+\pi_{2}(X+1)\quad.$ (61)

In the following, only the scalar variable $\pi_{2}$ needs to be varied. The
solution operator $G:I\rightarrow\mathcal{S}$ is sought, where $I=[0,1]$ is
fixed and $\mathcal{S}$ denotes the space of admissible deformations.

For each of the subnetworks, a hidden layer with 50 neurons is used. Their
architectures are thus given by $[1,50,50]$. Tanh is used as the activation
function according to the experience in neural FEM. The NN is trained over 40
epochs with L-BFGS on 10 000 training data set entries, and 1000 validation
data set entries.

The data set is built by selecting a single random collocation point for each
of the 100 realizations of $\pi_{2}$ ($P=1$) and 1000 grid points each on the
interval $[-1,1]$. No early stopping is used. The relative $L_{2}$ error over
the whole interval $I$ is shown in Fig. 25.

Figure 25: Distribution of the relative $L_{2}$ error over the parameter space
for displacements $u$ and strains $\varepsilon$.

The training of the mesh took approximately $118\mathrm{\leavevmode\nobreak\
s}$. The mean relative $L_{2}$ error for the displacements is $1.695\cdot
10^{-4}$ and $2.463\cdot 10^{-4}$ for the strains. For the calculation of the
complete test data set, the PIDeepONet took $0.033\mathrm{\leavevmode\nobreak\
s}$. This highlights the difference between the short run time of the
inference and the large computational effort for the training.

With L-BFGS and 1000 collocation points (to minimize the risk of overfitting)
the training of the DEM on 100 realizations took about 23 s. The higher
training effort of the operator model is profitable only for about 500
realizations of $\pi_{2}$ and more.

####  5.6.2 FNO and PINO

The FNO architecture proposed in [6] is applied for Example B2. Here, the
hyperparameter representing the hidden layer width $d_{v}=64$ (Section 4.2)
results in 549 569 NN parameters. Secondly, a smaller architecture with
$d_{v}=12$ is set up which results in 20 885 NN parameters. This is comparable
to the DeepONet architecture with 22 500 parameters.

Padding was considered as suggested in [6], since the considered example has
non-periodic boundary conditions in the input functions.

Gaussian Error Linear Unit (GELU) is used as the activation function,

$\displaystyle\operatorname{GELU}(x)=x\Phi(x)=\frac{x}{2}\left[1+\operatorname{erf}\left(\frac{x}{\sqrt{2}}\right)\right],$
(62)

where $\Phi(x)$ denotes the standard normal distribution.

Adam with decreasing learning rate (initial value 0.001, reduction factor 0.5
every 50 epochs) and weight decay $\lambda=10^{-4}$ is used as optimizer. The
other optimizer parameters are kept as the PyTorch default values. The
relative $L_{2}$ error is used for the loss function.

Similar to the DeepONet, a training data set with 1000 entries and a test data
set with 100 entries is used, 500 training epochs are carried out.

In the present work, the spatial derivatives in the loss value for the elastic
strain energy and body forces are approximated by a second order central
difference method instead of employing the autograd feature. This can
significantly reduce the computational effort, since the number of parameters
in the NN is usually much greater than the number of grid points. Exact
derivation methods are discussed in more detail in [13].

The results for an exemplary test load case are shown in Figs. 26 and 27. The
absolute errors in the displacements computed by FNO and PINO are relatively
similar, but the strains at the endpoints of the bar as computed by the pure
FNO show significant errors, rendering the solution practically unusable. The
mean errors in the displacements and strains on the whole test data set are
shown in Table 4. The inclusion of potential energy does not yield a
significant effect on the accuracy for the displacements. The unphysical
oscillations of the calculated strains at the edges of the computational
domain are eliminated (Fig. 27). The only drawback of PINO is the
discretization dependence of the numerical derivative. After training, the
error is no longer constant over different discretization levels which is
analyzed in Section 5.6.3.

(a) Displacements

(b) Absolute error

Figure 26: Displacements calculated with FNO and PINO (for an exemplary load
case of the test data set).

(a) Strains

(b) Absolute error

Figure 27: Strains calculated with FNO and PINO (for an exemplary load case of the test data set). | Mean ($d_{v}=64$) | Median ($d_{v}=64$) | Mean ($d_{v}=12$) | Median ($d_{v}=12$)
---|---|---|---|---
FNO, displacements | $2.438\cdot 10^{-3}$ | $1.447\cdot 10^{-3}$ | $1.102\cdot 10^{-2}$ | $6.428\cdot 10^{-3}$
PINO, displacements | $2.688\cdot 10^{-3}$ | $1.528\cdot 10^{-3}$ | $7.011\cdot 10^{-3}$ | $4.436\cdot 10^{-3}$
FNO, strains | $1.145\cdot 10^{-1}$ | $9.278\cdot 10^{-2}$ | $3.488\cdot 10^{-1}$ | $3.119\cdot 10^{-1}$
PINO, strains | $1.382\cdot 10^{-2}$ | $1.053\cdot 10^{-2}$ | $4.941\cdot 10^{-2}$ | $3.835\cdot 10^{-2}$
Table 4: Mean and median values for errors in displacements and strains from
FNO and PINO with $d_{v}=64$ and $d_{v}=12$.

The errors are concentrated at the edge, which indicates that the choice of
boundary conditions for the FNO might be unfavorable, although this approach
should be able to map non-periodic boundary conditions when padding data
arrays with zeros. For comparison, a new data set with periodic force fields
and a bar clamped on both sides has been modeled. For this purpose, the data
set for the 1D-Burgers problem is assumed [6]. The periodic initial conditions
are interpreted as force fields. The results that yield the largest relative
$L_{2}$ error among the test data set are shown in Fig. 28 for the
displacements and strains.

(a) Displacements

(b) Absolute error for displacements

(c) Strains

(d) Absolute error for strains

Figure 28: Displacements, strains and errors calculated with FNO for the
Burgers data set (load case with highest $L_{2}$ error within the test data
set).

The reason for the poor agreement at the boundaries cannot be attributed
solely to the non-periodic boundary conditions of the input functions,
although for the periodic boundary conditions, the median of the
$\epsilon_{\text{rel}}$ error is lower for the periodic boundary conditions
($2.396\cdot 10^{-2}$) than for the non-periodic ones ($9.278\cdot 10^{-2}$).
Large deviations at the boundaries of the periodic domain are still present
(Fig. 28). These errors can be reduced by regularization (for example by means
of an energy functional), as observed with the Neural FEM.

The average run times per epoch for computations on CPU and GPU are shown in
Table 5. The midrange mobile graphics in use accelerated the calculation by
factor 4 for the 32 Bit floating point accuracy. For a comparable acceleration
for FP64, instead of consumer graphics cards specific High Performance
Computing (HPC) accelerator cards are necessary. We skip a calculation of the
$\epsilon_{\mathrm{rel}}$ errors with double precision because no large
influence on the total error can be expected, according to the results of the
Neural FEM.

| CPU ($d_{v}=64$) | GPU ($d_{v}=64$) | CPU ($d_{v}=12$) | GPU ($d_{v}=12$)
---|---|---|---|---
FNO (FP32) | $9.754\mathrm{\leavevmode\nobreak\ s}$ | $2.447\mathrm{\leavevmode\nobreak\ s}$ | $2.737\mathrm{\leavevmode\nobreak\ s}$ | $0.765\mathrm{\leavevmode\nobreak\ s}$
PINO (FP32) | $9.926\mathrm{\leavevmode\nobreak\ s}$ | $2.491\mathrm{\leavevmode\nobreak\ s}$ | $2.874\mathrm{\leavevmode\nobreak\ s}$ | $0.837\mathrm{\leavevmode\nobreak\ s}$
FNO (FP64) | $18.120\mathrm{\leavevmode\nobreak\ s}$ | $14.694\mathrm{\leavevmode\nobreak\ s}$ | $5.163\mathrm{\leavevmode\nobreak\ s}$ | $4.458\mathrm{\leavevmode\nobreak\ s}$
PINO (FP64) | $17.041\mathrm{\leavevmode\nobreak\ s}$ | $14.763\mathrm{\leavevmode\nobreak\ s}$ | $5.076\mathrm{\leavevmode\nobreak\ s}$ | $4.491\mathrm{\leavevmode\nobreak\ s}$
Table 5: Average runtimes per epoch of FNO and PINO on FP32 and FP64 with
$d_{v}=64$ and $d_{v}=12$.

####  5.6.3 Zero-shot super resolution

A main feature of neural operators is the consistency of the numerical error
over different discretization levels. This is manifested by a nearly constant
progression of the error across the discretization level, as shown in Fig. 29.
Therefore, zero-shot super resolution becomes possible, which means that the
NN can be evaluated on a finer grid than the one that used for training. In
the present work, the NNs are trained on a dataset based on FEM solutions with
a grid with 1024 nodes, but can also be evaluated on finer discretization
levels with the same error. The only exception is the PINO which uses a finite
difference method in the optimization process in the current contribution. The
reference solution for the finer discretization task is obtained from an FEM
analysis with 8192 node points. This data set is cubically interpolated to all
other discretizations for comparison with the results of the NNs.

Figure 29: Relative $L_{2}$ error on a single instance of the test data set with varying discretization used for NN evaluation. Model | Error in $u$ | Error in $u_{,X}$ | Run time on CPU in s | Run time on GPU in s
---|---|---|---|---
FNO | $2.438\cdot 10^{-13}$ | $1.145\cdot 10^{-1}$ |  $81.289$ (estimated) | $20.392$
PINO | $2.688\cdot 10^{-13}$ | $1.382\cdot 10^{-2}$ |  $82.713$ (estimated) | $20.755$
DeepONet | $5.070\cdot 10^{-13}$ | $4.022\cdot 10^{-2}$ | $1.167$ | –
PIDeepONet | $2.783\cdot 10^{-13}$ | $2.800\cdot 10^{-3}$ | $10.117$ | –
Table 6: Comparison of the $\epsilon_{rel}$ errors in the test data set with
FP32.

##  6 Summary and Outlook

In this work, different NN methods have been analyzed and applied to examples
from elastostatics. Specifically, a 1D tensile bar with a hyperelastic
material (Example A, Section 5.1.1) and with a linear elastic material
(Example B, Section 5.1.2) has been investigated. Moreover, a plate made of a
Neo-Hookean material has been analyzed for two load cases, vertical loading
and uniaxial tension (Example C, Section 5.2).

##### Physics Informed Neural Networks (PINN)

In the basic form of classical PINN [7], the empirical risk is built from the
squared residuals of the differential operators. In various works [20, 19,
32], it has been shown that such a PINN is difficult or impossible to train
even for simple examples, so alternative forms of regularization have been
developed. The present work particularly studies the DEM, based on the
principle of minimal potential energy, and the cPINN, based on the game-
theory. The results for Example A using these three approaches (PINN, DEM,
cPINN) are compared in Fig. 30. The average accuracies are relatively similar,
but the comparatively long run time of the cPINN is disadvantageous.

Figure 30: Relative $L_{2}$ errors and run times for Example A with the Neural
FEM approaches.

According to [9], the relative errors can be reduced by up to 2 orders of
magnitude by the cPINN, which could not be demonstrated with Example A in this
work. The training of cPINN and DEM should converge in more cases than pure
PINN, i.e. it is more robust, as demonstrated with Example C. Overall, the
PINN performes best in Example A. However, the latter is not suitable to show
the training pathologies of the PINN. Those pathologies were demonstrated only
on Example C, where the training of the PINN fails, but the DEM can be applied
successfully. The error measured in the $L_{2}$ norm is relatively small. But
the absolute errors of the equivalent stresses for the 2D plate in the
vertical load case show deviations of up to $77\mathrm{\leavevmode\nobreak\
Nm}^{-2}$ at the restraint, for a maximum stress of about
$142\mathrm{\leavevmode\nobreak\ N}\mathrm{m}^{-2}$.

Our results with Example C underline the suggestion from [18] that PINN and
DEM as well as PIDeepONet in their present form are not able to resolve stress
concentrations. Further work is necessary to find and analyze alternative
approaches with improved accuracy and applicability for classical tasks in
solid mechanics like the investigation of critical areas in strength analysis.

###### DEM

In DEM, the convergence order of the integration method does not seem to be
significant after reaching a certain limit of accuracy. However, DEM holds the
risk of overfitting, which must be accounted for by early stopping or a
sufficient number of collocation points (support points). This topic has not
been addressed in the literature up to now.

An advantage from the use of the potential energy is the reduction of the
order of differentiation, which also decreases the numerical effort. The run
time is about a factor of six lower compared to the PINN.

All studies show an intense dependency of the result from the initialization
of the NN parameters. In extreme cases, the optimizations converges towards
different functions. The relative $L_{2}$ error of the DEM with trapezoidal
rule and 1000 collocation points results in the whole range from $4.898\cdot
10^{-6}$ up to $4.277\cdot^{-4}$ for displacements and $3.774\cdot 10^{-5}$ up
to $3.701\cdot 10^{-3}$ for strains. This indicates that the expensive
training has to be conducted several times, until an acceptable ML model is
found. Further improvement possibilities are Glorot initialization of the NN
parameters [31] and pretraining [33].

##### Neural Operator Methods

The Neural Operator methods have been applied to Example B. All models are
calculated with single precision floating point since the investigation of
DeepONet, PIDeepONet and the Neural FEM did not yield significant effects on
the accuracy of the results. Furthermore, the DeepONets have been optimized
with Full-Batch training.

The analysis shows that the DeepONet is significantly faster than the FNO,
even if calculated on a GPU. Both FNO and DeepONet can learn the solution
operator of the parametric PDE, but the achieved accuracies are not sufficient
in many cases.

##### Outlook

According to the results presented, the accuracy of NNs has to be improved for
a reliable engineering application in elastostatics.

One promising approach is a network architecture based on the Neural Attention
mechanism, that is claimed to improve the accuracy of PINNs about up to 2
orders of magnitude [19]. A similar suggestion is made with regard to the
Physics-Augmented Learning [32]. The sequential learning and the curriculum
learning are two further approaches that might improve the learning ability of
PINNs [20].

In future works, second order methods should be investigated for the
optimization of NN parameters during learning [34, 35]. The example of L-BFGS
shows, that the optimization with higher order methods can be immensely
accelerated. They also take profit from larger batches, which makes them more
efficient in terms of data parallelism.

## Abbreviations

  * •

ACGD Adaptive Competitive Gradient Descent

  * •

ANN Artificial Neural Network

  * •

BC Boundary Conditions

  * •

CGD Competitive Gradient Descent

  * •

cPINN competitive Physics Informed Neural Network

  * •

DEM Direct Energy Method

  * •

FCNN Fully Connected Neural Network

  * •

FEM Finite Element Method

  * •

FFT Fast Fourier Transformation

  * •

FNO Fourier Neural Operator

  * •

MC Monte Carlo

  * •

MSE Mean Square Error

  * •

NN Neural Network

  * •

PDE Partial Differential Equation

  * •

PINN Physics Informed Neural Network

  * •

PINO Physics Informed Neural Operator

  * •

PPDE Parametric Partial Differential Equation

  * •

SR Squared Residual

  * •

TF Transfer Learning

The authors cordially thank Mr. Emre Sahin for his contribution to the present
work.

## References

  * [1] Frederic E. Bock, Roland C. Aydin, Christian J. Cyron, Norbert Huber, Surya R. Kalidindi and Benjamin Klusemann  “A Review of the Application of Machine Learning and Data Mining Approaches in Continuum Materials Mechanics”  In _Frontiers in Materials_ 6, 2019  DOI: 10.3389/fmats.2019.00110
  * [2] T. Kirchdoerfer and M. Ortiz  “Data-driven computational mechanics”  In _Computer Methods in Applied Mechanics and Engineering_ 304, 2016, pp. 81–101  DOI: https://doi.org/10.1016/j.cma.2016.02.001
  * [3] Ronak Shoghi and Alexander Hartmaier  “Optimal Data-Generation Strategy for Machine Learning Yield Functions in Anisotropic Plasticity”  In _Frontiers in Materials_ 9, 2022, pp. 868248  DOI: 10.3389/fmats.2022.868248
  * [4] Adam Paszke et al.  “PyTorch: An Imperative Style, High-Performance Deep Learning Library”  In _Advances in Neural Information Processing Systems 32_ Curran Associates, Inc., 2019, pp. 8024–8035  URL: http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf
  * [5] Martín Abadi et al.  “TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems” Software available from tensorflow.org, 2015  URL: https://www.tensorflow.org/
  * [6] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart and Anima Anandkumar  “Fourier Neural Operator for Parametric Partial Differential Equations”  arXiv, 2020b  DOI: 10.48550/ARXIV.2010.08895
  * [7] M. Raissi, P. Perdikaris and G.E. Karniadakis  “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations”  In _Journal of Computational Physics_ 378, 2019, pp. 686–707  DOI: https://doi.org/10.1016/j.jcp.2018.10.045
  * [8] Vien Minh Nguyen-Thanh, Xiaoying Zhuang and Timon Rabczuk  “A deep energy method for finite deformation hyperelasticity”  In _European Journal of Mechanics - A/Solids_ 80, 2020, pp. 103874  DOI: https://doi.org/10.1016/j.euromechsol.2019.103874
  * [9] Qi Zeng, Yash Kothari, Spencer H. Bryngelson and Florian Schäfer  “Competitive Physics Informed Networks”  arXiv, 2022  DOI: 10.48550/ARXIV.2204.11144
  * [10] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang and George Karniadakis  “Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators”  In _Nature Machine Intelligence_ 3, 2021, pp. 218–229  DOI: 10.1038/s42256-021-00302-5
  * [11] Stefan Kollmannsberger  “Deep learning in computational mechanics : an introductory course”, Studies in computational intelligence ;  Cham: Springer, 2021, pp. 108  URL: https://link.springer.com/10.1007/978-3-030-76587-3
  * [12] Sifan Wang, Hanwen Wang and Paris Perdikaris  “Learning the solution operator of parametric partial differential equations with physics-informed DeepOnets”  arXiv, 2021b  DOI: 10.48550/ARXIV.2103.10974
  * [13] Zongyi Li, Hongkai Zheng, Nikola Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzadenesheli and Anima Anandkumar  “Physics-Informed Neural Operator for Learning Partial Differential Equations”  arXiv, 2021  DOI: 10.48550/ARXIV.2111.03794
  * [14] Diab W. Abueidda, Qiyue Lu and Seid Koric  “Meshless physics-informed deep learning method for three-dimensional solid mechanics”  In _International Journal for Numerical Methods in Engineering_ 122.23, 2021, pp. 7182–7201  DOI: https://doi.org/10.1002/nme.6828
  * [15] Kurt Hornik, Maxwell Stinchcombe and Halbert White  “Multilayer feedforward networks are universal approximators”  In _Neural Networks_ 2.5, 1989, pp. 359–366  DOI: https://doi.org/10.1016/0893-6080(89)90020-8
  * [16] Diederik P. Kingma and Jimmy Ba  “Adam: A Method for Stochastic Optimization”  arXiv, 2014  DOI: 10.48550/ARXIV.1412.6980
  * [17] Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zach DeVito, Zeming Lin, Alban Desmaison, Luca Antiga and Adam Lerer  “Automatic differentiation in PyTorch”, 2017
  * [18] Jan N. Fuhg and Nikolaos Bouklas  “The mixed Deep Energy Method for resolving concentration features in finite strain hyperelasticity”  In _Journal of Computational Physics_ 451, 2022, pp. 110839  DOI: https://doi.org/10.1016/j.jcp.2021.110839
  * [19] Sifan Wang, Yujun Teng and Paris Perdikaris  “Understanding and Mitigating Gradient Flow Pathologies in Physics-Informed Neural Networks”  In _SIAM Journal on Scientific Computing_ 43.5, 2021a, pp. A3055–A3081  DOI: 10.1137/20M1318043
  * [20] Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby and Michael W Mahoney  “Characterizing possible failure modes in physics-informed neural networks”  In _Advances in Neural Information Processing Systems_ 34 Curran Associates, Inc., 2021, pp. 26548–26560  URL: https://proceedings.neurips.cc/paper/2021/file/df438e5206f31600e6ae4af72f2725f1-Paper.pdf
  * [21] Zhangyong Liang, Huanhuan Gao and Tingting Li  “SEM: a shallow energy method for finite deformation hyperelasticity problems”  In _Acta Mechanica_ 233, 2022  DOI: 10.1007/s00707-022-03174-x
  * [22] Kaushik Bhattacharya, Bamdad Hosseini, Nikola B. Kovachki and Andrew M. Stuart  “Model Reduction And Neural Networks For Parametric PDEs”  In _The SMAI Journal of computational mathematics_ 7 Société de Mathématiques Appliquées et Industrielles, 2021, pp. 121–157  DOI: 10.5802/smai-jcm.74
  * [23] Yinhao Zhu and Nicholas Zabaras  “Bayesian deep convolutional encoder–decoder networks for surrogate modeling and uncertainty quantification”  In _Journal of Computational Physics_ 366, 2018, pp. 415–447  DOI: https://doi.org/10.1016/j.jcp.2018.04.018
  * [24] Weinan E and Bing Yu  “The Deep Ritz Method: A Deep Learning-Based Numerical Algorithm for Solving Variational Problems”  In _Communications in Mathematics and Statistics_ , 2018, pp. 1–12  DOI: https://doi.org/10.1007/s40304-018-0127-z
  * [25] Sharmila Karumuri, Rohit Tripathy, Ilias Bilionis and Jitesh Panchal  “Simulator-free solution of high-dimensional stochastic elliptic partial differential equations using deep neural networks”  In _Journal of Computational Physics_ 404, 2020, pp. 109120  DOI: https://doi.org/10.1016/j.jcp.2019.109120
  * [26] Florian Schäfer, Hongkai Zheng and Anima Anandkumar  “Implicit competitive regularization in GANs”  In _CoRR_ abs/1910.05852, 2019  arXiv: http://arxiv.org/abs/1910.05852
  * [27] Hongkai Zheng  “CGDs”  In _GitHub repository_ GitHub, https://github.com/devzhk/cgds-package, 2020
  * [28] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang and George Em Karniadakis  “Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators”  In _Nature Machine Intelligence_ 3.3, 2021, pp. 218–229  DOI: https://doi.org/10.1038/s42256-021-00302-5
  * [29] Sifan Wang, Hanwen Wang and Paris Perdikaris  “Learning the solution operator of parametric partial differential equations with physics-informed DeepONets”  In _Science Advances_ 7.40, 2021, pp. eabi8605  DOI: 10.1126/sciadv.abi8605
  * [30] Mohammad Amin Nabian and Hadi Meidani  “A deep learning solution approach for high-dimensional random differential equations”  In _Probabilistic Engineering Mechanics_ 57, 2019  DOI: 10.1016/j.probengmech.2019.05.001
  * [31] Xavier Glorot and Yoshua Bengio  “Understanding the difficulty of training deep feedforward neural networks”  In _International Conference on Artificial Intelligence and Statistics_ , 2010
  * [32] Ziming Liu, Yunyue Chen, Yuanqi Du and Max Tegmark  “Physics-Augmented Learning: A New Paradigm Beyond Physics-Informed Learning”  In _CoRR_ abs/2109.13901, 2021  arXiv: https://arxiv.org/abs/2109.13901
  * [33] Genki Yagawa and Atsuya Oishi  “Computational Mechanics with Neural Networks”, 2021  DOI: 10.1007/978-3-030-66111-3
  * [34] James Martens and Ilya Sutskever  “Training Deep and Recurrent Networks with Hessian-Free Optimization”  In _Neural Networks: Tricks of the Trade: Second Edition_ Berlin, Heidelberg: Springer Berlin Heidelberg, 2012, pp. 479–535  DOI: 10.1007/978-3-642-35289-8_27
  * [35] Zhewei Yao, Amir Gholami, Kurt Keutzer and Michael W. Mahoney  “PyHessian: Neural Networks Through the Lens of the Hessian”  In _2020 IEEE International Conference on Big Data (Big Data)_ , 2019, pp. 581–590