CANN-EUCLID: unsupervised constitutive artificial neural network model discovery from full-field data

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

CANN-EUCLID: unsupervised constitutive artificial neural network model discovery from full-field data

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

CANN–EUCLID: unsupervised constitutive artificial neural network model discovery from full-field data Benjamin Alheita,b , Siddhant Kumarb , Mathias Peirlincka a Department of BioMechanical Engineering, Faculty of Mechanical Engineering, Delft University of Technology,

arXiv:2606.14565v1 [cs.CE] 12 Jun 2026

b Department of Material Science Engineering, Faculty of Mechanical Engineering, Delft University of Technology,

Abstract Constitutive artificial neural networks (CANNs) provide an interpretable route to automated material model discovery, but they have so far been used predominantly in stress-supervised settings based on apparent stress–strain data from nominally homogeneous mechanical tests. Because each such test samples only a narrow loading path and provides homogenized rather than local stress information, robust constitutive discovery typically requires multiple complementary loading modes to constrain the multidimensional material response. This is particularly challenging for soft biological tissues, where repeated testing, damage, and sample-to-sample variability limit the amount of reliable mechanical information that can be extracted from a single specimen. Here, we combine CANNs with the stress-unsupervised full-field discovery framework EUCLID to identify sparse hyperelastic constitutive laws directly from displacement fields and reaction forces in a single heterogeneity-inducing loading case. The resulting CANN–EUCLID framework minimizes equilibrium imbalance while using sparsity-promoting regularization to select a compact set of active constitutive terms, without requiring local stress measurements or prescribing a specific constitutive law. We evaluate the approach on isotropic and anisotropic benchmark problems with prescribed ground-truth laws. When the ground-truth law is representable by the chosen CANN basis, the method recovers the correct active terms with near-exact accuracy, including exponential terms with parameters embedded inside nonlinear functions. When the exact law is not contained in the chosen CANN basis, the method retains shared terms and approximates missing contributions using available basis functions. Critically, generalization depends strongly on the deformation states sampled during discovery: exponential strain-stiffening terms can be recovered accurately when sufficiently probed, but can lead to large extrapolation errors when the stiffening regime lies outside the sampled deformation domain. Finally, through comparison of independent forward finite element validation simulations, we show that the behavior discovered by the CANN-EUCLID framework accurately replicates that of the ground truth. These results establish stress-unsupervised CANN discovery as a promising framework for interpretable, full-field constitutive model identification. Keywords: Constitutive artificial neural networks, EUCLID, Unsupervised constitutive discovery, Full-field inverse identification, Hyperelasticity, Sparse regression

