1.4pt

# Predicting and explaining nonlinear material response using deep Physically
Guided Neural Networks with Internal Variables

Javier Orera-Echeverria
University of Zaragoza
Mariano Esquillor, s.n. 50018, Zaragoza
<EMAIL_ADDRESS>
&Jacobo Ayensa-Jiménez
Aragon Institute of Engineering Research
University of Zaragoza
Mariano Esquillor, s.n. 50018, Zaragoza
<EMAIL_ADDRESS>
&Manuel Doblaré∗
Aragon Institute of Engineering Research
University of Zaragoza
Mariano Esquillor, s.n. 50018, Zaragoza
<EMAIL_ADDRESS>

###### Abstract

Nonlinear materials are often difficult to model with classical state model
theory because they have a complex and sometimes inaccurate physical and
mathematical description or we simply do not know how to describe such
materials in terms of relations between external and internal variables. In
many disciplines, Neural Network methods have arisen as powerful tools to
identify very complex and non-linear correlations. In this work, we use the
very recently developed concept of Physically Guided Neural Networks with
Internal Variables (PGNNIV) to discover constitutive laws using a model-free
approach and training solely with measured force-displacement data. PGNNIVs
make a particular use of the physics of the problem to enforce constraints on
specific hidden layers and are able to make predictions without internal
variable data. We demonstrate that PGNNIVs are capable of predicting both
internal and external variables under unseen load scenarios, regardless of the
nature of the material considered (linear, with hardening or softening
behavior and hyperelastic), unravelling the constitutive law of the material
hence explaining its nature altogether, placing the method in what is known as
eXplainable Artificial Intelligence (XAI).

_K_ eywords Nonlinear computational solid mechanics $\cdot$ Deep Neural
Network $\cdot$ Internal Variables $\cdot$ Explainable Artificial Intelligence
$\cdot$ Physics-Informed Machine Learning $\cdot$ Physically Guided Neural
Networks

##  1 Introduction

It is of common knowledge that our everyday life is being dramatically
challenged by Big Data and Artificial Intelligence (AI). According to the
International Data Corporation, the Global Datasphere (the summation of all
data, whether created, captured, or replicated) will grow from 33 ZB in 2018
to 175 ZB by 2025 [1]. This is mainly due to the explosion of social networks,
e-commerce and marketing, and the extension of the Internet of Things (IoT).
Among the top eight companies in terms of market capitalization, five are
based on data value leverage [2]. This huge amount of available information
justifies the prosperity of data science and data-based decision in fields as
diverse as sociology, economics, engineering and medicine.

As a response, Machine Learning (ML) methods have become today one of the main
tools in business, but also in science and technology. These methodologies
enable the extraction of information from data that would be intractable by
means of traditional methods [3]. They try to mimic the process of human
knowledge acquisition and structuring taking advantage of the aforementioned
advances in data generation, management, and storage, as well as huge
improvements in the performance of computers and algorithms [4]. In the
special case of Scientific ML [5], the natural adaptation of many supervised
as well as unsupervised ML algorithms to the vectorized representation that
most physical problems exhibit, makes the study of the convergence between
both of special importance.

Data-driven methods are used in many different physical disciplines such as
chemical and electrical processes [6], biology [7], spoken language
recognition [8] and a long etcetera. However, the link between classical
physical modelling and data-driven methods has not been quite clear so far,
since the physical description of most systems was built on the basis of
empirical knowledge rather than large data-bases. Nonetheless, the new era of
computation and Big Data has opened new perspectives where data can be
incorporated in this physical description in a consistent and comprehensive
way. Promising improvements can be therefore spotted, such as new forms of
empiricism that declare ”the end of theory” and the impending advancement of
data-driven methods over knowledge-driven science [9].

One of the most prominent strategies that has proven to be especially prolific
in recent years is the use of Artificial Neural Networks (ANNs). Since 1958,
when Rosenblatt developed the _perceptron_ [10], many works have been devoted
to ANNs, some of them demostrating their character as universal approximators
[11, 12, 13, 14, 15]. However, it has been only in the last decade of the XXth
century, thanks to the important progress of high performance computing
capabilities and the combination of back-propagation [16] with stochastic
gradient descent [17] algorithms, that ANNs have become a booming technology.
The progress has accelerated in the last decade with the advent of
Convolutional Neural Networks (CNNs) [18] and Recurrent Neural Networks (RNNs)
[19], in what is known nowadays as Deep Learning (DL) [20], and culminating
with the attention mechanism [21], transformer models and generative AI, whose
impact has gained nowadays great popularity thanks to Large Language Models
(LLMs) [22].

Leaving aside the progress of the DL as a research field, AI approaches have
changed the way we conceive science. On the one hand, AI has been used to
discover the hidden physical structure of the data and unravel the equations
of a system [23, 24, 25]. On the other hand, the tremendous predictive power
of AI has been blended with the scientific consistency of the explicit
mathematical representation of physical systems through the concept of _data-
driven_ models for simulation-based engineering and sciences (SBES) [26]. The
latter in turn may be done by the combination of raw data and physical
equations [27, 28], by enforcing a metriplectic structure to the model,
related with the fulfillment of thermodynamic laws [29] or by defining the
specific structure of the model [30, 31]. These novel _data-driven science_
approaches, coined as Scientific Machine Learning or Physics-Informed Machine
Learning (PIML) arise therefore with the main purpose of turning apparently
not-physically meaningful data-driven models, where approaches such as ANNs
have excelled, into physics-aware models.

However, the interplay between data and physical sciences has not been exempt
from setbacks of different nature. In fact, the use of complex DL models does
not fit well with the study of physical problems. Furthermore, in many
physical problems, many variables are involved, interacting in complex and
non-stationary ways. This requires huge amounts of data to get accurate
predictions using ANN techniques, sometimes in regions of the solution space
that are difficult to access or uncommon and therefore difficult to be
sampled. As a consequence, due to the bias–variance trade-off, poorly
extrapolation capacity is obtained out of the the usual data range for models
looking for recreating complex physics [32].

In addition, a physically based model is not only useful for making
predictions, but also to help in acquiring new knowledge by the interpretation
of its structure, parameters, and mathematical properties. Physical
interpretability is, in most cases, at least as important as predictive
performance. However, it is known that interpretability is one of the main
weaknesses of ANNs, as the acquired knowledge is encoded in the strength of
multiple connections, rather than stored at specific locations. That explains
the huge efforts that are currently being made towards “whitening” the “black-
box” character of ANN [33, 34] in what has been coined as eXplainable
Artificial Intelligence (XAI) [35, 36]. In the context of data-driven
simulation-based engineering and sciences (DDSBES) [37]. Two ways of
proceeding can be distinguished: building specific ANNs structures endowed
with the problem equations [38, 39], also known as the _inductive bias_
approach, and/or by regularizing the loss function using this same physical
information [40, 41, 42].

The approach to be followed depends decisively upon the data availability and
the way in which this data is used. If we follow a supervised approach, there
are two possibilities.

  * •

The first one assumes that we know the whole physics of the problem. In that
situation, supervised ML is used for the sake of computational requirements.
In that sense, ANNs act as Reduced Order Models (ROMs) that can be used as a
surrogate in problems involving optimization or control. Hybridizing DL with
physical information is a way of improving standard DL methods in terms of
data requirements, less expensive training or noise filtering at the
evaluation step, thanks to regularization [38, 40]. Therefore, we are
interested in the predictive capacity of the approach.

  * •

The second one assumes that we know some of the physics of the problem. In
that situation, we are interested in knowing the hidden physics that remains
unknown, expressed in terms of some model parameters [42] or functional
relations [43]. For that reason, we are rather interested in the explanatory
character of the method.

In other situations, we follow unsupervised approaches. This may be indeed due
to two different possibilities:

  * •

If we know the whole physics of the problem, the use of ANNs is merely
instrumental and is used as an alternative way of solving numerically some
system of Partial Differential Equations (PDEs) that require an important
computational effort [44].

  * •

When some of the physics of the problem are unknown, the intrinsic variability
of the data (for instance when measuring spatial or temporal fields) may be
exploited for an unsupervised discovery of some hidden constitutive models
[45, 46]. For that reason, we are rather interested in the explanatory
character of the method.

However, in the context of the IoT, where data quantity generally dominates
over data quality, the data availability and variability is key, and it is
difficult to guarantee whether, considering for instance the specific case of
computational solid mechanics, “the combination of geometry and loading
generates sufficiently diverse and heterogeneous strain states to train a
generalizable constitutive model with just a single experiment” [47].
Therefore, there are no other alternatives rather than introducing all the
control or measurable variables in the workflow, while maintaining the
desirable properties of the PIML approach, namely, its ability to get fast
predictions in real time (for optimization and control issues) together with
its explanatory capacity.

In this work, we demonstrate, how, in the context of computational solid
mechanics, Physically Guided Neural Networks with Internal Variables (PGNNIVs)
enable the compliance with these requirements particularly well. PGNNIVs
comprise ANNs that are able to incorporate some of the known physics of the
problem, expressed in terms of some measurable variables (for instance forces
and displacements) and some hidden ones (for instance stresses). Their
predictive character, improving many of the features of conventional ANN [38],
allow for fast and accurate predictions. In addition, PGNNIVs are able to
unravel the constitutive equations of different materials from unstructured
data, that is, uncontrolled test data obtained from system monitorization.

The content of this paper is structured as follows. First, the methodology is
described, including a brief overview of the state of the art, the use of
PGNNIVs in computational mechanics, the computational treatment of the
physical tensorial fields and operators, as well as the data-set generation
and training process. Then, the main results are presented, of both the
predictive and explanatory capacity of the method. Finally, conclusions are
summarized, together with the main limitations and a brief overview of future
work.

##  2 Methods

###  2.1 Brief overview of Physics Informed Machine Learning in computational
solid mechanics

PDEs are the standard way to describe physical systems under the continuum
setting thanks to their overarching capacity to model extremely different and
complex phenomena. However, analytic solutions are most times difficult or
even impossible to find. That is the reason why numerical methods have become
the universal tool to obtain approximate but accurate solutions to PDEs. These
methods consider a given discretization in space and time, which results in an
algebraic (in general non-linear) system that is then solved by means of
standard matrix manipulation.

In the last three decades, however, attempts to solve PDEs from a data-driven
point of view have been numerous. First tentatives [48, 49] generalized
earlier ideas for Ordinary Differential Equations (ODEs). Since then, many
different approaches have been proposed, from collocation methods [50, 51,
52], variational/energy approaches [53, 54, 55], loss regularization using
physical or domain knowledge [56, 57, 41, 58], to the most recent approaches
using automatic differentiation [42], that is nowadays known as Physics-
Informed Neural Networks (PINNs). Other works have extensively tried to
address this challenge using the stochastic representations of high-
dimensional parabolic PDEs [59, 60, 61].

In order to provide the data-driven models with a meaningful physical
character, remarkable efforts have been done in the recent years to embed
physical information into data-driven descriptions. The potential of solving
inverse problems with linear and non-linear behavior in solid mechanics, for
example, has been explored using DL methods [62], where the forward problem is
solved first to create a database, which is then used to train the ML
algorithms and determine the boundary conditions from assumed measurements.
Other approaches initially build a constitutive model into the framework by
enforcing constitutive constraints, and aim at calibrating the constitutive
parameters [63].

In this context, clearly differentiating between external and internal
variables becomes an important factor when approaching complex physical
problems, but this is always disregarded from a data-driven viewpoint.
External variables are those observable, measurable variables of the system,
that can be obtained directly from physical sensors such as position,
temperature or forces; internal variables are non-observable (not directly
measurable) variables, that integrate locally other observable magnitudes and
depend on the particular internal structure of the system [64]. This is very
important to consider in the ML framework since predictions associated to an
internal state model require explicitly the definition of the cloud of
experimental values that identifies the internal state model [65]. This
implies “measuring”, or better, assuming values for non-observable variables
[66, 67]. This is for example, the case of stresses in continuum mechanics
that can be determined _a priori_ only after making strong assumptions such as
their uniform distribution in the center section of a sample under uniform
tension.

An alternative is the use of PGNNIVs. This new methodological appraisal
permits to predict the values of the internal variables by mathematically
constraining some hidden layers of a deep ANN (whose structural topology is
appropriately predefined) by means of the fundamental laws of the continuum
mechanics such as conservation of energy or mechanical momenta that relate
internal non-measurable with external observable variables. With this, it is
possible to transform a pure ML based model into a physically-based one
without giving up the powerful tools of DL, including the implicit correlation
between observable data and the derived predictive capacity. This way, not
only the real internal variables of the problem are predicted, but also the
data needed to train the network decreases, convergence is reached faster,
data noise is better filtered and the extrapolation capacity out of the range
of the training data-set is improved, as recently demonstrated in [38]. In
line with the terminology and general framework used in PIML [68], PGNNIVs
showcase an intuitive interplay between an _inductive bias_ approach and a
_learning bias_ appraisal, where physical constrains are incorporated by means
of a physics-informed likelihood, i.e. additional terms in the loss function
(also known as _collocation_ or _regularization losses_).

###  2.2 Physically Guided Neural Networks with Internal Variables in
computational mechanics

####  2.2.1 Revisiting PGNNIVs

PGNNIVs are in essence a generalization of PINNS. In the latter, physical
equations constrain the values of output variables to belong to a certain
physical manifold that is built from the information provided by the data and
the specific form of the PDE considered. One of the shortcomings of PINNs is
that, in general, only simple ordinary PDEs that contain a few free parameters
and are closed-form can be considered. In contrast, PGNNIVs do not only apply
to scenarios where PDEs involve many parameters and have complex forms, but
also to those where a mathematical description is not available. In fact, this
new paradigm embodies a unique architecture where the values of the neuron
variables in some intermediate layers acquire physical meaning in an
unsupervised way, providing the network with an inherent explanatory capacity.

The main differentiating features that, up to the authors’ knowledge,
constitute a relevant contribution to the advances of data-driven physical
modelling and its particular application to nonlinear mechanics are two-
folded: first, physical constraints are applied in predefined internal layers
(PILs), in contrast to previous works on PINNs. Secondly, and even most
important, PGNNIVs are able to predict and explain the nature of the system
all at once, i.e. predictability of the variables as well as explainability of
the constitutive law are ensured and learned altogether. In all related works,
modeling assumptions on the constitutive law of the correspondent materials
are directly imposed, so that the material response complies with certain
constrains. On the contrary, PGNNIVs only enforce universal laws of the system
(e.g. balance equations) and no prior knowledge on the constitutive model is
incorporated.

