7. Material Identification

Finite element model updating, the virtual fields method, DIC full-field measurement, sensitivity and identifiability, hyperelastic soft materials, constitutive model discovery

1. Why Inverse Identification: From Standard Tests to Material Testing 2.0

Every forward simulation in this guide silently assumed that somebody, somewhere, measured the material. Given constitutive parameters $\boldsymbol{\theta}$, a geometry, and boundary conditions, the finite element method predicts the displacement field. This module is about the day you are that somebody: the material in front of you — say, the silicone skin you cast for a vision-based tactile sensor — has no datasheet worth trusting, and the parameters must come out of an experiment.

Material Testing 1.0: Engineer the Stress State Away

The classical approach never inverts anything. It engineers the experiment so that the stress state is statically determinate and homogeneous: in a dogbone uniaxial test, $\sigma = F/A$ everywhere in the gauge section, strain comes from an extensometer, and the stress–strain curve is read off directly. The price is that each such test probes exactly one loading mode — as the load rises it traces a one-dimensional family of stress states, every one of them proportional to the same uniaxial state, so no amount of pulling harder adds a genuinely new stress direction. An anisotropic law therefore does need a battery of separate tests — off-axis tension, shear rails, biaxial rigs — because a single direction never reaches the off-axis stiffnesses at all. A multi-parameter isotropic law is the weaker case, and the word to avoid is "needs": several parameters can enter one loading path with different amplitude dependence, so extra loading modes usually improve the conditioning and supply something to validate against, rather than being strictly necessary (§8 works this through for Mooney–Rivlin). And soft or small specimens rarely achieve homogeneity at all: grip effects, self-weight sag, bulging.

Material Testing 2.0: Invert the Philosophy

Material Testing 2.0 (Pierron & Grédiac 2021) turns the logic around. Deliberately design a heterogeneous test so that a single experiment activates many stress states at once; measure the full displacement field with a camera — $10^4$–$10^6$ data points per frame via digital image correlation — and extract the parameters simultaneously with an inverse method — all of them that the design actually excites separably, which is a real and checkable condition, not a formality (§8). Pierron's 2023 brief review surveys how far this program has come and what still blocks routine adoption.

Two inverse families dominate. Finite Element Model Updating (FEMU) iterates a forward FE model of the test until it reproduces the measurements; the idea dates back to Kavanagh & Clough (1971), who characterized elastic solids by matching FE predictions to measured response. The Virtual Fields Method (VFM) — introduced by Grédiac in 1989 and developed comprehensively in Pierron & Grédiac's 2012 book — bypasses the forward solve entirely by writing the principle of virtual work with the measured strain field.

For this site's audience the motivation is concrete: characterizing the silicone skin of a vision-based tactile sensor is precisely this problem, and the camera infrastructure is already sitting in the sensor. The module arc: measurement (DIC, §2) → problem statement and method map (§3) → FEMU (§4) → VFM linear (§5) → VFM nonlinear (§6) → soft materials (§7) → identifiability and test design (§8) → honest comparison (§9) → model discovery (§10) → lab workflow (§11).

2. The Data: Full-Field Measurement by Digital Image Correlation

Everything downstream is only as good as the measured field, so we start with the measurement. Digital image correlation (DIC) is the workhorse; the standard reference is the book by Sutton, Orteu & Schreier (2009), and Hild & Roux (2006) review the path from displacement measurement to elastic-property identification. The principle: apply a random speckle pattern to the specimen, photograph it before and during loading, and assume gray-level conservation between the reference image $f$ and the deformed image $g$:

Subset Matching

Around each point of a regular grid, take a subset — typically a 21–41 px square — and parametrize its motion with a shape function, usually affine (translation plus gradient terms). The subset's displacement is found by minimizing a correlation criterion; the zero-normalized sum of squared differences (ZNSSD) is standard because it is insensitive to uniform lighting offset and scale changes:

Subpixel accuracy comes from gray-level interpolation (bicubic splines): displacement resolution on the order of 0.01 px is routine, around 0.005 px under good conditions. The step size — the grid spacing, often 1/3 to 1/2 of the subset size — sets the data density.

From Displacements to Strains: the Virtual Strain Gauge

Identification needs strains, and strain is a derivative of noisy displacement data:

The strain measure written here is the infinitesimal one, and it is written only to make the point about differentiating noisy data. It is not objective at the silicone stretches of §7: for finite deformation the kinematic object is the deformation gradient $\mathbf{F} = \mathbf{I} + \nabla_0\mathbf{u}$ (§6), with $\mathbf{E} = (\mathbf{F}^{\mathsf{T}}\mathbf{F} - \mathbf{I})/2$ if a strain tensor is wanted. The noise scalings are properties of the differentiation and carry over to any of these measures. They are also the scalings of two different estimators, not universal bounds. Note the grid spacing $h$ in the first: without it, $\sigma_u/L_{\mathrm{VSG}}^{3/2}$ does not even have the dimensions of a strain.

Differentiation amplifies noise. In practice strain is therefore computed by fitting a local polynomial over a window of grid points — the virtual strain gauge (VSG) of size $L_{\mathrm{VSG}}$. The exponent is worth pinning down, because the two estimators diverge a lot at large windows. A crude two-point difference $\Delta u / L$ differences two independent noisy endpoints, giving $\sigma_{\varepsilon} \sim \sqrt{2}\,\sigma_u/L_{\mathrm{VSG}}$: it uses only the endpoints, so widening the window buys nothing but a longer baseline. A least-squares slope over all the points in the window also averages them, and for uncorrelated displacement noise on a grid of fixed spacing $h$ its exact standard deviation is $\sigma_u\sqrt{12/(h^2 w(w^2-1))}$ with $w = L_{\mathrm{VSG}}/h + 1$ points — which for $w \gg 1$ becomes $\sqrt{12h}\,\sigma_u/L_{\mathrm{VSG}}^{3/2}$, falling as $L_{\mathrm{VSG}}^{-3/2}$. Real DIC noise is spatially correlated, which eats most of that extra averaging gain and pushes the practical behaviour back toward $1/L_{\mathrm{VSG}}$; treat the pair as the two limiting estimator scalings rather than as bounds, and measure where your own setup lands. That window is the fundamental trade of the whole measurement chain: a large window averages the noise down but smears gradients, biasing the strain near holes and notches — exactly where a well-designed heterogeneous test puts its information. A small window follows the gradients faithfully but is noisy. Typical engineering strain noise floors sit in the $10^{-4}$–$10^{-3}$ range depending on smoothing — treat these as typical orders, and measure the actual floor for your setup from static image pairs.

