# FEM-based Real-Time Simulations of Large Deformations with Probabilistic
Deep Learning

Saurabh Deshpande<EMAIL_ADDRESS>
Department of Engineering; Faculty of Science, Technology and Medicine;
University of Luxembourg Jakub Lengiewicz<EMAIL_ADDRESS>
Department of Engineering; Faculty of Science, Technology and Medicine;
University of Luxembourg
and
Institute of Fundamental Technological Research, Polish Academy of Sciences
Stéphane P.A. Bordas<EMAIL_ADDRESS>
Department of Engineering; Faculty of Science, Technology and Medicine;
University of Luxembourg

###### Abstract

For many engineering applications, such as real-time simulations or control,
conventional solution techniques of the underlying nonlinear problems are
usually computationally too expensive. In this work we propose a highly
efficient deep-learning surrogate framework that is able to predict the
response of hyper-elastic bodies under load. The surrogate model takes a form
of special convolutional neural network architecture, so-called U-Net, which
is trained with force-displacement data obtained with the finite element
method. We propose deterministic- and probabilistic versions of the framework,
and study it for three benchmark problems. In particular, we check the
capabilities of the Maximum Likelihood and the Variational Bayes Inference
formulations to assessing the confidence intervals of solutions.

Keywords: Convolutional Neural Network, Bayesian Inference, Deep Learning,
Hyper-elastic, Finite Element Method, Real-Time Simulations, Variational
Inference

##  1 Introduction

Computer simulation is omnipresent in practical or industrial contexts, and it
has also become one of the most essential tools for scientific research in
general. This is especially true in engineering, where predictive
computational models of materials or whole systems are crucial in the process
of design, optimization or control. Some of these applications, like
biomedical simulations used to support or train surgeons (Courtecuisse et al.,
2014), (Wu et al., 2015), (Bui et al., 2018), (Cotin et al., 1999),
(Delingette et al., 1999), require computationally-efficient or even real-time
solutions. Finding convenient trade-offs between the accuracy and response
time of such computational models is currently an active area of research, and
also one of motivations for the research presented in this work.

Within approaches to reduce the computation time of complex engineering
simulations, we can find extremely simplified models and approximated
computational schemes, such as linear mass-spring systems with simplified
implicit dynamic, enabling, for instance, real time simulation of cloth
movements (Kang et al., 2000), muscle deformations (Nedel and Thalmann, 1998),
(Zhang et al., 2001) or whole surgical procedures (Patete et al., 2013).
However, they generally significantly compromise on accuracy, which makes
their predictive capabilities greatly limited.

When the accuracy is important, undoubtedly, the most general and widely used
methodology in engineering for solving boundary value problems is the finite
element method (FEM) (Zienkiewicz and Taylor, 1991). This accuracy, especially
when it comes to highly non-linear or history-dependent problems, may require
a significant computational effort. Although the advancements in hardware
development and software optimization enabled to some extent speeding up FEM
computations, for many applications it is still insufficient. This is why
various specialized FEM-based approaches have been developed to cut down the
prediction time, indispensably at the cost of sacrificing the accuracy, see,
e.g., (Marinkovic and Zehn, 2019). Below, we briefly characterize some of
them.

One notable way of speeding up FEM solutions is use of domain decomposition
methods, such as Finite Element- or Boundary Element Tearing and
Interconnecting (FETI/BETI) (Farhat and Roux, 1991), (Langer U., 2003). In
these methods, global domain is partitioned into a set of disconnected sub-
domains, which are computed in parallel on different processors. Lagrange
multipliers are introduced to enforce compatibility at the interface nodes of
different partitions. However, in many applications it is not possible to meet
real-time responses on the hardware available for industrial consumers. This
is due to a limited number of available cores and a significant communication
burden that deteriorates the overall time performance of the method.

Another way to speeding up FEM computations is to use model order reduction
(MOR) methods. The general concept of MOR, e.g., applied to a discretized FEM
formulation, is to ﬁnd a low dimension subspace in order to approximate the
full space at an acceptable loss of accuracy. This potentially enables to
controlling the trade-off between accuracy and computation time. Also MOR was
able to be adapted to work within the large-deformation regime, see, e.g.,
(Kerfriden et al., 2011), which was for instance applied to simulate and
control soft robotic arms (Goury and Duriez, 2018) or to reduce computational
costs in nonlinear fracture mechanics problems (Kerfriden et al., 2012).
However, the efficiency of MOR-based techniques deteriorates when large-
deformations are present.

We can distinguish yet another family of FEM-based approaches, in which the
expected speedups and approximation capabilities originate from underlying
Deep Neural Networks (DNNs) with Deep Learning (DL) techniques used to train
these networks. Generally, DL approaches make an important part of machine
learning techniques, and have allowed to solving highly complex problems that
had eluded scientists for decades. In particular, DL-based methods have also
been developed to efficiently solving problems in engineering (Goodfellow et
al., 2016). One of popular approaches utilises the idea of the so-called
Physic Informed Neural Networks (PINNs), in which both the data (either
synthetic or experimental) and the assumed governing Partial Differenntial
Equations (PDEs) are incorporated in the training phase; for early traces of
these see, e.g., (Lagaris et al., 1998), (Lagaris et al., 2000), (McFall and
Mahan, 2009), and for a recent study see, e.g., (Raissi et al., 2019). One of
benefits of this approach is that possibly much less data is needed for
training, which can be an important factor for many data-driven applications.
Note however, that even if the physics is not explicitly enforced in the
training phase (non-PINN case), it can still be recovered in the trained model
by implicitly following large amount of training data. In the case when all
training data are synthetically provided from FEM simulations, see, e.g.,
(Lorente et al., 2017), (Mendizabal et al., 2019), we will refer to it as the
_direct FEM-based approach_.

In this work we propose a framework that falls into the above mentioned class
of direct FEM-based DNN approaches. The framework is based on a particular DNN
architecture—the U-Net architecture (Ronneberger et al., 2015)—which in turn
can be viewed as a type of Convolutional Neural Networks (CNNs); further
explanations are provided later in this work. Originally, the U-Net
architecture has been developed for the purpose of biomedical image
segmentation, however, it turned out to be also suitable for other
applications. In particular, the present work is inspired by the recent
results of (Mendizabal et al., 2019), in which the authors demonstrate quite
accurate real-time non-linear force-displacement predictions done by U-Nets
trained on FEM-based data. It has been noticed, see (He and Xu, 2019), (Wang
et al., 2020), that this good accuracy is not totally accidental but can be
possibly linked to strong resemblances of U-Net architectures to multi grid
solution schemes (Scott, 2008). Such point of view places the U-Net approach
less as a brute-force black-box approximation and more as a suited solution
scheme, which makes this line of research very promising.

Another aspect that is studied in the present work is capabilities of U-Nets
to quantify uncertainties. This is motivated by the fact that in various real
life applications, such as surgical simulations (Bui et al., 2018) or
autonomous driving (McAllister et al., 2017), it is very important to produce
reliable uncertainty estimates in addition to the predictions. Their lack can
possibly lead to harmful consequences in these critical tasks. However,
deterministic neural networks are usually certain about their predictions, and
this overconfidence is especially evident when facing data away from the
training set (extrapolated region). Generally, the uncertainties can be
categorised as the ones associated to misfit of neural network models (the so
called, epistemic uncertainties) and the ones that refer to the noise in
training data (the so called, aleatoric uncertainties). In addition to that,
we can make another split, i.e., into uncertainties within and outside the
data region (interpolated and extrapolated regions, respectively). Within the
data region, the prediction uncertainties are caused both by the noisy data
and the the model misfit, which can be captured and quantified in many ways,
with the Maximum Likelihood Estimation (MLE) method as one of the most
straightforward (Kendall and Gal, 2017). However, for the extrapolated region
we have no data support, therefore, no direct quantification can be done
there. What we can only reasonably assume is that the uncertainty should
generally increase when moving away from the data region. To achieve this, in
this work we will extend the idea proposed in (Duerr et al., 2020)(Gal, 2016)
which relies on converting a neural network to its stochastic counterpart by
replacing discrete parameters by probability distributions and using a special
training technique that is based on Bayesian Inference.

As stated above, as a general methodology we focus on Bayesian Inference and
related Bayesian Neural Networks (BNNs). Indeed, Bayesian approaches have been
already considered in the context of FEM models; for instance in (Hauseux et
al., 2018) uncertainties are quantified for hyperelastic soft tissues by
incorporating stochastic parameters in a respective FEM model, while in
(Rappel et al., 2020) a thorough tutorial on using Bayesian Inference to solid
mechanics problems is provided. However, the above mentioned approaches use
Markov Chain Monte Carlo (MCMC) type of methods, which are computationally
intractable in the context of deep networks with millions of parameters. For
that reason, in this work we use one of the well-known methods of performing
approximate Bayesian inference—Variational Inference (VI) (Graves, 2011)—in
which an approximate distribution is used instead of the true Bayesian
posterior over model parameters.

To sum up, the scope of the present work is to study the applicability of
U-Net deep learning architectures to performing real-time predictions for
large-deformation problems with uncertainties, where synthetic training data
is provided by FEM simulations. The organisation of the paper is the
following. In Section 2 we present the general methodology applied to a
deterministic version of U-Net. In Section 3 we introduce the extension of the
framework to the Variational Bayesian Inference case. Then, in Section 4, an
extensive study of the proposed framework is performed, which is based on
several 2D and 3D benchmark examples. The conclusions and future research
directions are outlined in Section 5.

##  2 General FEM-based U-Net Methodology

The proposed approach can be divided into two main phases. The first phase
involves finite element simulations to prepare necessary datasets. The second
phase consists of building and training deterministic/probabilistic U-Net deep
neural networks, using the training datasets generated in the first phase. The
trained U-Nets are then used as surrogate models, and their predictions are
validated with the testing datasets.

Figure 1: Continuum problem is discretized by FEM mesh, training examples are
generated by applying random point loads on Neumann boundary, they are given
as input to the NN. Trained NN (blue mesh) then predicts the deformation for
new force. FEM mesh (red) is provided for comparision of the solution,
undeformed mesh is displayed in gray

###  2.1 FEM-Based Deep Learning Approach

As a problem to be solved, we consider a boundary value problem of deformation
of a hyperelastic solid, with a constant prescribed Dirichlet BC. Large body
deformations exhibit non-linear stress-strain behaviour when applied with
external forces, hence one needs to consider hyper-elastic constitutive laws
for simulating such systems.

Consider a boundary value problem in continuum mechanics in domain $\Omega$,
Dirichlet and Neumann boundary conditions are applied on $\Gamma_{D}$,
$\Gamma_{N}$ respectively. Neglecting the body forces, the virtual work
principle for nonlinear elastostatic equation reads

$\int_{\Omega}{\bf
P}(\mbox{\boldmath{$u$}})\cdot\nabla\delta\mbox{\boldmath{$u$}}\,dV-\int_{\Gamma_{t}}\mbox{\boldmath{$\bar{t}$}}\cdot\delta\mbox{\boldmath{$u$}}\,dS=0\quad\quad\forall\delta\mbox{\boldmath{$u$}},$
(1)

where $u$ and $\delta\mbox{\boldmath{$u$}}$ belong to appropriate functional
spaces, $\mbox{\boldmath{$u$}}=\bar{\mbox{\boldmath{$u$}}}$ and
$\delta\mbox{\boldmath{$u$}}={\bf 0}$ on $\Gamma_{u}$, and ${\bf
P}(\mbox{\boldmath{$u$}})$ is the first Piola-Kirchhoff stress tensor. The
required constitutive relationship will be defined through the (hyper-)elastic
strain energy potential $W(\mbox{\boldmath{$F$}})$ as