Classical Deep Neural Networks (DNNs) are often represented as _black boxes_
that theoretically can compute and learn any kind of function correlating the
input and output data [69]. In particular, they perform very well in areas of
science and technology where complex functions convey good approximations of
the governing phenomena. Although there exist some heuristic rules [70], these
_black boxes_ are usually trained via trial and error. Adding a physical
meaning to the hidden layers and constraining them by adding an extra term to
the cost function, has already proven to have significant advantages such as
less data requirements, higher accuracy and faster convergence in real
physical problems, as well as model unravelling capacities [38]. The basic
principles of PGNNIVs are briefly exposed in the next lines.

Let us consider a set of continuous Partial Differential Equations (PDEs) of
the form

$\displaystyle\mathcal{F}(u,v)$ $\displaystyle=f,\text{ in }\Omega,$ (1a)
$\displaystyle\mathcal{G}(u,v)$ $\displaystyle=g,\text{ in }\partial\Omega,$
(1b) $\displaystyle\mathcal{H}(u)$ $\displaystyle=v,\text{ in }\Omega,$ (1c)

where $u$ and $v$ are the unknown fields of the problem, $\mathcal{F}$ and
$\mathcal{H}$ are functionals representing the known and unknown physical
equations of the specific problem, $\mathcal{G}$ is a functional that
specifies the boundary conditions, and $f$ and $g$ are known fields.

The continuous problem has its analogous discretized representation in finite-
dimensional spaces in terms of vectorial functions $\boldsymbol{F}$,
$\boldsymbol{G}$ and $\boldsymbol{H}$ and nodal values $\boldsymbol{u}$,
$\boldsymbol{v}$, $\boldsymbol{f}$ and $\boldsymbol{g}$. Particularly,
$\boldsymbol{u}$ are the solution field nodal values and $\boldsymbol{v}$ are
the unknown internal field variables at the different nodes. The
discretization may be done using any discretization technique, such as the
Finite Element Method (FEM). Hence, Eqs. (1) become

$\displaystyle\boldsymbol{F}(\boldsymbol{u},\boldsymbol{v})=\boldsymbol{f},\text{
in }\Omega,$ (2a)
$\displaystyle\boldsymbol{G}(\boldsymbol{u},\boldsymbol{v})=\boldsymbol{g},\text{
in }\partial\Omega,$ (2b)
$\displaystyle\boldsymbol{H}(\boldsymbol{u})=\boldsymbol{v},\text{ in
}\Omega.$ (2c)

A PGNNIV may be defined for a problem of type (2) in the following terms:

$\displaystyle\boldsymbol{y}$ $\displaystyle=\mathsf{Y}(\boldsymbol{x}),$
$\displaystyle\boldsymbol{v}$ $\displaystyle=\mathsf{H}(\boldsymbol{u}),$
$\displaystyle\boldsymbol{x}$
$\displaystyle=\boldsymbol{I}(\boldsymbol{u},\boldsymbol{f},\boldsymbol{g}),$
$\displaystyle\boldsymbol{y}$
$\displaystyle=\boldsymbol{O}(\boldsymbol{u},\boldsymbol{f},\boldsymbol{g}),$
$\displaystyle\boldsymbol{R}(\boldsymbol{u},\boldsymbol{v},\boldsymbol{f},\boldsymbol{g})$
$\displaystyle=0,$

where:

  1. 1.

$\boldsymbol{u}$, $\boldsymbol{f}$ and $\boldsymbol{g}$ are the measurable
variables of the problem.

  2. 2.

$\boldsymbol{x}$ and $\boldsymbol{y}$ are the input and output variables
respectively and will be defined depending on which relation
$\boldsymbol{x}\mapsto\boldsymbol{y}$ wants to be predicted.

  3. 3.

$\boldsymbol{I}$ and $\boldsymbol{O}$ are functions that compute the input
$\boldsymbol{x}$ and the output $\boldsymbol{y}$ of the problem from the
measurable variables. In other words, functions $\boldsymbol{I}$ and
$\boldsymbol{O}$ define the data used as starting point to make predictions,
$\boldsymbol{x}$, and the data that we want to predict, that is,
$\boldsymbol{y}$.

  4. 4.

$\boldsymbol{R}$ are the physical constraints, related to the relations given
by $\boldsymbol{F}$ and $\boldsymbol{G}$.

  5. 5.

$\mathsf{Y}$ and $\mathsf{H}$ are DNN models:

     * •

$\mathsf{Y}$ is the predictive model, whose aim is to infer accurate values
for the output variables for a certain input set, that is, to surrogate the
relation $\boldsymbol{x}\mapsto\boldsymbol{y}$.

     * •

$\mathsf{H}$ is the explanatory model, whose objective is to unravel the
hidden physics of the relation $\boldsymbol{u}\mapsto\boldsymbol{v}$.

####  2.2.2 Adaptation to computational solid mechanics.

Our aim now is to reframe Eqs. (1) in the context of solid mechanics. To fix
ideas, although it is not difficult to adapt the methodology to other
constitutive models, we restrict this analysis to hyperelastic solids with
constant and known density $\rho$. First, we have to consider equilibrium
equations (momentum conservation) in the domain $\Omega$. In spatial
coordinates, equilibrium reads

$\mathrm{div}(\boldsymbol{\sigma})+\rho\boldsymbol{b}=\boldsymbol{0},$ (4)

where $\boldsymbol{\sigma}$ is Cauchy stress tensor, $\rho$ is the density,
$\boldsymbol{b}$ the spatial volumetric body force field and $\mathrm{div}$ is
the divergence operator in spatial coordinates. In material coordinates,
equilibrium reads

$\mathrm{DIV}(\boldsymbol{P})+\rho\boldsymbol{B}=\boldsymbol{0},$ (5)

where now
$\boldsymbol{P}=\mathrm{det}(\boldsymbol{F})\boldsymbol{\sigma}\boldsymbol{F}^{-\intercal}$
is the first Piola-Kirchhoff stress tensor,
$\boldsymbol{B}=\mathrm{det}(\boldsymbol{F})\boldsymbol{b}$ is the reference
volumetric body force field, and $\mathrm{DIV}$ is the divergence operator in
material coordinates.

If $\boldsymbol{\xi}=\chi(\boldsymbol{\Xi})$ is the motion function that
relates spatial ($\boldsymbol{\xi}$) and material ($\boldsymbol{\Xi}$)
coordinates, we define the deformation gradient tensor $\boldsymbol{F}$ as

$\boldsymbol{F}=\mathrm{GRAD}(\chi)=\mathrm{GRAD}\otimes\boldsymbol{\xi},$ (6)

where $\mathrm{GRAD}$ is the gradient operator in material coordinates. Eq.
(6) shows that $\boldsymbol{F}$ is a potential tensor field so
$\boldsymbol{F}$ must satisfy the following compatibility equation in each
connected component of $\Omega$:

$\mathrm{ROT}(\boldsymbol{F})=\boldsymbol{0},$ (7)

where $\mathrm{ROT}$ is the rotational in material coordinates.

For hyperelastic materials, the Cauchy stress tensor $\boldsymbol{\sigma}$ is
related to the deformation gradient tensor by means of the equation

$\boldsymbol{\sigma}=\frac{1}{\mathrm{det}(\boldsymbol{F})}\frac{\partial\Psi}{\partial\boldsymbol{F}}\cdot\boldsymbol{F}^{\intercal},$
(8)

where $\Psi$ is the strain energy function expressed as a function of the
deformation state given by $\boldsymbol{F}$,
$\Psi=\mathfrak{F}(\boldsymbol{F})$. Obtaining this particular function is the
subject of research of material sciences and many approaches are possible,
ranging from phenomenological descriptions to mechanistic and statistical
models.

Finally, Eqs. (4) or (5), (7) or (6), and (8) must be supplemented with
appropriate boundary conditions. We distinguish here between _essential_
boundary conditions and _natural_ boundary conditions. The former define the
motion of the solid at some boundary points $\Gamma_{E}\subset\partial\Omega$:

$\boldsymbol{\xi}=\bar{\boldsymbol{\xi}},$ (9)

whereas the latter define the traction vector at some other boundary points
$\Gamma_{N}\subset\partial\Omega$:

$\boldsymbol{\sigma}\cdot\boldsymbol{n}=\bar{\boldsymbol{t}},$ (10)

where $\boldsymbol{n}$ is the outwards normal vector (in the spatial
configuration) and $\bar{\boldsymbol{\xi}}$ and $\bar{\boldsymbol{t}}$ are
known values of the solid motion and spatial traction forces respectively. The
material analogous of Eq. (10) is

$\boldsymbol{P}\cdot\boldsymbol{N}=\bar{\boldsymbol{T}},$ (11)

where now $\boldsymbol{N}$ and $\bar{\boldsymbol{T}}$ are the material
analogous of $\boldsymbol{n}$ and $\bar{\boldsymbol{t}}$.

Let us assume that it is possible to measure (for instance using Digital Image
Correlation techniques [71]) the system response in terms of its motion given
by the map $\boldsymbol{\xi}=\chi(\boldsymbol{\Xi})$, that is,
$\boldsymbol{\xi}$ is a measurable variable. Let us also assume that we can
measure the volumetric loads,
$\boldsymbol{b}=\boldsymbol{b}(\boldsymbol{\xi})$, as well as the prescribed
traction forces $\bar{\boldsymbol{t}}$ (or, in an equivalent manner,
$\boldsymbol{B}=\boldsymbol{B}(\boldsymbol{\Xi})$ and $\bar{\boldsymbol{T}}$).
Therefore, using Eqs. (6) and (8) it is possible to express:

$\displaystyle\boldsymbol{F}$ $\displaystyle=\mathcal{A}(\boldsymbol{\xi}),$
(12a) $\displaystyle\boldsymbol{\sigma}$
$\displaystyle=\mathcal{B}(\Psi,\boldsymbol{\xi}),$ (12b)

for some appropriate differential operators $\mathcal{A}$ and $\mathcal{B}$.
Therefore, it is possible to recast hyperelastic solid mechanics as

$\displaystyle\mathrm{div}(\boldsymbol{\sigma}(\Psi,\boldsymbol{\xi}))$
$\displaystyle=-\rho\boldsymbol{b}\quad\mathrm{in}\quad\Omega,$ (13a)
$\displaystyle\boldsymbol{\xi}$
$\displaystyle=\bar{\boldsymbol{\xi}}\quad\mathrm{in}\quad\Gamma_{E},$ (13b)
$\displaystyle\boldsymbol{\sigma}\cdot\boldsymbol{n}$
$\displaystyle=\bar{\boldsymbol{t}}\quad\mathrm{in}\quad\Gamma_{N},$ (13c)
$\displaystyle\Psi$
$\displaystyle=\mathfrak{F}(\boldsymbol{F}(\boldsymbol{\xi}))\quad\mathrm{in}\quad\Omega.$
(13d)

where Eqs. (13a) and (13c) may be eventually substituted by their material
analogous. Groupping Eqs. (13b) and (13c), it is clear that we can express
computational solid mechanics for hyperelastic materials as

$\displaystyle\mathcal{F}(\boldsymbol{\xi},\boldsymbol{\sigma})=\boldsymbol{b},\text{
in }\Omega,$ (14a)
$\displaystyle\mathcal{G}(\boldsymbol{\xi},\boldsymbol{\sigma})=(\bar{\boldsymbol{\xi}},\bar{\boldsymbol{t}})\text{
in }\partial\Omega,$ (14b)
$\displaystyle\mathcal{H}(\boldsymbol{\xi})=\boldsymbol{\sigma},\text{ in
}\Omega,$ (14c)

that is in the form of Eqs. (1) with $u=\boldsymbol{\xi}$,
$v=\boldsymbol{\sigma}$, $f=\boldsymbol{b}$ and
$g=(\bar{\boldsymbol{\xi}},\bar{\boldsymbol{t}})$.

##### Particularization to small strains solid mechanics.

In small strains solid mechanics, it is common to work with the displacement
field $\boldsymbol{U}=\boldsymbol{\xi}-\boldsymbol{\Xi}$ and to define the
displacement gradient tensor $\boldsymbol{J}=\mathrm{GRAD}(\boldsymbol{U})$.
In that case, the constitutive equation is formulated as

$\boldsymbol{\sigma}=\mathfrak{G}(\boldsymbol{\varepsilon}),$ (15)

where $\boldsymbol{\varepsilon}$ is the Cauchy small strain tensor

$\boldsymbol{\varepsilon}=\mathrm{symgrad}(\boldsymbol{U})=\frac{1}{2}(\boldsymbol{J}+\boldsymbol{J}^{\intercal}),$
(16)

and $\mathfrak{G}$ is a tensor map. With these considerations, the equations
of the problem are now written as

$\displaystyle\mathcal{F}(\boldsymbol{U},\boldsymbol{\sigma})$
$\displaystyle=\boldsymbol{b},\text{ in }\Omega,$ (17a)
$\displaystyle\mathcal{G}(\boldsymbol{U},\boldsymbol{\sigma})$
$\displaystyle=(\bar{\boldsymbol{U}};\bar{\boldsymbol{t}}),\text{ in
}\partial\Omega,$ (17b) $\displaystyle\mathcal{H}(\boldsymbol{U})$
$\displaystyle=\boldsymbol{\sigma},\text{ in }\Omega,$ (17c)

again in the form of Eqs. (1) with $u=\boldsymbol{U}$,
$v=\boldsymbol{\sigma}$, $f=\boldsymbol{b}$ and
$g=(\bar{\boldsymbol{U}},\bar{\boldsymbol{t}})$.

Once we have discretized the problem, our aim is to predict a motion field
$\boldsymbol{\xi}$ (or a displacement field $\boldsymbol{U}$) from a
particular load case, expressed in terms of the volumetric loads and the
natural boundary conditions111It is important to recall that, as essential
boundary conditions can be measured as an output variable, it is not necessary
to include them as inputs of our problem., $\bar{\boldsymbol{t}}$, therefore
$\boldsymbol{I}(\boldsymbol{u},\boldsymbol{f},\boldsymbol{g})=(\boldsymbol{f},\boldsymbol{g})=(\boldsymbol{b},\bar{\boldsymbol{t}})$.
With these last remarks, the PGNNIV problem is stated for finite solid
mechanics as