2D, Stereo, and Volumetric Variants

Stereo-DIC uses two calibrated cameras to measure the 3D surface displacement field. It is the safe default for soft materials, because a single-camera 2D-DIC setup converts out-of-plane motion $\Delta z$ into fake in-plane strain of order $\Delta z / z$, with $z$ the standoff distance — and soft specimens bulge, sag and drift out of plane readily. Whether that matters is a quantitative question, not a property of the material class: compare $\Delta z / z$ against your target strain resolution. Digital Volume Correlation (DVC) extends the idea to volumetric images (µCT, confocal) and is the relevant variant for tissue.

Practicalities and Systematic Errors

• Speckle: feature size 3–5 px at the working magnification, roughly 50% coverage.
• Stereo calibration: use a calibration target, calibrate immediately before testing.
• Systematics: gray-level interpolation bias, lens distortion, and pattern degradation at large stretch on elastomers — speckle patterns crack, and a continuous paint film can locally stiffen a thin soft specimen (more in §11).
• Noise floor: acquire static image pairs and measure it; never assume it.

3. Problem Statement and the Map of Identification Methods

Formalize the problem. Given: a domain $\Omega$ (the imaged region); measured displacements $\mathbf{u}_{\mathrm{m}}$ at the DIC grid points, hence strains $\boldsymbol{\varepsilon}_{\mathrm{m}}$; and boundary information that is typically incomplete — a load cell gives the resultant $F$ on the grips, but never the traction distribution. Find $\boldsymbol{\theta} \in \mathbb{R}^p$ parametrizing the constitutive law:

Any solution must reconcile three ingredients: kinematic compatibility (supplied by the measurement), equilibrium ($\nabla\!\cdot\!\boldsymbol{\sigma} + \mathbf{b} = \mathbf{0}$, or its weak form), and the constitutive relation (the unknown). Identification methods differ in which residual $\mathcal{R}$ they minimize and whether a forward boundary-value problem is solved along the way.

The Five Methods

The authoritative overview by Avril et al. (2008) catalogues the full-field identification methods side by side. In compressed form:

What Avril et al. (2008) actually do deserves stating precisely, because the paper is the field's common reference point. All five methods are formulated side by side, and identifications are then run on full-field experimental data from four different experiments: (i) a tensile test, (ii) the Brazilian test, (iii) a shear-flexural test, and (iv) a biaxial test featuring a non-uniform damage field — hence non-uniform equivalent elastic properties. Tests (i)–(iii) identify uniform anisotropic elastic properties; tests (ii)–(iv) involve non-uniform strain fields in the region of interest. What the paper does not do is run all five methods on all four tests: the methods are formulated in common and then exercised where each is applicable, so the results are not a five-way head-to-head on a shared benchmark. (Consistent with that — and offered here as this module's reading rather than a claim about the paper's wording — RGM needs displacements and tractions over the boundary, which a DIC surface field plus a scalar load resultant does not supply; see the table above.) Read it for how the methods are formulated and exercised on real data, not as a single accuracy ranking.

The practical selection rule — offered as this module's synthesis, not as a headline verdict of that paper: choose on data type, computational budget, and how much of the boundary you actually know. Full coverage of a thin specimen and a linear-in-parameters law favor VFM; partial data, complex 3D geometry, or an exotic law favor FEMU; complete boundary knowledge is rare enough that RGM is the specialist's tool.

When the Unknown Is a Material Field: Elastography

Nothing in the formulation forces $\boldsymbol{\theta}$ to be a handful of numbers. The unknown can be a spatially varying field $\boldsymbol{\theta}(\mathbf{x})$ — a stiffness map rather than a stiffness value. Test (iv) of Avril et al. (2008) is exactly this situation: a non-uniform damage field means the equivalent elastic properties vary over the region, and the identification returns a property distribution. Pushed to its extreme, this is the defining problem of elastography: reconstruct a modulus map from a measured interior kinematic field, because stiffness contrast is the quantity of interest — the flagship application being medical imaging of tissue.

The same five residuals apply; only the parameter space explodes — from $p$ numbers to one or more unknowns per pixel or element. The problem then inherits everything Module 3 warned about: severe ill-posedness, regularization that trades resolution against noise amplification, and unobservable components (the scale indeterminacy above now holds per connected region unless force data or a reference stiffness pins it). FEMU generalizes by letting $\boldsymbol{\theta}$ vary element-wise, at the price of a sensitivity computation that grows with the field's dimension; the gap methods reconstruct contrast maps more directly from the measured field. Module 9 hands its configuration-side problems over to this section whenever the unknown turns out to be the material rather than the geometry.

4. FEMU: Finite Element Model Updating

FEMU is the oldest and most flexible approach (Kavanagh & Clough 1971): build a forward FE model of the test, compare its predictions to the measurements, and update $\boldsymbol{\theta}$ until they agree. Everything you know about forward FEM (Module 2) is reused; the new ingredients are a cost function, a sensitivity computation, and an update rule.

The Forward Model and Its Boundary Conditions

