1 Introduction↩︎

Finding suitable mathematical descriptions of mechanical material behavior is central to simulation-based engineering. In many material modeling frameworks, especially those grounded in thermodynamics [1]–[4], the material behavior is fully characterized by a set of constitutive functions, such as the strain-energy density potential or the dissipation density potential. The objective of mechanical material characterization is to infer these constitutive functions from experimental data, which typically include reaction forces and displacement fields of deformed specimens over time. Both the experimental observations and the underlying material behavior can naturally be interpreted as elements of infinite-dimensional function spaces. Consequently, material characterization can be formulated as the problem of identifying a mapping between such infinite-dimensional spaces. This perspective motivates the use of neural operator learning [5], which aims to learn mappings between infinite-dimensional function spaces. In this work, we take a first step in this direction by developing a neural operator architecture that maps full-field displacement data and reaction forces directly to the strain-energy density function of hyperelastic materials. Once trained, the neural operator characterizes a new material via a single forward pass, without solving any optimization problem which makes it significantly faster than optimization-based approaches. Our proposed framework is graphically illustrated in 1. Before presenting the framework in more detail, we briefly review recent advances in neural operator learning and machine learning-based material modeling.

Conventional neural networks approximate mappings between finite-dimensional spaces. For instance, a standard feed-forward neural network can be viewed as a mapping \(\text{NN} : \mathbb{R}^n \rightarrow \mathbb{R}^m\), where \(\mathbb{R}^n\) and \(\mathbb{R}^m\) are Euclidean spaces with dimensions \(n\) and \(m\), respectively. However, many problems in physics and engineering involve mappings between infinite-dimensional function spaces. These mappings are commonly called operators. A prominent example is the solution operator of a partial differential equation (PDE), which maps input functions, such as spatially varying parameters, geometric descriptions, boundary conditions, and initial conditions, to solution functions defined over space and time. Recognizing that conventional neural networks are inherently limited in approximating mappings between infinite-dimensional spaces, [5] proposed in their seminal work a neural network architecture specifically designed for this purpose. In contrast to standard networks, their approach aims to learn operators of the form \(\mathrm{NO} : \mathcal{U} \rightarrow \mathcal{W}\), where \(\mathcal{U}\) and \(\mathcal{W}\) denote infinite-dimensional function spaces. The proposed architecture consists of two coupled neural networks: one that encodes the input function in the space \(\mathcal{U}\), and another that represents the dependence on the evaluation point in the output space \(\mathcal{W}\). By combining these two components, the network effectively parameterizes a nonlinear operator between function spaces. [5] further showed that this architecture possesses a universal approximation property. It can approximate a broad class of continuous nonlinear operators between infinite-dimensional spaces to arbitrary accuracy, provided that the network complexity is sufficiently increased, see also [6] for earlier results on the approximation of functionals. The framework introduced by [5] was later adopted and further developed by [7], who popularized it under the name DeepONet. [7] established the terminology branch net and trunk net for the two aforementioned subnetworks in the architecture. DeepONets have since found widespread use across scientific and engineering domains, including phase-field modeling [8], porous media and geology [9], fluid dynamics [10], and robotics [11]. Other widely recognized architectures for approximating mappings between infinite-dimensional function spaces have been proposed by [12], who proposed neural operator architectures based on the composition of integral operators with nonlinear activation functions. In subsequent work, the same authors developed a variant based on Fourier operators [13], leading to the formulation of the Fourier Neural Operator (FNO) whose universal approximation abilities are discussed by [14]. An overview of different choices for the integral operators in the network is provided by [15]. FNOs and related architectures have since been applied across a wide range of problems in physics and engineering [16]–[19]. A fundamental challenge in learning the solution operator of a PDE is achieving independence from both the underlying geometry and its discretization [20]. A popular approach to address this are graph-based neural operators [12], [21], [22], which are adaptations of Graph Neural Networks (GNN) [23]–[25] that are discretization-convergent. The neural operators discussed above can be trained either in a supervised fashion using labeled input-output pairs, or in a physics-informed manner by incorporating the PDE residual directly into the training loss [26]–[29].

Owing to their universal approximation capabilities, neural operators have become increasingly popular in computational solid mechanics. For instance, [30] proposed a DeepONet for crack path prediction in brittle materials. [31] use a DeepONet to predict stresses in elastoplastic structures. [32] proposed and compared three DeepONet architectures for learning the solution operator of phase field fracture problems. [33] use DeepONets for predicting the stress field in small strain plasticity problems. [34] use FNOs to predict the deformation of hyperelastic microstructures on the unit domain, which is then used to compute the homogenized material response. [35] investigate the use of neural operators in solid mechanics and draw a comparison to the Finite Element Method (FEM). [36] propose physics-informed FNOs for linear isotropic elasticity on the unit domain. [37] use DeepONets to approximate the operator that maps heterogeneous microstructures to the deformation fields in small strain elasticity on the unit domain. [38] use FNOs for predicting the solution of micromechanics problems in linear isotropic or anisotropic elasticity on the unit domain and investigate their relation to methods based on the Fast Fourier Transform (FFT). [39] propose physics-informed reduced-order operator learning for homogenizing hyperelastic microstructures. [40] propose different neural operator architectures for learning the solution operator of geometrically linear and nonlinear elasticity, with a particular focus on arbitrary boundary conditions.

Most existing work on neural operators addresses forward problems and focuses on learning solution operators of PDEs, which map input quantities such as parameters, geometries, boundary conditions, and initial conditions to solution fields. In contrast, the use of neural operators for inverse problems has received comparatively little attention in the literature. Exceptions are the works by [41], who discuss the use of neural operators for PDE-constrained inverse problems, by [42], who leverage neural operators for solving inverse medium problems such as Magnetic Resonance Imaging (MRI), or by [43], who propose a special neural operator architecture for tackling the ill-posedness of inverse problems. In the field of computational solid mechanics, only a few works apply neural operator learning to inverse problems. For example, [44] use recurrent neural operators for learning plasticity models. [29] use physics-informed neural operators for the inference of thermoelectric properties from data. [45] leverage DeepONets for the inverse design of metamaterials. [46] use FNO-inspired architectures to learn the material response of elastoplastic materials and materials undergoing damage. [47] propose a time-resolution-independent architecture for learning the material response of anisotropic and temperature-dependent plasticity.

Prior work has explored the use of neural operators for learning constitutive models [44], [47]–[49]. However, in these works, the neural operators’ weights are trained to minimize the mismatch between the model predictions and the data for a specific material. For a new experimental dataset or a different material, the neural operators must be retrained. To the best of our knowledge, no existing study employs neural operators to directly map experimental fields to constitutive functions, enabling the trained neural operator to be repeatedly used for material model inference without retraining. Hence, we propose a neural operator framework that directly maps experimentally measured functions, such as displacement fields and reaction forces over time, to the material’s characteristic thermodynamic potentials. In this work, we focus on hyperelastic materials whose mechanical response is fully characterized by the strain-energy density function. Hence, our neural operators directly map from the spaces of the experimental data to the space of admissible strain-energy density functions. The latter space is constrained to include only mathematically and physically meaningful functions that fulfill fundamental physical laws, such as invariance to the observer, a stress-free undeformed configuration, the balance of angular momentum, and continuous differentiability, along with other desired properties such as material symmetry and polyconvexity. Imposing physical constraints on the material model reduces the design space, makes models less data-hungry, simplifies training, improves extrapolation, and prevents unphysical material responses [50]–[64]. Here, we specifically draw inspiration from two modeling approaches, Physics-Augmented Neural Networks (PANNs) [52], [54], [58] and Constitutive Artificial Neural Networks (CANNs) [51], [55]. PANNs are flexible neural network architectures for constitutive modeling that satisfy the fundamental physical constraints. We emphasize that physics-augmented should not be confused with physics-informed, which typically refers to loss functions informed by governing physical equations rather than labeled data pairs [26], [65]. CANNs, on the other hand, incorporate well-established constitutive modeling features into the neural network architecture, which reduces the number of weights and makes the networks more interpretable, generally at the cost of reduced flexibility.

Beyond physical admissibility, our neural operators offer a significant computational advantage over traditional methods for material characterization. Traditional approaches rely on optimization problems in which material model parameters are calibrated to minimize the discrepancy between model predictions and experimental data. This process can be computationally expensive, particularly when each model evaluation requires solving nonlinear equations with the Newton-Raphson method. In contrast, our proposed neural operators require only a single forward pass to infer the material model from experimental data, which is significantly faster than traditional optimization-based characterization methods. In this context, we note that our proposed neural operator framework shares methodological similarities with approaches that directly map experimental data to the parameters of a priori chosen material models [66]. However, in contrast to such works, our framework significantly reduces a priori assumptions about material behavior. Rather than fixing the functional form of the strain-energy density function, it can predict a variety of functional forms depending on the available experimental data. In addition, to the best of our knowledge, no learning framework has been presented so far that directly predicts the material model for given displacement fields in a discretization-independent manner.

We further note that the issue of computationally expensive inverse problems has been addressed in our recently proposed Material Fingerprinting framework [67]–[70], in which material models are selected from a precomputed database of material responses for a given experimental dataset. By replacing the continuous optimization problem with a fast pattern recognition algorithm, Material Fingerprinting enables significantly faster material characterization. While Material Fingerprinting precomputes and stores a database for use with a pattern recognition algorithm, our proposed neural operators can be interpreted as compressing the information contained in such a database into a neural operator that directly maps experimental data to constitutive functions.

Methods for mechanical material characterization can be broadly classified into those informed by labeled stress-strain data pairs from simple experiments with homogeneous deformation fields, such as uniaxial tension tests, and those informed by full-field displacement data and net reaction forces from complexly shaped geometries [71]–[81]. The former approach can yield ill-posed inverse problems, as different material models can produce identical stress-strain responses in the high-dimensional stress-strain space [82]–[84]. The latter approach offers the advantage that the space of deformations imposed on the material is enriched, which can reduce the ill-posedness of the inverse problem. Hence, our proposed neural operators follow this approach. However, the framework can be readily adapted to leverage stress-strain data pairs from experiments with homogeneous deformation fields instead.

We describe our proposed neural operator framework and its training in 2. In 3, we evaluate the predictive performance of the neural operators on unseen test data, noisy data, data with missing information, data from different spatial discretizations, and data from geometries of varying sizes. This is followed by a discussion of the ill-posedness of the infinite-dimensional learning problem in 4, where we show numerically that the restricted output space of the proposed neural operators provides sufficient regularization of the inverse problem. We conclude in 5 by summarizing the main findings and outlining directions for future work.

Notation: Parentheses \((\square)\) denote tuples, open intervals, or function dependencies. Square brackets \([\square]\) denote closed intervals or the order of computation in equations. Curly brackets \(\{\square\}\) denote sets. Subscripts and superscripts written in upright font do not denote indexing. Tensors and matrices may appear in compact or index notation. In compact notation, vectors (first-order tensors) and second-order tensors are denoted by bold letters, e.g., \(\boldsymbol{C}\). When appearing in index notation, the order of the tensor equals the number of indices, e.g., \(C_{ij}\). Unless stated otherwise, the Einstein convention for summation over repeated indices is used in equations appearing in index notation, e.g., \(a_i b_i = \sum_i a_i b_i\). Inner products are denoted by \(\cdot\), e.g., \(\boldsymbol{a}\cdot \boldsymbol{b}= a_i b_i\), \(\boldsymbol{A}\cdot \boldsymbol{B}= A_{ij} B_{ij}\), and outer products by \(\otimes\), e.g., \([\boldsymbol{a}\otimes \boldsymbol{b}]_{ij} = a_i b_j\). If no operation is indicated between two tensors, the juxtaposition implies tensor contraction, e.g., \([\boldsymbol{A}\boldsymbol{B}]_{ij}=A_{ik}B_{kj}\). The second-order identity is denoted by \(\boldsymbol{I}\). The transpose of a tensor is denoted by \(\square^{\top}\). The trace of a tensor is denoted by \(\text{tr}(\square)\), i.e., \(\text{tr}(\boldsymbol{C}) = C_{ii}\). As we present a purely numerical study, we do not specify any physical units throughout this work.

