Learning Physics-Consistent Material Behavior Without Prior Knowledge

inst1]Zhichao Han

[inst1]organization=Institute for Building Materials, ETH Zürich,
addressline=Laura-Hezner-Weg 7,

inst1]Mohit Pundir
inst2]Olga Finkfn1
inst1]David S. Kammercor1 fn1
[cor1]Corresponding author. E-mail: .
[fn1]Joint supervision.

[inst2]organization=Laboratory of Intelligent Maintenance and Operations Systems,
addressline=Station 18,

Accurately modeling the mechanical behavior of materials is crucial for numerous engineering applications. The quality of these models depends directly on the accuracy of the constitutive law that defines the stress-strain relation. Discovering these constitutive material laws remains a significant challenge, in particular when only material deformation data is available. To address this challenge, unsupervised machine learning methods have been proposed to learn the constitutive law from deformation data. However, existing approaches have several limitations: they either fail to ensure that the learned constitutive relations are consistent with physical principles, or they rely on a predefined library of constitutive relations or manually crafted input features. These dependencies require significant expertise and specialized domain knowledge. Here, we introduce a machine learning approach called uLED, which stands for unsupervised learning of hyperelastic constitutive relations in dynamic environments. uLED overcomes the limitations by using the input convex neural network (ICNN) as the surrogate constitutive model. Specifically, we improve the optimization strategy for training this special neural network. Our training approach allows ICNN to be trained end-to-end using direct strain invariants as input across various materials. Furthermore, we utilize the nodal force equilibrium condition at the internal domain as the training objective, which enables us to learn the constitutive relation solely from displacement recordings in dynamic settings. We validate the effectiveness of the proposed method on a diverse range of material laws. We demonstrate that it is robust to a significant level of noise and that it converges to the ground truth with increasing data resolution. We also show that the model can be effectively trained using a displacement field from a subdomain of the test specimen and that the learned constitutive relation from one material sample is transferable to other samples with different geometries. The developed methodology provides an effective tool for discovering constitutive relations. It is, due to its design based on dynamics, particularly suited for applications to strain-rate-dependent materials and situations where constitutive laws need to be inferred from in-situ measurements without access to global force data.

constitutive relation hyperelasticity dynamics unsupervised learning physics-consistent convergence

§ INTRODUCTION

Understanding the constitutive relation of materials is of fundamental importance for modeling and predicting their mechanical behavior. Accurately modeling the material behavior benefits numerous applications, including reliable prediction of landslides in nature [Cicoira et al., 2022], inverse design of mechanical materials in science and engineering [JH. Bastek, D. M. Kochmann, 2023], generating realistic scenarios in movies [Stomakhin et al., 2013], and many more. Since the constitutive relation cannot be measured directly in experiments, it is inferred indirectly through measurable and complementary quantities, such as macroscopic deformation and the corresponding mechanical response. To this end, conventional approaches typically assume a user-defined functional form (based on well-known material models) for the constitutive relation. Experimental data, such as displacement or strain data combined with the applied/reaction forces, are then used to infer the parameters within these models [Mihai et al., 2017]. Various techniques such as strain gauges [Harding and Welsh, 1983] and digital image correlation (DIC) [Rubino et al., 2019] are commonly used to collect such complementary data. Despite its usefulness, selecting an appropriate functional form for the constitutive relationship is challenging. It requires a high level of expertise and can often result in an oversimplified representation of the actual material behavior <cit.>.

The challenge to infer the constitutive model is further exacerbated in in-situ measurements where only one of the complementary data (either strain or boundary force but not both) is accessible. For example, measuring the biomechanical response (force measurement) in-situ of biological tissue to identify lesions and diseases is particularly difficult [YJ. Lim, D. Deo, T. P. Singh, D. B. Jones, S. De, 2009, Mazza et al., 2008].
In addition, understanding the material behavior of rate-dependent materials is crucial for applications such as simulating seismic site response [Miao et al., 2022] or analyzing the behavior of many polymer materials [Siviour and Jordan, 2016, J. E. Field, S. M. Walley, W. G. Proud, H. T. Goldrein, C. R. Siviour, 2004]. This requires consideration of dynamic deformation, where the constitutive relation depends on high strain rates. Consequently, determining constitutive relations without predefined functional forms and under various loading conditions (static or dynamic) remains a challenging but impactful task.

