# Automated discovery of generalized standard material models with EUCLID

Moritz Flaschel  Siddhant Kumar  Laura De Lorenzis  Department of Mechanical
and Process Engineering, ETH Zürich, 8092 Zürich, Switzerland  Department of
Materials Science and Engineering, Delft University of Technology, 2628 CD
Delft, The Netherlands

###### Abstract

We extend the scope of our recently developed approach for unsupervised
automated discovery of material laws (denoted as EUCLID) to the general case
of a material belonging to an unknown class of constitutive behavior. To this
end, we leverage the theory of generalized standard materials, which
encompasses a plethora of important constitutive classes including elasticity,
viscosity, plasticity and arbitrary combinations thereof. We show that, based
only on full-field kinematic measurements and net reaction forces, EUCLID is
able to automatically discover the two scalar thermodynamic potentials,
namely, the Helmholtz free energy and the dissipation potential, which
completely define the behavior of generalized standard materials. The a priori
enforced constraint of convexity on these potentials guarantees by
construction stability and thermodynamic consistency of the discovered model;
balance of linear momentum acts as a fundamental constraint to replace the
availability of stress-strain labeled pairs; sparsity promoting regularization
enables the automatic selection of a small subset from a possibly large number
of candidate model features and thus leads to a parsimonious, i.e., simple and
interpretable, model. Importantly, since model features go hand in hand with
the correspondingly active internal variables, sparse regression automatically
induces a parsimonious selection of the few internal variables needed for an
accurate but simple description of the material behavior. A fully automatic
procedure leads to the selection of the hyperparameter controlling the weight
of the sparsity promoting regularization term, in order to strike a user-
defined balance between model accuracy and simplicity. By testing the method
on synthetic data including artificial noise, we demonstrate that EUCLID is
able to automatically discover the true hidden material model from a large
catalog of constitutive classes, including elasticity, viscoelasticity,
elastoplasticity, viscoplasticity, isotropic and kinematic hardening.

###### keywords:

unsupervised learning, automated discovery, constitutive models, generalized
standard materials, interpretable models, sparse regression, inverse problems

##  Section 1 Introduction

The conventional strategy to mathematically describe the response of a given
material in solid mechanics entails the a priori choice of a material model
depending on a finite number of tunable parameters, and the calibration of
these parameters based on experimental observations. The calibrated material
model, which describes the relation between stresses and strains (and
potentially other material state variables), can then be deployed in numerical
(e.g. finite element) simulations to predict the mechanical response of
arbitrarily shaped components of the material under external influences. The
weakest aspect of this procedure lies in the a priori selection of a suitable
model, which is largely based on experience. An inappropriate initial choice
inevitably leads to failure to identify a set of material parameters which
enable the model to accurately represent the experimental data; in such case
the initial choice has to be modified and the procedure repeated, possibly
multiple times, resulting in a time-consuming and error-prone iterative
process.

Enabled by the tremendous amount of data made available by the recent advances
in experimental mechanics, data-driven approaches, often based on machine
learning tools, are increasingly being explored for applications in solid
mechanics, and especially in material modeling (Neggers et al. (2018); Montáns
et al. (2019); Liu et al. (2021); Peng et al. (2021); Kumar and Kochmann
(2022)). A popular strategy to mitigate the drawbacks of an inappropriate a
priori model selection is to choose a parametric material model ansatz that is
as general as possible, i.e., with a large number of tunable parameters.
Examples are splines (Sussman and Bathe, 2009), Gaussian processes (Frankel et
al., 2020), neural networks (Ghaboussi et al., 1991) or neural ordinary
differential equations (Tac et al., 2022). Due to the large number of
trainable parameters and functions that describe the material response, these
types of model ansatz are more flexible than conventional models and the risk
of not being able to well interpret the observed data is reduced. Combinations
of traditional and machine learning components can further increase the
performance of the resulting material models, as shown for example by González
et al. (2019a); Ibáñez et al. (2019); Fuhg et al. (2021, 2022). Recognizing
that constitutive models should not violate any well-known physical
requirements motivated the recent trend to supplement machine learning models
with physics-based constraints. These are enforced either by construction or
in a weak sense by introducing additional regularization terms to the loss
function that penalize the violation of physical laws. Examples of constraints
include convexity or polyconvexity in hyperelasticity (As’ad et al., 2022;
Klein et al., 2022a, b; Kalina et al., 2022; Tac et al., 2022; Thakolkaran et
al., 2022) or the second law of thermodynamics for dissipative materials (Masi
et al., 2021; Vlassis and Sun, 2021; Fuhg et al., 2022; Huang et al., 2022;
Masi and Stefanou, 2022). The mentioned works show that the enforcement of
physics constraints has the additional advantage of improving the
extrapolation power of machine-learning-based models. An alternative stream of
research avoids the fitting of a parametric material model altogether and aims
to solve boundary value problems in which the material model is replaced by
experimental data, resulting in so-called (material) model-free data-driven
approaches (Kirchdoerfer and Ortiz (2016); Ibañez et al. (2018)). Here the
modeling components are limited to the momentum balance equation and the
assumed kinematics, and thus are reduced to a minimum. The same problem can be
formulated in an inverse setting to recover stress fields from kinematic
measurements, see Leygue et al. (2018); Dalémat et al. (2019); Stainier et al.
(2019).

Most of the previously discussed methods have two weak points in common.
Firstly, both the training procedure for the machine learning models and the
construction of a material database for the model-free methods require a large
number of labeled stress-strain data pairs, which are not available from
standard experiments and can only be generated through computationally
expensive multiscale simulations (assuming the lower-scale geometry and
material behavior are known). The second weakness is the black box nature of
these methods, which makes it impossible to interpret their predictions with
mathematical or physical domain knowledge.

Against the outlined background and inspired by research on discovery of
ordinary or partial differential equations based on data in the physics
community (Schmidt and Lipson, 2009; Brunton et al., 2016), we recently
developed a new approach for automated discovery of material models, which we
denoted as EUCLID (Efficient Unsupervised Constitutive Law Identification and
Discovery). Originally developed for hyperelastic materials (Flaschel et al.,
2021) (see also the independent related work by Wang et al. (2021)) and
subsequently extended to elastoplasticity (Flaschel et al., 2022) and
viscoelasticity (Marino et al., in preparation), EUCLID solves both of the
aforementioned issues. In contrast to most machine-learning-based methods,
EUCLID relies solely on data that can be realistically acquired with common
experimental techniques, namely, full-field displacement data obtained, e.g.,
by Digital Image Correlation (DIC), and net reaction forces. In this sense,
EUCLID is related to parameter calibration methods based on full-field
displacement measurements such as the Virtual Fields Method (Grédiac, 1989;
Pierron and Grédiac, 2012), see also Hild and Roux (2006); Avril et al.
(2008); Roux and Hild (2020) for an overview. However, instead of calibrating
the parameters of an a priori assumed material model, EUCLID automatically
selects the model that leads to the best interpretation of the data starting
from a very large set of candidates. Leveraging sparse regression (Frank and
Friedman, 1993; Tibshirani, 1996; Brunton et al., 2016), it is ensured that a
parsimonious and interpretable material model, i.e., one with a small number
of terms and parameters, is discovered. Later works on EUCLID include a study
of the problem from a Bayesian perspective for quantifying the uncertainty in
the discovered models (Joshi et al., 2022) and an unsupervised training
framework for artificial neural networks (Thakolkaran et al., 2022). For
related works that aim to identify material parameters of an a priori known
constitutive model or train uninterpretable black-box machine learning models
in an unsupervised fashion, the reader is referred to works by Man and
Furukawa (2011); Huang et al. (2020); Liu et al. (2020); Amores et al. (2022);
Anton and Wessels (2022).

Thus far, we developed EUCLID for a few different material classes, including
hyperelasticity (Flaschel et al., 2021), elastoplasticity (Flaschel et al.,
2022)) and viscoelasticity (Marino et al., in preparation). For each material
class, we formulate a general parametric model library within the class, from
which the model is automatically selected (and its unknown parameters
simultaneously calibrated) based on sparse regression. This strategy is
appropriate if the material class is known a priori, which may be realistic in
many practical scenarios. On the other hand, in more complex situations it may
well happen that the studied material is completely or largely unknown, so
that the constitutive material class it belongs to is open to question. Thus,
the objective of this work is to extend the scope of EUCLID to the challenging
task of automatically discovering the material class along with the particular
model within this class, as usual on the basis of full-field displacement and
net force data. To pursue this goal, we leverage the theoretical framework of
generalized standard materials.

The theory of generalized standard materials (sometimes denoted as standard
dissipative materials) provides a general material modeling framework rooted
in the fundamental results by Onsager (1931a, b) on the thermodynamics of
irreversible processes. Onsager departs from a classical description of the
kinetics of irreversible processes by a linear system of equations and proves
the symmetry of its coefficient matrix, expressed by the so-called Onsager’s
reciprocal relations. This property in turn enables the definition of a
quadratic dissipation potential from which the linear kinetic evolution
equations can be obtained, and the establishment of a variational principle,
recently also denoted as Onsager’s variational principle by the soft matter
community (Doi, 2011). Onsager’s thermodynamic framework is extended to
material modeling by Ziegler (1957, 1958, 1972) and Halphen and Nguyen (1975)
(among many others), with Halphen and Nguyen (1975) introducing the notion of
generalized standard materials. Following these studies, it emerges that
constitutive relations for many types of materials can be deduced from a
single pair of scalar functions (thermodynamic potentials) fully
characterizing the material behavior: the free energy and the dissipation
potential. Unlike in the work by Onsager (1931a, b), Ziegler (1957) and
Halphen and Nguyen (1975) consider thermodynamic potentials that are not
necessarily quadratic, and hence may lead to non-linear relations between
thermodynamic forces and state variable rates, but convex, which is shown to
guarantee the fulfillment of thermodynamic consistency by construction. This
generality enables the description of a large variety of different material
classes with possibly non-linear constitutive equations, including elasticity,
plasticity, viscosity, damage and more (see Steinmann and Runesson (2021) for
a recent catalogue in the one-dimensional case).

Due to the general description of the material behavior and the automatic
fulfillment of physical laws, the theory of generalized standard materials
constitutes an attractive framework for developing data-driven and machine-
learning-based methods. It is hence surprising that thus far little effort has
been made in this direction. One exception is the work by Yu et al. (2021),
who use neural networks to parameterize the system matrices, the energy and
the external forces occurring in a generalized form of Onsager’s reciprocal
relations for modeling of the non-linear kinetics of fluids. Another exception
is the recent work by Huang et al. (2022), who propose Variational Onsager
Neural Networks to parameterize the free energy potential and the dissipation
potential for thermodynamic processes including phase transformations and
diffusion. The method is also applied to learn the constitutive relations of a
solid material; however, the provided benchmark test is restricted to one-
dimensional viscoelasticity, with one scalar internal variable and a quadratic
dissipation potential. Thus, more complicated material behaviors with several
internal variables and non-quadratic dissipation potentials (such as e.g. in
plasticity) remain untackled. A further stream of research (González et al.,
2019b; Hernández et al., 2021; Zhang et al., 2022) combines data-driven
methods and machine learning with the GENERIC framework, proposed by Grmela
and Öttinger (1997) and related to the theory of generalized standard
materials (Hütter and Svendsen, 2011; Mielke, 2011). Their objective is to
identify two thermodynamic potentials and additionally (unlike with the theory
of generalized standard materials) two system matrices that govern the
evolution of the state variables. The training, however, is supervised and
applications to solid material modeling are thus far restricted to quadratic
thermodynamic potentials, e.g., linear viscoelasticity (González et al.,
2019a), while once again more sophisticated material behavior remains
untackled.

In the following, we provide an introduction to the theory of generalized
standard materials and the resulting constitutive equations (Section 2).
Subsequently, we integrate this theory in the EUCLID framework (Section 3),
which we test on numerical benchmarks in Section 4. We finally draw
conclusions in Section 5.

##  Section 2 Generalized standard materials

In this section, we concisely introduce the fundamental notions on generalized
standard material models, their properties, and the numerical solution of the
related constitutive equations. Many more details on the models can be
retrieved e.g. in the early papers by Ziegler (1957, 1958) and Halphen and
Nguyen (1975) or in the recent presentation by Steinmann and Runesson (2021).

###  2.1 Thermodynamic potentials

We consider a material in three-dimensional space, undergoing infinitesimal
strains under isothermal conditions without phase transitions. The state of
the material is assumed to be uniquely described by a set of $n+1$ state
variables $\\{\bm{\varepsilon},\bm{\alpha}_{1},\dots,\bm{\alpha}_{n}\\}$,
where $\bm{\varepsilon}$ is the (observable) infinitesimal strain tensor and
$\bm{\alpha}_{i}$, $i\in 1,...,n$, are internal (i.e., non-observable) state
variables, typically describing the current state in relation to irreversible
phenomena taking place in the material. Although $\bm{\alpha}_{i}$ may be
tensors of any order, they will be denoted by boldface letters in the
subsequent derivations. For notational compactness we will also denote the set
of internal state variables as
$\bm{\alpha}=\\{\bm{\alpha}_{1},\dots,\bm{\alpha}_{n}\\}$.

By definition, generalized standard materials are materials for which the two
following scalar functions, also called thermodynamic potentials, exist:

  * 1.

The Helmholtz free energy density potential
$\psi\left(\bm{\varepsilon},\bm{\alpha}\right)$ (in the following also simply
denoted as free energy potential), that depends on the state variables and is
assumed to be strictly convex, non-negative
($\psi\left(\bm{\varepsilon},\bm{\alpha}\right)\geq 0$), zero at the origin
($\psi\left(\bf 0,\bf 0\right)=0$) and smooth, i.e., continuously
differentiable.

  * 2.

The dissipation rate density potential
$\pi\left(\dot{\bm{\varepsilon}},\dot{\bm{\alpha}}\right)$ (or dissipation
potential in short), that depends on the time derivatives of the state
variables (denoted by superposed dots), is convex, non-negative
($\pi\left(\dot{\bm{\varepsilon}},\dot{\bm{\alpha}}\right)\geq 0$), zero at
the origin ($\pi\left(\bf 0,\bf 0\right)=0$) and continuous. It may be weakly
convex, i.e., locally linear, and non-smooth, i.e., not continuously
differentiable.

For generalized standard materials, the specific choice of the functional form
of the two thermodynamic potentials, along with the knowledge of the values of
the contained parameters, completely characterizes the material behaviour. The
properties of the free energy potential guarantee material stability and
uniquely defined partial derivatives. Those of the dissipation potential will
be seen to guarantee by construction the fulfillment of the second law of
thermodynamics. Hence, a material model that fits within the framework of
generalized standard materials is by construction stable and thermodynamically
consistent.

Note that, in general, the thermodynamic potentials could also depend on
spatial gradients of the state variables, e.g.,
$\psi\left(\bm{\varepsilon},\nabla\bm{\varepsilon},\bm{\alpha},\nabla\bm{\alpha}\right)$
(see Nguyen (2010); Miehe (2011)), or the dissipation potential could also
depend on the state variables, i.e.,
$\pi\left(\bm{\varepsilon},\bm{\alpha},\dot{\bm{\varepsilon}},\dot{\bm{\alpha}}\right)$
(see Kumar and Lopez-Pamies (2016)). In this initial investigation, we will
consider the ”classical” dependencies and will not be concerned with such
cases. However, they would certainly be of interest for further generalization
of the method proposed in this paper, e.g., to softening material behavior.

