# The mixed Deep Energy Method for resolving concentration features in finite
strain hyperelasticity

Jan Niklas Fuhg
Sibley School of Mechanical and Aerospace Engineering
Cornell University, New York, USA
<EMAIL_ADDRESS>
&Nikolaos Bouklas
Sibley School of Mechanical and Aerospace Engineering
Cornell University, New York, USA
<EMAIL_ADDRESS>

###### Abstract

The introduction of Physics-informed Neural Networks (PINNs) has led to an
increased interest in deep neural networks as universal approximators of PDEs
in the solid mechanics community. Recently, the Deep Energy Method (DEM) has
been proposed. DEM is based on energy minimization principles, contrary to
PINN which is based on the residual of the PDEs. A significant advantage of
DEM, is that it requires the approximation of lower order derivatives compared
to formulations that are based on strong form residuals. However both DEM and
classical PINN formulations struggle to resolve fine features of the stress
and displacement fields, for example concentration features in solid mechanics
applications. We propose an extension to the Deep Energy Method (DEM) to
resolve these features for finite strain hyperelasticity. The developed
framework termed mixed Deep Energy Method (mDEM) introduces stress measures as
an additional output of the NN to the recently introduced pure displacement
formulation. Using this approach, Neumann boundary conditions are approximated
more accurately and the accuracy around spatial features which are typically
responsible for high concentrations is increased. In order to make the
proposed approach more versatile, we introduce a numerical integration scheme
based on Delaunay integration, which enables the mDEM framework to be used for
random training point position sets commonly needed for computational domains
with stress concentrations, i.e. domains with holes, notches, etc. We
highlight the advantages of the proposed approach while showing the
shortcomings of classical PINN and DEM formulations. The method is offering
comparable results to Finite-Element Method (FEM) on the forward calculation
of challenging computational experiments involving domains with fine geometric
features and concentrated loads, but additionally offers unique capabilities
for the solution of inverse problems and parameter estimation in the context
of hyperelasticity.

_Keywords_ Physics-informed neural networks $\cdot$ Deep energy method $\cdot$
Mixed deep energy method $\cdot$ Solid mechanics $\cdot$ Hyperelasticity
$\cdot$ Finite Strain

##  1 Introduction

Solid mechanics studies the motion and deformation of solid materials when
subjected to external forces. In the last years we have seen a rapid growth in
the utilization of machine learning and data-driven techniques in
computational solid mechanics. For example, Kirchdoerfer and Ortiz (2016)
defined a computational data-driven approach that enforces constraints and
conservation laws directly on experimental data. This data-driven approach has
been extended to a wide range of directions such as inelastic problems
(Eggersmann et al., 2019) and also extended to a non-local setting (González
et al., 2019). Machine learning has also been used to establish data-driven
elastic and inelastic constitutive laws (Ibañez et al., 2017; Huang et al.,
2020; Fuhg et al., 2021). These ideas make use of existing infrastructures
such as existing Finite-Element Method (FEM) solvers to solve the accompanying
partial differential equations (PDEs). Additionally, machine learning
techniques have been employed for the development of intrusive and non-
intrusive Reduced Order Modeling (ROM) schemes for the accelerated solutions
of PDEs (Kadeethum et al., 2021; Hernandez et al., 2021) following the
offline-online paradigm, in which synthetic training data is often constructed
from FEM solvers.

As an alternative approach, Neural Networks (NN) have recently gained a lot of
attention in the solid mechanics community, and more broadly in the physic-
based modeling community, for finding solutions of boundary value problems
from scratch, due to their abilities as universal approximators (Hornik et
al., 1989). While the idea is not new and has already been proposed over 20
years ago (Lagaris et al., 1998), it has seen remarkable success recently due
to advances in GPU hardware and neural network technologies ( e.g. deep
learning libraries, stochastic gradient descent algorithms, etc). In this
context Physics-informed Neural Networks (PINNS), first introduced in Raissi
et al. (2019), have emerged as a deep learning algorithm for solving PDEs,
based on automatic differentiation. The method has several advantages compared
to classical numerical techniques with the major ones being that the algorithm
is meshfree and is significantly easier to implement than traditional methods
such as FEM while at the same time allows for the incorporation of constraints
due to experimental observations and solution of inverse problems for
parameter estimation. As PINN is an optimization-based solution technique, it
also has the potential to be utilized in highly nonlinear problems where
traditional solution techniques are not sufficient or necessarily robust. PINN
has been extended into different directions such as variational problems
(VPINN) (Kharazmi et al., 2021) or a subdomain-version (XPINN) which defines
local neural networks in order to allow parallelization of the computation and
also to resolve the solution in heterogeneous domains (Jagtap and Karniadakis,
2020).

From these ideas different approaches have been developed to utilize PINNs for
applications in solid mechanics. Haghighat et al. (2020) trained PINN for
small strain elasticity and elastoplasticity problems, Kadeethum et al. (2020)
used the approach to solve the coupled Biot’s equations, and Rao et al. (2020)
explored a mixed formulation for regular physics-informed neural networks with
applications to small strain elastodynamics problems. Abueidda et al. (2020)
utilized the residual formulation of PINNs to study linear elasticity,
hyperelasticity, and plasticity. Recently, the Deep Energy Method (DEM) was
introduced, which defines the loss function not based on residual equations
but based on energy minimization principles (Weinan and Yu, 2018; Samaniego et
al., 2020; Nguyen-Thanh et al., 2020). Compared to classical PINN, this method
has advantages because it only requires first order differentiation through
the neural network. On the other hand it relies on accurate numerical
integration techniques. With this technique Nguyen-Thanh et al. (2020) solve
two and three dimensional finite-strain hyperelastic problems. However as will
be shown in this work, both the PINN and DEM struggle to resolve displacement
and stress concentrations. In this context we propose the mixed Deep Energy
Method (mDEM) which, by requiring both displacements and stress outputs from
the NN architecture, is able to resolve concentration features on comparable
levels to FEM. Hence, the presented method aims to improve the forward
capabilities of the physics-informed machine learning paradigm for solid
mechanics applications as a first step towards efficient and accurate inverse
problem calculations. Furthermore we provide some extension to the initial DEM
formulation by offering a more reliable integration scheme based on Delaunay
triangulation.