Machine learning (ML), particularly neural network-based approaches, offers a promising alternative for acquiring constitutive relations. These methods replace user-defined function forms with neural networks, reducing the reliance on domain expertise and minimizing user interference in the learning of constitutive laws. Neural network-based methods have been successfully applied to learn various constitutive relationships, such as elasticity [As'ad et al., 2022], viscoelasticity [Chen, 2021] and plasticity [Li et al., 2022]. However, traditional  ML methods typically depend on stress data for training, making them challenging to apply to experimental data where  stress is not directly measurable.
Instead, unsupervised methods have been developed to learn constitutive relations without using the direct stress data, such as the different variants of EUCLID [Flaschel et al., 2021, Flaschel et al., 2022, Joshi et al., 2022, Thakolkaran et al., 2022], and the work independently developed in [Wang et al., 2021].

In addition, we note an increasing interest in computer graphics for learning the constitutive relation to simulate or predict the motions of deformable objectives [Du, 2023]. An exemplary work is the NCLaw [Ma et al., 2023] that trains an MLP [Bishop, 2006] to approximate the constitutive relation. These unsupervised methods bypass the need for stress data by training the machine learning model to satisfy the linear momentum equation for a given displacement field.

Comparison between existing exemplary unsupervised and proposed ML approaches for learning material models considering three desired properties.

colspec=Q[m,c,6cm] Q[m,c,4cm]  Q[m,c,7cm] Q[m,c,3cm] ,row1=3.5ex, row2-5=5ex, row6=9ex, row7=5ex
Methods     No need for stress data for training      No need for constructing candidate constitutive relation library or hand-crafting input feature     Enforce physics-consistency

EUCLID variants [Flaschel et al., 2021, Flaschel et al., 2022, Joshi et al., 2022, Thakolkaran et al., 2022]

Yes     No     Yes

NCLaw [Ma et al., 2023]

Yes     Yes     No

Our approach     Yes     Yes     Yes

Despite their great potential, unsupervised ML approaches also have limitations, as summarized in Table <ref>. Some methods depend on prior knowledge of the constitutive relation to design a library of candidate constitutive relations [Flaschel et al., 2021, Flaschel et al., 2022, Joshi et al., 2022, Wang et al., 2021]. This reliance on a predefined candidate library restricts the ML model to approximating only a limited set of possible constitutive relations.  Consequently, if the actual material behavior deviates significantly from this pre-assumed set of possible constitutive relations, the model may fail to accurately capture the true constitutive relation.
NN-EUCLID [Thakolkaran et al., 2022] overcomes this limitation by using an input convex neural network (ICNN) [Amos et al., 2017], specifically  designed to approximate convex functions. This approach, previously employed in a supervised context by As'ad et al., 2022 serves as the surrogate constitutive model. However, NN-EUCLID requires additional feature transformation of the strain invariants, using hand-crafted features (e.g., mathematical terms specifically relevant for potential constitutive functions) rather than the direct strain invariants as  inputs. This is because ICNN requires some learnable parameters to be non-negative. This constraint makes it challenging to train ICNN end-to-end solely with strain invariants.
In contrast, NCLaw [Ma et al., 2023], which was developed primarily for computer graphics applications, bypass the need for preparing a candidate constitutive relation library or hand-crafting input features by using a standard MLP as the surrogate constitutive model. This  allows for direct training solely with strain invariants. However, this approach can result in learned constitutive relations that may look realistic but violate physical constraints. For instance, the learned energy density might not be convex with respect to strain invariants in hyperelasticity.

Further, while current methods circumvent the need for full-field stress data, they still rely on full-field displacement data along with information on boundaries (e.g., the reaction forces) to avoid learning trivial solutions for the constitutive relation (e.g., always outputting zero stress for any input strain in cases without external body force density) [Flaschel et al., 2021, Flaschel et al., 2022, Joshi et al., 2022, Wang et al., 2021, Thakolkaran et al., 2022]. These two complementary data sources are not readily available in some situations, such as in  in-situ experiments on living biological materials. To summarize, inferring a physics-consistent material model becomes challenging when there is no prior knowledge about the constitutive relation and  no access to complementary data such as stresses or boundary forces.

In this work, we aim to combine the advantages of both EUCLID and NCLaw (see the last entry in Table <ref>). We improve the training strategy for ICNN to satisfy the constraint that certain  parameters must be non-negative. Our proposed training strategy enables efficient  backpropagation of gradients under this constraint. Consequently, we can rely solely on strain invariants (without hand-crafting input features) to learn the constitutive relation for various materials. Moreover, we consider the nodal force equilibrium condition for the internal domain as the training objective. With this objective function, the proposed framework only requires using recorded temporal evolution of displacements to learn the constitutive relation, without any information about stress or boundary forces. Hence, we propose a new approach, named uLED, which leverages the governing equations in dynamics to learn material behavior. uLED is trained end-to-end, generalizes to different constitutive relations, and eliminates the reliance on prior knowledge of materials. Through extensive numerical simulations, we demonstrate the model's capability to uncover accurate strain-stress relations for various materials. Moreover, the learned constitutive relation converges to the ground-truth constitutive relation as we refine the data resolution used for training uLED.

In this work, we propose a new approach, named uLED, to learn material behavior by leveraging the governing equation in dynamics and hence eliminating the need to know the boundary forces. To achieve this property, uLED uses vanilla strain invariants as input to a neural network without hand-crafting the input feature, thus alleviating the reliance on prior knowledge of constitutive relations. Further, uLED learns the constitutive relation in a completely unsupervised way, i.e., only relying on the measured local displacements information without global properties, and guarantees the learned relationship is physics-consistent (e.g., convexity of strain energy). Through extensive numerical simulations, we demonstrate the model's capability to uncover the accurate strain-stress relation for various materials. We will show that the learned constitutive relation converges to the ground-truth constitutive relation when we refine the data resolution used for training uLED.

§ METHODOLOGY: UNSUPERVISED LEARNING OF HYPERELASTICITY IN DYNAMICS (ULED)

Illustration of the reference configuration $\domain$ and  the deformation over time. The displacement $\superscriptT{t}\vu$ of an arbitrary material point at each time $t$ is the difference between  its position $\superscriptT{t}\vx$ at time $t$ and its original position $\mX$ in the reference configuration $\domain$, i.e., $\superscriptT{t}\vu = \superscriptT{t}\vx - \mX$.

§.§ Considered mechanical system

We consider a deformable solid with its reference configuration $\domain$ in $\mathbb{R}^d$ of dimension $d$, undergoing dynamic deformation over time and occupying the domain $\superscriptT{t}\domain$ at time $t$. The solid is subject to applied displacements on the Dirichlet boundary $\partial \domain_{D}$ and applied tractions on the Neumann boundary $\partial \domain_{T}$.

An arbitrary material point located at position $\mX$ in the reference domain $\domain$ is positioned at $\superscriptT{t}\vx$ in the current domain $\superscriptT{t}\domain$. The displacement of this material point at time $t$ is given by $\superscriptT{t}\vu$, where $\superscriptT{t}\vu =  \superscriptT{t}\vx -  \mX$ (see Fig. <ref>). The deformation gradient at time $t$ is defined as $\superscriptT{t}\mF =  \mI + \nabla_\mX  \superscriptT{t}\vu \in \mathbb{R}^{d\times d}$, or, in index notation, $\superscriptT{t}F_{ij} = \partial \superscriptT{t}u_{i}/\partial X_j \equiv \superscriptT{t}u_{i,j}$ ($i,j = 1, \ldots, d$).

The stress in the solid is described by the first Piola-Kirchhoff stress tensor $\superscriptT{t}\mP = \cP(\superscriptT{t}\mF) \in \mathbb{R}^{d\times d}$, where $\cP$ represents the constitutive law relating the stress to the material's deformation. The solid must satisfy the conservation of linear momentum at time $t$, which, for the reference configuration, is given by:

\begin{equation}
\label{eq:governing_equ_strong_form}
\forall i = 1, \ldots, d: \quad \rho \superscriptT{t}{}{\ddot{u}{}_i} - \superscriptT{t}P_{ij, j} - \superscriptT{t}B_i = 0 ~ \text{ in } \Omega ~,
\end{equation}

where $\superscriptT{t}B_i$ is the $i$-th dimension of the body force density $\superscriptT{t}\mB \in \mathbb{R}^{d}$ at time $t$. In addition, the stresses must  balance with the applied tractions on $\partial \domain_{T}$.
Unlike previous methods [Thakolkaran et al., 2022, Flaschel et al., 2021, Wang et al., 2021] that consider static equilibrium and require knowledge of boundary forces to learn the constitutive relation, we consider the system under  dynamic deformation (Eq. <ref>). This approach  enables us to learn the constitutive relation without any information about boundary forces, as we will demonstrate.

§.§ Objective and proposed method

The objective of our proposed method is to infer the constitutive law that relates material stress to strain, i.e., $\cP(\mF)$, solely based on observations of the temporal evolution of the displacement field $\superscriptT{t}\vu$. Concretely, we aim to construct a neural-network-based surrogate model $\cPapprox(\mF; \theta)$ to approximate the constitutive material law, i.e., $\cPapprox(\mF; \theta) \approx \cP(\mF)$, where $\cPapprox(\mF;\theta)$ is parameterized by $\theta$. For a given experimentally measured displacement field $\superscriptT{t}\vugiven(\mX)$, we aim for $\cPapprox(\mF;\theta)$ to satisfy the conservation of linear momentum (Eq. <ref>) for all times $t \leq T$. We denote the predicted first Piola-Kirchhoff stress by $\superscriptT{t}\widehat{\mP}$, $\superscriptT{t}\widehat{\mP} = \cPapprox(\superscriptT{t}\mF; \theta)\in \mathbb{R}^{d\times d}$.  In contrast to prior work, we limit the spatial region of interest $\subdomain$ to a subdomain of $\domain$, i.e., $\subdomain \subseteq \domain$, which leads to the following strong form:

\begin{equation}
% \forall i = 1, \ldots, d: \quad \rho\superscriptT{t}{}{\ddot{u}_i} -  \superscriptT{t}\widehat{P}_{ij,j} - \superscriptT{t}B_i = 0~ \;  | \;\superscriptT{t}\vu(\mX) = \superscriptT{t}\vugiven(\mX),~ \forall \,\mX \in \subdomain  ~.
\forall i = 1, \ldots, d: \quad \rho\superscriptT{t}{}{\ddot{u}_i} -  \superscriptT{t}\widehat{P}_{ij,j} - \superscriptT{t}B_i = 0~ \text{ in } \subdomain  ~.
\label{eq:strong_opt_prob}
\end{equation}

Furthermore, since  we consider dynamic deformation where the term $\rho\superscriptT{t}{}{\ddot{u}_i}$ is nonzero, we do not need to know boundary forces to avoid learning the trivial solution[The trivial solution occurs  when the learned stress is a constant that only depends on the body force density, assuming the body force density is uniform throughout solid. In the case of static equilibrium (where $\rho\superscriptT{t}{}{\ddot{u}_i}$ is zero), boundary forces are  necessary to regularize the learned constitutive relation and  avoid such trivial solutions.].

To solve this problem, we convert it into its weak form , which can then be solved as an optimization problem. The weak form of Eq. <ref> is given by

\begin{equation}
\label{eq:governing_equ_weak_form}
\superscriptT{t}\functional_i(\theta; v_i) \equiv \int_{\subdomain}v_i \rho \superscriptT{t}{}{\ddot{{u}}_i} \dx + \int_{\subdomain}  v_{i,j}\superscriptT{t}\widehat{P}_{ij}  \dx - \int_{\subdomain}v_i \superscriptT{t}B_i \dx = 0~,
% \; | \;\superscriptT{t}\vu(\mX) = \superscriptT{t}\vugiven(\mX), ~\forall \,\mX \in \subdomain \quad
%\argmin_\theta  \int_{\domain}v_i \rho\ddot{{u}}_i(t) \dx + \int_{\domain}\sum_j \cPapprox_{i,j}(t; \theta)\nabla_j v_{i}  \dx - \int_{\domain}v_i B_i(t) \dx - \int_{\partial\domain_{T_i}}v_i T_i(t) \diff A = 0 \quad \forall \,\mX \in \subdomain
\end{equation}

where $v_i$ is an arbitrary test function in the functional space $V = \{v_i \mid v_i = 0 \text{ on } \subdomainboundary \}$. It is important to note that $\superscriptT{t}\functional_i(\theta; v_i)$ depends on the learnable parameters $\theta$ because $\superscriptT{t}\widehat{P}_{ij}$ is the output of the neural-network-based surrogated constitutive model parameterized by $\theta$. Given that the displacements are known everywhere, including on $\subdomainboundary$, we are free to choose this $V$. This choice ensures that the unknown tractions from the boundary do not affect the conservation of linear momentum.
Thus, the optimization problem is formulated as follows:

\begin{equation}
\label{eq:optimization_problem}
\text{Find $\theta$ such that } \superscriptT{t}\functional_i(\theta; v_i) = 0 \text{ for all } v_i \, .
\end{equation}

Finally, we note that the constitutive relation is learned to satisfy the principle of energy conservation through Eq. <ref> without requiring stress data. Specifically, we calculate the kinetic energy $\int_{\subdomain}v_i \rho \superscriptT{t}{}{\ddot{{u}}_i} \dx$ based on nodal accelerations and use it to determine the expected internal strain energy $\int_{\subdomain}  v_{i,j}\superscriptT{t}{P}_{ij}  \dx$ (assuming the body force is zero). Consequently, with the surrogate model $\cPapprox(\mF; \theta)$ satisfying Eq. <ref>, the predicted internal strain energy $\int_{\subdomain} v_{i,j}\superscriptT{t}\widehat{P}_{ij} \dx$  closely approximates its target value $\int_{\subdomain}  v_{i,j}\superscriptT{t}{P}_{ij}  \dx$.

The concept and architecture of the proposed uLED method. $\cPapprox(\mF; \theta)$ is the neural-network based surrogate model designed to approximate the constitutive material law, i.e., $\cPapprox(\mF; \theta) \approx \cP(\mF)$, where $\cPapprox(\mF;\theta)$ is parameterized by $\theta$. (a) The neural network architecture of the surrogate model $\cPapprox$. From the measured displacement field $\vu$, we compute the deformation gradient $\mF$, the Green-Lagrangian strain $\mE$, and strain invariants $I_1$ and $I_2$. The neural network (ICNN [Amos et al., 2017]) takes strain invariants as input and transforms them through  a series of latent representations $\vz_1, \vz_2, \ldots, \vz_5$ via multiple layers. The output at the last layer corresponds to the predicted energy density $\widehat{W}$, i.e., $\widehat{W}=\vz_5$. The predicted first Piola-Kirchhoff stress $\widehat{\mP}$ is computed by taking the gradient of the energy density. The learnable parameters $\theta$ in $\cPapprox(\mF; \theta)$ include the non-negative learnable layer-wise parameters
$\theta^{\text{(cvx)}}$ and normal learnable parameters $\theta^{\text{(fc)}}$ without the non-negative restriction. (b) Parameter optimization of the surrogate constitutive model $\cPapprox(\mF; \theta)$. The unique features (marked by the blue color) in uLED include the proposed training strategy which ensures the non-negative restriction of non-negative parameters $\theta^{\text{(cvx)}}$ during the training and the proposed loss function on nodal force equilibrium of the internal domain (Eq. <ref>).

To compute the integral numerically, we approximate the weak form using FEM (see  <ref> for the details). For each internal node $a$, i.e., not on Neumann or Dirichlet boundaries, the approximated weak form of governing equation is given by:

\begin{equation}
\label{eq:momentum_balance_every_node_internal}
\begin{split}
&\forall t, \ \forall i:\quad \quad \sum_{a=1}^{n_n} v_i^a \superscriptT{t}f_i^{a} = 0 \, , \\
&\superscriptT{t}f_i^{a} = \sum_b \left(\int_{\domain}N^a \rho N^b \dx \right) \ddot{u}^b_i(t) + \int_{\domain}\sum_j \nabla_j N^a \widehat{P}_{i,j}(t) \dx - \int_{\domain}N^a B_i(t) \dx
\end{split}
\end{equation}

where $N^a$ is the shape function associated at node $a$.
Since the test function $v_i^a$ can have arbitrary values for internal nodes, the $f_i^a$ in  <ref> must be zero for these internal nodes, i.e., $\forall a \text{ not on Neumann or Dirichlet boundaries}, \forall i\in d, \forall t:\, \superscriptT{t}f_i^{a} = 0 $. We use this property to train our neural network to approximate the unknown material model and minimize the absolute value of $\superscriptT{t}f_i^{a}$ for internal nodes at different time steps. This way, we circumvent the need for any information about the stress data or any information about the boundary force. We note that full-field observed displacement data is not necessary, i.e., we can select a part of the material to measure the dynamics and consider Eq. <ref> for this part. This property is very useful for real cases where we can only measure the dynamics of a part of a material.

§.§ Neural network design for $\cPapprox(\mF; \theta)$

We design uLED to learn the parameterized constitutive relation $\cPapprox(\mF; \theta)$  in three steps (see Fig. <ref>):

* Input processing: Transform the measured displacements into strain invariants

* Neural network module: Approximate the energy density function using the input convex neural network

* Output prediction: Extract the predicted energy density and compute the first Piola-Kirchhoff stress.

Input processing.
As mentioned in Sec. <ref>, the proposed method only relies on the recorded temporal evolution of nodal displacements at the given data resolution. The velocity and the acceleration of the displacements are computed from the displacements at consecutive time steps. Starting from the nodal displacements, we can compute the deformation gradient by its definition (see Sec. <ref>). While it is generally convenient to use the Green–Lagrange strain tensor $\mE = (\mF^T\mF - \mI)/2$ or equivalently the right Cauchy–Green deformation tensor $\mC = \mF^T\mF$ to  represent deformation without rigid body movement, we use strain invariants as the input for the neural network. In the 2D case, the strain invariants are:  $ I_1 = \tr(\mC), I_2 = \frac{1}{2} (\tr(\mC)^2 - \tr(\mC^2))$. The advantage of using strain invariants is that they ensure the derived constitutive relation possesses  desired properties [Satorras et al., 2021], such as  rotation-invariant predicted energy density   and  rotation-equivalent predicted stress  with the input deformation gradient.

Unlike previous approaches [Thakolkaran et al., 2022, As'ad et al., 2022], our inputs to the neural network are solely strain invariants, without any additional transformation or feature engineering. This greatly reduces the need for domain knowledge of different materials.

Neural network module.
We use the Input Convex Neural Network (ICNN) [Amos et al., 2017] to approximate the ground-truth energy density as a function of strain invariants. ICNN takes the strain invariants as input and outputs the energy density, i.e., $\widehat{W} = \text{ICNN}(I_1, I_2; \theta)$, where $\theta$ are the learnable parameters. The reason for choosing ICNN over other types  of neural networks is that ICNN guarantees the output is a convex function of the input. Convexity in the strain invariants ensures local material stability [Klein et al., 2022, Thakolkaran et al., 2022]. However, achieving  the desired convexity with ICNN comes with the restriction that part of its parameters ($\theta^{\text{(cvx)}}$ in the blue box of Fig. <ref>) must be non-negative. This restriction makes training ICNN non-trivial.
Specifically, we propose a new training strategy that is different from previous methods [Thakolkaran et al., 2022, As'ad et al., 2022] to enable the gradients efficiently backpropagated under the non-negative constraint. We introduce the proposed optimization strategy in Sec. <ref> and the detailed implementation of ICNN is explained in <ref>.

Stress prediction.
For hyperelastic materials, the first Piola-Kirchhoff stress $\mP$ can be computed as the gradient of the elastic energy density $W$, given by $\mP = \frac{\partial W}{\partial \mF}$, where $W = \mathcal{W}(\mF)$ depends on the deformation gradient $\mF$.  Equivalently, $\mP = \mF\mS$, where $\mS = \frac{\partial W}{\partial \mE}$ is the second Piola-Kirchhoff stress associated with the Green–Lagrange strain tensor $\mE$. The predicted energy density and the predicted first Piola-Kirchhoff stress should be zero in the absence of deformation, i.e., $W =0$ and $\mP = \mathbf{0}$ when $\mF = \mI$ or equivalently $\mE = \mathbf{0}$. We enforce this property, we use the predicted energy density and stress in the undeformed  state  as a reference to correct the predictions for any other state. Specifically, we denote the deformation gradient with no deformation by $\undeformedF$ and the Green strain with no deformation by $\undeformedE$ (i.e., $\undeformedF = \mI$, $\undeformedE = \mathbf{0}$), and the strain invariants of $\undeformedE$ by $I_1^{\undeformedSymbol}, I_2^{\undeformedSymbol}$  (i.e., $I_1^{\undeformedSymbol} = 0$ and $I_2^{\undeformedSymbol} = 0$). The predicted energy density for $\undeformedE$ is denoted by $\widehat{W}^{\undeformedSymbol}$, where $\widehat{W}^{\undeformedSymbol} = \text{ICNN}(I_1^{\undeformedSymbol}, I_2^{\undeformedSymbol}; \theta)$. The corresponding predicted first Piola-Kirchhoff stress $\widehat{\mP}^{\undeformedSymbol}$ for $\undeformedE$ is computed by $\widehat{\mP}^{\undeformedSymbol} = \partial \widehat{W}^{\undeformedSymbol} / \partial \undeformedF$. For any general state characterized by the principal invariants $I_1$ and $I_2$, the predicted energy density $\widehat{W}$ is computed by $\widehat{W} = \text{ICNN}(I_1, I_2; \theta) - \widehat{W}^{\undeformedSymbol}$. The corresponding predicted first Piola-Kirchhoff stress $\widehat{\mP}$ is computed by $\widehat{\mP} = \partial \widehat{W} / \partial \mF - \widehat{\mP}^{\undeformedSymbol}$. This ensures that the predicted energy density and the first Piola-Kirchhoff stress are zero in the absence of deformation.

§.§ Learning the optimal model

After defining the surrogate model for the constitutive relation in Sec. <ref>, we need to find the optimal values for parameters $\theta$ in $\cPapprox(\mF; \theta)$ to solve the optimization problem defined in Eq. <ref>. This problem presents two major challenges: 1) how to transform  the optimization problem in Eq. <ref> into a numerically computable objective function and 2) how to ensure  the non-negative requirements for parameters $\theta^{\text{(cvx)}}$. In the following, we outline our approach to address these challenges.

§.§.§ Training objective

To numerically solve the optimization problem defined in Eq. <ref>, we use the finite element method (FEM) to approximate $\functional(\theta; v_i)$ as defined in Eq. <ref> as $\functionalApprox(\theta; v_i)$. Suppose the region of interest $\subdomain$ is discretized by a set of elements $\mathcal{C}_e$ with $n_e = |\mathcal{C}_n|$ elements and a set of nodes $\mathcal{C}_n$ with $n_n = |\mathcal{C}_n|$ internal nodes. We denote the shape function associated with node $a\in \mathcal{C}_n$ by $N^a$ and the support elements on which $N^a$ is non-zero by $\text{supp}(N^a)$, i.e., $\text{supp}(N^a) = \{\text{el} \in \mathcal{C}_e :\mathrm{for~any~} \mX \in \text{el},  N^a(\mX) \neq 0\}$.
Thus, we get the approximated functional $\functionalApprox(\theta; v_i)$ to  approximate $\functional(\theta; v_i)$, i.e., $\superscriptT{t}\functionalApprox(\theta; v_i) \approx \superscriptT{t}\functional(\theta; v_i)$, as follows:

\begin{equation}
\label{eq:momentum_balance_every_node_internal}
\begin{split}
&\forall i=1, 2, \ldots, d:\quad \quad \superscriptT{t}\functionalApprox_i(\theta) \equiv \sum_{a=1}^{n_n} v_i^a \superscriptT{t}f_i^{a}(\theta)\, , \\
&\superscriptT{t}f_i^{a}(\theta) = \sum_{\text{el} \in \text{supp}(N^a)} \left[\int_{\text{el}}N^a \rho N^b \dx  \superscriptT{t}{}{\ddot{u}{}^b_i} + \int_{\text{el}} N^a_{,j}  \superscriptT{t}\widehat{P}_{ij} \dx - \int_{\text{el}}N^a \superscriptT{t}B_i \dx\right]
\end{split}
\end{equation}

where $\prescript{t}{}f_i^{a}(\theta)$ is the $i$-th component of the nodal force at node $a$ at time $t$, with $\widehat{P}_{ij}$ being  the components of the predicted stress from the surrogate constitutive model $\cPapprox(\mF, \theta)$. The computation of integrals in Eq. <ref> is discussed in <ref>.
Since the test function $v_i^a$ can have arbitrary values for internal nodes, $\superscriptT{t}f_i^{a}(\theta)$ in Eq. <ref> must be zero for these internal nodes at all times $t$, i.e., $\forall a \notin \partial\subdomain, \forall i\in d, \forall t:\, \superscriptT{t}f_i^{a} = 0 $. Therefore, we optimize the parameters $\theta$ to minimize the absolute value of $\superscriptT{t}f_i^{a}$ for internal nodes at all times $t$:

\begin{equation}
\label{eq:objective_function}
\argmin_{\theta = \{\theta^{\text{(cvx)}}, \theta^{\text{(fc)}}\}} \mathcal{L}(\theta) = \frac{1}{n_n}\frac{1}{d}\frac{1}{T}\sum_{a=1}^{n_n}\sum_{i=1}^d\sum_{t=1}^T |\superscriptT{t}f_i^{a}(\theta)| \; ,
\end{equation}

where $\superscriptT{t}f_i^{a}(\theta)$ is defined in Eq. <ref>. As indicated by the objective loss function (Eq. <ref>), the only required data for training the neural network are the nodal displacements $\{\superscriptT{t}\vu^{a}\}$ and accelerations $\{\superscriptT{t}\ddot{\vu}^{a}\}$ of the internal nodes $\{a \mid a \in \subdomain, a \notin \subdomainboundary\}$ at several time steps $t=1, 2, \ldots, T$. This approach circumvents the need for any information about the stress data or   boundary forces. Additionally, the region of interest $\subdomain$ can be a subset of $\domain$, indicating that full-field displacement observation is not necessary, and we can select a part of the material to infer the constitutive law. This property is particularly  useful in cases where we can only measure the dynamics of a part of a material.

§.§.§ Training strategy to satisfy the non-negative restriction

We use gradient-based methods to efficiently find the optimal values for $\theta$. However, compared to training a standard neural network, we have special constraints on some parameters (the non-negative parameters $\theta^{\text{(cvx)}}$) in ICNN. We introduce two additional steps specifically for training the ICNN to satisfy these constraints on the non-negative weights.
First, we convert any negative  weights in $\theta^{\text{(cvx)}}$ to their absolute values. This ensures  that $\theta^{\text{(cvx)}}$ is non-negative from the outset, even if they were initialized with random negative values.
Second, after each gradient descent step, we project the weights $\theta^{\text{(cvx)}}$ onto the non-negative space. Specifically, if a parameter value is negative, we set  it to 0; if the value is non-negative, we leave it unchanged.
These two steps ensure that the non-negative weights $\theta^{\text{(cvx)}}$ in ICNN, thereby guaranteeing that the predicted energy density is a convex function of the input. The pseudo-code for the training procedure is summarized in <ref>.

It is worth noting that our approach to enforcing the non-negative property differs from that of [Thakolkaran et al., 2022],  who apply a Softplus [Nair and Hinton, 2010] operation to the weights $\theta^{\text{(cvx)}}$ before multiplying the latent features.
We introduce the new approach because we observe that the previous approach often leads to the vanishing gradient problem, which makes it difficult to train the neural network.

We conjecture that this issue necessitates feature engineering in Thakolkaran et al., 2022, as the hand-crafted features help alleviate the training challenges.

§.§ Hyperparameter Settings

In our study, we primarily focus on the generalizability and transferability of the proposed uLED model, rather than accuracy. Consequently, we do not modify the architecture of the ICNN; instead, we adhere to the same architecture specified in the NN-EUCLID paper [Thakolkaran et al., 2022].  To ensure consistency across various experiments, we maintain consistent hyperparameter settings across different experiments. This includes using a uniform configuration for the ICNN in all experiments, including the number of layers, layer width, and the activation function. This standardized approach helps to assess the model's robustness across different testing scenarios.
Further details on the ICNN configuration are provided in Sec. <ref>.
For training  hyperparameters such as the learning rate, we adjust only the number  of training epochs to ensure convergence across  different experiments. Comprehensive  information about these training hyperparameters is available  in <ref>.

§.§ Key differences between proposed uLED and NN-EUCLID

The work of Thakolkaran et al., 2022 is the most relevant existing method to our proposed uLED approach. Both methods  use ICNN to guarantee local material stability and learn the strain-stress constitutive relation without stress data. However, the proposed uLED differs from [Thakolkaran et al., 2022] in three significant aspects.

* Improved Optimization Strategy: As mentioned in Sec. <ref>, we improve the optimization strategy in the training of  ICNN within uLED, enabling more effective backpropagation of  gradients.

* Input Data: We solely use strain invariants as input to ICNN, whereas  [Thakolkaran et al., 2022] requires hand-crafted transformations on the strain invariants based on the prior knowledge of the constitutive relation.

* Dynamic Formulation: We formulate the model in terms of dynamics and do not require  any information about the boundary forces. This setting makes our method applicable to in-situ measurements and scenarios where the constitutive relation depends on the strain rate.

§ VALIDATION AND EVALUATION OF THE PROPOSED FRAMEWORK

Illustration of the experimental configuration. (a) We simulate the transient motion of a 2D plate that is fixed at its left boundary and loaded with a temporally evolving traction on its right boundary. We record the displacements $\superscriptT{t}\vu^{a}$ of sampled material points (indicated by the yellow dots) from the FEM simulation. These sampled displacements $\superscriptT{t}\vu^{a}$ are used as the input for uLED. (b) The proposed uLED method is trained using   the sampled displacements $\superscriptT{t}\vu^{a}$. (c) The model predicts the energy density value and the first Piola-Kirchhoff stress tensor for an arbitrary element given the deformation gradient.

To showcase uLED's performance and versatility, we evaluate it with synthetic data from numerical simulations. We generate simulation results for a given material model, which constitutes the ground-truth constitutive relation. Then, we train uLED from the generated data and validate whether the learned constitutive relation approximates the ground-truth constitutive relations with sufficient accuracy.

To this end, we use FEM to generate a temporally evolving displacement field of a 2D plate as shown in Fig. <ref>. We fix the left boundary of the plate and apply transient tractions on the right boundary. We do not consider body forces in the simulation, i.e., $\superscriptT{t}\mB \equiv \mathbf{0}$.

Without loss of generality, we use dimensionless units for all simulations. Details of data generation are outlined in <ref>. The FEM simulation represents the underlying physics process, from which we sample a set of material points and record nodal displacements at different time steps. Nodal velocities and accelerations can be computed from the displacements. The sampled nodal displacements and accelerations are then used to train uLED.
For each constitutive relation, the sampled nodal displacements and accelerations are collected from 1,500 simulation time steps. To reduce training time, we construct a dataset by selecting every 7-th time step, resulting in 214 time steps. The first 80% of these 214 time steps are used for training the model, and the remaining 20% are used for validation. For each material, we conduct three independent experiments, randomly initializing the learnable parameters in uLED and training the model with the same training data. In each experiment, the validation data is used to select the best-trained model from the entire training period. Details to construct the uLED input from FEM simulation are presented in <ref>.

In real-world applications, the recorded data's resolution is determined by the measurement equipment, while the true physical process occurs typically at smaller scales. We reproduce this situation with our synthetic data by carefully choosing and distinguishing between the data resolution used during the generation and as the uLED input. During data generation using FEM simulations, we aim to accurately represent the relevant underlying physics and, hence, use a high resolution, referred to as $\frac{1}{h_{gen}}$ ($h_{gen}$ is the element size of the mesh). We use $h_{gen}=0.001$ consistently throughout this paper. To imitate the downsampling that generally occurs through the resolution of the measurement equipment, we use a lower resolution $\frac{1}{h}$ as input for uLED. Specifically, we set $\frac{1}{h} = \frac{1}{2 h_{gen}}$ throughout this paper, except in Sec. <ref>, where we explore the influence of data resolution on uLED.

Comparison between the predicted and ground-truth energy density and stress fields for a Neo-Hookean material. Data is shown for two arbitrarily chosen time steps (a) $t=7$, and (b) $t=1407$. On the left side, the predicted energy density ($\widehat{W}$) is compared to the ground-truth value $W$) via their normalized difference. On the right side, the magnitude of the first Piola-Kirchhoff stress fields ($|\widehat{\mP}|$) is compared to its ground-truth values ($|\mP|$). The predicted values shown are the mean of three independent experiments.