${\bf P}(\mbox{\boldmath{$F$}})=\frac{\partial
W(\mbox{\boldmath{$F$}})}{\partial\mbox{\boldmath{$F$}}}$ (2)

where
$\mbox{\boldmath{$F$}}=\mbox{\boldmath{$I$}}+\nabla\mbox{\boldmath{$u$}}$ is
the deformation gradient tensor. (The particular form of $W$, used in this
work, will be presented in Section 4.1.1.)

After standard FE discretization, the problem expressed by Eq. (1) will take
the form of system of non-linear equations

$\mathbf{R}(\mathbf{u})=\mathbf{f}_{\text{int}}({\bf
u})-\mathbf{f}_{\text{ext}}={\bf 0},$ (3)

which expresses the balance between external and internal nodal forces. By
solving the system of equations (3) (e.g., with the Newton-Raphson method) for
a given external force vector $\mathbf{f}_{\text{ext}}={\bf f}_{i}$ we obtain
a solution in the form of nodal displacements $\mathbf{u}_{i}$.

External forces can be applied to a selected region on the surface described
by $\Gamma_{t}$. For the current framework, we consider a single FE
discretization of a given domain $\Omega$, and apply large number of forces
with a prescribed range of magnitudes. As described in Figure 1, we apply
random loads vectors ${\bf f}_{i}$ on the nodes present in $\Gamma_{t}$ to
generate non-linear nodal displacements ${\bf u}_{i}$. The dataset
$\mathcal{D}={({\bf f}_{i},{\bf u}_{i})}_{i=1}^{N}$ consisting pairs of nodal
force and displacement vectors is treated as an input to train DNNs (the input
of the neural network is a vector of nodal forces, and the predicted output is
a vector of nodal displacements).

###  2.2 U-Net deep neural network architecture

U-Nets (Ronneberger et al., 2015) compose a sub-class of convolutional neural
networks. They owe their name to a specific U-shape of the architecture
diagram, e.g., see Figure 2, which is an effect of applying cascades of max
pooling operations, followed by cascades of upsampling operations.

At its input, the layer $\mbox{\boldmath{$d$}}^{0}$, the U-Net network,
$\mathcal{U}$, accepts the vector of external forces, ${\bf f}$, in the
original mesh format $n_{x}\times n_{y}\times n_{z}\times 3$, where $n_{x}$,
$n_{y}$, $n_{z}$ are the dimensions of the structured 3D mesh (for 2D problems
the format is $n_{x}\times n_{y}\times 2$). The output displacements, in the
layer $\mbox{\boldmath{$d$}}^{L}$, are in the analogous mesh format. The model
parameters $\theta$ are all weights (kernel $k$ and biases $b$) of the neural
network (see below for details).

For 2D problems, we use architecture as described in Figure 2. To a given
input mesh, first we add double zero padding in each spatial direction (This
is done to avoid the loss of information of corner nodes). There are two
layers in each encoding and decoding phases. To the padded input, we apply two
convolutions with batch normalisation followed by rectified linear unit
activation (ReLU), see Eq. (4). In encoding phase, at each step, $2\times 2$
Max Pooling layers are applied which decrease spatial dimensions by half, and
we multiply the number of channels by 2 (c=128 for first level). In decoding
steps, $2\times 2$ Up-sampling layers are applied which increase the spatial
dimension by two and no of channels are decreased by half, till we get back
the input dimensions. At each level, encoding and decoding outputs are
concatenated together. At the end $1\times 1$ convolution is applied with
linear activation to get the output.

Figure 2: A schematic of exemplary U-Net architecture for 2D domains,
($n_{\text{x}},n_{\text{y}}$) stand for number of nodes in x, y direction of
2D domain. First step contains ’$c$’ channels.

In order to better understand the idea of convolution operator, the $3\times
3\times c$ operator in 2D on the first U-Net level can be imagined as a filter
window that is applied to a $n_{x}\times n_{y}\times c$ mesh. For a two
dimensional domain, input tensor is of the dimension $n_{x}\times n_{y}\times
2$, which is identical to FEM mesh. Here $2$ stands for number of channels of
input convolutions, $n_{x},n_{y}$ stand for number of nodes in respective $x$,
$y$ direction. For the maximum leverage of CNNs, we keep $x$ & $y$-dofs in
separate channels. Convolutional filter is operated on local region and then
is slided along spatial directions $x,y$ with a stride of 1 as illustrated in
Figure 3. For a $3$D domain, $z$ channel is stacked at the back of $x,y$ and
local convolutional filter is operated in similar way as in Figure 3.

Figure 3: $x$ and $y$ dofs are stored in different channels, $3\times 3$
convolutional filter (gray) acts locally along the channel direction and it
slides along with the step of 1 in both horizontal and vertical directions
(red).

An example of a non-restrictive convolution operation (in 2D), between
subsequent U-Net layers $n$ and $n+1$, for the filter size
$3\times{}3\times{}c^{n}$ reads

$d^{n+1}_{p,r,\beta}=\mathcal{A}\left(b^{n+1}_{p,r,\beta}+\sum_{i=p-1,}^{p+1}\sum_{j=r-1,}^{r+1}\sum_{\alpha=1}^{c^{n}}d^{n}_{i,j,a}k^{n+1}_{i,j,\alpha,p,r,\beta}\right),$
(4)

where $d^{n}_{i,j,\alpha}$ and $d^{n+1}_{p,r,\beta}$ are neural network nodes
at layer $n$ and $n+1$, respectively, the weights
$k^{n+1}_{i,j,\alpha,p,r,\beta}$ are parameters of the convolution operator,
the weights $b^{n+1}_{p,r,\beta}$ are biases at a layer $n+1$, and
$\mathcal{A}(\cdot)$ is an activation function  (ReLU). The co-ordinates $i,j$
and $p,r$ are related to spacial configurations of nodes in a 2D grid, and in
the input and output layers, respectively, they directly refer to the
underlying structured FEM mesh. In our case, we add zero pad in each dimension
of input before applying the convolution, to ensure the same size of the input
and output. The indices $1\leq{}\alpha\leq{c^{n}}$ and
$1\leq{}\beta\leq{c^{n+1}}$ represent the channel number. Note that the number
of channels in subsequent layers do not need to be equal, i.e., in general
$c^{n}\neq{}c^{n+1}$. Note also that in the first and in the last layer, the
number of channels correspond to the spacial dimension of problem (2D or 3D).

Max-pooling operation is responsible for reducing spacial dimensions of its
input, channel dimensions are unaffected by it. It reads as follows

$d^{n+1}_{p,r,\alpha}=\max\limits_{\begin{subarray}{c}2p-1\leqslant i\leqslant
2p\\\\[1.39998pt] {2r-1\leqslant j\leqslant
2r}\end{subarray}}d_{i,j,\alpha}^{n}$ (5)

Upsampling operator can be seen as reverse of maxpooling, it increases the
spatial dimensions of the input by desirable factor without affecting the
channel dimension. As showed in Figure 2, in the decoder phase, outputs of up-
sampling are concatenated with respective layers from encoder phase of the
U-Net (in case of symmetric U-Nets, the $n$-th layer is concatenated with the
$(L-n-2)$-th layer, where $L$ is the index of output layer in a U-Net). The
up-sampling with concatenation read

$d_{i,j,\alpha}^{n+1}=\begin{cases}d_{{\lfloor}i/2{\rfloor}+1,{\lfloor}j/2{\rfloor}+1,\alpha}^{n}~{}~{},&1\leqslant\alpha\leqslant
c^{n}\\\\[10.00002pt]
d_{i,j,(\alpha-c^{n})}^{L-n-2}~{}~{}&c^{n}+1\leqslant\alpha\leqslant
c^{n}+c^{L-n-2}\end{cases}$ (6)

The final $1\times 1$ convolution operation reads

$d^{L}_{i,j,\beta}=b^{L}_{i,j,\beta}+\sum_{\alpha=1}^{c^{L-1}}k^{L}_{i,j,\alpha,\beta}d^{L-1}_{i,j,\alpha},\quad\quad\beta=1,2.$
(7)

For the 3D version of U-Nets, the operations given by Equations (4)-(7) are
straightforwardly extended by one additional dimension. This results in adding
one index to nodes’ and biases’ specifications, $d$ and $b$, respectively, and
two indexes to $3\times{}3$ convolution weights, $k$.

We denote by $\mathcal{U}({\bf f},\mbox{\boldmath{$\theta$}})={\bf d}^{L}$ the
prediction of a deterministic U-Net for a given input ${\bf f}$ and a given
set of neural network parameters

$\mbox{\boldmath{$\theta$}}=\bigcup_{n=1}^{L}\\{{\bf k}^{n},{\bf b}^{n}\\}.$
(8)

$\mathcal{U}({\bf f},\mbox{\boldmath{$\theta$}})$ is defined recursively (the
forward propagation), starting from the input layer
$\mbox{\boldmath{$d$}}^{0}={\bf f}$, then subsequently applying appropriate
transformations given by one of Eqs. (4)-(6), and finally applying the
transformation given by Eq. (7).

For a given training dataset $\mathcal{D}=\\{({\bf f}_{1},{\bf
u}_{1}),...,({\bf f}_{N},{\bf u}_{N})\\}$, the U-Net is trained by minimizing
the following mean squared error loss function

$\mathcal{L}_{\text{det}}(\mathcal{D},\mbox{\boldmath{$\theta$}})=\frac{1}{N}\sum_{i=1}^{N}\lVert\mathcal{U}({\bf
f}_{i},\bm{\theta})-{\bf u}_{i}\rVert_{2}^{2}$ (9)

which gives the optimal parameters

$\bm{\theta}^{*}=\text{arg
}\underset{\bm{\theta}}{\text{min}}~{}\mathcal{L}_{\text{det}}(\mathcal{D},\mbox{\boldmath{$\theta$}}).$
(10)

A particular training strategy, used in this work, will be introduced in
Section 4.1.2.

##  3 Probabilistic U-Net framework

###  3.1 Variational Bayesian Inference

Bayesian methods provide an approach to quantify uncertainty of prediction in
deep neural networks. In order to do so, in this framework, we replace the
deterministic weights by probability distributions, which helps in
incorporating uncertainty in neural networks (Blundell et al., 2015).
Originally this idea was only used to prevent overfitting, but was observed to
also increase variability of outputs in the extrapolated region. In order to
suitably control the level of introduced perturbations to parameters, we use
Bayesian Inference. This gives us a formal theoretical framework, allowing us
to apply suitable computational techniques (VI) to efficiently training
networks and predicting results. The input of a network remains the same, but
some of model parameters become probability distributions (stochastic), and
for that reason also the output of network becomes a distribution over
possible outputs. As per the standard Bayesian approach, we specify a prior
distribution $P(\bm{w})$ over parameters, we consider $\mathcal{D}=\\{({\bf
f}_{1},{\bf u}_{1}),...,({\bf f}_{N},{\bf u}_{N})\\}$ as the given training
dataset, then for a new vector ${\bf f}^{*}$, prediction ${\bf u}^{*}$ is
given by

$P({\bf u}^{*}|{\bf f}^{*},\mathcal{D})=\int P({\bf u}^{*}|{\bf
f}^{*},\bm{w})\cdot P(\bm{w}|\mathcal{D})\text{d}\bm{w}$ (11)

Bayesian inferencing involves the calculation of true parameter posterior
$P(\bm{w}|\mathcal{D})$ conditioned over the training data. Standard Markov
Chain Monte Carlo (MCMC) are more accurate methods, but they are
computationally very expensive and don’t scale well with the parameter size.
Hence variational inference (VI) is used to approximate true posterior
densities in bayesian neural networks (Graves, 2011). True posterior is
approximated by a variational posterior $q(w|\theta)$, parametrized by
$\theta$. Variational learning finds optimal parameters $\theta^{*}$ by
minimizing the Kullback-Leibler divergence (KL-divergence) between true and
variational posteriors, as per the following equation:

$\begin{split}\bm{\theta}^{*}&=\text{arg
}\underset{\bm{\theta}}{\text{min}}~{}\text{KL}\big{[}P(\bm{w}|{\bf f},{\bf
u})||q(\bm{w}|\bm{\theta})\big{]}\\\ &=\text{arg
}\underset{\bm{\theta}}{\text{min}}~{}\int\log\frac{q(\bm{w}|\bm{\theta})}{P(\bm{w})P(\mathcal{D}|\bm{w})}dw\\\
&=\text{arg
}\underset{\bm{\theta}}{\text{min}}~{}\text{KL}[q(\bm{w}|\bm{\theta})||P(\bm{w})]-\mathbb{E}_{q(\bm{w}|\bm{\theta})}[\log
P(\mathcal{D}|\bm{w})].\end{split}$ (12)

The resulting cost function is the loss function for training the neural
network. It consists of two parts, first is the prior dependent part
represented by KL-divergence term, it can be referred as model complexity
cost; it tells how close approximate posteriors are to priors. And later is
the data dependent part which can be referred as likelihood cost, it tells how
well network fits the data. Bayesian neural networks with prior distributions
are well known to induce regularisation effect (Neal, 1996); in particular,
using Gaussian priors is equivalent to weight decay ($L_{2}$ regularization)
(Vladimirova et al., 2019).

During the forward pass, weights are sampled from the variational posterior
$q(\bm{w}|\bm{\theta})$ but backpropagation can’t flow through random samples,
to avoid this reparameterization trick is used (Kingma et al., 2015). Sampling
of weights $w$ can be obtained by sampling a parameter free distribution i.e.
unit Gaussian, and shifting it by a mean $\bm{\mu}$ and scaling by a standard
deviation $\bm{\sigma}$. We parameterise the standard deviation point-wise as
$\bm{\sigma}=\log(1+\text{exp}(\bm{\rho}))$ to have $\bm{\sigma}$ always non-
negative. Thus the sample is
$\bm{w}=\bm{\mu}+\log(1+\text{exp}(\bm{\rho}))\odot\bm{\epsilon}$, where
$\odot$ is point-wise multiplication, $\bm{\epsilon}$ is drawn from
$\bm{\mathcal{N}}(\bm{0},\bm{I})$ and the variational posterior parameters are
$\bm{\theta}=(\bm{\mu},\bm{\rho})$ (Blundell et al., 2015).

###  3.2 Maximum likelihood estimation

In addition to Bayesian approach, we introduce Maximum Likelihood Estimation
(MLE) to cross-compare both methods about their capabilities to quantify
uncertainty in case of noisy data. MLE is a non-Bayesian approach, in the
sense that we do not place priors on parameters of the network, they are
treated as deterministic. It is a probabilistic model, in our case, the model
gives Gaussian distribution predictions. Loss function for MLE can be
recovered from the Equation (12) by removing KL divergence part. Thus MLE
model with nonconstant standard deviation can be fitted with a loss that
corresponds to the negative log-likelihood (NLL) of a normal probability
model, where both parameters (mean and standard deviation) depend on the
inputs and can be computed by an NN with double number of outputs yielding the
distribution as the output

$\bm{\theta}^{\text{MLE}}=\text{arg }\underset{\bm{\theta}}{\text{min}}-\log
P(\mathcal{D}|\bm{\theta}).$ (13)

###  3.3 Use of Empirical Bayes

Since NN parameters are latent variables of the model, it is very difficult to
make a proper choice of priors. If one sets the priors far from their true
values, then the posterior may be unduly affected by such choice. To overcome
this, we incorporate Empirical Bayes(EB) approach (Carlin and Louis, 1997),
where we use the observed data to estimate the prior hyper-parameters.
Otherwise it can be tough to do a good mean fit without a good prior knowledge
of network parameters. In training phase we update the prior mean, keeping the
prior standard deviation constant. Hence we optimise the loss by considering
gradients with respect to prior mean as well.

###  3.4 Loss functions for probabilistic U-Net

We modify the deterministic U-Net architecture by replacing some of the layers
with probabilistic layers, as a result output of the network is a probability
distribution itself. We choose Gaussian distributions to represent priors and
approximate posteriors of probabilistic layers, it is well known that use of
Gaussian prior induces L$2$-regularisation effect (Vladimirova et al., 2019).
Through various experiments for a completely probabilistic network, we
observed to have an overall bad fit to the data. Therefore we replace only
half of the deterministic layers with probabilistic ones in order to get a
good mean fit.

For the Bayesian CNN, we use loss function as given in Eq. (12). Expectations
of the Eq. (12) are approximated by $\mathcal{M}$ Monte Carlo samples drawn
from the approximate posterior $q(\bm{w}|\bm{\theta})$ as refereed in Eq. (14)

$\begin{split}\mathcal{L}_{\text{VB}}&=~{}\text{KL}[q(\bm{w}|\bm{\theta})||P(\bm{w})]-\mathbb{E}_{q(\bm{w}|\bm{\theta})}[\log
P(\mathcal{D}|\bm{w})]\\\ &\approx~{}\sum_{i=1}^{\mathcal{M}}\log
q(\bm{w}^{(i)}|\bm{\theta})-\log P(\bm{w}^{(i)})-\log
P(\mathcal{D}|\bm{w}^{(i)})\end{split}$ (14)

If we substitute Gaussian probability density functions, expression in the RHS
turns out to be as given in Eq. (15). Let us say we have ’$\mathcal{G}$’
Gaussian distributions for our Bayesian NN, every distribution is expressed in
terms its mean and standard deviation values, which are the parameters for
Bayesian CNN. Since standard deviations $\bm{\sigma}$ needs to be positive, we
first train the network on untransformed standard deviations $\bm{\rho}$ which
are later transformed to $\bm{\sigma}$ through soft-plus function. Also, for
the reasons discussed in Section 3.3, we involve prior means in the training
procedure as well. Hence parameters to be learned in the training procedure
are $\bm{\theta}=(\bm{\mu},\bm{\rho},\bm{\mu_{p}})$. $N,n$ are total number of
training examples and dof per problem respectively. Bayesian neural network
gives Gaussian distribution output, which is characterised by mean and
variance ($\mathcal{U}_{\mu}({\bf f},\bm{w}),\mathcal{U}_{\rho}({\bf
f},\bm{w})$). Finally the loss function for the variational bayes is given as
follows

$\begin{split}\mathcal{L}_{\text{VB}}(\mathcal{D},\bm{\theta})\approx~{}\sum_{i=1}^{\mathcal{M}}&\left[\sum_{j=1}^{\mathcal{G}}\left(\log\left(\sqrt{2\pi}~{}\sigma_{j}^{(i)}\right)-\frac{\left(w_{j}^{(i)}-\mu_{j}^{(i)}\right)^{2}}{2\left(\sigma_{j}^{(i)}\right)^{2}}-\log\left(\sqrt{2\pi}~{}\sigma_{p}\right)+\frac{\left(w_{j}^{(i)}-(\mu_{p})_{j}^{(i)}\right)^{2}}{2\sigma_{p}^{2}}\right)\right.\\\
&-\left.\sum_{k=1}^{N}\sum_{l=1}^{n}\left(\log\left(\sqrt{2\pi}~{}\mathcal{U}^{(l)}_{\sigma}({\bf
f}^{(k)},\bm{w}^{(i)})\right)-\frac{\left(\hat{u}^{(k)}_{l}-\mathcal{U}^{(l)}_{\mu}({\bf
f}^{(k)},\bm{w}^{(i)})\right)^{2}}{2\left(\mathcal{U}^{(l)}_{\sigma}({\bf
f}^{(k)},\bm{w}^{(i)})\right)^{2}}\right)\right]\end{split}$ (15)

where

$\begin{split}\mathcal{U}^{(l)}_{\sigma}({\bf
f}^{(k)},\bm{w}^{(i)})&=\log(1+\exp({\mathcal{U}^{(l)}_{\rho}({\bf
f}^{(k)},\bm{w}^{(i)})})),\\\
w_{j}^{(i)}&=\mu_{j}^{(i)}+\sigma_{j}^{(i)}\epsilon_{j}^{(i)},\hskip
43.36464pt\epsilon_{j}^{(i)}\mathcal{~{}}\sim~{}{\mathcal{N}}(0,1),\\\
\sigma_{j}^{(i)}&=\log(1+\exp({\rho_{j}^{(i)}})).\end{split}$ (16)

$\sigma_{p}$ stands for the standard deviation of prior, which is kept
constant in the training procedure. Where as $\bm{\mu_{p}}$, mean of the prior
is updated with every iteration during the back propagation step as mentioned
in Section 3.3. Once we obtain optimized parameters through the training
process, we replace true posterior with variational posterior
$q(\bm{w}|\bm{\theta})$ in Eq. (11)

$P({\bf u}^{*}|{\bf f}^{*},\mathcal{D})\approx\int P({\bf u}^{*}|{\bf
f}^{*},\bm{w})\cdot
q(\bm{w}|\bm{\theta})\text{d}\bm{w}\approx\frac{1}{T}\sum_{t=1}^{T}P({\bf
u}^{*}|{\bf f}^{*},\widetilde{\bm{w}_{t}})$ (17)

Resultant predictive distribution can be approximated by Monte Carlo
integration as given in Eq. (17) with $\widetilde{\bm{w}_{t}}\sim
q(\bm{w}|\bm{\theta})$. Hence the mean and variance of the prediction is
obtained by taking mean and variance of $T$ stochastic forwarded passes for
the same input.

In case of MLE, we do not place distributions over parameters, they are
discrete like in case of deterministic network as in Eq. (8), optimal
parameters of the network are computed by minimising the following loss
function

$\mathcal{L}_{\text{MLE}}(\mathcal{D,\bm{\theta}})\approx-\sum_{k=1}^{N}\sum_{l=1}^{n}\left[-\log\left(\sqrt{2\pi}~{}\mathcal{U}^{(l)}_{\sigma}({\bf
f}^{(k)},\bm{\theta})\right)-\frac{\left(\hat{u}^{(k)}_{l}-\mathcal{U}^{(l)}_{\mu}({\bf
f}^{(k)},\bm{\theta})\right)^{2}}{2\left(\mathcal{U}^{(l)}_{\sigma}({\bf
f}^{(k)},\bm{\theta})\right)^{2}}\right],$ (18)

where