$\displaystyle\boldsymbol{\xi}$
$\displaystyle=\mathsf{Y}(\bar{\boldsymbol{t}})$ (18a)
$\displaystyle\boldsymbol{\sigma}$
$\displaystyle=\mathsf{H}(\mathrm{KIN}(\boldsymbol{\xi}))$ (18b)
$\displaystyle\boldsymbol{x}$
$\displaystyle=(\bar{\boldsymbol{t}},\boldsymbol{b})$ (18c)
$\displaystyle\boldsymbol{y}$ $\displaystyle=\boldsymbol{\xi}$ (18d)
$\displaystyle\boldsymbol{R}(\boldsymbol{\xi},\boldsymbol{\sigma},\bar{\boldsymbol{t}})$
$\displaystyle=(\mathrm{div}(\boldsymbol{\sigma})-\rho\boldsymbol{b};\boldsymbol{\sigma}\cdot\boldsymbol{n}-\bar{\boldsymbol{t}};\boldsymbol{\xi}-\bar{\boldsymbol{\xi}}).$
(18e)

where $\mathrm{KIN}(\boldsymbol{\xi})$ is a selected kinematic descriptor of
the strain state, such as the deformation gradient
$\boldsymbol{F}=\mathrm{GRAD}(\boldsymbol{\xi})$, the right Cauchy - Green
deformation tensor $\boldsymbol{C}=\boldsymbol{F}^{\intercal}\boldsymbol{F}$,
or the Green - Lagrange strain tensor
$\boldsymbol{E}=\frac{1}{2}(\boldsymbol{C}-\boldsymbol{I})$, among others. For
small strains solid mechanics, the methodology simplifies to

$\displaystyle\boldsymbol{U}$ $\displaystyle=\mathsf{Y}(\bar{\boldsymbol{t}})$
(19a) $\displaystyle\boldsymbol{\sigma}$
$\displaystyle=\mathsf{H}(\mathrm{symgrad}(\boldsymbol{U}))$ (19b)
$\displaystyle\boldsymbol{x}$
$\displaystyle=(\bar{\boldsymbol{t}},\boldsymbol{b})$ (19c)
$\displaystyle\boldsymbol{y}$ $\displaystyle=\boldsymbol{U}$ (19d)
$\displaystyle\boldsymbol{R}(\boldsymbol{U},\boldsymbol{\sigma},\bar{\boldsymbol{t}})$
$\displaystyle=(\mathrm{div}(\boldsymbol{\sigma})-\rho\boldsymbol{b};\boldsymbol{\sigma}\cdot\boldsymbol{n}-\bar{\boldsymbol{t}};\boldsymbol{U}-\bar{\boldsymbol{U}}).$
(19e)

The appropriate structure and architecture of $\mathsf{H}$ and $\mathsf{Y}$
depend on the complexity of the material in hands, which is discussed later.

####  2.2.3 Case study: geometry and external forces.

For illustration purposes, the case study considered in this work consists in
a non-uniform biaxial test on a rectangular plate of height $L_{1}=16$ cm and
width $L_{2}=20$ cm, under plane stress. No volumetric loads are incorporated,
that is, $\boldsymbol{b}=\boldsymbol{0}$. We impose a certain arbitrary
compression load profile $p=p(s)$ (where $s$ is the coordinate along the right
and top contour). To accelerate computations we consider a load profile that
is symmetric with respect to the vertical and horizontal axis and acts
perpendicularly to the plate contour, as shown in Figure 1. The symmetry of
the problem allows therefore for the analysis of an equivalent problem by
extracting the upper-right portion of the plate and applying the corresponding
symmetry boundary conditions.

Figure 1: Dimensions and representation of the non-uniform biaxial test on a
2D plate. Load $p=p(s)$ acting perpendicularly to the contour is arbitrary,
provided that it is compatible with the symmetry of the problem. We locate the
origin $(x,y)=(0,0)$ at the geometrical center of the plate.

###  2.3 Data and operators representation

In this section we discuss how the different mathematical objects in our use
case problem (scalar, vectorial and tensorial fields and operators) are
represented. First, we describe how the different fields are encoded and
related to the measured data. Then, we explain how the different operators
involved are built. This includes both the known operators (equilibrium and
compatibility) as well as the unknown relationships comprising the predictive
and explanatory networks, $\mathsf{Y}$ and $\mathsf{H}$ respectively. Then, we
discuss how physical constrains are hardwired into the ANN, so that the built
PGNNIV is tailored towards the discovery of constitutive models that comply
with the physics of the solid mechanics problem, constraining therefore the
learning space and bypassing the parametrization of the constitutive law.

####  2.3.1 Data structures

The data that contains the nodal and element-wise variables (that is,
displacements and stresses/strains respectively) is stored in array
structures. Now we introduce the notation that will be used for referring to a
given tensor field, that is represented by an array $\mathtt{I}$ containing
the information of the tensorial field itself. The different dimensions of the
data are represented using the indexation $\mathtt{I}[I|J|K]$ where $I$ is a
multi-index associated with the discretization of the problem, $J$ with the
tensorial character and $K$ with the data instance. Thus, considering that we
have a data-set of size $N$ and a discretization of size $n_{x}\times n_{y}$,
the displacement field is represented by $\mathtt{U}[i,j|k|l]$ where
$i=1,\ldots,n_{x}$, $j=1,\ldots,n_{y}$, and where $k=1,2$ (2D problem) and
$l=1\ldots,N$, so that

$\mathtt{U}[i,j|k|l]=u_{k}(x_{i},y_{j})$

is the $k$ component of the displacement field evaluated at $(x_{i},y_{j})$
(that is, the node $(i,j)$) corresponding to the data $l$. Analogously, the
strain and stress fields are represented respectively by
$\mathtt{E}[i,j|k,l|m]$ and $\mathtt{S}[i,j|k,l|m]$ where
$i=1,\ldots,n_{x}-1$, $j=1,\ldots,n_{y}-1$, $k,l=1,2$ and $m=1\ldots,N$, such
that, for instance,

$\mathtt{E}[i,j|k,l|m]=E_{kl}\left(\frac{1}{2}(x_{i}+x_{i+1}),\frac{1}{2}(y_{j}+y_{j+1})\right)$

is the $k,l$ component of the strain tensor evaluated at the element $(i,j)$
corresponding to the data $m$.

Finally, the traction forces, $\bar{\boldsymbol{t}}^{\mathrm{top}}$ and
$\bar{\boldsymbol{t}}^{\mathrm{right}}$, which are treated as the inputs of
our problem (provided the volume forces are not considered), are represented
by $\mathtt{t}^{\mathrm{top}}[i|j|k]$ where $i=1,\ldots,n_{x}$, $j=1,2$, and
$k=1\ldots,N$, so that

$\mathtt{t}^{\mathrm{top}}[i|j|k]=t^{\mathrm{top}}_{j}(x_{i},L_{1}/2),$

and $\mathtt{t}^{\mathrm{right}}[i|j|k]$ where $i=1,\ldots,n_{y}$, $j=1,2$,
and $k=1\ldots,N$, so that

$\mathtt{t}^{\mathrm{right}}[i|j|k]=t^{\mathrm{right}}_{j}(L_{2}/2,y_{i}),$

both associated with the data instance $k$.

####  2.3.2 Operator construction

In this section we specify the details to build the predictive and explanatory
networks ($\mathsf{Y}$ and $\mathsf{H}$), as well as the constraint operator
$\boldsymbol{R}$. The definition of the complexity of a PGNNIV adapted to
solid mechanics stems from the architecture of the predictive and explanatory
networks, $\mathsf{Y}$ and $\mathsf{H}$, as they must be able to learn the
non-linearities between variables $\bar{\boldsymbol{t}}\mapsto\boldsymbol{U}$
or $\bar{\boldsymbol{t}}\mapsto\boldsymbol{\xi}$ and
$\boldsymbol{E}\mapsto\boldsymbol{P}$ (or
$\boldsymbol{\varepsilon}\mapsto\boldsymbol{\sigma}$), as explained in Section
2.2.2.

##### Predictive network.

The predictive network must be able to represent the data variability, so
typically it has an autoencoder-like structure. Its complexity is therefore
associated with the latent dimensionality and structure of the volumetric
loads and boundary conditions. Although more sophisticated approaches coming
from Manifold Learning theory are possible for analyzing data dimensionality
and structure [72], this is not the main interest of this work and therefore
is out of the scope of this particular study. Here we follow a much simpler
approach, where we build an ANN that is sufficiently accurate when predicting
the output $\boldsymbol{y}$ from an input $\boldsymbol{x}$.

Since we consider a biaxial quadratic load applied to the plate, the
complexity of the network depends on the data variability. For non-uniform
loads, the values of the elemental loads are the input of a DNN whose output
are the nodal displacements. For the uniform load, we use a much simpler
autoencoder-like DNN. It is important to note that a single, complex enough,
autoencoder-like DNN would be able to represent the data variability even in
the more complex scenario (the case with the largest latent space in the
process of data generation). However, we have decided to use two different
network architectures for illustrating this particular feature of PGNNIV: the
predictive network is associated to data variability, rather than data nature.
Therefore, we can adapt the network architecture to our problem
characteristics, aiming either at a better network performance (avoiding
overfitting) or at a lower computational cost.

In Appendix A, the particular architecture of the two predictive networks is
detailed, both for the non-uniform biaxial test and for the uniform biaxial
test. Anyway, although the $\mathsf{Y}$ network architecture was handcrafted,
it is expected that the more and variate data is available, the less relevant
the hand-engineering of the network becomes.

##### Explanatory network.

As in this work we restrict ourselves to the elastic regime, the input
variable is the given strain state at an arbitrary node
$\mathtt{E}[i,j|k,l|m]$ ($\mathsf{E}$ represents the Cauchy deformation tensor
for infinitesimal theory and the Green - Lagrange deformation tensor for
finite strains theory) and the output variable is the associated stress state
at that same element, $\mathtt{S}[i,j|k,l|m]$ (again, $\mathsf{S}$ represents
the Cauchy stress tensor or the first Piola-Kirchhoff tensor depending on
whether we are in the infinitesimal or finite strains theory). Note that under
the homogeneity assumption (and postulating that the stress state depends only
on the value of the deformation at the same point), the explanatory network is
a nonlinear map $\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$, due the symmetry of
both tensors, that may be expressed symbolically as

$\mathtt{S}[i,j|\cdot,\cdot|m]=\mathsf{H}\left(\mathtt{E}[i,j|\cdot,\cdot|m]\right).$

For non-local materials, given the described discretization, the explanatory
network could be in principle a map
$\mathbb{R}^{3n_{x}n_{y}}\rightarrow\mathbb{R}^{3n_{x}n_{y}}$ represented as

$\mathtt{S}[\cdot,\cdot|\cdot,\cdot|m]=\mathsf{H}\left(\mathtt{E}[\cdot,\cdot|\cdot,\cdot|m]\right).$

In particular, $\mathsf{H}$ has to be able to capture the highly nonlinear
dependencies that may exist between variables. This is in theory possible
thanks to the universal approximation theorem: by adding more internal layers
(also known as hidden layers) to the DNN model, we can provide the network
with the learning capability and complexity that a particular nonlinear
constitutive law might require. The implementation of the explanatory network
for homogeneous materials is efficiently implemented using a _convolutional_
filter to move across the domain element by element, but expanding the
features in a higher dimensional spaces, as illustrated in Fig. 2. We call
this type of architecture a moving Multilayer Perceptron (mMLP).

Figure 2: Representation of the explanatory network for the 2D plane stress
problem. On the left, the strain field computed from the displacement field
predicted by $\mathsf{Y}$ is represented on the plate. In an homogeneous
material, the DNN is fed with the strain state on each element, and the
correspondent stresses are obtained. The set of weights of this DNN moves
across the elements of the plate and updates after each iteration of the
optimization (see Figure 3 for the whole picture), acting as a
2D-convolutional filter. We call this type of architecture a moving Multilayer
Perceptron (mMLP).

For the very particular case of linear elasticity, we can explicitly
parameterize the explanatory network,
$\mathsf{H}(\boldsymbol{\sigma})=\boldsymbol{H}(\boldsymbol{\sigma};\boldsymbol{D})$
where $\boldsymbol{D}$ is the elastic tensor and
$\sigma_{ij}=D_{ijkl}\varepsilon_{kl}$. Using Voigt convention, the tensor
$\boldsymbol{D}$ is expressed, under plane stress conditions and in the most
general case, as a $3\times 3$ matrix.

$\boldsymbol{D}=\begin{pmatrix}\ d_{11}&d_{12}&d_{13}\\\ \
d_{12}&d_{22}&d_{23}\\\ \ d_{13}&d_{23}&d_{33}\\\ \end{pmatrix},$ (20)

where $d_{ij}$ are free model parameters, whereas for an isotropic elastic
material we have

$\boldsymbol{D}=\frac{E}{1-\nu^{2}}\begin{pmatrix}\ 1&\nu&0\\\ \ \nu&1&0\\\ \
0&0&1-\nu\\\ \end{pmatrix},$ (21)

where the elastic modulus $E$ and the Poisson ratio $\nu$ are the only free
fitting parameters learned during the training process.

An analogous reasoning holds for more complex parametric dependencies. It is
possible to express the explanatory network $\mathsf{H}$ as a parametric model
relating the strain and the stress states, that is

$\mathtt{S}[i,j|\cdot,\cdot|m]=\mathbf{H}\left(\mathtt{E}[i,j|\cdot,\cdot|m];\boldsymbol{\Lambda}\right).$

where $\boldsymbol{\Lambda}$ are some pre-defined fitting parameters that are
learned during the training step. In particular, an homogeneous material is
described as

$\mathtt{S}[i,j|\cdot,\cdot|m]=\mathbf{H}\left(\mathtt{E}[i,j|\cdot,\cdot|m];\boldsymbol{\Lambda}_{ij}\right).$

In this work, this approach is illustrated with different types of materials
ranging from the more simple case of a linear elastic material under
infinitesimal strain theory to an hyperelastic Ogden material under finite
strains theory.

##### Coupling the two networks using physical constraints.

The definition and subsequent formulation of the PGNNIV framework implies that
the loss function includes a term proportional to the quadratic error between
the predictions and true values of the output variable (minimization of the
maximum likelihood of the data given the parameters) and other penalty terms
related to some (physical) equations i.e. equilibrium constraints. Therefore
the different terms involved are:

  1. 1.

Loss term associated with the measurement of the displacement field:

$\mathrm{MSE}=\frac{1}{N}\sum_{i=1}^{N}||\bar{\boldsymbol{U}}^{(i)}-\mathsf{Y}(\boldsymbol{t}^{(i)})||^{2},$
(22)

where $\bar{\boldsymbol{U}}^{i}$ is the observed displacement corresponding to
sample $i$.

  2. 2.