###  2.2 State laws, complementary laws and constitutive laws

Under the present assumptions, the so-called intrinsic dissipation (or
mechanical dissipation) ${D}$ reads

${D}=\bm{\sigma}\colon\dot{\bm{\varepsilon}}-\dot{\psi}\geq 0,$ (1)

where $\bm{\sigma}$ is the Cauchy stress tensor and the inequality, also known
as dissipation inequality, is to be seen as a requirement posed to the
material model in order to satisfy the second principle of continuum
thermodynamics. Using the chain rule, the time derivative of the free energy
potential can be expressed as

$\dot{\psi}=\frac{\partial\psi}{\partial\bm{\varepsilon}}\colon\dot{\bm{\varepsilon}}+\frac{\partial\psi}{\partial\bm{\alpha}_{i}}\colon\dot{\bm{\alpha}}_{i},$
(2)

where we have adopted the Einstein convention for summation over repeating
indices. Let us now define the reversible or energetic stress
$\bm{\sigma}^{en}$ and the energetic driving forces
$\bm{\mathcal{A}}^{en}_{i}$ as follows

$\bm{\sigma}^{en}=\frac{\partial\psi}{\partial\bm{\varepsilon}}\left(\bm{\varepsilon},\bm{\alpha}\right),\qquad\bm{\mathcal{A}}^{en}_{i}=\frac{\partial\psi}{\partial\bm{\alpha}_{i}}\left(\bm{\varepsilon},\bm{\alpha}\right).$
(3)

In simple terms, $\bm{\sigma}^{en}$ represents the reversible part of the
stress, while $\bm{\mathcal{A}}^{en}_{i}$ is the energetic driving force
associated to the state variable $\bm{\alpha}_{i}$. Equations ((3)) are known
as state laws and define the thermodynamic forces which are available at a
given state. By substituting (2) and (3) in (1) we obtain

${D}=(\bm{\sigma}-\bm{\sigma}^{en})\colon\dot{\bm{\varepsilon}}-\bm{\mathcal{A}}^{en}_{i}\colon\dot{\bm{\alpha}}_{i}=\bm{\sigma}^{dis}\colon\dot{\bm{\varepsilon}}+\bm{\mathcal{A}}^{dis}_{i}\colon\dot{\bm{\alpha}}_{i}\geq
0.$ (4)

where we have defined the irreversible or dissipative stress
$\bm{\sigma}^{dis}$ and the dissipative driving force
$\bm{\mathcal{A}}^{dis}_{i}$ as follows

$\bm{\sigma}^{dis}=\bm{\sigma}-\bm{\sigma}^{en},\qquad\bm{\mathcal{A}}^{dis}_{i}=-\bm{\mathcal{A}}^{en}_{i}.$
(5)

For the case of non-dissipative material behavior (${D}=0$), we obtain
$\bm{\sigma}^{dis}=\bf 0$ and $\bm{\mathcal{A}}^{dis}_{i}=\bf 0$, implying
that the entire constitutive behavior is described by the state laws. For the
general case ${D}\geq 0$, which is considered here, the state laws are not
sufficient and must be supplemented by the so-called complementary laws

$\bm{\sigma}^{dis}=\frac{\partial\pi}{\partial\dot{\bm{\varepsilon}}}\left(\dot{\bm{\varepsilon}},\dot{\bm{\alpha}}\right),\qquad\bm{\mathcal{A}}^{dis}_{i}=\frac{\partial\pi}{\partial\dot{\bm{\alpha}}_{i}}\left(\dot{\bm{\varepsilon}},\dot{\bm{\alpha}}\right),$
(6)

which describe the relations between the advancement of irreversible phenomena
and the corresponding thermodynamic forces. In the case of a non-smooth
dissipation potential, the partial derivatives above have to be substituted by
subderivatives (Rockafellar, 2015). In the following, for notational
simplicity we denote both partial derivatives and subderivatives with the same
symbol. We further introduce the set of dissipative driving forces as
$\bm{\mathcal{A}}^{dis}=\\{\bm{\mathcal{A}}^{dis}_{1},\dots,\bm{\mathcal{A}}^{dis}_{n}\\}$.
Substituting (6) into (4) yields

${D}=\frac{\partial\pi}{\partial\dot{\bm{\varepsilon}}}\colon\dot{\bm{\varepsilon}}+\frac{\partial\pi}{\partial\dot{\bm{\alpha}}_{i}}\colon\dot{\bm{\alpha}}_{i}.$
(7)

It is thus evident that, whatever the specific choice of $\pi$ with the
postulated properties, it is $D\geq 0$, i.e., the dissipation inequality is
satisfied by construction.

The combination of state laws and complementary laws is sufficient to
completely describe the material behavior, which is governed by the
constitutive laws

$\bm{\sigma}=\frac{\partial\psi}{\partial\bm{\varepsilon}}\left(\bm{\varepsilon},\bm{\alpha}\right)+\frac{\partial\pi}{\partial\dot{\bm{\varepsilon}}}\left(\dot{\bm{\varepsilon}},\dot{\bm{\alpha}}\right),\qquad\bm{0}=\frac{\partial\psi}{\partial\bm{\alpha}_{i}}\left(\bm{\varepsilon},\bm{\alpha}\right)+\frac{\partial\pi}{\partial\dot{\bm{\alpha}}_{i}}\left(\dot{\bm{\varepsilon}},\dot{\bm{\alpha}}\right).$
(8)

###  2.3 Dual dissipation potential and dual complementary laws

An equivalent alternative to the complementary laws is obtained by applying
the Legendre-Fenchel transformation111As the dissipation potential might not
be strictly convex and smooth, the conventional Legendre transformation is not
applicable, see Touchette (2007) for details. to the dissipation potential to
obtain the dual dissipation potential

$\pi^{*}\left(\bm{\sigma}^{dis},\bm{\mathcal{A}}^{dis}\right)=\max_{\dot{\bm{\varepsilon}},\mathbf{\dot{\bm{\alpha}}}}\left\\{\bm{\sigma}^{dis}\colon\dot{\bm{\varepsilon}}+\bm{\mathcal{A}}^{dis}_{i}\colon\dot{\bm{\alpha}}_{i}-\pi\left(\dot{\bm{\varepsilon}},\dot{\bm{\alpha}}\right)\right\\},$
(9)

from which we obtain

$\dot{\bm{\varepsilon}}=\frac{\partial\pi^{*}}{\partial\bm{\sigma}^{dis}}\left(\bm{\sigma}^{dis},\bm{\mathcal{A}}^{dis}\right),\qquad\dot{\bm{\alpha}}_{i}=\frac{\partial\pi^{*}}{\partial\bm{\mathcal{A}}^{dis}_{i}}\left(\bm{\sigma}^{dis},\bm{\mathcal{A}}^{dis}\right),$
(10)

known as dual complementary laws. These laws describe the advancement of
irreversible phenomena as a function of the corresponding thermodynamic
forces, i.e., they are basically evolution equations for the strain and for
the internal variables, which is why they are often more useful than (6). Once
again, should $\pi^{*}$ be non-smooth, the partial derivatives have to be
intended as sub-derivatives. Note that, due to a property of the Legendre-
Fenchel transform, convexity of $\pi$ is equivalent to that of $\pi^{*}$.

###  2.4 Strain-rate-independent dissipation potential

For solid materials, it is often assumed that the dissipation potential does
not depend on the strain rate, i.e., $\pi=\pi\left(\dot{\bm{\alpha}}\right)$
(see, e.g., Lemaitre and Chaboche (1994); Michel and Suquet (2016); Blühdorn
et al. (2022)). As a consequence, it is
$\pi^{*}=\pi^{*}\left(\bm{\mathcal{A}}^{dis}\right)$ and the constitutive laws
can be simplified to

$\bm{\sigma}=\frac{\partial\psi}{\partial\bm{\varepsilon}}\left(\bm{\varepsilon},\bm{\alpha}\right),\qquad\dot{\bm{\alpha}}_{i}=\left.\frac{\partial\pi^{*}}{\partial\bm{\mathcal{A}}^{dis}_{i}}\left(\bm{\mathcal{A}}^{dis}\right)\right\rvert_{\bm{\mathcal{A}}^{dis}_{i}=-\frac{\partial\psi}{\partial{\bm{\alpha}}_{i}}}.$
(11)

In this work, we consider a dissipation potential which is independent of the
strain rate. However, only minor modifications would be needed to extend the
framework developed in the following to the case of a strain-rate-dependent
dissipation potential. Note that the choice of a strain-rate-independent
dissipation potential does not imply a strain-rate-independent material model,
as rate dependency is kept in the evolution of the internal variables.

###  2.5 Numerical solution of the constitutive equations

Given an initial condition for the state variables, Equations (11) can be
solved either in strain or in stress control for the history of the state
variables and of the dependent variables. Here, we will consider strain
control, i.e., the strain history is given and the objective is to compute the
stress and the internal variables at each time. The conventional strategy to
solve such problems is to compute the partial derivatives (or subderivatives)
of the thermodynamic potentials analytically and to solve at each time the
resulting equations, which are non-linear in general, using e.g. the Newton-
Raphson method. If the dissipation potential is non-smooth, algorithms like
the return mapping in elastoplasticity (Simo and Hughes, 1998; Neto et al.,
2008) are commonly applied.

Alternative strategies for solving equations like (11) employ automatic
differentiation, see Korelc (2002); Rothe and Hartmann (2015); Blühdorn et al.
(2022). In this way, derivatives do not need to be evaluated by hand and the
process of implementing the constitutive equation solver for new models is
greatly simplified, which can be particularly beneficial during the
development stage of new material models. However, such automatic solution
strategies also have drawbacks. For non-smooth thermodynamic potentials,
automatic differentiation cannot be applied, thus the potentials must be
approximated by a continuously differentiable regularization (examples are
provided by Steinmann and Runesson (2021)), whose choice is user-dependent and
does not follow rigorous rules. The regularization also introduces additional
parameters and non-linearity to the problem. Further, while conventional
methods often take advantage of problem-specific simplifications that increase
the computational efficiency (see, e.g., the scalar internal variable
substitution during the viscoplastic corrector step in A), automatic methods
are blind to such possibilities, thus they tend to be computationally more
expensive. For these reasons, in this work we rely on the conventional methods
for solving the constitutive equations and for computing the consistent
tangent (see A and B for details).

##  Section 3 Automated discovery of generalized standard material models

In the following, we describe our novel approach to automatically discover
generalized standard (hence, stable and thermodynamically consistent) material
models based on full-field displacement and global force data. The backbone of
the approach, which we denote as EUCLID in line with our previous work
(Flaschel et al., 2021, 2022; Joshi et al., 2022; Thakolkaran et al., 2022),
can be outlined as follows: we start from a very wide material model space,
which we construct by exploiting the flexibility and a priori thermodynamic
consistency of the generalized standard material model framework; we constrain
this model space by enforcing the satisfaction of momentum balance (in weak
sense) on the source data; we exploit sparse regression to ensure parsimonity
of the discovered model. Through these three key ingredients, we end up
simultaneously performing model selection and identification of the unknown
material parameters. Unlike in our previous papers, our model selection here
is not restricted to taking place within a predefined category of material
behavior (such as elasticity or plasticity), but the category itself is one of
the outcomes of the approach.

###  3.1 Material model library

We start by constructing a parametric library of three-dimensional and
isotropic generalized standard material models, out of which the EUCLID
algorithm will later select a model that best fits the data while at the same
time being mathematically concise. In particular, we choose a combination of
generalized Maxwell viscoelastic models and viscoplastic models including
isotropic and kinematic hardening mechanisms, see Figure 1 for a rheological
illustration. Thus, the model library includes a large portion of the models
included in the catalogue by Steinmann and Runesson (2021)222In the catalogue
by Steinmann and Runesson (2021), the models are presented in one-dimensional
form.. For similar viscoelastic-viscoplastic modeling approaches, the reader
is referred to Giunta and Angela Pisano (2006); Miled et al. (2011); Miled
(2011); Sun and Zhu (2013).

Figure 1: Rheological representation of the material model library. For the
interpretation of the isotropic hardening element see Bröcker and Matzenmiller
(2012).

####  3.1.1 Thermodynamic potentials

The adopted model library for the free energy potential $\psi$ reads

$\displaystyle\psi\left(\bm{\varepsilon},\bm{\alpha}\right)$
$\displaystyle=\frac{1}{2}\left(2G\operatorname{dev}(\bm{\varepsilon}-\bm{\alpha}_{{I}})\colon\operatorname{dev}(\bm{\varepsilon}-\bm{\alpha}_{I})+K\operatorname{tr}^{2}(\bm{\varepsilon}-\bm{\alpha}_{{I}})\right)+\frac{1}{2}H^{iso}\alpha_{{II}}^{2}+\frac{1}{2}H^{kin}\bm{\alpha}_{{III}}\colon\bm{\alpha}_{{III}}$
(12)
$\displaystyle+\frac{1}{2}\sum_{j=1}^{n_{\text{MW}}}\left(2G_{j}\operatorname{dev}(\bm{\varepsilon}-\bm{\alpha}_{j}-\bm{\alpha}_{{I}})\colon\operatorname{dev}(\bm{\varepsilon}-\bm{\alpha}_{j}-\bm{\alpha}_{{I}})+K_{j}\operatorname{tr}^{2}(\bm{\varepsilon}-\bm{\alpha}_{j}-\bm{\alpha}_{{I}})\right),$

where
$\bm{\alpha}=\\{\bm{\alpha}_{1},\dots,\bm{\alpha}_{n_{\text{MW}}},\bm{\alpha}_{{I}},\alpha_{{II}},\bm{\alpha}_{{III}}\\}$
are the $n_{\text{MW}}+3$ internal variables. Here, $n_{\text{MW}}$ denotes
the number of viscoelastic Maxwell elements considered in the model library;
$\bm{\alpha}_{1},\dots,\bm{\alpha}_{n_{\text{MW}}}$ have the meaning of
viscous strains in the Maxwell elements; $\bm{\alpha}_{I}$ is the plastic
strain tensor; $\bm{\alpha}_{II}$ and $\bm{\alpha}_{III}$ are the two
hardening variables corresponding to isotropic and kinematic hardening; $G$
and $K$ are the long-term shear and bulk moduli; $G_{i}$ and $K_{i}$ are the
shear and bulk moduli of the $i^{th}$ viscoelastic Maxwell element; and
$H^{iso}$ and $H^{kin}$ are material parameters related to isotropic and
kinematic hardening. At this point, the following considerations are in order:

  * 1.

The actual number of internal variables needed for the accurate description of
the material behavior is not known a priori. Our proposed strategy here
consists in introducing in the library a number of such variables deemed to be
larger than the actual number of internal variables needed to describe the
material state. As will become evident later, unnecessary internal variables
will become automatically inactive through the sparsity regularization.

  * 2.

Compared to our previous investigations focusing on a single material class,
here we are limiting ourselves to a reduced library per class. E.g., unlike in
Flaschel et al. (2021), we are only considering linear elasticity within the
kinematically linear framework, and in contrast to Flaschel et al. (2022), we
fix the shape of the initial yield surface and introduce a less versatile
hardening description. The reason is that here the focus lies on the ability
of EUCLID to discriminate among different material classes, rather than on its
versatility within a given material class (which was extensively proved in our
previous papers). Obviously, further extensions are possible and may well be
addressed in future developments.

The model library for the dual dissipation potential $\pi^{*}$ reads