§.§ Showcase example: Learning the Neo-Hookean constitutive relation

In the first experiment, we use the Neo-Hookean material (see first entry in Table <ref> for its definition) as an example to thoroughly examine the effectiveness of uLED in learning the constitutive relation before we generalize the results to a set of other constitutive laws in subsequent sections. We directly use strain invariants $I_1$ and $I_2$ as the inputs to ICNN and evaluate the learned constitutive relation through two investigations. First, we randomly select two deformation states used in training and compare the predicted energy density and stress fields by the trained uLED model to the ground-truth energy density and stress fields computed with the Neo-Hookean constitutive relation. The difference between the predicted and ground-truth values is two or three magnitudes smaller than the scale of the ground-truth values (see Fig. <ref>). This indicates that the uLED model after training can effectively approximate the ground-truth constitutive relation. Additionally, it  demonstrates that when the energy density is  accurately learned, the stress is also  accurately captured. Therefore, we will focus exclusively  on the energy density in the subsequent evaluations.

Validation of the predicted energy density of the Neo-Hookean material. (a) The curves of the predicted and ground-truth energy density of the Neo-Hookean material for different deformation states (see Eq. <ref>) at a single material point. (b) The training curves from an arbitrary experiment of uLED for the Neo-Hookean material.

As second evaluation, similar to [Thakolkaran et al., 2022], we compare the predicted and ground-truth energy density values for multiple deformation states at a single material point. We construct three  types of deformation: uniaxial deformation (UD), biaxial deformation (BD), and shear deformation (SD), defined as follows:

\begin{equation}
\label{eq:F_path}
\mF^{\text{UD}} = \begin{bmatrix}
1 + \gamma & 0 \\
0 & 1 \\
\end{bmatrix},\,\,\,\,\,\,
\mF^{\text{BD}} = \begin{bmatrix}
1 + \gamma & 0 \\
0 & 1 + \gamma\\
\end{bmatrix},\,\,\,\,\,\,
\mF^{\text{SD}} = \begin{bmatrix}
1 & \gamma \\
0 & 1\\
\end{bmatrix}
\end{equation}