In section 2 the general finite-strain solid mechanics framework for
hyperelasticity is introduced, and specialized to a Neo-Hookean type strain
energy. A short introduction to neural networks, PINNs and DEMs specialized
for 2D hyperelasticity is provided in Section 3. The newly proposed mixed Deep
Energy Method formulation as well as an accompanying integration scheme based
on Delaunay triangulation are introduced in Section 4. In Section 5 we test
the approach on various computational experiments comparing against existing
PINN- and DEM-based implementations. The paper is concluded in Section 6.

##  2 Finite-strain hyperelasticity

Let a bounded domain be occupied by an elastic body
$\mathcal{B}\subset\mathbb{R}^{3}$. The body’s boundary $\Gamma$ consists of
two non-overlapping regions $\Gamma_{t}$ and $\Gamma_{u}$ such that
$\Gamma=\Gamma_{t}\cup\Gamma_{u}$ (Fig. 1).

Figure 1: Solid domain with boundary conditions

Let the referential and current position vectors of the body be denoted by
$\bm{X}$ and $\bm{x}$ respectively. The time-dependent motion between the two
positions is then given by

$\bm{x}=\bm{\varphi}(\bm{X},t)=\bm{X}+\bm{u}(\bm{X},t)$ (1)

with $\bm{\varphi}(\bm{X},t)$ denoting the motion of the body and $\bm{u}$ the
time-dependent displacement field. The gradient of the motion with respect to
the initial position yield the deformation gradient

$\bm{F}=\text{Grad}\bm{\varphi}(\bm{X}).$ (2)

The governing equation for quasi-static problems in solid continua is the
equilibrium equation, which in strong form reads

$\displaystyle\text{Div}\bm{P}+\bm{f}_{b}$
$\displaystyle=\bm{0},\quad\mathrm{in}\,\,\mathcal{B}$ (3)

with

$\displaystyle\bm{u}$ $\displaystyle=\tilde{\bm{u}}\quad$
$\displaystyle\text{on}\,\,\Gamma_{u},$ (4) $\displaystyle\bm{P}\cdot\bm{N}$
$\displaystyle=\tilde{\bm{t}}\quad$ $\displaystyle\text{on}\,\,\Gamma_{t},$

where $\bm{P}$ and $\bm{f}_{b}$ denote the first Piola-Kirchhoff stress tensor
and the body force vector respectively. The body is in general subjected to
two types of boundary conditions with the outward normal $\bm{N}$, the
prescribed displacement $\tilde{\bm{u}}$ on $\Gamma_{t}$ and the prescribed
surface traction $\tilde{\bm{t}}$ on and $\Gamma_{u}$ as denoted in eq. 4.

In the hyperelastic framework, the existence of the strain energy function
$\Psi$ is postulated, which enables the computation of the first Piola-
Kirchhoff stress tensor by

$\bm{P}=\frac{\partial\Psi}{\partial\bm{F}}.$ (5)

In this paper we define a Neo-Hookean type strain energy which allows for
compressibility as

$\Psi=\frac{1}{4}\lambda(\log J^{2}-1-2\log
J)+\frac{1}{2}\mu(\text{tr}(\bm{C})-2-2\log J),$ (6)

where $\lambda$ and $\mu$ are the Lamé constants and $J$ and $\bm{C}$ are
given by $J=\det\bm{F}$ and $\bm{C}=\bm{F}^{T}\bm{F}$. Considering only
conservative loads and neglecting inertial contributions, the potential energy
of the body is defined by

$\Pi(\bm{\varphi})=\int_{\mathcal{B}}\Psi\,dV-\int_{\mathcal{B}}\bm{f}_{b}\cdot\bm{\varphi}\,dV-\int_{\Gamma_{t}}\tilde{\bm{t}}\cdot\bm{\varphi}\,dA.$
(7)

The energy functional from eq. 7 can be minimized with regards to
$\bm{\varphi}$ which yields the deformation that fulfills static equilibrium.

##  3 Deep neural networks, Physics-informed Neural Networks and the Deep
Energy Method

In this section we provide a brief overview of the neural networks
formulation, Physics-informed Neural Networks and the Deep Energy Method.
Neural networks are in general composed of one input, one output and $n_{D}-1$
hidden layers. Let the weights and biases of the $k^{\text{th}}$ layer be
denoted by $\bm{W}^{k}$ and $\bm{b}^{k}$. Consider that the $k^{\text{th}}$
hidden layer transfers some output $\bm{x}^{k}$ to the $(k+1)^{\text{th}}$
layers which applies an affine transformation

$\mathcal{L}(\bm{x}^{k})=\bm{W}^{k+1}\bm{x}^{k}+\bm{b}^{k+1}.$ (8)

and some activation activation function $\sigma$ to it. We consider the
activation function of the last layer to be identity. Since, equation (8) is
applied in every layer of the network an input $\bm{x}$ yields a network
output of the form