Geometry comes from the specimen. Boundary conditions come in two flavors: (a) measured resultant loads applied to the modeled grips, or (b) measured displacements imposed on the boundary of the DIC region itself. Option (b) is elegant — it removes grip modeling entirely — but it injects DIC noise into the BCs and, crucially, removes the stiffness-scale information: a displacement-driven model produces the same displacement field for any global modulus scaling. A force residual must then be added to restore the scale (§3's indeterminacy, again).

Cost Function

The first term is FEMU-U (displacement residual, with FE displacements interpolated to the DIC points), the second FEMU-F (reaction-force residual); the standard formulation weights and combines both. Statistically optimal weighting is the inverse noise covariance (generalized least squares); in practice, diagonal weights normalized by the noise variance of each data type.

Sensitivities: the Crux

Minimization is Gauss–Newton, driven by the sensitivity matrix

That displayed update minimizes the displacement term of $J$ alone. To descend the combined cost above, stack the force channel into both the residual and the sensitivity — $\mathbf{r} = \big[\mathbf{u}_{\mathrm{m}} - \mathbf{u}(\boldsymbol{\theta});\ \sqrt{w_F}\,(F_{\mathrm{m}} - F(\boldsymbol{\theta}))\big]$ and $\mathbf{S} = \big[\partial\mathbf{u}/\partial\boldsymbol{\theta};\ \sqrt{w_F}\,\partial F/\partial\boldsymbol{\theta}\big]$ — and the same formula then applies with $\mathbf{W}_u$ replaced by the conforming block-diagonal $\widetilde{\mathbf{W}} = \mathrm{diag}(\mathbf{W}_u,\, 1)$. It cannot be applied unchanged: $\mathbf{W}_u$ is sized for the displacement block, so $\mathbf{S}^{\mathsf{T}}\mathbf{W}_u\mathbf{S}$ is not even defined once the force row is stacked in. Equivalently, whiten the displacement block with $\mathbf{W}_u^{1/2}$ — the $\sqrt{w_F}$ factors already whiten the force row — and use the unweighted augmented update. The distinction is not cosmetic here: under measured-displacement BCs the force residual is the only term carrying stiffness-scale information, so a displacement-only update cannot converge on the scale no matter how many iterations it runs.

How $\mathbf{S}$ is computed decides both cost and reliability:

• Forward finite differences: $(p+1)$ FE solves per iteration; step-size dependent, and the step-size error interacts with noise.
• Direct differentiation: differentiate the discrete equilibrium — $\mathbf{K}\mathbf{u} = \mathbf{f}$ when it is linear, the residual form $\mathbf{R}(\mathbf{u}, \boldsymbol{\theta}) = \mathbf{0}$ when it is not — to get one extra linear solve per parameter, reusing the factorization already on hand (collapsible below).
• Adjoint methods: the winner when parameters are many (the elastography limit of §3; Module 10 owns the machinery).
• Automatic differentiation: Kumar, Seidl, Granzow, Yang & Fuhg (2025) derive numerically-exact AD sensitivities for finite-strain elastoplastic FEMU and VFM — and argue that earlier method comparisons were limited by their reliance on approximate gradients.

When $\mathbf{S}^{\mathsf{T}}\mathbf{W}_u\mathbf{S}$ is ill-conditioned — near-collinear parameters, flat directions — add Levenberg–Marquardt damping. For elastic problems from a reasonable initial guess, convergence is fast in practice: on the order of a few to a few tens of iterations. Treat that as a rough expectation, not a specification.

Read the Residual Field, Not Just the Minimum

Strengths: arbitrary laws, arbitrary geometry, tolerates partial and sparse data. Weaknesses: the cost of repeated (possibly nonlinear) FE solves, local minima, and sensitivity to the boundary-condition assumptions baked into the forward model.

5. The Virtual Fields Method: Linear Case

The VFM (Grédiac 1989; comprehensive treatment in Pierron & Grédiac 2012; review in Grédiac, Pierron, Avril & Toussaint 2006) starts from a statement your forward-FEM training already contains — the principle of virtual work — and weaponizes the fact that DIC hands you the actual strain field everywhere.

Here is the mechanism. DIC supplies $\boldsymbol{\varepsilon}(\mathbf{x})$ at every point of the region. For linear elasticity, $\boldsymbol{\sigma}$ depends linearly on the stiffness parameters. Substitute both into the PVW: each chosen virtual field converts the equation into one scalar linear equation in the unknown stiffnesses. Collect $N_v \ge p$ fields and you have

solved directly — square when $N_v = p$, but invertible only if $\operatorname{rank}\mathbf{A} = p$; overdetermined least squares when $N_v > p$. No FE model, no iteration — milliseconds. The condition to check is $\operatorname{rank}\mathbf{A} = p$, and it is not implied by choosing $p$ functionally independent virtual fields: what enters $\mathbf{A}$ is each field contracted against the measured strain field, and two perfectly distinct fields can still produce dependent rows. The worked example below exhibits exactly that — two different virtual fields, both blind to the same stiffness. Check the rank (or the condition number, §8) of the assembled $\mathbf{A}$, never the independence of the fields you wrote down.

The Canonical Worked Example

A rectangular coupon of area $S$ and thickness $t$, clamped at $x = 0$, with a transverse resultant $F$ applied at $x = L$. Take the virtual field $\mathbf{u}^{*} = (0,\ x)$. Its virtual strains: $\gamma^{*}_{xy} = 1$, all normal components zero. Now watch the boundary terms do their work:

• On the clamped edge, $\mathbf{u}^{*} = \mathbf{0}$: the unknown clamp reactions do zero virtual work.
• On the loaded edge, $u^{*}_y = L$ is constant: the unknown traction distribution enters only through its resultant, and the external work is exactly $F \cdot L$.

The PVW collapses to $t \int_S \sigma_{xy}\, \mathrm{d}S = F L$, and with $\sigma_{xy} = Q_{ss}\,\gamma_{xy}$ the shear modulus follows in closed form from an area integral of the measured shear strain:

One virtual field, one parameter, no FE forward solve — the mechanical model (equilibrium, plane stress, the constitutive form) is still very much present. For the orthotropic plane-stress case, additional independent fields are needed — and it pays to check, component by component, which stiffness each one actually reaches. Take $\mathbf{u}^{*} = (xy,\ 0)$: it gives $\varepsilon^{*}_{xx} = y$ and $\gamma^{*}_{xy} = x$, so it probes $Q_{xx}$, $Q_{xy}$ and $Q_{ss}$ — but $\varepsilon^{*}_{yy} = 0$, and the $Q_{yy}$ term in the work equation below is multiplied by exactly that, so the field is blind to $Q_{yy}$. Check its boundary behavior too, because it is not innocent: at the loaded edge $\mathbf{u}^{*} = (Ly,\ 0)$, which is neither zero nor constant, so the external work it collects there is $\int T_x(y)\,Ly\,\mathrm{d}y$ — a first moment of the axial grip traction, which the measured transverse resultant $F$ does not determine. By the tightrope rule below, the field is usable with the boundary data you actually have only if you can assert $T_x \equiv 0$ pointwise on that edge — it stays perfectly kinematically admissible either way; what fails is that its external work is not recoverable from the measured resultant. If you cannot, envelope it: $\mathbf{u}^{*} = (x(L-x)\,y,\ 0)$ has the same reach ($\varepsilon^{*}_{xx} = (L-2x)y$, $\gamma^{*}_{xy} = x(L-x)$) and the same $Q_{yy}$-blindness, but vanishes at both $x = 0$ and $x = L$, so it carries no unmeasured boundary work at all. The obvious-looking companion $\mathbf{u}^{*} = (0,\ x^2)$ does not help: $\varepsilon^{*}_{yy} = 0$ there too, leaving a shear-only equation. Reaching $Q_{yy}$ requires a virtual field whose $y$-component varies with $y$ — for instance $\mathbf{u}^{*} = (0,\ x(L-x)\,y)$, with $\varepsilon^{*}_{yy} = x(L-x)$, which also vanishes on the clamped edge and on the loaded edge and so contributes a purely homogeneous equation (fine, provided some other field still does work against the measured $F$ — see the tightrope below). The general work equation reads:

The Admissibility Tightrope

Special, Optimized, and Piecewise Virtual Fields

Three refinements (Pierron & Grédiac 2012):

• Special virtual fields: choose the fields so that $\mathbf{A}$ becomes the identity — each equation isolates exactly one parameter.
• Optimized virtual fields: among all special fields, pick the one that minimizes the variance of each parameter with respect to strain noise. This both stabilizes the identification and yields predicted standard deviations for free — error bars as a by-product of the solve.
• Piecewise (element-wise) virtual fields: constructed on a mesh, they handle complex geometries where global polynomials cannot be made admissible.

6. Nonlinear VFM and Sensitivity-Based Virtual Fields

When stress is not linear in the fitted parameters — Ogden with its exponents $\alpha_p$ among the unknowns, most plasticity laws — the PVW no longer collapses to a linear system. Note what that condition does not say: finite strain is not the criterion, and "hyperelastic" is not a synonym for "parameter-nonlinear". At fixed measured kinematics the neo-Hookean and Mooney–Rivlin energies of §7 are linear in $\mu$ and in $(C_{10}, C_{01})$ respectively — and stay so after the plane-stress pressure elimination — so their VFM identification is still one direct linear solve, at finite stretch. What breaks linearity is where the parameter sits in the law, not how large the strains are. Where it does break, VFM's core advantage survives intact: the stress is computed pointwise from the measured strain (or strain history) through the candidate law, with no boundary-value problem solved anywhere. Define the cost as the sum, over load steps and virtual fields, of squared virtual-work residuals, and minimize over $\boldsymbol{\theta}$ with gradient-based optimization — the gradients are cheap and, at fixed virtual fields, exact for the discrete evaluation, because the residual is then explicit in $\boldsymbol{\theta}$: they come from parameter sensitivities $\partial\mathbf{P}/\partial\boldsymbol{\theta}$ — obtained by differentiating the local stress-update equations with respect to $\boldsymbol{\theta}$ at fixed measured kinematics, in closed form or by automatic differentiation through the stress update — rather than from repeated forward solves:

Note the argument $\mathbf{F}_{1:t}$, the measured deformation history up to step $t$, not $\mathbf{F}_t$ alone. For an elastoplastic law the stress at a point is not a function of the current deformation gradient: it also depends on the path taken to get there, carried as internal state, so the honest arguments are $\mathbf{P}(\mathbf{F}_t, \boldsymbol{\xi}_t;\ \boldsymbol{\theta})$ with $\boldsymbol{\xi}_t$ produced by the return map from $\mathbf{F}_{1:t}$. Only for memoryless hyperelasticity does this collapse to $\mathbf{P}(\mathbf{F}_t;\ \boldsymbol{\theta})$. Note also the $t$ on $\mathbf{T}_{0,t}$: the measured traction is a function of the load step, exactly as $\mathbf{F}_t$ is — without it every step would be matched against the same constant external work, which is not the virtual-work imbalance. Body forces are dropped here relative to the full statement of the PVW in §5; restore $-\int_{V_0} \mathbf{b}_{0,t} \cdot \mathbf{u}^{*(j)}\,\mathrm{d}V_0$ inside the bracket if self-weight is not negligible — on a soft, low-modulus specimen it often is not. Two details: the minus sign, because the bracket subtracts external from internal work; and $\mathbf{b}_0 = J\mathbf{b}$, the body force per unit reference volume ($\rho_0\mathbf{g}$ for self-weight), because this residual is a total-Lagrangian statement over $V_0$ while §5's $\mathbf{b}$ is per unit current volume.

For large deformations, work in the reference configuration: the deformation gradient $\mathbf{F} = \mathbf{I} + \nabla_0\mathbf{u}$ comes from the measured displacement field, the first Piola–Kirchhoff stress $\mathbf{P}(\mathbf{F};\boldsymbol{\theta})$ from the candidate law. The internal virtual work is $\int \mathbf{P} : \nabla_0 \mathbf{u}^{*}\,\mathrm{d}V_0$ — $\mathbf{P}$ is work-conjugate to the material gradient of the virtual displacement, which is why $\nabla_0\mathbf{u}^{*}$ and not some strain measure appears there. For elastoplasticity, the stress at each DIC point comes from running the return-mapping algorithm on that point's measured strain history — still no forward boundary-value problem solved anywhere. Two qualifications, both of which the rest of this module depends on. A mesh may perfectly well be present — to build the piecewise virtual fields of §5, or simply to integrate the work terms; what VFM never does is solve the physical forward boundary-value problem at every parameter evaluation. Auxiliary global systems can still appear — the stiffness-based virtual fields of Marek, Davis & Pierron (2017) are built by inverting a global constrained system assembled from a modified stiffness matrix — but that system manufactures test functions; it is not the state solve FEMU repeats for every candidate $\boldsymbol{\theta}$, and that distinction is what actually buys the cost gap. And the return map is only piecewise smooth in $\boldsymbol{\theta}$, because its active set changes at yield, so the parameter sensitivities of the elastoplastic residual are obtained by differentiating the return-map equations themselves — in closed form, or by AD through the branch the algorithm actually executes — which is precisely what Kumar et al. (2025) supply. Both routes are exact within a fixed active set, and neither repairs the switch itself: where the active set changes, the classical derivative need not exist, and AD does not manufacture one — it differentiates the branch it happened to run down. That case needs a one-sided or generalized (Clarke) derivative, a smoothed yield transition, or an optimizer that tolerates nonsmoothness. Note which derivative this is: not the consistent algorithmic tangent $\partial\Delta\boldsymbol{\sigma}/\partial\Delta\boldsymbol{\varepsilon}$ that a Newton solver needs to iterate on displacements, but $\partial\boldsymbol{\sigma}/\partial\boldsymbol{\theta}$ at fixed measured strain. The kinematics are data here, so the strain-side tangent is exactly the object VFM never has to assemble globally.

Sensitivity-Based Virtual Fields

The Achilles heel of nonlinear VFM is virtual-field choice: hand-picked polynomials that work for one geometry fail for another, and poor fields leave the cost flat in some parameters. Sensitivity-based virtual fields (Marek, Davis & Pierron 2017) automate the choice: build each virtual field $\mathbf{u}^{*(i)}$ from the sensitivity of the stress field to parameter $\theta_i$,

— the virtual strain field is taken proportional to $\partial\boldsymbol{\sigma}/\partial\theta_i$ evaluated along the measured strain path, then rendered kinematically admissible. Each residual becomes maximally responsive to "its" parameter. Note the price this exacts on the exactness claim above: the fields now depend on $\boldsymbol{\theta}$, so $\mathrm{d}J/\mathrm{d}\boldsymbol{\theta}$ strictly acquires $\partial\mathbf{u}^{*}/\partial\boldsymbol{\theta}$ terms too. The stress-sensitivity gradient is exact only while the fields are held frozen — so freeze them through an inner minimization and rebuild them in an outer loop, or differentiate the field-construction step as well. This construction became a de facto standard for nonlinear VFM. Nikolov, Zhu & Estrada (2026) push further with variation-matching SBVFs: procedural construction via Gâteaux derivatives plus analytically-derived, model-tailored fields, demonstrated on simulated noisy volumetric (3D) data — removing manual field selection altogether.

The Incompressibility Subtlety

For exactly incompressible materials the hydrostatic pressure $p$ is not determined by the kinematics: it is a Lagrange multiplier enforcing $J = 1$, and no camera can see it. For a nearly incompressible model the situation is formally different and practically the same. There $p$ is determined — by the volumetric strain, through whatever volumetric law you chose, $p = -K(J-1)$ for the simplest one — but a bounded pressure corresponds to a volumetric strain of order $p/K$, and that number deserves checking rather than waving through. With $K/G = 2(1+\nu)/\big(3(1-2\nu)\big) \approx 500$ at $\nu = 0.499$ (§7), a pressure of order the shear modulus, $p \sim G$, gives $p/K \sim 2\times10^{-3}$ — comparable to, not safely below, the $10^{-4}$–$10^{-3}$ strain noise floor of §2; only at $p \sim 0.1G$ does it drop to $2\times10^{-4}$, still inside that band. So the volumetric strain is marginal, not invisible: determined in principle, and in practice almost never recoverable with useful precision — but that is a number to compute for your own test, not a blanket exemption. Either way, two escapes:

• Plane stress: eliminate $p$ via $\sigma_{zz} = 0$. The criterion is that the through-thickness normal stress is negligible — which is about thinness, not about how much the surface moves. A thin inflated membrane qualifies, and that is exactly why membrane inflation is a standard VFM test (§7). What breaks the trick is a genuinely thick section, or a flat two-dimensional idealization applied to a surface that is curving out of plane — there the error is in the geometry model, not in the plane-stress assumption, and it imports an uncontrolled bias into the identification either way.
• 3D: choose divergence-free virtual fields, so the pressure term vanishes identically and only deviatoric parameters are probed:
  $$\int_V p\ \nabla \cdot \mathbf{u}^{*}\, \mathrm{d}V = 0 \quad \text{when } \nabla \cdot \mathbf{u}^{*} = 0$$

7. Inverse Hyperelastic Characterization of Soft Materials

Now the specifics that matter for a tactile-sensor skin. Elastomers — the silicones cast over vision-based tactile sensors — are nearly incompressible ($\nu \approx 0.499$) with shear moduli spanning roughly tens of kPa for soft platinum-cure silicones up to the order of 1 MPa for stiffer PDMS formulations. Treat those as typical orders to be measured, not quoted: batch, cure schedule, and aging move the real value.

Candidate Strain-Energy Functions

The Cauchy stress for the incompressible neo-Hookean model shows the pressure problem explicitly:

with $p$ determined by the boundary conditions, not the kinematics — hence the previous section's divergence-free and plane-stress tricks on the VFM side. On the FEMU side, the same near-incompressibility bites the forward model — but only where the kinematics are genuinely volume-constrained. In 3D, in plane strain, and in axisymmetric models, a standard displacement-based FE solve locks volumetrically as $\nu \to 0.5$, so the forward model inside the updating loop needs a mixed displacement–pressure (or comparable) formulation. Condensed plane stress is the exception, and it is the case a thin skin or membrane usually lands in: the through-thickness stretch is free to accommodate the volume change, so the condensed constitutive matrix keeps the finite prefactor $E/(1-\nu^2)$ and the effective in-plane Lamé parameter is $\lambda^{*} = E\nu/(1-\nu^2)$, which at $\nu = 0.5$ is $2E/3$ against a shear modulus $E/3$ — a ratio of two, nowhere near the $\lambda/\mu \to \infty$ that causes locking. There is simply nothing to lock. Diagnose which regime your forward model is in before reaching for a mixed formulation. An identification is only as good as the forward solves — or the pointwise stress evaluations — inside it.

Two Landmark Applications

Arteries, with VFM. Avril, Badel & Duprey (2010) identified anisotropic hyperelastic parameters of in vitro human arterial segments from inflation and extension tests with full-field optical measurement — and did so using the virtual fields method, not a forward-model updating loop. That is exactly why it belongs in this module: their formulation handles the heterogeneous strain and stress distributions of an arterial segment rather than requiring a homogeneous gauge section. It is the template for membrane-type identification: pressurize, measure the surface kinematics full-field, enforce virtual work with the measured strains.

Elastomeric foam, with VFM-GA. Yan, Tao, Franck & Henann (2025) calibrate isotropic hyperelastic models — including compressible forms with coupled dependencies on the deformation invariants, applied to open-cell elastomeric foam — via a VFM objective searched with a genetic algorithm, using DIC displacement fields (both homogeneous and inhomogeneous) synchronized with force data. The GA is deployed, in the authors' own framing, for its ability to navigate the intricate landscape of material-parameter space; separately, the framework places no restriction on the target isotropic model, accommodating coupled invariant dependencies and specialized functional forms. It also includes a material-stability screen, eliminating parameter sets that would imply non-physical behavior.

Practical Guidance for the Tactile-Skin Use Case

• Test the actual cast material. Batch, cure, and aging matter more than datasheets.
• Prefer a heterogeneous test spanning the strain states seen in service. Indentation is closest to sensor loading; membrane inflation gives clean biaxial states.
• Mind viscoelasticity. Quasi-static hyperelastic fits are rate-dependent for silicones: test at service-representative rates, or model the rate dependence.
• Validate by predicting an independent experiment, never by the fitting residual (§11).

8. Identifiability, Noise, and Test Design

Identification quality is decided before the experiment runs. This section is the module's control center: what makes parameters recoverable, how noise reaches them, and how to design the test so they arrive with small error bars.

Local Identifiability: the Fisher Matrix

Near the optimum, the Gauss–Newton (Fisher) matrix $\mathbf{S}^{\mathsf{T}}\mathbf{W}\mathbf{S}$ governs everything. Its condition number flags near-collinear parameter directions; its inverse, times the noise variance, approximates the parameter covariance:

The off-diagonal structure is read through the correlation matrix:

A practical footnote on parameter scaling: $E$ in MPa and $\nu$ dimensionless differ by orders of magnitude in raw units, and a condition number computed on the raw parameterization mostly measures that scaling, not the physics. Read conditioning on relative (scaled) parameters — the Identifiability Lab below does.

The Classic Pathologies

Mooney–Rivlin: $C_{10}$ and $C_{01}$ are nearly unidentifiable from uniaxial-dominated data, because $I_1$ and $I_2$ vary almost proportionally along uniaxial paths — the two terms' stress contributions are nearly collinear, and $|r_{12}| \to 1$. Separation is enormously improved by strain states spanning the $(I_1, I_2)$ plane: simultaneous uniaxial, shear (planar tension), and biaxial regions in one heterogeneous field. Note the word: improved, not enabled. The two coefficients are not structurally unidentifiable from uniaxial data — for incompressible uniaxial tension the nominal stress is $P = 2(\lambda - \lambda^{-2})(C_{10} + C_{01}/\lambda)$, so they enter with different $\lambda$-dependence and a wide enough stretch range separates them in principle. It is the conditioning that is hopeless over any realistic range, which is a quantitative statement about error bars, not a rank deficiency. Ogden: the $(\mu_p, \alpha_p)$ pairs trade off strongly at small stretch. Hyperelastic parameters are, as a rule, correlated — report the correlation matrix, not just the diagonal.

Structural non-identifiability is worse than conditioning: with displacement-controlled BCs and no force channel, a uniform scaling of all moduli leaves the displacement field unchanged. The cost is exactly flat along the stiffness-scale direction — no amount of data or noise reduction helps. Any one independent stiffness datum removes it — a load measurement is the usual one, a known traction or a region of known modulus works too (§3, §5, and the Lab below, all telling the same story).

Noise Propagation, and Two Caveats on the Covariance Formula

For linear VFM, $\hat{\boldsymbol{\theta}} = \mathbf{A}^{-1}\mathbf{b}$, so strain noise propagates through $\mathbf{A}^{-1}$ — and optimized virtual fields (§5) explicitly minimize each parameter's noise sensitivity, returning predicted standard deviations as a by-product (Pierron & Grédiac 2012). For any method, two caveats before you trust $\sigma^2(\mathbf{S}^{\mathsf{T}}\mathbf{S})^{-1}$:

1. DIC noise is spatially correlated after subset matching and strain smoothing. Treating every DIC point as an independent sample grossly understates the uncertainty — absurdly small confidence intervals from $10^5$ "independent" points are a red flag, not a triumph. Characterize the actual noise from static image pairs.
2. Covariance quantifies noise-driven scatter only. Model error — wrong constitutive form, violated plane-stress assumption, idealized BCs — produces bias that is invisible to it. The residual field's spatial structure is the diagnostic (§4), not the error bar — read against a characterized noise model, since correlated DIC noise and interpolation bias leave spatial structure of their own.

Test Design: the Material Testing 2.0 Core

Design the specimen and loading — holes, notches, off-axis fibers, multiaxial rigs — to maximize sensitivity coverage (Pierron & Grédiac 2021). Rank candidate designs by predicted parameter standard deviations, or by D-optimality:

Best practice simulates the entire chain: FE-generate synthetic deformed speckle images, run them through the actual DIC engine, then through the identification — so the error budget includes DIC systematics (interpolation bias, VSG smoothing), not just idealized white noise.

9. FEMU vs VFM: An Honest Comparison

Both methods answer the same question from the same data. The differences are structural, and each direction of superiority has a matching direction of weakness.

What the Evidence Actually Supports

Avril et al. (2008) is the field's common reference point: it formulates FEMU, CEGM, VFM, EGM and RGM side by side and exercises the applicable ones on full-field experimental data from four tests — tensile, Brazilian, shear-flexural, and a biaxial test with a non-uniform damage field. Use it for how the methods relate and where each is applicable, not as a single accuracy ranking. Kumar, Seidl, Granzow, Yang & Fuhg (2025) revisited the comparison for finite-strain elastoplasticity with numerically-exact automatic-differentiation sensitivities for both methods — arguing that prior comparisons were confounded by gradient approximation error — and assessed robustness under practical data corruptions on an even footing.

10. Beyond Parameters: Constitutive Model Discovery

The 2020s reframed the question. Instead of "fit $\boldsymbol{\theta}$ inside a model someone chose", ask: discover the model itself from the data.

EUCLID: Sparse Discovery on a Feature Library

EUCLID (Flaschel, Kumar & De Lorenzis 2021) writes the strain-energy density as a large linear combination of interpretable candidate features — invariant powers, logarithmic and volumetric terms:

It enforces weak-form equilibrium evaluated with many test functions plus global force balance — structurally, a VFM residual with many virtual fields — and applies sparsity-promoting regression so only a few physically meaningful terms survive:

Both equilibrium blocks are written out, because the prose above claims both and a residual labelled only "weak form" would not encode the second. The first block is the weak linear-momentum balance at the free degrees of freedom; the second is the reaction-force balance at the constrained ones, weighted by its own hyperparameter $\lambda_r$ (needed because the boundary contributes far fewer equations than the interior); the third promotes sparsity. One notation warning: the original paper calls the sparsity exponent $p$, renamed to $q$ here because $p$ is the parameter count everywhere else in this module.

Inputs: the displacement field plus global reaction forces — both realistically available from DIC-instrumented testing. No stress data anywhere: stress is never measured, only inferred, which is what "unsupervised" means here. Output: a closed-form, interpretable hyperelastic law with fitted coefficients. The original paper demonstrates recovery of five hyperelastic models of varying complexity on numerically generated data, including graceful behaviour when the true function is absent from the candidate library.

Neural but Interpretable: CANN-EUCLID

CANN-EUCLID (Alheit, Kumar & Peirlinck 2026) trains a Constitutive Artificial Neural Network — architecture constrained so every node maps to an interpretable term — through the EUCLID residual on full-field displacement and force data, with sparsity-promoting regularization selecting a compact set of active terms; benchmarks cover isotropic and anisotropic hyperelasticity with prescribed ground-truth laws. The paper documents a sobering failure mode in its own words: exponential strain-stiffening terms "can be recovered accurately when sufficiently probed, but can produce large extrapolation errors when the stiffening regime lies outside the sampled domain" — a general warning for all data-driven constitutive laws, and the reason §7 insisted on independent-experiment validation.

The Unifying Frame: Reduced vs All-at-Once

Römer, Hartmann, Tröger, Anton, Wessels, Flaschel & De Lorenzis (2025) organize the whole field along one axis. At one end, all-at-once formulations carry state fields and parameters as joint unknowns, combining the model residual and the data misfit in a single objective — inverse (as opposed to parametric) PINNs sit here. The reduced formulation is the limit in which the discrete model equation is imposed as a strict constraint, so only $\boldsymbol{\theta}$ remains unknown — conventional FEMU sits here. The placement worth noting is VFM's: the review recovers it not alongside FEMU but at the other end of the same spectrum, as the limit in which the data-related objective is enforced strictly — and it also formulates an all-at-once VFM variant, which is the clearest sign that these are labels for formulations, not for method families. (Parametric PINNs, by the same logic, behave like the reduced approach.) They fold in statistical uncertainty quantification, propose a hybrid frequentist–Bayesian identification procedure for complex material models, and validate several methods on benchmarks with synthetic and real data. (That the classical least-squares estimators correspond to maximum-a-posteriori estimates under Gaussian noise and flat priors is the standard statistical reading — offered here as this module's synthesis, not as a headline claim of the review.)