2 Neural Operators↩︎

2.1 Function spaces↩︎

The general objective of material characterization based on full field data is to infer the material model from available displacement and reaction force data by solving an inverse problem. Mathematically, the objective is to infer an element in the material model space \(\mathcal{W}\) for a given element in the data space \(\mathcal{D} = \mathcal{U} \times \mathcal{R}\), where \(\mathcal{U}\) is the space of measurable displacements and \(\mathcal{R}\) is the space of measurable reaction forces. Hence, in this work, we propose neural operators that take elements in \(\mathcal{D}\) as input and output an element in \(\mathcal{W}\). In the following, we introduce the function spaces \(\mathcal{U}\), \(\mathcal{R}\), and \(\mathcal{W}\) and discuss how they are constrained by the physical principles and modeling assumptions adopted in this work.

2.1.1 Material model space↩︎

Since we assume incompressible, isotropic hyperelasticity in this work, the mechanical behavior is solely characterized by the strain-energy density function \(W\), and the space of material models reduces to the space of admissible strain-energy density functions \(\mathcal{W}\). In this work, we restrict \(\mathcal{W}\) to functions that satisfy fundamental physical constraints while also exhibiting additional desirable properties. Accordingly, we refer to our neural operators as Physics-Augmented Neural Operators (PANO) or Constitutive Artificial Neural Operators (CANO), in alignment with the notion of Physics-Augmented Neural Networks (PANNs) [52], [54], [58], [60] and Constitutive Artificial Neural Networks (CANNs) [51], [55] for constitutive modeling. Constraining the output space of the neural operators to physically meaningful models ensures that the predicted material models do not violate physical constraints and can hence be readily deployed in forward simulations without convergence, stability, or robustness issues. Furthermore, the predicted material models are expected to generalize better to unseen deformations than their unconstrained counterparts. Finally, restricting the output space reduces the ill-posedness of the problem, as physically inadmissible models that could otherwise explain the data with similar accuracy are excluded a priori.

The deformation of a body is mathematically described by the mapping \(\boldsymbol{x}= \phi(\boldsymbol{X},t)\), where \(t\) denotes time, \(\boldsymbol{X}\) represents a material point in the undeformed reference configuration, and \(\boldsymbol{x}\) denotes its corresponding position in the deformed configuration. The deformation gradient is defined as the Jacobian of the deformation mapping \(\boldsymbol{F}= \nabla_{\boldsymbol{X}} \phi\). Its determinant \(J = \det \boldsymbol{F}\) characterizes local volume changes and must therefore remain strictly positive to ensure physically admissible deformations. For incompressible material behavior, it is \(J=1\). The deformation gradient describes both volume changes and shape changes. We separate these contributions with the Flory split \(\boldsymbol{F}=\boldsymbol{F}^{\text{vol}} \bar \boldsymbol{F}\), where \(\boldsymbol{F}^{\text{vol}} = J^\frac{1}{3} \boldsymbol{I}\) is the volumetric, shape-preserving part of the deformation gradient, and \(\bar \boldsymbol{F}= J^{-\frac{1}{3}} \boldsymbol{F}\) is the isochoric, volume-preserving part. In the incompressible case, it is \(\boldsymbol{F}^{\text{vol}}=\boldsymbol{I}\), and the isochoric deformation gradient therefore plays a central role in the constitutive model. Since the isochoric deformation gradient is not invariant to rigid body motions of the deformed configuration and hence does not constitute an objective measure of deformation, we introduce the isochoric right Cauchy-Green strain \(\bar \boldsymbol{C}=\bar \boldsymbol{F}^{\top}\bar \boldsymbol{F}\), and its principal invariants \(\bar I_1=\text{tr}(\bar \boldsymbol{C})\), \(\bar I_2=\frac{1}{2}[\text{tr}(\bar \boldsymbol{C})^2 - \text{tr}(\bar \boldsymbol{C}^2)]\), and \(\bar I_3 = 1\). We further introduce the lifted polyconvex invariants \(\bar I_1^* = \bar I_1 - 3\) and \(\bar I_2^* = \bar I_2^{3/2} - 3^{3/2}\), where the motivation for the exponent \(\frac{3}{2}\) is detailed later in this section. We denote a tuple of these invariants by \(\bar\boldsymbol{I}^*=(\bar I_1^*,\bar I_2^*)\). For the region of possible values taken by these invariants, \(\bar\boldsymbol{I}^*\in\mathcal{I}_{\text{adm}}\subset \mathbb{R}_{\geq 0}^2\), we refer to 7.

To model incompressible hyperelastic material behavior, we consider strain-energy density functions of the form \(W(\boldsymbol{F}) = \bar W(\boldsymbol{F}) - p[J-1]\), where \(\bar W\) is the isochoric strain-energy and \(p\) is a Lagrange multiplier enforcing the incompressibility constraint. The stress response derives directly from the strain-energy density \[\boldsymbol{P} = \frac{\partial W}{\partial \boldsymbol{F}} = \frac{\partial \bar{W}}{\partial \boldsymbol{F}} - p \boldsymbol{F}^{-\top}.\] To ensure that the material model and its resulting stress-strain response satisfy fundamental physical constraints while also exhibiting additional desirable properties, we assume that the isochoric strain-energy can be expressed as a function of the two previously defined invariants \(\bar W(\bar\boldsymbol{I}^*)\), and we constrain the space of isochoric strain-energy functions to \[\label{eq:isochoric32strain-energy32density32space} \mathcal{W} = \left\{ \bar W \in C^2\left(\mathcal{I}_{\text{adm}}; \mathbb{R}_{\geq 0}\right) \;\middle| \; \bar{W} \;\text{separable, monotone, convex}, \bar{W}(0, 0) = 0 \right\},\tag{1}\] where the domain of the functions \(\bar W\) is restricted by the admissible region \(\mathcal{I}_{\text{adm}}\) of the invariants, see 7.

By restricting the output space of the neural operators to the constrained space defined above, we ensure that the material model and its resulting stress-strain response fulfill the following properties by construction:

1. Thermodynamic consistency: By defining the stress response as \(\boldsymbol{P}={\partial W}/{\partial \boldsymbol{F}}\), the material response is thermodynamically consistent by construction [85].

2. Objectivity: Since we are confined to strain-energy density functions that can be expressed solely as functions of the right Cauchy-Green tensor, which is invariant under rigid body motions of the deformed configuration, the material behavior is objective by construction, and the principle of material frame indifference is fulfilled. Specifically, it is \(W(\boldsymbol{F}) = W(\boldsymbol{Q}\boldsymbol{F})\) for all rotation tensors in the special orthogonal group \(\boldsymbol{Q}\in \text{SO}(3)\).

3. Balance of angular momentum: The balance of angular momentum states that the Cauchy stress tensor must be symmetric \(\boldsymbol{\sigma}= \boldsymbol{\sigma}^{\top}\). The Cauchy stress is given by \[\boldsymbol{\sigma}= \frac{1}{J} \boldsymbol{P}\boldsymbol{F}^{\top} = \frac{1}{J} \left[ \frac{\partial \bar{W}}{\partial \bar I_1} \left[2 J^{-2/3} \boldsymbol{B}- \frac{2}{3} \bar I_1 \boldsymbol{I}\right] + \frac{\partial \bar{W}}{\partial \bar I_2^{3/2}} \frac{\partial \bar I_2^{3/2}}{\partial \bar I_2} \left[2 J^{-4/3} [I_1 \boldsymbol{B}- \boldsymbol{B}\boldsymbol{B}] - \frac{4}{3} \bar I_2 \boldsymbol{I}\right] \right] - p \boldsymbol{I},\] with \(\boldsymbol{B}=\boldsymbol{F}\boldsymbol{F}^{\top}\), and hence is symmetric by construction.

4. Material symmetry: Since the strain-energy density is expressed solely as a function of invariants of the right Cauchy-Green tensor, the resulting material response is isotropic by construction. It is \(W(\boldsymbol{F}) = W(\boldsymbol{F}\boldsymbol{Q})\) for all rotation tensors in the special orthogonal group \(\boldsymbol{Q}\in \text{SO}(3)\).

5. Vanishing strain-energy density in the undeformed configuration: No deformation \(\boldsymbol{F}=\boldsymbol{I}\) implies \(\bar I_1 = 3\), \(\bar I_2 = 3^{3/2}\), and \(J=1\). Therefore, the constraint \(\bar{W}(0, 0) = 0\) ensures that the strain-energy density \(W\) vanishes in the undeformed configuration. This constraint is not a mathematical necessity. However, it represents a natural and widely accepted normalization of the strain-energy density.

6. Stress-free undeformed configuration: The isochoric contribution to the stress is \[\bar \boldsymbol{P} = \frac{\partial \bar{W}}{\partial \boldsymbol{F}} = \frac{\partial \bar{W}}{\partial \bar I_1} \frac{\partial \bar I_1}{\partial \boldsymbol{F}} + \frac{\partial \bar{W}}{\partial \bar I_2^{3/2}} \frac{\partial \bar I_2^{3/2}}{\partial \bar I_2} \frac{\partial \bar I_2}{\partial \boldsymbol{F}},\] with \[\frac{\partial \bar I_1}{\partial \boldsymbol{F}} = 2 J^{-2/3} \boldsymbol{F}- \frac{2}{3} \boldsymbol{F}^{-\top}\bar I_1 \quad and \quad \frac{\partial \bar I_2}{\partial \boldsymbol{F}} = 2 J^{-4/3} [I_1 \boldsymbol{F}- \boldsymbol{F}\boldsymbol{C}] - \frac{4}{3} \boldsymbol{F}^{-\top}\bar I_2.\] In the undeformed configuration, \(\boldsymbol{F}=\boldsymbol{I}\) implies \({\partial \bar I_1}/{\partial \boldsymbol{F}} = \boldsymbol{0}\) and \({\partial \bar I_2}/{\partial \boldsymbol{F}} = \boldsymbol{0}\). This, in turn, implies vanishing stress \(\bar \boldsymbol{P}=\boldsymbol{0}\).

7. Polyconvexity: A simple way to construct polyconvex strain-energy functions is by taking positive sums of monotone and convex functions of polyconvex invariants. The first isochoric principal invariant \(\bar I_1\) is polyconvex, however, the second isochoric principal invariant \(\bar I_2\) is not polyconvex [86]. Hence, we follow [86] and more recently [87] and express the isochoric strain-energy density as monotone and convex functions of the modified set of polyconvex invariants \(\bar\boldsymbol{I}^*\). While polyconvexity is not a strict physical requirement, when combined with appropriate coercivity conditions, it provides a sufficient framework to guarantee the existence of minimizers for the associated forward problem [88].

8. Differentiability: By ensuring that the strain-energy density is at least twice differentiable, \(\bar W \in C^2\left(\mathcal{I}_{\text{adm}}; \mathbb{R}_{\geq 0}\right)\), we guarantee well-defined stress and stiffness tensor computations.

For a more detailed description of how to enforce physical admissibility in hyperelasticity, we refer the reader to the manual provided by [58].

2.1.2 Data space↩︎

