Predicting and explaining nonlinear material response using deep
  Physically Guided Neural Networks with Internal Variables
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).
Javier Orera-Echeverria, Jacobo Ayensa-Jiménez, Manuel Doblare
Predicting and explaining nonlinear material response using deep Physically Guided Neural Networks with Internal Variables

###### 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).

Nonlinear computational solid mechanics Deep Neural Network Internal Variables Explainable Artificial Intelligence Physics-Informed Machine Learning 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 XX\({}^{\text{th}}\) 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

### 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_).

### 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

\[\mathcal{F}(u,v) =f,\text{ in }\Omega, \tag{1a}\] \[\mathcal{G}(u,v) =g,\text{ in }\partial\Omega,\] (1b) \[\mathcal{H}(u) =v,\text{ in }\Omega, \tag{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 \(\mathbf{F}\), \(\mathbf{G}\) and \(\mathbf{H}\) and nodal values \(\mathbf{u}\), \(\mathbf{v}\), \(\mathbf{f}\) and \(\mathbf{g}\). Particularly, \(\mathbf{u}\) are the solution field nodal values and \(\mathbf{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

\[\mathbf{F}(\mathbf{u},\mathbf{v})=\mathbf{f},\text{ in }\Omega, \tag{2a}\] \[\mathbf{G}(\mathbf{u},\mathbf{v})=\mathbf{g},\text{ in }\partial\Omega,\] (2b) \[\mathbf{H}(\mathbf{u})=\mathbf{v},\text{ in }\Omega. \tag{2c}\]

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

\[\mathbf{y} =\mathsf{Y}(\mathbf{x}),\] \[\mathbf{v} =\mathsf{H}(\mathbf{u}),\] \[\mathbf{x} =\mathbf{I}(\mathbf{u},\mathbf{f},\mathbf{g}),\] \[\mathbf{y} =\mathbf{O}(\mathbf{u},\mathbf{f},\mathbf{g}),\] \[\mathbf{R}(\mathbf{u},\mathbf{v},\mathbf{f},\mathbf{g})=0,\]

where:

1. \(\mathbf{u}\), \(\mathbf{f}\) and \(\mathbf{g}\) are the measurable variables of the problem.
2. \(\mathbf{x}\) and \(\mathbf{y}\) are the input and output variables respectively and will be defined depending on which relation \(\mathbf{x}\mapsto\mathbf{y}\) wants to be predicted.
3. \(\mathbf{I}\) and \(\mathbf{O}\) are functions that compute the input \(\mathbf{x}\) and the output \(\mathbf{y}\) of the problem from the measurable variables. In other words, functions \(\mathbf{I}\) and \(\mathbf{O}\) define the data used as starting point to make predictions, \(\mathbf{x}\), and the data that we want to predict, that is, \(\mathbf{y}\).
4. \(\mathbf{R}\) are the physical constraints, related to the relations given by \(\mathbf{F}\) and \(\mathbf{G}\).
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 \(\mathbf{x}\mapsto\mathbf{y}\).
* \(\mathsf{H}\) is the **explanatory model**, whose objective is to unravel the hidden physics of the relation \(\mathbf{u}\mapsto\mathbf{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}(\mathbf{\sigma})+\rho\mathbf{b}=\mathbf{0}, \tag{4}\]

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

\[\mathrm{DIV}(\mathbf{P})+\rho\mathbf{B}=\mathbf{0}, \tag{5}\]

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

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

\[\mathbf{F}=\mathrm{GRAD}(\chi)=\mathrm{GRAD}\otimes\mathbf{\xi}, \tag{6}\]

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

\[\mathrm{ROT}(\mathbf{F})=\mathbf{0}, \tag{7}\]

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

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

\[\mathbf{\sigma}=\frac{1}{\det(\mathbf{F})}\frac{\partial\Psi}{\partial\mathbf{F}}\cdot\bm {F}^{\intercal}, \tag{8}\]
where \(\Psi\) is the strain energy function expressed as a function of the deformation state given by \(\mathbf{F}\), \(\Psi=\mathfrak{F}(\mathbf{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\):

\[\mathbf{\xi}=\bar{\mathbf{\xi}}, \tag{9}\]

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

\[\mathbf{\sigma}\cdot\mathbf{n}=\bar{\mathbf{t}}, \tag{10}\]

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

\[\mathbf{P}\cdot\mathbf{N}=\bar{\mathbf{T}}, \tag{11}\]

where now \(\mathbf{N}\) and \(\bar{\mathbf{T}}\) are the material analogous of \(\mathbf{n}\) and \(\bar{\mathbf{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 \(\mathbf{\xi}=\chi(\mathbf{\Xi})\), that is, \(\mathbf{\xi}\) is a measurable variable. Let us also assume that we can measure the volumetric loads, \(\mathbf{b}=\mathbf{b}(\mathbf{\xi})\), as well as the prescribed traction forces \(\bar{\mathbf{t}}\) (or, in an equivalent manner, \(\mathbf{B}=\mathbf{B}(\mathbf{\Xi})\) and \(\bar{\mathbf{T}}\)). Therefore, using Eqs. (6) and (8) it is possible to express:

\[\mathbf{F} =\mathcal{A}(\mathbf{\xi}), \tag{12a}\] \[\mathbf{\sigma} =\mathcal{B}(\Psi,\mathbf{\xi}), \tag{12b}\]

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

\[\mathrm{div}(\mathbf{\sigma}(\Psi,\mathbf{\xi})) =-\rho\mathbf{b}\quad\mathrm{in}\quad\Omega, \tag{13a}\] \[\mathbf{\xi} =\bar{\mathbf{\xi}}\quad\mathrm{in}\quad\Gamma_{E},\] (13b) \[\mathbf{\sigma}\cdot\mathbf{n} =\bar{\mathbf{t}}\quad\mathrm{in}\quad\Gamma_{N},\] (13c) \[\Psi =\mathfrak{F}(\mathbf{F}(\mathbf{\xi}))\quad\mathrm{in}\quad\Omega. \tag{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

\[\mathcal{F}(\mathbf{\xi},\mathbf{\sigma})=\mathbf{b},\;\mathrm{in}\;\Omega, \tag{14a}\] \[\mathcal{G}(\mathbf{\xi},\mathbf{\sigma})=(\bar{\mathbf{\xi}},\bar{\mathbf{t}}) \;\mathrm{in}\;\partial\Omega,\] (14b) \[\mathcal{H}(\mathbf{\xi})=\mathbf{\sigma},\;\mathrm{in}\;\Omega, \tag{14c}\]

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

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

\[\mathbf{\sigma}=\mathfrak{G}(\mathbf{\varepsilon}), \tag{15}\]

where \(\mathbf{\varepsilon}\) is the Cauchy small strain tensor

\[\mathbf{\varepsilon}=\mathrm{symgrad}(\mathbf{U})=\frac{1}{2}(\mathbf{J}+\mathbf{J}^{\intercal }), \tag{16}\]

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

\[\mathcal{F}(\mathbf{U},\mathbf{\sigma}) =\mathbf{b},\;\mathrm{in}\;\Omega, \tag{17a}\] \[\mathcal{G}(\mathbf{U},\mathbf{\sigma}) =(\bar{\mathbf{U}};\bar{\mathbf{t}}),\;\mathrm{in}\;\partial\Omega,\] (17b) \[\mathcal{H}(\mathbf{U}) =\mathbf{\sigma},\;\mathrm{in}\;\Omega, \tag{17c}\]

again in the form of Eqs. (1) with \(u=\mathbf{U}\), \(v=\mathbf{\sigma}\), \(f=\mathbf{b}\) and \(g=(\bar{\mathbf{U}},\bar{\mathbf{t}})\).
Once we have discretized the problem, our aim is to predict a motion field \(\mathbf{\xi}\) (or a displacement field \(\mathbf{U}\)) from a particular load case, expressed in terms of the volumetric loads and the natural boundary conditions1, \(\bar{\mathbf{t}}\), therefore \(\mathbf{I}(\mathbf{u},\mathbf{f},\mathbf{g})=(\mathbf{f},\mathbf{g})=(\mathbf{b},\bar{\mathbf{t}})\). With these last remarks, the PGNNIV problem is stated for finite solid mechanics as

Footnote 1: It 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.

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

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

\[\mathbf{U} =\mathsf{Y}(\bar{\mathbf{t}}) \tag{19a}\] \[\mathbf{\sigma} =\mathsf{H}(\mathrm{symgrad}(\mathbf{U}))\] (19b) \[\mathbf{x} =(\bar{\mathbf{t}},\mathbf{b})\] (19c) \[\mathbf{y} =\mathbf{U}\] (19d) \[\mathbf{R}(\mathbf{U},\mathbf{\sigma},\bar{\mathbf{t}}) =(\mathrm{div}(\mathbf{\sigma})-\rho\mathbf{b};\mathbf{\sigma}\cdot\mathbf{n}- \bar{\mathbf{t}};\mathbf{U}-\bar{\mathbf{U}}). \tag{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, \(\mathbf{b}=\mathbf{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.

### 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 \(\mathsf{I}\) containing the information of the tensorial field itself. The different dimensions of the data are represented using the indexation \(\mathsf{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 \(\mathsf{U}[i,j|k|l]\) where \(i=1,\dots,n_{x}\), \(j=1,\dots,n_{y}\), and where \(k=1,2\) (2D problem) and \(l=1\dots,N\), so that

\[\mathsf{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 \(\mathsf{E}[i,j|k,l|m]\) and \(\mathsf{S}[i,j|k,l|m]\) where \(i=1,\dots,n_{x}-1\), \(j=1,\dots,n_{y}-1\), \(k,l=1,2\) and \(m=1\dots,N\), such that, for instance,

\[\mathsf{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, \(\mathbf{\bar{t}}^{\rm top}\) and \(\mathbf{\bar{t}}^{\rm right}\), which are treated as the inputs of our problem (provided the volume forces are not considered), are represented by \(\mathtt{t}^{\rm top}[i|j|k]\) where \(i=1,\ldots,n_{x}\), \(j=1,2\), and \(k=1\ldots,N\), so that

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

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

\[\mathtt{t}^{\rm right}[i|j|k]=t_{j}^{\rm right}(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 (Y and H), as well as the constraint operator \(\mathbf{R}\). The definition of the complexity of a PGNNIV adapted to solid mechanics stems from the architecture of the predictive and explanatory networks, Y and H, as they must be able to learn the non-linearities between variables \(\mathbf{\bar{t}}\mapsto\mathbf{U}\) or \(\mathbf{\bar{t}}\mapsto\mathbf{\xi}\) and \(\mathbf{E}\mapsto\mathbf{P}\) (or \(\mathbf{\varepsilon}\mapsto\mathbf{\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 \(\mathbf{y}\) from an input \(\mathbf{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.

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.
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 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 \(\mathbb{E}[i,j|k,l|m]\) (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, \(\mathbb{S}[i,j|k,l|m]\) (again, 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}\to\mathbb{R}^{3}\), due the symmetry of both tensors, that may be expressed symbolically as

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

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

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

In particular, 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).

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

\[\mathbf{D}=\left(\begin{array}{ccc}d_{11}&d_{12}&d_{13}\\ d_{12}&d_{22}&d_{23}\\ d_{13}&d_{23}&d_{33}\end{array}\right), \tag{20}\]

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 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).
where \(d_{ij}\) are free model parameters, whereas for an isotropic elastic material we have

\[\mathbf{D}=\frac{E}{1-\nu^{2}}\left(\begin{array}{ccc}1&\nu&0\\ \nu&1&0\\ 0&0&1-\nu\end{array}\right), \tag{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

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

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

\[\mathsf{S}[i,j]\cdot,\cdot|m]=\mathbf{H}\left(\mathsf{E}[i,j]\cdot,\cdot|m]; \mathbf{\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. Loss term associated with the measurement of the displacement field: \[\mathrm{MSE}=\frac{1}{N}\sum_{i=1}^{N}||\bar{\mathbf{U}}^{(i)}-\mathsf{Y}(\mathbf{t}^ {(i)})||^{2},\] (22) where \(\bar{\mathbf{U}}^{i}\) is the observed displacement corresponding to sample \(i\).
2. Constraint associated with the equilibrium equation. \[\mathbf{\nabla}\cdot\mathbf{P}=\mathbf{0},\quad\mathrm{or}\quad\mathbf{\nabla}\cdot\mathbf{\sigma }=\mathbf{0}.\] (23)
3. Constraint associated with the compatibility in the domain. \[\mathbf{E}-\frac{1}{2}\left(\mathbf{F}^{\intercal}\mathbf{F}-\mathbf{I}\right)=\mathbf{0},\quad \mathrm{or}\quad\mathbf{\varepsilon}-\frac{1}{2}(\nabla\otimes\mathbf{U}+\mathbf{U} \otimes\nabla)=\mathbf{0}.\] (24)
4. Constraint associated with the equilibrium of the stresses in the boundary. \[\mathbf{P}\cdot\mathbf{N}-\mathbf{T}=\mathbf{0},\quad\mathrm{or}\quad\mathbf{\sigma}\cdot\mathbf{n}- \mathbf{t}=\mathbf{0},\;\mathrm{in}\;\Gamma_{N}.\] (25)
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 \(\mathbf{E}\) (or \(\mathbf{\varepsilon}\)) and \(\mathbf{U}\), it is directly embedded in the network architecture. Therefore, the loss is expressed as:

\[\mathrm{CF}=\mathrm{MSE}+\mathrm{PEN}, \tag{27}\]

where

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

and

\[\mathrm{PEN}=\frac{1}{N}\sum_{i=1}^{N}\left[p_{2}\|\mathbf{\nabla}\cdot\mathbf{ \sigma}^{(i)}\|^{2}+p_{3}\|\mathbf{\sigma}^{(i)}\cdot\mathbf{n}-\bar{\mathbf{t}}^{(i)}\|^ {2}+p_{4}\|U_{x}^{(i)}(x=0,y)+U_{y}^{(i)}(x,y=0)\|^{2}\right], \tag{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 \(\mathbf{E}-\frac{1}{2}\left(\mathbf{F}^{\intercal}\mathbf{F}-\mathbf{I}\right)\) (or \(\mathbf{\varepsilon}-\frac{1}{2}(\nabla\otimes\mathbf{U}+\mathbf{U}\otimes\nabla)\)) is identically \(\mathbf{0}\).

The ANN minimization problem reads therefore:

\[\min_{\mathbf{W}}\mathrm{CF}(\mathcal{E};\mathbf{W}), \tag{30}\]

where \(\mathbf{W}\) are the network parameters and \(\mathcal{E}=\{\mathbf{\tilde{t}}^{i},\mathbf{\bar{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 \(\mathbf{\tilde{t}}^{(i)}\), measured displacements \(\mathbf{\bar{U}}^{(i)}\), \(i=1,\ldots,N\);  Initialization of PGNNIV parameters, \(\mathbf{w}=\mathbf{w}^{0}\), \(j=0\); repeat for\(i=1,\ldots,n_{b}\)do \(\mathtt{U}^{(i)}\leftarrow\mathsf{Y}\left(\mathbf{\bar{t}}^{(i)};\mathbf{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)};\mathbf{w}\right)\);/* Explanatory network */  end for \(\mathrm{MSE}=\frac{1}{n_{b}}\sum_{i=1}^{n_{b}}\left[p_{1}||\mathbf{\bar{0}}^{(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\mathbf{N}-\mathbf{\bar{T}}^{(i)} ||^{2}+p_{4}\left(||\mathtt{U}^{(i)}_{x}(x=0,y)||^{2}+||\mathtt{U}^{(i)}_{y}(x,y=0)||^{2}\right)\right]\); \(\mathrm{CF}=\mathrm{MSE}+\mathrm{PEN}\); \(\mathbf{w}\leftarrow\mathbf{w}-\nabla_{w}\mathrm{CF}\); /* Stochastic gradient descent step */ \(j\gets j+1\); until\(j=M\); Output: Optimal parameters \(\mathbf{w}^{*}=\mathbf{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.
#### 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.

Figure 3: **Graphical representation of the designed 2D-planar stress PGNNIV.** All significant tensorial fields 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).
For instance, it is possible to define the discrete gradient filter 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

\[\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),\] \[\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),\] \[\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),\] \[\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}\)

\[\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),\] \[\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 \(\mathbf{V}\) is \(\mathsf{SGRAD}(\mathbf{V})=\frac{1}{2}\left(\mathsf{GRAD}\otimes\mathbf{V}+\mathbf{V} \otimes\mathsf{GRAD}\right),\) where GRAD is the discrete gradient operator. Therefore, for large strain we have

\[\mathsf{E}=\frac{1}{2}\left(\mathsf{GRAD}\otimes\mathbf{U}+\mathbf{U}\otimes\mathsf{ GRAD}+\left(\mathsf{GRAD}\otimes\mathbf{U}\right)\left(\mathbf{U}\otimes\mathsf{GRAD} \right)\right)\]

and for small strains

\[\mathsf{E}=\mathsf{SGRAD}(\mathbf{U})=\frac{1}{2}\left(\mathsf{GRAD}\otimes\mathbf{U}+ \mathbf{U}\otimes\mathsf{GRAD}\right).\]

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

\[\mathsf{DIV}(\mathbf{T})=\mathsf{GRAD}\cdot\mathbf{T},\]

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

### 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.
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(\mathbf{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). \tag{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, \(\mathbf{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).\]

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.

\begin{table}
\begin{tabular}{|l|l|l|} \hline
**Material** & \(K\) [Pa] & \(n\) [-] \\ \hline Softening & \(18.69\) & \(0.45\) \\ \hline Hardening & \(1.869\times 10^{12}\) & \(3.5\) \\ \hline \end{tabular}
\end{table}
Table 1: Parameter values for the softening and hardening materials.

\begin{table}
\begin{tabular}{|c|c|} \hline
**Parameter** & **Value** \\ \hline \(\mu_{1}\) & \(281\) Pa \\ \hline \(\mu_{2}\) & \(-280\) Pa \\ \hline \(\mu_{3}\) & \(0.31\) Pa \\ \hline \(\alpha_{1}\) & \(1.66\) \\ \hline \(\alpha_{2}\) & \(1.61\) \\ \hline \(\alpha_{3}\) & \(38.28\) \\ \hline \end{tabular}
\end{table}
Table 2: Parameter values for the Ogden hyperelastic material.
_Predicting and explaining nonlinear material response using Deep Physically Guided Neural Networks with Internal Variables_. J. Orera-Echeverria et al.

### 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 \(\mathrm{I}\) as:

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

where \(\hat{\mathrm{I}}\) is the predicted value and \(\mathrm{I}\) the value obtained using FEM. For instance, for the displacement field \(\mathbf{U}\) represented by the array \(\mathrm{U}\):

\[\mathrm{RE}(\mathrm{\text{\text{\text{\text{\text{\text{\text{\text{ \text{\text{\text{\text{\text{\
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}, \tag{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:

\[\mathbf{D}_{\rm aniso}=\begin{pmatrix}1099&330&0\\ 330&1099&0\\ 0&0&511\end{pmatrix},\quad\mathbf{D}_{\rm 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}(\mathbf{D}_{\rm aniso})=14.4\%\) and \(\varepsilon_{r}(\mathbf{D}_{\rm 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

\[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.

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_{\rm min};\varepsilon_{\rm 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.

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:

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

\begin{table}
\begin{tabular}{|c|c|} \hline
**Parameter, \(\lambda\)** & **Relative error, \(\epsilon_{r}(\lambda)\) (\%)** \\ \hline Anisotropic model & \\ \hline \(d_{11}\) & \(0.02\) \\ \(d_{12}\) & \(0.08\) \\ \(d_{13}\) & \(\infty\) \\ \(d_{22}\) & \(0.02\) \\ \(d_{23}\) & \(\infty\) \\ \(d_{33}\) & \(33.59\) \\ \hline \hline Isotropic model & \\ \hline \(E\) & \(0.20\) \\ \hline \(\nu\) & \(0.52\) \\ \hline \end{tabular}
\end{table}
Table 5: Predicted and real values of the model parameters.
where \(\hat{\sigma}_{xx}\) and \(\sigma_{xx}\) are the predicted and FEM stresses respectively, which result from strains \(\varepsilon\in[\varepsilon_{min};\varepsilon_{max}]\).

#### 3.2.2 Finite strains

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

\[\mathrm{RE}_{xx}(\mathsf{H};\lambda_{2}) =\sqrt{\frac{\int_{\lambda_{\min}}^{\lambda_{\max}}(\hat{P}_{xx} (\lambda,\lambda_{2})-P_{xx}(\lambda,\lambda_{2}))^{2}\,\mathrm{d}\lambda}{ \int_{\lambda_{\min}}^{\lambda_{\max}}P_{xx}^{2}(\lambda,\lambda_{2})\, \mathrm{d}\lambda}}, \tag{38}\] \[\mathrm{RE}_{yy}(\mathsf{H};\lambda_{2}) =\sqrt{\frac{\int_{\lambda_{\min}}^{\lambda_{\max}}(\hat{P}_{yy} (\lambda,\lambda_{2})-P_{yy}(\lambda,\lambda_{2}))^{2}\,\mathrm{d}\lambda}{ \int_{\lambda_{\min}}^{\lambda_{\max}}P_{yy}^{2}(\lambda,\lambda_{2})\, \mathrm{d}\lambda}}. \tag{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

\begin{table}
\begin{tabular}{|l|l|} \hline
**Type of material** & \(\mathrm{RE}(\mathsf{H})\) (\%) \\ \hline Linear & 3.02 \\ \hline Softening & 0.96 \\ \hline Hardening & 5.92 \\ \hline \end{tabular}
\end{table}
Table 6: Explanatory errors for the \(\mathsf{H}\) model subjected to uniform uniaxial test.

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.**
considered, the constitutive relation writes \[P_{xx} =\frac{1}{\sqrt{2E_{xx}+1}}\sum_{p=1}^{3}\mu_{k}\left[(2E_{xx}+1)^{ \alpha_{k}/2}-\left((2E_{xx}+1)(2E_{yy}+1)\right)^{-\alpha_{k}/2}\right],\] \[P_{yy} =\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.

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.

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

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 superfluous model parameters are not accurately fitted.
## 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

\begin{table}
\begin{tabular}{|c|c|c|c|c|} \hline \(Y\) stretch ratio & \multicolumn{2}{c|}{**Parametric model**} & \multicolumn{2}{c|}{**Non-parametric model**} \\ \hline  & \(\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})\) (\%) \\ \hline \(\lambda_{2}=1.02\) & 0.98 & 2.65 & 0.34 & 3.57 \\ \hline \(\lambda_{2}=1.05\) & 1.15 & 1.93 & 0.33 & 1.77 \\ \hline \(\lambda_{2}=1.08\) & 1.32 & 1.10 & 0.39 & 1.52 \\ \hline \end{tabular}
\end{table}
Table 8: Explanatory errors for the finite strains problem subjected to uniform biaxial test.

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}\).**
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 \(\varepsilon\mapsto\mathbf{\sigma}\) or \(\mathbf{E}\mapsto\mathbf{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 \(\mathbf{P}\) or \(\mathbf{\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. 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. Defining a suitable architecture for the Y and 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. 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.