where $\gamma$ is varied to generate multiple deformation states. The range of  $\gamma$ values is determined from the training data. Specifically, we compute the maximum and minimum  values of $\gamma$ that appeared in the training data to evaluate the interpolation ability. Additionally, we extend the range by 20%  to generate deformation gradients that are not present in the training data, allowing us to evaluate the extrapolation ability of uLED. Results show that the curve of the learned energy density function (see Fig. <ref>a) closely matches the corresponding ground-truth curve closely in both  the interpolation and extrapolation regions.

The presented results have demonstrated  that the uLED is able to effectively approximate the constitutive relation for the Neo-Hookean material. However, it is essential to assess the efficiency of our proposed optimization strategy (as described in Sec. <ref>) in training uLED. We examine the training loss (Eq. <ref> on the training data) and the validation loss (Eq. <ref> on the validation data) throughout the training process for the Neo-Hookean material. Both the training loss and validation loss decrease significantly after about 420 epochs (see Fig. <ref>b), demonstrating that the proposed optimization strategy efficiently trains  uLED to learn the underlying material model. This  also indicates that the 1,000 training epochs used  are sufficient to achieve   a good accuracy. Furthermore, the validation loss closely follows the training loss throughout  the training process, indicating that there was no overfitting during  these 1,000 epochs.

§.§ Transferability to unseen geometries

Results for transferability on the unseen material sample (i.e. geometry). We train uLED on the plate sample and test the trained uLED on a new sample with holes. These two material samples used for training and testing are governed by the same Neo-Hookean law. We compare the predicted energy density field $\widehat{W}$ to the ground-truth $W$ on this new sample. The predicted values shown are the mean of three independent experiments.

One important characteristic of the learned material model is its transferability. This means that the uLED model trained on one sample of material should be applicable to other samples of the same material. For instance, if we have several samples of the same material, we can train uLED based on the observed dynamics of one sample. After this training, we expect that uLED can be directly applied to the other samples without additional training.

To demonstrate the transferability of uLED, we analyze the predicted energy density on an unseen sample with a different geometry from the sample used for training.  Specifically, we train uLED on a Neo-Hookean material sample with the plate geometry shown in Fig. <ref> (a). After training, we apply the trained model to predict the energy density of a new sample governed by the same Neo-Hookean law, using the plate geometry shown in Fig. <ref>. Our results show that the predicted energy density field on this new sample closely matches the ground-truth,  with differences that are two or three orders of magnitude smaller. This experiment indicates that the  constitutive relation learned by uLED is independent of the specific material sample used for training. Instead, uLED effectively learns the underlying constitutive law that governs the material behavior.

§.§ Generalizability to learn constitutive relation for various materials

Next, we investigate the generalizability of uLED across  different types  of materials. We consider various hyperelastic constitutive relations, as listed in <Ref>, which cover a wide range of material behavior features. Importantly, we maintain  the same ICNN architecture (including the number of layers, layer widths, activation function, etc.) for all considered material models. Additionally, we highlight that we consistently use strain invariants $I_1$ and $I_2$ as the inputs to ICNN without any further feature engineering.

Hyperelatic constitutive models in 2D. $\lambda$ and $\mu$ are Lamé's first and second parameter, respectively. $K$ is the bulk modulus. $I_1$ and $I_2$ are the strain invariants. $J = \det(\mF)$ is the determinate of the deformation gradient. $\bar{I}_1 = J^{-1} I_1$ and $\bar{I}_2 = J^{-2} I_2$ are the equivalent invariants of the unimodular $\bar{\mC} = (\det(\mC))^{-1/2}\mC$.

colspec=Q[m,c,6cm] Q[m,c,10cm] Q[m,c,3cm],row1=3.5ex, row2-3=4ex, row4=6ex, row5=4ex, row6=7.5ex, row7=4ex,
cell43 = font=1.5,
cell62 = font=1.5

Model     Energy density     Free parameters

Neo-Hookean (NH) [A. Logg, KA. Mardal, G. N. Wells (Eds.), 2012]     $W = \frac{\mu}{2} (I_1 - 2) - \mu \ln(J) + \frac{\lambda}{2} \ln(J)^2$     N/A

St. Venant–Kirchhoff (StVK) [A. Logg, KA. Mardal, G. N. Wells (Eds.), 2012]     $W = \frac{\lambda}{2} \tr(\mE)^2 + \mu \tr(\mE^2 )$     N/A

Two parameter Mooney-Rivlin (MR) [Bower, 2009]     $W = C_{10} (\bar{I}_1 - 2) + C_{01} (\bar{I}_2 - 1) + 0.5K(J - 1)^2$    $C_{10} = \frac{7}{16} \mu$
$C_{01} = \frac{1}{16}\mu$

Gent [Gent, 1996]     $W = - \frac{\mu}{2} (J_m \ln( 1 - ( \bar{I}_1 - 2) / J_m )  ) + 0.5 K (J - 1)^2$     $J_m = 10$

Arruda-Boyce (AB) [Bower, 2009]     $W = \mu ( \frac{1}{2} (\bar{I}_1 - 2) + \frac{1}{20N} ({\bar{I}{}}_1^2  - 4) + \frac{11}{1050 N^2} ({\bar{I}{}}_1^3  - 8) + \frac{19}{7000N^3} ({\bar{I}{}}_1^4  - 16) + \frac{519}{673750 N^4} ({\bar{I}{}}_1^5  - 32) ) + 0.5 K (J - 1)^2$     $N = 10$

Fung [YC. Fung, 2013]     $W = \frac{\mu}{2b} ( b (\bar{I}_1 - 2) + e^{b (\bar{I}_1 - 2)} - 1 ) + 0.5 K (J - 1)^2$     $b = 1$

We evaluate the performance of uLED for each constitutive relation by comparing the difference between the predicted and ground-truth energy density values for multiple deformation states at a single material point (as described in detail in Sec. <ref>). To measure uLED's performance, we use the normalized mean absolute error (NMAE) between the predicted energy density and the ground-truth energy density. For each deformation state, $\superscriptT{i}\mF$, we compute the predicted energy density, $\superscriptT{i}\widehat{W}$, and the ground-truth energy density, $\superscriptT{i}W$, corresponding to the deformation gradient (as given by Eq. <ref>). The NMAE is computed as the mean absolute error $\average{|\widehat{W} - W|}$, i.e., $\average{|\widehat{W} - W|} = \frac{1}{\Upsilon}\sum_{i=1}^{\Upsilon}|\superscriptT{i}\widehat{W} - \superscriptT{i}W|$, divided by the mean value of the ground-truth energy density $\average{|W|}$. Mathematically, this is expressed as: $\average{|W|} = \frac{1}{\Upsilon}\sum_{i=1}^{\Upsilon}|\superscriptT{i}W|$. The results indicate that the MAE between the predicted and ground-truth energy density values is two orders of magnitude smaller than the average ground-truth energy density values for all materials considered (see Fig. <ref>). However, uLED exhibits a relatively large error in shear deformation when learning the Fung constitutive law, suggesting that the Fung model is more challenging to learn compared to other constitutive relations. We hypothesize that this difficulty may require adjustments to  the neural network architecture (e.g., the number of layers, activation functions, etc.) to improve the performance in this case. Additional analyses of  the learned constitutive relations can be found in <ref>. These results demonstrate  that uLED is capable of effectively learning various constitutive relations without the need for manually constructing input features for each specific material. For example, using $(J-1)^2$ or $\bar{I}_1$ as input for learning the Mooney-Rivlin material is unnecessary since all cases considered cases were based simply on the strain invariants $I_1$ and $I_2$ as inputs for uLED.
Additionally, the decrease in both training loss and validation loss during the training process indicates  that the proposed optimization strategy is effective in training uLED to learn various constitutive relations  using only strain invariants as  input.
Detailed training curves illustrating this process are presented in <ref>.

Validation of the predicted energy density for various materials. We compute the NMAE of the predicted energy density compared to the ground-truth energy density for different deformation states at a single material point. We report the mean (indicated by the position of circles) and the standard derivation (indicated by the error bars) of the normalized MAE from three independent experiments.

§.§ Generalization on the unseen geometry

Results for validation on an unseen geometry. The comparison between the predicted energy density ($\widehat{W}$) / the magnitude of the first Piola-Kirchhoff stress fields ($|\widehat{\mP}|$) by uLED on an unseen geometry and the ground-truth values ($W, |\mP|$) for the Neo-Hookean material.  The predicted values shown are the mean of three independent experiments.

To demonstrate the generalization ability of the learned constitutive relation by uLED, we visualize the energy density/stress fields on an unseen geometry that is different from the geometry used to generate training data.  In detail, we train uLED on the plate geometry as shown in Fig. <ref>-left or Fig. <ref>. After training, we apply the trained model to the new geometry, as shown in Fig. <ref>. Our results show that the predicted fields match the ground-truth excellently, as the difference is several magnitudes smaller, indicating that the learned constitutive relation does not depend on a specific used for training the model.

§ TOWARDS APPLICABILITY TO EXPERIMENTAL DATA

Experimental data is typically less  precise and more noisy compared to the finite element method (FEM) data that was artificially generated and utilized in Sec. <ref> to evaluate the capabilities of uLED. Imperfections in training data can arise from several factors,  including coarse data resolution in the measurement of the displacement field, limited  observation windows,
the noise introduced by the digitalization of instruments,

and  dissipative processes in the material due to internal friction. This section examines how these factors impact uLED's capacity to discover constitutive relations.

§.§ Effect of data resolution

We begin  by investigating the effect of reduced data resolution on the performance of uLED. Specifically, we run a FEM simulation on a fine mesh ($h_{gen} = 0.001$) and train uLED using  different levels of coarsened data resolution $1/h$ where $h = 2h_{gen}, 3h_{gen}, \ldots, 10h_{gen}$. For this investigation, we use the Neo-Hookean constitutive law as an example,  though the results would be applicable to any of the other  constitutive laws considered (Table <ref>).

Convergence evaluation regarding data resolution with Neo-Hookean material.
(a) The difference between the predicted and the ground-truth energy density fields of a randomly selected state, normalized by the ground-truth value. The predicted energy density field is the mean of three independent experiments. (b) The mean absolute error (MAE) of the predicted energy density for different deformation states at a single material point, normalized by the average value of the ground-truth energy density. The mean (the solid line) and standard deviation (in shadow) of this normalized MAE metric are computed from three independent experiments.

First, we examine the discretization error in the training data for given data resolutions. For this purpose, we compute the average values of the kinetic and internal strain energy terms.
The kinetic energy refers to the average value of the inertia term in Eq. <ref> over the training data, i.e.,

\begin{equation}
\label{eq:resolution_convergence_kinetic_term}
\frac{1}{n_n}\frac{1}{d}\frac{1}{T} \sum_{a=1}^{n_n}\sum_{i=1}^d\sum_{t=1}^T \abs{\sum_{\text{el} \in \text{supp}(N^a)} \int_{\text{el}} N^a \rho N^b \dx\, \superscriptT{t}{\ddot{u}{}^b_i} }~.
\end{equation}

The internal strain energy refers to the average value of the strain energy term in Eq. <ref> over the training data, i.e.,

\begin{equation}
\label{eq:resolution_convergence_internal_energy_term}
\frac{1}{n_n}\frac{1}{d}\frac{1}{T}\sum_{a=1}^{n_n}\sum_{i=1}^d\sum_{t=1}^T \abs{\sum_{\text{el} \in \text{supp}(N^a)} \int_{\text{el}}  N_{,j}^a {P}_{ij} \dx}~.
\end{equation}

Notably, the internal strain energy is calculated with the ground-truth stress tensor instead of the predicted stress tensor. The discrepancy between these terms indicates the discretization error in the training data due to the coarsened data resolution $h$.
Without any discretization error, these two terms should be equal because of the conservation of linear momentum (see Eq. <ref>; supposing the body force density is zero). However, for training data with reduced resolution, there is, as one would expect, a discrepancy (see Fig. <ref> a), which shows that the information available in the data is (logically) of lower quality. This has direct consequences on the quality of the learned energy density function.

As expected, the data resolution has direct consequences on the quality of the learned energy density function. First, we analyze the predicted energy density fields by uLED trained with different data resolutions. To this end, we randomly select a deformation state used in training (the one shown  in Fig. <ref>) and apply  uLED,  trained with different data resolutions, to predict the energy density for this deformation state. The difference between the predicted and the ground-truth energy density fields decreases with higher data resolutions (see Fig. <ref> a). Second, similar to Sec. <ref>, we compare the energy density values predicted by the uLED model, trained  with different data resolutions, against  the ground-truth values for multiple deformation states at a single material point. As shown in Fig. <ref> b, the normalized mean absolute error of the predicted energy density consistently decreases with higher data resolution $1/h$ for all considered deformation states. These results illustrate that uLED learns the constitutive laws as well as possible given the quality of the available data – even relatively low-resolution data may still provide good predictions of the energy density. More importantly, our analysis demonstrates that the learned constitutive relation converges to the ground truth as we refine the data resolution used for training uLED. This suggests that uLED is potentially applicable to data with relatively low resolutions, as measured from physical experiments, and that the results can be improved by refining the data acquisition.