We consider a displacement-controlled experiment on an incompressible hyperelastic material over the time interval \([0,T]\). The undeformed three-dimensional specimen occupies the domain \(\Omega^{\text{3D}}\). We denote the Dirichlet boundaries at which the displacements are prescribed in \(X_i\)-direction by \(\Gamma^{\text{3D}}_{i}\). While conducting the experiment, observable quantities such as the displacement fields on a two-dimensional surface \(\Omega^{\text{2D}}\) of the specimen and the reaction force function at the boundary can be measured. The observable displacement fields are elements of \[\mathcal{U} = \left\{ \boldsymbol{u}(\boldsymbol{X},t) : \Omega^{\text{2D}} \times [0,T] \rightarrow \mathbb{R}^2 \;\middle| \; \begin{align} &u_1(\boldsymbol{X},t) = \bar u_1(\boldsymbol{X},t) && \forall \boldsymbol{X}\in \Gamma^{\text{2D}}_{1},\;\forall t \in [0,T], \\ &u_2(\boldsymbol{X},t) = \bar u_2(\boldsymbol{X},t) && \forall \boldsymbol{X}\in \Gamma^{\text{2D}}_{2},\;\forall t \in [0,T], \\ &\boldsymbol{u}(\boldsymbol{X},0) = \boldsymbol{0} && \forall \boldsymbol{X}\in \Omega^{\text{2D}} \end{align} \right\},\] where \(\bar u_1(\boldsymbol{X},t)\), \(\bar u_2(\boldsymbol{X},t)\) are the displacements prescribed on the Dirichlet boundaries \(\Gamma^{\text{2D}}_{i} = \Omega^{\text{2D}} \cap \Gamma^{\text{3D}}_{i}\). The observable reaction force functions, which are assumed to be scalar-valued in this work, are elements of \[\mathcal{R} = \left\{ R(t) : [0,T] \rightarrow \mathbb{R}\;\middle| \;R(0) = 0 \right\}.\] We denote a tuple consisting of an observable displacement field and a reaction force function by \(\boldsymbol{d}=(\boldsymbol{u},R)\) and note that not every pair \(\boldsymbol{d}\in \mathcal{D} = \mathcal{U} \times \mathcal{R}\) is physically meaningful. The observable displacement field and reaction force function must correspond to the same solution of the forward boundary value problem for an admissible material model. Hence, in the following, we introduce the constrained space \(\mathcal{D}_{\text{BVP}}\subset\mathcal{D}\) that only contains physically meaningful measurements.

For a given material \(\bar W \in \mathcal{W}\), the displacement \(\boldsymbol{u}^{\text{3D}}_{\bar W}\) and pressure field \(p^{\text{3D}}_{\bar W}\) over the three-dimensional domain \(\Omega^{\text{3D}}\) are determined by the forward boundary value problem of incompressible hyperelasticity (see 6). The corresponding deformation gradient and stress fields are \[\boldsymbol{F}^{\text{3D}}_{\bar W} = \nabla_{\boldsymbol{X}} \boldsymbol{u}^{\text{3D}}_{\bar W} + \boldsymbol{I}, \quad \boldsymbol{P}^{\text{3D}}_{\bar W} = \left. \frac{\partial \bar W}{\partial \boldsymbol{F}}\right\rvert_{\boldsymbol{F}=\boldsymbol{F}^{\text{3D}}_{\bar W}} - p^{\text{3D}}_{\bar W} [\boldsymbol{F}^{\text{3D}}_{\bar W}]^{-\top}.\] During the experiment, the three-dimensional fields are not directly observable. Instead, only selected quantities can be observed. We denote the observable two-dimensional displacement field over the specimen surface by \(\boldsymbol{u}^{\text{2D}}_{\bar W}\), and the observable reaction force function by \[R_{\bar W} =\int_{\Gamma^{\text{3D}}_{R}} \boldsymbol{P}^{\text{3D}}_{\bar W} \boldsymbol{N}\cdot \boldsymbol{e}_R \mathrm{d}A,\] where \(\Gamma^{\text{3D}}_{R}\) denotes the boundary at which the reaction force is measured, \(\boldsymbol{N}\) is the corresponding outward unit normal vector, and \(\boldsymbol{e}_R\) is the direction of the measured reaction force. The observed displacement fields and reaction force functions are fully determined by the material \(\bar W\). We hence define the observation mapping \[O : \mathcal{W} \rightarrow \mathcal{D}_{\text{BVP}}, \quad O : \bar W \mapsto (\boldsymbol{u}^{\text{2D}}_{\bar W},R_{\bar W}).\] Consequently, the constrained space of physically meaningful observations is \(\mathcal{D}_{\text{BVP}}=O(\mathcal{W})\subset\mathcal{D}\). If we assume a material \(\bar W_{\text{sim}} \in \mathcal{W}\) and synthetically generate data \(\boldsymbol{d}_{\text{sim}}\in\mathcal{D}_{\text{BVP}}\), it is \(O(\bar W_{\text{sim}}) = \boldsymbol{d}_{\text{sim}}\). In contrast, due to measurement noise and erroneous modeling assumptions, an experimental measurement \(\boldsymbol{d}_{\text{obs}}\in\mathcal{D}\) may not lie in the space \(\mathcal{D}_{\text{BVP}}\). Instead, the experimental measurement suggests that the material model lies in \[\mathcal{W}_{\boldsymbol{d}_\text{obs},\varepsilon} = \left\{ \bar W\in\mathcal{W} : \left\| O(\bar W)-\boldsymbol{d}_{\text{obs}} \right\|_{\mathcal{D}} \le \varepsilon \right\},\] where \(\varepsilon=0\) for synthetic data and \(\varepsilon>0\) for experimental data.

2.2 Statistical conditional inverse operator↩︎

In general, \(\mathcal{W}_{\boldsymbol{d}_\text{obs},\varepsilon}\) does not contain a unique material model. Thus, the inverse problem we aim to tackle with our proposed neural operators should not be understood as the application of a classical single-valued inverse operator \(S=O^{-1}\). Therefore, in this section, rather than postulating a canonical selector from \(\mathcal{W}_{\boldsymbol{d}_\text{obs},\varepsilon}\), we interpret material model discovery as a supervised statistical learning problem. We seek to learn a conditional inverse operator \(S_{\mu}: \mathcal{D}_{\text{BVP}} \rightarrow \mathcal{W}\) that maps an observation to the Bayes-optimal predictor induced by the training distribution \(\mu\). Although the operator is learned using observations in \(\mathcal{D}_{\text{BVP}}\), it is expected to generalize to observations lying near \(\mathcal{D}_{\text{BVP}}\).

Let \(\mu\) be the probability distribution from which training material models are sampled. We denote a random draw from this distribution by \(\bar W_{\mu} \sim \mu\) and the corresponding observation by \(\boldsymbol{d}_{\mu} = O(\bar W_{\mu})\). Further, let \[\mathcal{H} = L^2(\mathcal{I}_{\mathrm{adm}},\nu_{\mathcal{I}_{\mathrm{adm}}}) \quad \text{with} \quad \|V\|_{\mathcal{H}}^2 = \int_{\mathcal{I}_{\mathrm{adm}}} |V(\bar\boldsymbol{I}^*)|^2 \, \mathrm{d}\nu_{\mathcal{I}_{\mathrm{adm}}}(\bar\boldsymbol{I}^*),\] be the space of square-integrable functions over the admissible invariants domain, where \(\nu_{\mathcal{I}_{\mathrm{adm}}}\) denotes the probability distribution used to sample invariant query points for evaluating constitutive models. The Bayes-optimal predictor \(S_{\mu}: \mathcal{D}_{\text{BVP}} \rightarrow \mathcal{H}\) associated with the squared-error loss is given by \[S_{\mu} \in \operatornamewithlimits{arg\ min}_{S: \mathcal{D}_{\text{BVP}} \rightarrow \mathcal{W} } \, \mathbb{E} \left( \left\| {S}(\boldsymbol{d}_{\mu})-\bar W_{\mu} \right\|_{\mathcal{H}}^{2} \right).\] Equivalently, for a fixed observation \(\boldsymbol{d}_{\text{obs}}\in\mathcal{D}_{\mathrm{BVP}}\), the Bayes-optimal prediction \(\bar W_{\mu}( \, \cdot \, ; \boldsymbol{d}_{\text{obs}}) = S_{\mu}(\boldsymbol{d}_{\text{obs}})\) with \(\bar W_{\mu}( \, \cdot \, ; \boldsymbol{d}_{\text{obs}}): \mathcal{I}_{\text{adm}} \to \mathbb{R}_{\geq 0}\) satisfies \[\bar W_{\mu}( \, \cdot \, ; \boldsymbol{d}_{\text{obs}}) \in \operatornamewithlimits{arg\ min}_{V\in\mathcal{W}} \, \mathbb{E} \left( \left. \left\| V-\bar W_\mu \right\|_{\mathcal{H}}^{2} \,\right|\, \boldsymbol{d}_{\mu} = \boldsymbol{d}_{\text{obs}} \right).\]

If no two distinct material models represented by the training distribution produce the same observation, that is \[\bar W_{1},\bar W_{2} \in \operatorname{supp}\mu, \quad O(\bar W_{1}) = O(\bar W_{2}) \Rightarrow \bar W_{1}=\bar W_{2},\] then the material model is identifiable within the training class. In this case, \(S_{\mu}(O(\bar W_{\mu}))=\bar W_{\mu}\) for \(\mu\)-almost every \(\bar W\). If the inverse is not identifiable, several material models may correspond to the same observation. In that case, \(S_{\mu}(\boldsymbol{d}_{\text{obs}})\) is the prediction jointly preferred by the training distribution \(\mu\), the admissible hypothesis class, and the squared-error loss. In an unrestricted linear function space, the squared-error optimal prediction is the conditional mean \[\bar W_{\mu}( \, \cdot \, ; \boldsymbol{d}_{\text{obs}}) = \mathbb{E} \left[ \left. \bar W_\mu \right| \boldsymbol{d}_{\mu} = \boldsymbol{d}_{\text{obs}} \right].\] Such a prediction should not, in general, be interpreted as a unique physical inverse or as a canonical member of \(\mathcal{W}_{\boldsymbol{d}_{\text{obs}},\varepsilon}\).

The neural operators considered below provide a parameterized approximation of the statistical predictor \(S_{\mu}\) \[S_{\boldsymbol{\theta}}: \mathcal{D}_{\mathrm{BVP}} \to \mathcal{H}, \quad S_{\boldsymbol{\theta}}(\boldsymbol{d}) = \bar W_{\boldsymbol{\theta}} (\, \cdot \, ; \boldsymbol{d}),\] where the parameters \(\boldsymbol{\theta}\) are learned from labeled training data \((\boldsymbol{d}_i,\bar W_i)\) as discussed subsequently.

2.3 Neural operator architectures↩︎

In the following, we build the neural operators for learning a mapping from elements in \(\mathcal{D}\) to elements in \(\mathcal{W}\). The architectures chosen in this work are inspired by the DeepONet architectures proposed by [5], [7]. They consist of a branch net and a trunk net, as illustrated in 1. The trunk net maps the deformation gradient \(\boldsymbol{F}\) to a latent space. The discretization-independent and noise-robust branch net maps the given displacement and reaction force data to another latent space. The inner product of the latent spaces yields the strain-energy density \(W\). After a trained neural operator is queried with a given set of displacement and reaction force data, the branch net can be discarded, and the \(W\) can be evaluated for any arbitrary \(\boldsymbol{F}\). In this work, we consider two different architectures, the PANO (2) and the CANO (3), which share the same branch net, but differ in the design of the trunk net. In the following, we describe the shared branch net architecture and the different trunk net architectures.