$\mathcal{U}^{(l)}_{\sigma}({\bf
f}^{(k)},\bm{\theta})=\log(1+\exp({\mathcal{U}^{(l)}_{\rho}({\bf
f}^{(k)},\bm{\theta}})).$ (19)

At last, optimal parameters of probabilistic U-Nets are computed by minimising
respective loss functions as

$\bm{\theta}_{\text{VB/MLE}}^{*}=\text{arg
}\underset{\bm{\theta}}{\text{min}}~{}\mathcal{L}_{\text{VB/MLE}}(\mathcal{D},\mbox{\boldmath{$\theta$}}).$
(20)

##  4 Results

###  4.1 The numerical experiment procedure

####  4.1.1 Generation of Training Data from Hyperelastic FEM Simulations

(a)

(b)

(c)

Figure 4: Data generation strategy for three benchmark examples (a) 2D beam,
(b) 2D L-shape and (c) 3D beam. The parts of top surfaces marked in red color
indicate the nodes at which random nodal forces were applied to generate
training datasets.

Two 2D and one 3D benchmark problems are considered in this work, as
schematically shown in Figure 4. The Neo-Hookean hyperelastic material model
is used, with the Young’s modulus $E=0.5~{}\text{KPa}$ and the Poisson’s ratio
$\nu=0.4$. We use the following version of Neo-Hookean strain energy potential

$W(\mbox{\boldmath{$F$}})=\frac{\mu}{2}(I_{c}-3-2\ln{J})+\frac{\lambda}{4}(J^{2}-1-2\ln{J}),$
(21)

where the invariants $J$ and $I_{\text{c}}$ are given in terms of deformation
gradient $F$ as

$J=\text{det}(\mbox{\boldmath{$F$}}),\quad
I_{\text{c}}=\text{tr}(\mbox{\boldmath{$F$}}^{T}\mbox{\boldmath{$F$}}),\quad\text{where}\quad\mbox{\boldmath{$F$}}=\mbox{\boldmath{$I$}}+\nabla\mbox{\boldmath{$u$}},$
(22)

while $\mu$ and $\lambda$ are Lame’s parameters, which can be expressed in
terms of the Young’s modulus, $E$, and the Poisson’s ratio, $\nu$, as

$\lambda=\frac{E\nu}{(1+\nu)(1-2\nu)},\quad\mu=\frac{E}{2(1+\nu)}.$ (23)

As introduced in Section 2.1, for a given discretized problem, the
training/testing dataset is constructed as follows. Within nodes occupying a
prescribed region of the boundary (in red color in Figure 4), a particular
family of external force distribution is considered. Each loading case consist
of a single excited node, while for remaining nodes the external forces are
${\bf 0}$. For a given training/testing example, a single node is chosen for
which the external force vector is generated randomly, component-wise, from a
uniform distribution within a given range of magnitude. The example is then
solved with FEM, and the entire vector ${\bf f}$ of prescribed nodal external
forces (including unloaded nodes) and the vector ${\bf u}$ of computed nodal
displacements are saved. The procedure is repeated for all $N+M$ examples,
which creates the training/testing dataset $D=\\{({\bf f}^{(1)},{\bf
u}^{(1)}),...,({\bf f}^{(N+M)},{\bf u}^{(N+M)})\\}$.

The finite element simulations have been performed with the AceGen/AceFem
framework (Korelc, 2002) (standard library displacement-based Neo-Hoohean
finite elements are used). The non-linear FE problems are solved with the
Newton-Raphson method, and adaptive load-stepping scheme is used to avoid
convergence issues for large load cases.A single quad/hexahedral FE mesh per
problem is only considered.

_Remark:_ For 2D/3D beam examples the structured FE mesh is used, which is
compatible with the U-Net architecture introduced in Section 2.2. In the
L-shaped example, the FE mesh is not structured, which makes it impossible to
directly transform it to a compatible node numbering, with a possible
consequence of accuracy drop, as explained in Section 4.2.1. To correct this,
a special zero-padding operation is applied to each ${\bf f}^{(i)}$ and ${\bf
u}^{(i)}$ before using the dataset for training/testing, see Fig. 5.

Figure 5: Extra nodes with zero force/displacement are added to convert the
data to a structured format. Node numbers are mapped accordingly to follow the
assumed order of U-Net architecture.

The datasets are randomly split into training sets, $N$ (95%), and testing
sets, $M$ (5%). The characteristics of FE meshes and datasets for all three
problems are provided in Table 1.

Problem | N.of DOFs | Force component range [N] | dataset size N+M
---|---|---|---
2D beam | 128 | -2.5 to 2.5 | $5700+300$
2D beam | 128 | -1 to 1 | $3800+200$
2D L-shape | 160 ($256^{*}$) | -1 to 1 | $3800+200$
3D beam | 12096 | -2 to 2 | $33688+1782$
Table 1: FE datasets. The number of DOFs with the asterisk refer to the zero-
padded 2D L-Shaped mesh.

####  4.1.2 Implementation and training of U-Nets

2D U-Net Architecture: We use 3 level U-Net architecture as in Figure 2 with
with c=128 channels in the first level, at each level we use apply two
convolutional operators with $3\times 3$ filters. For Bayesian U-Net case, we
replace half of the layers by probabilistic layers (one layer out of two at
each level is replaced by probabilistic layer).

3D U-Net Architecture: For 3D case we use 4 level U-Net architecture with c=64
channels in the first layer. In the encoding phase, spatial dimensions of 3D
mesh decrease by a factor of 2 and channel dimensions increase by 2. And in
decoding phase, it is vice-versa as similar in Figure 2. Again, at each level
we apply 2 convolutional operators with $3\times 3$ filters.

Additionally, we use batch-normalisation on each layer. Batch normalization is
a technique for training very deep neural networks that standardizes the
inputs to a layer for each mini-batch (Ioffe and Szegedy, 2015). This has the
effect of stabilizing the learning process and dramatically reducing the
number of training epochs required to train deep networks.

Training: Trainings of deterministic and probabilistic versions of U-Net are
carried out using Keras (Chollet et al., 2015) and Tensorflow-probability
(Dillon et al., 2017) libraries respectively. Network is trained by minimising
loss function for given training dataset, minimisation is performed using Adam
optimizer, a well-known adaptive stochastic gradient-descent algorithm. We set
the learning rate to $1\times 10^{-4}$ and set other optimizer parameters as
per recommendations (Kingma and Ba, 2017). For the Monte Carlo sampling of
loss function ($\mathcal{L}_{\text{VB}}$) in Eq. (15), we use Flipout
estimator (Wen et al., 2018) with its recommended parameter values (usually
around 50 Monte Carlo samples). All the implementations are done on Tesla
V100-SXM2 GPU, on HPC facilities of the University of Luxembourg (Varrette et
al., 2014) using batch size of 4 and 600/75 epochs for 2D/3D cases.

####  4.1.3 Validation Metrics for the testing phase

Test data $\\{({\bf f}_{1},{\bf u}_{1}),...,({\bf f}_{M},{\bf u}_{M})\\}$ is
generated in a similar way as we generate the training data as described in
Section 4.1.1. We use mean absolute error norm as the validation metric. For
$m^{th}$ test example, $\mathcal{U}({\bf f}_{m})$ is the deterministic network
prediction and ${\bf u}_{m}$ is the finite element solution. When $n$ is
number of dofs of the mesh, mean norm error is given as

$e_{m}=\frac{1}{n}\sum_{i=1}^{n}{|\mathcal{U}({\bf f}_{m})^{i}-{\bf
u}_{m}^{i}|}.$ (24)

To have a single validation metric over the entire test set, we compute the
average mean norm $\bar{e}$ and the corrected sample standard deviation
$\sigma(e)$ as follows (Note: It is standard deviation of averaged errors
across the test set, not the standard deviation of all errors.)

$\begin{split}\bar{e}=\frac{1}{M}\sum_{m=1}^{M}e_{m},\qquad\sigma(e)=\sqrt{\frac{1}{M-1}\sum_{m=1}^{M}\left(e_{m}-\bar{e}\right)^{2}}.\end{split}$
(25)

In case of Bayesian CNN, output of the network is a probability distribution,
so we sample over the output distribution by taking multiple forward passes as
described in Eq. (17). And mean of these samples
$\mathcal{U}_{\mu}(\bm{f}_{m})$, is taken as mean prediction of the Bayesian
CNN, and standard deviation of these samples
$\mathcal{U}_{\sigma}(\bm{f}_{m})$, gives us the confidence intervals of
predictions. Now the error norm for $m^{th}$ test example is given as

$e(\mathcal{U}_{\mu}({\bf f}_{m}),{\bf
u}_{m})=\frac{1}{n}\sum_{i=1}^{n}{|\mathcal{U}_{\mu}({\bf f}_{m})-{\bf
u}_{m}^{i}|}.$ (26)

Again the average error norm and the corrected sample standard deviation for
all test examples is computed as

$\begin{split}\bar{e}=\frac{1}{M}\sum_{i=1}^{M}e(\mathcal{U}_{\mu}({\bf
f}_{m}),{\bf
u_{m}}),\qquad\sigma(e)=\sqrt{\frac{1}{M-1}\sum_{m=1}^{M}\left(e(\mathcal{U}_{\mu}({\bf
f}_{m}),{\bf u}_{m})-\bar{e}\right)^{2}}.\end{split}$ (27)

###  4.2 Deterministic U-Nets

####  4.2.1 Advantages of the U-Net convolutional architecture

U-Nets vs. fully-connected NNs

CNNs take advantage of local spatial coherence in the input, enabling them to
have fewer weights as some parameters are shared. Nearby nodes of the FEM mesh
show strong local correlation and CNNs leverage this fact. If we had to use a
feed-forward nerual network, number of latent parameters of the network will
be far more in number as compared to CNN.

3D beam example in Figure 4c, it has 12096 dof, so consider a simplest feed-
forward network without any hidden layer and having 12096 units in output and
input layer respectively. Because of absence of hidden layer, this network
corresponds to the linear approximation of the force displacement
relationship. In order to have the best linearised approximation, we
initialise trainable parameters of feed forward network with the inverse of
the stiffness matrix. Table 2 shows that the simplest feed forward network
(which is an in-accurate assumption) takes more than $\times 1.5$ parameters
than the deterministic CNN. In order to do a better (a non-linear
approximation), one needs to use a multi layer feed forwarded network, which
will have much more number of parameters and hence proportionally training
time will be significantly high. Hence choice of CNN makes complete sense
particularly for finer mesh problems.

NN type | N. of trainable parameters |  $\bar{e}$ [m] |  $\sigma(e)$ [m]
---|---|---|---
Deterministic U-Net |  $~{}94.1$ E+6 |  $0.6$ E-3 |  $0.3$ E-3
Fully-connected |  $146.3$ E+6 |  $7.0$ E-3 |  $8.0$ E-3
Table 2: U-Net vs fully-connected feed forward network comparison for 3D Beam
example

Effect of DOF ordering

CNNs take advantage of local spatial coherence in the input, they force the
extraction of local features by using local convolutional kernels (Bengio and
Lecun, 1997). Topology of input plays a crucial role in CNNs, hence order of
node numbering of finite element mesh plays a very important role in training
procedure. Different FEM pre-processors have different ways of numbering
nodes. (Geuzaine and Remacle, 2009) is a popular FE mesh generator, it first
number corner nodes, then edge nodes followed by internal nodes. This strategy
is not optimal for our purpose. To fully leverage advantages of CNN, care has
to be taken to number nodes in spatial proximity one after another.

(a)

(b)

(c)

Figure 6: Different node numbering strategies. (a) Numbering assumed by the
TensorFlow (the preferred one; used in this work), (b) Gmsh preprocessor
numbering, and (c) random numbering.

Figure 6 shows different node ordering strategies. Proper ordering strategy as
in Figure 6a is essential for the efficient utilisation of convolutional
networks. For the 3D beam as shown in the above analysis, we train the network
with different ordering strategy and predictions are made on identical test
sets same as in above subsection. Table 3 shows that error metrics are nearly
9 times high when ordering is not optimal.

Ordering strategy |  $\bar{e}$ [mm] | $\sigma(e)~{}[\text{mm}]$
---|---|---
Preferred ordering as in Figure 6a | 0.6 | 0.3
Random ordering as in Figure 6c | 5.4 | 2.8
Table 3: Error metrics for randomly ordered dof are nearly 9 times more than
the optimal ordering

####  4.2.2 Prediction accuracy

We consider different datasets generated as per the description in Table 1.
All the examples have Young’s modulus $Y$ of 500 Pa and a Poisson’s ratio of
0.4. Table 4 gives validation error metrics for entire test sets of different
examples, where $M$ stands for number of test examples

Example | M |  $\bar{e}$ [m] | $\sigma(e)~{}[\text{m}]$ |  $t_{\text{CNN}}$ [ms] |  $t_{train}$ [min]
---|---|---|---|---|---
2D Beam | 300 | 0.3 E-3 | 0.2 E-3 | 1 | 131
2D L-Shaped | 200 | 0.8 E-3 | 0.4 E-3 | 1 | 64
3D Beam | 1782 | 0.6 E-3 | 0.3 E-3 | 9 | 1060
Table 4: Error metrics for 2D and 3D test sets for deterministic U-Net.

We perform sensitivity analysis of errors with increase in deformation
amplitude. In order to do so, we fit a linear line to the error point data
($e$) of all test examples. From Figure 7, all the three examples show a very
low sensitivity with increase of displacement magnitude. Slopes of linear fit
to 2D-beam, 2D-L shape and 3D beam cases are 0.0008,0.0005,0.001 respectively.

We take some of the test cases from each type of example and show the
comparison of FEM and U-Net solutions. We show the overlap of meshes, and also
compute node-wise $\ell_{2}$ norm of the error (each node has two/three dof
associated to it, we take $\ell_{2}$ norm for the error associated to each
node) and interpolate the resultant error field on deformed meshes.

Figure 7: We sort the error $e$ for test examples as per the increasing
displacement magnitude at the point of application of force. The regression
line of equations $y=0.0008\times x$ (2D-beam), $y=0.0005\times x$(2D-L
shape), $y=0.001\times x$ (3D beam) show low sensitivity of the deterministic
U-Nets to displacement magnitudes

(a)  (b)

Figure 8: (a) Deformation of 2D beam discretised with $16\times 4$ mesh is
computed using deterministic U-Net. FEM reference solution is in red, U-Net
solution is in blue (both are overlapping) gray represents the rest position
of the mesh. Magnitude of tip deformation is 0.95 m and relative $\ell_{2}$
norm is 0.6%. (b) For the same example, we take $\ell_{2}$ norm of error at
each node. Error counter is obtained by interpolating nodal error values.

(a)  (b)

Figure 9: (a) Deformation of 2D-L shape mesh is computed using deterministic
U-Net. FEM reference solution is in red, U-Net solution is in blue (both are
overlapping), gray represents the rest position of the mesh. Magnitude of tip
deformation is 2.39 m and relative $\ell_{2}$ norm of FEM and deterministic
CNN is 0.4%. (b) For the same test example as in Figure 9, we take $\ell_{2}$
norm of error at each node. Error counter is obtained by interpolating nodal
error values.

(a)

(b)

Figure 10: (a) Deformation of 3D beam discretised with $28\times 12\times 12$
mesh is computed using deterministic U-Net. FEM reference solution is in red,
U-Net solution is in blue, gray represents the rest position of the mesh.
Magnitude of tip deformation is 1.1m and relative $\ell_{2}$ norm is 0.6% (b)
For the same test example as in Figure 10, we take $\ell_{2}$ norm of error at
each node. Error counter is obtained by interpolating nodal error values.

(a)

(b)

Figure 11: (a) Another example for same 3D Beam example as in Figure 10 FEM
reference solution is in red, U-Net solution is in blue, gray represents the
rest position of the mesh. Magnitude of tip deformation is 0.26m and relative
$\ell_{2}$ norm is 1.6% (b) For the same test example as in Figure 11, we take
$\ell_{2}$ norm of error at each node. Error counter is obtained by
interpolating nodal error values. One can observe that error at the point of
load application is higher.

Effect of changing the distribution of applied forces

Deterministic network has been trained by using single point load examples
only. But in the prediction phase we apply multiple force, Figure 12 shows one
such example where random multiple forces are applied on the top edge, U-Net
is able to closely follow the FEM solution. Solution accuracy of multiple
point load cases can be improved by incorporating multiple point loads in the
training data.

(a)  (b)

Figure 12: (a) Multiple point loads are applied in the prediction phase. FEM
reference solution is in red, CNN solution is in blue gray represents the rest
position of the mesh. Magnitude of tip deformation is 0.55 m and relative
$\ell_{2}$ norm of FEM and NN error is 0.6%. (b) For the same test example as
in Figure 12, we take $\ell_{2}$ norm of error at each node. Error counter is
obtained by interpolating nodal error values.

Also when we do predictions for the point loads applied on nodes as we move in
the perpendicular direction of training zone, error values increase. Consider
the 2D-beam example. Training has been done by applying random point loads on
the top line. To illustrate a small example, we apply a vertical force of 1.5N
on each of the 4 nodes of the free edge of the beam. Figure 13 shows as we
move away from the training zone, the mean error($e$) and maximum error
observed for the example increases. Figure 13 also shows uncertainty of
prediction for the dof on which point load is applied. As one can infer,
Bayesian U-Net is not giving reliable uncertainty estimates, as we move away
from training dataset in spatial directions, this is one of the limitations of
the method.

(a)  (b)

Figure 13: Mean error($e$) and maximum error are increasing for 4 point loads
as we move away in the perpendicular direction from the training line. Green
stands for prediction uncertainty of node on which force is applied. Bayesian
U-Net is not giving reliable uncertainty estimates as we move away in
perpendicular direction.

(a)  (b)

Figure 14: For the case in Figure 13, F=1.5 N is applied (a) on the
$2^{\text{nd}}$ node as (b) On the $4^{\text{th}}$ node. FEM solution is in
red, U-Net solution is in blue, as we move away from training line, U-Net
predictions show more deviation from true FEM solution.

###  4.3 Probabilistic U-Nets

The most important thing for probabilistic U-Nets to work is use of Emperical
Bayes as described in Section 3.3, i.e. inclusion of prior distribution
parameters in the training procedure.

####  4.3.1 Prediction accuracy

2D Beam: First we consider the same dataset as used for training deterministic
U-Net in Section 4.2.2 (N=5700, M=300). The training data-set was created by
applying point loads on the top edge with magnitude of x,y dof of force is not
more than 2.5N (For all point loads, $|F_{x}|,|F_{y}|\leq 2.5$). Ultimate goal
of our Bayesian CNN framework is to have a good mean fit and at the same time
give the confidence intervals of predictions. Output of the network being
distribution, for the trained network, we make 300 stochastic forward passes
to get the statistics of the output. Error metrics for the for test data over
Bayesian CNN are given in Table 5 and Figure 15 gives the sensitivity of
errors ($e(\mathcal{U}_{\mu}({\bf f}_{m}),{\bf u}_{m})$) with increasing
displacement magnitudes. Overall Bayesian U-Nets tend to give higher error
values, a behaviour which is expected because of introduction of
distributions, which regularises the model (Vladimirova et al., 2019).

Example | $M$ | $\bar{e}~{}[\text{m}]$ | $\sigma(e)~{}[\text{m}]$ | $t_{\text{pred}}~{}[\text{s}]$ | $t_{\text{train}}~{}[\text{min}]$
---|---|---|---|---|---
2D Beam (VB) | 300 | 1.3 E-3 | 1.3 E-3 | 0.1 | 226
2D Beam (D) | 300 | 0.3 E-3 | 0.2 E-3 | 0.001 | 131
2D L-Shaped (VB) | 200 | 5.3 E-3 | 3.7 E-3 | 0.1 | 143
2D L-Shaped (D) | 200 | 0.8 E-3 | 0.4 E-3 | 0.001 | 64
Table 5: Error metrics for 2D test sets using Bayesian U-Nets Figure 15: We
sort the error $e$ for test examples as per the increasing displacement
magnitudes at the point of application of force. slope for beam = 0.005,
l-shape = 0.006 The regression line of equation $y=0.005\times x$ (2D beam),
$y=0.006\times x$ shows the low sensitivity of Bayesian U-Net to displacement
magnitude

We consider the same point load test case as in Figure 8, which we considered
in deterministic case. This time we make the prediction using Bayesian U-Net.
Figure 16 shows the comparison of FEM and NN error and uncertainty associated
with the prediction (we plot single standard deviation values associated to
the prediction of respective dof). One can see that both are strongly co-
related spatially.

(a)

(b)

Figure 16: Error and uncertainty are co-related, (a) gives error between mean
prediction of Bayesian CNN and FEM. (b) gives uncertainty associated with the
prediction.

Similar kind of analysis is done for the multiple point load case (Figure. 12)
which is showed in the deterministic case. We plot the nodal errors and
uncertainties (single standard deviation values), we can see that they closely
follow each other.

(a)

(b)

Figure 17: Error and uncertainty are co-related, (a) gives error between mean
prediction of Bayesian CNN and FEM. (b) gives uncertainty associated with the
prediction.

Out of distribution case: Let us consider a test example in which a force of 5
N is applied on the corner node, which is far away from the training range
(-2.5 to 2.5 N). Again we compare error and uncertainty associated with the
Bayesian U-Net prediction. In Figure 18, one can see that both are strongly
corelated, rather both values are close to each other across the spatial
dimensions of beam. Thus the uncertainty predictions can give us an idea about
the error of U-Net predictions, thus making it more reliable. Almost all the
times, two standard deviation is able to capture the error.

(a)

(b)

Figure 18: Error and uncertainty are co-related, (a) gives error between mean
prediction of Bayesian U-Net and FEM. (b) gives uncertainty associated with
the prediction.

2D Beam: We focus displacement prediction of single dof (y-direction dof of
the last node). We apply multiple vertical forces varying from -8N to 8N on
corner node of the beam and predict its displacements using the Bayesian
U-Net. Figure 19a gives the prediction of displacement magnitude, as one can
see prediction matches with test FEM solution within the training region,
outside the training region, Bayesian U-Net prediction deviates from the FEM
solution. For the reference, we provide solutions for same cases using
deterministic U-Net as well, they deviate from FEM solution outside training
range. Figure 19b gives confidence intervals associated with the prediction.
One can see that network has very less uncertainty i.e. it is confident in the
regions of training data (-2.5 to 2.5 N) but as one moves away, uncertainty of
the prediction increases. We can also see that 95% confidence is able to
capture the error of Bayesian U-Net predictions outside the training region,
for the reference errors of deterministic U-Nets are presented as well.

(a)

(b)

Figure 19: Range of vertical forces applied on the corner node, (19a) gives
prediction(mean prediction) using Bayesian CNN. (19b) gives uncertainty of
network, as one can see, as force values move out of training set, network is
more uncertain in its prediction. Also, 95% confidence interval is capturing
errors of bayesian as well as deterministic U-Net outside the training region

2D Beam : This time we consider second training dataset for 2D beam with
forces less than 1N on each dof. (For all point loads in training set,
$|F_{x}|,|F_{y}|\leq 1\text{N}$). Figure 20 shows the comparisons of
uncertainty intervals for different training datasets. As the range of force
in training set is increasing, uncertainty starts decreasing.

Figure 20: Uncertainty intervals different input force range. As the range of
forces in training data increases, uncertainty is reduced.

2D L-shape: We consider the same data-set as used in deterministic case (N =
4000, M = 200). The training data-set was created by applying point loads on
the top edge with magnitude of x,y dof of force is not more than 1N (For all
point loads, $|F_{x}|,|F_{y}|\leq 1$). As similar to the beam case, statistics
of the output it obtained by doing 300 stochastic passes on trained network
for the same input.

2D L-shape: We apply multiple forces in horizontal direction varying from -6N
to 6N on the inner corner of the L-shape. Figure 21a shows how non-linear
deformation changes with the force value applied. Figure 21b gives the
uncertainty associated with the prediction. As the force value goes outside
the training value range, uncertainty increases.

(a)

(b)

Figure 21: 2D L-Shape. Force applied on the in corner node, (21a) gives mean
prediction using Bayesian CNN. We have also plotted some of the test points
whose deformation is calculated using FEM. (21b) gives uncertainty of network,
as one can see, as force values move out of training set, network is more
uncertain in its prediction.

####  4.3.2 Data vs. model uncertainties

Use of Bayesian approach enables us to give credible intervals of the
predictions. As one can see in Figure 19b and Figure 20, uncertainty shows an
increasing trend in the extrapolated region. Non-Bayesian approach, such as
MLE fail to show that. MLE approach is capable of detecting the noises in the
data region, but they don’t give any additional information about the
extrapolated region. Let us consider the same dataset with force magnitude is
less than 1 N, now we randomly incorporate 20% noise when magnitude of applied
force is less than 0.7 N (For all point loads, $|F_{x}|,|F_{y}|\leq
0.7\text{N}$). Figure 22 shows uncertainty is only capturing noises in the
data.

(a)

(b)

Figure 22: 2D Beam. MLE approach to detect noises in the data. We incoporated
20% noise for the forces with magnitude less 0.7 N. As one can see in (a), we
are able to capture the noise in the data with MLE, but it fails to quantify
the uncertainty in extrapolated region, network is overconfident in
extrapolated region. Where as in (b) we train the model with Variational
Bayes, where network captures the data noise and uncertainty grows quickly in
extrapolated region as well.

In a similar way as in case of the beam example, we implement a MLE approach
on noisy data. Let us consider the same dataset, now we randomly incorporate
20% noise in the dof deformations when magnitude of force is less than 0.7 N
(For all point loads, $|F_{x}|,|F_{y}|\leq 0.7\text{N}$). Figure 23 shows
uncertainty is more only in the region where data has noise in it, in the
extrapolated region the network is more confident. Hence MLE approach is only
capable of capturing noises in the data, fails to quantify the uncertainty in
the extrapolated region.

(a)

(b)

Figure 23: 2D L-Shaped. MLE approach to detect noises in the data. We
incoporated 20% noise for the forces with magnitude less 0.7 N. As one can see
in (a), we are able to capture the noise in the data with MLE, but it fails to
quantify the uncertainty in extrapolated region,network is overconfident in
extrapolated region. Where as in (b) we train the model with Variational
Bayes, where network captures the data noise and uncertainty grows quickly in
extrapolated region as well.

###  4.4 Training and prediction times

Training Times

For any neural network, the training phase of the deep learning model is the
most resource-intensive task. Hence modern machine learning open source
libraries such as Tensorflow, Keras, PyTorch are optimized to work with
Graphical Processing Units (GPU). Compared to the user end computer with its
Central Processing Unit (CPU), a GPU has a parallel structure that offers
faster computing and increased efficiency. Table 6 gives GPU training times
for different datasets for both determinsitic and probabilistic U-Nets.

Example | $N$ | $t_{\text{train}}~{}[\text{min}]$
---|---|---
2D Beam (D) | 5700 | 131
2D Beam (VB) | 5700 | 226
2D L-Shaped (D) | 3800 | 64
2D L-Shaped (VB) | 3800 | 143
3D Beam (D) | 33688 | 1060
Table 6: Training times for U-Nets. D= Determinitic, VB = Variational Bayes

Consider 2D-Beam example, training time for Bayesian U-Net is nearly 1.75
times of deterministic network Table 6. This is because we used half of the
layers as probabilistic layers, doubling the number of parameters for them,
and we also included prior means in the training for those layers. Hence
number of trainable parameters for Bayesian U-Net is nearly 1.75 times of its
deterministic counterpart. Similar explanation can be given for the 2D-L shape
case.

Training time for 3D beam case is high, it can be reduced by opting alternate
topologies of U-Net. One way of doing so is keeping the constant number of
channels in each U-Net level.

Effect of number of channels

Training time of U-Net can be reduced by decreasing the number of parameters
of model, this can be achieved by decreasing number of channels in each stage.
We use a architecture with constant number of channels in level of U-Net
architecture. Table 7 shows a comparison of training time and error metrics
for different architectures for the 3D-Beam dataset. Number of parameters
increase with number of channels. For 128 channels one can notice that test
error is more as compared to 64 channels, indicating over-fitting of the
network.

N. of channels, $c$ |  $\bar{e}$ [m] | $\sigma(e)~{}[\text{m}]$ | $t_{\text{train}}~{}[\text{min}]$
---|---|---|---
16 | 1.6 E-3 | 0.9 E-3 | 272
32 | 1.1 E-3 | 0.7 E-3 | 293
64 | 0.8 E-3 | 0.5 E-3 | 348
128 | 3.5 E-3 | 3.3 E-3 | 646
Table 7: Constant number of channels ’c’ are used in each level of U-Net. Use
of 64 channels is optimum to achieve error and computational time trade off

Prediction Times

Since recent years, GPU cloud computing is getting easily accessible, one can
leverage GPU support over the Internet as well. All these factors make our
framework easily deployable to the user end. Yet some readers might be
interested in knowing the prediction time comparison on standard CPU instead
of GPU.

Type | dof |  $t\\_fem_{\text{CPU}}$ [s] |  $t_{\text{CPU}}$ [s] |  $t_{\text{GPU}}$ [s] | $\frac{t\\_\text{fem}_{\text{CPU}}}{t_{\text{CPU}}}$ | $\frac{t\\_fem_{\text{CPU}}}{t_{\text{GPU}}}$
---|---|---|---|---|---|---
2D Beam | 128 | 0.123 | 0.005 | 0.001 | 25 | 123
2D L-shape | 256 | 0.120 | 0.007 | 0.001 | 17 | 120
3D Beam | 12096 | 3.1 | 0.1 | 0.009 | 31 | 345
Table 8: Prediction times of deterministic U-Net on CPU, GPU. Under similar
computational resources on CPU, U-Net shows x30 speedup, which can be even
more depending on the boundary conditions. Where as GPU shows nearly x350
speedup.

Table 8 gives the comparison of prediction times for different examples. For
some of the force values in the training set, FEM solution took more than
$3$s. Hence under the identical computational resources, deterministic CNN
gave x30 speed up for given dataset. Another important point to mention here
is, both deterministic and Bayesian U-Net take almost same time for prediction
irrespective of the value of the applied force. Where as in case of FEM
solution, time changes with the value of force, because we use an iterative
solver and adaptive load-stepping scheme to avoid convergence issues for large
load cases. Hence even on the local, we can expect much more speedup than x30
if we go towards the higher range forces.

3D Beam: Prediction time for test examples using deterministic U-Net is only
$9$ ms. For the finite element solution, we use one of the very efficient
solvers, MKL Pardiso solver (Schenk and Gärtner, 2004). Some of the test
examples with $\ell_{2}$ norm of nodal force greater than 3N took as much as
3.1 secs. Thereby deterministic U-Net is observed to give nearly x350 speedup
on GPU.

Time of prediction of Bayesian U-Net is time of sampling over output
distribution, which is time of 300 stochastic forward passes in our case. In
case of above example, 300 forward passes for a single test example took 0.1
secs for both 2D beam and 2D-L shape examples. Compared to deterministic case,
average time for the prediction of single pass is less (0.3 ms) because of
efficient utilisation of batch prediction. Hence even Bayesian inference takes
a very little time in the prediction phase.

##  5 Conclusion

In this work we have proposed a Bayesian CNN framework which is capable of
accurately predicting non-linear deformations of hyper-elastic bodies in real
time. With the use of CNNs, we observed nearly x350 speedup in computational
time, thus enabling goal of real time accurate simulations. This has a
potential for problems requiring very fast simulations of objects undergoing
interactions, such as surgical simulators.

We gained a detailed understanding of efficient utilisation of U-Nets for real
time hyper-elastic simulations. We elaborated important factors such as node
number ordering, use of Empirical Bayes, thus allowing us to use state of the
art data driven methods for the solid mechanics applications.

We also proposed a probabilistic framework which is capable of giving reliable
uncertainty estimates of U-Net models in addition to the predictions. We
showed the correlation of U-Net errors (deviation from FEM solution) with
prediction uncertainty of Bayesian U-Nets. We also showed that uncertainty of
prediction rapidly increases in the extrapolated region. With our framework we
could give uncertainty of U-Net models as well as capture the uncertainty
present in the data. Thereby our framework goes a step further, to make real
time large deformation simulations more trustworthy, which is crucial in
engineering applications.

Besides showing promising results, we also demonstrated several important
limitations of the current framework. Firstly, the convolution operations that
are used in U-Nets strongly rely on underlying structured meshes. With
moderate effort this limitation can be possibly alleviated, e.g., by creating
a mapping between a given arbitrary unstructured mesh and a suitable
structured grid. Another way of generalising framework for irregular meshes is
to use newly developed graph convolutions (Hanocka et al., 2019), learning on
complex shapes and geometries is relatively easier with graph convolution
operations. Secondly, the current framework performs poorly in prediction of
multiple point loads, especially when point loads are applied on adjacent
nodes. Mechanical coupling can be one of the reasons for this, and one way to
alleviate this problem is to include multiple point load examples in the
training phase. However, this will obviously further increase the
computational cost of training. Thirdly, we observed that the proposed novel
technique to quantify uncertainties in extrapolated regions does not always
give reliable predictions. Bayeisan U-Nets failed to give reliable credible
intervals of predictions when we applied the force on the node which is
outside of training region (in spatial dimensions), which will be addressed in
the future work. Also they didn’t perform efficiently when data was missing in
the interpolated region, reliable uncertainty estimates in the missing
interpolated region can be obtained by using noise contrastive priors (NCPs)
(Hafner et al., 2020). We leave these challenging questions for future
research.

_Acknowledgements:_

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Sklodowska-Curie grant
agreement No. 764644. Jakub Lengiewicz would like to acknowledge the support
from EU Horizon 2020 Marie Sklodowska Curie Individual Fellowship _MOrPhEM_
under Grant 800150. This paper only contains the author’s views and the
Research Executive Agency and the Commission are not responsible for any use
that may be made of the information it contains.

## References

  * Bengio and Lecun (1997) Y. Bengio and Yann Lecun.  Convolutional networks for images, speech, and time-series.  _The handbook of brain theory and neural networks_ , 11 1997.
  * Blundell et al. (2015) Charles Blundell, Julien Cornebise, Koray Kavukcuoglu, and Daan Wierstra.  Weight uncertainty in neural network.  In Francis Bach and David Blei, editors, _Proceedings of the 32nd International Conference on Machine Learning_ , volume 37 of _Proceedings of Machine Learning Research_ , pages 1613–1622, Lille, France, 07–09 Jul 2015. PMLR.
  * Bui et al. (2018) Huu Phuoc Bui, Satyendra Tomar, Hadrien Courtecuisse, Stéphane Cotin, and Stéphane P. A. Bordas.  Real-time error control for surgical simulation.  _IEEE Transactions on Biomedical Engineering_ , 65(3):596–607, 2018.  doi: 10.1109/TBME.2017.2695587.
  * Carlin and Louis (1997) Bradley P. Carlin and Thomas A. Louis.  Bayes and empirical bayes methods for data analysis.  _Statistics and Computing_ , 1997.  doi: https://doi.org/10.1023/A:1018577817064.
  * Chollet et al. (2015) François Chollet et al.  Keras.  https://github.com/fchollet/keras, 2015.
  * Cotin et al. (1999) S. Cotin, H. Delingette, and N. Ayache.  Real-time elastic deformations of soft tissues for surgery simulation.  _IEEE Transactions on Visualization and Computer Graphics_ , 5(1):62–73, 1999.  doi: 10.1109/2945.764872.
  * Courtecuisse et al. (2014) Hadrien Courtecuisse, Jérémie Allard, Pierre Kerfriden, Stéphane P.A. Bordas, Stéphane Cotin, and Christian Duriez.  Real-time simulation of contact and cutting of heterogeneous soft-tissues.  _Medical Image Analysis_ , 18(2):394–410, 2014\.  ISSN 1361-8415.  doi: https://doi.org/10.1016/j.media.2013.11.001.  URL https://www.sciencedirect.com/science/article/pii/S1361841513001692.
  * Delingette et al. (1999) H. Delingette, S. Cotin, and N. Ayache.  A hybrid elastic model allowing real-time cutting, deformations and force-feedback for surgery training and simulation.  In _Proceedings Computer Animation 1999_ , pages 70–81, 1999.  doi: 10.1109/CA.1999.781200.
  * Dillon et al. (2017) Joshua V. Dillon, Ian Langmore, Dustin Tran, Eugene Brevdo, Srinivas Vasudevan, Dave Moore, Brian Patton, Alex Alemi, Matthew D. Hoffman, and Rif A. Saurous.  Tensorflow distributions.  _CoRR_ , abs/1711.10604, 2017.  URL http://arxiv.org/abs/1711.10604.
  * Duerr et al. (2020) O. Duerr, B. Sick, and E. Murina.  _Probabilistic Deep Learning: With Python, Keras and TensorFlow Probability_.  Manning Publications, 2020.  ISBN 9781617296079.  URL https://books.google.co.in/books?id=-bYCEAAAQBAJ.
  * Farhat and Roux (1991) Charbel Farhat and Francois-Xavier Roux.  A method of finite element tearing and interconnecting and its parallel solution algorithm.  _International Journal for Numerical Methods in Engineering_ , 32(6):1205–1227, 1991.  doi: https://doi.org/10.1002/nme.1620320604.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620320604.
  * Gal (2016) Yarin Gal.  Uncertainty in deep learning.  2016\.
  * Geuzaine and Remacle (2009) Christophe Geuzaine and Jean-François Remacle.  Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities.  _International Journal for Numerical Methods in Engineering_ , 79:1309 – 1331, 09 2009.  doi: 10.1002/nme.2579.
  * Goodfellow et al. (2016) Ian J. Goodfellow, Yoshua Bengio, and Aaron Courville.  _Deep Learning_.  MIT Press, Cambridge, MA, USA, 2016.  http://www.deeplearningbook.org.
  * Goury and Duriez (2018) Olivier Goury and Christian Duriez.  Fast, generic and reliable control and simulation of soft robots using model order reduction.  _IEEE Transactions on Robotics_ , 34(6):1565 – 1576, 2018.  doi: 10.1109/TRO.2018.2861900.  URL https://hal.archives-ouvertes.fr/hal-01834483.
  * Graves (2011) Alex Graves.  Practical variational inference for neural networks.  In J. Shawe-Taylor, R. Zemel, P. Bartlett, F. Pereira, and K. Q. Weinberger, editors, _Advances in Neural Information Processing Systems_ , volume 24. Curran Associates, Inc., 2011.  URL https://proceedings.neurips.cc/paper/2011/file/7eb3c8be3d411e8ebfab08eba5f49632-Paper.pdf.
  * Hafner et al. (2020) Danijar Hafner, Dustin Tran, Timothy Lillicrap, Alex Irpan, and James Davidson.  Noise contrastive priors for functional uncertainty.  In Ryan P. Adams and Vibhav Gogate, editors, _Proceedings of The 35th Uncertainty in Artificial Intelligence Conference_ , volume 115 of _Proceedings of Machine Learning Research_ , pages 905–914. PMLR, 22–25 Jul 2020.  URL https://proceedings.mlr.press/v115/hafner20a.html.
  * Hanocka et al. (2019) Rana Hanocka, Amir Hertz, Noa Fish, Raja Giryes, Shachar Fleishman, and Daniel Cohen-Or.  Meshcnn: A network with an edge.  38(4), 2019.  ISSN 0730-0301.  doi: 10.1145/3306346.3322959.  URL https://doi.org/10.1145/3306346.3322959.
  * Hauseux et al. (2018) Paul Hauseux, Jack S.Hale, Stéphane Cotin, and Stéphane P.A.Bordas.  Quantifying the uncertainty in a hyperelastic soft tissue model with stochastic parameters.  _Applied Mathematical Modelling_ , 62:86–102, 2018.  ISSN 0307-904X.  doi: https://doi.org/10.1016/j.apm.2018.04.021.  URL https://www.sciencedirect.com/science/article/pii/S0307904X18302063.
  * He and Xu (2019) Juncai He and Jinchao Xu.  Mgnet: A unified framework of multigrid and convolutional neural network.  _Science china mathematics_ , 62(7):1331–1354, 2019.
  * Ioffe and Szegedy (2015) Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  02 2015.
  * Kang et al. (2000) Young-Min Kang, Jeong-Hyeon Choi, Hwan-Gue Cho, and Chan-Jong Park.  Fast and stable animation of cloth with an approximated implicit method.  In _Proceedings Computer Graphics International 2000_ , pages 247–255, 2000.  doi: 10.1109/CGI.2000.852340.
  * Kendall and Gal (2017) Alex Kendall and Yarin Gal.  What uncertainties do we need in bayesian deep learning for computer vision?  In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , NIPS’17, page 5580–5590, Red Hook, NY, USA, 2017. Curran Associates Inc.  ISBN 9781510860964.
  * Kerfriden et al. (2011) P. Kerfriden, P. Gosselet, S. Adhikari, and S.P.A. Bordas.  Bridging proper orthogonal decomposition methods and augmented newton–krylov algorithms: An adaptive model order reduction for highly nonlinear mechanical problems.  _Computer Methods in Applied Mechanics and Engineering_ , 200(5):850–866, 2011.  ISSN 0045-7825.  doi: https://doi.org/10.1016/j.cma.2010.10.009.  URL https://www.sciencedirect.com/science/article/pii/S0045782510002872.
  * Kerfriden et al. (2012) Pierre Kerfriden, Olivier Goury, Timon Rabczuk, and Stéphane Bordas.  A partitioned model order reduction approach to rationalise computational expenses in multiscale fracture mechanics.  _Computer Methods Appl Mech Eng_ , pages 169–188, 08 2012.  doi: https://doi.org/10.1016/j.cma.2012.12.004.
  * Kingma and Ba (2017) Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2017.
  * Kingma et al. (2015) Durk P Kingma, Tim Salimans, and Max Welling.  Variational dropout and the local reparameterization trick.  In C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 28. Curran Associates, Inc., 2015.  URL https://proceedings.neurips.cc/paper/2015/file/bc7316929fe1545bf0b98d114ee3ecb8-Paper.pdf.
  * Korelc (2002) Joze Korelc.  Multi-language and multi-environment generation of nonlinear finite element codes.  _Engineering with Computers_ , 18:312–327, 01 2002.  doi: 10.1007/s003660200028.
  * Lagaris et al. (1998) I.E. Lagaris, A. Likas, and D.I. Fotiadis.  Artificial neural networks for solving ordinary and partial differential equations.  _IEEE Transactions on Neural Networks_ , 9(5):987–1000, 1998.  doi: 10.1109/72.712178.
  * Lagaris et al. (2000) I.E. Lagaris, A.C. Likas, and D.G. Papageorgiou.  Neural-network methods for boundary value problems with irregular boundaries.  _IEEE Transactions on Neural Networks_ , 11(5):1041–1049, 2000.  doi: 10.1109/72.870037.
  * Langer U. (2003) Steinbach O Langer U.  Boundary element tearing and interconnecting methods.  _Computing_ , 71:205–228, 2003.  doi: https://doi.org/10.1007/s00607-003-0018-2.
  * Lorente et al. (2017) D. Lorente, F. Martínez-Martínez, M.J. Rupérez, M.A. Lago, M. Martínez-Sober, P. Escandell-Montero, J.M. Martínez-Martínez, S. Martínez-Sanchis, A.J. Serrano-López, C. Monserrat, and J.D. Martín-Guerrero.  A framework for modelling the biomechanical behaviour of the human liver during breathing in real time using machine learning.  _Expert Systems with Applications_ , 71:342–357, 2017.  ISSN 0957-4174.  doi: https://doi.org/10.1016/j.eswa.2016.11.037.  URL https://www.sciencedirect.com/science/article/pii/S0957417416306728.
  * Marinkovic and Zehn (2019) Dragan Marinkovic and Manfred Zehn.  Survey of finite element method-based real-time simulations.  _Applied Sciences_ , 9(14), 2019.  ISSN 2076-3417.  doi: 10.3390/app9142775.  URL https://www.mdpi.com/2076-3417/9/14/2775.
  * McAllister et al. (2017) Rowan McAllister, Yarin Gal, Alex Kendall, Mark van der Wilk, Amar Shah, Roberto Cipolla, and Adrian Weller.  Concrete problems for autonomous vehicle safety: Advantages of bayesian deep learning.  In _Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence, IJCAI-17_ , pages 4745–4753, 2017.  doi: 10.24963/ijcai.2017/661.  URL https://doi.org/10.24963/ijcai.2017/661.
  * McFall and Mahan (2009) Kevin Stanley McFall and James Robert Mahan.  Artificial neural network method for solution of boundary value problems with exact satisfaction of arbitrary boundary conditions.  _IEEE transactions on neural networks_ , 20(8):1221—1233, August 2009.  ISSN 1045-9227.  doi: 10.1109/tnn.2009.2020735.  URL https://doi.org/10.1109/TNN.2009.2020735.
  * Mendizabal et al. (2019) Andrea Mendizabal, Pablo Márquez-Neila, and Stéphane Cotin.  Simulation of hyperelastic materials in real-time using deep learning.  _Medical Image Analysis_ , 59:101569, 10 2019.  doi: 10.1016/j.media.2019.101569.
  * Neal (1996) Radford M. Neal.  _Bayesian Learning for Neural Networks_.  Springer-Verlag, Berlin, Heidelberg, 1996.  ISBN 0387947248.
  * Nedel and Thalmann (1998) L. Nedel and D. Thalmann.  Real time muscle deformations using mass-spring systems.  _Proceedings. Computer Graphics International (Cat. No.98EX149)_ , pages 156–165, 1998.
  * Patete et al. (2013) Paolo Patete, Maria Ida Iacono, Maria Francesca Spadea, Giovanna Trecate, Daniele Vergnaghi, Luca Tommaso Mainardi, and Guido Baroni.  A multi-tissue mass-spring model for computer assisted breast surgery.  _Medical Engineering & Physics_, 35(1):47–53, 2013.  ISSN 1350-4533.  doi: https://doi.org/10.1016/j.medengphy.2012.03.008.  URL https://www.sciencedirect.com/science/article/pii/S1350453312000562.
  * Raissi et al. (2019) 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.  _Journal of Computational Physics_ , 378:686–707, 2019\.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2018.10.045.  URL https://www.sciencedirect.com/science/article/pii/S0021999118307125.
  * Rappel et al. (2020) H Rappel, L Beex, J S Hale, L Noels, and S P A Bordas.  A tutorial on bayesian inference to identify material parameters in solid mechanics.  _Archives of Computational Methods in Engineering_ , 27, 2020.  doi: https://doi.org/10.1007/s11831-018-09311-x.
  * Ronneberger et al. (2015) O. Ronneberger, P.Fischer, and T. Brox.  U-net: Convolutional networks for biomedical image segmentation.  In _Medical Image Computing and Computer-Assisted Intervention (MICCAI)_ , volume 9351 of _LNCS_ , pages 234–241. Springer, 2015.  URL http://lmb.informatik.uni-freiburg.de/Publications/2015/RFB15a.  (available on arXiv:1505.04597 [cs.CV]).
  * Schenk and Gärtner (2004) Olaf Schenk and Klaus Gärtner.  Solving unsymmetric sparse systems of linear equations with pardiso.  _Future Generation Computer Systems_ , 20(3):475–487, 2004.  ISSN 0167-739X.  doi: https://doi.org/10.1016/j.future.2003.07.011.  URL https://www.sciencedirect.com/science/article/pii/S0167739X03001882.  Selected numerical algorithms.
  * Scott (2008) Susanne C. BrennerL. Ridgway Scott.  _Finite Element Multigrid Methods_ , pages 155–173.  Springer New York, New York, NY, 2008.  ISBN 978-0-387-75934-0.  doi: 10.1007/978-0-387-75934-0˙7.  URL https://doi.org/10.1007/978-0-387-75934-0_7.
  * Varrette et al. (2014) Sébastien Varrette, Pascal Bouvry, Hyacinthe Cartiaux, and Fotis Georgatos.  Management of an academic hpc cluster: The ul experience.  07 2014.  doi: 10.1109/HPCSim.2014.6903792.
  * Vladimirova et al. (2019) Mariia Vladimirova, Jakob Verbeek, Pablo Mesejo, and Julyan Arbel.  Understanding priors in Bayesian neural networks at the unit level.  In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , pages 6458–6467. PMLR, 09–15 Jun 2019.  URL http://proceedings.mlr.press/v97/vladimirova19a.html.
  * Wang et al. (2020) Feng Wang, Alberto Eljarrat, Johannes Müller, Trond Henninen, Rolf Erni, and Christoph Koch.  Multi-resolution convolutional neural networks for inverse problems.  _Scientific Reports_ , 10:5730, 03 2020.  doi: 10.1038/s41598-020-62484-z.
  * Wen et al. (2018) Yeming Wen, Paul Vicol, Jimmy Ba, Dustin Tran, and Roger Grosse.  Flipout: Efficient pseudo-independent weight perturbations on mini-batches, 2018.
  * Wu et al. (2015) Jun Wu, Rüdiger Westermann, and Christian Dick.  A survey of physically based simulation of cuts in deformable bodies.  _Computer Graphics Forum_ , 34, 03 2015.  doi: 10.1111/cgf.12528.
  * Zhang et al. (2001) Yu Zhang, E.C. Prakash, and E. Sung.  Real-time physically-based facial expression animation using mass-spring system.  In _Proceedings. Computer Graphics International 2001_ , pages 347–350, 2001.  doi: 10.1109/CGI.2001.934696.
  * Zienkiewicz and Taylor (1991) OC Zienkiewicz and RL Taylor.  The finite element method, ; volume 2: solid and fluid mechanics, dynamics and non-linearity.  1991\.


# Cloud Removal from Satellite Images

Rutvik Chauhan  Computing Science(Multimedia)
University of Alberta
Edmonton, Canada
<EMAIL_ADDRESS>Antarpuneet Singh  Computing Science(Multimedia)
University of Alberta
Edmonton, Canada
<EMAIL_ADDRESS>Sujoy Saha  Computing Science(Multimedia)
University of Alberta
Edmonton, Canada
<EMAIL_ADDRESS>
###### Abstract

In this report, we have analyzed available cloud detection technique using
sentinel hub. We have also implemented spatial attention generative
adversarial network and improved quality of generated image compared to
previous solution [7].

##  I Introduction

Observing and monitoring the Earth’s surface from space is essential to handle
climate change, natural resources management, and disaster mitigation. Remote
Sensing (RS) imagery is critical to perform these challenging operations.
However, clouds are one of the common factors which can create a complete or
partial obstacle in the scene.

Satellites examine the Earth’s surface across visible and infrared bands to
capture images. The resulting multispectral image has more channels compared
to an RGB image, corresponding to a different wavelength. Approximately half
of the Earth’s surface is covered in opaque clouds with an additional 20%
being obstructed by cirrus or thin clouds. In some places, the clouds can
cover the Area of Interest (AOI) for several months in a row. An abundance of
research has already been done to remove the clouds and get the accurate
pixels in the satellite images. The approaches include Traditional and ML-
based solutions.

##  II Motivation

In today’s world, Satellite images depict the necessary information about the
Earth’s surface. This information can be crucial in various studies like
disaster management and prevention. Cloud can impede the visibility of these
image scenes. Thus, it is necessary to identify and replace the clouds with
accurate pixels to make the images efficient for further general or ML-based
studies. Researchers have done several types of research to detect and replace
the clouds from satellite images. These studies can be categorized into the
Traditional approaches and advanced ML-based approaches.

Traditional approaches can be classified into three main groups: 1) Multi-
Spectral, 2) Multi-Temporal, and 3) Inpainting. Multi-Spectral approaches can
be applied in partially cloud-affected images. Multi-Temporal approaches use
previously taken data to replace the cloud in the image. Lastly, Inpainting
approaches can use cloud-detection algorithms to identify the cloudy part of
an image so as to reconstruct it using the other non-cloudy ones.