$\hat{\bm{y}}(\bm{x})=(\mathcal{L}_{k}\circ\sigma\circ\mathcal{L}_{k-1}\circ\cdots\circ\sigma\circ\mathcal{L}_{1})(\bm{x})$
(9)

where $\circ$ is a composition operator. The goal of neural networks is to
find the optimal value set of trainable parameters
$\bm{\Theta}=\\{\bm{W}^{k},\bm{b}^{k}\\}_{k=1}^{n_{D}}$ such that the networks
provides the best fit for the input-output mapping. This is achieved by
following an optimization procedure defined over some loss function
$J(\bm{\Theta})$

$\bm{\Theta}^{\star}=\operatorname*{arg\,min}_{\bm{\Theta}}J(\bm{\Theta}).$
(10)

This problem is typically solved in an iterative manner by employing a
stochastic gradient-descent approach. A review on neural network optimization
methods is provided in Bottou et al. (2018). For more information on neural
networks we refer to Goodfellow et al. (2016).
PINNs define a global shape function for the field(s) of interest over the
computational domain and then require the stresses derived with automatic
differentiation to fulfill a residual based on the strong form, including the
equilibrium equation (3) as well as the corresponding boundary conditions (4).
The computational domain is spanned by a set of collocation points. Here we
will specialize the overview of PINN in the context of hyperelasticity. In
this case the unkown field that will be approximated through PINN is the
displacement field. In absence of body forces, the residual is of the form

$\bm{\mathcal{R}}=\text{Div}\bm{P}=\bm{0}.$ (11)

The collocation points hereby act as inputs to the neural network which yield
the outputs $\bm{z}(\bm{X},\bm{\Theta})$. The outputs
$\bm{z}(\bm{X},\bm{\Theta})$ can either be subjected to

$\hat{\bm{u}}(\bm{X},\bm{\Theta})=\bm{A}(\bm{X})+\bm{B}(\bm{X})\circ\bm{z}(\bm{X},\bm{\Theta})$
(12)

where $\bm{A}$ and $\bm{B}$ are a-priori chosen in a way such that the
displacement prediction $\hat{\bm{u}}$ fulfills the displacement boundary
conditions acting on the body. Here, $\bm{A}$ and $\bm{B}$ can either chosen
analytically or can be neural networks themselves. For a detailed discussion
refer to Rao et al. (2020). Alternately, we can choose

$\hat{\bm{u}}(\bm{X},\bm{\Theta})=\bm{z}(\bm{X},\bm{\Theta})$ (13)

where the boundary conditions need to be trained. From the displacements the
deformation gradient can be obtained by

$\displaystyle\hat{\bm{F}}(\bm{X},\bm{\Theta})$
$\displaystyle=\bm{I}+\text{Grad}\hat{\bm{u}}(\bm{X},\bm{\Theta}).$ (14)

Consequently, the stresses $\bm{P}(\hat{\bm{F}}(\bm{X},\bm{\Theta}))$ can be
computed following their definitions eq. 5.

Consider a set of domain training points
$\\{\bm{X}_{\mathcal{R}}^{i}\\}_{i=1}^{N_{\mathcal{R}}}$, displacement
boundary training points $\\{\bm{X}_{u}^{i}\\}_{i=1}^{N_{u}}$ and traction
boundary training points $\\{\bm{X}_{t}^{i}\\}_{i=1}^{N_{t}}$ where
$N_{\mathcal{R}}$, $N_{u}$ and $N_{t}$ denote the number of points,
respectively. The PINN loss function is then defined as

$\displaystyle\bm{\Theta}^{\star}$
$\displaystyle=\operatorname*{arg\,min}_{\bm{\Theta}}\qquad
W_{\mathcal{R}}\text{MSE}_{\mathcal{R}}\left(\\{\bm{X}_{\mathcal{R}}^{i}\\}_{i=1}^{N_{\mathcal{R}}},\bm{\Theta}\right)+W_{t}\text{MSE}_{t}\left(\\{\bm{X}_{t}^{i}\\}_{i=1}^{N_{t}},\bm{\Theta}\right)+\underbrace{W_{u}\text{MSE}_{u}\left(\\{\bm{X}_{u}^{i}\\}_{i=1}^{N_{u}},\bm{\Theta}\right)}_{\text{If
not a-priori fulfilled}}$ (15)

with some weight values $W_{\mathcal{R}}$, $W_{t}$ and $W_{u}$ and
specializing to a 2D case