The main purpose of the branch net is to encode the available displacement and reaction force data in a latent space (2 3). We assume the displacement fields over the specimen domain are available through full-field measurement techniques such as Digital Image Correlation. More specifically, we assume that the two-dimensional displacement vector \(\boldsymbol{u}\) is known at a total of \(N_{\boldsymbol{X}}\) points over the domain \(\Omega^{\text{2D}}\) and at a total of \(N_t\) time steps in the interval \([0,T]\), resulting in a displacement data array \(\boldsymbol{u}_{\text{data}} \in \mathbb{R}^{2 \times N_{\boldsymbol{X}} \times N_t}\). However, as we discuss at several points throughout this work, our neural operators are designed such that they can also take displacement fields sampled at arbitrary spatial locations in the domain \(\Omega^{\text{2D}}\) as input.

In the first layer of the branch net, the displacement field data are encoded in a lower-dimensional space. Specifically, we represent both components of the displacement field in the reduced basis \[\label{eq:Laplacian95interpolation} u_1(\boldsymbol{X},t) \approx u_1^{\text{lift}}(\boldsymbol{X},t) + \sum_{j=1}^{N_{\phi}} a_{1jt} \, \phi_{1j}(\boldsymbol{X}) \quad and \quad u_2(\boldsymbol{X},t) \approx u_2^{\text{lift}}(\boldsymbol{X},t) + \sum_{j=1}^{N_{\phi}} a_{2jt} \, \phi_{2j}(\boldsymbol{X}), \\\tag{2}\] where \(\phi_{1j}(\boldsymbol{X}) : \Omega^{\text{2D}} \rightarrow \mathbb{R}\) and \(\phi_{2j}(\boldsymbol{X}) : \Omega^{\text{2D}} \rightarrow \mathbb{R}\) are time-independent basis functions with non-local support on \(\Omega^{\text{2D}}\), and \(a_{1jt} \in \mathbb{R}\) and \(a_{2jt} \in \mathbb{R}\) are real-valued coefficients. The number of considered basis functions \(N_\phi\) for representing the displacement fields is chosen to be smaller than the number of available displacement measurements \(N_{\boldsymbol{X}}\). Hence, the first layer compresses the given displacement data array from \(\boldsymbol{u}_{\text{data}} \in \mathbb{R}^{2 \times N_{\boldsymbol{X}} \times N_t}\) to the encoded array \(\boldsymbol{a}_{\text{data}} \in \mathbb{R}^{2 \times N_\phi \times N_t}\). This lower-dimensional representation of the displacement field at each time step serves several purposes. First, the compressed displacement fields are easier for the branch net to process and reduce the number of its trainable weights. Second, because the basis functions have non-local support, the interpolation in 2 reduces the sensitivity of the branch net to local noise in the displacement field. Third, the branch net becomes independent of the spatial resolution because the interpolation in 2 can also be performed when displacements are provided at spatial locations other than the \(N_{\boldsymbol{X}}\) points discussed above. We finally note that the interpolation in 2 is fully deterministic and does not depend on any trainable weights.

In this work, we choose the basis functions \(\phi_{1j}(\boldsymbol{X})\) and \(\phi_{2j}(\boldsymbol{X})\) as the first \(N_\phi\) Laplacian eigenfunctions of the geometry \(\Omega^{\text{2D}}\). These functions are obtained by solving a generalized eigenvalue problem that is derived by discretizing the Laplace partial differential equation with a reference mesh on the geometry. A detailed description of the Laplacian eigenfunction computation is provided in 9. The Laplacian eigenfunctions vanish at the Dirichlet boundaries \(\Gamma^{\text{2D}}_{1}\), \(\Gamma^{\text{2D}}_{2}\), respectively. The fixed lifting fields \(u_1^{\text{lift}}(\boldsymbol{X},t)\) and \(u_2^{\text{lift}}(\boldsymbol{X},t)\), which fulfill the Dirichlet boundary conditions, ensure that the interpolated displacement fields fulfill the Dirichlet boundary conditions. For an illustration of the first Laplacian eigenfunctions \(\phi_{1j}(\boldsymbol{X})\), we refer to 2 3.

In addition to the displacement field, the scalar-valued vertical reaction force is assumed to be measured by a load cell during the experiment at the same \(N_t\) time steps, resulting in a force data array \(\boldsymbol{R}_{\text{data}} \in \mathbb{R}^{N_t}\). The force data array is normalized as \(\boldsymbol{R}_{\text{data}} / \|\boldsymbol{R}_{\text{data}}\|\). Both the displacement data and the normalized force data are mapped to a latent space of non-negative real valued numbers \(\mathbb{R}_{\geq 0}^{N_{\text{latent}}}\) using a neural network. We denote this neural network by \[\label{eq:latent95network} \text{FFNN}: \mathbb{R}^{2 \times N_{\phi} \times N_t} \times \mathbb{R}^{N_t} \rightarrow \mathbb{R}_{\geq 0}^{N_{b}}.\tag{3}\] Here, we chose a conventional, fully-connected feed-forward neural network that takes the flattened arrays \(\boldsymbol{a}_{\text{data}}\) and \(\boldsymbol{R}_{\text{data}} / \|\boldsymbol{R}_{\text{data}}\|\) as input. Enforcing the non-negativity of the output \(\boldsymbol{b}\in \mathbb{R}_{\geq 0}^{N_{b}}\) is straightforward by applying an element-wise, non-negative activation function \(\sigma:\mathbb{R}\rightarrow \mathbb{R}_{\geq 0}\) in the last layer, such as the ReLU function or smooth alternatives. In general, the dimension of the latent space \(N_b\) is much smaller than the dimension of the data \(N_b \ll 2 \, N_{\phi} \, N_t + N_t\). The dimension of the latent space can be considered a hyperparameter. Its meaning becomes clearer in the following, where we introduce the trunk net for describing the material model.

While the branch net encodes the information from available displacement and reaction force data, the trunk net (2 3) describes the mechanical material behavior. Specifically, the material behavior is governed by the isochoric strain-energy density function in the constrained space \(\mathcal{W}\). In this work, we approximate elements in this space using either Physics-Augmented Neural Networks (PANNs) [52], [54], [58] or Constitutive Artificial Neural Networks (CANNs) [51], [55] \[\text{PANN}: \mathbb{R}_{\geq 0}^{2} \rightarrow \mathbb{R}_{\geq 0}^{N_{b}} \quad or \quad \text{CANN}: \mathbb{R}_{\geq 0}^{2} \rightarrow \mathbb{R}_{\geq 0}^{N_{b}},\] which are designed to fulfill the desired physical constraints discussed above. A detailed description of these networks is provided in the subsequent [sec:trunk_pano] [sec:trunk_cano]. Note that, as opposed to the original PANN and CANN formulations, our networks do not map directly to the isochoric strain-energy density function but instead to a latent feature vector \(\boldsymbol{\psi}\in \mathbb{R}_{\geq 0}^{N_{b}}\), which shares the same dimension as the latent output space of the branch net \(\boldsymbol{b}\in \mathbb{R}_{\geq 0}^{N_{b}}\). Each entry in this vector can be interpreted as a physics-constraint modeling feature for the isochoric strain-energy density function. The isochoric strain-energy density function is then finally obtained through the inner product of the outputs of the branch and trunk nets \[\label{eq:inner95bt} \bar W = \boldsymbol{b}\cdot \boldsymbol{\psi}\cdot \|\boldsymbol{R}_{\text{data}}\|,\tag{4}\] where we use the norm of the reaction force data to appropriately scale the isochoric strain-energy density.

The composition of the branch and trunk nets forms our physics-aware neural operators for constitutive model discovery. The neural operators can be trained in a supervised manner. To this end, labeled data of the form \((\boldsymbol{u}_{\text{data}},\boldsymbol{R}_{\text{data}},\bar\boldsymbol{I}^*_{\text{data}},\bar \boldsymbol{W}_{\text{data}})_i\) with \(i=1,\dots,N_{\text{data}}\) are generated in a simulation-based data-mining process. To generate a single data point, an isochoric strain-energy density function \(\bar W\) is assumed and evaluated at a total of \(N_{\bar I}\) instances of invariant pairs \((\bar I_1, \bar I_2)\). This yields the data arrays \(\bar\boldsymbol{I}^*_{\text{data}} \in \mathbb{R}^{2 \times N_{\bar I}}\) and \(\bar \boldsymbol{W}_{\text{data}} \in \mathbb{R}^{N_{\bar I}}\). Next, the displacements and reaction forces are simulated for the standardized experiments, assuming the same isochoric strain-energy density function \(\bar W\). This provides the data arrays \(\boldsymbol{u}_{\text{data}}\) and \(\boldsymbol{R}_{\text{data}}\), where the displacement data \(\boldsymbol{u}_{\text{data}}\) is compressed into \(\boldsymbol{a}_{\text{data}}\) through the Laplacian encoding (9).

More information regarding training will be provided in 2.5. After training is completed, the neural operator is capable of predicting the isochoric strain-energy density function for a given set of experimentally measured displacements and reaction forces. This function can subsequently be evaluated for arbitrary input pairs of invariants. In this way, the neural operator learns a mapping between function spaces rather than between finite-dimensional Euclidean spaces. When experimental displacements and reaction forces are provided as input, the predicted isochoric strain-energy density function can be directly employed in forward simulations. Notably, for this purpose, it is not necessary to implement the entire neural operator. The branch net can be omitted, as only its outputs are required for the implementation of the material model.

2.3.1 PANO↩︎

In the following, we describe the trunk net architectures in detail. 2 shows an illustration of the PANN-inspired trunk architecture of the PANO. We assume that the dimension \(N_b\) of the latent space is even, and define the trunk architecture as \[\label{eq:features95pano} \begin{align} \psi_i(\boldsymbol{F}) = \begin{cases} \text{CMNN}_i(\bar I_1^*) - \text{CMNN}_i(0) & \text{if } i \leq \frac{N_b}{2} \\ \text{CMNN}_i(\bar I_2^*) - \text{CMNN}_i(0) & \text{if } i > \frac{N_b}{2} \end{cases}, \end{align}\tag{5}\] where \(\text{CMNN}_i : \mathbb{R}_{\geq 0} \rightarrow \mathbb{R}_{\geq 0}\) for \(i=1,\dots,N_b\) are neural networks that are convex and monotone with respect to their inputs. Such networks can be constructed through a slight modification of the input-convex neural networks proposed by [89]. For detailed information on how to impose the desired convexity and monotonicity constraints, we refer to [52]–[54], [90]. Due to their specifically designed architecture, the features \(\psi_i(\boldsymbol{F})\) are convex and monotone with respect to the normalized polyconvex invariants. Hence, the isochoric strain-energy density in 4 is polyconvex. The constants \(\text{CMNN}_i(0)\) ensure that the features vanish in the undeformed configuration \(\psi_i(\boldsymbol{I})=0\). We note that the convex, monotone neural networks \(\text{CMNN}_i\) have independent, trainable weights. The flexibility of the trunk net in fitting the constitutive material response is determined by the number of layers and neurons in \(\text{CMNN}_i\), as well as by the number of features \(N_b\).

2.3.2 CANO↩︎