Constraint associated with the equilibrium equation.

$\boldsymbol{\nabla}\cdot\boldsymbol{P}=\boldsymbol{0},\quad\mathrm{or}\quad\boldsymbol{\nabla}\cdot\boldsymbol{\sigma}=\boldsymbol{0}.$
(23)

  3. 3.

Constraint associated with the compatibility in the domain.

$\boldsymbol{E}-\frac{1}{2}\left(\boldsymbol{F}^{\intercal}\boldsymbol{F}-\boldsymbol{I}\right)=\boldsymbol{0},\quad\mathrm{or}\quad\boldsymbol{\varepsilon}-\frac{1}{2}(\nabla\otimes\boldsymbol{U}+\boldsymbol{U}\otimes\nabla)=\boldsymbol{0}.$
(24)

  4. 4.

Constraint associated with the equilibrium of the stresses in the boundary.

$\boldsymbol{P}\cdot\boldsymbol{N}-\boldsymbol{T}=\boldsymbol{0},\quad\mathrm{or}\quad\boldsymbol{\sigma}\cdot\boldsymbol{n}-\boldsymbol{t}=\boldsymbol{0},\text{
in }\Gamma_{N}.$ (25)

  5. 5.

Constraints associated with the compatibility of the displacements in the
boundary.

$U_{x}(x=0,y)=0,\quad U_{y}(x,y=0)=0.$ (26)

The global cost function (which turns out to be a _virtual_ physics-informed
likelihood in a Bayesian formulation or, equivalently, a regularized cost
function in the most common terminology) can be computed as a weighted sum of
$\mathrm{MSE}$ and $\mathrm{PEN}$, with $\mathrm{PEN}$ referring to the
physical terms, that is, Eqs. (23), (24), (25) and (26). As Eq. (24) may be
expressed as an explicit relation between $\boldsymbol{E}$ (or
$\boldsymbol{\varepsilon}$) and $\boldsymbol{U}$, it is directly embedded in
the network architecture. Therefore, the loss is expressed as:

$\mathrm{CF}=\mathrm{MSE}+\mathrm{PEN},$ (27)

where

$\mathrm{MSE}=\frac{1}{N}\sum_{i=1}^{N}\left[p_{1}\|\bar{\boldsymbol{U}}^{(i)}-\mathsf{Y}\left(\bar{\boldsymbol{t}}^{(i)}\right)\|^{2}\right],$
(28)

and

$\mathrm{PEN}=\frac{1}{N}\sum_{i=1}^{N}\left[p_{2}\|\boldsymbol{\nabla}\cdot\boldsymbol{\sigma}^{(i)}\|^{2}+p_{3}\|\boldsymbol{\sigma}^{(i)}\cdot\boldsymbol{n}-\bar{\boldsymbol{t}}^{(i)}\|^{2}+p_{4}\\\
\|U_{x}^{(i)}(x=0,y)+U_{y}^{(i)}(x,y=0)\|^{2}\right],$ (29)

where with the superscript $(i)$ we refer to the $i$-th piece of data and
$p_{j},j=1,2,3,4$ are penalty coefficients that account for the relative
importance of each term in the global CF (and may be seen as Lagrange
multipliers that softly enforce the constrains). Recall that no penalty for
the compatibility in the domain is included since
$\boldsymbol{E}-\frac{1}{2}\left(\boldsymbol{F}^{\intercal}\boldsymbol{F}-\boldsymbol{I}\right)$
(or
$\mathbf{\boldsymbol{\varepsilon}}-\frac{1}{2}\mathbf{(\nabla\otimes\boldsymbol{U}}+\boldsymbol{U}\otimes\nabla)$)
is identically 0.

The ANN minimization problem reads therefore:

$\min_{\boldsymbol{W}}\mathrm{CF}(\mathcal{E};\boldsymbol{W}),$ (30)

where $\boldsymbol{W}$ are the network parameters and
$\mathcal{E}=\\{\bar{\boldsymbol{t}}^{i},\bar{\boldsymbol{U}}^{i}|i=1,\cdots,N\\}$
is a given training data-set. By minimizing this function (and assuring that
not overfitting is observed by examining the predictions for test data) we
will obtain predictions of displacement, stresses and strains. For simplicity,
Algorithm 1 details a stochastic gradient descent version of the optimization,
even if in this work we always used the Adam optimizer.

From the theoretical point of view, Eq. (30) presents a complex constrained
optimization problem that has been widely studied in the context of applied
mathematics, i.e. Langrange multipliers. However, when NNs come into play
along with PDEs, the optimization becomes more involved as the complex nature
of the Pareto front, extensively studied in [73] for Physics-Informed Neural
Networks (PINNs), determines that the optimum is a state where an individual
loss cannot be further decreased without increasing at least one of the
others, and therefore the optimal set of weighting hyperparameters $p_{i}$
cannot be inferred in advance. This weighting hyperparameters $p_{i}$,
commonly referred as penalties, arise in a natural way if they are regarded as
real numbers scaling the covariance matrix of the variables’ _virtual_ maximum
likelihood probability distribution. This concept was introduced in [74] and
[75] in the context of state-space particle dynamics.

Input: PGNNIV architecture, batch size $n_{b}$, penalties $p_{k}$,
$k=1,2,3,4$, and number of iterations $M$;

Data: external forces $\bar{\boldsymbol{t}}^{(i)}$, measured displacements
$\bar{\boldsymbol{U}}^{(i)}$, $i=1,\ldots,N$;

Initialization of PGNNIV parameters, $\boldsymbol{w}=\boldsymbol{w}^{0}$,
$j=0$;

repeat

for  _$i=1,\ldots,n_{b}$_ do

$\mathtt{U}^{(i)}\leftarrow\mathsf{Y}\left(\bar{\mathtt{T}}^{(i)};\boldsymbol{w}\right)$;/*
Predictive network */

$\mathtt{E}^{(i)}=\mathsf{KIN}\left(\mathtt{U}^{(i)}\right)$;/* Green-Lagrange
or Cauchy strain tensor */

$\mathtt{S}^{(i)}\leftarrow\mathsf{H}\left(\mathtt{E}^{(i)};\boldsymbol{w}\right)$;/*
Explanatory network */

end for

$\mathrm{MSE}=\frac{1}{n_{b}}\sum_{i=1}^{n_{b}}\left[p_{1}||\bar{\mathtt{U}}^{(i)}-\mathtt{U}^{(i)}||^{2}\right]$;

$\mathrm{PEN}=\frac{1}{n_{b}}\sum_{i=1}^{n_{b}}\left[p_{2}||\mathsf{DIV}(\mathtt{S}^{(i)})||^{2}+p_{3}||\mathtt{S}^{(i)}\cdot\boldsymbol{N}-\bar{\mathtt{T}}^{(i)}||^{2}+p_{4}\left(\|\mathtt{U}_{x}^{(i)}(x=0,y)\|^{2}+\|\mathtt{U}_{y}^{(i)}(x,y=0)\|^{2}\right)\right]$;

$\mathrm{CF}=\mathrm{MSE}+\mathrm{PEN}$;

$\boldsymbol{w}\leftarrow\boldsymbol{w}-\nabla_{w}\mathrm{CF}$; /* Stochastic
gradient descent step */

$j\leftarrow j+1$;

until  _$j=M$_ ;

Output: Optimal parameters $\boldsymbol{w}^{*}=\boldsymbol{w}$ for
$\mathsf{Y}$ and $\mathsf{H}$;

Algorithm 1 PGNNIV learning algorithm

Figure 3 shows a graphical representation of the different structures involved
(tensorial fields) and the links between them (known and unknown operators)
for finite strains solid mechanics.

Figure 3: Graphical representation of the designed 2D-planar stress PGNNIV.
All significant tensorial fiels of the problem are represented: the input
variables (top and right tractions and volume forces, where the latter are
assumed to be null so removed formally from the input), the output variables
(displacement field at each nodal value), as well as the internal variables of
the problem (stress and strain fields, represented in Voigt notation).

####  2.3.3 Details about the discretization.

The discretization of both space and time domains lies on the basis of
numerical methods. In the particular case of solid mechanics under the
hypotheses considered here, time discretization turns out not to be relevant
for the overall computations, since loads are applied in a quasi-static way
and the sole discretization of the geometry provides a very good approximation
of how the continuum solid behaves.

Traditional FEMs follow a matrix-based approach to have algebraic systems
whose solution is an approximate solution. By subdividing the whole domain
into small parts (_finite elements_), PDEs governing the physical phenomena
occurring in the particular geometry can be approximated by means of
computable functions to generate algebraic systems even for complex
geometries. However, FEMs require exact knowledge of the properties of the
material and are usually time-consuming. On the contrary, PGNNIVs require no
information about the material properties since these are learned during the
training process of the network, and the calculation time for the forward
problem is reduced to seconds at prediction time in the online loop.

The discrete nature of methods such as FEM, which subdivide space in small
elements, very closely resembles that of PGNNIVs, which comprise a number of
discrete units (neurons) to represent field variables. Moreover, commonly used
differential operators in these methods are also subject of a suitable
description in the PGNNIV framework using convolutional filters.

For instance, it is possible to define the discrete gradient filter
$\mathtt{GRAD}$ acting on the nodes for obtaining values on the elements or
acting on the elements for obtaining values on the nodes. For example, if
$\mathtt{W}=\mathsf{GRAD}\otimes\mathtt{U}$, then

$\displaystyle\mathtt{W}[i,j|1,1|m]=\frac{1}{2h_{x}}\left(\Delta_{x}\mathtt{U}[i,j+1|1|m]+\Delta_{x}\mathtt{U}[i,j-1|1|m]\right),$
$\displaystyle\mathtt{W}[i,j|1,2|m]=\frac{1}{2h_{y}}\left(\Delta_{y}\mathtt{U}[i+1,j|1|m]+\Delta_{y}\mathtt{U}[i-1,j|1|m]\right),$
$\displaystyle\mathtt{W}[i,j|2,1|m]=\frac{1}{2h_{x}}\left(\Delta_{x}\mathtt{U}[i,j+1|2|m]+\Delta_{x}\mathtt{U}[i,j-1|2|m]\right),$
$\displaystyle\mathtt{W}[i,j|2,2|m]=\frac{1}{2h_{y}}\left(\Delta_{y}\mathtt{U}[i+1,j|2|m]+\Delta_{y}\mathtt{U}[i-1,j|2|m]\right),$

where
$\Delta_{x}\mathtt{U}[i,\cdot|1|m]=\mathtt{U}[i+1,\cdot|1|m]-\mathtt{U}[i-1,\cdot|1|m]$
and
$\Delta_{y}\mathtt{U}[\cdot,j|1|m]=\mathtt{U}[\cdot,j+1|1|m]-\mathtt{U}[\cdot,j-1|1|m]$.
Analogously, if $\mathtt{R}=\mathsf{GRAD}\cdot\mathtt{T}$

$\displaystyle\mathtt{R}[i,j|1|m]=\frac{1}{2h_{x}}\left(\Delta_{x}\mathtt{T}[i,j+1|1,1|m]+\Delta_{x}\mathtt{T}[i,j-1|1,1|m]\right)+\frac{1}{2h_{y}}\left(\Delta_{y}\mathtt{T}[i+1,j|1,2|m]+\Delta_{y}\mathtt{T}[i-1,j|1,2|m]\right),$
$\displaystyle\mathtt{R}[i,j|2|m]=\frac{1}{2h_{x}}\left(\Delta_{x}\mathtt{T}[i,j+1|2,1|m]+\Delta_{x}\mathtt{T}[i,j-1|2,1|m]\right)+\frac{1}{2h_{y}}\left(\Delta_{y}\mathtt{T}[i+1,j|2,2|m]+\Delta_{y}\mathtt{T}[i-1,j|2,2|m]\right),$

Now we can define the different discretized differentials. For instance the
2D-discretized symmetric gradient of a vector field $\boldsymbol{V}$ is
$\mathsf{SGRAD}(\boldsymbol{V})=\frac{1}{2}\left(\mathsf{GRAD}\otimes\boldsymbol{V}+\boldsymbol{V}\otimes\mathsf{GRAD}\right),$
where $\mathsf{GRAD}$ is the discrete gradient operator. Therefore, for large
strain we have

$\mathsf{E}=\frac{1}{2}\left(\mathsf{GRAD}\otimes\boldsymbol{U}+\boldsymbol{U}\otimes\mathsf{GRAD}+\left(\mathsf{GRAD}\otimes\boldsymbol{U}\right)\left(\boldsymbol{U}\otimes\mathsf{GRAD}\right)\right)$

and for small strains

$\mathsf{E}=\mathsf{SGRAD}(\boldsymbol{U})=\frac{1}{2}\left(\mathsf{GRAD}\otimes\boldsymbol{U}+\boldsymbol{U}\otimes\mathsf{GRAD}\right).$

Similarly, the 2D-discretized divergence of a tensor field $\boldsymbol{T}$ is
defined as:

$\mathsf{DIV}(\boldsymbol{T})=\mathsf{GRAD}\cdot\boldsymbol{T},$

so we can express the equilibrium equation both in infinitesimal as well as
finite strains theory as $\mathsf{DIV}(\boldsymbol{\sigma})=\boldsymbol{0}$ or
$\mathsf{DIV}(\boldsymbol{P})=\boldsymbol{0}$ respectively.

###  2.4 Data generation and training process

In this work we generate synthetic data although the methodology is the same
for real experimental data. Synthetic data generation is more practical, as it
allows for error validation by comparing to the real solutions, and
inexpensive, given that an accurate numerical solution is available via FEM
Analysis.

##### Small strains: linear, softening and hardening elastic materials.

For the creation of the linear, softening and hardening materials, we used
Matlab, in co-simulation with Abaqus CAE/6.14-2. Matlab was used to
automatically and iteratively generate an Abaqus input file that contained the
geometry and load profiles. Once each numerical simulation is completed, the
results are stored back into Matlab. For all the test cases, the geometry is
the same, whereas the variability in the data-set is achieved by randomly
changing the load profiles so that all the examples correspond with different
experiments. The load profiles are parabolic and are generated for both the
right and top contours.

We consider three elastic materials with different constitutive laws. On the
one hand, we choose an isotropic linear elastic material with elastic modulus
$E=1000$ Pa and Poisson’s coefficient $\nu=0.3$. On the other hand, we
considered two nonlinear materials, one with softening properties and another
one with hardening properties. In Abaqus, they were modeled as plastic
materials with no discharging effects caused by the removal of the load, and
strain ranges confined within very small values, that allowing for the
compliance of the nonlinear constitutive law with the infinitesimal strains
hypothesis. A relation of the type $\sigma=K\varepsilon^{n}$ was used, with
values of $K$ and $n$ specified in Table 1.