$\displaystyle\text{MSE}_{\mathcal{R}}\left(\\{\bm{X}_{\mathcal{R}}^{i}\\}_{i=1}^{N_{\mathcal{R}}},\bm{\Theta}\right)$
$\displaystyle=\frac{1}{N_{\mathcal{R}}}\sum_{i=1}^{N_{\mathcal{R}}}\absolutevalue{\frac{\partial}{\partial
X_{1}}P_{11}(\hat{\bm{F}}(\bm{X}_{\mathcal{R}}^{i},\bm{\Theta}))+\frac{\partial}{\partial
X_{2}}P_{12}(\hat{\bm{F}}(\bm{X}_{\mathcal{R}}^{i},\bm{\Theta}))}^{2}$ (16)
$\displaystyle+\frac{1}{N_{\mathcal{R}}}\sum_{i=1}^{N_{\mathcal{R}}}\absolutevalue{\frac{\partial}{\partial
X_{1}}P_{21}(\hat{\bm{F}}(\bm{X}_{\mathcal{R}}^{i},\bm{\Theta}))+\frac{\partial}{\partial
X_{2}}{P}_{22}(\hat{\bm{F}}(\bm{X}_{\mathcal{R}}^{i},\bm{\Theta}))}^{2},$
$\displaystyle\text{MSE}_{t}\left(\\{\bm{X}_{t}^{i}\\}_{i=1}^{N_{t}},\bm{\Theta}\right)$
$\displaystyle=\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}\absolutevalue{\bm{P}(\hat{\bm{F}}(\bm{X}_{t}^{i},\bm{\Theta}))\bm{N}_{i}-\overline{\bm{t}}(\bm{X}_{t}^{i})}^{2},$
$\displaystyle\text{MSE}_{u}\left(\\{\bm{X}_{u}^{i}\\}_{i=1}^{N_{u}},\bm{\Theta}\right)$
$\displaystyle=\frac{1}{N_{u}}\sum_{i=1}^{N_{u}}\absolutevalue{\hat{\bm{u}}(\bm{X}_{u}^{i},\bm{\Theta})-\bm{u}(\bm{X}_{u}^{i})}^{2}.$

The PINN concept is summarized in Figure 2.

Figure 2: PINN overview

Maintaining some of the architecture of PINNs, Nguyen-Thanh et al. (2020)
proposed an alternative approach. The neural network formulation is still used
to solve quasi-static equilibrium problems by utilizing the network as a
global shape function for the displacement over the body. However instead of
defining the residual of the governing PDE as the loss function, the authors
define the loss as the potential energy of the solid body, as defined in
equation 7. By following that, the minimization of the loss function over some
collocation points that span the computational domain and its boundary,
through stationarity of the potential energy in nonlinear problems, leads to
an approximate solution. The advantage of this approach is that only first
order automatic differentiation is needed for the solution procedure, instead
of second order as required in PINN. This, can potentially lead to quicker
convergence and higher accuracy. As a downside the method requires an
effective integration procedure over the domain spanned by the collocation
points. In essence DEM uses the same concept as PINNs for the input and the
output to the neural network, i.e. collocation point positions as input and
the output is defined by either eq. 12 or 13. In addition to the deformation
gradient (eq. 14), DEM requires a formulation for the motion which is given by

$\displaystyle\hat{\bm{\varphi}}(\bm{X},\bm{\Theta})$
$\displaystyle=\bm{X}+\hat{\bm{u}}(\bm{X},\bm{\Theta}).$ (17)

Consider a set of domain training points
$\\{\bm{X}_{\Pi}^{i}\\}_{i=1}^{N_{\Pi}}$ and boundary training points
$\\{\bm{X}_{u}^{i}\\}_{i=1}^{N_{u}}$ where $N_{\Pi}$ and $N_{u}$ denote the
number of points, respectively. Hence, the loss function can be written as

$\displaystyle\bm{\Theta}^{\star}$
$\displaystyle=\operatorname*{arg\,min}_{\bm{\Theta}}\qquad$
$\displaystyle\Pi(\hat{\bm{\varphi}}(\\{\bm{X}_{\Pi}^{i}\\}_{i=1}^{N_{\Pi}},\bm{\Theta}))+\underbrace{W_{u}\text{MSE}_{u}\left(\\{\bm{X}_{u}^{i}\\}_{i=1}^{N_{u}},\bm{\Theta}\right)}_{\text{If
not a-priori fulfilled}}$ (18)

with a weight value $W_{u}$ and where

$\displaystyle\Pi(\hat{\bm{\varphi}}(\\{\bm{X}_{\Pi}^{i}\\}_{i=1}^{N_{\Pi}},\bm{\Theta}))$
$\displaystyle=\int_{\mathcal{B}}\Psi(\hat{\bm{F}}(\\{\bm{X}_{\Pi}^{i}\\}_{i=1}^{N_{\Pi}},\bm{\Theta}))\,dV-\int_{\mathcal{B}}\bm{f}_{b}\cdot\hat{\bm{\varphi}}(\\{\bm{X}_{\Pi}^{i}\\}_{i=1}^{N_{\Pi}},\bm{\Theta})\,dV-\int_{\Gamma_{t}}\tilde{\bm{t}}\cdot\hat{\bm{\varphi}}(\\{\bm{X}_{\Pi}^{i}\\}_{i=1}^{N_{\Pi}},\bm{\Theta})\,dA$
(19)
$\displaystyle\text{MSE}_{u}\left(\\{\bm{X}_{u}^{i}\\}_{i=1}^{N_{u}},\bm{\Theta}\right)$
$\displaystyle=\frac{1}{N_{u}}\sum_{i=1}^{N_{u}}\absolutevalue{\hat{\bm{u}}(\bm{X}_{u}^{i},\bm{\Theta})-\bm{u}(\bm{X}_{u}^{i})}^{2}.$

Figure 3: DEM overview

The process is summarized in Figure 3.

##  4 Mixed Deep Energy Method (MDEM)

For solid mechanics and multiphysics applications, problematic for the DEM and
PINN formulations is that they average out stress and displacement
concentrations as will be shown in Section 5. In order to circumvent these
issues we propose the Mixed Deep Energy Method (mDEM) which introduces
stresses as additional collocation point outputs, which in turn have to
fulfill the traction boundary conditions. Furthermore to reduce the
computation time, training a different neural network for each individual
output, as e.g. proposed in Haghighat et al. (2020), has been explicitly
avoided in this work. Hence, in a 2D setting, in addition to the two outputs
$\bm{z}$ from DEM we also obtain four more outputs $\bm{Z}$ corresponding to
the first Piola-Kirchhoff stress components. Equivalently these outputs can
individually be trained to directly fulfill the corresponding Dirichlet and
Neumann boundary conditions respectively, or can be subjected to