The PANN-inspired trunk architecture leaves the functional form of the latent features \(\boldsymbol{\psi}\in \mathbb{R}_{\geq 0}^{N_b}\) flexible and determines it during training. As an alternative, we investigate a CANN-inspired trunk architecture in which \(\boldsymbol{\psi}\) is composed of a priori chosen, well-established constitutive features. Specifically, we consider the first six features of a generalized Mooney–Rivlin material model, which can be derived from a Taylor expansion of the isochoric strain-energy density function about the undeformed configuration (8) \[\label{eq:features95cano} \boldsymbol{\psi}(\boldsymbol{F}) = [ \bar I_1^*, ~ \bar I_2^*, ~ [\bar I_1^*]^2, ~ [\bar I_2^*]^2, ~ [\bar I_1^*]^3, ~ [\bar I_2^*]^3 ]^{\top}.\tag{6}\] Since the polynomial features are convex and monotone for non-negative inputs, the resulting isochoric strain-energy density in 4 is polyconvex. We note that \(\psi_i(\boldsymbol{I})=0\). The CANO architecture is illustrated in 3. The trunk net of the CANO is less flexible in fitting the constitutive material response than that of the PANO. However, this reduced flexibility decreases the number of trainable weights and makes the material model description more interpretable.

2.4 Experimental design↩︎

In this work, we assume an experimental setup with a fixed geometry, prescribed Dirichlet boundary conditions, and a fixed number of time steps. The experiment is designed to be reproducible for different materials across laboratories to enable a neural operator that is trained on this setup to be used repeatedly for constitutive model discovery. While we consider a relatively simple specimen geometry here, future work may explore specimen geometries optimized for material characterization [91], [92].

As we will show later, a neural operator trained on the reference geometry can be used for inference on geometries with different side lengths and thicknesses. Furthermore, we will demonstrate that, although the discretization is assumed to be fixed during training data generation, the neural operators can be applied to different discretizations of the geometry and can accommodate situations in which only partial displacement field information is available over the specimen surface. Spatial discretization independence is important because we cannot always guarantee that the full displacement fields are captured at all time steps throughout an experiment. In contrast, temporal discretization independence is generally less critical, since displacements and reaction forces can typically be measured at high sampling frequencies [68], [78] such that we can construct data at arbitrary time points through interpolation.

Here, we consider a thin plate with a central hole, as illustrated in 4. In the undeformed configuration, the plate occupies the domain \(\Omega^{\text{3D}}\). The length and height of the specimen are \(2\), and the thickness is \(0.01\). The plate is clamped at the top and bottom boundaries and is deformed under displacement control over \(N_t\) time steps in constant increments until a maximum displacement of \(\bar u_2 = -2.0\) is reached. Due to the symmetry of the problem, only one eighth of the plate is considered in the simulations, resulting in a computational domain with length and height equal to 1 and a thickness of 0.005. On the symmetry-reduced domain, the displacement in the \(x_1\)-direction is prescribed at the bottom and right boundaries, denoted by \(\Gamma^{\text{3D}}_{1}\). The displacement in the \(x_2\)-direction is prescribed at the top and bottom boundaries, denoted by \(\Gamma^{\text{3D}}_{2}\). The displacement in the \(x_3\)-direction is fixed at the back surface. We assume that the two-dimensional displacement field on the surface of the specimen, denoted by \(\Omega^{\text{2D}}\), can be measured using image-based measurement techniques such as Digital Image Correlation (DIC), and the scalar-valued vertical reaction force can be measured using a load cell. In actual experiments, the symmetry condition may be slightly violated. In such cases, the displacement field over \(\Omega^{\text{2D}}\) can be obtained by averaging the displacement fields of the four quadrants of the surface.

We note that, for the considered boundary value problem under displacement control and near plane stress conditions, certain rescalings of the geometry and boundary conditions leave the solution invariant up to a scaling factor. Specifically, scaling the plate thickness by a factor \(s_{3}\) scales the reaction force by the same factor while leaving the displacements on the measurement surface unchanged, provided that the thickness remains sufficiently small compared to the length and height, so that the plane stress assumption remains valid. Furthermore, if the length, height, and hole radius are scaled by a factor \(s_{1}\), and the prescribed displacement is scaled by the same factor while keeping the thickness fixed, both the displacement field and the reaction force scale by \(s_{1}\). These scaling relations will be exploited later to enable a neural operator trained on a reference geometry to perform inference on data generated from geometries with different side lengths and thicknesses.

2.5 Training↩︎

To generate training data, we discretize the three-dimensional geometry and perform finite element simulations to solve the displacement-controlled boundary value problem for different strain-energy density functions. The finite element simulations are performed using FEniCSx. From the simulation results, we extract the two-dimensional displacement fields on the surface of the specimen \(\Omega^{\text{2D}}\) at \(N_{\boldsymbol{X}}\) points. For details, we refer to 6.

To obtain a meaningful dataset that covers a large portion of the constitutive modeling space, we consider the Taylor expansion of the isochoric strain-energy density function about the undeformed configuration (see 8 for details). Since classical material models are generally restricted to low-order polynomial terms, we restrict ourselves to terms up to third order \[\bar W(\bar I_1^*,\bar I_2^*) = C_{10} \bar I_1^* + C_{01} \bar I_2^* + C_{20} [\bar I_1^*]^2 + C_{02} [\bar I_2^*]^2 + C_{30} [\bar I_1^*]^3 + C_{03} [\bar I_2^*]^3.\] We use Latin hypercube sampling to sample parameters in the six-dimensional space \([0,1]^6\) to create parameter samples \(\bar C_{ij} \in [0,1]\). For inputs greater than one, the higher-order polynomial features exhibit a greater influence on the total model. The maximum invariant \(\bar I_1^*\) reached in the simulations is about 3, and the maximum invariant \(\bar I_2^*\) is about 5. Therefore, we scale the parameters obtained through Latin hypercube sampling by \(C_{10}=\bar C_{10}/3\), \(C_{01}=\bar C_{10}/5\), \(C_{20}=\bar C_{20}/3^2\), \(C_{02}=\bar C_{02}/5^2\), \(C_{30}=\bar C_{30}/3^3\), \(C_{03}=\bar C_{03}/5^3\). We note that, owing to the reaction force normalization (see 2 3), the proposed neural operators are invariant under scalar scaling of the strain-energy density. Consequently, only the relative magnitudes of the constants \(C_{ij}\) are of interest, and it suffices to restrict the Latin hypercube sampling to the unit interval \([0,1]^6\).

The finite element simulations provide labeled data pairs of strain-energy density functions and the associated displacement and reaction force data. To train the neural operators, we sample the invariants uniformly from the domain \(\mathcal{I}_{\text{sample}} = \{ (\bar I_1 - 3, \bar I_2^{3/2} - 3^{3/2}) \in \mathcal{I}_{\text{adm}} \, | \, \bar I_2^{3/2} - 3^{3/2} \leq I_1 \}\). This sampling domain is chosen heuristically to approximately cover the range of invariants encountered in the training dataset. For an illustration of the invariants sampling, the reader is referred to 10 in 7. In our numerical experiments, alternative sampling domains did not significantly affect the results.

We obtain labeled data of the form \((\boldsymbol{u}_{\text{data}},\boldsymbol{R}_{\text{data}},\bar\boldsymbol{I}^*_{\text{data}},\bar \boldsymbol{W}_{\text{data}})_i\) with \(i=1,\dots,N_{\text{data}}\), \(\boldsymbol{u}_{\text{data}} \in \mathbb{R}^{2 \times N_{\boldsymbol{X}} \times N_t}\), \(\boldsymbol{R}_{\text{data}} \in \mathbb{R}^{N_t}\), \(\bar\boldsymbol{I}^*_{\text{data}} \in \mathbb{R}^{2 \times N_{\bar I}}\), and \(\bar \boldsymbol{W}_{\text{data}} \in \mathbb{R}^{N_{\bar I}}\). In this work, we run \(N_{\text{data}}=3000\) simulations. \(80\%\) of the data is used for training (\(N_{\text{data}}^{\text{train}}=0.8N_{\text{data}}\)), \(10\%\) for validation (\(N_{\text{data}}^{\text{val}}=0.1N_{\text{data}}\)), and \(10\%\) for testing (\(N_{\text{data}}^{\text{test}}=0.1N_{\text{data}}\)). The loss function is defined as \[\text{Loss}(\boldsymbol{\theta}) = \frac{1}{N_{\text{data}}^{\text{train}}} \sum_i^{N_{\text{data}}^{\text{train}}} \text{MSE}_i(\boldsymbol{\theta}) \quad \text{with} \quad \text{MSE}_i(\boldsymbol{\theta}) = \frac{1}{N_{\bar I}} \left\| \bar \boldsymbol{W}_{\text{NO}}\left((\boldsymbol{u}_{\text{data}},\boldsymbol{R}_{\text{data}},\bar\boldsymbol{I}^*_{\text{data}})_i;\boldsymbol{\theta}\right) - (\bar \boldsymbol{W}_{\text{data}})_i \right\|^2,\] where \(\bar \boldsymbol{W}_{\text{NO}}\) is the neural operator prediction, and \(\boldsymbol{\theta}\) represents all the trainable weights. The mean squared error \(\text{MSE}_i\) will be used in the following to assess the prediction accuracy of the neural operator for one individual training or testing sample \(i\). The hyperparameters of the training procedure are provided in 10.

3 Results↩︎

In the following, we present the results obtained with the trained PANO and CANO on the previously described training dataset, comprising 3,000 simulations generated by Latin hypercube sampling of the parameter domain of the cubic Taylor expansion model. We then discuss the discretization independence and noise robustness of the neural operators. Furthermore, we demonstrate that the neural operators can infer material models from full-field data even when the geometry side lengths and thickness differ from those assumed during the generation of the training data.

3.1 PANO↩︎

We use the Adam optimizer to train the PANO on \(80\%\) of the simulation data. \(10\%\) of the simulation data is used for validation during the training process. The corresponding validation loss is used to select the best model over all training epochs. Another \(10\%\) of the simulation data is held aside and is entirely unseen during training. These testing data are used in the following to assess the accuracy and reliability of the trained networks. For a detailed description of the hyperparameters of the neural operator architecture and the training process, we refer to 10.

5 (a) shows a histogram of the MSE computed for all simulations in the unseen testing dataset. The histogram indicates that the majority of samples exhibit low MSE values, while only a few outliers attain larger errors. We highlight the minimum, median, and maximum MSE values observed in the test dataset in the histogram. 5 (b) 5 (c) 5 (d) show the corresponding true and discovered isochoric strain-energy densities for the minimum, median, and maximum MSE cases, respectively. Excellent agreement is observed for the minimum and median MSE cases. Only the prediction corresponding to the maximum MSE exhibits noticeable deviations from the true isochoric strain-energy density. However, as shown by the histogram in 5 (a), this case represents an exceptional outlier, which is not representative of the overall predictive performance of PANO. The second worst MSE is substantially smaller than the worst MSE, and shows a significantly improved agreement between the true and discovered models.

3.2 CANO↩︎

Next, we use the Adam optimizer to train the CANO on \(80\%\) of the simulation data while using \(10\%\) of the data for validation and \(10\%\) of the data for testing. We refer to 10 for a detailed description of the hyperparameters. 6 (a) shows a histogram of the MSE computed from the unseen testing dataset, which shows that the majority of samples exhibit low MSE values. 6 (b) 6 (c) 6 (d) show the true and discovered models corresponding to the minimum, median, and maximum MSE. The agreement for the minimum and median MSE cases is excellent, while the prediction corresponding to the maximum MSE exhibits small deviations from the true isochoric strain-energy density.

A comparison of the PANO and CANO results reveals that the MSE distribution for CANO is shifted towards lower values, and the minimum and median MSE values are closer to one another than those obtained with PANO. This indicates that CANO predicts the underlying material models more consistently and accurately than PANO. Importantly, the maximum MSE achieved by CANO is lower than that of PANO. This improvement is reflected in the much closer agreement between the true and discovered isochoric strain-energy densities for the worst-performing case shown in 6 (d). These findings can be attributed to the constitutive structure that is explicitly embedded in the CANO architecture. By incorporating constitutive features directly into the network design, CANO is better equipped to learn accurate constitutive representations from the training data and to extrapolate to unseen data.