Whichever end of the spectrum you work at, the closing discipline is the same: the provenance of the identification — which data, which residual, which regularization — determines what the "identified material" actually means. A modulus is never just a number; it is a number plus the experiment and the residual that produced it.

11. A Practical Identification Workflow for a Silicone Skin

The module condensed into a recipe you can run in the lab. Each step references the section that justifies it.

1. Design (§8). Pick a heterogeneous test spanning the service strain states: open-hole tension or planar-tension-with-notch for in-plane characterization; membrane inflation or instrumented indentation with stereo-DIC when out-of-plane response matters. Before cutting silicone, run the synthetic pipeline: FE-simulate the candidate test with plausible parameters, deform synthetic speckle images, run them through the actual DIC software, run the identification. This reveals identifiability problems and realistic error bars while they are still free to fix.
2. Speckle (§2). Airbrushed or spray pattern, feature size 3–5 px at the chosen magnification. On soft thin samples use flexible ink sparingly — a continuous paint film can locally stiffen the very surface being measured. Acquire static image pairs and measure the actual displacement/strain noise floor.
3. Test (§7). Synchronized load cell + camera(s); quasi-static rate matched to service conditions — silicones are mildly viscoelastic, and a rate mismatch biases hyperelastic fits. Record the true unloaded reference frame. For stereo, calibrate immediately before testing.
4. DIC (§2). Choose subset, step, and VSG from the noise-floor study, not from software defaults. Export displacements, strains, and coordinates in the specimen frame.
5. Identify (§5–§8). Start with the simplest credible model (neo-Hookean, one parameter) and add parameters only when the data can support them — watch the correlation matrix $r_{ij}$. Run VFM first where the specimen justifies a reduced model — thin sheet, membrane, plane stress — because its work integral needs $\mathbf{F}$ throughout the analyzed volume, and surface DIC on a thick indented block does not supply that; for genuinely 3D tests, go to FEMU (or measure volumetrically with DVC). Where VFM applies it costs seconds, with a deviatoric or plane-stress treatment for incompressibility; then cross-check with FEMU. Report $\hat{\boldsymbol{\theta}}$ with covariance and with residual-field maps.
6. Validate (§7, §9). Predict an experiment not used in fitting — e.g. identify from tension, predict indentation force–depth. Agreement there, not the fitting residual, is the evidence that the material model transfers to the sensor simulation.
7. Report. Images, calibration, DIC settings, noise floor, virtual fields or FEMU configuration — so the identification is reproducible.