$\displaystyle\hat{\bm{u}}(\bm{X},\bm{\Theta})$
$\displaystyle=\bm{A}(\bm{X})+\bm{B}(\bm{X})\circ\bm{z}(\bm{X},\bm{\Theta})$
(20) $\displaystyle\hat{\bm{P}}(\bm{X},\bm{\Theta})$
$\displaystyle=\bm{C}(\bm{X})+\bm{D}(\bm{X})\circ\bm{Z}(\bm{X},\bm{\Theta})$

such that these conditions are fulfilled a-priori. In order to connect the
displacements and stresses, mDEM requires training of the constitutive
behavior by enforcing

$\hat{\bm{P}}(\bm{X},\bm{\Theta})=\bm{P}(\hat{\bm{F}}(\bm{X},\bm{\Theta}))$
(21)

where the latter term is derived from the displacement output and the
hyperleastic constitutive behavior of equation (5).

Therefore we consider a set of domain training points
$\\{\bm{X}_{\Pi}^{i}\\}_{i=1}^{N_{\Pi}}$, displacement boundary condition
points $\\{\bm{X}_{u}^{i}\\}_{i=1}^{N_{u}}$, and traction boundary condition
points $\\{\bm{X}_{t}^{i}\\}_{i=1}^{N_{t}}$, where $N_{\Pi}$, $N_{u}$ and
$N_{t}$ denote the number of points. The loss function can the be written as

$\displaystyle\bm{\Theta}^{\star}$
$\displaystyle=\operatorname*{arg\,min}_{\bm{\Theta}}\qquad$
$\displaystyle\Pi(\hat{\bm{\varphi}}(\\{\bm{X}_{\Pi}^{i}\\}_{i=1}^{N_{\Pi}},\bm{\Theta}))+W_{P}\text{MSE}_{P}\left(\\{\bm{X}_{\Pi}^{i}\\}_{i=1}^{N_{\Pi}},\bm{\Theta}\right)$
(22)
$\displaystyle+\underbrace{W_{u}\text{MSE}_{u}\left(\\{\bm{X}_{u}^{i}\\}_{i=1}^{N_{u}},\bm{\Theta}\right)}_{\text{If
not a-priori
fulfilled}}+\underbrace{W_{t}\text{MSE}_{t}\left(\\{\bm{X}_{t}^{i}\\}_{i=1}^{N_{t}},\bm{\Theta}\right)}_{\text{If
not a-priori fulfilled}}$

with weight values $W_{u}$,$W_{t}$, $W_{P}$ and where

$\displaystyle\Pi(\hat{\bm{\varphi}}(\\{\bm{X}_{\Pi}^{i}\\}_{i=1}^{N_{\Pi}},\bm{\Theta}))$
$\displaystyle=\int_{\mathcal{B}}\Psi(\hat{\bm{F}}(\\{\bm{X}_{\Pi}^{i}\\}_{i=1}^{N_{\Pi}},\bm{\Theta}))\,dV-\int_{\mathcal{B}}\bm{f}_{b}\cdot\hat{\bm{\varphi}}(\\{\bm{X}_{\Pi}^{i}\\}_{i=1}^{N_{\Pi}},\bm{\Theta})\,dV-\int_{\Gamma_{t}}\tilde{\bm{t}}\cdot\hat{\bm{\varphi}}(\\{\bm{X}_{\Pi}^{i}\\}_{i=1}^{N_{\Pi}},\bm{\Theta})\,dA,$
(23)
$\displaystyle\text{MSE}_{P}\left(\\{\bm{X}_{\Pi}^{i}\\}_{i=1}^{N_{\Pi}},\bm{\Theta}\right)$
$\displaystyle=\frac{1}{N_{\Pi}}\sum_{i=1}^{N_{\Pi}}\absolutevalue{\hat{\bm{P}}(\bm{X}_{\Pi}^{i},\bm{\Theta})-\bm{P}(\hat{\bm{F}}(\bm{X}_{\Pi}^{i},\bm{\Theta}))}^{2},$
$\displaystyle\text{MSE}_{u}\left(\\{\bm{X}_{u}^{i}\\}_{i=1}^{N_{u}},\bm{\Theta}\right)$
$\displaystyle=\frac{1}{N_{u}}\sum_{i=1}^{N_{u}}\absolutevalue{\hat{\bm{u}}(\bm{X}_{u}^{i},\bm{\Theta})-\bm{u}(\bm{X}_{u}^{i})}^{2},$
$\displaystyle\text{MSE}_{t}\left(\\{\bm{X}_{t}^{i}\\}_{i=1}^{N_{t}},\bm{\Theta}\right)$
$\displaystyle=\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}\absolutevalue{\hat{\bm{P}}(\bm{X}_{t}^{i},\bm{\Theta})\bm{N}_{i}-\overline{\bm{t}}(\bm{X}_{t}^{i})}^{2}.$

The process is visualized in Figure 4.

Figure 4: DEM overview

###  4.1 Delauney integration

