Hyperreduction accelerates finite element model updating for the discovery of neural materials

Hagen Holthusen; Matthias Knipper; Kevin Linka

cs.CEOct 22, 2026

Hyperreduction accelerates finite element model updating for the discovery of neural materials

Authors: Hagen Holthusen, Matthias Knipper, Kevin Linka

Organizations: Computational Mechanics in Medicine, Institute of Applied Medical Engineering, RWTH Aachen University, Pauwelsstraße 20, 52074 Aachen, Germany

Abstract

Finite element model updating is among the most general approaches to discover constitutive models from structural experiments: since the displacement field is computed rather than measured, it accommodates arbitrary geometries, boundary conditions, and sparse or noisy data. However, every loss evaluation requires a nonlinear finite element solve per load step, and training neural materials requires many such evaluations. Intrusive model order reduction lowers this cost while retaining boundary conditions and variational structure. How can a reduced model remain accurate while the network weights change at every training step? We combine proper orthogonal decomposition and energy-conserving sampling and weighting with adjoint sensitivities of the reduced model, and adapt both the trial space and the sampled elements along the optimization path. The constitutive law is a polyconvex neural network: a logarithmic normalization term, together with a zero-anchored input convex architecture and a growth term, renders the energy polyconvex, coercive, and non-negative by construction, independent of the trained weights. We demonstrate the approach on two experimental datasets and one synthetic three-dimensional example. For a two-phase specimen, training time decreases from nine hours to 41 minutes compared with a full-order identification, without loss of accuracy. In the three-dimensional example, six modes replace 48 672 degrees of freedom, 116 of 15 625 elements suffice to evaluate the internal forces, and training takes 13 minutes. In all cases, the identified material reproduces the measured reaction force and predicts the displacement field from the force alone. Hyperreduction makes the discovery of neural materials from structural data affordable while retaining the advantages of finite element model updating.

Figures & tables

Category	Parameter	Case 1	Case 2	Case 3
Data	specimen	MNIST sample 22	TT1 [ 25 ]	synthetic block
	load steps $N_{\mathrm{t}}$	$21$	$51$ (DIC at $11$ )	$10$
	loss	force	force; force + displ.	force
	normalization of the force	target, $\bar{f}$	target and model, $\bar{f}$	target, $\bar{f}$
	$\bar{f}$	$15.94$ N	$236.9$ N	$85.61$ N
Discretization	element	hexahedron	hexahedron, Q1/P0	hexahedron

Table 1: Settings of the inverse identification for all three cases. In each case, the full-order solves (bootstrap, canary and final evaluation) use the same solver settings as the reduced runs. The reduction parameters of Case 1 define the base configuration, around which the study in Fig. 7 varies one parameter at a time.

Figure 1: Case 1, sample 22 of the Mechanical MNIST Challenge [ 92 , 93 ] in the reference configuration. The inclusion is shown in green and the matrix in grey. The plate of $40\times 40\times 2$ mm is discretized by $50\times 50\times 1$ hexahedral elements, whose nodes on the observed face coincide with the DIC points. On the top and bottom edges, the in-plane displacements $u_{x}(x,t)$ and $u_{y}(x,t)$ measured by DIC are prescribed for $N_{\mathrm{t}}=21$ load steps, and the out-of-plane displacement is fixed. The reaction force on the top edge is the only quantity that enters the loss.

Figure 2: Case 1, training history of the hyperreduced identification. (a) The loss decreases by about four orders of magnitude within the first

150

epochs and reaches

\mathcal{L}=2.50\cdot 10^{-4}

, compared with

4.09\cdot 10^{-4}

of the full-order reference. (b) The gradient norm drops from about

10^{2}

1.34\cdot 10^{-2}

; rebuilds cause no visible disturbance. (c) Wall time per epoch: reduced epochs take

1.2

1.6

s, canary epochs (orange markers) between

23

and

72

s, and canary epochs with a rebuild up to

114

s. (d) The cumulative wall-clock time of

0.69

h is

13

times smaller than the

8.95

h of the full-order reference (dashed). Grey dashed lines mark canary epochs and black dashed lines rebuilds of the reduced model.