In parallel with the traditional approaches, advanced machine learning-based
algorithms like deep Convolutional Neural Networks (CNNs), Generative
Adversarial Networks (GANs) have proven their efficiency to restore the
missing information. However, the performance of the model can further be
improved by ensembling more neural network algorithms.

##  III Related Work

For cloud removal from satellite images, many different approaches have been
proposed over the years. They can be broadly classified into the following
categories: Traditional approaches which involve detection of clouds from
images and using image inpainting techniques to fill the cloudy areas with the
help of cloud-free images taken over different time intervals. In [1], a
similar approach for removal of clouds from satellite images is presented. It
involves usage of a threshold value to detect cloudy pixels from satellite
images and replacing those pixels with the data from another images of the
same location. Due to variations of patterns observed in clouds, selection of
one definite threshold is a challenge.

Generative Adversarial Network(GAN) based approaches: In [2], authors have
proposed a cyclic generative adversarial network, that eliminates the need of
a paired training dataset and Synthetic Aperture Radar (SAR) images. They
introduced a set of generator and discriminator which compete with each other
until the generator can produce realistic cloud free images. Due to lack of
paired cloudy-cloud free images they didn’t have any quantitative results.

Combination of multiple models based on deep learning techniques: In [3],
authors have introduced two new paired datasets which they have created using
real-world images available from European Space Agency satellite Sentinel-2
and a novel Spatiotemporal Generator Networks (STGAN) which utilizes RGB and
IR data to determine relationships across multiple images over a location.