In the DEM paper (Nguyen-Thanh et al., 2020), the authors explore different
numerical integration techniques for the evaluation of the integral including
Trapezoidal and Simpson’s rule as well as a Monte Carlo integration version.
Even though the presented techniques prove to be sufficient for the
investigated numerical examples, they are difficult to employ when the sample
points need to represent solid continua with fine geometric features, which
commonly lead to stress concentrations. Hence, in this paper we employ a form
of Delaunay-integration to resolve this issue. The computational domain is
tesselated based on its respective sample positions resulting in $n_{e}$
triangles or tetrahedra. Given a triangle $j$ in two dimensions, with the
three points denoted by $P_{i}={x_{i},y_{i}}$ and a corresponding function
value $f_{i}$ with $i=1,2,3$ for the field of interest given at each point.
Then define $J_{j}$ as the area of the triangle and let

$\overline{f_{\bigtriangleup_{j}}}=\frac{f_{1}+f_{2}+f_{3}}{3}$ (24)

denote the mean function value of an element. Following eq. 7, the potential
energy of the domain is then approximated by

$\displaystyle\Pi(\bm{\varphi})$
$\displaystyle=\int_{\mathcal{B}}\Psi\,dV-\int_{\mathcal{B}}\bm{f}_{b}\cdot\bm{\varphi}\,dV-\int_{\Gamma_{t}}\tilde{\bm{t}}\cdot\bm{\varphi}\,dA$
(25)
$\displaystyle=\sum_{i=1}^{n_{e}}\int_{\bigtriangleup_{i}}\Psi_{\bigtriangleup_{i}}d\bigtriangleup_{i}-\sum_{i=1}^{n_{e}}\int_{\bigtriangleup_{i}}\left[\bm{f}_{b}\cdot\bm{\varphi}\right]_{\bigtriangleup_{i}}d\bigtriangleup_{i}-\int_{\Gamma_{t}}\tilde{\bm{t}}\cdot\bm{\varphi}\,dA$
$\displaystyle=\sum_{i=1}^{n_{e}}J_{i}\overline{\Psi_{\bigtriangleup_{i}}}-\sum_{i=1}^{n_{e}}J_{i}\overline{\left[\bm{f}_{b}\cdot\bm{\varphi}\right]_{\bigtriangleup_{i}}}-\sum_{i=1}^{N_{b}}w_{i}\left[\tilde{\bm{t}}_{i}\cdot\bm{\varphi}_{i}\right]$

where $w_{i}$ represents some weights to a common numerical integration
technique such as Simpson’s rule. An example of a constraint Voronoi mesh for
a plate with a "C"-shaped hole (red zone) is displayed in Figure 5.

(a) a

(b) b

Figure 5: Constraint Voronoi tesselation for integration of defected area.

##  5 Numerical examples

In this section we compare the solutions of mDEM to DEM and PINN as well as to
a Finite-Element Method (FEM) solution which will be taken as ground truth.
The FEM results were obtained using the Fenics framework (Alnæs et al., 2015).
The deep learning formulations were implemented111Codes will be made public
after acceptance of this paper. in Pytorch (Paszke et al., 2019) and the
network parameters were optimized using the Adam optimizer (Kingma and Ba,
2014) and the limited Broyden–Fletcher–Goldfarb–Shanno (LFBGS) algorithm
consecutively. In the following, we use analytical formulations to fulfill the
displacement boundary conditions of the applications a-priori, while the
traction boundaries are learned from data. All neural networks are composed of
6 hidden layers with 60 neurons each. Furthermore all networks are trained
with an Adam and LFBGS learning rate of $1e-3$.

###  5.1 Uniaxial loading

In a first example we study the displacement and stress results from PINN,
DEM, mDEM and PINN for a simple uniaxial loading problem as seen in Figure
6(a). A $1\text{m}\times 1\text{m}$ block is subjected to a uniform line load
on the right-hand side while the movement on the left-hand side is restricted.
The FEM solution is generated from a $100\times 100$ elements uniform mesh.
The displacement boundary conditions of PINN, DEM and mDEM are a-priori
fulfilled by setting

$\hat{\bm{u}}(\bm{X},\bm{\theta})=\bm{X}\circ\bm{z}(\bm{X},\bm{\theta}).$ (26)

All the methods are trained using a grid of $200\times 200$ training points
over the computational domain (see Figure 6(b)). The boundary traction
boundary conditions of PINN and mDEM are trained with an additional 5000
training points per edge. Even though the loading is uniaxial, as the left
side of the domain is clamped we are expecting stress concentrations around
the top and bottom of that side.

(a)

(b)

Figure 6: Uniaxial loading problem and training point positions.

After training the networks, the displacement magnitudes of DEM, PINN and mDEM
follow the expected pattern providing a relatively good match to the FEM
result, as displayed in Figure 7. However the $P_{11}$ and $P_{12}$ stress
components approximated by the methods show significant differences, see
Figure 8. It is clear that mDEM is able to resolve the stresses that arise
near the clamping whereas DEM and PINN fail in doing so.

(a) DEM displacement magnitude

(b) PINN displacement magnitude

(c) mDEM displacement magnitude

(d) FEM displacement magnitude

Figure 7: Uniaxial loading problem displacement magnitude.

(a) DEM $P_{11}$

(b) PINN $P_{11}$

(c) mDEM $P_{11}$

(d) FEM $P_{11}$

(e) DEM $P_{12}$

(f) PINN $P_{12}$

(g) mDEM $P_{12}$

(h) FEM $P_{12}$

Figure 8: Uniaxial loading problem. $P_{11}$ (a-d) and $P_{12}$ (e-h) stress
components.

###  5.2 Localized traction boundary

The first example considers the $1m\times 1m$ block as depicted in Figure 9. A
load is applied locally on the right-hand side of the block. The same number
of domain and boundary training points are used as in the previous example
(see Section 5.1 and Figure 6(b)). The displacement boundary conditions of
PINN, DEM and mDEM are a-priori fulfilled by setting