The data-set comprises $N=10^{3}$ FEM-simulations for the linear material and
$N=10^{4}$ FEM-simulations for the hardening and softening materials.

Material |  $K$ [Pa] |  $n$ [-]
---|---|---
Softening | $18.69$ | $0.45$
Hardening | $1.869\times 10^{12}$ | $3.5$
Table 1: Parameter values for the softening and hardening materials.

##### Finite strains: Ogden-like hyperelastic material.

For the finite strains case, an incompressible Ogden-like hyperelastic
material of order 3 is used for the data-set generation. The incompressible
Ogden hyperelastic material of order $m$ is defined in terms of its strain-
energy density function [76]:

$\Psi(\boldsymbol{C})=\Psi(\lambda_{1},\lambda_{2},\lambda_{3})=\sum_{i=p}^{m}\frac{\mu_{p}}{\alpha_{p}}\left(\lambda_{1}^{\alpha_{p}}+\lambda_{2}^{\alpha_{p}}+\lambda_{3}^{\alpha_{p}}-3\right).$
(31)

In Table 2 we report the material parameters used for the data generation. We
produce $N=10^{4}$ examples corresponding to uniform biaxial tests, where
$\lambda_{1},\lambda_{2}\in[1;1.10]$. For an incompressible membrane under
biaxial deformation, assuming a plane stress state, the solution of the
problem using an Ogden’s strain-energy function has analytical solution [77].
The displacement fields corresponding to uniform biaxial deformations are

$U_{x}(x,y)=\lambda_{1}x,\quad U_{y}(x,y)=\lambda_{2}y,\quad
U_{z}(x,y)=\frac{1}{\lambda_{1}\lambda_{2}}x,$

so, the non-vanishing components of the Green-Lagrange deformation tensor,
$\boldsymbol{E}$, are

$E_{xx}=\frac{1}{2}(\lambda_{1}^{2}-1),\quad
E_{yy}=\frac{1}{2}(\lambda_{2}^{2}-1),\quad
E_{zz}=\frac{1}{2}((\lambda_{1}\lambda_{2})^{-2}-1).$

As we are assuming plane stress, that is
$\sigma_{zz}=\sigma_{xz}=\sigma_{yz}=0$, the non-vanishing components of the
first order Piola-Kirchhoff stress tensor are

$P_{xx}=\frac{1}{\lambda_{1}}\sum_{k=1}^{3}\mu_{k}\left(\lambda_{1}^{\alpha_{k}}-(\lambda_{1}\lambda_{2})^{-\alpha_{k}}\right),\quad
P_{yy}=\frac{1}{\lambda_{2}}\sum_{k=1}^{3}\mu_{k}\left(\lambda_{2}^{\alpha_{k}}-(\lambda_{1}\lambda_{2})^{-\alpha_{k}}\right).$

Parameter | Value
---|---
$\mu_{1}$ |  $281$ Pa
$\mu_{2}$ |  $-280$ Pa
$\mu_{3}$ |  $0.31$ Pa
$\alpha_{1}$ | $1.66$
$\alpha_{2}$ | $1.61$
$\alpha_{3}$ | $38.28$
Table 2: Parameter values for the Ogden hyperelastic material.

##### Training process.

For the evaluation of the methodology, we have trained four PGNNIVs
corresponding to four cases:

  * •

Linear material with parametric explanatory network.

  * •

Linear, softening and hardening materials with non-parametric explanatory
network.

  * •

Ogden-like material with parametric explanatory network.

  * •

Ogden-like material with non-parametric explanatory network.

For all the data-sets considered, there is a number of hyperparameters that
have been tuned for obtaining the proceeding results with the different
networks, namely the learning rate $\beta$ and the four penalty coefficients
$p_{i}$, $i=1,\ldots,4.$. The specific values of these hyperparameters are
reported together with the different network topologies in Appendix A.2.

##  3 Results

When used as forward solvers, PGNNIVs can either predict measurable variables
if force-displacement data is available, for example, through Digital Image
Correlation (DIC) techniques, or explain the internal state of the solid if
this is needed for a certain application. In this section, we validate the
performance of PGNNIVs acting as a forward-solver against standard FEM
solutions for the plate using the different materials described in Section
2.4, and also as a method for constitutive equation discovery.

###  3.1 Predictive capacity

####  3.1.1 Infinitesimal strains case

We first evaluate the prediction capacity of PGNNIVs for the different tested
materials under random parabolic loads. For a quantitative evaluation of the
predictive capacity of the PGNNIV, we define the Relative Error
($\mathrm{RE}$) of an array field $\mathtt{I}$ as:

$\mathrm{RE}(\mathtt{I})=\frac{\sum_{I,J,K}\left(\hat{\mathtt{I}}[I|J|K]-\mathtt{I}[I|J|K]\right)^{2}}{\sum_{I,J,K}\mathtt{I}[I|J|K]^{2}},$
(32)

where $\hat{\mathtt{I}}$ is the predicted value and $\mathtt{I}$ the value
obtained using FEM. For instance, for the displacement field $\boldsymbol{U}$
represented by the array $\mathtt{U}$:

$\mathrm{RE}(\mathtt{U})=\sqrt{\frac{\sum_{i=1}^{n_{x}}\sum_{j=1}^{n_{y}}\sum_{k=1}^{2}\sum_{l=1}^{N}\left(\hat{\mathtt{U}}[i,j|k|l]-\mathtt{U}[i,j|k|l]\right)^{2}}{\sum_{i=1}^{n_{x}}\sum_{j=1}^{n_{y}}\sum_{k=1}^{2}\sum_{l=1}^{N}\mathtt{U}[i,j|k|l]^{2}}},$
(33)

or using standard field notation:

$\mathrm{RE}(\boldsymbol{U})=\sqrt{\frac{\sum_{i=1}^{n_{x}}\sum_{j=1}^{n_{y}}\sum_{k=1}^{2}\sum_{l=1}^{N}\left(\hat{U}^{(l)}_{k}(x_{i},y_{j})-U^{(l)}_{k}(x_{i},y_{j})\right)^{2}}{\sum_{i=1}^{n_{x}}\sum_{j=1}^{n_{y}}\sum_{k=1}^{2}\sum_{l=1}^{N}(U^{(l)}_{k}(x_{i},y_{j}))^{2}}},$
(34)

where $U_{k}^{(l)}(x,y)$ represents the value of $U_{k}(x,y)$ corresponding to
the data sample $l$.

Table 3 shows the Relative Error (RE) obtained by the PGNNIV for the
softening, hardening and linear materials. We observe predictive errors below
$10\%$ for all the cases.

Type of material |  $\mathrm{RE}(\boldsymbol{U})$ (%) |  $\mathrm{RE}(\boldsymbol{\sigma})$ (%) |  $\mathrm{RE}(\boldsymbol{\varepsilon})$ (%)
---|---|---|---
Linear | 0.53 | 3.21 | 3.61
Softening | 0.92 | 2.07 | 4.55
Hardening | 1.17 | 4.48 | 5.51
Table 3: PGNNIV predictive errors for the linear, softening and hardening
materials.

####  3.1.2 Finite strains case

We evaluate now the Ogden material under finite strains. Table 4 shows the
Relative Error (RE) obtained by the PGNNIV for the Ogden material, resulting
in that case in errors below $2\%$ (due to the lower data-set variability).

|  $\mathrm{RE}(\boldsymbol{U})$ (%) |  $\mathrm{RE}(\boldsymbol{P})$ (%) |  $\mathrm{RE}(\boldsymbol{E})$ (%)
---|---|---|---
Ogden material | 0.34 | 1.56 | 0.37
Table 4: PGNNIV predictive errors for the hyperelastic material.

###  3.2 Explanatory capacity

The explanatory capacity of the method is evaluated here in two levels:
parametric identification and state model unveiling.

####  3.2.1 Infinitesimal strains

##### Parameter identification.

First, we train and evaluate a linearized version of the PGNNIV (where
$\mathsf{Y}$ and $\mathsf{H}$ are replaced with weight matrices with no
activation functions) with the data-set corresponding to the isotropic linear
elastic material parametrized with elastic modulus $E$ and Poisson ratio
$\nu$. It is possible to evaluate the relative error when learning an elastic
tensor $\boldsymbol{D}$ by means of

$\epsilon_{r}(\boldsymbol{D})=\frac{\|\hat{\boldsymbol{D}}-\boldsymbol{D}\|_{\mathrm{F}}}{\|\boldsymbol{D}\|_{\mathrm{F}}},$
(35)

where $\hat{\boldsymbol{D}}$ is the predicted value, $\boldsymbol{D}$ the real
one, and $\|\cdot\|_{\mathrm{F}}$ depicts the Fröbenius norm,
$\|\boldsymbol{A}\|_{F}=\sqrt{\boldsymbol{A}\mathrel{\mathop{\mathchar
58\relax}}\boldsymbol{A}}$. For assessing the identifiability of the
structural parameters , it is possible to evaluate the prediction of the model
parameters associated with Eqs. (20) and (21) by computing the relative errors
for a given parameter $\lambda$, that are defined as:

$\epsilon_{r}(\lambda)=\frac{|\hat{\lambda}-\lambda|}{\lambda},$ (36)

where $\hat{\lambda}$ and $\lambda$ are the predicted and real values
respectively. The learned anisotropic and isotropic tensors, up to a precision
of $1$ Pa are:

$\boldsymbol{D}_{\mathrm{aniso}}=\begin{pmatrix}1099&330&0\\\ 330&1099&0\\\
0&0&511\\\
\end{pmatrix},\quad\boldsymbol{D}_{\mathrm{iso}}=\begin{pmatrix}1098&331&0\\\
331&1098&0\\\ 0&0&767\\\ \end{pmatrix}.$

The errors of the elastic tensor when using the anisotropic or the isotropic
elastic material model are respectively
$\varepsilon_{r}(\boldsymbol{D}_{\mathrm{aniso}})=14.4\%$ and
$\varepsilon_{r}(\boldsymbol{D}_{\mathrm{iso}})=0.2\%$, i.e when the assumed
hypotheses are true, the explanatory capacity of the network increases.
Notwithstanding, the general anisotropic model has a certain explanatory
capacity, as we may detect the supplementary structural symmetries in the
resulting elastic tensor, that is, $d_{13},d_{23}\ll
d_{11},d_{12},d_{22},d_{33}$, $d_{11}\simeq d_{22}$. The last symmetry
condition $d_{33}=d_{11}-d_{12}$ is not fulfilled, as the data-set is not very
rich in large pure shear-stress states, so the model is not able to detect
this symmetry.

Moving to specific structural parameter identification, Table 5 describes how
the model is able to predict accurately the different elastic parameters. As
commented before, the worst prediction is observed for the anisotropic model
and the parameter that correlates shear stresses and strains. If, using the
anisotropic model, we were interested in finding the value of the isotropic
model parameters $E$ and $\nu$, it is possible to compute $\nu$ using the
relations

$\nu=\frac{d_{12}}{d_{11}}=\frac{d_{12}}{d_{22}},$

and then to compute $E$ using

$\quad
E=d_{11}(1-\nu^{2})=d_{22}(1-\nu^{2})=d_{12}\frac{1-\nu^{2}}{\nu}=d_{33}(1+\nu).$

Of course, we obtain a different accuracy for each of these expressions. Using
the anisotropic model, we obtain values of $0.1\%$ and $0.1\%$ for
$\epsilon_{r}(\nu)$ and values of $0.04\%$, $0.04\%$, $0.4\%$ and $34\%$ for
$\epsilon_{r}(E)$, in agreement with the previous observations.

Parameter, $\lambda$ |  Relative error, $\epsilon_{r}(\lambda)$ (%)
---|---
Anisotropic model |
$d_{11}$ | $0.02$
$d_{12}$ | $0.08$
$d_{13}$ | $\infty$
$d_{22}$ | $0.02$
$d_{23}$ | $\infty$
$d_{33}$ | $33.59$
Isotropic model |
$E$ | $0.20$
$\nu$ | $0.52$
Table 5: Predicted and real values of the model parameters.

##### State model discovery.

While the predictive capacity of PGNNIVs does not necessarily surpass that of
a classical (unconstrained) NN, the significant improvement is visible when
assessing the explanatory capacity of the PGNNIV, which can be evaluated by
its ability to learn the material constitutive law. We perform a virtual
uniaxial test using the explanatory network, which corresponds to the
functional representation of $\sigma_{xx}=\mathsf{H}_{1}(\varepsilon,0,0)$ for
$\varepsilon\in[\varepsilon_{\mathrm{min}};\varepsilon_{\mathrm{max}}]$ and we
compare the PGNNIV predictions with a virtual uniaxial test produced with FEM,
as described in Section 2.4. Results are shown in Figure 4.

(a) Softening material.

(b) Linear material.

(c) Hardening material.

Figure 4: PGNNIV prediction versus FEM solution of the uniaxial test curve for
the different data-sets. We observe good agreement between FEM solution
(continuous line) and PGNNIV prediction (dashed line), for the softening (a),
linear (b) and hardening (c) materials.

The explanatory error is quantified as the normalized area confined between
the real uniaxial test curve and the PGNNIV-predicted one in Figure 4. It is
expressed as:

$\mathrm{RE}(\mathsf{H})=\sqrt{\frac{\int_{\varepsilon_{\mathrm{min}}}^{\varepsilon_{\mathrm{max}}}(\hat{\sigma}_{xx}(\varepsilon)-\sigma_{xx}(\varepsilon))^{2}\,\mathrm{d}\varepsilon}{\int_{\varepsilon_{\mathrm{min}}}^{\varepsilon_{\mathrm{max}}}\sigma_{xx}^{2}(\varepsilon)\,\mathrm{d}\varepsilon}},$
(37)

where $\hat{\sigma}_{xx}$ and $\sigma_{xx}$ are the predicted and FEM stresses
respectively, which result from strains
$\varepsilon\in[\varepsilon_{min};\varepsilon_{max}]$.

Type of material |  $\mathrm{RE}(\mathsf{H})$ (%)
---|---
Linear | 3.02
Softening | 0.96
Hardening | 5.92
Table 6: Explanatory errors for the $\mathsf{H}$ model subjected to uniform
uniaxial test.

####  3.2.2 Finite strains