$\pi^{*}\left(\bm{\mathcal{A}}^{dis}\right)=\pi^{*}_{\text{VE}}\left(\bm{\mathcal{A}}^{dis}_{1},\dots,\bm{\mathcal{A}}^{dis}_{n_{\text{MW}}}\right)+\pi^{*}_{\text{VP}}\left(\bm{\mathcal{A}}^{dis}_{{I}},\mathcal{A}^{dis}_{{II}},\bm{\mathcal{A}}^{dis}_{{III}}\right),$
(13)

with a viscoelastic contribution

$\pi^{*}_{\text{VE}}\left(\bm{\mathcal{A}}^{dis}_{1},\dots,\bm{\mathcal{A}}^{dis}_{n_{\text{MW}}}\right)=\frac{1}{2}\sum_{j=1}^{n_{\text{MW}}}\left(\frac{1}{2G_{j}g_{j}}\operatorname{dev}(\bm{\mathcal{A}}^{dis}_{j})\colon\operatorname{dev}(\bm{\mathcal{A}}^{dis}_{j})+\frac{1}{9K_{j}k_{j}}\operatorname{tr}^{2}(\bm{\mathcal{A}}^{dis}_{j})\right),$
(14)

and a viscoplastic contribution

$\pi^{*}_{\text{VP}}\left(\bm{\mathcal{A}}^{dis}_{{I}},\mathcal{A}^{dis}_{{II}},\bm{\mathcal{A}}^{dis}_{{III}}\right)=\begin{cases}0&\text{if}~{}\sqrt{\frac{3}{2}}\left\|\operatorname{dev}\left(\bm{\mathcal{A}}^{dis}_{{I}}+\bm{\mathcal{A}}^{dis}_{{III}}\right)\right\|\leq\sigma_{y}-\mathcal{A}^{dis}_{{II}}\\\
\frac{1}{2}\frac{1}{\eta_{p}}\left(\sqrt{\frac{3}{2}}\left\|\operatorname{dev}\left(\bm{\mathcal{A}}^{dis}_{{I}}+\bm{\mathcal{A}}^{dis}_{{III}}\right)\right\|-\sigma_{y}+\mathcal{A}^{dis}_{{II}}\right)^{2}&\text{otherwise},\end{cases}$
(15)

where
$\bm{\mathcal{A}}^{dis}=\\{\bm{\mathcal{A}}^{dis}_{1},\dots,\bm{\mathcal{A}}^{dis}_{n_{\text{MW}}},\bm{\mathcal{A}}^{dis}_{{I}},\mathcal{A}^{dis}_{{II}},\bm{\mathcal{A}}^{dis}_{{III}}\\}$
are the dissipative driving forces associated to the internal variables with
the same subscripts, $g_{i}$ and $k_{i}$ are the shear and bulk relaxation
times of the $i^{th}$ viscoelastic Maxwell element, $\sigma_{0}$ is the yield
stress, and $\eta_{p}$ is the viscosity associated to the viscoplastic
behavior.

####  3.1.2 Constitutive equations

The constitutive equations for the chosen model library are obtained by
substituting the thermodynamic potentials of Section 3.1.1 in (11). The stress
is derived as

$\displaystyle\bm{\sigma}=\frac{\partial\psi}{\partial\bm{\varepsilon}}$
$\displaystyle=2G\operatorname{dev}(\bm{\varepsilon}-\bm{\alpha}_{{I}})+K\operatorname{tr}(\bm{\varepsilon}-\bm{\alpha}_{{I}})\mathbf{I}$
(16)
$\displaystyle+\sum_{j=1}^{n_{\text{MW}}}\left(2G_{j}\operatorname{dev}(\bm{\varepsilon}-\bm{\alpha}_{j}-\bm{\alpha}_{{I}})+K_{j}\operatorname{tr}(\bm{\varepsilon}-\bm{\alpha}_{j}-\bm{\alpha}_{{I}})\mathbf{I}\right).$

Differentiation of the dual dissipation potential with respect to the
viscoelastic driving forces leads to the viscoelastic evolution laws

$\dot{\bm{\alpha}}_{j}=\frac{\partial\pi^{*}}{\partial\bm{\mathcal{A}}^{dis}_{j}}=\frac{1}{2G_{j}g_{j}}\operatorname{dev}(\bm{\mathcal{A}}^{dis}_{j})+\frac{1}{9K_{j}k_{j}}\operatorname{tr}(\bm{\mathcal{A}}^{dis}_{j})\mathbf{I}.$
(17)

These laws need to be evaluated at

$\bm{\mathcal{A}}^{dis}_{j}=-\frac{\partial\psi}{\partial\bm{\alpha}_{j}}=\left(2G_{j}\operatorname{dev}(\bm{\varepsilon}-\bm{\alpha}_{j}-\bm{\alpha}_{{I}})+K_{j}\operatorname{tr}(\bm{\varepsilon}-\bm{\alpha}_{j}-\bm{\alpha}_{{I}})\mathbf{I}\right),$
(18)

which leads to

$\dot{\bm{\alpha}}_{j}=\frac{1}{g_{j}}\operatorname{dev}(\bm{\varepsilon}-\bm{\alpha}_{j}-\bm{\alpha}_{{I}})+\frac{1}{3k_{j}}\operatorname{tr}(\bm{\varepsilon}-\bm{\alpha}_{j}-\bm{\alpha}_{{I}}).$
(19)

Differentiation of the dual dissipation potential with respect to the
viscoplastic driving forces leads to the viscoplastic evolution laws

$\displaystyle\dot{\bm{\alpha}}_{{I}}=\dot{\bm{\alpha}}_{{III}}=\frac{\partial\pi^{*}}{\partial\bm{\mathcal{A}}^{dis}_{{I}}}$
$\displaystyle=\begin{cases}0&\text{if}~{}f\leq 0\\\
\frac{1}{\eta_{p}}f\sqrt{\frac{3}{2}}\frac{\operatorname{dev}\left(\bm{\mathcal{A}}^{dis}_{{I}}+\bm{\mathcal{A}}^{dis}_{{III}}\right)}{\left\|\operatorname{dev}\left(\bm{\mathcal{A}}^{dis}_{{I}}+\bm{\mathcal{A}}^{dis}_{{III}}\right)\right\|}&\text{otherwise},\end{cases}$
(20)
$\displaystyle\dot{{\alpha}}_{{II}}=\frac{\partial\pi^{*}}{\partial{\mathcal{A}}^{dis}_{{II}}}$
$\displaystyle=\begin{cases}0&\text{if}~{}f\leq 0\\\
\frac{1}{\eta_{p}}f&\text{otherwise},\end{cases}$

where we defined

$f=\sqrt{\frac{3}{2}}\left\|\operatorname{dev}\left(\bm{\mathcal{A}}^{dis}_{{I}}+\bm{\mathcal{A}}^{dis}_{{III}}\right)\right\|-\sigma_{y}+\mathcal{A}^{dis}_{{II}}.$
(21)

We notice that $\dot{\bm{\alpha}}_{{I}}=\dot{\bm{\alpha}}_{{III}}$. These
evolution laws need to be evaluated at
$\bm{\mathcal{A}}^{dis}_{{I}}=-\frac{\partial\psi}{\partial\bm{\alpha}_{{I}}}=\bm{\sigma}$,
${\mathcal{A}}^{dis}_{{II}}=-\frac{\partial\psi}{\partial{\alpha}_{{II}}}=-H^{iso}\alpha_{{II}}$
and
$\bm{\mathcal{A}}^{dis}_{{III}}=-\frac{\partial\psi}{\partial\bm{\alpha}_{{III}}}=-H^{kin}\bm{\alpha}_{{III}}$,
leading to

$\displaystyle\dot{\bm{\alpha}}_{{I}}=\dot{\bm{\alpha}}_{{III}}$
$\displaystyle=\begin{cases}0&\text{if}~{}f\leq 0\\\
\frac{1}{\eta_{p}}f\sqrt{\frac{3}{2}}\frac{\operatorname{dev}\left(\bm{\sigma}-H^{kin}\bm{\alpha}_{{III}}\right)}{\left\|\operatorname{dev}\left(\bm{\sigma}-H^{kin}\bm{\alpha}_{{III}}\right)\right\|}&\text{otherwise},\end{cases}$
(22) $\displaystyle\dot{{\alpha}}_{{II}}$
$\displaystyle=\begin{cases}0&\text{if}~{}f\leq 0\\\
\frac{1}{\eta_{p}}f&\text{otherwise},\end{cases}$

with

$f=\sqrt{\frac{3}{2}}\left\|\operatorname{dev}\left(\bm{\sigma}-H^{kin}\bm{\alpha}_{{III}}\right)\right\|-\sigma_{y}-H^{iso}\alpha_{{II}}.$
(23)

For solving the incremental constitutive equations at a material point, we
choose an implicit Euler discretization in time. Given the state variables at
the previous time step, the strain at the current time step and the time step
size, a viscoelastic predictor / viscoplastic corrector algorithm (see A) is
applied to compute the stress and the internal variables at the current time
step. The consistent tangent modulus is given in B.

####  3.1.3 Material parameters

All tunable material parameters in the material model library are collected in
the material parameter vector

$\bm{\theta}=\left[G,K,G_{1},\dots,G_{n_{\text{MW}}},\frac{1}{g_{1}},\dots,\frac{1}{g_{n_{\text{MW}}}},K_{1},\dots,K_{n_{\text{MW}}},\frac{1}{k_{1}},\dots,\frac{1}{k_{n_{\text{MW}}}},\frac{1}{\sigma_{0}},\eta_{p},H^{iso},H^{kin}\right]^{T}.$
(24)

To ensure the convexity of the thermodynamic potentials, all parameters in
$\bm{\theta}$ are assumed to be non-negative.

Note that some parameters appear in the material parameter vector through
their reciprocals. This facilitates the simplification of the model in special
cases for which these parameters tend to infinity (see also Section 3.4). For
example, if the yield stress tends to infinity, the model can be considered as
purely viscoelastic, and if a relaxation time tends to infinity, the viscous
effect of the corresponding Maxwell element can be ignored. In practice, it is
more convenient to work with vanishing quantities than with diverging ones.
For the numerical treatment, the parameters in reciprocal form are assumed to
be strictly positive.

###  3.2 Available data

By conducting a mechanical test on a specimen made of the material under
consideration, two-dimensional full-field displacements and global reaction
force data can be acquired to serve as input for EUCLID. The point-wise
displacement data are given in the form
$\\{\mathbf{u}^{a,t}:a=1,\dots,n_{n};t=1,\dots,n_{t}\\}$, where $n_{n}$ is the
number of points in space and $n_{t}$ is the number of time steps at which
displacements are known. The reaction force data are given in the form
$\\{\hat{R}^{\beta,t}:\beta=1,\dots,n_{\beta};t=1,\dots,n_{t}\\}$, where
$n_{\beta}$ is the number of measured reaction force components.

To interpolate the given displacement data, we create a finite element mesh
with shape functions $\\{N^{a}(\mathbf{X}):a=1,\dots,n_{n}\\}$, such that each
node of the mesh corresponds to one of the points where the displacement is
known. The displacement field at each time is then given as

$\mathbf{u}^{t}(\mathbf{X})=\sum_{a=1}^{n_{n}}N^{a}(\mathbf{X})\mathbf{u}^{a,t},\quad\forall\
t=1,\dots,n_{t},$ (25)

which is differentiated in space to obtain the infinitesimal strain field at
each time step

$\bm{\varepsilon}^{t}(\mathbf{X})=\nabla^{\text{sym}}\mathbf{u}^{t}(\mathbf{X})=\sum_{a=1}^{n_{n}}\frac{1}{2}\left[\mathbf{u}^{a,t}\otimes\nabla
N^{a}(\mathbf{X})+\nabla
N^{a}(\mathbf{X})\otimes\mathbf{u}^{a,t}\right],\quad\forall\
t=1,\dots,n_{t}.$ (26)

###  3.3 Cost function based on the weak linear momentum balance

Given full-field kinematic data and net reaction forces, EUCLID employs the
balance of linear momentum as a physical constraint to determine the optimal
values of the unknown parameters in the model library. In our previous
investigations, we observed this constraint to be effective in contrasting the
strong ill-posedness induced by the lack of labeled stress-strain data pairs.
In the following, we introduce a cost function based on the sum of squared
unbalanced forces, thus establishing the weak linear momentum balance in the
interior and at the boundary of the test specimen (see Flaschel et al. (2021,
2022) for more information). Later, this cost function will be minimized along
with a sparsity promoting penalty term to select a parsimonious model while
determining its unknown parameters.

The weak formulation of the linear momentum balance under the assumption of
vanishing inertia and body forces at a generic time step $t$ reads

$\int_{\Omega}\bm{\sigma}^{t}(\Delta
t,\bm{\varepsilon}^{t},\bm{\alpha}^{t-1},\bm{\theta})\colon\nabla\mathbf{v}\;\\!\mathrm{d}A-\int_{\partial\Omega}\hat{\mathbf{t}}^{t}\cdot\mathbf{v}\;\\!\mathrm{d}S=0,\quad\forall\
\text{admissible}\ \mathbf{v},\quad\forall\ t=1,\dots,n_{t},$ (27)

where $\Omega$ and $\partial\Omega$ denote the specimen domain and boundary,
respectively, and $\hat{\mathbf{t}}$ is the surface traction on
$\partial\Omega$. The stress at the considered time step $\bm{\sigma}^{t}$
depends on the time step size $\Delta t$, the current strain
$\bm{\varepsilon}^{t}$, the internal variables at the previous time step
$\bm{\alpha}^{t-1}$ and the material parameters $\bm{\theta}$ through the
stress update algorithm (see A). The weak form must be satisfied for any
admissible (i.e., sufficiently smooth) test functions $\mathbf{v}$; we
discretize them with the same shape functions adopted for the displacement
data interpolation

$\mathbf{v}(\mathbf{X})=\sum_{a=1}^{n_{n}}N^{a}(\mathbf{X})\,\mathbf{v}^{a},$
(28)

and we assume them to be constant in time. This leads to the discretized weak
linear momentum balance

$\sum_{a=1}^{n_{n}}\mathbf{v}^{a}\cdot\left[\underbrace{\int_{\Omega}\bm{\sigma}^{t}(\Delta
t,\bm{\varepsilon}^{t},\bm{\alpha}^{t-1},\bm{\theta})\nabla
N^{a}(\mathbf{X})\;\\!\mathrm{d}A}_{\mathbf{F}^{a,t}(\bm{\theta})}-\int_{\partial\Omega}{\hat{\mathbf{t}}}^{t}N^{a}(\mathbf{X})\;\\!\mathrm{d}S\right]=0,\quad\forall\
\text{admissible}\ \mathbf{v},\quad\forall\ t=1,\dots,n_{t},$ (29)

where $\mathbf{F}^{a,t}$ denotes the internal force at node $a$ of the finite
element mesh for time step $t$. Since the internal forces depend on the stress
at the current time step, which in turn depends on the internal variables at
the previous time step, they can only be computed (for a given set of material
parameters $\bm{\theta}$) by iterating through all load steps and solving the
constitutive equations at each step. To this end, the internal variables at
the initial time step are assumed to be zero.