Figure 3: Case 1, adaptivity of the reduction. (a) The snapshot matrix reaches its capacity of

N_{\mathrm{snap}}^{\max}=30

directly after the two bootstrap epochs, while the dimension of the trial space grows stepwise from

N_{\mathrm{r}}=9

N_{\mathrm{r}}=13

. (b) Training errors of POD and ECSW at each canary epoch; squares mark the POD error of the updated snapshots with the previous dimension, which exceeds

\varepsilon_{\mathrm{POD}}=10^{-4}

at four canary epochs and triggers the growth of

N_{\mathrm{r}}

. (c) The number of active ECSW elements grows from

186

356

, whereas the padded capacity remains at

384

, such that no change of the ECSW rule requires a recompilation. (d) Split of the wall-clock time: the

19

canary epochs consume

41\,\%

of the total time, almost as much as the

979

reduced epochs.

Panel of Figure 4 — Figure 4: Case 1, identified material evaluated by a full-order solve, compared with the measurement. (a) The machine curve is reproduced with a root-mean-square error of $0.249$ N, and (b) the force parity yields $R^{2}=0.9970$ . (c) The displacements, which are not part of the loss, agree with the DIC measurement with $R^{2}_{x}=0.909$ and $R^{2}_{y}=0.967$ . (d)–(e) Contour lines of the measured (black) and discovered (red) transverse and loading displacement nearly coincide away from the inclusion. (f) The relative displacement error remains below $5\,\%$ in most of the plate and concentrates along the upper flank of the inclusion (green outline) with up to $15\,\%$ . (g)–(h) Deformed configuration at the last load step: the inclusion moves almost rigidly and remains almost undeformed with $F_{yy}\approx 1.0$ , whereas the matrix above and below the inclusion is stretched up to $F_{yy}=1.87$ . Since only the reaction force enters the loss, this contrast is learned from the structural response alone.

Figure 5: Case 1, trial space at the end of the training. Top left: normalized singular values $\hat{\sigma}_{k}$ of the snapshot matrix and projection error $\mathrm{error}_{\mathrm{r}}(k)$ . After a decay of four orders of magnitude within eight modes, the spectrum forms a plateau, in which the projection error crosses $\varepsilon_{\mathrm{POD}}=10^{-4}$ at $N_{\mathrm{r}}=13$ . Remaining panels: magnitude (color) and direction (arrows) of the first five modes. Mode 1 describes the stretching in loading direction, mode 2 a lateral motion of the upper part, mode 3 the deformation around the upper part of the inclusion, mode 4 the region around its lower tip, and mode 5 boundary layers along the clamped edges.

Figure 6: Case 1, final ECSW rule on the deformed mesh. The remaining mesh is shown as wireframe and the inclusion in green. The $356$ active elements (about $14\,\%$ of the mesh) are colored by their weight; the $28$ padded entries of the set of $384$ carry zero weight and add no element. The element rows along both clamped edges are selected completely with weights close to one, and further elements concentrate in and around the inclusion, whereas the matrix away from the inclusion is represented by few elements with weights of up to about $220$ .

Figure 7: Case 1, one-at-a-time variation of the reduction parameters around the base configuration of Table 1 . Shown are the force error of the identified material (blue, left axis) and the speed-up per epoch against the full-order reference (red, right axis); the dotted line marks the force error of the full-order reference. (a) Canary interval

N_{\mathrm{canary}}

: only the largest interval deteriorates the accuracy. (b)–(c) Tolerances of the model selection: no influence in the tested range. (d) Snapshot update

\varepsilon_{\mathrm{proj}}

: stricter tolerances improve the accuracy at constant cost. (e) POD tolerance

\varepsilon_{\mathrm{POD}}

and (f) ECSW tolerance

\varepsilon_{\mathrm{ECSW}}

: the speed-up varies between

8.6

and

60

, while the force error remains close to the base configuration.

Figure 8: Case 1, comparison of all configurations of the parameter study. (a) Force parity: all reduced runs (squares) and the full-order reference (circles) reproduce the measured force nearly identically. (b) Force error against the speed-up per epoch: all reduced runs except