3.3 Noisy data↩︎

To assess the predictive performance of the neural operators under noisy measurements, we select a representative simulation from the test dataset and perturb the corresponding data with random noise. Specifically, we consider the simulation corresponding to the median MSE obtained by CANO. To emulate experimental measurements, uncorrelated Gaussian noise with a standard deviation \(\sigma\) is added to the displacement fields of all load steps. Following the discussion in [76], realistic noise levels are represented by \(\sigma=10^{-4}\) and \(\sigma=10^{-3}\). In addition, we investigate the extreme case of \(\sigma=10^{-2}\). Considering that the maximum absolute displacement in the first load step is \(0.2\), this corresponds to an exceptionally high noise level and therefore represents a challenging test case for the neural operators.

7 compares the true and discovered material models for all considered noise levels with the corresponding noise-free prediction. For the realistic noise levels of \(\sigma=10^{-4}\) and \(\sigma=10^{-3}\), the discovered material models exhibit only minor deviations from the ground truth, indicating that the neural operator is robust to measurement noise. Even for the extreme and rather unrealistic noise level of \(\sigma=10^{-2}\), the predicted material model remains reasonably close to the ground truth, demonstrating a remarkable degree of robustness under severely perturbed input data. The observed robustness to measurement noise can likely be attributed to the first layer of the branch net, where the displacement field is represented using Laplacian eigenfunctions. Owing to their non-local support, the Laplacian eigenfunctions inherently smooth the displacement field, thereby attenuating high-frequency noise while preserving the underlying deformation patterns.

3.4 Discretization independence and missing data↩︎

In the following, we investigate the robustness of the neural operators to incomplete data, varying spatial discretizations, and geometries with different dimensions. Such robustness is essential for practical applications, as it enables a trained neural operator to be applied to experimental data acquired on different measurement grids and from specimens with geometries that differ from those used during training.

First, we investigate the effect of incomplete measurement data on the predictions of the neural operators. This investigation is motivated by practical applications, where full-field measurements, such as those obtained from DIC, may contain missing or invalid data due to occlusions, poor image texture, or correlation failures. As a representative example, we consider the simulation corresponding to the median MSE obtained by CANO. From the full set of spatial measurement points, we randomly select a subset that is provided as input to the neural operator, while the remaining points are discarded. 8 (a) shows the spatial distribution of the measurement points for the complete dataset, whereas 8 (b) 8 (c) 8 (d) show the randomly selected subsets for progressively decreasing numbers of measurement points. The corresponding displacement fields are projected onto the Laplacian eigenfunction basis and subsequently serve as input to the branch net of the neural operators. The resulting true and discovered material models are shown in 8 (e) 8 (f) 8 (g) 8 (h). The results demonstrate that the predictions of the neural operators remain remarkably robust, exhibiting only minor changes even when a substantial fraction of the measurement data is omitted.

In practical applications, it may not always be feasible to fabricate specimens with the same geometric dimensions as those used to generate the training data. Nevertheless, the trained neural operators can also be applied to specimens of different sizes. The only requirements are that the geometry remains a thin square plate with a central hole, that the ratio of the hole radius to the specimen side length is preserved, and that the same amount of total stretch is applied. The absolute specimen dimensions, however, remain arbitrary, and both the side length and the thickness may differ from those considered during training. To demonstrate invariance with respect to specimen dimensions, we consider a specimen whose side length and thickness are scaled by factors of \(s_{1}=2.5\) and \(s_{3}=1.5\), respectively, relative to the geometry used for training. As a representative example, we select the material model corresponding to the median MSE obtained by CANO. The deformation of the scaled specimen is then simulated using a spatial discretization that deliberately differs from that employed during the generation of the training data. 9 (a) shows the spatial locations at which displacements are evaluated on the surface of the new specimen geometry. The resulting simulation data are then used as input to the neural operator. To this end, the displacement data are scaled by \(s_{1}\), and the reaction force is scaled by \(s_{1}\) and \(s_{3}\) prior to being passed to the branch net. 9 (b) shows the resulting discovered material model in comparison with the true material model. The good agreement indicates that the trained neural operators are not restricted to specific geometry dimensions or spatial discretizations.

4 Ill-posedness of the infinite-dimensional inverse problem↩︎

4.1 General remarks↩︎

It is well known that inverse problems, including those considered in this work, can be ill-posed. Different material modeling functions in \(\mathcal{W}\) may be equally consistent with a given set of observed data in \(\mathcal{U}\) and \(\mathcal{R}\) [83], [84]. In the following, we discuss the ill-posed and infinite-dimensional nature of the problem and present numerical experiments suggesting that our proposed neural operator architecture regularizes the inverse problem and reduces its ill-posedness.

A particularly challenging aspect of the learning problem considered in this work is that the functions in the image space \(\mathcal{W}\) are defined over the unbounded domain \(\mathcal{I}_{\text{adm}}\). This stands in contrast to operator learning for forward problems, where the quantities of interest are typically PDE solutions defined over a bounded domain. The deformations contained in any experimental dataset used to train our neural operators are necessarily bounded. Regardless of how the standardized experiment is designed, there will always exist deformations that lie outside the range encountered during training. This renders the unique identification of a function in \(\mathcal{W}\) for a given experimental dataset impossible. We note, however, that this characteristic of the inverse problem is not specific to our operator learning framework, but is inherent to material characterization frameworks in general, including traditional optimization-based and machine learning-based approaches.

Even if the domain of the functions in \(\mathcal{W}\) is assumed to be bounded, uniquely identifying a function in \(\mathcal{W}\) from a given experimental dataset remains challenging because of the infinite-dimensional nature of the problem. The experimental observations provide only finite-dimensional information about the underlying material model, whereas the sought function belongs to an infinite-dimensional space. Consequently, multiple functions in \(\mathcal{W}\) may be consistent with the same set of observations, even when their domains are restricted to a bounded set \(\mathcal{I}_{\text{bounded}} \subset \mathcal{I}_{\text{adm}}\).

In practice, additional assumptions or regularization are required to restrict the admissible solution space and obtain meaningful reconstructions. Hence, in our study, we restrict the design space of the constitutive functions. As described previously, the space \(\mathcal{W}\) is constrained to only include physically admissible functions. This reduces the space of admissible constitutive behaviors and thereby mitigates the ill-posedness of the problem. Furthermore, while neural operators are designed to address infinite-dimensional learning problems, their practical implementations are always restricted to finite-dimensional latent representations, such as the latent feature vector \(\boldsymbol{\psi}\) in our PANO and CANO architectures. The principal difference between PANO and CANO lies in the construction of these latent features. In CANO, the features are prescribed a priori, whereas in PANO they are learned during training. The universal approximation properties of neural operators have been studied theoretically [5]. However, these results are established in the asymptotic regime of infinitely many parameters and therefore provide only theoretical guarantees regarding the expressive power of the architectures. They do not directly address the approximation capabilities or regularization properties of practical implementations with finite-dimensional latent representations and finite network capacity.

The fact that the practical implementations of the PANO and CANO operate on finite-dimensional feature vectors \(\boldsymbol{\psi}\) provides an additional source of regularization. In the following, we present a numerical investigation to demonstrate that the proposed neural operator architectures regularize the learning problem sufficiently to render it well-posed.

4.2 Material model identifiability through the weak equilibrium equations↩︎

As an attempt to assess the well-posedness of the learning problem within a constrained class of material models, we investigate the sensitivity of the weak equilibrium equations and the measured reaction forces to perturbations in the material model. We restrict our analysis to the class of material models that can be represented by our CANO architecture and comment on the PANO architecture subsequently.

We consider CANO’s constitutive modeling features \(\psi_i(\boldsymbol{F})\) (see 6 ) obtained from the cubic Taylor expansion of the isochoric strain-energy density (see 8). Within this constrained class of material models, the constitutive behavior is fully characterized by the material parameters \(C_{ij}\), which we collect in the parameter vector \(\boldsymbol{c}\in \mathbb{R}^6\). Consequently, the strain-energy density and the first Piola-Kirchhoff stress can be expressed as \[W(\boldsymbol{F}) = c_i \psi_i(\boldsymbol{F}) - p [J-1], \quad \boldsymbol{P}(\boldsymbol{F}) = c_i \frac{\partial \psi_i}{\partial \boldsymbol{F}} - p \boldsymbol{F}^{-\top},\]

Our goal is to prove the identifiability of the unknown parameters \(\boldsymbol{c}\) through the weak equilibrium equations for the given displacement field and reaction force data. We consider the reduced plane stress problem and assume that the displacement fields \(\boldsymbol{u}_t(\boldsymbol{X}): \Omega^{\text{2D}}\rightarrow\mathbb{R}^2\) and reaction forces \(R_t\) are known from measurements at time steps \(t\). Under the incompressibility condition, the deformation gradient computes to \[\boldsymbol{F}_t(\boldsymbol{X}) =\begin{bmatrix} \boldsymbol{I}+ \nabla_{\boldsymbol{X}} \boldsymbol{u}_t(\boldsymbol{X}) & 0 \\ 0 & \frac{1}{\det(\boldsymbol{I}+ \nabla \boldsymbol{u}_t(\boldsymbol{X}))} \\ \end{bmatrix}.\]

Under the assumption of plane stress conditions, the first Piola-Kirchhoff stress components \(P_{13}\), \(P_{31}\), \(P_{23}\), \(P_{32}\), and \(P_{33}\) vanish identically. From \(P_{33}=0\), we compute the pressure term \(p={c_i \, {\partial \psi_i}/{\partial F_{33}}} \cdot {1}/{F^{-1}_{33}}\), and the stress at time step \(t\) is computed to \[\boldsymbol{P}(\boldsymbol{F}_t) = c_i \left[ \left.\frac{\partial \psi_i}{\partial \boldsymbol{F}}\right|_{\boldsymbol{F}=\boldsymbol{F}_t} - \left.{\frac{\partial \psi_i}{\partial F_{33}}}\right|_{\boldsymbol{F}=\boldsymbol{F}_t} \frac{\boldsymbol{F}_t^{-\top}}{F^{-1}_{t33}}\right].\]