§.§ Learning constitutive relation with a limited observation window

An important challenge for inference methods is that experimental data, particularly for in-situ observations, is often limited to a part of the (test) specimen. The formulation of uLED offers  a significant  advantage in that it does not require a full displacement field to learn a constitutive law. To evaluate this capability of uLED, we consider a scenario  where the data is only available from a subdomain of the material. Specifically, we train uLED using the deformation data from  only a quarter of the full domain (marked by the dashed box in Fig. <ref>). For evaluation, we use the trained model to predict the energy density across  the entire domain (see Fig. <ref>). The results show that, compared to the reference case with training on the full domain,  the predicted energy density exhibits a larger error (see the normalized errors in Fig. <ref> and Fig. <ref>, noting that the ground-truth energy density fields are identical). This outcome is expected because the training data is more limited in this scenario, making the model trained with  partial observation more likely to make inaccurate  predictions for the deformation gradients not covered by the training data. We also note that the normalized error is large in areas where the ground-truth scale is small (the denominator in the normalization is small), similar to what was observed in the reference case with training on the full domain (see Fig. <ref>). This indicates that while uLED can function with partial observation data, the accuracy improves with more comprehensive training data.

Results for validation with partial observation on the Neo-Hookean material case. The uLED model is trained using the displacement fields from the area marked with the dashed box only (excluding  the border of the box). The model predicts the energy density for the full domain. The predicted values represent  the mean of three independent experiments. For comparison purposes, the selected deformation states are identical to those in the reference case with training on the full domain (Fig. <ref>)

§.§ Effect of noisy data

The quality of the learned model is directly related to the quality of the data. Since the inference of material constitutive laws  is aimed at  applications involving  experimentally acquired data, which is unavoidably noisy, it is crucial to evaluate the effect of such imperfections on the learning capabilities. One specific example of noisy experimental data originates from DIC measurements, which depend on the imaging device's pixel accuracy. In this study, we analyze the degree to which noise in the input measurements influences the performance of uLED, using the Neo-Hookean constitutive law as an example.

Convergence evaluation regarding noise level with Neo-Hookean material at data resolution $h_{gen}/h=1/2$.
(a) The difference between the predicted and the ground-truth energy density fields of a randomly selected state, normalized by the ground-truth value. The predicted energy density field is the mean of three independent experiments. (b) The mean absolute error (MAE) of the predicted energy density for different deformation states at a single material point, normalized by the average value of the ground-truth energy density. The mean and standard derivation of this normalized MAE are computed from three independent experiments.

We generate the noisy data as follows: First, we construct the clean data with a resolution of $1/h$. Similar to Thakolkaran et al., 2022, we add random noise from a zero-mean normal distribution with a standard deviation $\sigma$ to the clean displacements $\vu$ to obtain the noisy displacements $\widetilde{\vu}$. In real measurements,  noise is introduced to $\vu$ due to inaccurate measurements of $\vx$  (current position) and $\mX$ (reference position). Such measurement inaccuracies scale with the data resolution, meaning that  higher data resolution results in  smaller errors. Therefore, we also add noise to our clean data that scales according to the data resolution. For our analysis, we consider four different noise levels: $\sigma =10^{-8}h, 10^{-7}h, 10^{-6}h, 10^{-5}h$. Subsequently, from these noisy displacements, we compute the noisy velocities $\dot{\widetilde{\vu}}$ and noisy accelerations $\ddot{\widetilde{\vu}}$ (see <ref> for details). This approach provides us with noisy simulation data at the data resolution $1/h$.

We quantitatively compute the values of the inertia term and the internal energy term in the momentum conservation equation (Eq. <ref>) with different noisy levels. The kinetic energy refers to the average value of the inertia term in Eq. <ref> over the training data, i.e.,

\begin{equation}
\label{eq:noise_kinetic_term}
\frac{1}{n_n}\frac{1}{d}\frac{1}{T}\sum_{a=1}^{n_n}\sum_{i=1}^d\sum_{t=1}^T \abs{\sum_{\text{el} \in \text{supp}(N^a)} \int_{\text{el}} N^a \rho N^b \dx\,  \superscriptT{t}{\ddot{\widetilde{u}}{}^b_i} }~.
\end{equation}

The inertia term is computed with the noisy nodal accelerations $\ddot{\widetilde{\vu}}$. The internal strain energy refers to the average value of the strain energy term in Eq. <ref> over the training data, i.e.,

\begin{equation}
\label{eq:noise_internal_energy_term}
\frac{1}{n_n}\frac{1}{d}\frac{1}{T}\sum_{a=1}^{n_n}\sum_{i=1}^d\sum_{t=1}^T \abs{\sum_{\text{el} \in \text{supp}(N^a)} \int_{\text{el}}  N_{,j}^a \widetilde{P}_{ij} \dx}~.
\end{equation}

The stress tensor $\widetilde{\mP}$ in this term is computed by applying the ground-truth Neo-Hookean constitutive relation on the noisy deformation gradients. The discrepancy between these two terms indicates to which extent the recorded data violates the conservation of linear momentum due to the introduced noise.
These two terms should be equal (assuming the body force density is zero) without any discretization error and noise.

The introduced noise directly impacts the quality of the learned energy density function.  As shown in Fig. <ref> a, we analyze the energy density fields predicted by uLED,  trained with data affected by different noise levels $\sigma/h$ (here referred to as noisy-data-uLED). We randomly select a deformation state used in the training (the reference case in Fig. <ref>) and use the noisy-data-uLED model to predict the energy density for this deformation state. The difference between the predicted and the ground-truth energy density fields increases with higher noise levels ($\sigma/h$). Additionally, we compare the energy density values predicted by noisy-data-uLED and the ground-truth values for multiple deformation states at a single material point. As shown in Fig. <ref> b, the NMAE of the predicted energy density consistently decreases with decreasing noise levels $\sigma/h$ across all considered deformations. This behavior occurs because the introduced noise prevents the constitutive relation from satisfying Eq. <ref>. Consequently, the learned constitutive relation deviates from the ground truth when optimizing Eq. <ref>.

Therefore, we conclude that uLED's performance is quite robust to a certain level of noise. More importantly, the learned constitutive relation converges to the ground-truth constitutive relation as the noise in the measured displacements is reduced. The results suggest that uLED is potentially applicable to experimental data and that its accuracy may be further improved by increasing the precision of displacement measurements.

§.§ Learn constitutive relation in a dissipative setting

In this experiment, we evaluate whether uLED can be applied to cases where the deformation rate also affects material behavior. To achieve this, we consider the same dynamic setting described  in Sec. <ref>, but with the inclusion of energy dissipation. Consequently, the weak form of the governing equation now reads:

\begin{equation}
\label{eq:governing_equ_weak_form_with_dissipation}
\int_{\subdomain}v_i \rho \superscriptT{t}{}{\ddot{{u}}_i} \dx + \int_{\subdomain}  v_{i,j}\superscriptT{t}{P}_{(e)ij}  \dx - \underbrace{\int_{\subdomain}  v_{i,j}\superscriptT{t}{P}_{(d)ij}  \dx}_{\text{dissipative energy}} - \int_{\subdomain}v_i \superscriptT{t}B_i \dx = 0~.
\end{equation}

In Eq. <ref>, the internal strain energy $\int_{\subdomain}  v_{i,j}\superscriptT{t}{P}_{(e)ij}  \dx$ is determined by the elastic strain-stress relation $\cP_{e}$ with ${\mathbf{P}}_{(e)} = \cP_{e}(\mF)$. We use St. Venant-Kirchhoff constitutive law for $\cP_{e}(\mF)$.  The dissipative energy $\int_{\subdomain}  v_{i,j}\superscriptT{t}{P}_{(d)ij}  \dx$ is determined by the dissipative process $\cP_{d}$ with ${\mathbf{P}}_{(d)} =\cP_{d}(\dot{\mF})$.
In this experiment, we define the dissipation function as $\cP_{d}(\dot{\mF}) = \alpha_d\cdot g(\dot{\mF})$, where $\dot{\mF} = \odv*{(\mI + \nabla_\mX  \vu)}{t}= \nabla_\mX\dot{\vu} $ and $\alpha_d$ is the dissipation coefficient. We choose $g$ to be a linear function of $\dot{\mF}$.
Similar to previous experiments, we set the body force to zero.

We keep the ICNN module unchanged to learn the $\cP_{e}$ but add another neural network, a simple multilayer perceptron (MLP) [Bishop, 2006], to learn $\cP_{d}$. Importantly, we do not assume any form for $\cP_{d}$ when extending uLED to learn it. While the dissipation function is linearly related to the deformation gradient rate in the simulation, we do not provide this information to uLED. We only know that the dissipation energy depends on the deformation gradient rate $\dot{\mF}$ and  should be zero when the deformation gradient rate is zero. We enforce this constraint by using the same approach as in Sec. <ref>. We denote the predicted $\widehat{\mP}_d$ as $\widehat{\mP}_d^{\undeformedSymbol} = \text{MLP}(\dot{\mF}^{\undeformedSymbol}; \theta_{mlp})$ by the MLP when the deformation gradient rate is zero (i.e., $\dot{\mF}^{\undeformedSymbol} = \mathbf{0}$). The predicted $\widehat{\mP}_d$ for any given state with $\dot{\mF}$ is $\widehat{\mP}_d(\dot{\mF}) = \text{MLP}(\dot{\mF}; \theta_{mlp}) - \widehat{\mP}_d^{\undeformedSymbol}$.
The only change in the input data compared to previous experiments (without energy dissipation) is the inclusion of the deformation gradient rate $\dot{\mF}$.

To evaluate the learning performance, we compare the evolution of energies in this dissipative system, as shown in Fig. <ref>.
We observe that both the elastic energy and the accumulated dissipative energy predicted by uLED  closely align with the ground truth. This indicates that the elastic strain-stress relationship $\cP_e$ and the energy dissipation function $\cP_d$ are effectively learned. Additionally, we examine the curve of the predicted total energy, which is the sum of the elastic, dissipative, and kinetic energies, and find it to be in close agreement with the ground truth. These results suggest that uLED is potentially applicable for learning the constitutive relations in dissipative settings.

In detail, the elastic strain-stress relation $\cP_e$ determines the elastic potential energy at time $t_T$ by $E_{elastic}(t_T)= \int_{\domain} \superscriptT{t_T}W_e \dx$, where $\superscriptT{t_T}W_e$ is the elastic energy density at time $t_T$ determined by the deformation $\superscriptT{t_T}\mF$ and $\cP_e$. The damping term $\mP_d=\cP_{d}(\dot{\mE})$ contributes to the dissipated energy accumulated over time, which is referred to as $E_{damping}(t_T) = \int_{0}^{t_T}\int_{\domain} \frac{1}{2} \superscriptT{t}\mP_{d}: \nabla \superscriptT{t}\vu \dx \diff t $.
The kinetic energy in the system at time $t_T$ is computed by $E_{kinetic}(t_T) = \int_{\domain} \frac{1}{2}\rho \superscriptT{t_T}{\dot{\vu}}^2 \dx  $.
The total energy is the summation of the kinetic energy, the elastic potential energy and the dissipated energy, i.e., $E_{total} = E_{kinetic} + E_{elastic} + E_{damping}$.
The predicted elastic potential energy $\widehat{E}_{elastic}$ and the predicted dissipated energy $\widehat{E}_{damping}$ are computed by using the learned constitutive relation instead of the ground-truth relation. The predicted total energy $\widehat{E}_{total}$ are computed as $\widehat{E}_{total} = E_{kinetic} + \widehat{E}_{elastic} + \widehat{E}_{damping}$.
Fig. <ref> demonstrates that the curve of $E_{elastic}$ matches the curve of $\widehat{E}_{elastic}$ closely and the curve of $E_{damping}$ matches the curve of $\widehat{E}_{damping}$ closely, indicating the learned $\widehat{W}_e$ (or $\widehat{\mP}_{e}$) approximates the ground-truth $W_e$ (or $\mP_{e}$) well, and the learned $\widehat{\mP}_{d}$ approximates the ground-truth $\mP_{d}$ well. And the closeness of curves of ground-truth total energy $E_{total}$ and predicted total energy $\widehat{E}_{total}$  further verifies the good performance of uLED.

Results for validation in a dissipative setting.
The figure shows the evolution of energies in this dissipative system over time. The plotted predicted energies are the mean computed from three independent experiments. The standard derivation is vanishingly small.

§ DISCUSSION

The proposed uLED model is a powerful tool for uncovering the constitutive relation of materials.
Its key advantage lies in the ability to infer the underlying constitutive relation solely from recorded displacements over time, without requiring prior knowledge about the constitutive relation or measurements of the boundary forces. These features make uLED particularly suitable for two challenging scenarios: 1) learning the constitutive relation of novel and unconventional materials, and 2) learning the constitutive relation in in-situ situations where only displacements can be conveniently measured.

The advantage of uLED stems from its hybrid approach, which combines deep learning with physical laws. On the one hand, we leverage the universal approximation capability of neural networks to approximate constitutive relations without requiring any prior knowledge of the material model. On the other hand, by considering the force equilibrium for internal nodes as the training objective, we only need to record the displacements over time for training the neural network. Stress data and boundary forces are not included in the training objective, thus eliminating the  need for  this additional information.