\varepsilon_{\mathrm{proj}}=10^{-2}

reach an error at or below that of the full-order reference (dotted lines) at speed-ups between

8.6

and

60

. (c) Median epoch times: reduced epochs require

0.5

4

s, canary epochs between

25

and

85

s. (d) Final dimension of the trial space

N_{\mathrm{r}}

and size of the padded ECSW set

|\tilde{\mathcal{E}}|

, compared with the

n_{\mathrm{el}}=2\,500

elements of the mesh; only the tolerances of POD and ECSW change these sizes considerably.

Figure 9: Case 2, region of interest of the perforated specimen TT1 of Abbasi et al. [25] in the reference configuration. The membrane of $160\times 50\times 2.5$ mm with six elliptical holes is discretized by $1\,518$ hexahedral elements with one element across the thickness. On the top and bottom edges, B-spline fits of the measured in-plane displacements $u_{x}(x,t)$ and $u_{y}(x,t)$ are prescribed for $N_{\mathrm{t}}=51$ load steps, and the out-of-plane displacement is fixed. The reaction force on the bottom edge enters both losses, the DIC displacements at $11$ load steps only the second one.

Figure 10: Case 2, training history of the force-based identification on TT1. (a) The loss decreases by three orders of magnitude within

100

epochs and reaches

\mathcal{L}=6.0\cdot 10^{-5}

. (b) The gradient norm shows short peaks at canary epochs, when the ECSW rule is updated, which decay within a few epochs. (c) Reduced epochs take about

11

s, canary epochs between

123

and

169

s. (d) The training requires

2.73

h in total. Grey dashed lines mark canary epochs and black dashed lines rebuilds.

Figure 11: Case 2, adaptivity of the reduction in the force-based identification. (a) The snapshot matrix is full from the first bootstrap epoch on, since one epoch contains

51

states, and the trial space grows once from

N_{\mathrm{r}}=26

27

. (b) After epoch

300

, the projection error of the snapshots drops to about

10^{-11}

, as the snapshots of the converging material become nearly linearly dependent. (c) ECSW requires

613

801

of the

1\,518

elements, stored in a padded set of

1\,024

entries. (d) The reduced epochs dominate the wall-clock time with

79\,\%

Figure 12: Case 2, identification from force and displacements. Top: training history; the total loss (solid) is dominated by the displacement term, which levels off at about

1.1\cdot 10^{-1}

, while the force term (dashed blue) falls below

10^{-3}

. The jump of the displacement term at the first canary epoch shows that the outdated trial space had underestimated the displacement misfit. Bottom: adaptivity of the reduction; the trial space grows from

N_{\mathrm{r}}=25

27

at epochs

52

and

102

, and ECSW requires up to

928

active elements. The total wall-clock time of

2.78

h is nearly identical to the force-based run.

Panel of Figure 13 — Figure 13: Case 2, material identified from the force alone, evaluated by a full-order solve on the training specimen TT1 (top) and on the unseen specimen TT2 (bottom). Top: (a)–(b) The machine curve of TT1 is reproduced with a root-mean-square error of $1.80$ N and $R^{2}=0.9992$ . (c) The displacements, which are not part of the loss, agree with $R^{2}_{x}=0.933$ and $R^{2}_{y}=0.997$ . (d)–(e) Measured (black) and discovered (red) contour lines nearly coincide in loading direction and deviate slightly in transverse direction. (f) The relative displacement error remains below $6\,\%$ and concentrates around the largest hole and at the right edge. Bottom: For TT2, the material overestimates the force at large displacements by about $10\,\%$ (root-mean-square error of $16.4$ N), in line with the full-order results of [ 40 ] , and predicts the displacement field with $R^{2}_{x}=0.892$ and $R^{2}_{y}=0.997$ . The corresponding results for the material identified from force and displacements are shown in Fig. 1 .

Figure 14: Case 2, trial space at the end of the force-based identification. Top left: normalized singular values

\hat{\sigma}_{k}

and projection error

\mathrm{error}_{\mathrm{r}}(k)

; the spectrum decays slowly over almost twenty modes, and the dimension

N_{\mathrm{r}}=27