Next, we show that the parameters \(\boldsymbol{c}\) can be uniquely identified from the given data through the weak equilibrium equations. This procedure is closely aligned with weak formulation-based inverse parameter identification methods from full-field measurements, such as the Virtual Fields Method (VFM) [71], [72], the Equilibrium Gap Method (EGM) [93], or EUCLID [76], [94]. We discretize the weak form of the linear momentum balance by introducing linear shape functions \(N_k(\mathbf{X})\) for the discretization of the test functions. Denoting basis vectors by \(\boldsymbol{e}_j\), this leads to the following expression for the internal nodal forces \[f^{\text{int}}_{tjk} = \int_\Omega^{\text{2D}} \, \boldsymbol{P}(\boldsymbol{F}_t) \colon \nabla (N_k \boldsymbol{e}_j) \, \mathrm{d}\boldsymbol{X},\] where \(t\) denotes the time step, \(j\) the spatial dimension, and \(k\) the node in the finite element mesh. By substituting the stress, we obtain \[f^{\text{int}}_{tjk} = c_i \int_\Omega^{\text{2D}} \, \left[ \left.\frac{\partial \psi_i}{\partial \boldsymbol{F}}\right|_{\boldsymbol{F}=\boldsymbol{F}_t} - \left.{\frac{\partial \psi_i}{\partial F_{33}}}\right|_{\boldsymbol{F}=\boldsymbol{F}_t} \frac{\boldsymbol{F}_t^{-\top}}{F^{-1}_{t33}}\right] \colon \nabla (N_k \boldsymbol{e}_j) \, \mathrm{d}\boldsymbol{X},\] where all integrals are evaluated using standard numerical quadrature. All internal forces are collected in the force vector \(\boldsymbol{f}^{\text{int}}\in\mathbb{R}^{2N_tN_{\boldsymbol{X}}}\) which is written as \[\boldsymbol{f}^{\text{int}} = \boldsymbol{A}^{\text{int}} \boldsymbol{c},\] where the matrix \(\boldsymbol{A}^{\text{int}} \in \mathbb{R}^{2N_tN_{\boldsymbol{X}} \times 6}\) collects contributions arising from the integrals above. For a detailed description of the assembly of this matrix, we refer to [76]. At all free degrees of freedom, the internal nodal forces should vanish. We introduce the selection matrix \(\boldsymbol{S}^{\text{free}}\) that selects the appropriate entries from the internal nodal force vector \(\boldsymbol{f}^{\text{free}}=\boldsymbol{S}^{\text{free}}\boldsymbol{f}^{\text{int}}=\boldsymbol{S}^{\text{free}}\boldsymbol{A}^{\text{int}} \boldsymbol{c}= \boldsymbol{0}\). Further, the internal nodal forces at the boundary must sum up to the measured reaction forces. We introduce the matrix \(\boldsymbol{S}^{\text{react}}\) such that \(\boldsymbol{S}^{\text{react}}\boldsymbol{f}^{\text{int}}=\boldsymbol{S}^{\text{react}}\boldsymbol{A}^{\text{int}} \boldsymbol{c}= \boldsymbol{R}^{\text{react}}\), where \(\boldsymbol{R}^{\text{react}}\) collects the reaction forces of all load steps. We concatenate the equations and obtain \[\label{eq:equilibrium95gap} \boldsymbol{A}\boldsymbol{c}= \boldsymbol{R}, \quad \text{with} \quad \boldsymbol{A}= \begin{bmatrix} \boldsymbol{S}^{\text{free}}\boldsymbol{A}^{\text{int}} \\ \boldsymbol{S}^{\text{react}}\boldsymbol{A}^{\text{int}} \\ \end{bmatrix}, ~ \boldsymbol{R}= \begin{bmatrix} \boldsymbol{0} \\ \boldsymbol{R}^{\text{react}} \\ \end{bmatrix}.\tag{7}\] The system above can be used to identify the parameters \(\mathbf{c}\) in a least-squares sense. The system of equations clarifies why the reaction forces are needed for the inverse problem. The equations corresponding to the free degrees of freedom are homogeneous. They may constrain the relative combination of material parameters, but they do not determine their absolute value. The measured reaction force provides this scale information.