$\hat{\bm{u}}(\bm{X},\bm{\theta})=\bm{X}\circ\bm{z}(\bm{X},\bm{\theta}).$ (27)

Figure 10 displays the resulting displacement magnitudes of the four different
numerical techniques. From the FEM solution it can be seen that there is a
displacement concentration around the locally applied traction boundary
condition. Comparing DEM, PINN and mDEM results shows that only mDEM is able
to resolve this feature accurately. Figure 12(a) shows the error evolution
over the training iterations for PINN. It can be seen that even though the
error converges sufficiently, PINN is not able to recover the fine features of
the FEM solution, which might be due to the stiffness of problem, i.e.
requiring an accurate traction boundary fit on a first order automatic
differentiation level while needing to resolve the residual at a second order
level. Similarly, the same problems with DEM and PINN can be observed when
looking at the resulting stress components $P_{11}$ and $P_{12}$ (Figure 11).
In contrast, mDEM is able to achieve proficient stress fields compared to FEM.

Figure 9: Localized traction problem

(a) DEM

(b) PINN

(c) mDEM

(d) FEM

Figure 10: Localized Uniaxial loading problem displacement magnitude.

(a) DEM $P_{11}$

(b) PINN $P_{11}$

(c) mDEM $P_{11}$

(d) FEM $P_{11}$

(e) DEM $P_{12}$

(f) PINN $P_{12}$

(g) mDEM $P_{12}$

(h) FEM $P_{12}$

Figure 11: Localized uniaxial loading problem stress components $P_{11}$ (a-d)
and $P_{12}$ (e-f).

###  5.3 Beam with a circular hole

In the last example we look at a 2D beam with a circular hole, see Figure
13(a). On each edge $5000$ training points for the traction boundary were
sampled, whereas the domain training points are arranged in a point-grid of
$300\times 150$ samples where points lying in the hole were removed. The final
training point positions are displayed in Figure 13(b).

(a) Local traction boundary

(b) Beam with a circular hole

Figure 12: Normalized loss value convergence over iteration steps of PINN.

The displacement boundary conditions of PINN, DEM and mDEM are a-priori
fulfilled by defining

$\displaystyle\bm{u}_{1}(\bm{X},\bm{\theta})$
$\displaystyle=XYz_{1}(\bm{X},\bm{\theta}),$ (28)
$\displaystyle\bm{u}_{2}(\bm{X},\bm{\theta})$
$\displaystyle=Yz_{2}(\bm{X},\bm{\theta}).$

The displacement solutions in the $x-$ and $y-$ directions of all investigated
techniques are shown in Figure 14 respectively. DEM and mDEM are in close
agreement to the FEM solutions, whereas PINN does not appear to find the
correct solution. This is supported by the evolution of the loss values over
the training procedure as reported in Figure 12(b) which highlights a quick
saturation of the error value pointing towards stagnation in local minima.
When looking at the expected stress concentration features, here displayed for
the $P_{11}$ and $P_{12}$ components in Figure 15, it can be seen that even
though DEM is able to follow the overall stress field in an accurate manner,
only mDEM resolves the stress concentrations on comparable levels to FEM.

(a) Beam with defection

(b) Training point distribution

Figure 13: Problem setting and training point positions

(a) DEM $u_{x}$

(b) PINN $u_{x}$

(c) mDEM $u_{x}$

(d) FEM $u_{x}$

(e) DEM $u_{y}$

(f) PINN $u_{y}$

(g) mDEM $u_{y}$

(h) FEM $u_{y}$

Figure 14: Defection problem displacement components $u_{x}$ (a-d) and $u_{y}$
(e-f) .

(a) DEM $P_{11}$

(b) PINN $P_{11}$

(c) mDEM $P_{11}$

(d) FEM $P_{11}$

(e) DEM $P_{12}$

(f) PINN $P_{12}$

(g) mDEM $P_{12}$

(h) FEM $P_{12}$

Figure 15: Defection problem stress components $P_{11}$ (a-d) and $P_{12}$
(e-f).

##  6 Conclusion and outlook

This paper proposes an extension to the formulation of the Deep Energy Method
(DEM) called mixed Deep Energy Method (mDEM). It enhances the original
approach which uses a neural network as a global shape function of the
displacement over the computational domain by additionally defining the stress
components of the first Piola-Kirchhoff stress tensor as outputs of the
networks. This proves to be useful in combating the shortcomings of DEM with
regards to resolving concentration features. In order to make mDEM more
versatile we introduce a numerical integration approach based on Delaunay
tesselation that does not require the training point positions to be arranged
in grid-like fashion.
We test the proposed approach on three computational experiments which exhibit
stress concentrations and compare the results to DEM, a Physics-Informed
Neural Network and FEM solutions. It can be seen that mDEM is able to yield
comparable results to FEM in all investigated problems, whereas DEM and PINN
prove to not be reliably able to resolve local stress and displacement
features. In future works we plan to apply the presented formulation to three
dimensional applications, time-dependent computational experiments as well as
problems governed by non-elastic constitutive laws as well as for inverse
problems focusing on parameter estimation.