The learned constitutive relation is ensured to be physics-consistent by three key design features. First, the hybrid approach inherently imposes physical constraints on the learned constitutive relation; the predicted stress naturally satisfies the conservation of linear momentum, which  is our training objective. Second, we exclusively  use strain invariants ($I_1, I_2$ for the 2D case) as  inputs to the neural network. This ensures that the learned energy density is rotation-invariant and that the learned stress is rotation-equivalent. Third, we employ the Input Convex Neural Network (ICNN) to approximate the underlying constitutive relation. With ICNN, we guarantee that the learned energy density is a convex function of strain invariants, thereby further ensuring local material stability.

While the ICNN ensures the desired local material stability, it requires  that some of its learnable parameters remain  non-negative. Training the ICNN efficiently while maintaining these non-negative parameters  poses a challenge. Previous  optimization strategies struggled to backpropagate gradients effectively under these constraints. To address this challenge, we developed  a new optimization strategy specifically for training the ICNN. Our results demonstrate the  effectiveness of this approach, successfully maintaining the non-negativity of parameters while ensuring efficient training.

When validating the proposed uLED method, we maintain a consistent  neural network architecture (e.g., the number of layers, the number of neurons in each layer, the activation function, etc.) across  all constitutive relations considered in the paper. Despite  this uniform  setup, we have demonstrated that uLED effectively learns different constitutive relations. However, we believe that performance can be further improved by tuning the neural network architecture for each specific constitutive relation and selecting the optimal architecture accordingly.

The main limitation of the proposed uLED method is that the constitutive relation is learned implicitly, meaning that a  a black-box neural network encodes the underlying constitutive relation, and we do not obtain an explicit form of this function. This issue can  potentially be addressed by applying symbolic regression (e.g., [Petersen et al., 2021]) to learn the explicit form of the material model. Specifically, for a given strain, the trained neural network can output the corresponding stress. Then, one can employ any  symbolic regression method to discover the explicit equation relating  the input strain to  the corresponding stress. Notably, the field of symbolic regression is rapidly advancing, and the  quality of the discovered symbolic equation will significantly depend  on the chosen symbolic regression algorithm.

Although this work presents the application of uLED in 2-dimensional space, the approach can easily be extended to 3-dimensional space. The only difference for 3-dimensional space is that the strain invariants will have three values instead of two. We can use these three strain invariants as inputs to the ICNN in uLED to learn the constitutive relation in 3-dimensional space. In the future, we aim to extend the framework to encompass a broader spectrum of constitutive relations of solids, such as viscoelasticity and plasticity, by incorporating strain history and hardening phenomena through tailored neural network architectures [Chen, 2021, Li et al., 2022].
While  we employ the finite element method to solve the optimization problem numerically (Eq. <ref>), uLED could also be developed using  other numerical methods. Since some materials, such as granular materials, are better described by numerical frameworks other than FEM (e.g., material point method [Yue et al., 2018, Daviet and Bertails-Descoubes, 2016]), we could consider the most suitable numerical method to solve the optimization problem in different scenarios. This flexibility in integrating different numerical methods enables uLED to be adapted for learning a wide range of materials.

§ CONCLUSIONS

In this research, we propose an unsupervised machine learning approach called uLED to learn the constitutive laws of hyperelastic materials. The key features of our approach are that it does not require any information about the stress labels or boundary forces as input. These advantages are achieved by leveraging dynamic, rather than static, observations of the displacement field. We have demonstrated that the proposed method can be applied to learn a diverse range of materials solely using  strain invariants, without the need for  hand-crafted input features for the neural network. This significantly  reduces the reliance on domain knowledge and expertise in potential constitutive relations. The constitutive relation learned by our approach is inherently  physics-consistent. We have also demonstrated that partial observations of the displacement field are sufficient to learn the constitutive relation.
Moreover, once uLED is trained based on the observed dynamics of one sample, the trained model can be directly applied to other samples of the same material with different geometries. This is because the constitutive relation learned by uLED is independent of the specific material sample used for training. uLED demonstrates excellent generalizability, effectively learning constitutive relations for various materials.

Furthermore, we have demonstrated that the learned constitutive relation is robust to significant levels of noise and converges to the ground truth as  data resolution increases. Finally, the proposed approach, which builds on observations of dynamic displacement fields, is particularly powerful for in-situ applications. This  includes  inferring the material model of biological tissues and strain-rate-dependent materials, providing an opportunity to discover accurate constitutive relations in complex real-life scenarios. In the future, we aim to extend uLED to real applications in 3-dimensional space.

§ ACKNOWLEDGEMENTS

This project has been funded by ETH Grant no. ETH-12 21-1.

§ CODE AVAILABILITY

The FEM simulation is implemented with FEniCS [Alnaes et al., 2015] (version 2019).
The implementation of the proposed uLED method is based on PyTorch [Paszke et al., 2019].
The source code is available on Gitlab: <https://gitlab.ethz.ch/cmbm-public/papers-supp-info/2024/learning-physics-consistent-material-behavior-without-prior-knowledge>.

[Cicoira et al., 2022]
authorA. Cicoira, authorL. Blatny,
authorX. Li, authorB. Trottet,
authorJ. Gaume,
titleTowards a predictive multi-phase model for alpine
mass movements and process cascades,
journalEngineering Geology volume310
(year2022) pages106866.

[JH. Bastek, D. M. Kochmann, 2023]
authorJH. Bastek, D. M. Kochmann,
titleInverse design of nonlinear mechanical metamaterials
via video denoising diffusion models,
journalNature Machine Intelligence volume5
(year2023) pages1466–1475.

[Stomakhin et al., 2013]
authorA. Stomakhin, authorC. Schroeder,
authorL. Chai, authorJ. Teran,
authorA. Selle,
titleA material point method for snow simulation,
journalACM Transactions on Graphics (TOG)
volume32 (year2013) pages1–10.

[Mihai et al., 2017]
authorL. A. Mihai, authorS. Budday,
authorG. A. Holzapfel, authorE. Kuhl,
authorA. Goriely,
titleA family of hyperelastic models for human brain
journalJournal of the Mechanics and Physics of Solids
volume106 (year2017) pages60–79.

[Harding and Welsh, 1983]
authorJ. Harding, authorL. M. Welsh,
titleA tensile testing technique for fibre-reinforced
composites at impact rates of strain,
journalJournal of materials science volume18
(year1983) pages1810–1826.

[Rubino et al., 2019]
authorV. Rubino, authorA. Rosakis,
authorN. Lapusta,
titleFull-field ultrahigh-speed quantification of dynamic
shear ruptures using digital image correlation,
journalExperimental Mechanics volume59
(year2019) pages551–582.

[Balmforth et al., 2014]
authorN. J. Balmforth, authorI. A. Frigaard,
authorG. Ovarlez,
titleYielding to stress: recent developments in
viscoplastic fluid mechanics,
journalAnnual review of fluid mechanics
volume46 (year2014) pages121–146.

[Kolymbas, 2012]
authorD. Kolymbas, titleConstitutive modelling of
granular materials, publisherSpringer Science & Business Media,

[YJ. Lim, D. Deo, T. P. Singh, D. B. Jones, S. De, 2009]
authorYJ. Lim, D. Deo, T. P. Singh, D. B. Jones, S. De,
titleIn situ measurement and modeling of biomechanical
response of human cadaveric soft tissues for physics-based surgical
journalSurgical endoscopy volume23
(year2009) pages1298–1307.

[Mazza et al., 2008]
authorE. Mazza, authorP. Grau,
authorM. Hollenstein, authorM. Bajka,
titleConstitutive modeling of human liver based on in vivo
in: booktitleMedical Image Computing and Computer-Assisted
Intervention–MICCAI 2008: 11th International Conference,
organizationSpringer, year2008, pp.

[Miao et al., 2022]
authorY. Miao, authorH. He, authorH. Liu,
authorS. Wang,
titleReproducing ground response using in-situ soil
dynamic parameters,
journalEarthquake Engineering & Structural Dynamics
volume51 (year2022) pages2449–2465.

[Siviour and Jordan, 2016]
authorC. R. Siviour, authorJ. L. Jordan,
titleHigh strain rate mechanics of polymers: a review,
journalJournal of Dynamic Behavior of Materials
volume2 (year2016) pages15–32.

[J. E. Field, S. M. Walley, W. G. Proud, H. T. Goldrein, C. R. Siviour, 2004]
authorJ. E. Field, S. M. Walley, W. G. Proud, H. T. Goldrein,
C. R. Siviour,
titleReview of experimental techniques for high rate
deformation and shock studies,
journalInternational journal of impact engineering
volume30 (year2004) pages725–775.