In the following, we use the system of equations above to assess the well-posedness of the inverse problem by studying the properties of \(\boldsymbol{A}\). We consider two parameter vectors \(\boldsymbol{c}\) and \(\boldsymbol{c}+\Delta \boldsymbol{c}\). They fulfill the equilibrium equations if \(\boldsymbol{A}\boldsymbol{c}= \boldsymbol{R}\) and \(\boldsymbol{A}[\boldsymbol{c}+\Delta \boldsymbol{c}] = \boldsymbol{R}\). Subtracting the two corresponding systems gives \[\boldsymbol{A}\Delta \boldsymbol{c}= \boldsymbol{0}.\] Thus, the first-level non-uniqueness within the restricted material class is \[\ker\boldsymbol{A} = \left\{ \Delta \boldsymbol{c}\in\mathbb{R}^6: \boldsymbol{A}\Delta \boldsymbol{c}= \boldsymbol{0} \right\}.\] We compute the singular value decomposition \[\boldsymbol{A}= \boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^{\top}.\] The numerical rank is defined as \[\operatorname{rank}_{\tau}(\boldsymbol{A}) = \#\left\{ i: \frac{\sigma_i}{\sigma_{\max}}>\tau \right\},\] where \(\tau\) is a prescribed relative tolerance, \(\sigma_i=\Sigma_{ii}\) (no Einstein summation) are the singular values of \(\boldsymbol{A}\), and \(\sigma_{\max}\) is the largest singular value. The numerical nullity is then \[\dim\ker_{\tau}\boldsymbol{A} = 6-\operatorname{rank}_{\tau}(\boldsymbol{A}).\] If \(\boldsymbol{A}\) possesses right singular vectors associated with zero or sufficiently small singular values, these vectors represent material parameter perturbations that produce no observable change in the measured displacement and reaction force data. Such directions correspond to non-identifiable parameter combinations, indicating that the inverse problem is ill-posed. Conversely, if all singular values exceed the prescribed relative tolerance, every material parameter perturbation produces a detectable change in the available observations. In this case, no non-identifiable directions exist, and the inverse problem is well-posed.

We assemble the matrix \(\boldsymbol{A}\) and compute its singular values for representative datasets in the training data. For all cases considered, the smallest singular value is at least \(15\%\) of the largest singular value, and is therefore well above the numerical tolerance \[\frac{\sigma_{\min}}{\sigma_{\max}} \gtrsim 0.15 \gg \tau.\] Thus, \(\boldsymbol{A}\) has full numerical rank, which implies that its numerical nullity is zero \[\operatorname{rank}_{\tau}(\boldsymbol{A})=6, \qquad \dim\ker_{\tau} \boldsymbol{A}=0.\] Consequently, no non-identifiable material parameter directions are detected within the prescribed numerical tolerance, and our standardized experiment exhibits no detectable first-level non-uniqueness. We note that this conclusion is restricted to the chosen material class with six parameters and to the fixed-displacement first-level analysis. It does not prove the uniqueness of the full constitutive-function inverse problem. Non-uniqueness may reappear when the material class is enriched by introducing coupling terms, higher-order terms, or more general constitutive features, such as those in the trunk net of the PANO architecture. For the latter case, the numerical rank and nullity of the matrix \(\boldsymbol{A}\) depend on the feature functions \(\psi_i(\boldsymbol{F})\) learned during PANO training.

5 Conclusions↩︎

We presented a novel operator learning framework for directly inferring constitutive functions from experimentally measurable quantities, namely displacement fields and reaction force functions. The proposed network architectures, PANO and CANO, successfully learn the mapping from experimental data to constitutive functions using the training dataset and demonstrate accurate predictions on unseen test data. Importantly, the projection onto Laplacian eigenfunctions enables the trained neural operators to perform robust inference from noisy measurements, incomplete data, varying spatial discretizations, and specimens with different dimensions. At the same time, the discovered constitutive functions can be evaluated at arbitrary deformation states, which makes the predictions independent of the discretization used during training.

Several directions for future research emerge from the present work. The proposed framework was developed and validated for a thin plate with a central hole. Extending the approach to varying specimen geometries represents a natural next step and would further enhance its applicability to experimental settings. The present framework employs Laplacian eigenfunctions to encode the displacement fields. Future work may systematically investigate alternative basis functions and assess their influence on the accuracy, robustness, and generalization capabilities of the neural operators. Although our proposed approach is independent of the spatial discretization, it is not invariant with respect to the temporal discretization. This is a reasonable assumption, as experimental data can generally be acquired at high frequencies in time. Nevertheless, future work could achieve time-discretization independence by representing both the displacement fields and the reaction force functions via an appropriate temporal basis. In future work, we further aim to broaden the applicability of neural operators across different experimental configurations by making them independent of specimen geometry and topology. Furthermore, the present study focused on data obtained from experiments on complex specimen geometries. Extending the framework to simpler experimental configurations with homogeneous deformation states is expected to be straightforward and may facilitate applications in standard material characterization procedures. Finally, an important direction for future research concerns the design of the constitutive model space used during training data generation. We expect that the generality and predictive capabilities of the proposed framework can be further improved by enriching the training database with a broader range of constitutive models. In particular, future work will consider more complex material classes, including anisotropic and dissipative constitutive behavior.

In summary, our work demonstrates that the inverse discovery of constitutive functions from experimental data can be significantly accelerated by leveraging the information contained in previously generated simulation data. In conventional approaches, such as Finite Element Model Updating, solving an inverse problem requires a large number of forward simulations, which are typically discarded once the optimization has converged. We envision a paradigm in which these simulation datasets are instead retained and used to train general neural operators that can subsequently solve inverse problems directly from experimental measurements to eliminate the need for repeated computationally expensive forward simulations.

Code and data availability↩︎

The code and data will be publicly available upon publication at https://github.com/mflaschel.

Acknowledgments↩︎

The authors gratefully acknowledge the use of an existing finite element method implementation in FEniCSx developed by Miguel Angel Moreno-Mateos and Simon Wiesheier. The authors acknowledge support from the European Research Council (ERC) Grant 101141626 DISCOVER. Funded by the European Union. Views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

6 Forward boundary value problem↩︎

In the following, we detail the forward boundary value problem that is solved numerically to generate the simulation data used for training the neural operators. We consider an incompressible hyperelastic material that occupies the domain \(\Omega^{\text{3D}}\) in the undeformed configuration and deforms over the time interval \([0,T]\). We neglect the moment of inertia and body forces throughout this work. Further, we assume a specimen under full displacement control, such that the Neumann boundary conditions are homogeneous. The boundary value problem then states \[\begin{align} &\boldsymbol{F}(\boldsymbol{X},t) = \nabla_{\boldsymbol{X}} \boldsymbol{u}(\boldsymbol{X},t) + \boldsymbol{I}\quad \text{in } \Omega^{\text{3D}} \times [0,T],\\ &J(\boldsymbol{X},t) = \det(\boldsymbol{F}(\boldsymbol{X},t)) = 1, \quad \text{in } \Omega^{\text{3D}} \times [0,T],\\ &\boldsymbol{P}(\boldsymbol{X},t) = \frac{\partial \bar W}{\partial \boldsymbol{F}} - p(\boldsymbol{X},t) \boldsymbol{F}^{-\top}(\boldsymbol{X},t), \quad \text{in } \Omega^{\text{3D}} \times [0,T],\\ &\nabla_{\boldsymbol{X}} \cdot \boldsymbol{P}(\boldsymbol{X},t) = \boldsymbol{0}, \quad \text{in } \Omega^{\text{3D}} \times [0,T],\\ &u_1(\boldsymbol{X},t) = \bar u_1(\boldsymbol{X},t), \quad \text{on } \Gamma^{\text{3D}}_1 \times [0,T],\\ &u_2(\boldsymbol{X},t) = \bar u_2(\boldsymbol{X},t), \quad \text{on } \Gamma^{\text{3D}}_2 \times [0,T],\\ &P_{1j}(\boldsymbol{X},t) N_j(\boldsymbol{X},t) = 0, \quad \text{on } \Gamma^{\text{3D}}_{N1} \times [0,T],\\ &P_{2j}(\boldsymbol{X},t) N_j(\boldsymbol{X},t) = 0, \quad \text{on } \Gamma^{\text{3D}}_{N2} \times [0,T],\\ \end{align}\] where \(\Gamma^{\text{3D}}_i\) denote the Dirichlet boundaries at which the displacements \(\bar u_i\) are prescribed, and \(\Gamma^{\text{3D}}_{Ni}\) are the homogeneous Neumann boundaries with outward unit normal vector \(\boldsymbol{N}\). Note that the full boundary is composed of the Dirichlet and Neumann boundaries \(\Gamma^{\text{3D}} = \Gamma^{\text{3D}}_i \cup \Gamma^{\text{3D}}_{Ni}\). For generating the simulation data, we use an existing finite element method implementation. We refer to [68] for details. The simulations are performed on a tetrahedral mesh. The displacements on the surface, \(\Omega^{\text{2D}}\), can then be extracted on the corresponding triangular surface mesh, whose nodes coincide with a subset of the tetrahedral mesh nodes.

7 Admissible range and sampling of the invariants↩︎

The admissible region of the isochoric principal invariants has been studied by [95], see also the recent works by [87], [96], [97], who found that the admissible region is bounded by the uniaxial tension (UT) and equibiaxial tension (BT) curves in the invariants plane. Assuming incompressibility and uniaxial tension, the deformation gradient is \(\boldsymbol{F}=\text{diag}(\lambda,\lambda^{-1/2},\lambda^{-1/2})\) and consequently \(\boldsymbol{C}=\text{diag}(\lambda^2,\lambda^{-1},\lambda^{-1})\). The invariants compute to \[\label{eq:curve95UT} \bar I_1 = \lambda^2 + 2\lambda^{-1}, \quad \bar I_2 = 2 \lambda + \lambda^{-2}. \quad \text{(uniaxial tension/compression)}\tag{8}\] For incompressible equibiaxial tension, we have \(\boldsymbol{F}=\text{diag}(\lambda,\lambda,\lambda^{-2})\) and \(\boldsymbol{C}=\text{diag}(\lambda^2,\lambda^2,\lambda^{-4})\). The invariants compute to \[\label{eq:curve95EBT} \bar I_1 = \lambda^{-4} + 2\lambda^{2}, \quad \bar I_2 = 2 \lambda^{-2} + \lambda^4. \quad \text{(equibiaxial tension/compression)}\tag{9}\] It is easy to see that the possible deformations of uniaxial tension and equibiaxial tension, and hence also the corresponding invariants, are equivalent if we substitute \(\lambda=\hat{\lambda}^{-1/2}\) into the equibiaxial deformation. Hence, 8 9 describe the same curves in the invariants plane, which are parametrized by \(\lambda\). Specifically, the curve for uniaxial tension with \(\lambda \geq 1\) equals the curve for equibiaxial tension with \(\lambda \leq 1\), and vise versa. As discussed by the aforementioned works, these curves bound the admissible region of the invariants [87], [95]–[97].

The mapping \((\bar I_1, \bar I_2) \mapsto (\bar I_1 - 3, \bar I_2^{3/2} - 3^{3/2})\) preserves the topology of the admissible domain, and hence the curve defined by \[\bar I_1 = \lambda^2 + 2\lambda^{-1} - 3, \quad \bar I_2^{3/2} = [2 \lambda + \lambda^{-2}]^{3/2} - 3^{3/2},\] describes the bounds of the admissible region of the polyconvex invariants set \((\bar I_1 - 3, \bar I_2^{3/2} - 3^{3/2})\). We denote the admissible region by \(\mathcal{I}_{\text{adm}}\) and refer to 10 for an illustration.

The generation of training data for the neural operators requires sampling the invariants. In this work, we sample invariants in the domain \(\mathcal{I}_{\text{sample}} = \{ (\bar I_1^*, \bar I_2^*) \in \mathcal{I}_{\text{adm}} \, | \, \bar I_2^* \leq \bar I_1^* + 3 \}\). This sampling domain is chosen heuristically to capture the range of invariants observed in a representative simulation (see 10). We note that in our numerical experiments, different sampling strategies did not significantly affect the training results. The primary purpose of sampling the invariants is to provide a quantifiable measure of the mismatch between the true and predicted isochoric strain-energy density functions. Since these functions are defined on unbounded domains, there is no universal way to numerically quantify the discrepancy between them.

8 Taylor expansion of the strain-energy density↩︎

We recall \(I_1^*=\bar I_1 - 3\) and \(I_2^*=\bar I_2^{3/2} - 3^{3/2}\), and write the isochoric strain-energy density function as \(\bar W(I_1^*,I_2^*)\). Any general, sufficiently smooth function \(\bar W(I_1^*,I_2^*)\) can be approximated through a Taylor expansion about the undeformed state \((I_1^*,I_2^*)=(0,0)\) \[\begin{align} \bar W(I_1^*,I_2^*) =\;& \bar W(0,0) + I_1^*\, \bar W_{,I_1^*}(0,0) + I_2^*\, \bar W_{,I_2^*}(0,0) \\ &+ \frac{1}{2}\left[ [I_1^*]^2 \bar W_{,I_1^* I_1^*}(0,0) + 2 I_1^* I_2^* \bar W_{,I_1^* I_2^*}(0,0) + [I_2^*]^2 \bar W_{,I_2^* I_2^*}(0,0) \right] \\ &+ \frac{1}{6}\left[ [I_1^*]^3 \bar W_{,I_1^* I_1^* I_1^*}(0,0) + 3 [I_1^*]^2 I_2^* \bar W_{,I_1^* I_1^* I_2^*}(0,0) + 3 I_1^* [I_2^*]^2 \bar W_{,I_1^* I_2^* I_2^*}(0,0) + [I_2^*]^3 \bar W_{,I_2^* I_2^* I_2^*}(0,0) \right] \\ &+ \mathcal{O}\!\left(\|(I_1^*,I_2^*)\|^4\right). \end{align}\] Here, we used the comma convention for denoting partial derivatives. We note that \(\bar W(0,0) = 0\), and by introducing the respective constants \(C_{ij}\), this can be written as \[\begin{align} \bar W(I_1^*,I_2^*) =\;& C_{10} \, I_1^* + C_{01} \, I_2^* + C_{20} [I_1^*]^2 + C_{11} I_1^* I_2^* + C_{02} [I_2^*]^2 \\ &+ C_{30} [I_1^*]^3 + C_{21} [I_1^*]^2 I_2^* + C_{12} I_1^* [I_2^*]^2 + C_{03} [I_2^*]^3 + \mathcal{O}\!\left(\|(I_1^*,I_2^*)\|^4\right). \end{align}\] This is a modified version of the generalized Mooney-Rivlin model, which depends directly on the invariants \(I_1^*\) and \(I_2^*\) instead of \(I_1\) and \(I_2\). As higher-order contributions are generally negligible for most practical applications, only terms up to third order are considered for training data generation in the present work. Further, since we restrict our attention to separable strain-energy densities, we neglect coupling terms and consider functions of the form \[\bar W(I_1^*,I_2^*) = C_{10} \, I_1^* + C_{01} \, I_2^* + C_{20} [I_1^*]^2 + C_{02} [I_2^*]^2 + C_{30} [I_1^*]^3 + C_{03} [I_2^*]^3,\] in which we set \(C_{ij} \geq 0\).

9 Laplacian displacement field encoding↩︎

To construct a discretization-independent neural operator, we map the input displacement fields to a finite basis in the first layer of the branch net. We consider an interpolation with basis functions that have non-local support on the domain \(\Omega^{\text{2D}}\) such that the interpolation function fulfills the Dirichlet boundary conditions of the displacement fields. We additively split the displacement fields \[\label{eq:Laplacian95lifting} u_1(\boldsymbol{X},t) = u_1^{\text{lift}}(\boldsymbol{X},t) + u_1^{\text{hom}}(\boldsymbol{X},t), \quad u_2(\boldsymbol{X},t) = u_2^{\text{lift}}(\boldsymbol{X},t) + u_2^{\text{hom}}(\boldsymbol{X},t), \\\tag{10}\] into a priori chosen lifting fields \(u_i^{\text{lift}}(\boldsymbol{X},t)\), which fulfill the Dirichlet boundary conditions, and the remaining fields \(u_i^{\text{hom}}(\boldsymbol{X},t)\), which vanish on the Dirichlet boundaries. Specifically, for the geometry and loading conditions described in 2.4, we choose \(u_1^{\text{lift}}(\boldsymbol{X},t) = 0\) and \(u_2^{\text{lift}}(\boldsymbol{X},t) = \bar u_2(t) [1 - X_2]\), which fulfill the boundary conditions shown in 4.

It remains to interpolate the functions \(u_i^{\text{hom}}(\boldsymbol{X},t)\) \[\label{eq:Laplacian95interpolation95hom} u_1^{\text{hom}}(\boldsymbol{X},t) \approx \sum_{j=1}^{N_{\phi}} a_{1jt} \, \phi_{1j}(\boldsymbol{X}), \quad u_2^{\text{hom}}(\boldsymbol{X},t) \approx \sum_{j=1}^{N_{\phi}} a_{2jt} \, \phi_{2j}(\boldsymbol{X}), \\\tag{11}\] where \(a_{ijt}\) are interpolation coefficients and \(\phi_{ij}(\boldsymbol{X})\) are basis functions that vanish on the respective Dirichlet boundaries. We note that due to the lifting of the displacement fields, the basis functions are assumed to be independent of time.

While different choices of basis functions are possible and are expected to yield comparable neural operator accuracy, we focus here on a specific set of basis functions that proved particularly effective in our numerical studies. The basis functions are chosen as the eigenfunctions of the Laplace operator \[- \nabla \phi_{ij}(\boldsymbol{X}) = \lambda_{ij} \phi_{ij}(\boldsymbol{X}), \quad \text{(no Einstein summation)}\] subject to homogeneous Dirichlet conditions on the displacement-constrained boundaries and homogeneous Neumann conditions on the remaining parts of the boundary. To generate the basis functions, we discretize the eigenvalue problem above on a triangular finite element mesh. This yields the generalized eigenvalue problem \[\boldsymbol{K}\boldsymbol{\phi}^h_{ij} = \lambda_{ij} \boldsymbol{M}\boldsymbol{\phi}^h_{ij}, \quad \text{(no Einstein summation)}\] where \(\boldsymbol{\phi}^h_{ij} \in \mathbb{R}^{N_X}\) denotes the vector of nodal values of the basis function on the chosen reference mesh, and \(\boldsymbol{K}\) and \(\boldsymbol{M}\) denote the associated finite element stiffness and mass matrices. Since \(u_1\) and \(u_2\) generally do not share the same Dirichlet boundaries, we impose the homogeneous Dirichlet boundary conditions separately for each displacement component, yielding two generalized eigenvalue problems \[\boldsymbol{K}_1^{\text{free}} \boldsymbol{\phi}^h_{1j} = \lambda_{1j} \boldsymbol{M}_1^{\text{free}} \boldsymbol{\phi}^h_{1j}, \quad \boldsymbol{K}_2^{\text{free}} \boldsymbol{\phi}^h_{1j} = \lambda_{2j} \boldsymbol{M}_2^{\text{free}} \boldsymbol{\phi}^h_{2j}. \quad \text{(no Einstein summation)}\] Solving the eigenvalue problems finally yields the basis functions \(\boldsymbol{\phi}^h_{1j}\) and \(\boldsymbol{\phi}^h_{2j}\), from which we retain the first \(N_\phi\) eigenfunctions, corresponding to the lowest-frequency spatial modes. The first five eigenfunctions are illustrated in 11. These functions are used to interpolate the displacement fields, as shown in 11 . Since \(N_\phi < N_X\), in general, this reduces the dimensionality of the displacement data. For our training dataset, we have \(N_X=499\) and \(N_\phi=100\). We further expect that the interpolation improves the neural operators’ robustness to spatial noise, as small local oscillations in the displacement fields have only a minor effect on the interpolation coefficients.

10 Hyperparameters↩︎

The hyperparameters of the PANO and CANO architectures, together with the corresponding training settings, are summarized in 1. Most hyperparameters are chosen heuristically. However, for training the CANO, the number of neurons in the branch net, the initial learning rate, the weight decay, and the batch size, are optimized using a random search with the Python library optuna. Similarly, for training the PANO, the number of neurons in the branch and trunk nets, the number of latent features, the initial learning rate, the scheduler parameter T_0, the weight decay, and the batch size, are optimized with optuna. To select the optimal parameters, we only use the validation data, not the testing data. Finally, for the PANO, we observed that restarting the training process three times for the selected hyperparameters further reduced the validation loss.

References↩︎