We denote the set of all nodal degrees of freedom in the finite element mesh
as $\mathcal{D}=\\{(a,i):a=1,\dots,n_{n};i=1,2\\}$. We further partition this
set into two subsets: $\mathcal{D}^{\text{free}}$, containing the free degrees
of freedom, i.e., those that do not correspond to any of the constrained
boundaries, and $\mathcal{D}^{\text{disp}}$, containing the displacement-
constrained degrees of freedom, such that
$\mathcal{D}^{\text{free}}\cup\mathcal{D}^{\text{disp}}=\mathcal{D}$ and
$\mathcal{D}^{\text{free}}\cap\mathcal{D}^{\text{disp}}=\emptyset$. Further,
we denote the subset of $\mathcal{D}^{\text{disp}}$ corresponding to the
reaction force $\hat{R}^{\beta,t}$ as $\mathcal{D}^{\text{disp},\beta}$.

At the free degrees of freedom, the traction force $\hat{\mathbf{t}}$ vanishes
and we obtain from (29)

$F^{a,t}_{i}(\bm{\theta})=0,\quad\forall\
(a,i)\in\mathcal{D}^{\text{free}},\quad\forall\ t=1,\dots,n_{t},$ (30)

whereas at the displacement-constrained boundaries the sum of the nodal
internal forces must equal the global external reaction force

$\sum_{(a,i)\in\mathcal{D}^{\text{disp},\beta}}F^{a,t}_{i}(\bm{\theta})=\sum_{(a,i)\in\mathcal{D}^{\text{disp},\beta}}\int_{\partial\Omega}{\hat{t}}^{t}_{i}N^{a}(\mathbf{X})\;\\!\mathrm{d}S=\hat{R}^{\beta,t},\quad\forall\
\beta=1,\dots,n_{\beta},\quad\forall\ t=1,\dots,n_{t}.$ (31)

Global measures of the violation of linear momentum balance in the interior
and at the boundary can thus be obtained by the following sums of squared
residuals

$C^{\text{free}}(\bm{\theta})=\sum_{t=1}^{n_{t}}\sum_{(a,i)\in\mathcal{D}^{\text{free}}}\left|F^{a,t}_{i}(\bm{\theta})\right|^{2},$
(32)

and

$C^{\text{disp}}(\bm{\theta})=\sum_{t=1}^{n_{t}}\sum_{\beta=1}^{n_{\beta}}\left|\hat{R}^{\beta,t}-\sum_{(a,i)\in\mathcal{D}^{\text{disp},\beta}}F^{a,t}_{i}(\bm{\theta})\right|^{2},$
(33)

which are finally combined in the cost function

$C(\bm{\theta})=C^{\text{free}}(\bm{\theta})+\lambda_{r}C^{\text{disp}}(\bm{\theta}).$
(34)

Here $\lambda_{r}$ is a weighting parameter that scales the interior and
boundary contributions to the cost function. As there are less reaction force
measurements (see (31)) than free degrees of freedom (see (30)), the weighting
parameter should be chosen sufficiently larger than one ($\lambda_{r}>>1$).
Following previous works (Flaschel et al., 2021, 2022), we choose
$\lambda_{r}=100$ and keep it constant throughout this work. Based on our
experience, the choice of $\lambda_{r}$ is not crucial for the success of the
method. Finally, note that in the quantification of the cost we are summing
upon all time steps.

###  3.4 Sparsity promoting regularization

Material models with a large number of material parameters are known to be
less interpretable and poor at extrapolation to unseen strains (Hartmann,
2001). For this reason, our objective with EUCLID is to obtain a parsimonious
model, i.e., a model with as few features as needed to accurately interpret
the available data. Minimizing the cost function in (34) leads in general to a
dense parameter vector $\bm{\theta}$, i.e., a vector with many non-zero
entries. This, in turn, implies that the resulting material model contains
nearly all terms of the thermodynamic potentials contained in the initially
chosen library (see Section 3.1). Thus, nearly all internal variables are
active and the calibrated model is as complex as the entire model space
adopted at the outset, which goes against the objective of discovering an
interpretable model containing only a few terms. In order to obtain a sparse
parameter vector with many zero entries and hence simple expressions for the
thermodynamic potentials, we apply sparse regression, i.e., we minimize (34)
supplemented with a sparsity promoting regularization term

$\bm{\theta}^{\text{opt}}=\operatornamewithlimits{arg\
min}_{\bm{\theta}\geq\bm{\theta}^{\text{min}}}\left(C(\bm{\theta})+\lambda_{p}\|\bm{\theta}\|_{p}^{p}\right),\qquad\text{where}\qquad\|\bm{\theta}\|_{p}=\left(\sum_{i\geq
3}|\theta_{i}|^{p}\right)^{1/p}.$ (35)

The sparsity promoting regularization term was originally proposed by Frank
and Friedman (1993) and Tibshirani (1996) and later applied to problems in
dynamics by Brunton et al. (2016) and to problems in continuum mechanics by
Flaschel et al. (2021, 2022). This term takes high values for dense solution
vectors and small values for sparse solution vectors, whereby the weighting
parameter $\lambda_{p}$ influences its impact on the minimization problem. We
will show later how the value of $\lambda_{p}$ can be automatically selected
to achieve a user-defined compromise between model accuracy and sparsity (see
Section 3.5). The parameter $p$ influences the curvature of the regularization
term. For $p\rightarrow 0$, this term converges to the $L_{0}$-(pseudo)-norm
of $\bm{\theta}$, which counts the number of non-zero entries in the vector.
However, this choice turns the optimization problem into a combinatorial
subset selection problem which becomes computationally intractable for a large
number of features. For this reason, $L_{0}$ regularization is often relaxed
to the convex $L_{1}$ regularization ($p=1$), also known as LASSO (least
absolute shrinkage and selection operator) (Tibshirani, 1996). As the $L_{0}$
regularization, the LASSO promotes sparse solution vectors with many zero
entries, while shrinking the remaining non-zero parameters. Within $[0,1]$, a
smaller value of $p$ reduces the shrinkage of the non-zero parameters at the
cost of computational complexity, as the degree of non-convexity in the
penalty term increases with decreasing $p$. In this work, we thus choose $p=1$
to limit the computational complexity, and keep it constant throughout the
numerical examples. Based on our experience, the choice of $0<p\leq 1$ is not
crucial for the success of the method (Flaschel et al., 2021, 2022).

As different material parameters have different units, the summation in the
regularization term is unphysical and only fulfills a numerical purpose. While
a non-dimensionalization strategy would be recommendable (and is planned for
future studies), here we adopt the simplified approach to compute the
regularization term based on the numerical values of the parameters for fixed
units (chosen as s for time, mm for length and kN for force) and treat the
result as dimensionless, such that the unit of the weighting factor
$\lambda_{p}$ is assumed to be $\text{kN}^{2}$. With this specific choice of
units, the numerical values for the material parameters studied in this work
appear within a reasonable range of orders of magnitude, such that no scaling
issues were observed during numerical testing. Note that we exclude
$\theta_{1}$ and $\theta_{2}$ from the sparsity promoting regularization as we
assume that the elasticity constants $G$, $K$ are always non-zero for a solid
material. The parameter constraints $\bm{\theta}>\bm{\theta}^{\text{min}}$ in
(35) enforce that all parameters are larger than zero, and all parameters that
appear in the minimization problem with their reciprocals (i.e., $1/g_{i}$,
$1/k_{i}$, $1/\sigma_{0}$) are larger than a lower bound, here chosen as
$\theta^{\text{lb}}=10^{-6}$ (see also Section 3.1.3).

###  3.5 Optimization strategy and hyperparameter selection

In the following, we describe the strategies adopted for the solution of the
optimization problem in (35) and for the selection of the hyperparameter
$\lambda_{p}$ (for a summary of all other (hyper-)parameters see C). For the
optimization, we rely on the trust region reflective Newton algorithm
implemented in the Matlab built-in function lsqnonlin, which has proven to be
a suitable choice in terms of convergence and efficiency in previous work
(Flaschel et al., 2022).

To obtain an initial estimate of the solution, the problem is first solved
without regularization, i.e., $\lambda_{p}=0$, and for a reduced set of
parameters, i.e., the linear elastic constants $G$ and $K$, while constraining
the other parameters to be equal to the corresponding values in
$\bm{\theta}^{\text{min}}$. This preconditioning step is computationally
efficient due to the small number of trainable parameters, and at the same
time it provides a good estimate about the order of magnitude of the elastic
parameters. Subsequently, the optimization problem is solved (still with
$\lambda_{p}=0$) for the entire set of unknown parameters. In general, the
function to be minimized in (35) is non-convex and the minimization problem
admits multiple local minima. In order to increase the chance of finding (or
approaching as closely as possible) the global optimum, we run the
optimization in parallel for multiple randomly generated initial guesses
($n_{g}=24$ in this work) and, from all generated solutions, we select the one
corresponding to the lowest cost value. This solution serves as initial guess
for the regularized optimization problem.

The value of the hyperparameter $\lambda_{p}$ is of paramount importance for
the achievement of a good compromise between accuracy and sparsity of the
discovered material model. As follows, we propose an automated algorithm for
its selection. We solve (35) for multiple different values of $\lambda_{p}$
which are successively increased by a factor of two, starting from an initial
small value. In this work, we adopt $\lambda_{p}\in\\{10^{-4}\cdot 2^{j}\
\text{kN}^{2}:j=0,\dots,23\\}$. For every value of $\lambda_{p}$, (35)
delivers a solution for the material parameters and an associated value of the
cost function, which we store in the set

$\mathcal{S}=\\{(\bm{\theta}^{\text{opt}}_{j},C_{j}):j=1,\dots,n_{\lambda}\\},$
(36)

with $n_{\lambda}=24$ as the number of considered $\lambda_{p}$ values. Figure
2 shows an exemplary plot of the cost function and regularization term values
for the different choices of $\lambda_{p}$, which follows a consistently
encountered trend. As $\lambda_{p}$ increases, the cost function value
increases and the regularization term value decreases, implying that the
discovered model becomes decreasingly accurate and increasingly simple. It is
thus evident that this plot carries the information needed to strike a balance
between the conflicting objectives of accuracy and parsimony. Importantly, the
initial range of $\lambda_{p}$ values is characterized by a low rate of
variation of both cost and regularization term, so that at least visually it
is not difficult to identify a threshold value of $\lambda_{p}$ beyond which
the cost start increasing (and the regularization term decreasing) at a much
faster rate. Rather than defining a threshold on $\lambda_{p}$, we prefer to
define one on the cost, $C^{\text{th}}$, with the objective to limit the
tolerated error on the fitting accuracy of the discovered models. In this
work, we choose $C^{\text{th}}$ to be slightly larger than the smallest cost
$C^{\text{min}}$ in the solution set $\mathcal{S}$ but not smaller than
$10^{-5}\ \text{kN}^{2}$, i.e., $C^{\text{th}}=\max\\{10^{-5}\
\text{kN}^{2},1.1\ C^{\text{min}}\\}$. Once this threshold is chosen, all
solutions with $C_{j}\geq C^{\text{th}}$ are considered irrelevant. The
remaining relevant solutions are

$\mathcal{S}^{\text{th}}=\\{(\bm{\theta}^{\text{opt}}_{j},C_{j})\in\mathcal{S}:C_{j}<C^{\text{th}}\\},$
(37)

from which we select the sparsest solution, i.e., the solution in
$\mathcal{S}^{\text{th}}$ with the smallest value of
$\|\bm{\theta}^{\text{opt}}_{j}\|_{p}^{p}$. In all cases, this automated
selection algorithm leads to a solution with both low cost and low
regularization term (see Figure 2), thus striking a user-defined compromise
between model accuracy and sparsity.

Figure 2: Pareto analysis for the automated selection of the hyperparameter
$\lambda_{p}$. The left figure shows the cost function values (in
$\text{kN}^{2}$) and the regularization term values as functions of the
different choices of $\lambda_{p}$ (in $\text{kN}^{2}$). A magnification in
the region around the automatically selected hyperparameter is shown in the
right figure. The threshold for the cost function value and the solution
obtained from the automated selection algorithm are highlighted with blue and
black dashed lines, respectively. The solutions were obtained from the data
corresponding to model LEVP with noise level $\sigma=0.5\mu\text{m}$ (see
Section 4)

###  3.6 Material classification

In a final postprocessing step, after having selected a suitable value for
$\lambda_{p}$, we set all components of $\bm{\theta}$ that lie below a
threshold $\theta^{\text{th}}$ (chosen as $\theta^{\text{th}}=10^{-4}$ in this
work) to their minimal value, i.e., to the corresponding value in
$\bm{\theta}^{\text{min}}$. Further, if the product of the modulus and the
reciprocal of a relaxation time of a Maxwell element lies below the threshold,
both the modulus and the reciprocal of the relaxation time are set to their
minimal values. In this way, unnecessary terms in the thermodynamic potentials
are removed and the discovered material can be classified into different
material classes. If $G$ or $K$ are larger than the threshold, the material is
classified as elastic. If one of the reciprocal relaxation times $1/g_{i}$ or
$1/k_{i}$ is larger than the threshold, the material is classified as
viscoelastic. If, instead, any of the reciprocal relaxation times $1/g_{i}$ or
$1/k_{i}$ is below the threshold, we assume $g_{i}\rightarrow\infty$ or
$k_{i}\rightarrow\infty$, respectively. If the reciprocal of the yield stress
is below the threshold, neither plasticity nor viscoplasticity is active, and
we assume $\sigma_{0}\rightarrow\infty$. If, however, the reciprocal of the
yield stress is larger than the threshold, the material is classified as
plastic and the value of $\eta_{p}$ decides whether the material is further
classified as viscoplastic. Finally, isotropic or kinematic hardening are
active if respectively $H^{iso}$ or $H^{kin}$ are larger than the threshold.

##  Section 4 Numerical benchmarks

To assess the performance of the developed method, we test it on synthetic
data generated by finite element simulations. While testing on real
experimental data will be an important future objective, the deployment of
computational data has the obvious advantage that the ground truth, i.e., the
true material model underlying the data, is exactly known. We select five
different benchmark material models:

  * 1.

E: Elastic

  * 2.

VE: Viscoelastic

  * 3.

VEEP: Viscoelastic (only shear) & elastoplastic with isotropic hardening

  * 4.

EVP: Elastic & viscoplastic with kinematic hardening

  * 5.

VEVP: Viscoelastic & viscoplastic with mixed hardening

Table 1 illustrates more in detail the features active within each benchmark
material model, while the corresponding material parameters are given in Table
2.

###  4.1 Data generation

Figure 3: Displacement-controlled loading conditions of a square domain with
two elliptic holes. The dimensions of the specimen and the loading parameter
$\delta$ are in mm and the loading rate $\dot{\delta}$ is in
$\text{mm}/\text{s}$.