In [4], SpaceEye, a first of its kind model which utilizes both multi-modal
and multi-temporal satellite images. SpaceEye takes multi spectral and SAR
images as input, and uses a combination of a coarse and refinement network to
patch the missing cloud pixels in an image. The predictions generated by this
model are validated using a discriminator network which further improves
quality of SpaceEye. Authors also stated the applications of SpaceEye by
conducting case-studies.

In [5], authors have been able to generate the entire image instead of only
cloudy pixels by using a hierarchical combination of 3 neural networks -
conditional Generative Adversarial Network (cGAN), Convolutional Long Short-
Term Memory (ConvLSTM) and U-shaped Convolutional Neural Network. This method
has shown encouraging results when compared with state-of-the-art techniques.

##  IV Dataset

We have the RICE dataset, which was proposed in [6]. The RICE1 dataset
contains 500 data samples with each sample having a cloudy image and a
cloudless image under 512×512 resolution. The dataset is collected by Google
Earth, and the cloudy/cloudless images are obtained by setting the cloud layer
whether to display. Figure 1 shows some examples of RICE1 dataset.

Figure 1: Upper row shows cloudy images and lower row shows cloud-free images
from RICE1 dataset

##  V Implementation

###  V-A Cloud detection using Sentinelhub and s2cloudless