## References

  * Abueidda et al. (2020) Abueidda, D.W., Lu, Q., Koric, S., 2020.  Deep learning collocation method for solid mechanics: Linear elasticity, hyperelasticity, and plasticity as examples.  arXiv preprint arXiv:2012.01547 .
  * Alnæs et al. (2015) Alnæs, M.S., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Richardson, C., Ring, J., Rognes, M.E., Wells, G.N., 2015\.  The fenics project version 1.5.  Archive of Numerical Software 3.  doi:10.11588/ans.2015.100.20553.
  * Bottou et al. (2018) Bottou, L., Curtis, F.E., Nocedal, J., 2018.  Optimization methods for large-scale machine learning.  Siam Review 60, 223–311.
  * Eggersmann et al. (2019) Eggersmann, R., Kirchdoerfer, T., Reese, S., Stainier, L., Ortiz, M., 2019.  Model-free data-driven inelasticity.  Computer Methods in Applied Mechanics and Engineering 350, 81–99.
  * Fuhg et al. (2021) Fuhg, J.N., Boehm, C., Bouklas, N., Fau, A., Wriggers, P., Marino, M., 2021\.  Model-data-driven constitutive responses: application to a multiscale computational framework.  arXiv:2104.02650.
  * González et al. (2019) González, D., Chinesta, F., Cueto, E., 2019.  Thermodynamically consistent data-driven computational mechanics.  Continuum Mechanics and Thermodynamics 31, 239–253.
  * Goodfellow et al. (2016) Goodfellow, I., Bengio, Y., Courville, A., Bengio, Y., 2016\.  Deep learning. volume 1.  MIT press Cambridge.
  * Haghighat et al. (2020) Haghighat, E., Raissi, M., Moure, A., Gomez, H., Juanes, R., 2020.  A deep learning framework for solution and discovery in solid mechanics.  arXiv preprint arXiv:2003.02751 .
  * Hernandez et al. (2021) Hernandez, Q., Badías, A., González, D., Chinesta, F., Cueto, E., 2021.  Deep learning of thermodynamics-aware reduced-order models from data.  Computer Methods in Applied Mechanics and Engineering 379, 113763.
  * Hornik et al. (1989) Hornik, K., Stinchcombe, M., White, H., 1989.  Multilayer feedforward networks are universal approximators.  Neural networks 2, 359–366.
  * Huang et al. (2020) Huang, D., Fuhg, J.N., Weißenfels, C., Wriggers, P., 2020\.  A machine learning based plasticity model using proper orthogonal decomposition.  Computer Methods in Applied Mechanics and Engineering 365, 113008.
  * Ibañez et al. (2017) Ibañez, R., Borzacchiello, D., Aguado, J.V., Abisset-Chavanne, E., Cueto, E., Ladeveze, P., Chinesta, F., 2017.  Data-driven non-linear elasticity: constitutive manifold construction and problem discretization.  Computational Mechanics 60, 813–826.
  * Jagtap and Karniadakis (2020) Jagtap, A.D., Karniadakis, G.E., 2020\.  Extended physics-informed neural networks (xpinns): A generalized space-time domain decomposition based deep learning framework for nonlinear partial differential equations.  Communications in Computational Physics 28, 2002–2041.
  * Kadeethum et al. (2021) Kadeethum, T., Ballarin, F., Bouklas, N., 2021.  Non-intrusive reduced order modeling of poroelasticity of heterogeneous media based on a discontinuous galerkin approximation.  arXiv preprint arXiv:2101.11810 .
  * Kadeethum et al. (2020) Kadeethum, T., Jørgensen, T.M., Nick, H.M., 2020.  Physics-informed neural networks for solving nonlinear diffusivity and biot’s equations.  PloS one 15, e0232683.
  * Kharazmi et al. (2021) Kharazmi, E., Zhang, Z., Karniadakis, G.E., 2021.  hp-vpinns: Variational physics-informed neural networks with domain decomposition.  Computer Methods in Applied Mechanics and Engineering 374, 113547.
  * Kingma and Ba (2014) Kingma, D.P., Ba, J., 2014.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980 .
  * Kirchdoerfer and Ortiz (2016) Kirchdoerfer, T., Ortiz, M., 2016\.  Data-driven computational mechanics.  Computer Methods in Applied Mechanics and Engineering 304, 81–101.
  * Lagaris et al. (1998) Lagaris, I.E., Likas, A., Fotiadis, D.I., 1998.  Artificial neural networks for solving ordinary and partial differential equations.  IEEE transactions on neural networks 9, 987–1000.
  * Nguyen-Thanh et al. (2020) Nguyen-Thanh, V.M., Zhuang, X., Rabczuk, T., 2020.  A deep energy method for finite deformation hyperelasticity.  European Journal of Mechanics-A/Solids 80, 103874.
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., Chintala, S., 2019.  Pytorch: An imperative style, high-performance deep learning library, in: Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., Garnett, R. (Eds.), Advances in Neural Information Processing Systems 32. Curran Associates, Inc., pp. 8024–8035.  URL: http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf.
  * Raissi et al. (2019) Raissi, M., Perdikaris, P., Karniadakis, G.E., 2019.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics 378, 686–707.
  * Rao et al. (2020) Rao, C., Sun, H., Liu, Y., 2020.  Physics informed deep learning for computational elastodynamics without labeled data.  arXiv preprint arXiv:2006.08472 .
  * Samaniego et al. (2020) Samaniego, E., Anitescu, C., Goswami, S., Nguyen-Thanh, V.M., Guo, H., Hamdia, K., Zhuang, X., Rabczuk, T., 2020\.  An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications.  Computer Methods in Applied Mechanics and Engineering 362, 112790.
  * Weinan and Yu (2018) Weinan, E., Yu, B., 2018.  The deep ritz method: a deep learning-based numerical algorithm for solving variational problems.  Communications in Mathematics and Statistics 6, 1–12.