For each benchmark material model, a two-dimensional finite element simulation
under plane-strain conditions is conducted. To obtain a heterogeneous strain
field and hence a rich data set during the experiment, complex geometries are
preferred over simple geometries. We therefore adopt a square specimen with
two elliptic holes as depicted in Figure 3 (left), which has been already
shown to produce a diverse strain field upon deformation (see Flaschel et al.
(2022)), and we load it under displacement-controlled tension and compression.
The loading parameter $\delta$ determines the displacement at the upper
boundary of the specimen. As the model library considered in this work
includes a few rate-dependent models, different rates of the loading parameter
$\dot{\delta}$ are included during data generation as shown in Figure 3
(right). In particular, we consider a number of $n_{k}=20$ loading phases,
indexed by ${k}=1,\dots,n_{k}$, out of which $n_{k}/2$ under tension and
$n_{k}/2$ under compression. Each loading phase consists of a loading period,
during which $\dot{\delta}$ is constant, and a subsequent relaxation period,
during which no further deformation is applied at the boundary
($\dot{\delta}=0$). During each loading period, the same absolute value of
deformation $|\delta|=0.5\text{mm}$ is applied over different loading time
intervals $t^{\text{load}}_{k}$, which are chosen evenly spaced on a
logarithmic scale between $t^{\text{load}}_{1}=10^{-2}\text{s}$ and
$t^{\text{load}}_{20}=10^{2}\text{s}$, i.e.,
$t^{\text{load}}_{k}\in\\{10^{(-2+\frac{4}{n_{k}-1}i)}\text{s}:i=0,\dots,(n_{k}-1)\\}$.
The time intervals of the relaxation periods are assumed all equal to
$t^{\text{relax}}_{k}=10^{2}\text{s}$. The entire problem is discretized in
time with a total number of $n_{t}=800$ time steps, choosing the same number
of time steps for each loading and relaxation period.

From the finite element simulations, the nodal displacements and net reaction
force at the upper specimen boundary are extracted. To mimic real experiments,
we add to the data artificial Gaussian noise as follows

$u^{a,t}_{i}=u^{\text{fem},a,t}_{i}+u^{\text{noise},a,t}_{i},\quad
u^{\text{noise},a,t}_{i}\sim\mathcal{N}(0,\sigma)\quad\forall\quad
a\in\\{1,\dots,n_{n}\\},~{}i\in\\{1,2\\},~{}t\in\\{1,\dots,n_{t}\\},$ (38)

where $u^{\text{fem},a,t}_{i}$ is the $i^{\text{th}}$ displacement component
of node $a$ at the $t^{\text{th}}$ time step obtained from the finite element
simulations, and $u^{\text{noise},a,t}_{i}$ denotes the added noise, sampled
from a Gaussian distribution with zero mean and standard deviation $\sigma$.
Considering modern DIC setups, a reasonable upper limit for the noise is
$\sigma=0.1~{}\mu\text{m}$ (Pierron et al., 2010; Marek et al., 2019). For our
benchmarks, we test different noise levels
$\sigma\in\\{0\mu\text{m},0.1\mu\text{m},0.3\mu\text{m},0.5\mu\text{m}\\}$.
After applying artificial noise to the data, temporal Savitzky–Golay denoising
(Savitzky and Golay, 1964) with quadratic polynomial order and a moving-window
length of $10$ time steps is applied to each loading period and relaxation
period.

###  4.2 Discovered material models

EUCLID is applied to the full-field displacement and reaction force data
obtained from the five different benchmark models with different noise levels
to discover the hidden material behavior. Table 1 illustrates the discovered
material classes in comparison to the true ones, revealing that EUCLID is
indeed able to automatically discriminate among different classes of material
behavior. In other words, purely based on displacement and force data on one
single test, EUCLID is able to discern whether the hidden material behavior is
elastic, viscoelastic, plastic or viscoplastic, and whether isotropic or
kinematic hardening takes place. Note that the automatic selection of the few
terms in the thermodynamic potentials needed to describe the material behavior
implies the automatic selection of the correspondingly active internal
variables. Only for one benchmark, namely, the VEEP model with noise level
$\sigma=0.3\mu\text{m}$, EUCLID fails to predict the true material class as it
discovers a Maxwell element for the bulk material response, which is not
contained in the true material model. The discovered material model is
nevertheless parsimonious, as the effects of rate-dependent plasticity and
kinematic hardening are correctly excluded from the material response.

Table 1: True and discovered material types for different noise levels $\sigma$ (in mm). Benchmarks |  Elasticity |  Viscoelasticity |  # Maxwell (shear) |  # Maxwell (bulk) |  Plasticity |  Viscoplasticity |  Isotropic Hardening |  Kinematic Hardening
---|---|---|---|---|---|---|---|---
E | Truth | on | off | 0 | 0 | off | off | off | off
| $\sigma=0$ | on | off | 0 | 0 | off | off | off | off
| $\sigma=10^{-4}$ | on | off | 0 | 0 | off | off | off | off
| $\sigma=3\cdot 10^{-4}$ | on | off | 0 | 0 | off | off | off | off
| $\sigma=5\cdot 10^{-4}$ | on | off | 0 | 0 | off | off | off | off
VE | Truth | on | on | 1 | 1 | off | off | off | off
| $\sigma=0$ | on | on | 1 | 1 | off | off | off | off
| $\sigma=10^{-4}$ | on | on | 1 | 1 | off | off | off | off
| $\sigma=3\cdot 10^{-4}$ | on | on | 1 | 1 | off | off | off | off
| $\sigma=5\cdot 10^{-4}$ | on | on | 1 | 1 | off | off | off | off
VEEP | Truth | on | on | 1 | 0 | on | off | on | off
| $\sigma=0$ | on | on | 1 | 0 | on | off | on | off
| $\sigma=10^{-4}$ | on | on | 1 | 0 | on | off | on | off
| $\sigma=3\cdot 10^{-4}$ | on | on | 1 | 1 | on | off | on | off
| $\sigma=5\cdot 10^{-4}$ | on | on | 1 | 0 | on | off | on | off
EVP | Truth | on | off | 0 | 0 | on | on | off | on
| $\sigma=0$ | on | off | 0 | 0 | on | on | off | on
| $\sigma=10^{-4}$ | on | off | 0 | 0 | on | on | off | on
| $\sigma=3\cdot 10^{-4}$ | on | off | 0 | 0 | on | on | off | on
| $\sigma=5\cdot 10^{-4}$ | on | off | 0 | 0 | on | on | off | on
VEVP | Truth | on | on | 1 | 1 | on | on | on | on
| $\sigma=0$ | on | on | 1 | 1 | on | on | on | on
| $\sigma=10^{-4}$ | on | on | 1 | 1 | on | on | on | on
| $\sigma=3\cdot 10^{-4}$ | on | on | 1 | 1 | on | on | on | on
| $\sigma=5\cdot 10^{-4}$ | on | on | 1 | 1 | on | on | on | on

Besides and simultaneously to identifying the material class, EUCLID computes
the unknown material parameters. Table 2 reports the comparison between the
true and the calibrated material parameters for the different benchmarks,
revealing a very good agreement. As expected, the discrepancy between true and
discovered parameters increases for increasing noise level.

Table 2: True and discovered material parameters for different noise levels
$\sigma$ (in mm).

Benchmarks |  $G$ $\left[\frac{\text{kN}}{\text{mm}^{2}}\right]$ |  $K$ $\left[\frac{\text{kN}}{\text{mm}^{2}}\right]$ |  $G_{1}$ $\left[\frac{\text{kN}}{\text{mm}^{2}}\right]$ |  $g_{1}$ $\left[\text{s}\right]$ |  $K_{1}$ $\left[\frac{\text{kN}}{\text{mm}^{2}}\right]$ |  $k_{1}$ $\left[\text{s}\right]$ |  $\sigma_{0}$ $\left[\frac{\text{kN}}{\text{mm}^{2}}\right]$ |  $\eta_{p}$ $\left[\frac{\text{kN}\ \text{s}}{\text{mm}^{2}}\right]$ |  $H^{iso}$ $\left[\frac{\text{kN}}{\text{mm}^{2}}\right]$ |  $H^{kin}$ $\left[\frac{\text{kN}}{\text{mm}^{2}}\right]$
---|---|---|---|---|---|---|---|---|---|---
E | Truth | 0.6000 | 1.3000 | 0 | $\rightarrow\infty$ | 0 | $\rightarrow\infty$ | $\rightarrow\infty$ | 0 | 0 | 0
| $\sigma=0$ | 0.6000 | 1.3000 | 0 | $\rightarrow\infty$ | 0 | $\rightarrow\infty$ | $\rightarrow\infty$ | 0 | 0 | 0
| $\sigma=10^{-4}$ | 0.6000 | 1.3000 | 0 | $\rightarrow\infty$ | 0 | $\rightarrow\infty$ | $\rightarrow\infty$ | 0 | 0 | 0
| $\sigma=3\cdot 10^{-4}$ | 0.5998 | 1.3004 | 0 | $\rightarrow\infty$ | 0 | $\rightarrow\infty$ | $\rightarrow\infty$ | 0 | 0 | 0
| $\sigma=5\cdot 10^{-4}$ | 0.5998 | 1.3002 | 0 | $\rightarrow\infty$ | 0 | $\rightarrow\infty$ | $\rightarrow\infty$ | 0 | 0 | 0
VE | Truth | 0.6000 | 1.3000 | 0.3500 | 110.0000 | 0.4000 | 15.0000 | $\rightarrow\infty$ | 0 | 0 | 0
| $\sigma=0$ | 0.6000 | 1.2999 | 0.3500 | 109.9953 | 0.3999 | 15.0047 | $\rightarrow\infty$ | 0 | 0 | 0
| $\sigma=10^{-4}$ | 0.6017 | 1.2930 | 0.3517 | 109.6126 | 0.3929 | 15.5063 | $\rightarrow\infty$ | 0 | 0 | 0
| $\sigma=3\cdot 10^{-4}$ | 0.6042 | 1.2674 | 0.3544 | 109.7273 | 0.3671 | 16.5166 | $\rightarrow\infty$ | 0 | 0 | 0
| $\sigma=5\cdot 10^{-4}$ | 0.6094 | 1.2702 | 0.3597 | 108.2126 | 0.3693 | 15.1302 | $\rightarrow\infty$ | 0 | 0 | 0
VEEP | Truth | 0.6000 | 1.3000 | 0.3500 | 110.0000 | 0 | $\rightarrow\infty$ | 0.0300 | 0 | 0.0300 | 0
| $\sigma=0$ | 0.6000 | 1.3000 | 0.3500 | 110.0008 | 0 | $\rightarrow\infty$ | 0.0300 | 0 | 0.0300 | 0
| $\sigma=10^{-4}$ | 0.6008 | 1.2997 | 0.3506 | 109.5033 | 0 | $\rightarrow\infty$ | 0.0303 | 0 | 0.0276 | 0
| $\sigma=3\cdot 10^{-4}$ | 0.5935 | 1.3055 | 0.3435 | 114.3310 | 0.0090 | 31.5015 | 0.0314 | 0 | 0.0213 | 0
| $\sigma=5\cdot 10^{-4}$ | 0.5908 | 1.2963 | 0.3409 | 116.6722 | 0 | $\rightarrow\infty$ | 0.0329 | 0 | 0.0123 | 0
EVP | Truth | 0.6000 | 1.3000 | 0 | $\rightarrow\infty$ | 0 | $\rightarrow\infty$ | 0.0300 | 0.0400 | 0 | 0.0100
| $\sigma=0$ | 0.6000 | 1.3000 | 0 | $\rightarrow\infty$ | 0 | $\rightarrow\infty$ | 0.0300 | 0.0399 | 0 | 0.0100
| $\sigma=10^{-4}$ | 0.6013 | 1.2997 | 0 | $\rightarrow\infty$ | 0 | $\rightarrow\infty$ | 0.0303 | 0.0360 | 0 | 0.0099
| $\sigma=3\cdot 10^{-4}$ | 0.6015 | 1.2991 | 0 | $\rightarrow\infty$ | 0 | $\rightarrow\infty$ | 0.0307 | 0.0288 | 0 | 0.0099
| $\sigma=5\cdot 10^{-4}$ | 0.6020 | 1.2965 | 0 | $\rightarrow\infty$ | 0 | $\rightarrow\infty$ | 0.0313 | 0.0229 | 0 | 0.0096
VEVP | Truth | 0.6000 | 1.3000 | 0.3500 | 110.0000 | 0.4000 | 15.0000 | 0.0300 | 0.0400 | 0.0300 | 0.0100
| $\sigma=0$ | 0.6000 | 1.3000 | 0.3500 | 110.0096 | 0.4000 | 15.0007 | 0.0300 | 0.0399 | 0.0300 | 0.0100
| $\sigma=10^{-4}$ | 0.5996 | 1.3011 | 0.3494 | 110.5654 | 0.4020 | 15.1755 | 0.0304 | 0.0321 | 0.0276 | 0.0098
| $\sigma=3\cdot 10^{-4}$ | 0.5970 | 1.2930 | 0.3466 | 113.9556 | 0.3965 | 14.6939 | 0.0316 | 0.0158 | 0.0201 | 0.0095
| $\sigma=5\cdot 10^{-4}$ | 0.5967 | 1.3001 | 0.3468 | 114.2066 | 0.4026 | 14.7275 | 0.0320 | 0.0160 | 0.0181 | 0.0094

Finally, in order to test if the discovered material models describe a similar
material response as the true ones under unseen loading histories (and
especially for the cases with noise, featuring some discrepancy between true
and discovered material parameters), we compute the stress response of the
material along two different arbitrarily chosen deformation paths, see Figure
4. We select these strain paths as uniaxial tension (UT) and simple shear
(SS), given by

$\bm{\varepsilon}^{\text{UT}}=\begin{bmatrix}\epsilon&0&0\\\ 0&0&0\\\ 0&0&0\\\
\end{bmatrix},\qquad\bm{\varepsilon}^{\text{SS}}=\begin{bmatrix}0&\epsilon&0\\\
\epsilon&0&0\\\ 0&0&0\\\ \end{bmatrix},$ (39)

where $\epsilon$ is a scalar deformation parameter that is let to increase to
a value of $0.1$ during ten loading phases and to subsequently decrease to
zero during additional ten loading phases. The loading and relaxation periods
of each loading phase, as well as the loading and relaxation times, are chosen
in the same way as described for the loading parameter $\delta$ in Section
4.1. Figure 4 demonstrates the good agreement between true and discovered
material responses, with discrepancies that, as expected, increase for
increasing noise level.

Figure 4: True and discovered material responses under uniaxial tension (UT)
and simple shear (SS) for different noise levels $\sigma$ (in mm). The stress
components $\sigma_{11}$ and $\sigma_{12}$ are in
$\frac{\text{kN}}{\text{mm}^{2}}$.

##  Section 5 Conclusions

We extended the scope of our recently developed EUCLID approach for
unsupervised automated discovery of material laws to the highly general case
of a material belonging to an unknown class of constitutive behavior. To this
end, we leveraged the theory of generalized standard materials, which
encompasses a large number of important constitutive classes including
elasticity, viscosity, plasticity and arbitrary combinations thereof. We
showed that, based solely on full-field displacements and net reaction forces,
EUCLID is able to automatically discover the two scalar thermodynamic
potentials, namely, the Helmholtz free energy and the dissipation potential,
which completely define the behavior of generalized standard materials.
Stability and thermodynamic consistency of the discovered model are guaranteed
by construction by the a priori requirement of convexity for the two
potentials; balance of linear momentum (enforced locally in the interior and
globally at each constrained side of the domain) acts as a fundamental
constraint to counteract the ill-posedness of the inverse problem due to the
unavailability of stress-strain labeled pairs; sparsity promoting
regularization leads to the automatic selection of a small subset from a
possibly large number of candidate model features in the initially chosen
large model library and thus leads to a parsimonious, i.e., simple and
interpretable, model. Importantly, since model features go hand in hand with
the corresponding internal variables, sparse regression automatically induces
a parsimonious selection of the few internal variables needed for an accurate
but simple description of the material behavior. We propose a fully automatic
procedure for the selection of the hyperparameter controlling the weight of
the sparsity promoting regularization term, in order to strike a user-defined
balance between model accuracy and simplicity.