, set by

\varepsilon_{\mathrm{POD}}=10^{-6}

earlier in the training, is never decreased. Remaining panels: magnitude (color) and direction (arrows) of the first five modes; mode 1 describes the stretching in loading direction, while the higher modes resolve the motion of individual ligaments between the holes.

Figure 15: Case 2, final ECSW rule of the force-based identification on the deformed mesh of TT1. The $801$ active elements ( $53\,\%$ of the mesh) are colored by their weight, the remaining mesh is shown as wireframe. They line both clamped edges and the boundaries of all holes, where the deformation localizes. With weights of at most $17$ , the hyperreduction is considerably less sparse than for the plate of Case 1.

Panel of Figure 16 — Figure 16: Case 3, cube with two spherical inclusions in the reference configuration. The cube of edge length $10$ mm contains two spheres of radii $R_{1}=1.5$ mm and $R_{2}=2.0$ mm (green) and is discretized by $25^{3}$ hexahedral elements. Top: three-dimensional view; the bottom face is clamped, and the top face is displaced by $\bar{u}_{y}(t)$ up to $3$ mm in loading direction $y$ in $N_{\mathrm{t}}=10$ load steps, with the transverse displacements fixed. Bottom: front view (x–z), top view (x–y) and side view (y–z) with dimensions in mm. The inclusions are ten times stiffer than the matrix.

Figure 17: Case 3, training history of the hyperreduced identification. (a) The loss decreases from about

60

\mathcal{L}=4.2\cdot 10^{-5}

, with a steep drop around epoch

190

. (b) The gradient norm decreases from about

5\cdot 10^{2}

1.1\cdot 10^{-2}

. (c) Reduced epochs take

0.17

0.3

s, canary epochs about

40

s and canary epochs with a rebuild about

85

s, i.e., a reduced epoch is more than two orders of magnitude cheaper than a canary epoch. (d) The training requires

0.21

h in total, most of which is spent in the full-order epochs.

Figure 18: Case 3, adaptivity of the reduction. (a) With ten load steps per epoch, the snapshot matrix reaches its capacity of

30

only at epoch

252

; the trial space grows from

N_{\mathrm{r}}=4

6

. (b) At epochs

102

and

302

, the POD error with the previous dimension exceeds

\varepsilon_{\mathrm{POD}}=10^{-4}

(squares). (c) At epoch

252

, the active ECSW set exceeds the padded capacity of

96

elements, which is then raised to

192

; this rebuild changes no dimension of the trial space. (d) Bootstrap and canary epochs account for

85\,\%

of the wall-clock time, the

391

reduced epochs for

15\,\%

Panel of Figure 19 — Figure 19: Case 3, identified material compared with the synthetic reference. (a)–(b) The force is reproduced with a root-mean-square error of $0.556$ N and $R^{2}=0.9995$ . (c) All three displacement components, which are not part of the loss, agree with $R^{2}\geq 0.990$ . (d)–(e) Reference (black) and discovered (red) contour lines of the transverse and loading displacement in the mid-plane $y=5$ mm nearly coincide, also around both inclusions. (f) The relative displacement error remains below $1.5\,\%$ . (g)–(h) Deformed configuration at the last load step: the displacement in loading direction increases from the clamped bottom face to the top face, which is displaced by $3$ mm, and $F_{yy}$ varies between $1.044$ and $1.61$ on the surface, with the largest stretches in the vicinity of the clamped faces.

Figure 20: Case 3, energy density of the identified (red) and the reference Neo-Hookean material Eq. 44 (black) at a single material point. Top: matrix; bottom: inclusion; columns: uniaxial tension, simple shear and equibiaxial tension, each over the range the phase experiences in the reference solve, with the deformation gradients Eq. 46 ; in tension, the unprescribed diagonal components follow from vanishing normal stresses. The matrix law is reproduced in uniaxial tension (root-mean-square error of

0.004

MPa), which dominates the response of the block, whereas it is too stiff in simple shear (

0.048

MPa) and in equibiaxial tension (

0.44

MPa), which the experiment hardly activates. The stiffer and less deformed inclusion is identified less accurately in all three load cases (

0.013

0.073

and

0.069

MPa).