We have implemented a process for detecting clouds from sentinelheub2
satellite images using Sentinelhub and s2cloudless. The method receives
dynamic satellite data from the Sentinel2 server using API requests. Then it
generates the cloud mask of the particular satellite image using
S2PixelCloudDetector. We used 10 out of 13 bands to create the cloud mask.

Figure 2: Sentinel2 image Figure 3: Sentinel2 image after cloud detection

We have also checked the detection results for only the RGB bands of the
satellite image of sentinelhub. This time, we have taken a series of satellite
images of different timestamps to detect the clouds. It can also identify the
clouds efficiently using only the RGB bands.

Figure 4: Sentinel2 image after cloud detection

###  V-B Cloud Removal using GAN

We have used [7] as reference to implement cloud removal from the satellite
images. In this paper, author’s have proposed a spatial attention generative
adversarial networks. We have taken reference of spatial attentive block and
spatial attentive residual block to create the architecture of generator. The
Figure 5 and Figure 6 show the block.

Figure 5: Spatial Attentive Block Figure 6: Spatial Attentive Residual Block

We have implemented 2 new approaches to improve results of cloud removal using
spatial attention generative adversarial network. In the first approach we
have changed the architecture of the generator network. This new architecture
first extracts the features from image and then generate spatial attention
map. This generated spatial map is given as input to extract features and
generate new spatial map from that features and in the last given to
convolution layer, which give reconstructed cloud-free image as output. Using
this architecture, in the result, the PSNR is slightly increased and SSIM is
similar to previous solution. For the second approach, we have used
8-neighbour pixels to generate spatial attention map. This approach shown a
significant increase in PSNR and slight improvement in SSIM. Using 8-neighbour
it has made significant change in the genrated attention map which resulted in
better reconstruction of cloud-free images from cloudy images.