1. Introduction Data-driven constitutive modeling has emerged as a powerful alternative to classical phenomenological constitutive model development in solid mechanics. Rather than prescribing a constitutive ansatz a priori and calibrating a corresponding set of parameters, machine-learning-based constitutive models aim to infer material behavior directly from mechanical testing data while preserving, to varying degrees, the physical structure required for robust prediction and finite element implementation [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. Within this broader landscape, hyperelasticity has proven to be a particularly fertile setting for constitutive machine learning because the existence of a strain-energy density provides a natural route to enforce thermodynamic consistency and to derive stresses by differentiation [6, 16, 17]. However, the central

challenge remains the same: the model must be expressive enough to capture complex nonlinear behavior, while remaining interpretable enough to support insight, trust, and generalization. Constitutive artificial neural networks (CANNs) have recently emerged as a prominent neural constitutive model class for hyperelasticity because they combine nonlinear expressivity with term-level interpretability [18, 19]. Like other neural constitutive models [20, 21, 22], CANNs can represent a broad class of nonlinear responses without requiring commitment to a predefined constitutive form. Unlike conventional black-box approaches, however, they decompose the constitutive response into interpretable basis terms, enabling direct insight into the contributions of invariants, stretches, and nonlinear functions [19]. This makes CANNs particularly attractive for automated model discovery, where predictive accuracy and mechanistic understanding are essential. Driven by these advantages, CANNs have seen rapid adoption in biomechanics and soft-tissue modeling, including applications to brain [23, 24, 25], skin [26], artificial and real meat [27], arteries [28, 29], and cardiac tissue [30, 31], and have been extended to path-dependent behavior such as viscoelasticity [32, 33], growth and remodeling [34], and Bayesian discovery [35]. Their growing practical relevance is further reflected by integration into commercial finite element workflows via universal Abaqus subroutines [24, 28, 36], positioning CANNs as a leading paradigm for interpretable constitutive machine learning. Despite this progress, CANNs have so far predominantly been trained in a stress-supervised setting, in which the network is calibrated using apparent stress–strain pairs obtained from mechanical tests such as uniaxial tension, biaxial tension, and triaxial shear [23, 26, 30]. In this learning paradigm, the experiments are either designed to approximate homogeneous deformation states or interpreted through homogenized force– displacement measures, yielding representative stress–strain data rather than directly measured local stress fields. While combining multiple loading paths improves identifiability and reduces the risk of overfitting to a single experiment, each test samples only a narrow subset of the admissible deformation space. As a result, stress-supervised calibration provides only sparse coverage of the multidimensional stress–strain response of the underlying material. This sparsity is problematic for constitutive discovery: a model may appear accurate along the measured loading paths while remaining weakly constrained elsewhere in the relevant deformation domain. This limitation is particularly pronounced for soft biological tissues, whose behavior is highly nonlinear, often anisotropic, and sensitive to the range and diversity of deformation states used during calibration [37, 38, 39, 40, 30]. In practice, obtaining dense local stress data over a sufficiently rich multidimensional deformation space is infeasible for two distinct reasons. Mechanically, even advanced multiaxial tests probe only restricted deformation manifolds, and stress fields are generally not directly measurable in heterogeneous specimens. Biologically, the missing deformation-space coverage cannot simply be recovered by repeatedly loading the same specimen, which may induce preconditioning, stress softening, or damage, nor by pooling nominally matched specimens, which is confounded by inter-sample and intra-sample microstructural variability [41, 42, 43]. In contrast, full-field deformation measurements are increasingly accessible through techniques such as digital image correlation and magnetic resonance or ultrasound-based strain imaging [44, 45, 46, 47], motivating an alternative approach in which constitutive laws are identified directly from full-field kinematics and global equilibrium rather than stress–strain pairs. In solid mechanics, this idea underlies a long line of stress-unsupervised inverse-identification methods based on full-field measurements. Classical examples include the virtual fields method (VFM), which identifies constitutive parameters through the principle of virtual work [48, 49, 44, 50], and the equilibrium gap method (EGM), which leverages violations of equilibrium as a basis for identification [51]. More recently, Efficient Unsupervised Constitutive Law Identification and Discovery (EUCLID) was introduced to identify sparse interpretable constitutive models by minimizing momentum imbalance while selecting a small number of active terms from an a priori library through sparsity-inducing Lp regularization [52]. However, its original formulation represents the constitutive response as a linear combination of library terms, which limits its ability to identify models in which material parameters enter nonlinear functions. This limitation is especially relevant for soft biological tissues, where exponential-type strain-energy functions are ubiquitous [53, 54, 55, 56, 42]. NN-EUCLID addresses part of this limitation by enabling unsupervised training of input-convex neural-network-based constitutive models [20]. This increases expressivity, but the resulting models do not retain the sparse term-level interpretability that makes CANNs attractive for constitutive 2

discovery. Thus, despite substantial advances in the EUCLID family of methods for hyperelasticity, plasticity [57, 58], viscoelasticity [59], generalized standard materials [60], and Bayesian identification [61], a stressunsupervised EUCLID framework for discovering sparse, interpretable constitutive models with nonlinear parameter dependence remains lacking. Here, we address this gap by combining CANNs [19] with the stress-unsupervised discovery framework of NN-EUCLID [20] and the sparsity-promoting Lp regularization strategy of EUCLID [52]. The resulting framework discovers CANN-based hyperelastic laws directly from displacement fields and reaction forces, without requiring local stress measurements or a predefined constitutive ansatz. It thereby brings together the best of all worlds: stress-unsupervised full-field training, sparse term-level interpretability, and nonlinearin-parameter expressivity. The remainder of this article is structured as follows. We present the overall framework for stress-unsupervised training of CANNs in Sec. 2. In Sec. 3, we assess the framework using synthetic full-field data generated from finite element simulations with prescribed ground-truth constitutive models for isotropic and anisotropic hyperelasticity. These benchmarks test whether the method can recover the correct constitutive terms when the ground-truth law is contained in the CANN basis, how it approximates ground-truth behavior when required terms are absent from the basis, and under which conditions the discovered models generalize to unseen deformation states. In Sec. 4, we discuss the implications of the results particularly in the context of soft biological tissues. Finally, we provide concluding remarks in Sec. 5. 2. Methodology 2.1. CANNs The architecture of CANNs is presented in detail in its originating work [18] with architectural variations proposed in following works [19, 29, 30, 62, 28, 24, 25, 23, 26, 27, 35, 63]. Here, without loss of generality, we provide a brief description of a particular class of CANN used in this work. We note that the contributions in this work are not limited to this specific architecture and can be similarly generalized to other CANN architectures. Details regarding intrinsic enforcement of objectivity, material symmetry, thermodynamic consistency, polyconvexity, etc. are available in the literature [19, 6, 16, 21, 20]. We note that the architecture of CANNs is designed so that these principles are guaranteed in the discovered behavior [19]. 2.1.1. Kinematics

We define the motion of a body φ as a unique map from points X ∈ Ω0 ⊂ 3 to x ∈ Ω ⊂ 3 at some time t, i.e. x = φ (X, t), where Ω0 is the initial undeformed domain and Ω is the current deformed domain. We define the deformation gradient F and right Cauchy-Green tensor C, respectively, as follows: F :=

∂φ , ∂X

C : = FTF .

(1)

To satisfy objectivity and material symmetry in the case of hyperelasticity, a strain energy density function is typically defined in terms of invariants of C. A typical choice of invariants in the case of anisotropy [56] is i 1h 2 I1 = tr (C) , I2 = tr (C) − tr C 2 , I3 = det (C) , I4,ij = Ai · CAj , I5,ij = Ai · C 2 Aj , (2) 2

where Ai , i = 1, ..., ns is the ith unit structural vector in the reference configuration, which represents the direction of a microstructural anisotropic constituent, e.g. the axis of a collagen fiber family or the normal of packed myofiber sheets1 . To capture the significantly different isochoric (volume-preserving) and volumetric 1 When there is only one structural vector, n

s = 1, we write I4 as opposed to I4,ij for the sake of brevity.

responses commonly observed in hyperelastic materials, we multiplicatively decompose the deformation gradient into an isochoric F̃ and volumetric part F̄ , i.e., F = F̃ F̄ ,

F̃ := J −1/3 F ,

F̄ := J 1/3 I ,

J = det (F ) ,

(3)

where J is the (volumetric) Jacobian and ˜ • and ¯• denote the isochoric and volumetric parts of •, respectively. The isochoric right Cauchy-Green tensor is then defined as (4)

C̃ := F̃ T F̃ , and the associated isochoric invariants are given by 2 1 ˜ ˜ tr C̃ − tr C̃ 2 I1 = tr C̃ , I2 = , 2

I˜4,ij = Ai · C̃Aj ,

I˜5,ij = Ai · C̃ 2 Aj .

(5)

2.1.2. Pseudo-invariants and their requirements Here, we briefly discuss the inputs to CANNs as this aids in their definition that follows in Sec. 2.1.3. The (isochoric) deformation invariants from Eq. (5) are not directly input into a CANN, but rather, a set of nK pseudo-invariants Ki , i = 1, . . . , nK . The specific form of these pseudo-invariants must be chosen carefully to satisfy the following requirements: (i) Zero strain energy density at F = I. CANNs are designed such that they evaluate to 0 when all of the inputs are 0. Consequently, to obtain a material model that results in Ψ (I) = 0, it is required that Ki F =I = 0, i = 1, . . . , nK . (ii) Zero stress at F = I. The first Piola-Kirchhoff stress is given by n

P =

K ∂Ψ ∂Ki ∂Ψ X = , ∂F ∂Ki ∂F i

(6)

which, in the special case that a given Ki is defined solely in terms of isochoric invariant I˜i , expands to P =

nK X ∂Ψ ∂Ki ∂ I˜i . ∂Ki ∂ I˜i ∂F

(7)

Ii ∂Ψ ∂Ki To ensure that P (I) = 0, it is required that at least one of the factors ∂K , ∂ I˜ , or ∂∂F evaluate to 0 when i i

∂Ψ F = I. In general, CANNs do not ensure that ∂K = 0, hence when using isochoric invariants either i F =I ˜ ˜ ∂Ki ∂ I˜i , or must be zero at F = I. For I˜1 and I˜2 this is trivial as ∂ I1 = ∂ I2 = 0. However, this ∂ I˜i

∂F

∂F F =I

∂J is not the case for I˜4,ij and I˜5,ij . Similarly, for pseudoinvariants Ki that reflect J, we have ∂F ̸= 0. F =I ∂Ki In these cases, the expression for the pseudo-invariants Ki must be chosen such that ∂ I˜ F =I = 0. i

(iii) Polyconvexity. Polyconvexity of the strain energy function [64] is the de facto sufficient condition used to guarantee the existence of a solution to a finite strain elasticity boundary value problem (BVP). Hence, the CANN should ultimately result in a polyconvex strain energy function. As is detailed in the following subsection (Sec. 2.1.3), the CANN architecture generates a sum of interpretable terms by composing various activation functions that are applied to the pseudo-invariant inputs. To ensure that the overall discovered behavior is polyconvex, one must ensure that each resulting term in the sum is polyconvex. As such, the pseudo-invariants must be chosen appropriately. We show that the choices for the activation functions and pseudo-invariants used in this work result in a polyconvex strain energy density function in Appendix A. 4

In this work, without loss of generality, we choose the following pseudo-invariants as inputs to the CANN: D E2  ˜ 3/2 I4,ij − 1 if i = j , 2 3/2 ˜ ˜ K1 = I1 − 3 , K2 = I2 −3 , K3 = [J − 1] , K4,ij = (8) I˜2 otherwise . 4,ij

However, other choices are also possible, see e.g. [25, 65]. Moreover, for the sake of polyconvexity, we use D E2 h i2 3/2 I˜2 as opposed to I˜2 [66] and I˜4,ij − 1 as opposed to I˜4,ij − 1 [67] (see Appendix A). 2.1.3. CANNs map invariants to strain energy density CANNs provide a structured neural architecture for representing the strain-energy density as a sum of interpretable terms. The computational graph for achieving this is illustrated in Fig. 1. In contrast to the densely connected layers typical of many neural network architectures, CANNs utilize a sparse, branching tree-like architecture. An additional contrast to conventional neural network architectures is that CANNs layer 1

layer 2

layer 3

layer 4

layer 5

Figure 1: Constitutive Artificial Neural Network (CANN) architecture. CANNs define a map from the deformation gradient F and structural vectors Ai to the strain energy density Ψ. This is done by defining a branching tree which automatically generates a large sum of interpretable terms – i.e. it results in nested for-loops in code. First, a set of nK scalar invariants is obtained as defined in Eq. (8). These are then operated on by a set of functions fj , j = 1, . . . , nf . The outputs of these functions are then multiplied by weights ϕijk and fed into a following set of functions gk , k = 1, . . . , ng . Finally, the outputs of functions gk k = 1, . . . , ng are multiplied by a second set of weights θijk and summed to provide the final strain energy density.

typically have a fixed depth of five layers, with each layer corresponding to a level of the branching tree. The first layer of the CANN takes the deformation gradient F and structural vectors Ai , i = 1, . . . , ns as input to the network. A set of pseudo-invariants, as discussed in Sec. 2.1.2, is calculated in the second layer, resulting in nK nodes. In the third layer, a set of functions f1 (◦) , . . . , fnf (◦) are applied to each pseudo-invariant from the second layer, leading to nK × nf nodes. The output of each node in the third layer is replicated ng times, each replicated output is multiplied by a unique weight ϕijk (i ∈ [1, nK ], 5

j ∈ [1, nf ], j ∈ [1, ng ]) and fed into a corresponding function in the fourth layer g1 (◦) , . . . , gng (◦), leading to nK × nf × ng nodes. The output of each node in the fourth layer is multiplied by a unique weight θijk and summed to provide the strain energy density as output in the fifth layer. The indices of the weights ϕijk and θijk associated with selected edges have been indicated in Fig. 12 . The mathematical description of the process illustrated visually in Fig. 1 and textually in the previous paragraph is concisely represented by the following equation: n

nK X f g X X , Ψ F , A1 , . . . , Ans := θijk gk ϕijk fj Ki F , A1 , . . . , Ans {z } | i j k

(9)

Ψijk

which clearly yields a sum of nK × nf × ng terms, or “sub-strain energy density functions” Ψijk . To create a concretized instance of a CANN from this abstract description, one must choose the set of pseudoinvariants to use Ki , i = 1, . . . , nK , and functions f1 (◦) , . . . , fnf (◦) and g1 (◦) , . . . , gng (◦). In this work, we choose to use the pseudo-invariants listed in Eq. (8). The typical choice for the functions f1 (◦) , . . . , fnf (◦) [36], which is also used in this work, is fj (◦) = ◦j . (10) For the functions g1 (◦) , . . . , gng (◦), we choose ( ◦ if k = 1 , gk (◦) = exp (◦) − 1 if k = 2 .

(11)

In Appendix A, we show that these activation-function choices and pseudo-invariants yield a polyconvex strain-energy density, provided that the CANN weights are non-negative. Accordingly, throughout training we constrain all trainable weights θijk and ϕijk to satisfy θijk , ϕijk ≥ 0. 2.2. CANNs meet unsupervised discovery Minimizing imbalance of equilibrium Following prior EUCLID works [52, 20], as a starting point, we assume that one has a dataset of displacement values measured at a finite number of material points nn that correspond with nodes in a FE mesh, and a finite number of times steps nt , i.e. U = uI,t ∈ 3 | I ∈ [1, nn ] , t ∈ [1, nt ] , (12)

and, if the material is anisotropic, we also assume that one has access to an approximation of the structural vector fields which we denote collectively as A. Additionally, we assume that one has access to the aggregated reaction forces on nβ boundaries at those time steps, i.e. R = Rβ,t ∈ 3 | β ∈ [1, nβ ] , t ∈ [1, nt ] . (13)

The goal is to discover the constitutive model and parameters that satisfy the weak form of equilibrium balance, i.e. to find θ = {(ϕijk , θijk ) | i ∈ [1, nK ], j ∈ [1, nf ], j ∈ [1, ng ]} that satisfies Z Z P (F , A; θ) : ∇v dΩ0 − t · v dΓN = 0 , ∀v , (14) Ω0

ΓN

2 When referring to specific components of ϕ ijk and θijk , we use commas between indices to disambiguate cases where one or more indices are larger than nine. For example, when writing ϕ1213 the values of specific indices are unclear (i = 1, j = 21, k = 3, or i = 12, j = 1, k = 3, etc.), so instead we write, e.g., ϕ1,21,3 which clarifies that i = 1, j = 21, and k = 3.

where v is an arbitrary test function that is from the same function space as u and we have assumed that body forces are negligible. To obtain a continuous field for u from the discrete measurements in U, we use basis functions N I associated with each node I that define a piecewise polynomial function space over the domain, i.e. nn X t u ≈ N I uI,t , (15) I

where •t denotes the value of a field at time t. This allows for obtaining an approximation for the deformation gradient over the domain as follows F t = I + ∇ut ≈ I +

nn X

uI,t ⊗ ∇N I .

(16)

Since u and v are from the same function space, we choose the same basis functions to provide an approximation for the test function v, i.e. nn X N I vI . (17) v≈ I

Substituting Eqs. (15) and (17) into Eq. (14) results in the following discretized version of the weak form of momentum balance: Z Z nn X I t I t I vi P F , A; θ ij ∇Nj dΩ0 − ti N dΓN = 0 , ∀viI . (18) I

Ω0

ΓN

Due to the arbitrariness of viI , Eq. (18) holds iff the expression in parentheses is zero for all (I, i). Accordingly, we interpret this as an unbalanced force at each node riI,t , i.e. Z Z I,t t I tti N I dΓN = 0 , (19) ri := P F , A; θ ij ∇Nj dΩ0 − ΓN Ω0 {z } | | {z } I,t rint,i

I,t rext,i

I,t I,t where we further identify internal and external unbalanced forces denoted by rint,i and rext,i , respectively, for node I in direction i at time t. Eq. (19) provides a basis for an objective function that can be minimized with respect to material parameters to approximately satisfy balance of equilibrium, namely X h I,t i2 L= ri . (20) I,i,t

However, in practice, one does not typically know the traction t on the constrained surfaces and so riI,t cannot be determined for the nodes on those surfaces. To circumvent this issue we note that although one does not have access to the traction, one does have access to the reaction force for the constrained surfaces, i.e. Eq. (13). Additionally, since the internal and external forces should balance pointwise, their sums should also balance over a given boundary, i.e. X I,t X I,t Riβ,t = rext,i = rint,i , (21) I∈D β

I∈D β

where D denotes the set of nodes on boundary β. Accordingly, instead of minimizing Eq. (20), we additively decompose the objective function into a component for the free and fixed nodes, respectively, 1 X X h I,t i2 rint,i , (22) Lint = nt nn t (I,i)∈D free  2 X X 1 I,t  Rβ,t − Lext = rint,i , (23) nt nβ β β

β,t

(I,i)∈D

where Dfree denotes all nodes that are not constrained. Sparsity promoting regularization The goal of model discovery is not only to identify a model that minimizes the imbalance of equilibrium, but also to select, among admissible fits, a compact constitutive representation with a small number of interpretable and physically meaningful terms. This requires balancing accuracy against sparsity, and thus predictability against interpretability. The de facto approach for achieving this in the context of constitutive modeling is to penalize the number of non-zero terms in a discovered model by adding the Lp -norm of the model coefficients to the loss equation, i.e. for a model of the form h (x) =

np X

(24)

αi hi (x) ,

the Lp -norm of its coefficients is given by ||α||p =

" np X

|αi |

#1/p

(25)

where np is the number of parameters. In particular this penalizes the number of non-zero coefficients for 0 ≤ p ≤ 1. Using such a regularization term to promote sparsity was first suggested in [68] where the authors generalized the use of L2 (Lp with p = 2) regularization originally proposed in [69] not for promoting sparsity but for stabilizing regression problems with a high degree of collinearity. This was first done in the context of constitutive model discovery in the original EUCLID work [52] and has later been adopted extensively in the training of CANNs [70, 62, 30, 23]. We note that the stress given by a CANN (the strain energy density for which is given in Eq. (9)) is given by nf ng nK X X X ∂Ψ X ∂Ψi,k,l ′ ′ ∂Ki P = θi,k,l ϕi,k,l f2,l f1,k = = , (26) ∂F ∂F | {z } | {z ∂F} i i,k,l

≡αi in Eq. (24)

≡hi in Eq. (24)

where • denotes the derivative of a function and we identify that θi,k,l ϕi,k,l is equivalent to αi in Eq. (24) ′ ′ ∂Ki and f2,l f1,k ∂F is equivalent to hi in Eq. (24). As such, we propose that the appropriate regularization term to use is nf ng nK X X X 1 p Lp = |θi,k,l ϕi,k,l | . (27) nK nf ng i ′

Hence, the final loss equation is given by

L = λint Lint + λext Lext + λp Lp

(28)

where Lint penalizes the imbalance of equilibrium on unconstrained nodes (Eq. (22)), Lext penalizes the imbalance of reaction forces on constrained nodes (Eq. (23)), Lp penalizes the number of non-zero terms in the final model, and λint , λext , and λp weight the relative contributions of Lint , Lext , and Lp , respectively, to the overall loss. With the loss Eq. (28) in hand, the overall CANN-EUCLID training procedure is summarized as follows (see Fig. 2): the structural vector and deformation gradient fields (Fig. 2 (a)), which are assumed to be measured data, act as inputs to a CANN (Fig. 2 (b)). The CANN maps these fields to a strain energy density, from which a trial stress field is obtained (Fig. 2 (c)) by differentiating the strain energy density with respect to the input F . The imbalance of equilibrium for the trial stress field is calculated and added to the loss (Fig. 2 (d)), in addition to the imbalance of the reaction forces (Fig. 2 (e)) and the sparsity-promoting Lp -norm of the model parameters (Fig. 2 (f)). If the loss is not below some tolerance ϵ, then the weights 8

(e) reaction forces

(a)

structural vector fields

(b) CANN during training

(f) (c) trial stress field

(d) internal force balance

deformation gradient field

external force balance

sparsity reg.

(g) update weights no (h) interpretable discovered behaviour sparsifiedCANN

relative contributions of terms

yes

Figure 2: Schematic summary of the CANN-EUCLID framework. (a) the structural vector and deformation gradient fields, which are assumed to be measured data, act as inputs to a CANN (b). The CANN maps these fields to a strain energy density, from which a trial stress field is obtained (c) by differentiating the strain energy density with respect to the input F . The imbalance of equilibrium for the trial stress field is calculated and added to the loss (d), in addition to the imbalance of the reaction forces (e) and the sparsity-promoting Lp -norm of the model parameters (f). If the loss is not below some tolerance ϵ, then the weights are updated using a gradient-based optimization algorithm (g). This process repeats until the loss is below ϵ, at which point we conclude that a sparse CANN has been identified that satisfies the internal balance of equilibrium and reaction forces within an acceptable tolerance (h).

are updated using a gradient-based optimization algorithm (Fig. 2 (g)). This process repeats until the loss is below ϵ, at which point we conclude that a sparse CANN has been identified that satisfies the internal balance of equilibrium and reaction forces within an acceptable tolerance (Fig. 2 (h)). In practice, direct optimization of Eq. (28) with sparsity-promoting regularization can be sensitive to initialization and can bias the final parameter values away from the primary equilibrium objective. We therefore use a three-stage training procedure consisting of pre-regularization, sparsity-driven subset selection, and post-regularization polishing. The full optimization protocol, together with all training and data-generation hyperparameters, is reported in Appendix B. 2.3. Evaluation of the CANN-EUCLID framework 2.3.1. Synthetic data generation We verify the success of the unsupervised training of CANNs using canonical benchmarks for unsupervised constitutive behavior discovery, as established in [52, 20, 21], for both isotropic and anisotropic hyperelasticity. To this end, we generate synthetic displacement field and reaction force data by applying biaxial loading to a plate with a central hole, as illustrated in Fig. 3 (a). The amount of deformation is parameterized by δ, which is linearly ramped from zero to one over ten time steps (see Appendix B for details). These data are obtained from FE simulations [11] with chosen ground truth constitutive models. In the isotropic cases

Figure 3: Schematic illustration of geometries used in training and validation. The training geometry (a) consists of a plate with a hole in the bottom left corner. It is subjected to vertical displacement on the top face, horizontal displacement on the right face, and roller boundary conditions on the left and bottom faces. The amount of deformation is parameterized by δ which is linearly ramped from zero to one over ten time steps. The validation geometry (b) consists of a plate with elliptical holes. It is subject to a unit of vertical displacement on the top face while the bottom face remains fixed. All other boundaries are traction free.

the ground truth constitutive models considered 3 are as follows: NeoHookean [71]: h i 2 Ψ = 0.5 I˜1 − 3 + 1.5 [J − 1] ,

Demiray [53]: i i h h 2 Ψ = 0.1 exp 5 I˜1 − 3 − 1 + 1.5 [J − 1] ,

Isihara [72]: h i h i h i2 2 Ψ = 0.5 I˜1 − 3 + I˜2 − 3 + I˜1 − 3 + 1.5 [J − 1] ,

Arruda-Boyce [73]: p p sinh βc 2 Ψ = 2.5 Nc βc λc − Nc log − cAB + 1.5 [J − 1] , βc Gent-Thomas [74]: h i 2 Ψ = 0.5 I˜1 − 3 + log I˜2 /3 + 1.5 [J − 1] , Ogden [75]: 2 1.3 1.3 Ψ = λ1.3 1 + λ2 + λ3 − 3 + 1.5 [J − 1] .

(29a)

(29b)

(29c)

(29d)

(29e) (29f)

3 We replace the original material parameters by arbitrary benchmark values and write the models directly to avoid introducing notation used only once.

√ Here, λc = I˜1 /3, βc = L−1 λc / Nc , and L−1 denotes the inverse Langevin function. The constants are set to Nc = 28, corresponding to the number of polymer chain segments, and cAB ≈ 3.7910, where the latter offsets the energy density to zero at F = I, since the Arruda–Boyce formulation does not vanish in the undeformed configuration. In the anisotropic cases, the considered ground truth behaviors are as follows: Anisotropic Neohookean (AN) [76]: h i D E2 2 Ψ = I˜1 − 3 + I˜4 − 1 + 1.5 [J − 1] ,

Holzapfel–Gasser–Ogden (HGO) [37]: D h i E2 2 Ψ = I˜1 − 3 + 0.25 exp 2 I˜4 − 1 − 1 + 1.5 [J − 1] , Meaney [77]: h i h 2 i −1 2 Ψ = 0.5 I˜1 − 3 + 0.5 I˜4 + 2I˜4 − 3 + 1.5 [J − 1] , Merodio–Ogden (MO)[78]: h i h i2 2 Ψ = 0.5 I˜1 − 3 + 0.5 I˜5 − 1 + 1.5 [J − 1] ,

Humphrey–Yin (HY)[54]: " # q 2 ! i i h h 2 Ψ = 0.1 exp 5 I˜1 − 3 − 1 + exp 8 I˜4 − 1 − 1 + 2.5 [J − 1] ,

Gasser–Ogden–Holzapfel (GOH)[38]:   *  +2  h i 1 ˜1 ˜4 I I 2 + − 1  − 1 + 2.5 [J − 1] . Ψ = 0.5 I˜1 − 3 + exp 3 6 6 2

(30a)

(30b)

(30c)

(30d)

(30e)

(30f)

Moreover, the structural vector is taken to be homogeneously aligned in the positive y-direction, i.e. A = T [0 1 0] . This benchmark set is chosen deliberately to include both ground-truth laws that are exactly representable by the chosen CANN basis and laws containing invariants or functional forms outside that basis. This enables separate assessment of exact term recovery and out-of-basis approximation. For further details on synthetic data generation, we refer the interested reader to Tab. B.1 in Appendix B. 2.3.2. Error metric definitions During comparison of ground truth and discovered behavior, we quantify model performance using the coefficient of determination, commonly referred to as the R2 -score. This metric is defined as 2 Pn true pred n y − y i i i=1 1 X true R2 := 1 − Pn , ȳ := y . (31) true − ȳ)2 n i=1 i i=1 (yi

Here, n denotes the number of data points, and y represents the quantity of interest, such as the strain energy density or a component of the stress tensor. The values yitrue and yipred correspond to the true and predicted responses, respectively, while ȳ denotes the mean of the true values.

The coefficient of determination is a global metric in the sense that it summarizes model performance over an entire dataset with a single scalar value. However, it is also useful to assess performance locally, i.e. to quantify the prediction error for an individual input. For this purpose, we use the normalized error, defined as y true − y pred , (32) ϵnorm := median(|y true |) 11

where median (•) denotes the median4 of some set of values • in the training dataset. When the value of interest is a tensor, e.g. the Kirchhoff stress, then we take the median of the norm of that tensor in the training dataset. 3. Results The results are organized to progress from interpretable, path-wise comparisons to stringent tests of generalization and practical predictive performance. We first consider sparse term recovery and path-wise agreement in Sec. 3.1, where the discovered constitutive models are compared with the ground truth along canonical loading paths; this reveals whether the CANN identifies the correct active terms and reproduces the associated stress response in mechanically interpretable settings. We then examine generalization in multi-dimensional deformation space in Sec. 3.2, motivated by the fact that one-dimensional loading paths probe only a small subset of the admissible deformation states; this analysis therefore assesses accuracy over broader domains and distinguishes performance on deformation states that are seen versus unseen during training. Finally, we perform an independent forward finite element validation in Sec. 3.3 on an unseen geometry and loading configuration to evaluate how the discovered models perform in a practically relevant setting. Collectively, these subsections assess interpretability, generalization, and predictive utility. 3.1. Sparse term recovery and path-wise agreement The stress responses of the discovered models are compared to the isotropic and anisotropic ground truth models in Figs. 4 and 5, respectively, for various loading paths. The considered loading paths include uniaxial tension (UT), confined compression (CC), biaxial tension (BT), and simple shear (SS), which are defined, respectively, as  1+γ UT  0 F (γ) = 0  1+γ BT  0 F (γ) = 0

0 √1 1+γ

 0 0 ,

√1 1+γ

0 1+γ 0

0 0

1 [1+γ]2

 1

 0 0 F CC (γ) =  0 1 0 , 0 0 1   1 γ 0 SS F (γ) = 0 1 0 , 0 0 1 1+γ



,

(33a)

(33b)

and are presented in rows 1–4, respectively, in the isotropic cases (Fig. 4), with the Neohookean, Demiray, Isihara, Arruda–Boyce, Gent–Thomas, and Ogden models being presented in columns 1–6, respectively. For the anisotropic benchmarks in Fig. 5, the loading paths are defined by the same deformation gradients in Eqs. (33); the material orientation is varied by choosing the structural vector relative to the loading axes. T T Specifically, UT is evaluated twice, with A = [1 0 0] and A = [0 1 0] , to probe loading parallel (row 1) and transverse (row 2) to the structural direction, respectively. The remaining rows show CC (row T T 3) and BT (row 4) with A = [1 0 0] , and SS (row 5) with A = [0 1 0] .

4 The median is used instead of the mean to define a robust normalization scale, since several benchmark models contain exponential terms that generate large outlying values for selected inputs.

Neohookean

P11 F CC

Arruda-Boyce Gent-Thomas 2

Ogden 2

R = 1.00 8.2

R = 1.00 2.1

R = 0.83 1.2

R = 0.98 0.6

R2 = 1.00 −2.2

R2 = 1.00 −4.1

R2 = 1.00 −8.3

R2 = 1.00 −3.9

R2 = 0.92 −3.1

R2 = 0.90 −2

R2 = 1.00 0.5

R2 = 1.00 908740

R2 = 0.93 16.4

R2 = 1.00 1.2

R2 = 0.92 0.7

R2 = 0.84

R2 = 1.00 0.7

R2 = 1.00 7

R2 = 1.00 4.9

R2 = 1.00 1.8

R2 = 0.99 1.1

R2 = 0.94

P11 F BT

Isihara 2

R = 1.00 164.8

P12 F SS

Demiray 2

R = 1.00 0.8

P11 F UT

θ1,1,1 φ1,1,1 (I˜1 − 3) θ1,1,2 exp(φ1,1,2 (I˜1 − 3)1 ) 1

θ1,2,1 φ1,2,1 (I˜1 − 3)2 θ1,2,2 exp(φ1,2,2 (I˜1 − 3)2 )

3/2 θ2,1,1 φ2,1,1 (I˜2 − 33/2 )1

3/2 θ2,1,2 exp(φ2,1,2 (I˜2 − 33/2 )1 ) 3/2

θ2,2,1 φ2,2,1 (I˜2

− 33/2 )2 3/2

θ2,2,2 exp(φ2,2,2 (I˜2

− 33/2 )2 )

0.3

0.5

1 2 1

θ3,1,1 φ3,1,1 ((J − 1) )

θ3,1,2 exp(φ3,1,2 ((J − 1)2 )1 ) θ3,2,1 φ3,2,1 ((J − 1)2 )2

θ3,2,2 exp(φ3,2,2 ((J − 1)2 )2 )

Figure 4: Comparison of discovered isotropic models first Piola-Kirchhoff stress against ground truth. Relevant components of the first Piola-Kirchhoff stress for uniaxial tension (UT), confined compression (CC), biaxial tension (BT), and simple shear (SS), as defined in Eqs. 33a and 33b, are shown in rows 1–4, respectively, while the Neohookean, Demiray, Isihara, Arruda-Boyce, Gent-Thomas, and Ogden models are shown in columns 1–6, respectively. The ground truth response is shown as a dashed black line, whereas the discovered response is presented as a stacked area plot, with colours indicating the contribution of individual terms. The coefficient of determination (R2 -score) is shown in the top left of each plot.

HGO 2

P11 F UT

GOH 2

R = 1.00 6.1

R = 0.92 17.8

R = 0.98 4.5

R2 = 1.00 1

R2 = 0.91 1.5

R2 = 0.93 1

R2 = 0.27 9.1

R2 = 1.00 1

R2 = 1.00 −1.7

R2 = 1.00 −2.2

R2 = 0.49 −3.6

R2 = 0.95 −1.7

R2 = 0.99 −2.7

R2 = 1.00 −2.2

R2 = 1.00 1.7

R2 = 0.83 1.8

R2 = 0.97 1.6

R2 = 0.96 5.3

R2 = 0.86 0.8

R2 = 1.00 1.1

R2 = 1.00 0.9

R2 = 0.81 1.3

R2 = 0.88 1.1

R2 = 0.94 1.5

R2 = 0.99 0.8

P11 F CC

HY 2

R = 0.96 3.9

P11 F BT

MO 2

R = 1.00 40.8

P12 F SS

Meaney 2

R = 1.00 4.5

P11 F UT

0.0

0.5 0.0

0.5 0.0 3/2

θ1,1,1 φ1,1,1 (I˜1 − 3)1 θ1,1,2 exp(φ1,1,2 (I˜1 − 3)1 )

θ2,2,1 φ2,2,1 (I˜2

θ1,2,2 exp(φ1,2,2 (I˜1 − 3)2 )

θ3,1,2 exp(φ3,1,2 ((J − 1)2 )1 )

θ1,2,1 φ1,2,1 (I˜1

− 3)2

3/2 θ2,1,1 φ2,1,1 (I˜2 − 33/2 )1

− 33/2 )2

0.5 0.0

3/2 θ2,2,2 exp(φ2,2,2 (I˜2 − 33/2 )2 ) θ3,1,1 φ3,1,1 ((J − 1)2 )1

θ3,2,1 φ3,2,1 ((J − 1)2 )2

0.5 0.0

0.5

θ3,2,2 exp(φ3,2,2 ((J − 1)2 )2 )

θ4,1,1 φ4,1,1 (< I4,11 − 1 >2 )1

θ4,1,2 exp(φ4,1,2 (< I4,11 − 1 >2 )1 ) θ4,2,1 φ4,2,1 (< I4,11 − 1 >2 )2

θ4,2,2 exp(φ4,2,2 (< I4,11 − 1 >2 )2 )

3/2 θ2,1,2 exp(φ2,1,2 (I˜2 − 33/2 )1 )

Figure 5: Comparison of discovered anisotropic models stress against ground truth. The first Piola–Kirchhoff stress for uniaxial tension (UT) with A = [1 0 0]T , UT with A = [0 1 0]T , confined compression (CC) with A = [1 0 0]T , biaxial tension (BT) with A = [1 0 0]T , and simple shear (SS) with A = [0 1 0]T , are shown in rows 1–5, respectively, while the Anisotropic Neohookean (AN), Holzapfel–Gasser–Ogden (HGO), Meaney, Merodio–Ogden (MO), Humphrey–Yin (HY), and Gasser–Ogden–Holzapfel (GOH) models are shown in columns 1–6, respectively. The ground-truth behavior is shown as a single dashed black line, whereas the discovered behavior is presented as a stacked area plot with colours indicating the contribution of individual terms as specified in the legend. Additionally, the coefficient of determination (R2 -score) between the true and discovered behavior for each model and loading path is shown in the top left corner of each corresponding plot.

Overall, the discovered models exhibit excellent agreement with the ground truth (shown with a dashed line), with coefficients of determination R2 exceeding 0.95 in the majority of cases (Figs. 4 and 5). Ground-truth laws that are representable by the chosen CANN basis are near-exactly recovered. In the isotropic setting this is observed for the Neohookean and Demiray models where the isochoric terms θ1,1,1 ϕ1,1,1 [I˜1 −3] and θ1,1,2 exp(ϕ1,1,2 [I˜1 −3]) are correctly identified, respectively, and the volumetric term θ3,1,1 ϕ3,1,1 [J−1]2 is correctly identified for both models (Fig. 4 columns 1 and 2, respectively). In the anisotropic case, this was demonstrated for the anisotropic Neohookean (AN) and Holzapfel-Gasser-Ogden (HGO) ground truth models where the CANN correctly identifies the anisotropic terms θ4,1,1 ϕ4,1,1 ⟨I˜4 −1⟩2 and 2 θ4,1,2 exp ϕ4,1,2 ⟨I˜4 −1⟩ , respectively, as well as the isotropic terms θ1,1,1 ϕ1,1,1 [I˜1 −3] and θ3,1,1 ϕ3,1,1 [J−1]2 for both models (Fig. 5 columns 1 and 2, respectively). The recovery of exponential terms is particularly important, because the corresponding material parameter appears inside a nonlinear function; this goes beyond the standard EUCLID formulation, which identifies linear coefficients of fixed library terms but does not directly calibrate parameters embedded inside nonlinear functions [52]. Shared terms are recovered even when the full ground-truth law is not representable by the chosen CANN basis. For example, in the isotropic setting, the Isihara model and chosen CANN architecture both contain the terms [I˜1 −3], [I˜1 −3]2 , and [J−1]2 , but the Isihara model also contains the term [I˜2 −3] which is not contained in the CANN basis. Nevertheless, the terms [I˜1 −3], [I˜1 −3]2 , and [J−1]2 are correctly identified by the CANN. A similar pattern is observed in the anisotropic setting: in the Meaney, Merodio–Ogden (MO), and Gasser–Ogden–Holzapfel (GOH) models, CANN-EUCLID correctly identifies the term θ1,1,1 ϕ1,1,1 [I˜1 −3], while in the Humphrey–Yin (HY) model it correctly identifies the term θ1,1,2 exp(ϕ1,1,2 [I˜1 −3]). This occurs despite the fact that these models contain additional terms that are not present in the chosen CANN basis. The only exception is the MO model: although [I˜1 −3] is the only I˜1 -based term in the ground-truth law and is correctly identified, CANN-EUCLID also activates the additional term θ1,1,2 exp(ϕ1,1,2 [I˜1 −3]) in the discovered MO behavior. Apart from this exception, all I˜1 - and J -based terms that are present in the CANN basis but absent from the Meaney, MO, GOH, and HY ground-truth models are correctly omitted where appropriate. Ground truth terms that are absent from the CANN basis are typically approximated by other terms. This is exemplified by the Isihara model, which contains the non-polyconvex term [I˜2 −3]. In this case, i h 3/2 3/2 ˜ CANN-EUCLID instead identifies the polyconvex term θ2,1,1 ϕ2,1,1 I2 −3 , resulting in comparable pathwise behavior (Fig. 4, column 3). Moreover, complex functional forms, such as logarithmic terms, hyperbolic sine functions, and principal stretches raised to non-integer powers, are approximated by simpler terms contained in the chosen CANN basis. This is demonstrated by CANN-EUCLID’s ability to approximate the response of the Gent–Thomas, Arruda–Boyce, and Ogden models with just two terms: θ1,1,1 ϕ1,1,1 [I˜1 −3] and θ3,1,1 ϕ3,1,1 [J−1]2 (Fig. 4, columns 4–6). Similarly, all anisotropic I˜4 - and I˜5 -based terms in the Meaney, MO, HY, and GOH models are absent from the chosen CANN basis. Nevertheless, CANN-EUCLID approximates the anisotropic stress contributions with good accuracy in loading cases where these terms are active, i.e., i h T UT with A = [1 0 0] , BT, and SS. For the Meaney model, the term I˜4 2 +2I˜4 −1 −3 is approximated by 2 2 θ4,1,1 ϕ4,1,1 ⟨I˜4 −1⟩ (Fig. 5, column 3). For the MO model, the term [I˜5 −1] is approximated by a combination 2 2 4 of θ4,1,1 ϕ4,1,1 ⟨I˜4 −1⟩ , θ4,1,2 exp ϕ4,1,2 ⟨I˜4 −1⟩ , and θ4,2,1 ϕ4,2,1 ⟨I˜4 −1⟩ (Fig. 5, column 4). For the HY model, the D√ E2 term exp 8 I˜4 −1 −1 is approximated by θ4,1,2 exp ϕ4,1,2 ⟨I˜4 −1⟩2 and θ4,2,1 ϕ4,2,1 ⟨I˜4 −1⟩4 (Fig. 5, column 5).

Finally, for the GOH model, the term 61 exp 3 I61 + I24 −1 ˜

4 −1 is approximated by θ4,2,1 ϕ4,2,1 ⟨I˜4 −1⟩ (Fig.

column 6). 3.2. Generalization in multi-dimensional deformation space Observing the discovered behavior along distinct loading paths as shown in Figs. 4 and 5 is useful for elucidating which terms contribute to the behavior and for providing a cursory comparison between the ground-truth and discovered behavior. However, evaluating the performance of the discovered behavior in this manner has at least two notable limitations: 15

Neohookean

log(λ̃2 )

Isotropic

0.5

BT BT–2:1

Demiray

Isihara

BT BT–2:1

0.3

SS −0.4

−0.3

0.0

0.5

log(λ̃1 )

HGO

log(F̃22 )

Anisotropic

0.9

0.0

UT 0.5

1.0

−0.8

BT–2:1

UT 1

log(λ̃1 )

GOH

0.5 UT

−0.8

2 0

log(λ̃1 )

MO 0.6

BT UT

BT–2:1 BT

BT–2:1

0.5

0.0 0.0

UT 0.5

Meaney

0.4

0.5 BTBT–2:1

log(λ̃1 )

0.9

BT UT

log(F̃11 )

−0.5

1.0 0.0

Ogden

BT 0.4 BT–2:1

log(λ̃1 )

BT–2:1 BT–2:1

−0.4

Gent-Thomas

BT BT–2:1

log(λ̃1 )

0.3

BT–2:1

Arruda-Boyce 0.4

−0.1 0.5−0.25 0.00

0.25

log(F̃11 ) 100

−0.1

−0.5

0.0

log(F̃11 )

0.5−0.25

200

−0.2

−0.1 0.00

0.25

log(F̃11 ) 300

−0.25 0.00

0.25

0.0

log(F̃11 )

400

0.5

log(F̃11 ) 500

No. material points

Figure 6: Deformation states seen in training data. In the isotropic cases (top), the strain energy density function and principal stresses can be written in terms of the principal stretches. The frequency of the first and second isochoric principal stretches, λ̃1 and λ̃2 , in the training data is shown for the Neohookean, Demiray, Isihara, Arruda–Boyce, Gent–Thomas, and Ogden models in columns 1–6, respectively. In the anisotropic cases (bottom), the strain energy density function and principal stresses cannot, in general, be written solely in terms of the principal stretches. We therefore plot the frequency of the F̃11 and F̃22 components of the isochoric deformation gradient. The anisotropic Neohookean (AN), Holzapfel–Gasser–Ogden (HGO), Meaney, Merodio–Ogden (MO), Humphrey–Yin (HY), and Gasser–Ogden–Holzapfel (GOH) models are shown in columns 1–6, respectively. The paths corresponding to uniaxial tension (UT), equibiaxial tension (BT), and simple shear (SS), as defined in Eqs. (33), are indicated by dashed orange lines; the simple-shear path is not shown in the anisotropic case because it degenerates to a point in this projection. The line corresponding to biaxial tension with a 2:1 ratio, i.e., the loading applied to the training sample in Fig. 3, is indicated by the annotation “BT–2:1”. The figure illustrates that a single heterogeneous test in the unsupervised paradigm induces deformation states associated with multiple conventional supervised loading paths. It also shows that the resulting dataset is imbalanced, with most deformation states clustering around the BT–2:1 line.

(i) One-dimensional lines through a six-dimensional deformation space. General three-dimensional hyperelastic behavior depends on six independent variables, i.e., the six independent components of the left or right Cauchy–Green tensor. The loading paths described in Eqs. (33) and used for Figs. 4 and 5 define one-dimensional lines through this deformation space. These lines necessarily omit large regions of the six-dimensional domain that describes the possible deformation state at a material point. As such, evaluating the behavior on these one-dimensional paths alone is not sufficient, as it excludes large regions of deformation space that the discovered models may encounter during deployment, e.g., in a finite element solver. This is illustrated in Fig. 6, where uniaxial tension (UT), biaxial tension (BT), and simple shear (SS) are shown as dashed orange lines overlaid on frequency plots of the deformation states observed in the training data for each ground-truth model. Since the full sixdimensional space cannot be visualized directly, these paths are projected onto two-dimensional manifolds corresponding to the isochoric principal stretches in the isotropic cases and to two components of the deformation gradient in the anisotropic cases. These projected paths occupy only small regions of the sampled deformation space, highlighting domains in which the accuracy of the discovered models is not assessed by path-wise comparisons alone.

(ii) Lack of discrimination between predictions on seen and unseen deformation states. In supervised model calibration or discovery using, for example, uniaxial test data, the seen deformation states can be defined naturally from the maximum stretch reached during training. Generalization can then be assessed by evaluating the discovered model at stretches beyond those observed during training. In the context of the loading paths used in Figs. 4 and 5, this would amount to choosing a value of γ corresponding to the maximum stretch seen during training. However, in the unsupervised case, the deformation states in the training data are diverse and heterogeneous. It is therefore not immediately clear which range of γ values should be considered seen during training and which should be considered unseen. To address these points, we assess the accuracy of the discovered behavior on higher-dimensional domains of the possible deformation space that are seen and unseen during training in Fig. 7. Since effectively visualizing a three- or six-dimensional domain is intractable on a two-dimensional page or screen, we choose to assess the behavior on two-dimensional manifolds through the six-dimensional deformation space. The normalized error, as defined in Eq. (32), in the first, second, and third principal Kirchhoff stresses, are shown in rows 1–3, respectively, for the isotropic cases (top) and anisotropic cases (bottom) of Fig. 7. In the case of isotropic behavior, we use the manifold defined by isochoric deformations as done in Fig. 6; i.e. we use λ̃1 and λ̃2 and set λ̃3 =1/[λ̃1 λ̃2 ]. In the anisotropic case, we restrict ourselves to visualizing the F̃11 and F̃22 components of the isochoric deformation gradient whilst noting that this is only a slice through the six-dimensional domain of independent variables. Moreover, the sampled deformation gradients are chosen T to be diagonal with F̃33 =1/[F̃11 F̃22 ] and the direction of the structural vector is taken to be A = [0 1 0] . The sampled domains are chosen to evaluate the performance and generalizability of the models on seen and unseen deformation states. The domains are chosen by obtaining the convex hull of the deformation states used during training as shown in Fig. 6. The sampling domain is then obtained by simply scaling the points on the perimeter of the convex hull by a factor cgen (1 < cgen ). We make use of log values since this places the undeformed location at the origin, which allows for increasing the sampling domain by simply scaling the points without translating the undeformed location5 . The larger the value of cgen the larger the domain of unseen deformation states; in this work, we choose cgen = 1.5. The convex hull of the seen deformation states is shown as a grey line in Fig. 7 and the associated sampled domains are clear from the colored areas. Additionally, the coefficients of determination on the seen Rs2 and unseen Ru2 data are shown in the top-right corner of each normalized error plot. To aid in the analysis and discussion of Fig. 7, the median and maximum of the stress prediction errors in the seen and unseen deformation domains are presented as a bar plot for each model in Fig. 8. Exceptional performance on seen and unseen deformation states is obtained when the ground truth models are recovered near-exactly. This is demonstrated by the Neohookean and Demiray models in the isotropic case and the anisotropic Neohookean (AN) and Holzapfel–Gasser–Ogden (HGO) models in the anisotropic cases, all of which were shown to be captured exactly in Sec. 3.1, Figs. 4 and 5. Here, these models all achieve R2 -scores of 1.00 for all three principal stresses on seen and unseen deformations (Fig. 7, columns 1 and 2, Anisotropic and Isotropic). Additionally, the highest normalized error between all of them on seen deformation states is below 1.1 × 10−3 (Fig. 8). However, despite the near-exact recovery of the ground truth models, and median errors on unseen deformations being below 5.4 × 10−4 (Fig. 8), the Demiray and HGO models, i.e. the models containing exponentials, still result in maximum normalized errors of 69 and 4.4 × 103 on the unseen deformation domain, respectively. This highlights the sensitivity of exponential-based models and provides a first indication of the challenges surrounding generalizability in the presence of such terms. The median error in the seen deformation domain is low across all models. This holds true even 5 Scaling compressive stretches in the log space results in more compressive stretches (log(λ̃ ) < 0 ⇒ c gen log(λ̃2 ) < log(λ̃2 )), 2

and similarly for tensile stresses (0 < log(λ̃1 ) ⇒ log(λ̃1 ) < cgen log(λ̃1 )). If the log of the stretches was not used, scaling would result in all stretches becoming more tensile (λ̃1 < cgen λ̃1 ).

Neohookean

Demiray

−0.5

0.3

0.4

−0.2

−0.4

−0.5

−1.0

−0.9

0.3

0.4

0.5

−0.2

−0.4

−0.5

−1.0

−0.9

0.3

0.4

0.5

−0.2

−0.4

−0.5

−1.0

−0.9

log(λ̃2 )

τ2 τ3

−0.5 0

0.0

0.5

log(λ̃1 )

0.0

0.5

log(λ̃1 )

HGO

1.0

log(λ̃1 )

Meaney

0.6

1.3

0.8

0.0

−0.1

1.2

0.6

1.3

0.8

0.0

−0.1

1.2

0.6

1.3

0.8

0.0

−0.1

log(λ̃1 )

GOH

0.7

−0.1

−0.3

0.7

−0.1

−0.3

0.7

−0.1

−0.3

log(F̃22 )

1.2

log(F̃22 )

τ1

Ogden

0.3

0.6

τ2

Gent-Thomas

log(λ̃2 )

Isotropic

−0.5

Anisotropic

Arruda-Boyce

0.5

0.6

τ3

Isihara

log(λ̃2 )

τ1

0.6

0.0

log(F̃11 )

0.5

0.00

0.25

log(F̃11 ) 0.2

0.50

0.0

0.5

0.00

0.25

log(F̃11 )

0.4

0.6

0.50

0.0

0.5

log(F̃11 ) 0.8

log(F̃11 ) 1.0

normalized error

Figure 7: Accuracy of discovered CANNs evaluated in isochoric deformation space. The normalized error of the first, second, and third principal Kirchhoff stresses discovered by the CANN as compared to the ground truth are displayed in rows 1–3, respectively, for both the isotropic (top) and anisotropic (bottom) cases, with results for different material models grouped in columns. Plots in the isotropic cases (top) use the log of the first and second isochoric principal h i stretches for the x– and y–axes respectively; the third isochoric principal stretch is implicitly defined through λ̃3 = 1/ λ̃1 λ̃2 . In the anisotropic

cases (bottom), the F̃11 and F̃22 components of the isotropic deformation gradient are used for the axes. Moreover, the sampled h i deformation gradients are chosen to be diagonal with F̃33 = 1/ F̃11 F̃22 and the direction of the structural vector is taken

to be A = [0 1 0]T . The convex hull of the training data (visible in Fig. 6) is indicated by a grey line and the evaluation domain has been extended beyond the convex hull to assess the generalizability of the discovered CANNs. The R2 -scores for 2 are reported in the top-right corner of each plot. the seen data Rs2 and unseen data Ru

Figure 8: Statistics for CANNs stress prediction accuracy on seen and unseen deformations. The median and maximum of the normalized errors in the Kirchhoff principal stresses presented in Fig. 7 are shown for the seen and unseen domains of deformation for the isotropic and anisotropic cases.

in cases where the ground truth contains a pseudo-invariant not used in the CANN basis, e.g. the Isihara and Merodio–Ogden (MO) models, which have median errors of 4.3 × 10−3 and 1.8 × 10−2 , respectively (Fig. 8). It also holds when the ground truth has a substantially different functional form than the activation functions in the CANN, e.g. the Ogden and the Meaney models, with median errors of 6.2 × 10−2 and 4.9 × 10−2 , respectively (Fig. 8). Moreover, regions of the seen domain with more training data points tend to exhibit lower errors than sparsely sampled regions, such as those close to the boundary of the seen domain. For example, the errors for the HY and GOH models in Fig. 7 are much larger close to the convex hull of the seen deformation domain; these regions correspond to deformation states reached by fewer material points in Fig. 6. The low number of data points in these deformation states gives them a smaller relative contribution to the loss function, which likely contributes to the larger normalized errors observed in these regions. This suggests a practical avenue for increasing the accuracy of the discovered model in deformation regimes of interest: design heterogeneous tests that drive more material points into those regimes. Generalization performance is highly dependent on the presence or absence of exponential terms. This is most clear for the HY and GOH models, both of which contain exponential terms that are not contained in the chosen CANN basis. In these cases, the median errors in the unseen domain are 0.41 and 2.2, respectively. By contrast, the median errors in the unseen domain for the Meaney and MO cases are 0.11 and 0.043, respectively. These results indicate that reliable extrapolation cannot be assumed when the sampled deformation states do not sufficiently probe a potential strain-stiffening regime that governs the material response outside the training domain. 3.3. Independent forward finite element validation The solution to the finite strain elasticity BVP inherently minimizes the total stored elastic energy. As such, many of the deformation states present in Fig. 7 may be unlikely to arise in a forward FE simulation because of their excessively large associated strain energy density values. Hence, evaluating the accuracy of the discovered models by sampling deformation space alone, as done in Sec. 3.2, may give a more pessimistic impression of their practical accuracy. To estimate the accuracy that can be expected when using the discovered material models in a forward FE simulation, we perform an independent validation 19

Neohookean

Isihara

Arruda-Boyce

Gent-Thomas

Ogden

0.050 0.025

1.00

100

25 0 0

cumulative density

Demiray

0.075

predicted stress

norm. error

0.100

R2 = 1.00

true stress

R2 = 1.00

100

true stress

R2 = 1.00

true stress

R2 = 1.00

true stress

R2 = 0.98

true stress

R2 = 0.97 10

true stress

0.75 0.50 0.25 0.00 10−3 10−2 10−1

10−3 10−2 10−1

norm. error

Figure 9: Comparison of ground truth and discovered stress fields in validation simulations for isotropic behavior. The normalized error in the Kirchhoff stress defined in Eq. (32) is shown on the deformed configuration in the top row, parity plots of principal Kirchhoff stresses are shown in the middle row, and cumulative density plots of the normalized error are shown in the bottom row along with the median and maximum values. Results for the six ground truth models are shown column-wise.

simulation and compare the stress field obtained with the ground-truth behavior against that obtained with the discovered behavior. The loading and geometric parameters of the validation case are chosen to represent a challenging scenario that induces extreme deformations, as illustrated in Fig. 3: the geometry contains elliptical holes that induce stress concentrations and is stretched by 100 % of its height. Additionally, in the T anisotropic cases, the structural vector is chosen to be in the positive y-direction, i.e. A = [0 1 0] . The results of these simulations are displayed in Figs. 9 and 10. Ground-truth models are grouped by column, normalized stress errors are shown spatially on the deformed geometries in the first row, the three principal Kirchhoff stresses are displayed collectively in parity plots in the second row, and the cumulative distribution of normalized stress errors is shown in the third row. Forward validation errors are substantially lower than errors obtained from uniform sampling of the multidimensional deformation space. This is particularly clear for the two least accurate isotropic cases, namely the Gent–Thomas and Ogden models, for which the median normalized stress errors are only 6.0 × 10−3 and 7.6 × 10−3 , respectively (Fig. 9, bottom row, columns 5 and 6). These errors are more than an order of magnitude lower than the median errors on the unseen cgen ∈ [1.0, 1.5] multidimensional deformation states from Sec. 3.2, which were 0.29 and 0.26 for the Gent–Thomas and Ogden models, respectively, in Fig. 8. Similarly, the validation simulations for the anisotropic Meaney, HY, and GOH models yield median errors that are approximately one order of magnitude lower than those obtained from direct sampling of the deformation domain. Median errors of 0.011, 0.036, and 0.037 are observed for the validation simulations of the Meaney, HY, and GOH models, respectively (Fig. 10, columns 3, 5, and 6, respectively), compared to median errors of 0.11, 0.41, and 2.2 on unseen cgen ∈ [1.0, 1.5] multi-dimensional 20

Meaney

GOH

0.050 0.025

200

1.00

600

1000

400 400

200

100

500

100

cumulative density

HGO

0.075

predicted stress

norm. error

0.100

200

200 R2 = 1.00

true stress

200

0 0

R2 = 1.00

1000

true stress

R2 = 0.99

100

true stress

R2 = 0.94

250

true stress

R2 = 0.96 500

true stress

R2 = 0.92 0

true stress

500

0.75 0.50 0.25 0.00 10−3 10−2 10−1

10−3 10−2 10−1

norm. error

Figure 10: Comparison of ground truth and discovered stress fields in validation simulations for anisotropic behavior. The normalized error, as defined in Eq. (32), for the final loading step is shown on the deformed geometry in the first row, parity plots for predicted and true principal Kirchhoff stresses are shown in the second row together with R2 -scores, with frequency indicated by colour, and cumulative density plots of the number of material points versus normalized error are shown in the third row. Results for the anisotropic Neohookean (AN), Holzapfel–Gasser–Ogden (HGO), Meaney, Merodio– Ogden (MO), Humphrey–Yin (HY), and Gasser–Ogden–Holzapfel (GOH) models are shown in columns 1–6, respectively.

deformations in Fig. 7 (Fig. 8). This discrepancy arises because the statistics in Fig. 8 are computed from uniform sampling of the deformation domain and therefore do not reflect the likelihood of individual deformation states occurring in a forward simulation. By contrast, the error statistics obtained from the validation simulations naturally account for how frequently different deformation states arise in the solution, making the resulting median errors more representative of practical FE predictions. Within the validation simulations, median errors remain low, while maximum errors reveal localized or dispersed model-specific discrepancies. In the isotropic cases, this is evident from the predicted-to-true stress parity plots, in which the predicted principal stresses lie close to the identity line for all six models (Fig. 9, middle row). The coefficient of determination is R2 = 1.00 for the Neohookean, Demiray, Isihara, and Arruda–Boyce cases, and remains high at R2 = 0.98 and R2 = 0.97 for the Gent– Thomas and Ogden cases, respectively (Fig. 9, middle row). Moreover, the median normalized stress errors remain small across all isotropic models: below 8.3 × 10−4 for the Neohookean, Demiray, Isihara, and Arruda–Boyce models, and below 7.6 × 10−3 for the more challenging Gent–Thomas and Ogden models (Fig. 9, bottom row). Thus, even for the two least accurate isotropic validation cases, the median error remains below 1 %. The spatial error distributions further show that the larger discrepancies in the isotropic cases are localized rather than widespread. For the Neohookean, Demiray, Isihara, and Arruda–Boyce models, the error is nearly uniformly close to zero, whereas for the Gent–Thomas and Ogden models, the larger stress errors are confined to relatively small regions near the stress concentrations around the holes (Fig. 9, top row). In the anisotropic cases, the near-exactly recovered anisotropic Neohookean (AN) and Holzapfel–Gasser–Ogden (HGO) models again display extremely low median normalized errors of 4.4 × 10−9 21

and 8.2 × 10−8 , respectively (Fig. 10 columns 1 and 2, respectively). The remaining anisotropic models also demonstrate low median errors of 0.011, 0.028, 0.036, and 0.037 for the Meaney, Merodio–Ogden (MO), Humphrey–Yin (HY), and Gasser–Ogden–Holzapfel (GOH) models, respectively (Fig. 10 columns 3–6, respectively). However, these cases also exhibit large maximum normalized errors, with values of 2.6, 20, 17, and 24 for the Meaney, MO, HY, and GOH models, respectively. Again, larger errors tend to arise in the presence of exponential terms: the HY and GOH models contain exponential terms in the ground truth, while in the MO case the exponential term θ4,1,2 exp ϕ4,1,2 ⟨I˜4 −1⟩2 is included in the discovered CANN. Moreover, these errors tend to be more spatially dispersed in the presence of exponential terms: the maximum error in the Meaney model is highly localized, whereas the larger errors in the MO, HY, and GOH cases are more spread out. 4. Discussion The central result of this study is that sparse, interpretable CANN-based constitutive laws can be discovered directly from full-field kinematic data and reaction forces. Across the isotropic and anisotropic benchmarks, CANN–EUCLID recovered the correct constitutive structure when the ground-truth law was contained in the chosen basis, produced useful approximations when it was not, and revealed clear limits of extrapolation when relevant deformation regimes were insufficiently sampled. These results support stress-unsupervised full-field discovery as a promising route for constitutive model identification, while also highlighting the importance of basis design, deformation-state coverage, and forward validation. When the ground-truth constitutive law was contained in the chosen CANN basis, the framework recovered the correct active terms with near-exact accuracy. This was observed for the isotropic Neohookean and Demiray models and for the anisotropic Neohookean and Holzapfel–Gasser–Ogden models. In these cases, the discovered CANNs selected the relevant isochoric, volumetric, and anisotropic terms while suppressing inactive alternatives. Importantly, this included exponential terms whose material parameters appear inside nonlinear functions. This capability goes beyond classical sparse linear-library discovery and is particularly relevant for soft biological tissues, where exponential strain-stiffening terms are common. Moreover, for these exactly representable cases, near-exact term recovery also translated into excellent generalization over the evaluated deformation domains. The framework also produced useful constitutive approximations when the exact ground-truth law was not representable by the selected CANN basis. In these cases, CANN–EUCLID generally retained shared terms that were present in both the ground truth and the basis, and compensated for missing contributions using the closest available basis functions. This behavior was observed for the Isihara, Arruda–Boyce, Gent– Thomas, and Ogden models in the isotropic setting, and for the Meaney, Merodio–Ogden, Humphrey– Yin, and Gasser–Ogden–Holzapfel models in the anisotropic setting. These results indicate that exact representability is not strictly required for useful model discovery. However, they also show that the quality of the discovered approximation depends on the expressiveness of the basis and on the deformation states sampled during training. A central finding of this study is that constitutive recovery and constitutive generalization must be evaluated separately. The CANN architecture used here retains important physical structure, including objectivity, material symmetry, zero energy and zero stress in the reference configuration, and polyconvexity for the selected pseudo-invariants and activation functions. Nevertheless, these physical constraints do not by themselves guarantee accurate extrapolation outside the deformation states represented in the training data. This distinction was clear in the out-of-basis benchmarks. Acceptable generalization was obtained for the Isihara, Arruda–Boyce, Meaney, and Merodio–Ogden models, whereas poorer generalization was observed for the Gent–Thomas, Ogden, Humphrey–Yin, and Gasser–Ogden–Holzapfel models. Even within a single model, generalization was not uniform: for several benchmarks, including Gent–Thomas, Ogden, Meaney, and Merodio–Ogden, the discovered laws performed well in some unseen regions of deformation space but poorly in others. This highlights the importance of evaluating discovered constitutive laws over multidimensional deformation domains rather than along a small number of canonical loading paths alone. 22

Generalization was particularly sensitive in the presence of exponential-type terms. When the ground-truth response contained an exponential contribution compatible with the CANN basis, as in the Demiray and Holzapfel–Gasser–Ogden benchmarks, the corresponding term could be recovered accurately. However, even small discrepancies in exponential terms produced large errors outside the sampled deformation domain because of their rapid growth. Conversely, when exponential terms were introduced by the discovered approximation but were not part of the true response, extrapolation also deteriorated, as observed most clearly in the anisotropic benchmarks. This does not imply that exponential terms should be avoided. Rather, it shows that exponential terms require careful interpretation, sufficient probing of the stiffening regime during training, and explicit generalization checks over relevant deformation domains. The results further highlight the importance of the training deformation-state distribution. The heterogeneous test we considered here generated a much richer set of deformation states than a conventional homogeneous test, but this set was still imbalanced, with many material points clustered around the dominant biaxial-tension-with-2:1-ratio loading direction. The discovered models were most accurate in densely sampled regions and less reliable near the boundaries of the observed domain or in sparsely populated regions. This suggests that robust constitutive discovery requires not only an expressive and physically constrained model class, but also experimental designs that deliberately populate the deformation regimes in which the model is expected to be used. Recent specimen- and topology-optimization approaches for one-shot material identification provide a natural route toward this goal. From this perspective, the present findings provide a strong motivation for the stress-unsupervised full-field paradigm. The challenge of generalization is not unique to unsupervised discovery. If anything, it is more restrictive in conventional supervised calibration, where data are typically obtained from a sparse collection of one-dimensional or nominally homogeneous loading paths. By contrast, a single heterogeneous full-field experiment can expose the material model to a substantially broader and more application-relevant portion of deformation space. Since accurate constitutive predictions cannot generally be expected far outside the observed domain, improving discovery requires expanding and balancing the domain of seen deformations. This is precisely where unsupervised full-field discovery offers a fundamental advantage over traditional stress-supervised testing. Independent forward finite element validation provided a more application-relevant perspective on these generalization errors. Although uniform sampling of the multidimensional deformation space revealed large errors in some unseen regions, the median stress errors in forward simulations remained low for most benchmarks. Larger discrepancies were typically localized near stress concentrations and occurred in deformation states that were less frequently realized by the mechanical boundary value problem. Thus, broad deformation-space error maps are essential for diagnosing the limitations of a discovered law, but they may be conservative relative to practical finite element use, where mechanically realized deformation states occupy only a subset of the full admissible domain. Several extensions follow naturally from this work. First, the framework should be evaluated with broader CANN bases, including generalized invariants, logarithmic terms, principal-stretch-based terms, and other features known to be important in rubber-like materials and soft tissues. Second, CANN–EUCLID should be coupled to optimized experimental designs that produce broader and more balanced coverage of applicationrelevant deformation states. Third, the framework should be extended toward realistic biological settings with spatially varying fiber architectures. More broadly, the long-term opportunity is to move from discovering a single homogeneous constitutive law toward discovering spatially informed, physics-constrained constitutive representations directly from rich experimental fields. Such developments would make unsupervised CANN discovery a practical route toward interpretable material laws for complex and intrinsically heterogeneous soft tissues. 5. Conclusions In this work, we introduced CANN–EUCLID, a stress-unsupervised framework for discovering sparse hyperelastic constitutive laws from full-field displacement data and reaction forces. By combining the term-level 23

interpretability of constitutive artificial neural networks (CANNs) with the equilibrium-based training strategy of EUCLID, the framework enables constitutive model discovery without local stress measurements, prescribed stress–strain pairs, or a fixed constitutive ansatz. Across isotropic and anisotropic benchmark problems, the results show that CANN–EUCLID can recover compact, interpretable constitutive models directly from a single heterogeneous loading experiment. These findings establish CANN–EUCLID as a promising framework for interpretable, stress-unsupervised constitutive model discovery from full-field data. Its main advantage is not simply that it avoids stress measurements, but that it allows heterogeneous experiments to be used directly as information-rich discovery problems. This is especially attractive for soft biological tissues, where local stresses are difficult to measure, repeated testing can alter the specimen, and pooling data across samples is complicated by biological variability. In such settings, a single well-designed heterogeneous experiment can provide a substantially richer basis for constitutive discovery than a collection of nominally homogeneous loading paths. Data availability The unsupervised training framework developed in this work will be released as a module in the COMMET codebase https://github.com/COMMET-code after publication. In addition to the unsupervised training of CANNs, the module provides support for training of other neural constitutive models including ICNNs [79, 20], ICKANs [80, 21], and custom user-defined neural constitutive models.

Appendix A. Polyconvexity of chosen CANN architecture A strain energy density function Ψ (F ) is said to be polyconvex if there exists a convex function W : 3×3 × 3×3 × + → 0+ such that + +

Ψ (F ) = W (F , adj (F ) , det (F )) ,

(A.1)

where 3×3 denotes the set of second order 3 × 3 tensors with positive determinant, + denotes the set + of positive real values, 0+ denotes the set of non-negative real values, and adj (•) denotes the adjoint of •. This condition has been shown in [64] to be sufficient for guaranteeing a solution to a finite strain elastic boundary value problem with suitable boundary conditions while also allowing for the existence of bifurcations which are physically observed in e.g. elastic buckling. Additionally, this condition can be verified for a given strain energy function with relative ease in comparison to the weaker sufficient condition of quasiconvexity [81, 66]. For these reasons, the use of polyconvex strain energy functions has arguably become the standard approach for guaranteeing a solution to a finite strain elastic boundary value problem in the community.

To keep with this trend, we choose the architecture and pseudo-invariants of the CANN such that the resulting discovered models are guaranteed to be polyconvex. Note that the architectural choices listed in Sec. 2.1.3 lead to a model consisting of the summation of the following terms: h ij θ1,j,1 ϕ1,j,1 I˜1 − 3 , h 3/2 ij θ2,j,1 ϕ2,j,1 I˜2 − 33/2 , 2j

θ3,j,1 ϕ3,j,1 [J − 1] , D E2j θ4,j,1 ϕ4,j,1 I˜4 − 1 , ij h θ1,j,2 exp ϕ1,j,2 I˜1 − 3 −1 , h 3/2 ij θ2,j,2 exp ϕ2,j,2 I˜2 − 33/2 −1 , h i 2j θ3,j,2 exp ϕ3,j,2 [J − 1] −1 , D E2j ˜ θ4,j,2 exp ϕ4,j,2 I4 − 1 −1 .

(A.2) (A.3) (A.4) (A.5) (A.6) (A.7) (A.8) (A.9)

Hence, the discovered model is guaranteed to be polyconvex if θijk , ϕijk ≥ 0 and j ∈ N, where N is the set of natural numbers, for the following reasons: (i) Eqs. (A.2) and (A.3) are proven to be polyconvex in [66] (see Lemma C.6 6 ). (ii) Eq. (A.4) is convex in det (F ) for j ∈ N. Hence, it is also polyconvex. (iii) It has been shown that I˜4 is polyconvex in [66] Lemma C.3, 27 . Additionally, the function h (•) = 2j ⟨• − 1⟩ is convex and monotonically non-decreasing in • for j ∈ N. Hence, Eq. (A.5) is polyconvex. (iv) Since exp (•) is convex and monotonically increasing in • and Eqs. (A.2)–(A.5) are polyconvex, it follows that Eqs. (A.6)–(A.9) are also polyconvex. (v) The summation of polyconvex functions is polyconvex. 3/2 ||F ||2 ||adj(F )||3 = I˜1 and det(F )2 = I˜2 as written in [66] Lemma C.6. det(F )2/3 m T 7 Note that tr(F F A⊗A) ˜ = I4 when m = 1 as written in [66] Lemma C.3. 1/3 det(F T F )

6 Note that

Appendix B. Multistage training details and hyperparameters Despite the utility of Lp regularization, it also introduces some challenges. Such challenges include: (i) The final selection of terms being highly sensitive to model initialization [70]. (ii) The “optimized” parameters being suboptimal for the primary objective (i.e. minimizing the imbalance of equilibrium) due to the fact that the Lp -norm must be minimized simultaneously. To address these challenges, we take a three-stage training approach: (Stage 1) Pre-regularization. The sensitivity to model initialization arises from the fact that terms closer to zero are penalized more strongly by Lp regularization than terms further from zero. Consequently, terms close to zero by happenstance after initialization are arbitrarily penalized more than others. To reduce this sensitivity, we first train without imposing Lp regularization, i.e. we set λp = 0. Doing so allows each term to naturally adjust its contribution to the overall behavior from its initialized state before any term is penalized. (Stage 2) Regularization and subset selection. We then apply Lp regularization to induce sparsity for a chosen number of epochs. At the end of this stage, we select all terms that have a non-negligible contribution to the overall behavior. This is determined by evaluating the strain energy density for each deformation gradient in the dataset and keeping track of the proportion that is contributed by each term. Terms that have an average proportional contribution above a small chosen cut-off threshold are selected and all other terms are set to zero. (Stage 3) Post-regularization polishing. After selecting the subset of active terms, we again train without imposing Lp regularization, i.e. we set λp = 0 again. Hence, the parameters of terms in the resulting subset are truly optimized for minimizing the imbalance of equilibrium. Taken together, these three stages reduce sensitivity to model initialization and ensure that the parameter values of the selected model truly minimize the imbalance of equilibrium. Further details and values for the training and validation data generation along with model and training hyperparameters are presented in Tab. B.1.

Table B.1: Parameters and hyperparameters used for the data generation, validation, and model training.

Parameter Training specimen: Number of nodes Number of reaction force constraints Number of data snapshots Loading parameter Validation specimen: Number of nodes CANN hyperparameters: Number of power functions Number of activation functions Activation functions

Notation

Value

nn nβ nt δ

1 464 6 10 {0.1 × t : t = 1, . . . , nt }

20 082

nf ng −

2 2 •, exp (•)

Pseudo-invariants (isotropic cases)

−

Pseudo-invariants (anisotropic cases)

−

3/2 2 I˜1 − 3, I˜2 − 33/2 , [J − 1] D E2 3/2 2 I˜1 − 3, I˜2 − 33/2 , [J − 1] , I˜4 − 1

λint λext λp − − −

1 1 0 Adam 4000 0.025

λint λext λp p − − − −

1 1 0.001 0.25 Adam 4000 0.025 1 × 10−4

λint λext λp − − −

1 1 0 Adam 4000 0.005

Training hyperparameters: Pre-regularization: Internal force balance weight External force balance weight Regularization weight Optimizer Epochs Learning rate Regularization and subset selection: Internal force balance weight External force balance weight Regularization weight Regularization exponent Optimizer Epochs Learning rate Subset cut-off threshold Post-regularization polishing: Internal force balance weight External force balance weight Regularization weight Optimizer Epochs Learning rate

References [1] J. Ghaboussi, J. H. Garrett, X. Wu, Knowledge-based modeling of material behavior with neural networks, Journal of Engineering Mechanics 117 132–153. doi:10.1061/(ASCE)0733-9399(1991)117: 1(132). URL https://ascelibrary.org/doi/10.1061/%28ASCE%290733-9399%281991%29117%3A1%28132% 29 [2] F. Masi, I. Stefanou, P. Vannucci, V. Maffi-Berthier, Thermodynamics-based artificial neural networks for constitutive modeling, Journal of the Mechanics and Physics of Solids 147 (2021) 104277. doi: 10.1016/j.jmps.2020.104277. URL https://doi.org/10.1016/j.jmps.2020.104277 [3] N. N. Vlassis, W. Sun, Sobolev training of thermodynamic-informed neural networks for interpretable elasto-plasticity models with level set hardening, Computer Methods in Applied Mechanics and Engineering 377 (2021) 113695. doi:10.1016/j.cma.2021.113695. URL https://doi.org/10.1016/j.cma.2021.113695 [4] V. Tac, V. D. Sree, M. K. Rausch, A. B. Tepole, Data-driven modeling of the mechanical behavior of anisotropic soft biological tissue, Engineering with Computers 38 (2022) 4167–4182. doi:10.1007/ s00366-022-01733-3. URL http://dx.doi.org/10.1007/s00366-022-01733-3 [5] V. Tac, F. Sahli Costabal, A. B. Tepole, Data-driven tissue mechanics with polyconvex neural ordinary differential equations, Computer Methods in Applied Mechanics and Engineering 398 (2022) 115248. doi:10.1016/j.cma.2022.115248. URL http://dx.doi.org/10.1016/j.cma.2022.115248 [6] L. Linden, D. K. Klein, K. A. Kalina, J. Brummund, O. Weeger, M. Kästner, Neural networks meet hyperelasticity: A guide to enforcing physics, Journal of the Mechanics and Physics of Solids 179 (2023) 105363. doi:10.1016/j.jmps.2023.105363. URL http://dx.doi.org/10.1016/j.jmps.2023.105363 [7] J. Dornheim, L. Morand, H. J. Nallani, D. Helm, Neural networks for constitutive modeling: From universal function approximators to advanced models and the integration of physics, Archives of Computational Methods in Engineering 31 (2024) 1097–1127. doi:10.1007/s11831-023-10009-y. URL https://doi.org/10.1007/s11831-023-10009-y [8] J. N. Fuhg, G. Anantha Padmanabha, N. Bouklas, B. Bahmani, W. Sun, N. N. Vlassis, M. Flaschel, P. Carrara, L. De Lorenzis, A review on data-driven constitutive laws for solids, Archives of Computational Methods in Engineering (2024) 1–43. [9] V. Tac, E. Kuhl, A. B. Tepole, Data-driven continuum damage mechanics with built-in physics, Extreme Mechanics Letters 71 (2024) 102220. doi:10.1016/j.eml.2024.102220. URL https://doi.org/10.1016/j.eml.2024.102220 [10] V. Taç, M. K. Rausch, I. Bilionis, F. Sahli Costabal, A. B. Tepole, Generative hyperelasticity with physics-informed probabilistic diffusion fields, Engineering with Computers 41 (2025) 51–69. doi: 10.1007/s00366-024-01984-2. URL https://doi.org/10.1007/s00366-024-01984-2 [11] B. Alheit, M. Peirlinck, S. Kumar, Commet: Orders-of-magnitude speed-up in finite element method via batch-vectorized neural constitutive updates, Computer Methods in Applied Mechanics and Engineering 452 (2026) 118728. doi:10.1016/j.cma.2026.118728. URL https://doi.org/10.1016/j.cma.2026.118728 28

[12] A. A. Jadoon, K. A. Kalina, M. K. Rausch, R. Jones, J. N. Fu hg, Inverse design of anisotropic microstructures using physics-augmented neural networks, Journal of the Mechanics and Physics of Solids 203 (2025) 106161. doi:10.1016/j.jmps.2025.106161. URL https://doi.org/10.1016/j.jmps.2025.106161 [13] M. Flaschel, D. Martonová, C. Veil, E. Kuhl, Material fingerprinting: A shortcut to material model discovery without solving optimization problems, Computer Methods in Applied Mechanics and Engineering 450 (2026) 118573. doi:10.1016/j.cma.2025.118573. URL https://doi.org/10.1016/j.cma.2025.118573 [14] D. Martonová, E. Kuhl, M. Flaschel, Material fingerprinting for rapid discovery of hyperelastic models: First experimental validation, Journal of the Mechanics and Physics of Solids 208 (2026) 106463. doi: 10.1016/j.jmps.2025.106463. URL https://doi.org/10.1016/j.jmps.2025.106463 [15] M. Flaschel, M. A. Moreno-Mateos, S. Wiesheier, P. Steinmann, E. Kuhl, Unsupervised material fingerprinting: Ultra-fast hyperelastic model discovery from full-field experimental measurements (1 2026). arXiv:2601.14965, doi:10.48550/arXiv.2601.14965. URL http://arxiv.org/abs/2601.14965 [16] D. K. Klein, M. Fernández, R. J. Martin, P. Neff, O. Weeger, Polyconvex anisotropic hyperelasticity with neural networks, Journal of the Mechanics and Physics of Solids 159 (2022) 104703. doi:10.1016/ j.jmps.2021.104703. URL https://doi.org/10.1016/j.jmps.2021.104703 [17] E. Magaña, S. Pezzuto, F. Sahli Costabal, Ensemble learning of the atrial fibre orientation with physicsinformed neural networks, The Journal of Physiology n/a (9 2025). doi:10.1113/jp288001. URL https://doi.org/10.1113/jp288001 [18] K. Linka, M. Hillgärtner, K. P. Abdolazizi, R. C. Aydin, M. Itskov, C. J. Cyron, Constitutive artificial neural networks: A fast and general approach to predictive data-driven constitutive modeling by deep learning, Journal of Computational Physics 429 (2021) 110010. doi:10.1016/j.jcp.2020.110010. URL http://dx.doi.org/10.1016/j.jcp.2020.110010 [19] K. Linka, E. Kuhl, A new family of constitutive artificial neural networks towards automated model discovery, Computer Methods in Applied Mechanics and Engineering 403 (2023) 115731. doi:10.1016/ j.cma.2022.115731. URL https://doi.org/10.1016/j.cma.2022.115731 [20] P. Thakolkaran, A. Joshi, Y. Zheng, M. Flaschel, L. De Lorenzis, S. Kumar, Nn-euclid: Deep-learning hyperelasticity without stress data, Journal of the Mechanics and Physics of Solids 169 (2022) 105076. doi:10.1016/j.jmps.2022.105076. URL http://dx.doi.org/10.1016/j.jmps.2022.105076 [21] P. Thakolkaran, Y. Guo, S. Saini, M. Peirlinck, B. Alheit, S. Kumar, Can kan cans? input-convex kolmogorov-arnold networks (kans) as hyperelastic constitutive artificial neural networks (cans), Computer Methods in Applied Mechanics and Engineering 443 (2025) 118089. [22] K. P. Abdolazizi, R. C. Aydin, C. J. Cyron, K. Linka, Constitutive kolmogorov–arnold networks (ckans): Combining accuracy and interpretability in data-driven material modeling, Journal of the Mechanics and Physics of Solids 203 (2025) 106212. doi:10.1016/j.jmps.2025.106212. URL http://dx.doi.org/10.1016/j.jmps.2025.106212 [23] K. Linka, S. R. St. Pierre, E. Kuhl, Automated model discovery for human brain using constitutive artificial neural networks, Acta Biomaterialia 160 (2023) 134–151. doi:10.1016/j.actbio.2023.01. 055. URL https://doi.org/10.1016/j.actbio.2023.01.055 29

[24] M. Peirlinck, K. Linka, J. A. Hurtado, E. Kuhl, On automated model discovery and a universal material subroutine for hyperelastic materials, Computer Methods in Applied Mechanics and Engineering 418 (2024) 116534. doi:10.1016/j.cma.2023.116534. URL https://doi.org/10.1016/j.cma.2023.116534 [25] S. R. St. Pierre, K. Linka, E. Kuhl, Principal-stretch-based constitutive neural networks autonomously discover a subclass of ogden models for human brain tissue, Brain Multiphysics 4 (2023) 100066. doi: 10.1016/j.brain.2023.100066. URL https://doi.org/10.1016/j.brain.2023.100066 [26] K. Linka, A. Buganza Tepole, G. A. Holzapfel, E. Kuhl, Automated model discovery for skin: Discovering the best model, data, and experiment, Computer Methods in Applied Mechanics and Engineering 410 (2023) 116007. doi:10.1016/j.cma.2023.116007. URL https://doi.org/10.1016/j.cma.2023.116007 [27] S. R. St. Pierre, D. Rajasekharan, E. C. Darwin, K. Linka, M. E. Levenston, E. Kuhl, Discovering the mechanics of artificial and real meat, Computer Methods in Applied Mechanics and Engineering 415 (2023) 116236. doi:10.1016/j.cma.2023.116236. URL https://doi.org/10.1016/j.cma.2023.116236 [28] M. Peirlinck, K. Linka, J. A. Hurtado, G. A. Holzapfel, E. Kuhl, Democratizing biomedical simulation through automated model discovery and a universal material subroutine, Computational mechanics (2024) 1–21. [29] T. Vervenne, M. Peirlinck, N. Famaey, E. Kuhl, Constitutive neural networks for main pulmonary arteries: discovering the undiscovered, Biomechanics and Modeling in Mechanobiology 24 (2025) 615– 634. doi:10.1007/s10237-025-01930-1. URL https://doi.org/10.1007/s10237-025-01930-1 [30] D. Martonová, M. Peirlinck, K. Linka, G. A. Holzapfel, S. Leyendecker, E. Kuhl, Automated model discovery for human cardiac tissue: Discovering the best model and parameters, Computer Methods in Applied Mechanics and Engineering 428 (2024) 117078. doi:10.1016/j.cma.2024.117078. URL https://doi.org/10.1016/j.cma.2024.117078 [31] M. Peirlinck, K. Linka, E. Kuhl, Atrial constitutive neural networks (4 2025). arXiv:2504.02748, doi:10.48550/arXiv.2504.02748. URL http://arxiv.org/abs/2504.02748 [32] L. M. Wang, K. Linka, E. Kuhl, Automated model discovery for muscle using constitutive recurrent neural networks, Journal of the Mechanical Behavior of Biomedical Materials 145 (2023) 106021. doi: 10.1016/j.jmbbm.2023.106021. URL https://doi.org/10.1016/j.jmbbm.2023.106021 [33] K. P. Abdolazizi, K. Linka, C. J. Cyron, Viscoelastic constitutive artificial neural networks (vcanns) –a framework for data-driven anisotropic nonlinear finite viscoelasticity, Journal of Computational Physics 499 (2024) 112704. doi:10.1016/j.jcp.2023.112704. URL https://doi.org/10.1016/j.jcp.2023.112704 [34] H. Holthusen, T. Brepols, K. Linka, E. Kuhl, Automated model discovery for tensional homeostasis: Constitutive machine learning in growth and remodeling, Computers in Biology and Medicine 186 (2025) 109691. doi:10.1016/j.compbiomed.2025.109691. URL https://doi.org/10.1016/j.compbiomed.2025.109691 [35] K. Linka, G. A. Holzapfel, E. Kuhl, Discovering uncertainty: Bayesian constitutive artificial neural networks, Computer Methods in Applied Mechanics and Engineering 433 (2025) 117517. doi:10. 1016/j.cma.2024.117517. URL https://doi.org/10.1016/j.cma.2024.117517 30

[36] M. Peirlinck, J. A. Hurtado, M. K. Rausch, A. B. Tepole, E. Kuhl, A universal material model subroutine for soft matter systems, Engineering with Computers (9 2024). doi:10.1007/s00366-024-02031-w. URL https://doi.org/10.1007/s00366-024-02031-w [37] G. A. Holzapfel, T. C. Gasser, R. W. Ogden, A new constitutive framework for arterial wall mechanics and a comparative study of material models, Journal of Elasticity 61 (2000) 1–48. doi:10.1023/a: 1010835316564. URL https://doi.org/10.1023/a:1010835316564 [38] T. C. Gasser, R. W. Ogden, G. A. Holzapfel, Hyperelastic modelling of arterial layers with distributed collagen fibre orientations, Journal of The Royal Society Interface 3 (2006) 15–35. doi:10.1098/rsif. 2005.0073. URL http://dx.doi.org/10.1098/rsif.2005.0073 [39] G. A. Holzapfel, R. W. Ogden, Constitutive modelling of passive myocardium: a structurally based framework for material characterization, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 367 (2009) 3445–3475. doi:10.1098/rsta.2009.0091. URL https://doi.org/10.1098/rsta.2009.0091 [40] R. Avazmohammadi, J. S. Soares, D. S. Li, S. S. Raut, R. C. Gorman, M. S. Sacks, A contemporary look at biomechanical models of myocardium, Annual Review of Biomedical Engineering 21 (2019) 417–442. doi:10.1146/annurev-bioeng-062117-121129. URL https://doi.org/10.1146/annurev-bioeng-062117-121129 [41] A. Aggarwal, L. T. Hudson, D. W. Laurence, C.-H. Lee, S. Pant, A bayesian constitutive model selection framework for biaxial mechanical testing of planar soft tissues: Application to porcine aortic valves, Journal of the Mechanical Behavior of Biomedical Materials 138 (2023) 105657. doi:10.1016/j.jmbbm. 2023.105657. URL http://dx.doi.org/10.1016/j.jmbbm.2023.105657 [42] R. P. Krijnen, A. Joshi, S. Kumar, M. Peirlinck, Unsupervised full-field bayesian inference of orthotropic hyperelasticity from a single biaxial test: a myocardial case study, Computer Methods in Applied Mechanics and Engineering 459 (2026) 119034. doi:10.1016/j.cma.2026.119034. URL https://doi.org/10.1016/j.cma.2026.119034 [43] N. Famaey, H. Fehervary, Y. Lafon, A. Akyildiz, S. Dreesen, K. Bruyère-Garnier, J.-M. Allain, M. Alloisio, A. Aparici-Gil, C. Catalano, F. Chassagne, S. Chokhandre, K. Crevits, H. Crielaard, E. Cunnane, C. Cunnane, K. De Leener, A. Desai, R. Driessen, A. Erdemir, M. Eskandari, S. Evans, C. Gasser, M. Gebhardt, B. Glasmacher, G. A. Holzapfel, M. Isasi, L. Jennings, S. Kurz, S. LealMarin, P. Lecomte, A. Morch, J. Mulvihill, F. Nemavhola, T. Pandelani, S. Pasta, E. Peña, B. Pierrat, H.-L. Ploeg, S. Polzer, M. Rausch, D. Schwarz, H. Screen, S. Sherifova, G. Sommer, S. Wang, D. Walsh, D. Yadav, T. Marchal, L. Geris, Community challenge towards consensus on characterization of biological tissue: C4bio’s first findings, Journal of Biomechanics 194 (2026) 113021. doi:10.1016/j.jbiomech.2025.113021. URL https://doi.org/10.1016/j.jbiomech.2025.113021 [44] F. Pierron, M. Grédiac, The virtual fields method: extracting constitutive mechanical parameters from full-field deformation measurements, Springer Science \& Business Media, 2012. [45] F. O. Kolawole, M. Peirlinck, T. E. Cork, M. Levenston, E. Kuhl, D. B. Ennis, Validating mri-derived myocardial stiffness estimates using in vitro synthetic heart models, Annals of Biomedical Engineering 51 (7) (2023) 1574–1587. doi:10.1007/s10439-023-03164-7. URL http://dx.doi.org/10.1007/s10439-023-03164-7

[46] S. Meng, A. A. K. Yousefi, S. Avril, Machine-learning-based virtual fields method: Application to anisotropic hyperelasticity, Computer Methods in Applied Mechanics and Engineering 434 (2025) 117580. doi:10.1016/j.cma.2024.117580. URL https://doi.org/10.1016/j.cma.2024.117580 [47] X. Navy, Z. Sheng, K. Kim, J. M. Cormack, Three-dimensional tissue strain measurement using a row–column array during biaxial testing of excised ventricular porcine myocardium, Ultrasound in Medicine & Biology 51 (2025) 1622–1626. doi:10.1016/j.ultrasmedbio.2025.05.007. URL https://doi.org/10.1016/j.ultrasmedbio.2025.05.007 [48] M. Grédiac, Principe des travaux virtuels et identification, Comptes rendus de l'Acad\'emie dessciences. S\'erie 2, M\'ecanique, Physique, Chimie, Sciences de l'univers, Sciences de la Terre 309 (1) (1989) 1–5. [49] M. Grédiac, F. Pierron, S. Avril, E. Toussaint, The virtual fields method for extracting constitutive parameters from full-field measurements: a review, Strain 42 (2006) 233–253. doi:10.1111/j.1475-1305. 2006.tb01504.x. URL https://doi.org/10.1111/j.1475-1305.2006.tb01504.x [50] S. Avril, Recent advances in the virtual fields method for evaluating and identifying tissue biomechanical properties and constitutive parameters, Journal of Biomechanical Engineering 148 (5 2026). doi: 10.1115/1.4071211. URL https://doi.org/10.1115/1.4071211 [51] D. Claire, F. Hild, S. Roux, A finite element formulation to identify damage fields: the equilibrium gap method, International Journal for Numerical Methods in Engineering 61 (2004) 189–208. doi: 10.1002/nme.1057. URL https://doi.org/10.1002/nme.1057 [52] M. Flaschel, S. Kumar, L. De Lorenzis, Unsupervised discovery of interpretable hyperelastic constitutive laws, Computer Methods in Applied Mechanics and Engineering 381 (2021) 113852. doi:10.1016/j. cma.2021.113852. URL http://dx.doi.org/10.1016/j.cma.2021.113852 [53] H. Demiray, A note on the elasticity of soft biological tissues, Journal of Biomechanics 5 (3) (1972) 309–311. doi:10.1016/0021-9290(72)90047-4. URL http://dx.doi.org/10.1016/0021-9290(72)90047-4 [54] J. D. Humphrey, F. C. P. Yin, On constitutive relations and finite deformations of passive cardiac tissue: I. a pseudostrain-energy function, Journal of Biomechanical Engineering 109 (1987) 298–304. doi:10.1115/1.3138684. URL https://doi.org/10.1115/1.3138684 [55] J. M. Guccione, A. D. McCulloch, L. K. Waldman, Passive material properties of intact ventricular myocardium determined from a cylindrical model, Journal of Biomechanical Engineering 113 (1991) 42–55. doi:10.1115/1.2894084. URL https://doi.org/10.1115/1.2894084 [56] G. A. Holzapfel, Nonlinear solid mechanics: a continuum approach for engineering, Wiley, 2000. [57] M. Flaschel, S. Kumar, L. De Lorenzis, Discovering plasticity models without stress data, npj Computational Materials 8 (2022) 1–10. doi:10.1038/s41524-022-00752-4. URL http://dx.doi.org/10.1038/s41524-022-00752-4 [58] H. Xu, M. Flaschel, L. De Lorenzis, Discovering non-associated pressure-sensitive plasticity models with euclid, Advanced Modeling and Simulation in Engineering Sciences 12 (2025) 1. doi:10.1186/ s40323-024-00281-3. URL https://doi.org/10.1186/s40323-024-00281-3 32

[59] E. Marino, M. Flaschel, S. Kumar, L. De Lorenzis, Automated identification of linear viscoelastic constitutive laws with euclid, Mechanics of Materials 181 (2023) 104643. doi:10.1016/j.mechmat. 2023.104643. URL http://dx.doi.org/10.1016/j.mechmat.2023.104643 [60] M. Flaschel, S. Kumar, L. De Lorenzis, Automated discovery of generalized standard material models with euclid, Computer Methods in Applied Mechanics and Engineering 405 (2023) 115867. doi:10. 1016/j.cma.2022.115867. URL https://www.sciencedirect.com/science/article/pii/S0045782522008234 [61] A. Joshi, P. Thakolkaran, Y. Zheng, M. Escande, M. Flaschel, L. De Lorenzis, S. Kumar, Bayesianeuclid: Discovering hyperelastic material laws with uncertainties, Computer Methods in Applied Mechanics and Engineering 398 (2022) 115225. doi:10.1016/j.cma.2022.115225. URL http://dx.doi.org/10.1016/j.cma.2022.115225 [62] D. Martonová, S. Leyendecker, G. A. Holzapfel, E. Kuhl, Discovering dispersion: How robust is automated model discovery for human myocardial tissue?, Biomechanics and Modeling in Mechanobiology 24 (2025) 2023–2037. doi:10.1007/s10237-025-02005-x. URL https://doi.org/10.1007/s10237-025-02005-x [63] D. Martonová, A. Goriely, E. Kuhl, Generalized invariants meet constitutive neural networks: A novel framework for hyperelastic materials, Journal of the Mechanics and Physics of Solids 206 (2026) 106352. doi:10.1016/j.jmps.2025.106352. URL https://doi.org/10.1016/j.jmps.2025.106352 [64] J. M. Ball, Convexity conditions and existence theorems in nonlinear elasticity, Archive for Rational Mechanics and Analysis 63 (1976) 337–403. doi:10.1007/bf00279992. URL https://doi.org/10.1007/bf00279992 [65] D. Martonová, A. Goriely, E. Kuhl, Generalized invariants meet constitutive neural networks: A novel framework for hyperelastic materials (8 2025). [66] J. Schröder, P. Neff, Invariant formulation of hyperelastic transverse isotropy based on polyconvex free energy functions, International Journal of Solids and Structures 40 (2003) 401–445. doi:10.1016/ s0020-7683(02)00458-4. URL http://dx.doi.org/10.1016/s0020-7683(02)00458-4 [67] D. Balzani, P. Neff, J. Schröder, G. Holzapfel, A polyconvex framework for soft biological tissues. adjustment to experimental data, International Journal of Solids and Structures 43 (2006) 6052–6070. doi:10.1016/j.ijsolstr.2005.07.048. URL http://dx.doi.org/10.1016/j.ijsolstr.2005.07.048 [68] l. E. Frank, J. H. Friedman, A statistical view of some chemometrics regression tools, Technometrics 35 (1993) 109–135. doi:10.1080/00401706.1993.10485033. URL https://doi.org/10.1080/00401706.1993.10485033 [69] A. E. Hoerl, R. W. Kennard, Ridge regression: Biased estimation for nonorthogonal problems, Technometrics 12 (1970) 55–67. doi:10.1080/00401706.1970.10488634. URL https://doi.org/10.1080/00401706.1970.10488634 [70] J. A. McCulloch, S. R. St. Pierre, K. Linka, E. Kuhl, On sparse regression, lp-regularization, and automated model discovery, International Journal for Numerical Methods in Engineering 125 (14) (2024) e7481.

[71] R. S. Rivlin, Large elastic deformations of isotropic materials iv. further developments of the general theory, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 241 (1948) 379–397. doi:10.1098/rsta.1948.0024. URL https://doi.org/10.1098/rsta.1948.0024 [72] A. Isihara, N. Hashitsume, M. Tatibana, Statistical theory of rubber-like elasticity. iv. (two-dimensional stretching), The Journal of Chemical Physics 19 (1951) 1508–1512. doi:10.1063/1.1748111. URL https://doi.org/10.1063/1.1748111 [73] E. M. Arruda, M. C. Boyce, A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials, Journal of the Mechanics and Physics of Solids 41 (1993) 389–412. doi: 10.1016/0022-5096(93)90013-6. URL https://doi.org/10.1016/0022-5096(93)90013-6 [74] A. N. Gent, A. G. Thomas, Forms for the stored (strain) energy function for vulcanized rubber, Journal of Polymer Science 28 (1958) 625–628. doi:10.1002/pol.1958.1202811814. URL https://doi.org/10.1002/pol.1958.1202811814 [75] R. W. Ogden, Large deformation isotropic elasticity –on the correlation of theory and experiment for incompressible rubberlike solids, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 326 (1972) 565–584. doi:10.1098/rspa.1972.0026. URL https://doi.org/10.1098/rspa.1972.0026 [76] X. Ning, Q. Zhu, Y. Lanir, S. S. Margulies, A transversely isotropic viscoelastic constitutive equation for brainstem undergoing finite deformation, Journal of Biomechanical Engineering 128 (2006) 925–933. doi:10.1115/1.2354208. URL https://doi.org/10.1115/1.2354208 [77] D. F. Meaney, Relationship between structural modeling and hyperelastic material behavior: application to cns white matter, Biomechanics and Modeling in Mechanobiology 1 (2003) 279–293. doi:10.1007/ s10237-002-0020-1. URL https://doi.org/10.1007/s10237-002-0020-1 [78] J. Merodio, R. Ogden, Mechanical response of fiber-reinforced incompressible non-linearly elastic solids, International Journal of Non-Linear Mechanics 40 (2005) 213–227. doi:10.1016/j.ijnonlinmec. 2004.05.003. URL https://doi.org/10.1016/j.ijnonlinmec.2004.05.003 [79] B. Amos, L. Xu, J. Z. Kolter, Input convex neural networks, in: Proceedings of the 34th International Conference on Machine Learning, PMLR, 2017, pp. 146–155. URL https://proceedings.mlr.press/v70/amos17b.html [80] Z. Liu, Y. Wang, S. Vaidya, F. Ruehle, J. Halverson, M. Soljacić, T. Y. Hou, M. Tegmark, Kan: Kolmogorov-arnold networks (2 2025). arXiv:2404.19756, doi:10.48550/arXiv.2404.19756. URL http://arxiv.org/abs/2404.19756 [81] J. Morrey, B. Charles, Quasi-convexity and the lower semicontinuity of multiple integrals, Pacific Journal of Mathematics 2 (1952) 25–53.

Record · ID 271840 · SHA-256 e2a5e5c00a7a60e7
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.