We explore now the explanatory capacity for the finite strains case. If
$\lambda\in[\lambda_{\mathrm{min}};\lambda_{\mathrm{max}}]$ is the
longitudinal stretch, $\lambda_{1}$, the relative explanatory errors for a
given transversal stretch $\lambda_{2}$ are defined as:

$\displaystyle\mathrm{RE}_{xx}(\mathsf{H};\lambda_{2})$
$\displaystyle=\sqrt{\frac{\int_{\lambda_{\mathrm{min}}}^{\lambda_{\mathrm{max}}}(\hat{P}_{xx}(\lambda,\lambda_{2})-P_{xx}(\lambda,\lambda_{2}))^{2}\,\mathrm{d}\lambda}{\int_{\lambda_{\mathrm{min}}}^{\lambda_{\mathrm{max}}}P_{xx}^{2}(\lambda,\lambda_{2})\,\mathrm{d}\lambda}},$
(38) $\displaystyle\mathrm{RE}_{yy}(\mathsf{H};\lambda_{2})$
$\displaystyle=\sqrt{\frac{\int_{\lambda_{\mathrm{min}}}^{\lambda_{\mathrm{max}}}(\hat{P}_{yy}(\lambda,\lambda_{2})-P_{yy}(\lambda,\lambda_{2}))^{2}\,\mathrm{d}\lambda}{\int_{\lambda_{\mathrm{min}}}^{\lambda_{\mathrm{max}}}P_{yy}^{2}(\lambda,\lambda_{2})\,\mathrm{d}\lambda}}.$
(39)

Note that, as the roles of $x$ and $y$ are symmetrical in the considered
biaxial test (which is uniform, meaning that $P_{xx}$ on the top contour has
the same vale as $P_{yy}$ on the right contour), the indicated errors are
sufficient for illustrating the explanatory capacity of the method. For
structural parameter identification, the formula used for error quantification
is Eq. (36).

##### Parameter identification.

As explained previously, we first explicitly state the parametric shape of the
constitutive equation, that is, we prescribe the material to be Ogden-like.
Under these assumptions and for the uniform biaxial test considered, the
constitutive relation writes

$\displaystyle P_{xx}$
$\displaystyle=\frac{1}{\sqrt{2E_{xx}+1}}\sum_{p=1}^{3}\mu_{k}\left[\left(2E_{xx}+1\right)^{\alpha_{k}/2}-\left((2E_{xx}+1)(2E_{yy}+1)\right)^{-\alpha_{k}/2}\right],$
$\displaystyle P_{yy}$
$\displaystyle=\frac{1}{\sqrt{2E_{yy}+1}}\sum_{p=1}^{3}\mu_{k}\left[\left(2E_{yy}+1\right)^{\alpha_{k}/2}-\left((2E_{xx}+1)(2E_{yy}+1)\right)^{-\alpha_{k}/2}\right].$

Therefore, the parameters $\alpha_{k}$, $\mu_{k}$, $k=1,2,3$ are, in
principle, the ones that ought to be learned by the explanatory network. We
obtain values for the parameters of $\mu_{1}=276\,\mathrm{Pa}$,
$\mu_{2}=-277\,\mathrm{Pa}$, $\mu_{3}=0.31\,\mathrm{Pa}$, $\alpha_{1}=1.53$,
$\alpha_{2}=1.47$ and $\alpha_{3}=38.32$. The relative errors for the
different parameters are shown in Table 7. It is important to note that there
are some parameters that are more accurately predicted than others, i.e. they
are superfluous. This fact relies on the capacity of each of the parameters
for explaining the material response, as observed in Fig. 5, where we compare
the theoretical constitutive relation with the one obtained using the learned
parameters. The explanatory errors are reported in Table 8, adding evidence of
the explanatory power of the method despite the discrepancies in some
parameters.

Parameter, $\lambda$ |  Relative error, $\epsilon_{r}(\lambda)$ (%)
---|---
$\mu_{1}$ | $1.7$
$\mu_{2}$ | $-1.1$
$\mu_{3}$ | $0.5$
$\alpha_{1}$ | $7.7$
$\alpha_{2}$ | $8.6$
$\alpha_{3}$ | $0.1$
Table 7: Predicted and real values of the model parameters for the Ogden
material.

(a) Component $xx$.

(b) Component $yy$.

Figure 5: Parametric PGNNIV prediction versus analytic solution of the uniform
biaxial test curve for the Ogden material. We observe good agreement between
analytical solution (continuous line) and PGNNIV prediction (dashed line), for
the different values of $\lambda_{2}$. This indicates that the network has a
good explicability capacity even though some superflous model parameters are
not accurately fitted.

##### State model discovery.

We now evaluate the model discovered by the PGNNIV with a virtual biaxial test
using the explanatory network. The functional representation is now
$(P_{xx},P_{yy},P_{xy})=\mathsf{H}(E_{xx},E_{yy},0)$ for
$E_{xx}\in[E_{\mathrm{min}};E_{\mathrm{max}}]$ and we compare the PGNNIV
predictions with the model used for the data generation.

The results are shown in Fig. 6 for three different values of $\lambda_{2}$,
and the errors, computed according to Eqs.(39) are displayed in Table 8.

(a) Component $xx$.

(b) Component $yy$.

Figure 6: Non-parametric PGNNIV prediction versus analytic solution of the biaxial test curve for the Ogden material. We observe good agreement between FEM solution (continuous line) and PGNNIV prediction (dashed line), for the different values of $\lambda_{2}$. $Y$ stretch ratio | Parametric model | Non-parametric model
---|---|---
|  $\mathrm{RE}_{xx}(\mathsf{H};\lambda_{2})$ (%) |  $\mathrm{RE}_{yy}(\mathsf{H};\lambda_{2})$ (%) |  $\mathrm{RE}_{xx}(\mathsf{H};\lambda_{2})$ (%) |  $\mathrm{RE}_{yy}(\mathsf{H};\lambda_{2})$ (%)
$\lambda_{2}=1.02$ | 0.98 | 2.65 | 0.34 | 3.57
$\lambda_{2}=1.05$ | 1.15 | 1.93 | 0.33 | 1.77
$\lambda_{2}=1.08$ | 1.32 | 1.10 | 0.39 | 1.52
Table 8: Explanatory errors for the finite strains problem subjected to
uniform biaxial test.

##  4 Discussion, conclusions and future work

Throughout this work we have presented the mathematical foundations of PGNNIVs
in the field of computational solid mechanics, which demonstrates to be a
particularly interesting niche for the use of such methodology. We have
demonstrated that PGNNIVs have both predictive and explanatory capacity:

  * •

Predictive capacity: PGNNIVs are able to accurately predict the solid response
to new external stimuli in real time, something fundamental for optimization,
control and probabilistic problems. They are also able to predict not only the
solid response in terms of displacement field, but also the deformation and
stress fields, without the need of any extra post-processing. This has been
demonstrated in Section 3.1, where we have obtained relative errors always
below $10\%$. Controlling non-primary fields is sometimes important in
engineering problems, as high stresses cause damage, plasticity or structural
failure. As the explanatory network, once trained, encodes all the information
about the material properties, it can be used for the prediction of stresses
directly from the displacement fields, if necessary.

  * •

Explanatory capacity: PGNNIVs are able to unveil hidden state model equations,
that is, the constitutive equations of computational solid mechanics. First,
for parameter identification and fitting, PGNNIV are able to identify inherent
material symmetries (such as isotropy) and also to predict the value of the
structural model parameters with high accuracy. In the latter, PGNNIVs are in
a certain sense an alternative to conventional least-square minimization
problems [78] (e.g. using standard methods such as Levenberg-Marquardt
algorithm), but making the use of software and hardware tools associated with
ANN technology: Graphical Processor Units (GPUs) and Tensor Processor Units
(TPUs), distributed and cloud computation, scalability, transfer and federate
learning strategies among others. In addition, PGNNIVs address the more
challenging problem of model-free unravelling of nonlinear materials
constitutive laws. In Section 3.2, we have demonstrated the explanatory
capacity of PGNNIVs both for parameter identification and state model
discovery with many examples (linear and nonlinear materials both in the
infinitesimal and finite strains framework). The relative error when
predicting structural parameters is always below $2\%$ except if the data-set
does not contain information about the material response to stimuli associated
with a certain parameter or the parameter itself has not a direct impact in
the explanatory capacity (superfluous parameters). Besides, the explanatory
relative error is below $10\%$ for all the cases analyzed.

One important characteristic of PGNNIVs related to this double capacity is the
fact that it is possible to decouple two sources of variability in the data
obtained from a mechanical system that can be measured using any kind of
sensor, namely, the stimuli variability and the variability related to system
response.

  * •

The stimuli variability is stored in the predictive network, which acts as an
autoencoder. The encoder is able to map data that lives in a space of
dimension $D$, to a latent space of dimension $d\ll D$. The size of the latent
space is therefore informative about the data variability, and the values of
the latent variables are a compressed representation of data. In addition to
theoretical considerations, this fact has important practical consequences: if
we know the sources of variability of our system, it is possible to design the
predictive network accordingly.

  * •

The physical interpretable knowledge, that is, the constitutive equation of
the material, is delocalized, spread and diluted in the weights of the
predictive network, but is also encoded in the relation
$\boldsymbol{\varepsilon}\mapsto\boldsymbol{\sigma}$ or
$\boldsymbol{E}\mapsto\boldsymbol{P}$ that is learned by the explanatory
network, in a much more structured way. Indeed, if it is intended to identify
and adjust some structural parameters, this is possible. Otherwise, the
constitutive model as a whole may be also unravelled using the expressiveness
of ANN methods.

In that sense, PGNNIVs are knowledge generators from data as most of ML
techniques are, with the difference that for this particular method, the
physical knowledge is directly distilled in a separated component. This
particularity is what makes the difference between PINNs and PGNNIVs. In
PINNs, data and mathematical physics models are seamlessly integrated for
solving parametric PDEs of a given problem [42, 68], but, by construction, the
information cannot be extrapolated to other situations. That means that, in
the context of solid mechanics, the network trained for a given problem using
PINNs cannot be used for predicting the response of the system under different
volume loads or boundary conditions, which greatly weakens its ability as a
predictive method. PGNNIVs overcome this difficulty precisely by distilling
the physical information from the intrinsic variability of the stimuli.

There is another paramount characteristic of PGNNIVs for computational solid
mechanics that has been largely discussed and has explained their emergence.
There is no need to have access to the values of internal variables, that are,
strictly speaking, non-measurable as they are mathematical constructs coming
from a scientific theory. In that sense, even if the bases for
thermodynamically-appropriated ANN for constitutive equations in solid
mechanics have been investigated [79], it is important to recall that stress
fields, such as $\boldsymbol{P}$ or $\boldsymbol{\sigma}$ are not accessible
without the need of extra hypothesis (geometry or load specific
configurations). This fundamental issue should not be overlooked and many
recent works have worked for this purpose. Efficient Unsupervised Constitutive
Law Identification and Discovery (EUCLID) method is one of the most acclaimed
efforts in that direction, either using sparse identification [80, 45, 81],
clustering [82], Bayesian methods [83] or ANN [47]. However, EUCLID paradigm
relies on the fact that the geometry and loads are appropriate enough to
ensure strain-stress fields variability in a single specimen. When the
geometry or the data acquisition capabilities do not satisfy this requirement,
the only possibility is to take action on the data-set or the network, or at
least to track the different load conditions and incorporate them to the
computational pipeline.

Finally, among the different PIML methods, PGNNIVs are more transparent than
other approaches that have demonstrated to be very performant for computing
the evolution of dynamical systems by incorporating thermodynamical
constraints. Structure Preserving Neural Networks enforce first and second
laws of thermodynamics as regularization term [84]. Even if in the cited work
the GENERIC structure allows for a split between reversible and non-reversible
(dissipative) components, the physical information is again diluted in all the
network weights, rather than in some specific components. A dimensionality
reduction of the dynamic data was also explored [85]. This may be interpreted
also as an information reduction to distil physical knowledge, although the
interpretability is still dark. In PGNNIVs, however, interpretable physical
information (that is, knowledge) is located in specific ANN components.

Nevertheless, the presented methodology still has some limitations and there
exists room for exploration in several directions:

  1. 1.

The data requirements for the problem in hands are high. In this work, we have
generated data synthetically, but in reality sensors usually collect noisy
data from experimental tests and there exist also important limitations
concerning the size of the data-sets. A probabilistic (Bayesian) viewpoint
will enable a new interpretation of PGNNIVs in the _small data_ regime,
although this methodology is rather thought for systems where intensive data
mining is possible, in which data quantity prevails over data quality.

  2. 2.

Defining a suitable architecture for the $\mathsf{Y}$ and $\mathsf{H}$
networks is not a simple task and requires an iterative process, including the
tuning of many hyperparameters. In addition, PGNNIVs require extra hyper-
parameters, i.e. penalty coefficients $p_{i}$, making the process even more
involved and time consuming. We have presented some insights into to the
complexity and architecture of both predictive and explanatory networks,
related with both their prediction and explanation character, but either a
great intuition for network design, or time-intensive trial-and-error
iterative testing are needed.

  3. 3.

In this work, we have made used of relatively coarse discretizations, but
finer meshes will result in more expensive training processes due to the
exponentially larger number of parameters required. More powerful
computational strategies (distributed computing, parallelization) as well as
more advanced hardware (GPUs and TPUs) will enable the acceleration of the
PGNNIVs training processes although the problem of dealing with
multidimensional and unstructured meshes still remains open.

Many challenges lie ahead of the development of a more general PGNNIV
framework. Next lines of research will address the formulation of PGNNIVs
under finite strain assumptions using a much more theoretical basis, such as
the presented in some recent works [79, 86, 87], leveraging its predictive and
explanatory power. Furthermore, extensions of the 2D planar stress
architecture to general 3D problems with more complex geometries and load
scenarios as well as to more complex constitutive laws that might depend on
time (visco-elasticity), or heterogeneous conditions, still pose major
challenges for the future.

Finally, the usage of real data from sensors, for example, through Digital
Image or Volume Correlation tests (DIC, DVC) [71] or even more advanced
methods such as Finite Element Model Updating (FEMU) [88] or Virtual Fields
Methods (VFM) [89] will put to the test the applicability of PGNNIVs to real
scientific problems in the field of engineering.

In conclusion, we have demonstrated that in the context of computational solid
mechanics, PGNNIVs are a family of ANN for accurately predicting measurable
and non-measurable variables such as displacement and stress fields, in real
time, and are also able to describe or unravel the constitutive model with
high accuracy for different linear and nonlinear (hyper-)elastic materials.
Even if this work is preliminary, the ingredients that it comprises correspond
to a general approach and the methodology can be applied to cases of
scientific interest with the necessary adaptations of the network
architectures.