To train the GAN model, we have used 400 images from RICE dataset from which
320 images are used to train the model and 80 images are used to validate the
model. We trained the model for 100 epochs. Table 1 shows the comparison
between 3 solutions to remove clouds from satellite images in terms of PSNR
and SSIM. Figure 7, Figure 8 and Figure 9 shows comparison between result of
previous solution and our proposed solution.

| PSNR | SSIM
---|---|---
Previous Solution | 26.4239 | 0.9545
Our Solution | 29.4300 | 0.9670
TABLE I: Comparison of PSNR & SSIM Figure 7: Ground Truth Figure 8: Our
solution Figure 9: Previous Solution

##  VI Future Work and Conclusion

In this project, we have implemented previously available cloud detection
technique using sentinel hub tool and cloud removal techniques using GAN and
tried to improve the available algorithms. Using GAN, we have successfully
generated the cloud-free images from cloudy satellite images and improved the
reconstructed image quality. For the future work, we can improve the model
architecture using U-net like architecture and integrating spatial and
temporal information available of satellite images by training the model using
multiple cloudy images.

## References

  * [1] Das, Sanju, Purnendu Das, and Bishwa Ranjan Roy. ”Cloud Detection and Cloud Removal of Satellite Image—A Case Study.” Trends in Communication, Cloud, and Big Data (2020): 53-63.
  * [2] Singh, Praveer, and Nikos Komodakis. ”Cloud-gan: Cloud removal for sentinel-2 imagery using a cyclic consistent generative adversarial networks.” IGARSS 2018-2018 IEEE International Geoscience and Remote Sensing Symposium. IEEE, 2018.
  * [3] Uzkent, Burak Uzkent, et al. ”Cloud Removal in Satellite Images Using Spatiotemporal Generative Networks.” (2019).
  * [4] Zhao, Mingmin, Peder A. Olsen, and Ranveer Chandra. ”Seeing Through Clouds in Satellite Images.” arXiv preprint arXiv:2106.08408 (2021).
  * [5] Sebastianelli, Alessandro, et al. ”Sentinel-1 and Sentinel-2 Spatio-Temporal Data Fusion for Clouds Removal.” arXiv e-prints (2021): arXiv-2106.
  * [6] Lin, Daoyu, et al. ”A remote sensing image dataset for cloud removal.” arXiv preprint arXiv:1901.00600 (2019).
  * [7] Pan, Heng. ”Cloud Removal for Remote Sensing Imagery via Spatial Attention Generative Adversarial Network.” arXiv preprint arXiv:2009.13015 (2020).