Figure 21: Case 3, trial space at the end of the training. Top left: normalized singular values $\hat{\sigma}_{k}$ and projection error $\mathrm{error}_{\mathrm{r}}(k)$ ; in contrast to the experimental cases, the spectrum decays steadily, and $\varepsilon_{\mathrm{POD}}=10^{-4}$ is reached with $N_{\mathrm{r}}=6$ modes. Remaining panels: magnitude of the first five modes on the deformed block; mode 1 describes the uniaxial stretching, the higher modes the bulging of the lateral faces and the deformation at the edges of the clamped faces.

Figure 22: Case 3, final ECSW rule. The mesh is shown as wireframe and the inclusions in green. Only $116$ of the $15\,625$ elements carry a positive weight. They cluster inside and around both inclusions and at the corners of the clamped faces, whereas the matrix is represented by a few elements with weights of up to about $1\,200$ . The $76$ padded entries carry zero weight and add no element to the figure.

Appendix figures & tables1 asset

Supplementary material from the paper’s appendix.

Appendix

Panel of Figure 1 — Figure 1: Case 2, material identified from force and displacements, evaluated by a full-order solve on the training specimen TT1 (top) and on the unseen specimen TT2 (bottom). Top: (a)–(b) The force is reproduced less accurately than in Fig. 13 , with a root-mean-square error of $6.47$ N and $R^{2}=0.990$ ; the material is too soft at small and too stiff at large displacements. (c)–(f) In turn, the displacement field agrees with $R^{2}_{x}=0.991$ and $R^{2}_{y}=0.998$ , and the relative error falls below $2\,\%$ in most of the specimen. The trade-off between both channels indicates that force and field are not fully consistent within the present model. Bottom: For TT2, the force at large displacements is overestimated by about $10\,\%$ (root-mean-square error of $16.5$ N), as for the force-based material, whereas the displacement field benefits from the additional data: $R^{2}_{x}$ increases from $0.892$ to $0.998$ , and the relative error drops below $5\,\%$ almost everywhere.

Finite Element Method

Finite Element Method (FEM) approximates physical fields by solving equations over small mesh elements.

See all 31 papers

May 22, 2026

Automatic-Differentiation-Enabled Finite Element Model Updating: Learning Constitutive Neural Networks from Structural Deformation Data

The identification of constitutive neural network models from heterogeneous full-field deformation data provides a robust alternative to traditional calibration methods based on homogeneous stress-strain experiments, particularly given the high dimensionality of trainable parameters. This work presents ADE-FEMU (Automatic-Differentiation-Enabled Finite Element Model Updating), an end-to-end differentiable framework that integrates a constitutive neural network model within a JAX-FEM nonlinear solver and identifies its parameters through gradient-based minimization of a measurement-mismatch loss, allowing for training on structural data with limited kinetic information. Thus, it can be considered as an extension of finite element model updating (FEMU) for hyperelastic materials, implemented within a GPU-accelerated, fully automatically differentiable finite element framework Newton tangent stiffness and loss gradients are computed automatically using forward- and reverse-mode automatic differentiation throughout the entire pipeline, thereby removing the need for analytic adjoints or offline surrogate models. ADE-FEMU is demonstrated for two architectures: a grey-box Constitutive Artificial Neural Network (CANN), a polyconvex, fully connected model with high flexibility, and a white-box CANN, an expert-system network with phenomenologically interpretable strain-energy terms. Focusing on incompressible isotropic hyperelasticity, ADE-FEMU is evaluated on three open experimental datasets: (1) full digital image correlation (DIC) of a perforated tensile specimen, (2) a reduced-data scenario with a one-dimensional stretch profile and global force-displacement curve, and (3) a heterogeneous matrix-inclusion system in which both phases constitutive laws are identified and generalized to twenty-two previously unseen samples.

Matthias Knipper, Chenyi Ji, Malte Brand +1

Computational Mechanics in Medicine, Applied Medical Engineering, RWTH Aachen University Pauwelsstraße 20, 52074 Aachen, Germany · Institute for Continuum and Material Mechanics, Hamburg University of Technology, Eißendorfer Straße 42, 21073 Hamburg, Germany