The present study paves the way for a number of interesting extensions. A
possible important direction would be a further generalization of the model
library for the thermodynamic potentials to include additional phenomena such
as anisotropy, heterogeneity, softening (which would require the consideration
of gradient terms), coupled thermo- or electromechanical effects, to name but
a few. Another important aim of future research would be the enhancement of
the computational efficiency, as solving the constitutive equations becomes
more demanding for an increasing number of unknown parameters in the inverse
problem. To ease the implementation of the constitutive equation solver for a
wide library of thermodynamic potentials, relying on automatic differentiation
would be beneficial; however, as discussed in this paper, automatic
differentiation requires a regularization of the potentials and currently
implies a decrease in the overall computational efficiency, which could be
addressed in the future. Last but not least, testing of EUCLID on experimental
data would be an important future goal.

## Acknowledgments

MF and LDL would like to acknowledge funding by SNF through grant N.
$200021\\_204316$ “Unsupervised data-driven discovery of material laws”.

##  Appendix A Solving the constitutive equations

We solve the constitutive equations (see Section 3.1.2) in strain control.
Given the state variables and dependent variables at the previous time step
$(t-1)$ and the strain $\bm{\varepsilon}^{t}$ at the current time step $t$,
the objective is to compute first the internal variables $\bm{\alpha}^{t}$
using Equations (19), (22) and then the stress $\bm{\sigma}^{t}$ using
Equation (16). To this end, Equations (19), (22), (16) are discretized in time
using an implicit Euler scheme, i.e., $\dot{\square}^{t}=\frac{1}{\Delta
t}\left(\square^{t}-\square^{t-1}\right)$, where
$\Delta\square=\square^{t}-\square^{t-1}$ is used to denote differences of
variables between the previous and current time step. The discretized
equations are then solved with a viscoelastic predictor / viscoplastic
corrector algorithm (see Simo and Hughes (1998); Neto et al. (2008) for
similar algorithms).

###  A.1 Viscoelastic predictor step

For the viscoelastic predictor step (also known as trial step), we assume a
purely viscoelastic behavior, i.e., $f^{t,\text{trial}}\leq 0$. Under this
assumption, the viscoplastic internal variables in Equations (22) do not
evolve, e.g., $\bm{\alpha}_{{I}}^{t,\text{trial}}=\bm{\alpha}_{{I}}^{t-1}$.
The viscoelastic internal variables are updated using the time-discretized
Equation (19)

$\frac{1}{\Delta
t}\left(\bm{\alpha}_{j}^{t,\text{trial}}-\bm{\alpha}_{j}^{t-1}\right)=\frac{1}{g_{j}}\operatorname{dev}(\bm{\varepsilon}^{t}-\bm{\alpha}_{j}^{t,\text{trial}}-\bm{\alpha}_{{I}}^{t-1})+\frac{1}{k_{j}}\operatorname{vol}(\bm{\varepsilon}^{t}-\bm{\alpha}_{j}^{t,\text{trial}}-\bm{\alpha}_{{I}}^{t-1}).$
(40)

Splitting this equation into its deviatoric and volumetric contributions, we
obtain explicit expressions for the update of the viscoelastic internal
variables

$\displaystyle\operatorname{dev}\bm{\alpha}_{j}^{t,\text{trial}}$
$\displaystyle=\frac{1}{\left(1+\frac{\Delta
t}{g_{j}}\right)}\left(\operatorname{dev}\bm{\alpha}_{j}^{t-1}+\frac{\Delta
t}{g_{j}}\operatorname{dev}\left(\bm{\varepsilon}^{t}-\bm{\alpha}_{{I}}^{t-1}\right)\right),$
(41) $\displaystyle\operatorname{vol}\bm{\alpha}_{j}^{t,\text{trial}}$
$\displaystyle=\frac{1}{\left(1+\frac{\Delta
t}{k_{j}}\right)}\left(\operatorname{vol}\bm{\alpha}_{j}^{t-1}+\frac{\Delta
t}{k_{j}}\operatorname{vol}\left(\bm{\varepsilon}^{t}-\bm{\alpha}_{{I}}^{t-1}\right)\right).$

Knowing all internal variables at the current time step $t$, the trial stress
$\bm{\sigma}^{t,\text{trial}}$ is computed using Equation (16). After some
algebra, we arrive at the explicit expressions

$\displaystyle\operatorname{dev}\bm{\sigma}^{t,\text{trial}}$
$\displaystyle=2\bar{G}\operatorname{dev}\left(\bm{\varepsilon}^{t}-\bm{\alpha}_{{I}}^{t-1}\right)-\sum_{j=1}^{n_{\text{MW}}}2\frac{G_{j}}{1+\frac{\Delta
t}{g_{j}}}\operatorname{dev}\bm{\alpha}_{j}^{t-1},$ (42)
$\displaystyle\operatorname{vol}\bm{\sigma}^{t,\text{trial}}$
$\displaystyle=3\bar{K}\operatorname{vol}\left(\bm{\varepsilon}^{t}-\bm{\alpha}_{{I}}^{t-1}\right)-\sum_{j=1}^{n_{\text{MW}}}3\frac{K_{j}}{1+\frac{\Delta
t}{k_{j}}}\operatorname{vol}\bm{\alpha}_{j}^{t-1},$

where we defined

$\displaystyle\bar{G}$
$\displaystyle=G+\sum_{j=1}^{n_{\text{MW}}}G_{j}\left(1-\frac{\frac{\Delta
t}{g_{j}}}{1+\frac{\Delta t}{g_{j}}}\right),$ (43) $\displaystyle\bar{K}$
$\displaystyle=K+\sum_{j=1}^{n_{\text{MW}}}K_{j}\left(1-\frac{\frac{\Delta
t}{k_{j}}}{1+\frac{\Delta t}{k_{j}}}\right),$

We can determine whether the viscoelastic predictor is admissible by computing
$f^{t,\text{trial}}$ dependent on $\bm{\sigma}^{t,\text{trial}}$. If
$f^{t,\text{trial}}>0$, the viscoelastic predictor is not admissible and the
viscoplastic corrector step needs to be applied.

###  A.2 Viscoplastic corrector step

The computation of the viscoplastic internal variables in the viscoplastic
corrector step can be greatly simplified by taking advantage of some
characteristics of the problem. To this end, we define a scalar internal
variable $\gamma=\sqrt{\frac{3}{2}}\alpha_{{II}}$. After discretization in
time and noting that $f^{\text{{t,trial}}}>0$, the evolution laws in (22) can
then be written as

$\Delta{\bm{\alpha}}_{{I}}=\Delta{\bm{\alpha}}_{{III}}=\Delta\gamma\frac{\operatorname{dev}\left(\bm{\sigma}^{t}-H^{kin}\bm{\alpha}_{{III}}^{t}\right)}{\left\|\operatorname{dev}\left(\bm{\sigma}^{t}-H^{kin}\bm{\alpha}_{{III}}^{t}\right)\right\|},\qquad\Delta{{\alpha}}_{{II}}=\sqrt{\frac{2}{3}}\Delta\gamma.$
(44)

It will be shown that the internal variable update $\Delta\gamma$ can be
computed explicitly dependent on the variables at the previous time step, and
afterwards, the update of the other viscoplastic internal variables (i.e.,
$\Delta\bm{\alpha}_{{I}}=\Delta\bm{\alpha}_{{III}}$ and $\Delta\alpha_{{II}}$)
can be computed dependent on $\Delta\gamma$. The explicit expression for
$\Delta\gamma$ and the relation between $\Delta\gamma$ and the other internal
variable updates are derived in the following. The availability of these
relations leads to a computationally efficient implementation of the
viscoplastic corrector step.

With reference to (42), the deviatoric and volumetric parts of the stress at
the current time step are

$\displaystyle\operatorname{dev}\bm{\sigma}^{t}$
$\displaystyle=2\bar{G}\operatorname{dev}\left(\bm{\varepsilon}^{t}-\left(\bm{\alpha}_{{I}}^{t,\text{trial}}+\Delta\bm{\alpha}_{{I}}\right)\right)-\sum_{j=1}^{n_{\text{MW}}}2\frac{G_{j}}{1+\frac{\Delta
t}{g_{j}}}\operatorname{dev}\bm{\alpha}_{j}^{t-1},$ (45)
$\displaystyle\operatorname{vol}\bm{\sigma}^{t}$
$\displaystyle=3\bar{K}\operatorname{vol}\left(\bm{\varepsilon}^{t}-\left(\bm{\alpha}_{{I}}^{t,\text{trial}}+\Delta\bm{\alpha}_{{I}}\right)\right)-\sum_{j=1}^{n_{\text{MW}}}3\frac{K_{j}}{1+\frac{\Delta
t}{k_{j}}}\operatorname{dev}\bm{\alpha}_{j}^{t-1},$

where we made use of
$\bm{\alpha}_{{I}}^{t}=\bm{\alpha}_{{I}}^{t,\text{trial}}+\Delta\bm{\alpha}_{{I}}$.
Once again referring to (42), we write (45) as

$\displaystyle\operatorname{dev}\bm{\sigma}^{t}$
$\displaystyle=\operatorname{dev}\bm{\sigma}^{t,\text{trial}}-2\bar{G}\operatorname{dev}\Delta\bm{\alpha}_{{I}},$
(46) $\displaystyle\operatorname{vol}\bm{\sigma}^{t}$
$\displaystyle=\operatorname{vol}\bm{\sigma}^{t,\text{trial}}-3\bar{K}\operatorname{vol}\Delta\bm{\alpha}_{{I}},$

Knowing that the volumetric part of $\Delta\bm{\alpha}_{{I}}$ vanishes (see
(44)), we observe that
$\operatorname{vol}\bm{\sigma}^{t}=\operatorname{vol}\bm{\sigma}^{t,\text{trial}}$.
Therefore, we can focus on the stress deviator in the following. Subtracting
the term $\operatorname{dev}\left(H^{kin}\bm{\alpha}_{{III}}^{t}\right)$ from
both sides of the first equation in (46) and using
$\Delta\bm{\alpha}_{{I}}=\Delta\bm{\alpha}_{{III}}$ gives

$\operatorname{dev}\left(\bm{\sigma}^{t}-H^{kin}\bm{\alpha}_{{III}}^{t}\right)=\operatorname{dev}\left(\bm{\sigma}^{t,\text{trial}}-H^{kin}\bm{\alpha}_{{III}}^{t,\text{trial}}\right)-(2\bar{G}+H^{kin})\operatorname{dev}\Delta\bm{\alpha}_{{I}}.$
(47)

After inserting the time-discretized evolution law for
$\Delta\bm{\alpha}_{{I}}$ (see (44)), we obtain

$\operatorname{dev}\left(\bm{\sigma}^{t}-H^{kin}\bm{\alpha}_{{III}}^{t}\right)=\operatorname{dev}\left(\bm{\sigma}^{t,\text{trial}}-H^{kin}\bm{\alpha}_{{III}}^{t,\text{trial}}\right)-(2\bar{G}+H^{kin})\Delta\gamma\frac{\operatorname{dev}\left(\bm{\sigma}^{t}-H^{kin}\bm{\alpha}_{{III}}^{t}\right)}{\left\|\operatorname{dev}\left(\bm{\sigma}^{t}-H^{kin}\bm{\alpha}_{{III}}^{t}\right)\right\|},$
(48)

which can be rearranged to yield

$\left(1+\frac{(2\bar{G}+H^{kin})\Delta\gamma}{\left\|\operatorname{dev}\left(\bm{\sigma}^{t}-H^{kin}\bm{\alpha}_{{III}}^{t}\right)\right\|}\right)\operatorname{dev}\left(\bm{\sigma}^{t}-H^{kin}\bm{\alpha}_{{III}}^{t}\right)=\operatorname{dev}\left(\bm{\sigma}^{t,\text{trial}}-H^{kin}\bm{\alpha}_{{III}}^{t,\text{trial}}\right).$
(49)

This implies that the tensors
$\operatorname{dev}\left(\bm{\sigma}^{t}-H^{kin}\bm{\alpha}_{{III}}^{t}\right)$
and
$\operatorname{dev}\left(\bm{\sigma}^{t,\text{trial}}-H^{kin}\bm{\alpha}_{{III}}^{t,\text{trial}}\right)$
only differ by a scalar factor and hence are collinear, i.e.

$\frac{\operatorname{dev}\left(\bm{\sigma}^{t}-H^{kin}\bm{\alpha}_{{III}}^{t}\right)}{\left\|\operatorname{dev}\left(\bm{\sigma}^{t}-H^{kin}\bm{\alpha}_{{III}}^{t}\right)\right\|}=\frac{\operatorname{dev}\left(\bm{\sigma}^{t,\text{trial}}-H^{kin}\bm{\alpha}_{{III}}^{t,\text{trial}}\right)}{\left\|\operatorname{dev}\left(\bm{\sigma}^{t,\text{trial}}-H^{kin}\bm{\alpha}_{{III}}^{t,\text{trial}}\right)\right\|}.$
(50)

The importance of this result is best understood when revisiting (44). With
this new result, the update of the internal variables
$\Delta\bm{\alpha}_{{I}}$ and $\Delta\bm{\alpha}_{{III}}$ can be written as a
function of $\Delta\gamma$ and of known variables from the previous time step.
Hence, it only remains to compute $\Delta\gamma$, as derived in the following.

Inserting (50) in (48), we obtain

$\operatorname{dev}\left(\bm{\sigma}^{t}-H^{kin}\bm{\alpha}_{{III}}^{t}\right)=\left(1-\frac{(2\bar{G}+H^{kin})\Delta\gamma}{\left\|\operatorname{dev}\left(\bm{\sigma}^{t,\text{trial}}-H^{kin}\bm{\alpha}_{{III}}^{t,\text{trial}}\right)\right\|}\right)\operatorname{dev}\left(\bm{\sigma}^{t,\text{trial}}-H^{kin}\bm{\alpha}_{{III}}^{t,\text{trial}}\right).$
(51)

Using this result, we can write

$\displaystyle\Delta\gamma$ $\displaystyle=\sqrt{\frac{3}{2}}\frac{\Delta
t}{\eta_{p}}f^{t},$ (52) $\displaystyle=\sqrt{\frac{3}{2}}\frac{\Delta
t}{\eta_{p}}\left(\sqrt{\frac{3}{2}}\left\|\operatorname{dev}\left(\bm{\sigma}^{t}-H^{kin}\bm{\alpha}_{{III}}^{t}\right)\right\|-\sigma_{y}-H^{iso}\alpha_{{II}}^{t,\text{trial}}-H^{iso}\sqrt{\frac{2}{3}}\Delta\gamma\right),$
$\displaystyle=\sqrt{\frac{3}{2}}\frac{\Delta
t}{\eta_{p}}\left(\sqrt{\frac{3}{2}}\left\|\operatorname{dev}\left(\bm{\sigma}^{t,\text{trial}}-H^{kin}\bm{\alpha}_{{III}}^{t,\text{trial}}\right)\right\|-\sqrt{\frac{3}{2}}(2\bar{G}+H^{kin})\Delta\gamma-\sigma_{y}-H^{iso}\alpha_{{II}}^{t,\text{trial}}-H^{iso}\sqrt{\frac{2}{3}}\Delta\gamma\right),$
$\displaystyle=\sqrt{\frac{3}{2}}\frac{\Delta
t}{\eta_{p}}\left(f^{t,\text{trial}}-\sqrt{\frac{3}{2}}(2\bar{G}+H^{kin})\Delta\gamma-H^{iso}\sqrt{\frac{2}{3}}\Delta\gamma\right),$