## References

  * [1] David Reinsel-John Gantz-John Rydning, John Reinsel, and John Gantz.  The digitization of the world from edge to core.  Framingham: International Data Corporation, 16:1–28, 2018.
  * [2] George Morahan.  The 20 most valuable companies in market capitalisation, 2022.  Accessed: 2022-09-30.
  * [3] Christopher M Bishop and Nasser M Nasrabadi.  Pattern recognition and machine learning, volume 4.  Springer, 2006.
  * [4] Luigi Atzori, Antonio Iera, and Giacomo Morabito.  The internet of things: A survey.  Computer networks, 54(15):2787–2805, 2010.
  * [5] Jeyan Thiyagalingam, Mallikarjun Shankar, Geoffrey Fox, and Tony Hey.  Scientific machine learning benchmarks.  Nature Reviews Physics, 4(6):413–420, 2022.
  * [6] Adnan Nuhic, Tarik Terzimehic, Thomas Soczka-Guth, Michael Buchholz, and Klaus Dietmayer.  Health diagnosis and remaining useful life prognostics of lithium-ion batteries using data-driven methods.  Journal of power sources, 239:680–688, 2013.
  * [7] Li C Xue, Drena Dobbs, Alexandre MJJ Bonvin, and Vasant Honavar.  Computational prediction of protein interfaces: A review of data driven methods.  FEBS letters, 589(23):3516–3526, 2015.
  * [8] Oliver Lemon and Olivier Pietquin.  Data-driven methods for adaptive spoken dialogue systems: Computational learning for conversational interfaces.  Springer Science & Business Media, 2012.
  * [9] Rob Kitchin.  Big data, new epistemologies and paradigm shifts.  Big data & society, 1(1):2053951714528481, 2014.
  * [10] Frank Rosenblatt.  The perceptron: A probabilistic model for information storage and organization in the brain.  Psychological Review, 65(6):386, 1958.
  * [11] George Cybenko.  Approximations by superpositions of a sigmoidal function.  Mathematics of Control, Signals and Systems, 2:183–192, 1989.
  * [12] Kurt Hornik.  Approximation capabilities of multilayer feedforward networks.  Neural Networks, 4(2):251–257, 1991.
  * [13] Allan Pinkus.  Approximation theory of the mlp model in neural networks.  Acta Numerica, 8:143–195, 1999.
  * [14] Zhou Lu, Hongming Pu, Feicheng Wang, Zhiqiang Hu, and Liwei Wang.  The expressive power of neural networks: A view from the width.  In Advances in Neural Information Processing Systems, pages 6231–6239, 2017.
  * [15] Boris Hanin.  Universal function approximation by deep neural nets with bounded width and relu activations.  Mathematics, 7(10):992, 2019.
  * [16] Yves Chauvin and David E Rumelhart.  Backpropagation: theory, architectures, and applications.  Psychology press, 1995.
  * [17] Shun-ichi Amari.  Backpropagation and stochastic gradient descent method.  Neurocomputing, 5(4-5):185–196, 1993.
  * [18] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  Advances in neural information processing systems, 25, 2012.
  * [19] Alex Graves, Abdel-rahman Mohamed, and Geoffrey Hinton.  Speech recognition with deep recurrent neural networks.  In 2013 IEEE international conference on acoustics, speech and signal processing, pages 6645–6649. Ieee, 2013.
  * [20] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [21] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  Advances in neural information processing systems, 30, 2017.
  * [22] Timm Teubner, Christoph M Flath, Christof Weinhardt, Wil van der Aalst, and Oliver Hinz.  Welcome to the era of chatgpt et al. the prospects of large language models.  Business & Information Systems Engineering, 65(2):95–101, 2023\.
  * [23] Michael Schmidt and Hod Lipson.  Distilling free-form natural laws from experimental data.  Science, 324(5923):81–85, 2009.
  * [24] Steven L Brunton, Joshua L Proctor, and J Nathan Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the National Academy of Sciences, 113(15):3932–3937, 2016.
  * [25] Silviu-Marian Udrescu and Max Tegmark.  Ai feynman: a physics-inspired method for symbolic regression.  Science Advances, 6(16):eaay2631, 2020.
  * [26] Francisco Chinesta and Elias Cueto.  Empowering engineering with data, machine learning and artificial intelligence: a short introductive review.  Advanced Modeling and Simulation in Engineering Sciences, 9(1):21, 2022.
  * [27] Trenton Kirchdoerfer and Michael Ortiz.  Data-driven computational mechanics.  Computer Methods in Applied Mechanics and Engineering, 304:81–101, 2016.
  * [28] Jacobo Ayensa-Jiménez, Mohamed H Doweidar, Jose A Sanz-Herrera, and Manuel Doblaré.  An unsupervised data completion method for physically-based data-driven models.  Computer Methods in Applied Mechanics and Engineering, 344:120–143, 2019.
  * [29] Elias Cueto and Francisco Chinesta.  Thermodynamics of learning physical phenomena.  Archives of Computational Methods in Engineering, pages 1–14, 2023\.
  * [30] Samuel Greydanus, Misko Dzamba, and Jason Yosinski.  Hamiltonian neural networks.  Advances in neural information processing systems, 32, 2019.
  * [31] Pengzhan Jin, Zhen Zhang, Aiqing Zhu, Yifa Tang, and George Em Karniadakis.  Sympnets: Intrinsic structure-preserving symplectic networks for identifying hamiltonian systems.  Neural Networks, 132:166–179, 2020.
  * [32] Xue Ying.  An overview of overfitting and its solutions.  In Journal of physics: Conference series, volume 1168, page 022022\. IOP Publishing, 2019.
  * [33] Aravindh Mahendran and Andrea Vedaldi.  Understanding deep image representations by inverting them.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 5188–5196, 2015.
  * [34] Ravid Shwartz-Ziv and Naftali Tishby.  Opening the black box of deep neural networks via information.  arXiv preprint arXiv:1703.00810, 2017.
  * [35] Amina Adadi and Mohammed Berrada.  Peeking inside the black-box: a survey on explainable artificial intelligence (xai).  IEEE access, 6:52138–52160, 2018.
  * [36] Feiyu Xu, Hans Uszkoreit, Yangzhou Du, Wei Fan, Dongyan Zhao, and Jun Zhu.  Explainable ai: A brief survey on history, research areas, approaches and challenges.  In Natural Language Processing and Chinese Computing: 8th CCF International Conference, NLPCC 2019, Dunhuang, China, October 9–14, 2019, Proceedings, Part II 8, pages 563–574. Springer, 2019.
  * [37] Ruben Ibañez, Domenico Borzacchiello, Jose Vicente Aguado, Emmanuelle Abisset-Chavanne, Elias Cueto, Pierre Ladeveze, and Francisco Chinesta.  Data-driven non-linear elasticity: constitutive manifold construction and problem discretization.  Computational Mechanics, 60:813–826, 2017.
  * [38] Jacobo Ayensa-Jimenez, Mohamed H Doweidar, Jose A Sanz-Herrera, and Manuel Doblare.  Prediction and identification of physical systems by means of physically-guided neural networks with meaningful internal layers.  Computer Methods in Applied Mechanics and Engineering, 381:113816, 2021.
  * [39] Filippo Masi, Ioannis Stefanou, Paolo Vannucci, and Victor Maffi-Berthier.  Thermodynamics-based artificial neural networks for constitutive modeling.  Journal of the Mechanics and Physics of Solids, 147:104277, 2021\.
  * [40] Anuj Karpatne, Gowtham Atluri, James H Faghmous, Michael Steinbach, Arindam Banerjee, Auroop Ganguly, Shashi Shekhar, Nagiza Samatova, and Vipin Kumar.  Theory-guided data science: A new paradigm for scientific discovery from data.  IEEE Transactions on knowledge and data engineering, 29(10):2318–2331, 2017.
  * [41] Nikhil Muralidhar, Mohammad Raihanul Islam, Manish Marwah, Anuj Karpatne, and Naren Ramakrishnan.  Incorporating prior domain knowledge into deep neural networks.  In 2018 IEEE international conference on big data (big data), pages 36–45. IEEE, 2018.
  * [42] Maziar Raissi, Paris Perdikaris, and George 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.
  * [43] Jacobo Ayensa-Jiménez, Mohamed H Doweidar, Jose A Sanz-Herrera, and Manuel Doblare.  Understanding glioblastoma invasion using physically-guided neural networks with internal variables.  PLoS Computational Biology, 18(4):e1010019, 2022.
  * [44] Isaac E Lagaris, Aristidis C Likas, and Dimitris G Papageorgiou.  Neural-network methods for boundary value problems with irregular boundaries.  IEEE Transactions on Neural Networks, 11(5):1041–1049, 2000.
  * [45] Moritz Flaschel, Siddhant Kumar, and Laura De Lorenzis.  Discovering plasticity models without stress data.  npj Computational Materials, 8(1):91, 2022.
  * [46] Alexandre M Tartakovsky, Carlos Ortiz Marrero, Paris Perdikaris, Guzel D Tartakovsky, and David Barajas-Solano.  Learning parameters and constitutive relationships with physics informed deep neural networks.  arXiv preprint arXiv:1808.03398, 2018.
  * [47] Prakash Thakolkaran, Akshay Joshi, Yiwen Zheng, Moritz Flaschel, Laura De Lorenzis, and Siddhant Kumar.  Nn-euclid: Deep-learning hyperelasticity without stress data.  Journal of the Mechanics and Physics of Solids, 169:105076, 2022\.
  * [48] B Ph van Milligen, V Tribaldos, and JA Jiménez.  Neural network differential equation and plasma equilibrium solver.  Physical Review Letters, 75(20):3594, 1995.
  * [49] Isaac E Lagaris, Aristidis Likas, and Dimitrios I Fotiadis.  Artificial neural networks for solving ordinary and partial differential equations.  IEEE Transactions on Neural Networks, 9(5):987–1000, 1998.
  * [50] Keith Rudd and Silvia Ferrari.  A constrained integration (cint) approach to solving partial differential equations using artificial neural networks.  Neurocomputing, 155:277–285, 2015.
  * [51] Justin Sirignano and Konstantinos Spiliopoulos.  Dgm: A deep learning algorithm for solving partial differential equations.  Journal of Computational Physics, 375:1339–1364, 2018.
  * [52] Diab W Abueidda, Qiyue Lu, and Seid Koric.  Deep learning collocation method for solid mechanics: Linear elasticity, hyperelasticity, and plasticity as examples.  arXiv e-prints, pages arXiv–2012, 2020.
  * [53] Bing Yu et al.  The deep ritz method: a deep learning-based numerical algorithm for solving variational problems.  Communications in Mathematics and Statistics, 6(1):1–12, 2018.
  * [54] Esteban Samaniego, Cosmin Anitescu, Somdatta Goswami, Vien Minh Nguyen-Thanh, Hongwei Guo, Khader Hamdia, X Zhuang, and T Rabczuk.  An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications.  Computer Methods in Applied Mechanics and Engineering, 362:112790, 2020.
  * [55] Vien Minh Nguyen-Thanh, Xiaoying Zhuang, and Timon Rabczuk.  A deep energy method for finite deformation hyperelasticity.  European Journal of Mechanics-A/Solids, 80:103874, 2020.
  * [56] Anuj Karpatne, William Watkins, Jordan Read, and Vipin Kumar.  Physics-guided neural networks (pgnn): An application in lake temperature modeling.  arXiv preprint arXiv:1710.11431, 2, 2017.
  * [57] Russell Stewart and Stefano Ermon.  Label-free supervision of neural networks with physics and domain knowledge.  In Thirty-First AAAI Conference on Artificial Intelligence, pages 2576–2582, 2017.
  * [58] Jim Magiera, Deep Ray, Jan S Hesthaven, and Christian Rohde.  Constraint-aware neural networks for riemann problems.  Journal of Computational Physics, 409:109345, 2020.
  * [59] Weinan E, Jiequn Han, and Arnulf Jentzen.  Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations.  Communications in Mathematics and Statistics, 5(4):349–380, 2017\.
  * [60] Jiequn Han, Arnulf Jentzen, and Weinan E.  Solving high-dimensional partial differential equations using deep learning.  Proceedings of the National Academy of Sciences, 115(34):8505–8510, 2018.
  * [61] Côme Huré, Huyên Pham, and Xavier Warin.  Deep backward schemes for high-dimensional nonlinear pdes.  Mathematics of Computation, 89(324):1547–1579, 2020.
  * [62] Hamid Reza Tamaddon-Jahromi, Neeraj Kavan Chakshu, Igor Sazonov, Llion M Evans, Hywel Thomas, and Perumal Nithiarasu.  Data-driven inverse modelling through neural network (deep learning) and computational heat transfer.  Computer Methods in Applied Mechanics and Engineering, 369:113217, 2020.
  * [63] Craig M Hamel, Kevin N Long, and Sharlotte LB Kramer.  Calibrating constitutive models with full-field data via physics informed neural networks.  Strain, page e12431, 2022.
  * [64] Jacobo Ayensa-Jiménez, Mohamed H Doweidar, Jose A Sanz-Herrera, and Manuel Doblaré.  On the application of physically-guided neural networks with internal variables to continuum problems.  arXiv preprint arXiv:2011.11376, 2020.
  * [65] Konstantinos Karapiperis, Laurent Stainier, Michael Ortiz, and José E Andrade.  Data-driven multiscale modeling in mechanics.  Journal of the Mechanics and Physics of Solids, 147:104239, 2021\.
  * [66] Chengping Rao, Hao Sun, and Yang Liu.  Physics-informed deep learning for computational elastodynamics without labeled data.  Journal of Engineering Mechanics, 147(8):04021043, 2021.
  * [67] BVSS Bharadwaja, Mohammad Amin Nabian, Bharatkumar Sharma, Sanjay Choudhry, and Alankar Alankar.  Physics-informed machine learning and uncertainty quantification for mechanics of heterogeneous materials.  Integrating Materials and Manufacturing Innovation, pages 1–21, 2022.
  * [68] George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang.  Physics-informed machine learning.  Nature Reviews Physics, 3(6):422–440, 2021.
  * [69] Michael A Nielsen.  Neural networks and deep learning, volume 25.  Determination press San Francisco, CA, 2015.
  * [70] Steven Walczak and Narciso Cerpa.  Heuristic principles for the design of artificial neural networks.  Information and software technology, 41(2):107–117, 1999.
  * [71] TC Chu, WF Ranson, and Michael A Sutton.  Applications of digital-image-correlation techniques to experimental mechanics.  Experimental mechanics, 25(3):232–244, 1985.
  * [72] Yunqian Ma and Yun Fu.  Manifold learning theory and applications, volume 434.  CRC press Boca Raton, 2012.
  * [73] Franz M Rohrhofer, Stefan Posch, and Bernhard C Geiger.  On the pareto front of physics-informed neural networks.  arXiv preprint arXiv:2105.00862, 2021.
  * [74] Sebastian Kaltenbach and Phaedon-Stelios Koutsourelakis.  Incorporating physical constraints in a deep probabilistic machine learning framework for coarse-graining dynamical systems.  Journal of Computational Physics, 419:109673, 2020.
  * [75] Sebastian Kaltenbach and Phaedon-Stelios Koutsourelakis.  Physics-aware, deep probabilistic modeling of multiscale dynamics in the small data regime.  arXiv preprint arXiv:2102.04269, 2021.
  * [76] Raymond William Ogden.  Large deformation isotropic elasticity–on the correlation of theory and experiment for incompressible rubberlike solids.  Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 326(1567):565–584, 1972.
  * [77] Gerhard A Holzapfel.  Nonlinear solid mechanics: a continuum approach for engineering science, 2002.
  * [78] John E Dennis Jr and Robert B Schnabel.  Numerical methods for unconstrained optimization and nonlinear equations.  SIAM, 1996.
  * [79] Kevin Linka and Ellen Kuhl.  A new family of constitutive artificial neural networks towards automated model discovery.  Computer Methods in Applied Mechanics and Engineering, 403:115731, 2023.
  * [80] Moritz Flaschel, Siddhant Kumar, and Laura De Lorenzis.  Unsupervised discovery of interpretable hyperelastic constitutive laws.  Computer Methods in Applied Mechanics and Engineering, 381:113852, 2021.
  * [81] Moritz Flaschel, Siddhant Kumar, and Laura De Lorenzis.  Automated discovery of generalized standard material models with euclid.  Computer Methods in Applied Mechanics and Engineering, 405:115867, 2023.
  * [82] Enzo Marino, Moritz Flaschel, Siddhant Kumar, and Laura De Lorenzis.  Automated identification of linear viscoelastic constitutive laws with euclid.  Mechanics of Materials, 181:104643, 2023.
  * [83] Akshay Joshi, Prakash Thakolkaran, Yiwen Zheng, Maxime Escande, Moritz Flaschel, Laura De Lorenzis, and Siddhant Kumar.  Bayesian-euclid: Discovering hyperelastic material laws with uncertainties.  Computer Methods in Applied Mechanics and Engineering, 398:115225, 2022.
  * [84] Quercus Hernández, Alberto Badías, David González, Francisco Chinesta, and Elías Cueto.  Structure-preserving neural networks.  Journal of Computational Physics, 426:109950, 2021.
  * [85] Quercus Hernandez, Alberto Badias, David Gonzalez, Francisco Chinesta, and Elias Cueto.  Deep learning of thermodynamics-aware reduced-order models from data.  Computer Methods in Applied Mechanics and Engineering, 379:113763, 2021.
  * [86] Kevin Linka, Sarah R St Pierre, and Ellen Kuhl.  Automated model discovery for human brain using constitutive artificial neural networks.  Acta Biomaterialia, 160:134–151, 2023.
  * [87] Lennart Linden, Dominik K Klein, Karl A Kalina, Jörg Brummund, Oliver Weeger, and Markus Kästner.  Neural networks meet hyperelasticity: A guide to enforcing physics.  Journal of the Mechanics and Physics of Solids, page 105363, 2023\.
  * [88] John E Mottershead, Michael Link, and Michael I Friswell.  The sensitivity method in finite element model updating: A tutorial.  Mechanical systems and signal processing, 25(7):2275–2296, 2011\.
  * [89] Michel Grediac, Fabrice Pierron, Stéphane Avril, and Evelyne Toussaint.  The virtual fields method for extracting constitutive parameters from full-field measurements: a review.  Strain, 42(4):233–253, 2006.