Flashcards

Seventeen cards: the two method families, the PVW machinery, the identifiability results, and the numbers worth knowing cold.

References

• Grédiac, Michel (1989). Principe des travaux virtuels et identification. Comptes Rendus de l’Académie des Sciences, Série II, 309, 1–5. The paper that introduces the virtual fields method.
• Kavanagh, Kenneth T.; Clough, Ray W. (1971). Finite element applications in the characterization of elastic solids. International Journal of Solids and Structures, 7(1):11–23. doi:10.1016/0020-7683(71)90015-1
• Pierron, F.; Grédiac, M. (2021). Towards Material Testing 2.0. A review of test design for identification of constitutive parameters from full-field measurements. Strain, 57(1):e12370 (online 2020). doi:10.1111/str.12370
• Pierron, Fabrice (2023). Material Testing 2.0: A brief review. Strain, 59(3):e12434. doi:10.1111/str.12434
• Sutton, Michael A.; Orteu, Jean-José; Schreier, Hubert (2009). Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts, Theory and Applications. Springer US (book). doi:10.1007/978-0-387-78747-3
• Hild, François; Roux, Stéphane (2006). Digital Image Correlation: from Displacement Measurement to Identification of Elastic Properties – a Review. Strain, 42(2):69–80. doi:10.1111/j.1475-1305.2006.00258.x
• Avril, Stéphane; Bonnet, Marc; Bretelle, Anne-Sophie; Grédiac, Michel; Hild, François; Ienny, Patrick; Latourte, Félix; Lemosse, Didier; Pagano, Stéphane; Pagnacco, Emmanuel; Pierron, Fabrice (2008). Overview of Identification Methods of Mechanical Parameters Based on Full-field Measurements. Experimental Mechanics, 48(4):381–402. doi:10.1007/s11340-008-9148-y
• Grédiac, Michel; Pierron, Fabrice; Avril, Stéphane; Toussaint, Evelyne (2006). The Virtual Fields Method for Extracting Constitutive Parameters From Full-Field Measurements: a Review. Strain, 42(4):233–253. doi:10.1111/j.1475-1305.2006.00283.x
• Pierron, Fabrice; Grédiac, Michel (2012). The Virtual Fields Method: Extracting Constitutive Mechanical Parameters from Full-field Deformation Measurements. Springer New York (book). doi:10.1007/978-1-4614-1824-5
• Marek, Aleksander; Davis, Frances M.; Pierron, Fabrice (2017). Sensitivity-based virtual fields for the non-linear virtual fields method. Computational Mechanics, 60(3):409–431. doi:10.1007/s00466-017-1411-6
• Nikolov, Denislav P.; Zhu, Zhiren; Estrada, Jonathan B. (2026). Variation-matching sensitivity-based virtual fields for hyperelastic material model calibration. Strain, 62(1):e70031 (2026); arXiv:2509.03363 (2025). doi:10.1111/str.70031
• Avril, Stéphane; Badel, Pierre; Duprey, Ambroise (2010). Anisotropic and hyperelastic identification of in vitro human arteries from full-field optical measurements. Journal of Biomechanics, 43(15):2978–2985. doi:10.1016/j.jbiomech.2010.07.004
• Yan, Zicheng; Tao, Jialiang; Franck, Christian; Henann, David L. (2025). A Virtual Fields Method-Genetic Algorithm (VFM-GA) calibration framework for isotropic hyperelastic constitutive models with application to an elastomeric foam material. arXiv:2510.07683. doi:10.48550/arXiv.2510.07683
• Kumar, Sanjeev; Seidl, D. Thomas; Granzow, Brian N.; Yang, Jin; Fuhg, Jan N. (2025). A comparative study of calibration techniques for finite strain elastoplasticity: Numerically-exact sensitivities for FEMU and VFM. arXiv:2503.19782. doi:10.48550/arXiv.2503.19782
• Flaschel, Moritz; Kumar, Siddhant; De Lorenzis, Laura (2021). Unsupervised discovery of interpretable hyperelastic constitutive laws. Computer Methods in Applied Mechanics and Engineering, 381:113852. doi:10.1016/j.cma.2021.113852
• Alheit, Benjamin; Kumar, Siddhant; Peirlinck, Mathias (2026). CANN-EUCLID: unsupervised constitutive artificial neural network model discovery from full-field data. arXiv:2606.14565. doi:10.48550/arXiv.2606.14565
• Römer, Ulrich; Hartmann, Stefan; Tröger, Jendrik-Alexander; Anton, David; Wessels, Henning; Flaschel, Moritz; De Lorenzis, Laura (2025). Reduced and All-at-Once Approaches for Model Calibration and Discovery in Computational Solid Mechanics. Applied Mechanics Reviews, 77(4) (2025); arXiv:2404.16980 (2024).