which is solved for $\Delta\gamma$ to obtain

$\Delta\gamma=\frac{f^{t,\text{trial}}}{\sqrt{\frac{2}{3}}\left(\frac{\eta_{p}}{\Delta
t}+H^{iso}\right)+\sqrt{\frac{3}{2}}(2\bar{G}+H^{kin})}.$ (53)

After computing $\Delta\gamma$, the other internal variables (see (44) and
(50)) and the stress (see (45)) can be computed.

##  Appendix B Consistent tangent

Implementing the material model library in forward finite element simulations
requires the computation of the consistent tangent

$\mathbb{C}^{t}=\frac{\partial\bm{\sigma}^{t}}{\partial\bm{\varepsilon}^{t}}.$
(54)

We consider two different cases. If the viscoelastic predictor in A.1 is
admissible, it is

$\mathbb{C}^{t}=\mathbb{C}^{t,\text{trial}}=\frac{\partial\bm{\sigma}^{t,\text{trial}}}{\partial\bm{\varepsilon}^{t}}=2\bar{G}\mathbb{I}^{\text{dev}}+3\bar{K}\mathbb{I}^{\text{vol}}.$
(55)

Otherwise, we obtain

$\displaystyle\mathbb{C}^{t}$
$\displaystyle=\left(1-\frac{2\bar{G}\Delta\gamma}{\left\|\operatorname{dev}\left(\bm{\sigma}^{t,\text{trial}}-H^{kin}\bm{\alpha}_{{III}}^{t,\text{trial}}\right)\right\|}\right)2\bar{G}\mathbb{I}^{\text{dev}}+3\bar{K}\mathbb{I}^{\text{vol}}$
(56)
$\displaystyle-\frac{2\bar{G}\Delta\gamma\left(\sigma_{y}+H^{iso}\alpha_{{II}}^{t,\text{trial}}\right)}{f^{t,\text{trial}}\left\|\operatorname{dev}\left(\bm{\sigma}^{t,\text{trial}}-H^{kin}\bm{\alpha}_{{III}}^{t,\text{trial}}\right)\right\|^{3}}\mathbb{C}^{t,\text{trial}}\colon\left(\operatorname{dev}\left(\bm{\sigma}^{t,\text{trial}}-H^{kin}\bm{\alpha}_{{III}}^{t,\text{trial}}\right)\otimes\operatorname{dev}\left(\bm{\sigma}^{t,\text{trial}}-H^{kin}\bm{\alpha}_{{III}}^{t,\text{trial}}\right)\right).$

Note that the consistent tangent is required for the solution of the forward
finite element problem but not for the inverse problem, i.e., not for EUCLID.

##  Appendix C Numerical settings

Table 3 provides a list of parameters and hyperparameters used during the data
generation via finite element simulations and during the inverse discovery
process via EUCLID.

Table 3: (Hyper-)parameters for the finite element simulations and EUCLID. Parameter | Notation | Value | Unit
---|---|---|---
Number of nodes in mesh | $n_{n}$ | $2,179$ | -
Number of measured reaction force components | $n_{\beta}$ | $2$ | -
Number of time steps | $n_{t}$ | $800$ | -
Number of loading phases | $n_{k}$ | $20$ | -
Applied loading per period | $\delta$ | $\pm 0.5$ | mm
Time per loading period | $t^{\text{load}}_{k}$ | $\\{10^{(-2+\frac{4}{n_{k}-1}i)}:i=0,\dots,(n_{k}-1)\\}$ | s
Time per relaxation period | $t^{\text{relax}}_{k}$ | $10^{2}$ | s
Displacement noise standard deviation | $\sigma$ | $\left\\{0,10^{-4},3\cdot 10^{-4},5\cdot 10^{-4}\right\\}$ | mm
Moving-window length for denoising | - | $10$ | -
Coefficient for reaction force balance | $\lambda_{r}$ | 100 | -
Number of random initial guesses | $n_{g}$ | $24$ | -
Exponent for $L_{p}$ regularization | $p$ | $1$ | -
Number of choices of $\lambda_{p}$ | $n_{\lambda}$ | $24$ | -
Coefficient for $L_{p}$ regularization | $\lambda_{p}$ | $\\{10^{-4}\cdot 2^{i}:i=0,\dots,23\\}$ | $\text{kN}^{2}$
Lower bound for $\bm{\theta}$ | $\theta^{\text{min}}_{i}$ | $\\{0,\theta^{\text{lb}}\\}$ | various
Lower bound for $\bm{\theta}$ (reciprocals) | $\theta^{\text{lb}}$ | $10^{-6}$ | various
Threshold for $C$ | $C^{\text{th}}$ | $\max\\{10^{-5},1.1\ C^{\text{min}}\\}$ | $\text{kN}^{2}$
Threshold for $\bm{\theta}$ | $\theta^{\text{th}}$ | $10^{-4}$ | various

## Code and data availability

Codes and data are publicly available at https://euclid-code.github.io/.