##  Appendix A Details about network architectures and hyperparameters.

###  A.1 Network architectures used

As explained in Section 2.3.2, we use different combinations of architectures
for both predictive and explanatory networks, $\mathsf{Y}$ and $\mathsf{H}$,
depending on the problem. Table 9 shows the architecture details for the
unique PGNNIV trained on the linear, softening and hardening materials in
small strains theory when the parametric structure of the constitutive
equation is not prescribed. For the anisotropic and isotropic parametric
discovery, the explanatory network is substituted by the direct parametric
expression. The $\mathsf{Y}$ network layer includes empty elements (null
value) corresponding to the most general natural boundary conditions where
shear loads are applied on the contour. Therefore, the input layer comprises
$2\times(n_{x}-1)\times(n_{y}-1)=2\times(10+8)=36$ neurons. The number of
neurons of the output layer corresponds to $2\times n_{x}\times n_{y}=2\times
11\times 9=198$.

Network | Class of architecture | Detailed layers (input and output neurons included)
---|---|---
Predictive network $\mathsf{Y}$ | MLP | 36/50/tanh/50/tanh/50/tanh/50/tanh/198
Explanatory network $\mathsf{H}$ | mMLP | 3/100/tanh/100/tanh/100/tanh/100/tanh/100/tanh/100/tanh/100/tanh/3
Table 9: Architecture of the PGNNIV for the softening, hardening and linear
materials.

Table 10 shows the implementation details for the unique PGNNIV trained on the
Ogden material in finite strain theory when the parametric structure of the
constitutive equation is not prescribed. For the Ogden constitutive model
parametric discovery, the explanatory network is substituted by the direct
parametric expression.

Network | Class of architecture | Detailed layers (input and output neurons included)
---|---|---
Predictive network $\mathsf{Y}$ | autoencoder | 18/9/tanh/3/tanh/3/tanh/9/tanh/198
Explanatory network $\mathsf{H}$ | mMLP | 3/60/ReLU/60/ReLU/60/ReLU/60/ReLU/60/ReLU/3
Table 10: Architecture of the PGNNIV for the Ogden material.

###  A.2 Hyperparameters

For reproducibility purposes, Table 11 shows a summary of the hyperparameters
that were found to give optimal results for the optimization process.

| Loss function coefficients | Learning algorithm
---|---|---
| $p_{1}$ | $p_{2}$ | $p_{3}$ | $p_{4}$ |  $\alpha$ (Learn. rate) |  $n_{b}$ (batch size)
Linear (model-free) | $1\times 10^{3}$ | $1\times 10^{-3}$ | $1\times 10^{-3}$ | $1\times 10^{3}$ | $5\times 10^{-5}$ | $100$
Linear (anisotropic model) | $1\times 10^{3}$ | $1\times 10^{-3}$ | $1\times 10^{-3}$ | $1\times 10^{3}$ | $5\times 10^{-5}$ | $100$
Linear (isotropic model) | $1\times 10^{3}$ | $1\times 10^{-3}$ | $1\times 10^{-3}$ | $1\times 10^{3}$ | $5\times 10^{-5}$ | $100$
Softening (model-free) | $1\times 10^{3}$ | $1\times 10^{-3}$ | $1\times 10^{-3}$ | $1\times 10^{3}$ | $5\times 10^{-5}$ | $100$
Hardening (model-free) | $1\times 10^{4}$ | $1\times 10^{-8}$ | $1\times 10^{-8}$ | $1\times 10^{4}$ | $5\times 10^{-4}$ | $100$
Ogden (model-free) | $1.0\times 10^{0}$ | $1\times 10^{-3}$ | $1\times 10^{-3}$ | $1\times 10^{-3}$ | $1\times 10^{-3}$ | $128$
Ogden (model-based) | $1.0\times 10^{0}$ | $1\times 10^{-3}$ | $1\times 10^{-3}$ | $1\times 10^{-3}$ | $1\times 10^{-3}$ | $128$
Table 11: Value of the hyperparameters for all the training processes.

##  Appendix B Some field solution predictions of the PGNNIV for the linear,
softening, hardening and Ogden materials.

As an illustration of the predictive capacity of the methodology, in this
Appendix we illustrate some solution fields associated to the linear,
softening and hardening material problems, as well as for the Ogden-like
material.

###  B.1 Infinitesimal theory.

####  B.1.1 Linear material.

In Fig. 7 the displacement field (FEM solution versus PGNNIV prediction) is
represented for a test example. In Fig. 8 the components of the strain tensor
are illustrated and in Fig. 9 the components of the stress tensor. We observe
high similarity between the ground truth and predictive fields despite the
coarse discretization.

(a) Real horizontal component ($u_{x}$).

(b) Predicted horizontal component ($u_{x}$).

(c) Real vertical component ($u_{y}$).

(d) Predicted vertical component ($u_{y}$).

Figure 7: PGNNIV prediction versus FEM solution of the components of the
displacement field for a single test-set example of the linear material.

(a) Real normal component ($\varepsilon_{xx}$).

(b) Predicted normal component ($\varepsilon_{xx}$).

(c) Real normal component ($\varepsilon_{yy}$).

(d) Predicted normal component ($\varepsilon_{yy}$).

(e) Real shear component ($\varepsilon_{xy}$).

(f) Predicted shear component ($\varepsilon_{xy}$).

Figure 8: PGNNIV prediction versus FEM solution of the components of the
strain fields for a single test-set example of the linear material.

(a) Real normal component ($\sigma_{xx}$).

(b) Predicted normal component ($\sigma_{xx}$).

(c) Real normal component ($\sigma_{yy}$).

(d) Predicted normal component ($\sigma_{yy}$).

(e) Predicted normal component ($\sigma_{xy}$).

(f) Predicted shear component ($\sigma_{xy}$).

Figure 9: PGNNIV prediction versus FEM solution of the components of the
stress fields for a single test-set example of the linear material.

####  B.1.2 Softening material.

In Fig. 10 the displacement field (FEM solution versus PGNNIV prediction) is
represented for a test example. In Fig. 11 the components of the strain tensor
are illustrated and in Fig. 12 the components of the stress tensor. We observe
high similarity between the ground truth and predictive fields despite the
coarse discretization.

(a) Real horizontal component ($u_{x}$).

(b) Predicted horizontal component ($u_{x}$).

(c) Real vertical component ($u_{y}$).

(d) Predicted vertical component ($u_{y}$).

Figure 10: PGNNIV prediction versus FEM solution of the components of the
displacement field for a single test-set example of the softening material.

(a) Real normal component ($\varepsilon_{xx}$).

(b) Predicted normal component ($\varepsilon_{xx}$).

(c) Real normal component ($\varepsilon_{yy}$).

(d) Predicted normal component ($\varepsilon_{yy}$).

(e) Real shear component ($\varepsilon_{xy}$).

(f) Predicted shear component ($\varepsilon_{xy}$).

Figure 11: PGNNIV prediction versus FEM solution of the components of the
strain fields for a single test-set example of the softening material.

(a) Real normal component ($\sigma_{xx}$).

(b) Predicted normal component ($\sigma_{xx}$).

(c) Real normal component ($\sigma_{yy}$).

(d) Predicted normal component ($\sigma_{yy}$).

(e) Real shear component ($\sigma_{xy}$).

(f) Predicted shear component ($\sigma_{xy}$).

Figure 12: PGNNIV prediction versus FEM solution of the components of the
stress fields for a single test-set example of the softening material.

####  B.1.3 Hardening material.

In Fig. 13 the displacement field (FEM solution versus PGNNIV prediction) is
represented for a test example. In Fig. 14 the components of the strain tensor
are illustrated and in Fig. 15 the components of the stress tensor. We observe
high similarity between the ground truth and predictive fields despite the
coarse discretization.

(a) Real horizontal component ($u_{x}$).

(b) Predicted horizontal component ($u_{x}$).

(c) Real vertical component ($u_{y}$).

(d) Predicted vertical component ($u_{y}$).

Figure 13: PGNNIV prediction versus FEM solution of the components of the
displacement field for a single test-set example of the hardening material.

(a) Real normal component ($\varepsilon_{xx}$).

(b) Predicted normal component ($\varepsilon_{xx}$).

(c) Real normal component ($\varepsilon_{yy}$).

(d) Predicted normal component ($\varepsilon_{yy}$).

(e) Real shear component ($\varepsilon_{xy}$).

(f) Predicted shear component ($\varepsilon_{xy}$).

Figure 14: PGNNIV prediction versus FEM solution of the components of the
strain fields for a single test-set example of the hardening material.

(a) Real normal component ($\sigma_{xx}$).

(b) Predicted normal component ($\sigma_{xx}$).

(c) Real normal component ($\sigma_{yy}$).

(d) Predicted normal component ($\sigma_{yy}$).

(e) Real shear component ($\sigma_{xy}$).

(f) Predicted shear component ($\sigma_{xy}$).

Figure 15: PGNNIV prediction versus FEM solution of the components of the
stress fields for a single test-set example of the hardening material.

###  B.2 Finite strains.

In Fig. 16 the displacement field (analytic solution versus PGNNIV prediction)
is represented for a test example. In Fig. 17 the components of the strain
tensor are illustrated and in Fig. 18 the components of the stress tensor. We
observe high similarity between the ground truth and predictive fields despite
the coarse discretization.

(a) Real horizontal component ($u_{x}$).

(b) Predicted horizontal component ($u_{x}$).

(c) Real vertical component ($u_{y}$).

(d) Predicted vertical component ($u_{y}$).

Figure 16: Nonlinear PGNNIV prediction versus FEM solution of the components
of the displacement field for a single test-set example of the Ogden material.

(a) Real normal component ($E_{xx}$).

(b) Predicted normal component ($E_{xx}$).

(c) Real normal component ($E_{yy}$).

(d) Predicted normal component ($E_{yy}$).

(e) Real shear component ($E_{xy}$).

(f) Predicted shear component ($E_{xy}$).

Figure 17: PGNNIV prediction versus FEM solution of the components of the
strain fields for a single test-set example of the Ogden material.

(a) Real normal component ($P_{xx}$).

(b) Predicted normal component ($P_{xx}$).

(c) Real normal component ($P_{yy}$).

(d) Predicted normal component ($P_{yy}$).

(e) Real shear component ($P_{xy}$).

(f) Predicted hear component ($P_{xy}$).

Figure 18: PGNNIV prediction versus FEM solution of the components of the
stress fields for a single test-set example of the Ogden material.