[As'ad et al., 2022]
authorF. As'ad, authorP. Avery,
authorC. Farhat,
titleA mechanics-informed artificial neural network
approach in data-driven constitutive modeling,
journalInternational Journal for Numerical Methods in
Engineering volume123 (year2022)

[Chen, 2021]
authorG. Chen,
titleRecurrent neural networks (rnns) learn the
constitutive law of viscoelasticity,
journalComputational Mechanics volume67
(year2021) pages1009–1019.

[Li et al., 2022]
authorX. Li, authorY. Cao, authorM. Li,
authorY. Yang, authorC. Schroeder,
authorC. Jiang,
titlePlasticitynet: Learning to simulate metal, sand, and
snow for optimization time integration,
journalProceedings of the 36th International Conference on
Neural Information Processing Systems volume35
(year2022) pages27783–27796.

[Flaschel et al., 2021]
authorM. Flaschel, authorS. Kumar,
authorL. De Lorenzis,
titleUnsupervised discovery of interpretable hyperelastic
constitutive laws,
journalComputer Methods in Applied Mechanics and
Engineering volume381 (year2021)

[Flaschel et al., 2022]
authorM. Flaschel, authorS. Kumar,
authorL. De Lorenzis,
titleDiscovering plasticity models without stress data,
journalnpj Computational Materials volume8
(year2022) pages91.

[Joshi et al., 2022]
authorA. Joshi, authorP. Thakolkaran,
authorY. Zheng, authorM. Escande,
authorM. Flaschel, authorL. De Lorenzis,
authorS. Kumar,
titleBayesian-EUCLID: Discovering hyperelastic material
laws with uncertainties,
journalComputer Methods in Applied Mechanics and
Engineering volume398 (year2022)

[Thakolkaran et al., 2022]
authorP. Thakolkaran, authorA. Joshi,
authorY. Zheng, authorM. Flaschel,
authorL. De Lorenzis, authorS. Kumar,
titleNN-EUCLID: Deep-learning hyperelasticity without
stress data,
journalJournal of the Mechanics and Physics of Solids
volume169 (year2022) pages105076.

[Wang et al., 2021]
authorZ. Wang, authorJ. B. Estrada,
authorE. M. Arruda, authorK. Garikipati,
titleInference of deformation mechanisms and constitutive
response of soft material surrogates of biological tissue by full-field
characterization and data-driven variational system identification,
journalJournal of the Mechanics and Physics of Solids
volume153 (year2021) pages104474.

[Du, 2023]
authorT. Du,
titleDeep learning for physics simulation,
in: booktitleACM SIGGRAPH 2023 Courses, SIGGRAPH '23,
publisherAssociation for Computing Machinery,
addressNew York, NY, USA, year2023.

[Ma et al., 2023]
authorP. Ma, authorP. Y. Chen,
authorB. Deng, authorJ. B. Tenenbaum,
authorT. Du, authorC. Gan,
authorW. Matusik,
titleLearning neural constitutive laws from motion
observations for generalizable PDE dynamics,
in: booktitleProceedings of the 40th International
Conference on Machine Learning, volume volume202,
publisherPMLR, year2023, pp.

[Bishop, 2006]
authorC. M. Bishop,
titlePattern recognition and machine learning,
journalSpringer google schola volume2
(year2006) pages5–43.

[Amos et al., 2017]
authorB. Amos, authorL. Xu, authorJ. Z.
titleInput convex neural networks,
in: booktitleProceedings of the 34th International
Conference on Machine Learning, organizationJMLR,
year2017, pp. pages146–155.

[Satorras et al., 2021]
authorV. G. Satorras, authorE. Hoogeboom,
authorM. Welling,
titleE(n) equivariant graph neural networks,
in: editorM. Meila, editorT. Zhang (Eds.),
booktitleProceedings of the 38th International Conference on
Machine Learning, volume volume139, publisherPMLR,
year2021, pp. pages9323–9332.

[Klein et al., 2022]
authorD. K. Klein, authorM. Fernández,
authorR. J. Martin, authorP. Neff,
authorO. Weeger,
titlePolyconvex anisotropic hyperelasticity with neural
journalJournal of the Mechanics and Physics of Solids
volume159 (year2022) pages104703.

[Nair and Hinton, 2010]
authorV. Nair, authorG. E. Hinton,
titleRectified linear units improve restricted boltzmann
in: booktitleProceedings of the 27th international
conference on machine learning, publisherOmnipress,
year2010, pp. pages807–814.

[A. Logg, KA. Mardal, G. N. Wells (Eds.), 2012]
authorA. Logg, KA. Mardal, G. N. Wells (Eds.),
titleAutomated solution of differential equations by the finite
element method: The FEniCS book, volume volume84,
publisherSpringer, year2012.

[Bower, 2009]
authorA. F. Bower, titleApplied mechanics of solids,
publisherCRC press, year2009.

[Gent, 1996]
authorA. N. Gent,
titleA new constitutive relation for rubber,
journalRubber chemistry and technology
volume69 (year1996) pages59–61.

[YC. Fung, 2013]
authorYC. Fung, titleBiomechanics: mechanical
properties of living tissues, publisherSpringer Science &
Business Media, year2013.

[Petersen et al., 2021]
authorB. K. Petersen, authorM. Landajuela,
authorT. N. Mundhenk, authorC. P. Santiago,
authorS. K. Kim, authorJ. T. Kim,
titleDeep symbolic regression: Recovering mathematical
expressions from data via risk-seeking policy gradients,
in: booktitleProceedings of the 9th International
Conference on Learning Representations, year2021.

[Yue et al., 2018]
authorY. Yue, authorB. Smith, authorP. Y.
Chen, authorM. Chantharayukhonthorn,
authorK. Kamrin, authorE. Grinspun,
titleHybrid grains: Adaptive coupling of discrete and
continuum simulations of granular media,
journalACM Transactions on Graphics (TOG)
volume37 (year2018) pages1–19.

[Daviet and Bertails-Descoubes, 2016]
authorG. Daviet, authorF. Bertails-Descoubes,
titleA semi-implicit material point method for the
continuum simulation of granular materials,
journalACM Transactions on Graphics (TOG)
volume35 (year2016) pages1–13.

[Alnaes et al., 2015]
authorM. S. Alnaes, authorJ. Blechta,
authorJ. Hake, authorA. Johansson,
authorB. Kehlet, authorA. Logg,
authorC. N. Richardson, authorJ. Ring,
authorM. E. Rognes, authorG. N. Wells,
titleThe FEniCS project version 1.5,
journalArchive of Numerical Software volume3

[Paszke et al., 2019]
authorA. Paszke, authorS. Gross,
authorF. Massa, authorA. Lerer,
authorJ. Bradbury, authorG. Chanan,
authorT. Killeen, authorZ. Lin,
authorN. Gimelshein, authorL. Antiga, et al.,
titlePytorch: An imperative style, high-performance deep
learning library,
in: booktitleProceedings of the 33rd International
Conference on Neural Information Processing Systems,
volume volume32, year2019, pp.

[Clevert et al., 2016]
authorD. Clevert, authorT. Unterthiner,
authorS. Hochreiter,
titleFast and accurate deep network learning by
exponential linear units (ELUs),
in: booktitleProceedings of the 4th International
Conference on Learning Representations, year2016.

[Glorot and Bengio, 2010]
authorX. Glorot, authorY. Bengio,
titleUnderstanding the difficulty of training deep
feedforward neural networks,
in: editorY. W. Teh, editorM. Titterington
(Eds.), booktitleProceedings of the 13th international conference
on artificial intelligence and statistics, volume volume9,
publisherPMLR, year2010, pp.

[Kingma and Ba, 2015]
authorD. P. Kingma, authorJ. Ba,
titleAdam: A method for stochastic optimization,
in: booktitleProceedings of the 3rd International
Conference on Learning Representations, year2015.

§ ICNN MODEL CONFIGURATION

To ensure that the learned energy density is a convex function of strain invariants, we use the ICNN [Amos et al., 2017] as the neural network model.
The layer-wise propagation rule of ICNN is:

\begin{equation*}
\vz_{i+1} = g_i(\vz_i W_i^{(z)} + \vx W_i^{(x)} + b_i), \quad \text{for } i = 0, 1, \ldots
\end{equation*}

where $\vz_{i+1}$ is the intermediate latent variable at layer $i+1$ ($\vz_0 \equiv 0$). $\vx$ is the input (strain invariants in this paper).  $W_i^{(z)}$ are the non-negative learnable weights at the $i$-th layer ($W_0^{(z)}\equiv 0$). $W_i^{(x)}$ and $b_i$ are the standard  learnable weights and biases without the non-negative restriction at the $i$-th layer. For convenience, we refer to the non-negative weights $W_i^{(z)}$ as $\theta_{i}^{\text{(cvx)}}$ and the standard weights $W_i^{(x)}$ and bias $b_i$ at the $i$-th layer as $\theta_{i}^{\text{(fc)}}$. The terms $\theta^{\text{(cvx)}}$ and $\theta^{\text{(fc)}}$ refer to all non-negative and standard  learnable parameters, respectively.

In this work, we use a 5-layer ICNN to learn the constitutive relation of all considered materials. The neural network architecture is shown in the following table:

Layer weights and bias configuration.

colspec=Q[m,c,3cm] Q[m,c,12cm],row1=4ex, row2-6=5ex
layer     parameter configuration

layer 0     $W_0^{(z)} \equiv 0$, $\quad W_0^{(x)}$ shape: [2, 128], $\quad b_0$ shape: [128]

layer 1     $W_1^{(z)}$ shape: [128, 128], $W_1^{(x)}$ shape: [2, 128], $b_1$ shape: [128]

layer 2     $W_2^{(z)}$ shape: [128, 128], $W_2^{(x)}$ shape: [2, 128], $b_2$ shape: [128]

layer 3     $W_3^{(z)}$ shape: [128, 128], $W_3^{(x)}$ shape: [2, 128], $b_3$ shape: [128]

layer 4     $W_4^{(z)}$ shape: [128, 1], $W_4^{(x)}$ shape: [2, 1], $b_4$ shape: [1]

Additionally, ICNN requires the non-linear activation functions $g_i$ to be convex and non-decreasing. In this work, we use the Exponential Linear Unit (ELU) [Clevert et al., 2016] as the activation function for all layers.

§ NUMERICAL COMPUTATION OF THE OPTIMIZATION PROBLEM (EQ. <REF>)

As the proposed method is developed within Finite Element Method (FEM) scheme, we start by giving a brief introduction to the essential background of FEM for continuum mechanics. We refer to <cit.> for details.

We represent vectors as column vectors. We use $\mX\in \mathbb{R}^d$ to denote the position of a material point in the reference configuration and $\superscriptT{t}\vx \in \mathbb{R}^d$ in the current configuration at time $t$, where $d$ ($d$=2 or 3) is the spatial dimension. The governing equation of linear momentum conservation is given by Eq. <ref>, i.e.,

\begin{equation}
\forall i = 1, \ldots, d: \quad \rho\superscriptT{t}{}{\ddot{u}_i} -  \superscriptT{t}\widehat{P}_{ij,j} - \superscriptT{t}B_i = 0~ \text{ in } \subdomain  ~.
\end{equation}

where $\superscriptT{t}\vu$ is the displacement field that is computed from the current position $\superscriptT{t}\vx$ of each material point to its original position $\mX$ in the reference domain, $\superscriptT{t}\vu = \superscriptT{t}\vx- \mX$. $\superscriptT{t}{\ddot{\vu}}$ is the acceleartion field.  $\superscriptT{t}\mF \in \mathbb{R}^{d\times d}$ is the deformation gradient at time $t$ and is defined as $\superscriptT{t}\mF =  \mI + \nabla_\mX  \superscriptT{t}\vu$, or, in index notation $\superscriptT{t}F_{ij} = \partial \superscriptT{t}u_{i}/\partial X_j \equiv \superscriptT{t}u_{i,j}$ ($i,j = 1, \ldots, d$). We denote by $\mP = \cP(\superscriptT{t}\mF) \in \mathbb{R}^{d\times d}$ the first Piola-Kirchhoff stress, and $\superscriptT{t}\mB \in \mathbb{R}^{d}$ the body force density.

The equivalent weak form of Eq. <ref> is:

\begin{equation}
% \label{eq:governing_equ_weak_form}
\int_{\subdomain}v_i \rho \superscriptT{t}{\ddot{{u}}_i} \dx + \int_{\subdomain} v_{i,j}\superscriptT{t}{P}_{ij} \dx = \int_{\subdomain}v_i \superscriptT{t}B_i \dx + \int_{\partial\domain_{T_i}}v_i \superscriptT{t}T_i \diff A
\end{equation}

where $v_i$ is the $i$-th dimension of an arbitrary test function $\vv$ in the functional space $V = \{v_i \mid v_i = 0 \text{ on } \partial \domain_{D_i}\}$ where $\partial \domain_{D_i}$ is the Dirichlet boundary at dimension $i$. $\domain$ is the domain of reference configuration. $\superscriptT{t}T_i$ is the traction force on the Neumann boundary at dimension $i$ denoted by $\partial \domain_{T_i}$.

To solve the optimization problem defined in Eq. <ref>, we use the FEM to approximate the $\superscriptT{t}\functional_i(\theta; v_i)$, as defined in Eq. <ref>. In this work, we use triangular elements to discretize the reference domain and use linear basis functions $\{N^a(\mX)\}$, where $a \in \mathcal{C}_n$ nodes (with $\mathcal{C}_n$ being  the set of internal nodes), to approximate the trial and the test function, as indicated by Eq. <ref>:

\begin{equation}
\label{eq:approximate_u_v}
\superscriptT{t}\vu(\mX) \approx \sum_{a=1}^{n_n}N^a(\mX)\superscriptT{t}\vu^{a}\, , \quad
\vv(\mX) \approx \sum_{a=1}^{n_n}N^a(\mX)\vv^a
\end{equation}

where $\superscriptT{t}\vu^{a}$ and $\vv^a$ are the values of the trial function and test function of node $a$ at time t.

By approximating the $\vu$ and $\vv$ in Eq. <ref> using Eq. <ref>, the balance equation for these internal nodes reduces to (Eq. <ref>)

\begin{equation}
\label{eq:appendix_momentum_balance_internal_node}
\begin{split}
&\forall t, \ \forall i:\quad \quad \sum_{a=1}^{n_n} v_i^a \superscriptT{t}f_i^{a}(\theta) = 0 \, , \\
&\superscriptT{t}f_i^{a}(\theta) = \sum_{\text{el} \in \text{supp}(N^a)} \left[\int_{\text{el}}N^a \rho N^b \dx \cdot \superscriptT{t}{\ddot{u}}_i^{b} + \int_{\text{el}}  N^a_{,j} \superscriptT{t}\widehat{P}_{ij} \dx - \int_{\text{el}}N^a \superscriptT{t}B_i \dx \right]
\end{split}
\end{equation}

where $\text{supp}(N^a)$ is the set of elements for which $N^a$ is non-zero.
Since the test function can have arbitrary values for internal nodes, the function $f_i^a(\theta)$ in Eq. <ref> must be zero for these internal nodes, i.e.,

\begin{equation}
\forall i\in d, \forall t:\quad \superscriptT{t}f_i^{a}(\theta) = 0
\end{equation}

Therefore, we derive the optimization objective as shown in Eq. <ref>. In the following, we explain how to numerically compute each term in $\superscriptT{t}f_i^{a}(\theta)$.
Since  we use linear basis functions in this paper, the term $\int_{\text{el}}N^a \rho N^b \dx$ can be computed analytically. Additionally, $N^a_{,j}$ and $\superscriptT{t}{\mF}$ are constant within each element because the basis functions are linear. Consequently, the $\superscriptT{t}\widehat{\mP} = \cPapprox(\superscriptT{t}\mF; \theta)$ is constant within each element. Therefore, the term $\int_{\text{el}} N^a_{,j} \superscriptT{t}\widehat{P}_{ij} \dx$ can be computed by multiplying  $ N^a_{,j}$, $\superscriptT{t}\widehat{P}_{ij}$, and the volume of the element. We do not consider body forces, so the term $\int_{\text{el}}N^a \superscriptT{t}B_i \dx$ is zero for each element.

§ DETAILS IN ULED TRAINING

We randomly initialize the parameters in uLED using the Xavier initialization algorithm [Glorot and Bengio, 2010]. We use the Adam optimizer [Kingma and Ba, 2015] with a learning rate of 5E-4 to optimize all learnable parameters. We split the input along the time dimension to mini-batches, with  the batch size  chosen to avoid exceeding the GPU memory capacity. The number of training epochs is set to 1,800 for the Fung material and 1,000 for all other materials to ensure the proposed ML model achieves good accuracy. (see Fig. <ref>). To ensure  $\theta^{\text{(cvx)}}$ remains non-negative, we project the negative values in $\theta^{\text{(cvx)}}$ to zero after each gradient descent step.
The pseudo-code for the training procedure is summarized in Algorithm <ref>.
The experiments are conducted on a server equipped with four RTX 3090 GPUs.

uLED training procedure

Prepare time-invariant inputs: $\rho, \{N^a\}_{a=1,\ldots,n_n}, \{N_{,j}^a\}_{a=1,\ldots,n_n}$
Prepare displacements at data resolution $1/h$: $\{\superscriptT{t}\vu^{a}\}_{a=1,\ldots,n_n}^{t=1,\ldots,T}$, $\{\superscriptT{t}\ddot{\vu}^{a}\}_{a=1,\ldots,n_n}^{t=1,\ldots,T}$
Randomly initialize $\theta^{\text{(cvx)}}$ as $\theta^{\text{(cvx)}, \text{now}}$,  $\theta^{\text{(fc)}}$ as $\theta^{\text{(fc)}, \text{now}}$
$\theta^{\text{(cvx)}, \text{now}} \gets$  ($\theta^{\text{(cvx)}, \text{now}}$)
converge or reach max training epochs

Split $\{\superscriptT{t}\vu^{a}\}_{a=1,\ldots,n_n}^{t=1,\ldots,T}$, $\{\superscriptT{t}\ddot{\vu}^{a}\}_{a=1,\ldots,n_n}^{t=1,\ldots,T}$ to mini-batches along time steps
*[l]If necessary, split along spatial domain as well

each batch
*[l]Compute the batch loss
Compute $\mathcal{L}(\theta^{\text{(cvx)}}, \theta^{\text{(fc)}})$ by Eq. <ref>

*[l]Optimize parameters, using Adam as the optimizer

$\theta^{\text{(fc)}, \text{new}} \gets$  Optimizer.step $\left(\nabla_{\theta^{\text{(fc)}}} \mathcal{L}(\theta^{\text{(cvx)}}, \theta^{\text{(fc)}}),  \theta^{\text{(fc)}, \text{now}} \right)$

$\theta^{\text{(cvx)}, \text{new}} \gets$  Optimizer.step $\left(\nabla_{\theta^{\text{(cvx)}}} \mathcal{L}(\theta^{\text{(cvx)}}, \theta^{\text{(fc)}}), \theta^{\text{(cvx)}, \text{now}} \right)$

$\theta^{\text{(cvx)}, \text{new}} \gets$  ($\theta^{\text{(cvx)}, \text{new}}$)

*[l]Update the current estimation of $\theta^{\text{(cvx)}}$ and $\theta^{\text{(fc)}}$

$\theta^{\text{(cvx)}, \text{now}} \gets \theta^{\text{(cvx)}, \text{new}}$
$\theta^{\text{(fc)}, \text{now}} \gets \theta^{\text{(fc)}, \text{new}}$

each element $\theta_{i,j}$ in $\theta$
$\theta_{i,j}$ < 0
$\theta_{i,j}$ = $-\theta_{i,j}$

each element $\theta_{i,j}$ in $\theta$
$\theta_{i,j}$ < 0
$\theta_{i,j}$ = 0

§ FEM SIMULATION SETTINGS

We use dimensionless units for all simulations, as the
machine-learning algorithm is not designed for any specific scale. For  the simulations, we set Young's module to 10,000, the Possion's ratio to 0.3, and the density $\rho$ to 1.0. Correspondingly, this results in $\lambda = 5769.2$ (Lamé's first parameter), $\mu = 3846.2$ (Lamé's second parameter), and $K = 8333.3$ (bulk modulus).

All simulations run 1,500 time steps with a time step size of 0.002.
We use the implicit Newmark method to update the dynamics for the FEM simulation. We use triangle mesh for the FEM simulation, with element size of $h_{gen} = 0.001$. The velocity and acceleration are updated as follows:

\begin{equation}
\label{eq:newmark_update}
\begin{split}
&\superscriptT{t+1}\ddot{\vu} = \frac{\superscriptT{t+1}\vu - \superscriptT{t}\vu - \Delta t \cdot \superscriptT{t}\dot{\vu}}{\beta \Delta t^2} - \frac{1 - 2\beta}{2 \beta}\superscriptT{t}\ddot{\vu}\\
& \superscriptT{t+1}\dot{\vu} =  \superscriptT{t}\dot{\vu} + \Delta t ((1-\gamma) \superscriptT{t}\ddot{\vu} + \gamma\cdot \superscriptT{t+1}\ddot{\vu})
\end{split}
\end{equation}

We denote  the displacements at the next time step $t+1$ as $\superscriptT{t+1}\vu$ and the displacements at the current time step $t$ as $\superscriptT{t}\vu$  (similarly  for velocities and accelerations). We set $\beta = \frac{1}{4}$ and $\gamma = \frac{1}{2}$ to ensure  the simulation unconditionally stable.

In Sec. <ref> where we generate noisy data, we impose noise
on the sampled nodal displacements at the resolution $1/h$ using  the following equation:

\begin{equation*}
\superscriptT{t}{\widetilde{u}}_{i}^{a} \leftarrow \superscriptT{t}{u}_{i}^{a} + \delta u
\end{equation*}

where $\superscriptT{t}{\widetilde{u}}_{i}^{a}$ is the $i$-th dimension of the noisy displacement of node $a$ at time $t$, and $\delta u \sim \mathcal{N}(0, \sigma^2)$ is a random number sampled independently from a zero-mean normal distribution with standard derivation $\sigma$. After obtaining  the noisy sampled nodal displacements $\{\superscriptT{t}\widetilde{\vu}\}_{t=1,\ldots,T}$, we generate the noisy velocities and accelerations by using the noisy displacements in the Newmark update scheme (Eq. <ref>).

§ DETAILS TO CONSTRUCT THE ULED INPUT FROM FEM SIMULATION

The FEM simulation generates nodal dynamics, including displacements and accelerations at data resolution of $1/h_{gen}$. To construct the input for training uLED at the coarser data resolution $1/h$ ($1/h < 1/h_{gen}$), we create a coarse mesh by selecting  every $k$-th node of the fine mesh (at resolution $1/h_{gen}$) along each spatial dimension, as illustrated  in Fig. <ref>. We vary  $k$ from $1$ to $10$ to create different levels of coarse-grained meshes. When $k$ is $1$, the new mesh is identical to the fine mesh, resulting in no discretization error. Conversely, when $k$ is $10$,  a significant  discretization error is introduced. The nodal displacements and accelerations in each coarse-grained mesh are directly interpolated from the recorded nodal displacements and accelerations in the fine mesh (see Fig. <ref>).
Throughout the main manuscript (except in Sec. <ref>, where we study the effect of different data resolutions), we set $h_{gen}/h  = 1/2$, i.e., $k=2$.

Illustration of the mesh coarse-graining and nodal dynamics interpolation. (a) The fine mesh with a mesh size of 0.05. (b) A coarse mesh created by taking every third node ($k=3$) along each dimension. We only simulate the nodal dynamics (displacements, velocities, accelerations) on the fine mesh. The nodal dynamics on the coarse-grained mesh are directly interpolated from the corresponding nodes on the fine data resolution. For example, the node marked by `a' with position (0.3, 0.15) on the coarse mesh has the same dynamics as the node marked by `$\alpha$' at the same position on the fine mesh. Note that this figure is for illustration purposes. We use a smaller mesh size ($h_{gen}=0.001$) for the fine data resolution in this paper.

When constructing the input for uLED, we sample every 7-th time step from the 1500 time steps for each simulation, resulting in  214 time steps for training and validation. Among these 214 time steps, we use the first 80% for training and the remaining 20%  for validation. The best-trained model is selected based on the objective function (mean absolute error of nodal force balance, Eq. <ref>) evaluated on the validation dataset. Notably, we do not use any information about the ground-truth stress labels.

§ EVALUATION METRICS

We consider the following criteria to evaluate the performance of the trained model:

1. Mean absolute error of momentum conservation (note the body force is zero)

\begin{equation*}
\frac{1}{n_n}\frac{1}{d}\frac{1}{T}\sum_{a=1}^{n_n}\sum_{i=1}^d\sum_{t=1}^T\abs{\sum_b \left(\int_{\domain}N^a \rho N^b \dx \right) \superscriptT{t}{\ddot{u}}^b_i + \int_{\domain}\sum_j  N^a_{,j} \superscriptT{t}{\widehat{P}}_{ij} \dx}
\end{equation*}

The momentum conservation error is used as the training objective. We also directly use loss to refer to this metric. Further, in Sec. <ref> and Sec. <ref>, we use the ground-truth ${\mP}$ instead of the predicted $\widehat{\mP}$ in the this metric to measure the mesh discretization error and the the noisy level.

2. Mean absolute error of the predicted energy density

\begin{equation*}
\frac{1}{n_n}\frac{1}{d}\frac{1}{T}\sum_{a=1}^{n_n}\sum_{i=1}^d\sum_{t=1}^T\abs{\superscriptT{t}{\widehat{W}}_{i}^{a} - \superscriptT{t}{W}_{i}^{a}}
\end{equation*}

This metric describes the average discrepancy between the predicted energy density and the ground-truth energy density at each dimension for all internal nodes computed from the deformation states at different time steps.

3. Mean absolute error of the predicted first Piola-Kirchhoff stress

\begin{equation*}
\frac{1}{n_n}\frac{1}{d}\frac{1}{T}\sum_{a=1}^{n_n}\sum_{i=1}^d\sum_{t=1}^T\abs{\superscriptT{t}{\widehat{P}}_{i}^{a} - \superscriptT{t}{P}_{i}^{a}}
\end{equation*}

This metric describes the average discrepancy between the predicted stress and the ground-truth stress at each dimension for all internal nodes computed from the deformation states at different time steps.

§ VISUALIZATION OF THE LEARNED ENERGY DENSITY  AND STRESS FUNCTIONS OF DIFFERENT MATERIALS

In Sec. <ref>, we  analyzed the curves of the predicted and ground-truth energy densities for the Neo-Hookean material across  different deformation states (see Eq. <ref>) at a single material point. Here, we extend this comparison to other materials (see Table <ref>), by comparing the curves of the predicted energy density and stress to the ground truth, as shown in Fig. <ref>. The results align with the analysis in Sec. <ref>, including  the relatively large error observed for the Fung material in Fig. <ref> .

Validation of the predicted energy density and stress of various materials.
(a) The curves of the predicted and  ground-truth energy density for different deformation states (see Eq. <ref>) at a single material point. (b) The curves of the predicted and ground-truth stress for different deformation states at a single material point. The interpolation region indicates the range of  $\gamma$  covered by the training data. The extrapolation region indicates the range of  $\gamma$  unseen in the training data. The predicted values represent the mean of three independent experiments. The standard deviation is not visible because it is too small.

§ TRAINING CURVES

To examine the efficiency of the optimization strategy for various material constitutive laws, we analyze the training loss and validation loss during the training process. For all considered materials, both the training loss and validation loss decrease significantly during the training process (see Fig. <ref>), demonstrating that the proposed optimization strategy efficiently trains uLED to learn the material models. Notably, for all material laws except St. Venant-Kirchhoff, the training curves exhibit a plateau where the loss remains almost constant (e.g., epochs $<300$ for the Mooney-Rivlin case) before a substantial decrease occurs within relatively few epochs. We believe this plateau arises from local minima in the optimization process,  influenced by the landscape of the constitutive relation functions. The St. Venant-Kirchhoff model does not display this plateau because its constitutive relation does not have a complex landscape, making it relatively straightforward to learn. Additionally, we observe that the Fung law requires more epochs to train uLED to achieve good accuracy. This indicates that the Fung constitutive relation is relatively difficult to learn, as evidenced by the relatively large error observed in Figure <ref>. Finally, we note that the validation and training loss are similar throughout the training process for all considered material laws, indicating that there is generally no overfitting in training uLED.

Training curves of uLED for different materials. The training loss and validation loss are computed by Eq. <ref> using the training data and validation data. The Fung material is trained with 1,800 epochs and all other materials are trained with 1,000 epochs. Each experiment of training uLED is repeated three times with randomly initialized parameters in ICNN – one arbitrary example is shown.

Comparison between the predicted and ground-truth energy density along different deformation states (Eq. <ref>). The interpolation region indicates the range of the $\gamma$ covered by training data. The extrapolation region indicates the range of the $\gamma$ that is unseen in training data. The predicted value is the mean of three independent experiments. The standard deviation is not visible because it is too small.

The predicted values and the ground-truth values of the magnitude of the first Piola-Kirchhoff stress along different deformation states (Eq. <ref>). The interpolation region indicates the range of the $\gamma$ covered by training data. The extrapolation region indicates the range of the $\gamma$ that is unseen in training data. The predicted values are the mean out of three independent experiments, while the standard derivation is too small to be visible.

§ VISUALIZATION OF THE PREDICTED ENERGY DENSITY/STRESS WITH DIFFERENT COARSE DATA RESOLUTIONS

The comparison between predicted energy density with different mesh sizes and the ground-truth energy density. The input state is randomly selected at a time step. The predicted energy density fields are the mean across three independent experiments.

The comparison between predicted first Piola-Kirchhoff stress with different mesh sizes and the ground-truth values. The input state is randomly selected at a time step. The plotted values are the mean magnitude of the first Piola-Kirchhoff stress across three independent experiments.

§ VISUALIZATION OF THE PREDICTED ENERGY DENSITY/STRESS WITH DIFFERENT NOISE LEVELS

The comparison between predicted energy density with different levels of noise and the ground-truth energy density. The input state is randomly selected at a time step. The predicted energy density fields are the mean across three independent experiments.

The comparison between predicted first Piola-Kirchhoff stress with different levels of noise and the ground-truth values. The input state is randomly selected at a time step. The plotted values are the mean magnitude of the first Piola-Kirchhoff stress across three independent experiments.