## References

  * Amores et al. (2022) Amores, V.J., Montáns, F.J., Cueto, E., Chinesta, F., 2022\.  Crossing Scales: Data-Driven Determination of the Micro-scale Behavior of Polymers From Non-homogeneous Tests at the Continuum-Scale.  Frontiers in Materials 9, 879614\.  URL: https://www.frontiersin.org/articles/10.3389/fmats.2022.879614/full, doi:10.3389/fmats.2022.879614.
  * Anton and Wessels (2022) Anton, D., Wessels, H., 2022\.  Identification of Material Parameters from Full-Field Displacement Data Using Physics-Informed Neural Networks .
  * As’ad et al. (2022) As’ad, F., Avery, P., Farhat, C., 2022.  A Mechanics-Informed Artificial Neural Network Approach in Data-Driven Constitutive Modeling, in: AIAA SCITECH 2022 Forum, American Institute of Aeronautics and Astronautics, San Diego, CA & Virtual.  URL: https://arc.aiaa.org/doi/10.2514/6.2022-0100, doi:10.2514/6.2022-0100.
  * Avril et al. (2008) Avril, S., Bonnet, M., Bretelle, A.S., Grédiac, M., Hild, F., Ienny, P., Latourte, F., Lemosse, D., Pagano, S., Pagnacco, E., Pierron, F., 2008.  Overview of Identification Methods of Mechanical Parameters Based on Full-field Measurements.  Experimental Mechanics 48, 381–402.  URL: http://link.springer.com/10.1007/s11340-008-9148-y, doi:10.1007/s11340-008-9148-y.
  * Blühdorn et al. (2022) Blühdorn, J., Gauger, N.R., Kabel, M., 2022.  AutoMat: automatic differentiation for generalized standard materials on GPUs.  Computational Mechanics 69, 589–613.  URL: https://link.springer.com/10.1007/s00466-021-02105-2, doi:10.1007/s00466-021-02105-2.
  * Brunton et al. (2016) Brunton, S.L., Proctor, J.L., Kutz, J.N., 2016.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the National Academy of Sciences 113, 3932–3937.  URL: http://www.pnas.org/lookup/doi/10.1073/pnas.1517384113, doi:10.1073/pnas.1517384113.
  * Bröcker and Matzenmiller (2012) Bröcker, C., Matzenmiller, A., 2012\.  Thermoviscoplasticity deduced from enhanced rheological models.  PAMM 12, 327–328.  URL: https://onlinelibrary.wiley.com/doi/10.1002/pamm.201210152, doi:10.1002/pamm.201210152.
  * Dalémat et al. (2019) Dalémat, M., Coret, M., Leygue, A., Verron, E., 2019\.  Measuring stress field without constitutive equation.  Mechanics of Materials 136, 103087\.  doi:10.1016/j.mechmat.2019.103087.
  * Doi (2011) Doi, M., 2011.  Onsager’s variational principle in soft matter.  Journal of Physics: Condensed Matter 23, 284118.  URL: https://iopscience.iop.org/article/10.1088/0953-8984/23/28/284118, doi:10.1088/0953-8984/23/28/284118.
  * Flaschel et al. (2021) Flaschel, M., Kumar, S., De Lorenzis, L., 2021.  Unsupervised discovery of interpretable hyperelastic constitutive laws.  Computer Methods in Applied Mechanics and Engineering 381, 113852.  doi:10.1016/j.cma.2021.113852.
  * Flaschel et al. (2022) Flaschel, M., Kumar, S., De Lorenzis, L., 2022.  Discovering plasticity models without stress data.  npj Computational Materials 8, 91\.  URL: https://www.nature.com/articles/s41524-022-00752-4, doi:10.1038/s41524-022-00752-4.
  * Frank and Friedman (1993) Frank, l.E., Friedman, J.H., 1993\.  A Statistical View of Some Chemometrics Regression Tools.  Technometrics 35, 109–135.  URL: http://www.tandfonline.com/doi/abs/10.1080/00401706.1993.10485033, doi:10.1080/00401706.1993.10485033.
  * Frankel et al. (2020) Frankel, A.L., Jones, R.E., Swiler, L.P., 2020.  TENSOR BASIS GAUSSIAN PROCESS MODELS OF HYPERELASTIC MATERIALS.  Journal of Machine Learning for Modeling and Computing 1, 1–17.  URL: http://www.dl.begellhouse.com/journals/558048804a15188a,583c4e56625ba94e,651a2e6b0260f708.html, doi:10.1615/JMachLearnModelComput.2020033325.
  * Fuhg et al. (2021) Fuhg, J.N., Boehm, C., Bouklas, N., Fau, A., Wriggers, P., Marino, M., 2021\.  Model-data-driven constitutive responses: application to a multiscale computational framework.  arXiv:2104.02650 [cs, math, stat] URL: http://arxiv.org/abs/2104.02650. arXiv: 2104.02650.
  * Fuhg et al. (2022) Fuhg, J.N., Fau, A., Bouklas, N., Marino, M., 2022\.  Elasto-plasticity with convex model-data-driven yield functions , 27.
  * Ghaboussi et al. (1991) Ghaboussi, J., Garrett, J.H., Wu, X., 1991.  Knowledge‐Based Modeling of Material Behavior with Neural Networks.  Journal of Engineering Mechanics 117, 132–153.  URL: http://ascelibrary.org/doi/10.1061/%28ASCE%290733-9399%281991%29117%3A1%28132%29, doi:10.1061/(ASCE)0733-9399(1991)117:1(132).
  * Giunta and Angela Pisano (2006) Giunta, M., Angela Pisano, A., 2006\.  One‐Dimensional Visco‐Elastoplastic Constitutive Model for Asphalt Concrete.  Multidiscipline Modeling in Materials and Structures 2, 247–264.  URL: https://www.emerald.com/insight/content/doi/10.1163/157361106776240761/full/html, doi:10.1163/157361106776240761.
  * González et al. (2019a) González, D., Chinesta, F., Cueto, E., 2019a.  Learning Corrections for Hyperelastic Models From Data.  Frontiers in Materials 6, 14\.  URL: https://www.frontiersin.org/article/10.3389/fmats.2019.00014/full, doi:10.3389/fmats.2019.00014.
  * González et al. (2019b) González, D., Chinesta, F., Cueto, E., 2019b.  Thermodynamically consistent data-driven computational mechanics.  Continuum Mechanics and Thermodynamics 31, 239–253.  URL: http://link.springer.com/10.1007/s00161-018-0677-z, doi:10.1007/s00161-018-0677-z.
  * Grmela and Öttinger (1997) Grmela, M., Öttinger, H.C., 1997\.  Dynamics and thermodynamics of complex fluids. I. Development of a general formalism.  Physical Review E 56, 6620–6632.  URL: https://link.aps.org/doi/10.1103/PhysRevE.56.6620, doi:10.1103/PhysRevE.56.6620.
  * Grédiac (1989) Grédiac, M., 1989.  Principle of virtual work and identification.  Comptes Rendus de L Academie des Sciences Serie Ii , 1–5.
  * Halphen and Nguyen (1975) Halphen, B., Nguyen, Q.S., 1975\.  Sur les matériaux standard généralisés , 26\.
  * Hartmann (2001) Hartmann, S., 2001.  Numerical studies on the identification of the material parameters of Rivlin’s hyperelasticity using tension-torsion tests.  Acta Mechanica 148, 129–155.  URL: http://link.springer.com/10.1007/BF01183674, doi:10.1007/BF01183674.
  * Hernández et al. (2021) Hernández, Q., Badías, A., González, D., Chinesta, F., Cueto, E., 2021.  Structure-preserving neural networks.  Journal of Computational Physics 426, 109950.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0021999120307245, doi:10.1016/j.jcp.2020.109950.
  * Hild and Roux (2006) Hild, F., Roux, S., 2006.  Digital Image Correlation: from Displacement Measurement to Identification of Elastic Properties - a Review.  Strain 42, 69–80.  URL: http://doi.wiley.com/10.1111/j.1475-1305.2006.00258.x, doi:10.1111/j.1475-1305.2006.00258.x.
  * Huang et al. (2020) Huang, D.Z., Xu, K., Farhat, C., Darve, E., 2020\.  Learning constitutive relations from indirect observations using deep neural networks.  Journal of Computational Physics 416, 109491.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0021999120302655, doi:10.1016/j.jcp.2020.109491.
  * Huang et al. (2022) Huang, S., He, Z., Chem, B., Reina, C., 2022.  Variational Onsager Neural Networks (VONNs): A thermodynamics-based variational learning strategy for non-equilibrium PDEs.  Journal of the Mechanics and Physics of Solids 163, 104856.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0022509622000692, doi:10.1016/j.jmps.2022.104856.
  * Hütter and Svendsen (2011) Hütter, M., Svendsen, B., 2011\.  On the Formulation of Continuum Thermodynamic Models for Solids as General Equations for Non-equilibrium Reversible-Irreversible Coupling , 12.
  * Ibañez et al. (2018) Ibañez, R., Abisset-Chavanne, E., Aguado, J.V., Gonzalez, D., Cueto, E., Chinesta, F., 2018\.  A Manifold Learning Approach to Data-Driven Computational Elasticity and Inelasticity.  Archives of Computational Methods in Engineering 25, 47–57.  URL: http://link.springer.com/10.1007/s11831-016-9197-9, doi:10.1007/s11831-016-9197-9.
  * Ibáñez et al. (2019) Ibáñez, R., Abisset-Chavanne, E., González, D., Duval, J.L., Cueto, E., Chinesta, F., 2019\.  Hybrid constitutive modeling: data-driven learning of corrections to plasticity models.  International Journal of Material Forming 12, 717–725.  URL: http://link.springer.com/10.1007/s12289-018-1448-x, doi:10.1007/s12289-018-1448-x.
  * Joshi et al. (2022) Joshi, A., Thakolkaran, P., Zheng, Y., Escande, M., Flaschel, M., De Lorenzis, L., Kumar, S., 2022.  Bayesian-EUCLID: Discovering hyperelastic material laws with uncertainties.  Computer Methods in Applied Mechanics and Engineering 398, 115225.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782522003681, doi:10.1016/j.cma.2022.115225.
  * Kalina et al. (2022) Kalina, K.A., Linden, L., Brummund, J., Metsch, P., Kästner, M., 2022.  Automated constitutive modeling of isotropic hyperelasticity based on artificial neural networks.  Computational Mechanics 69, 213–232.  URL: https://link.springer.com/10.1007/s00466-021-02090-6, doi:10.1007/s00466-021-02090-6.
  * Kirchdoerfer and Ortiz (2016) Kirchdoerfer, T., Ortiz, M., 2016\.  Data-driven computational mechanics.  Computer Methods in Applied Mechanics and Engineering 304, 81–101.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782516300238, doi:10.1016/j.cma.2016.02.001.
  * Klein et al. (2022a) Klein, D.K., Fernández, M., Martin, R.J., Neff, P., Weeger, O., 2022a.  Polyconvex anisotropic hyperelasticity with neural networks.  Journal of the Mechanics and Physics of Solids 159, 104703.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0022509621003215, doi:10.1016/j.jmps.2021.104703.
  * Klein et al. (2022b) Klein, D.K., Ortigosa, R., Martínez-Frutos, J., Weeger, O., 2022b.  Finite electro-elasticity with physics-augmented neural networks.  Computer Methods in Applied Mechanics and Engineering 400, 115501.  URL: https://linkinghub.elsevier.com/retrieve/pii/S004578252200514X, doi:10.1016/j.cma.2022.115501.
  * Korelc (2002) Korelc, J., 2002.  Multi-language and Multi-environment Generation of Nonlinear Finite Element Codes.  Engineering with Computers 18, 312–327.  URL: http://link.springer.com/10.1007/s003660200028, doi:10.1007/s003660200028.
  * Kumar and Lopez-Pamies (2016) Kumar, A., Lopez-Pamies, O., 2016\.  On the two-potential constitutive modeling of rubber viscoelastic materials.  Comptes Rendus Mécanique 344, 102–112.  URL: https://linkinghub.elsevier.com/retrieve/pii/S1631072115001448, doi:10.1016/j.crme.2015.11.004.
  * Kumar and Kochmann (2022) Kumar, S., Kochmann, D.M., 2022\.  What Machine Learning Can Do for Computational Solid Mechanics, in: Aldakheel, F., Hudobivnik, B., Soleimani, M., Wessels, H., Weißenfels, C., Marino, M. (Eds.), Current Trends and Open Problems in Computational Mechanics. Springer International Publishing, Cham, pp. 275–285.  URL: https://link.springer.com/10.1007/978-3-030-87312-7_27, doi:10.1007/978-3-030-87312-7_27.
  * Lemaitre and Chaboche (1994) Lemaitre, J., Chaboche, J.L., 1994\.  Mechanics of solid materials.  Cambridge university press.
  * Leygue et al. (2018) Leygue, A., Coret, M., Réthoré, J., Stainier, L., Verron, E., 2018.  Data-based derivation of material response.  Computer Methods in Applied Mechanics and Engineering 331, 184–196.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782517307156, doi:10.1016/j.cma.2017.11.013.
  * Liu et al. (2020) Liu, X., Tao, F., Du, H., Yu, W., Xu, K., 2020.  Learning Nonlinear Constitutive Laws Using Neural Network Models Based on Indirectly Measurable Data.  Journal of Applied Mechanics 87, 081003\.  URL: https://asmedigitalcollection.asme.org/appliedmechanics/article/doi/10.1115/1.4047036/1083320/Learning-Nonlinear-Constitutive-Laws-Using-Neural, doi:10.1115/1.4047036.
  * Liu et al. (2021) Liu, X., Tian, S., Tao, F., Du, H., Yu, W., 2021\.  Machine learning-assisted modeling of composite materials and structures: a review, in: AIAA Scitech 2021 Forum, American Institute of Aeronautics and Astronautics, VIRTUAL EVENT.  URL: https://arc.aiaa.org/doi/10.2514/6.2021-2023, doi:10.2514/6.2021-2023.
  * Man and Furukawa (2011) Man, H., Furukawa, T., 2011\.  Neural network constitutive modelling for non-linear characterization of anisotropic materials.  International Journal for Numerical Methods in Engineering 85, 939–957.  URL: https://onlinelibrary.wiley.com/doi/10.1002/nme.2999, doi:10.1002/nme.2999.
  * Marek et al. (2019) Marek, A., Davis, F.M., Rossi, M., Pierron, F., 2019\.  Extension of the sensitivity-based virtual fields to large deformation anisotropic plasticity.  International Journal of Material Forming 12, 457–476.  URL: http://link.springer.com/10.1007/s12289-018-1428-1, doi:10.1007/s12289-018-1428-1.
  * Masi and Stefanou (2022) Masi, F., Stefanou, I., 2022\.  Evolution TANN and the discovery of the internal variables and evolution equations in solid mechanics.  URL: http://arxiv.org/abs/2209.13269. arXiv:2209.13269 [cond-mat].
  * Masi et al. (2021) Masi, F., Stefanou, I., Vannucci, P., Maffi-Berthier, V., 2021\.  Thermodynamics-based Artificial Neural Networks for constitutive modeling.  Journal of the Mechanics and Physics of Solids 147, 104277.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0022509620304841, doi:10.1016/j.jmps.2020.104277.
  * Michel and Suquet (2016) Michel, J.C., Suquet, P., 2016\.  A model-reduction approach in micromechanics of materials preserving the variational structure of constitutive relations.  Journal of the Mechanics and Physics of Solids 90, 254–285.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0022509616300928, doi:10.1016/j.jmps.2016.02.005.
  * Miehe (2011) Miehe, C., 2011.  A multi-field incremental variational framework for gradient-extended standard dissipative solids.  J. Mech. Phys. Solids , 26.
  * Mielke (2011) Mielke, A., 2011.  Formulation of thermoelastic dissipative material behavior using GENERIC.  Continuum Mechanics and Thermodynamics 23, 233–256.  URL: http://link.springer.com/10.1007/s00161-010-0179-0, doi:10.1007/s00161-010-0179-0.
  * Miled (2011) Miled, B., 2011.  Coupled viscoelastic-viscoplastic modeling of homogeneous and reinforced thermoplastic polymers.  Ph.D. thesis.
  * Miled et al. (2011) Miled, B., Doghri, I., Delannay, L., 2011.  Coupled viscoelastic–viscoplastic modeling of homogeneous and isotropic polymers: Numerical algorithm and analytical solutions.  Computer Methods in Applied Mechanics and Engineering 200, 3381–3394.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782511002702, doi:10.1016/j.cma.2011.08.015.
  * Montáns et al. (2019) Montáns, F.J., Chinesta, F., Gómez-Bombarelli, R., Kutz, J.N., 2019\.  Data-driven modeling and learning in science and engineering.  Comptes Rendus Mécanique 347, 845–855.  URL: https://linkinghub.elsevier.com/retrieve/pii/S1631072119301809, doi:10.1016/j.crme.2019.11.009.
  * Neggers et al. (2018) Neggers, J., Allix, O., Hild, F., Roux, S., 2018.  Big Data in Experimental Mechanics and Model Order Reduction: Today’s Challenges and Tomorrow’s Opportunities.  Archives of Computational Methods in Engineering 25, 143–164.  URL: http://link.springer.com/10.1007/s11831-017-9234-3, doi:10.1007/s11831-017-9234-3.
  * Neto et al. (2008) Neto, E.d.S., Peric, D., Owen, D., 2008.  Computational methods for plasticity.  John Wiley & Sons.
  * Nguyen (2010) Nguyen, Q.S., 2010.  On standard dissipative gradient models.  Annals of Solid and Structural Mechanics 1, 79–86.  URL: http://link.springer.com/10.1007/s12356-010-0006-0, doi:10.1007/s12356-010-0006-0.
  * Onsager (1931a) Onsager, L., 1931a.  Reciprocal Relations in Irreversible Processes. I.  Physical Review 37, 405–426.  URL: https://link.aps.org/doi/10.1103/PhysRev.37.405, doi:10.1103/PhysRev.37.405.
  * Onsager (1931b) Onsager, L., 1931b.  Reciprocal Relations in Irreversible Processes. II.  Physical Review 38, 2265–2279.  URL: https://link.aps.org/doi/10.1103/PhysRev.38.2265, doi:10.1103/PhysRev.38.2265.
  * Peng et al. (2021) Peng, G.C.Y., Alber, M., Buganza Tepole, A., Cannon, W.R., De, S., Dura-Bernal, S., Garikipati, K., Karniadakis, G., Lytton, W.W., Perdikaris, P., Petzold, L., Kuhl, E., 2021\.  Multiscale Modeling Meets Machine Learning: What Can We Learn?  Archives of Computational Methods in Engineering 28, 1017–1037.  URL: https://link.springer.com/10.1007/s11831-020-09405-5, doi:10.1007/s11831-020-09405-5.
  * Pierron et al. (2010) Pierron, F., Avril, S., The Tran, V., 2010.  Extension of the virtual fields method to elasto-plastic material identification with cyclic loads and kinematic hardening.  International Journal of Solids and Structures 47, 2993–3010.
  * Pierron and Grédiac (2012) Pierron, F., Grédiac, M., 2012\.  The Virtual Fields Method.  Springer New York, New York, NY.  URL: http://link.springer.com/10.1007/978-1-4614-1824-5, doi:10.1007/978-1-4614-1824-5.
  * Rockafellar (2015) Rockafellar, R.T., 2015.  Convex Analysis.  Princeton University Press.  URL: https://doi.org/10.1515/9781400873173, doi:10.1515/9781400873173.
  * Rothe and Hartmann (2015) Rothe, S., Hartmann, S., 2015\.  Automatic differentiation for stress and consistent tangent computation.  Archive of Applied Mechanics 85, 1103–1125.  URL: http://link.springer.com/10.1007/s00419-014-0939-6, doi:10.1007/s00419-014-0939-6.
  * Roux and Hild (2020) Roux, S., Hild, F., 2020.  Optimal procedure for the identification of constitutive parameters from experimentally measured displacement fields.  International Journal of Solids and Structures 184, 14–23.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0020768318304542, doi:10.1016/j.ijsolstr.2018.11.008.
  * Savitzky and Golay (1964) Savitzky, A., Golay, M.J.E., 1964\.  Smoothing and Differentiation of Data by Simplified Least Squares Procedures.  Analytical Chemistry 36, 1627–1639.  URL: https://pubs.acs.org/doi/abs/10.1021/ac60214a047, doi:10.1021/ac60214a047.
  * Schmidt and Lipson (2009) Schmidt, M., Lipson, H., 2009\.  Distilling Free-Form Natural Laws from Experimental Data.  Science 324, 81–85.  URL: https://www.sciencemag.org/lookup/doi/10.1126/science.1165893, doi:10.1126/science.1165893.
  * Simo and Hughes (1998) Simo, J.C., Hughes, T.J.R., 1998\.  Computational inelasticity.  Number v. 7 in Interdisciplinary applied mathematics, Springer, New York.
  * Stainier et al. (2019) Stainier, L., Leygue, A., Ortiz, M., 2019.  Model-free data-driven methods in mechanics: material data identification and solvers.  Computational Mechanics 64, 381–393.  URL: http://link.springer.com/10.1007/s00466-019-01731-1, doi:10.1007/s00466-019-01731-1.
  * Steinmann and Runesson (2021) Steinmann, P., Runesson, K., 2021\.  The Catalogue of Computational Material Models: Basic Geometrically Linear Models in 1D.  Springer International Publishing, Cham.  URL: http://link.springer.com/10.1007/978-3-030-63684-5, doi:10.1007/978-3-030-63684-5.
  * Sun and Zhu (2013) Sun, L., Zhu, Y., 2013.  A serial two-stage viscoelastic–viscoplastic constitutive model with thermodynamical consistency for characterizing time-dependent deformation behavior of asphalt concrete mixtures.  Construction and Building Materials 40, 584–595.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0950061812007489, doi:10.1016/j.conbuildmat.2012.10.004.
  * Sussman and Bathe (2009) Sussman, T., Bathe, K.J., 2009\.  A model of incompressible isotropic hyperelastic material behavior using spline interpolations of tension-compression test data.  Communications in Numerical Methods in Engineering 25, 53–63.  URL: http://doi.wiley.com/10.1002/cnm.1105, doi:10.1002/cnm.1105.
  * Tac et al. (2022) Tac, V., Sahli Costabal, F., Tepole, A.B., 2022.  Data-driven tissue mechanics with polyconvex neural ordinary differential equations.  Computer Methods in Applied Mechanics and Engineering 398, 115248.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782522003838, doi:10.1016/j.cma.2022.115248.
  * Thakolkaran et al. (2022) Thakolkaran, P., Joshi, A., Zheng, Y., Flaschel, M., De Lorenzis, L., Kumar, S., 2022\.  NN-EUCLID: Deep-learning hyperelasticity without stress data.  Journal of the Mechanics and Physics of Solids 169, 105076.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0022509622002538, doi:10.1016/j.jmps.2022.105076.
  * Tibshirani (1996) Tibshirani, R., 1996.  Regression Shrinkage and Selection via the Lasso.  Journal of the Royal Statistical Society: Series B (Methodological) 58, 267–288.  URL: http://doi.wiley.com/10.1111/j.2517-6161.1996.tb02080.x, doi:10.1111/j.2517-6161.1996.tb02080.x.
  * Touchette (2007) Touchette, H., 2007.  Legendre-Fenchel transforms in a nutshell .
  * Vlassis and Sun (2021) Vlassis, N.N., Sun, W., 2021\.  Sobolev training of thermodynamic-informed neural networks for interpretable elasto-plasticity models with level set hardening.  Computer Methods in Applied Mechanics and Engineering 377, 113695.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782521000311, doi:10.1016/j.cma.2021.113695.
  * Wang et al. (2021) Wang, Z., Estrada, J., Arruda, E., Garikipati, K., 2021\.  Inference of deformation mechanisms and constitutive response of soft material surrogates of biological tissue by full-field characterization and data-driven variational system identification.  Journal of the Mechanics and Physics of Solids 153, 104474.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0022509621001459, doi:10.1016/j.jmps.2021.104474.
  * Yu et al. (2021) Yu, H., Tian, X., E, W., Li, Q., 2021.  OnsagerNet: Learning stable and interpretable dynamics using a generalized Onsager principle.  Physical Review Fluids 6, 114402\.  URL: https://link.aps.org/doi/10.1103/PhysRevFluids.6.114402, doi:10.1103/PhysRevFluids.6.114402.
  * Zhang et al. (2022) Zhang, Z., Shin, Y., Karniadakis, G.E., 2022.  GFINNs: GENERIC Formalism Informed Neural Networks for Deterministic and Stochastic Dynamical Systems.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 380, 20210207\.  URL: http://arxiv.org/abs/2109.00092, doi:10.1098/rsta.2021.0207. arXiv:2109.00092 [cs, math].
  * Ziegler (1957) Ziegler, H., 1957.  Thermodynamik und rheologische Probleme.  Ingenieur-Archiv 25, 58–70.  URL: http://link.springer.com/10.1007/BF00536645, doi:10.1007/BF00536645.
  * Ziegler (1958) Ziegler, H., 1958.  An attempt to generalize Onsager’s principle, and its significance for rheological problems.  ZAMP Zeitschrift für Angewandte Mathematik und Physik 9, 748–763.  URL: http://link.springer.com/10.1007/BF02424793, doi:10.1007/BF02424793.
  * Ziegler (1972) Ziegler, H., 1972.  Systems with internal parameters obeying the orthogonality condition.  Zeitschrift für angewandte Mathematik und Physik ZAMP 23, 553–566.  URL: http://link.springer.com/10.1007/BF01593978, doi:10.1007/BF01593978.