Automatic Differentiation Hyperelasticity Materials Science Physics-Informed ML Finite Element Method

Ask PDF

Mar 27, 2026

Hyperelastic constitutive model discovery with differentiable finite elements and structure-preserving neural networks

The discovery of constitutive laws from experimentally accessible measurements is a central problem in nonlinear computational mechanics. Many data-driven constitutive identification approaches rely either on paired strain-stress data or on full-field displacement measurements, both of which are difficult to obtain in realistic three-dimensional settings. We present a differentiable finite element framework for the discovery of hyperelastic material laws from partial observations, including boundary-only displacement measurements and global reaction forces. The method embeds the nonlinear finite element equilibrium problem directly into the learning loop, so that candidate strain-energy densities are assessed through the deformation fields and reactions they induce. This formulation enforces mechanical equilibrium as a constraint and allows the loss function to be evaluated only at observed locations. To ensure physical admissibility and promote numerical solvability throughout training, the constitutive response is represented by Hyperelastic Neural Networks, a structure-preserving neural class that enforces residual energy and stress-free conditions, frame indifference, isotropic material symmetry, polyconvexity, coercivity, and controlled volumetric growth by construction. The resulting PDE-constrained learning problem is solved using a quasi-Newton strategy combined with continuation and solver-aware backtracking. Numerical experiments in two- and three-dimensional finite elasticity demonstrate accurate recovery of hyperelastic isotropic responses from boundary-only data, robustness to measurement noise, and generalization across geometries, loading conditions, and boundary conditions.

Francesco Regazzoni

MOX, Department of Mathematics, Politecnico di Milano, Milan, Italy

Structure-Preserving Neural Networks Hyperelasticity Physics-Informed ML Differentiable Physics

Ask PDF

Jun 29, 2026

Implementation of Hyperelastic Physics-Augmented Neural Networks in the Explicit Finite Element Codes Simcenter Radioss and OpenRadioss with Applications to Impact Events

Data-driven material modeling techniques have gained significant attention due to their ability to capture complex constitutive behaviors beyond the limitations of classical material models. Physics-augmented neural networks (PANNs), which embed physical constraints directly into their architecture, combine the flexibility of machine learning with the reliability required for engineering simulations. This work presents an approach to integrate such network architectures into the explicit finite element solvers Simcenter Radioss and OpenRadioss (Siemens). A framework for transferring pretrained network architectures and their parameters to a standalone user material routine is developed. Networks are trained using PyTorch, though the procedure can be adapted to other frameworks such as TensorFlow, enabling the use of PANNs within existing finite element technology without requiring specialized solvers. Particular emphasis is placed on computational efficiency. The influence of network architecture on simulation performance is investigated, and strategies for reducing evaluation costs while preserving accuracy are discussed. Specifically, replacing the SoftPlus activation function with SQuarePlus is shown to reduce computational cost. A publicly available GitHub repository automates the generation of Fortran user material routines, requiring only the specification of the network architecture and trained parameters. An example impact simulation demonstrates that the generated PANN user material reproduces the nonlinear behavior characteristic of hyperelastic materials under large strains, providing a practical route toward machine-learning-based constitutive models in explicit finite element simulations.

Lukas Maurer, Sascha Eisenträger, Marian Bulla +1

Institute of Materials, Technologies and Mechanics Otto von Guericke University Magdeburg Magdeburg, D-39106 · Siemens Digital Industries Am Kabellager 9 Köln, D-51063

Neural Surrogate Modeling Efficient Neural Network Inference Hyperelasticity Physics-Informed ML Finite Element Method

Ask PDF

Abstract

Finite Element Method

Explore similar work

Automatic-Differentiation-Enabled Finite Element Model Updating: Learning Constitutive Neural Networks from Structural Deformation Data

Hyperelastic constitutive model discovery with differentiable finite elements and structure-preserving neural networks

Implementation of Hyperelastic Physics-Augmented Neural Networks in the Explicit Finite Element Codes Simcenter Radioss and OpenRadioss with Applications to Impact Events