9. Reference Configuration & Prestress

The unknown reference configuration, inverse motion and Shield duality, the Eshelby stress, Govindjee–Mihalic inverse elastostatic finite elements, Sellier fixed-point unloading with relaxation and continuation, MULF prestressing, load-free vs stress-free, joint nonuniqueness of geometry and material

Contents

1. The Configuration You Never Measured 2. Finite Kinematics, Just In Time 3. Four Configurations, Not Two 4. Shield Duality and the Eshelby Stress 5. Inverse Elastostatic Finite Elements 6. Fixed-Point Unloading: Relaxation, Acceleration, Continuation 7. Prestressing Instead of Un-Deforming 8. Uniqueness, Conditioning, and Static Determinacy 9. Practice: Choosing a Method, and the Pitfalls That Bite 10. Applications and the 2023–2026 Frontier Interactive: Unloaded-Shape Laboratory Interactive: Joint Identifiability Explorer Interactive: Unloading a Patient-Specific Aorta Flashcards

1. The Configuration You Never Measured

The constitutive kinematics of finite strain are written from a reference configuration $\mathcal{B}_0$: strain is measured from it, and the strain-energy function is a function of the gradient taken with respect to it. Equilibrium itself can be posed either way — materially on $\mathcal{B}_0$ or spatially on $\mathcal{B}$ — and the element shape functions are defined on a parent element and mapped onto whichever domain you integrate over. That flexibility is what §4 and §5 exploit. But the constitutive dependence on $\mathcal{B}_0$ is not negotiable, and in most real problems you never observe $\mathcal{B}_0$.

A CT scan of an aorta is taken at roughly 80 mmHg diastolic pressure, so the segmented mesh is already stretched. A turbine blade is designed for its hot, spun-up shape, but the machinist needs the cold shape. A silicone tactile skin is designed for the shape it should hold on the sensor body under gravity and membrane pretension, but the mold must be cut for the un-sagged shape. In all three cases the deformed geometry is data and the reference geometry is the unknown.

Core Problem

Given the deformed configuration $\mathcal{B}$ of a hyperelastic body — imaged, scanned, or specified by design — together with the loads acting on it and its constitutive law $W$, find the reference configuration $\mathcal{B}_0$ from which it came, and the stress it currently carries. In the guide's observation model $y = \mathcal{H}(u, p, f, X_0) + b + \eta$, the material parameters $p$ and the loads $f$ are assumed known, the deformed geometry is the measurement, and the unknown is $X_0$ itself.

	In this module
Measured	The deformed surface or volume mesh $\mathcal{B}$ — segmented from CT/MRI, scanned, or given as a design target. Carries voxel-scale geometric noise.
Unknown	The reference configuration $\mathcal{B}_0$, equivalently the material coordinates $\mathbf{X}(\mathbf{x})$ of every point of $\mathcal{B}$ — or, in the prestressing route, a prestretch field $\mathbf{F}_{\text{pre}}(\mathbf{x})$ on $\mathcal{B}$.
Assumed known	The constitutive law $W$ and its parameters, the loads (pressure level, body force), the boundary conditions, and — for anisotropic tissue — the fibre architecture.
Unobservable	Residual stress in the unloaded state (§3); the split between stiffness and reference geometry from a single loaded state (§8); any reference configuration reachable through a different branch of a multi-valued forward map (§8).

The Cost of Pretending the Image Is Stress-Free

This is not a refinement. Ramella et al. (2024) ran TAVI and TEVAR device-deployment simulations with and without arterial pre-stress recovered by an inverse elastostatic solve, and compared. Including pre-stress raises the predicted contact pressure by 48% and 55%, the aortic wall stresses by 162% and 157%, and the aortic wall strains by 77% and 21% relative to treating the CT geometry as the reference, for the TAVI and TEVAR models respectively. (Read the direction carefully: those are increases measured against the CT-as-reference result. Expressed the other way — how far the CT-referenced model falls short of the pre-stressed one — the same numbers are 32% / 36%, 62% / 61% and 44% / 17%. A shortfall can never exceed 100%; an increase can.) Treating the imaged geometry as stress-free is a first-order modelling error.

Three Problems Sharing One Forward Operator

The whole module hangs on keeping three distinct inverse problems apart. They share the forward operator and differ only in which slot of it is unknown.

\begin{aligned}\text{forward:}\quad & (\mathcal{B}_0,\; W,\; \text{loads}) \;\longmapsto\; \boldsymbol{\varphi},\ \mathcal{B}\\[2pt] \text{inverse motion:}\quad & (\mathcal{B},\; W,\; \text{loads}) \;\longmapsto\; \boldsymbol{\varphi}^{-1},\ \mathcal{B}_0\\[2pt] \text{prestressing:}\quad & (\mathcal{B},\; W,\; \text{loads}) \;\longmapsto\; \mathbf{F}_{\text{pre}}\ \text{on}\ \mathcal{B}\\[2pt] \text{elastography:}\quad & (\mathcal{B}_0,\; \mathbf{u}^{\text{meas}},\; \text{loads}) \;\longmapsto\; W\ \text{or}\ \mu(\mathbf{x})\end{aligned}

Sections 4–6 and 8–10 treat the inverse motion problem; Section 7 treats prestressing. The fourth line — the unknown is a material field, not a configuration — is a different inverse family and has its own module.

Pointer — elastography lives in Module 7

Distributed-parameter material identification (recover $\mu(\mathbf{x})$ from a measured displacement field) is the canonical subject of Module 7, and this module does not re-derive it. Three results from there you will need when Section 8 discusses nonuniqueness: the adjoint method delivers the full gradient of the misfit for one forward plus one adjoint solve, at a cost independent of the number of parameters (Oberai, Gokhale & Feijóo, 2003) — the mechanics of which are Module 10's subject; the direct reconstruction has no unique solution unless stiffness is known a priori on a sufficient portion of the boundary, and identical strain fields follow from different modulus distributions under displacement boundary conditions (Barbone & Bamber, 2002); and with $N$ independent measured displacement fields the residual freedom is two arbitrary functions at $N=1$, at most four arbitrary constants at $N=2$, and one constant at $N=4$ (Barbone & Gokhale, 2004). Doyley (2012) surveys the model-based family. The structural lesson transfers directly: a raw strain image is not a modulus image, and an unloaded geometry recovered from one loaded state is not an identified material model.

The organising picture. One forward operator, four inverse families distinguished only by which slot is unknown. The reference configuration carries a question mark because in the problems of this module it is never observed — only the current configuration is.

Warning — two unrelated things called “inverse FEM”

In the structural-health-monitoring and aerospace literature, iFEM means the Tessler–Spangler inverse Finite Element Method for shape sensing from discrete surface strain gauges (Modules 4–6): geometry and material are known, the unknown is the displacement field. This module's inverse FEM is inverse elastostatics: the displacement data are complete, and the unknown is the reference geometry. Search for “inverse finite element method” and you will get plate, shell, hull and riser shape-sensing papers, not Govindjee–Mihalic. The two families are not variants of each other; they do not even share an unknown. Say this out loud once, or your own literature search will mislead you for a week.

2. Finite Kinematics, Just In Time

Modules 1–6 could live comfortably in small strain. This one cannot: the entire problem is that a 15% stretch makes the reference and current configurations genuinely different domains. Here is the minimum kinematic apparatus, stated in the two dual descriptions you will need.

The motion $\boldsymbol{\varphi}: \mathcal{B}_0 \to \mathcal{B}$ maps material points $\mathbf{X}$ to spatial points $\mathbf{x} = \boldsymbol{\varphi}(\mathbf{X})$. Its gradient is the deformation gradient $\mathbf{F}$, and $J = \det\mathbf{F} > 0$ is the volume ratio. If $\boldsymbol{\varphi}$ is a diffeomorphism it has an inverse motion $\boldsymbol{\Phi} = \boldsymbol{\varphi}^{-1}: \mathcal{B}\to\mathcal{B}_0$, $\mathbf{X} = \boldsymbol{\Phi}(\mathbf{x})$, whose gradient is the inverse deformation gradient $\mathbf{f}$.

\mathbf{F} = \frac{\partial \boldsymbol{\varphi}}{\partial \mathbf{X}},\qquad \mathbf{f} = \frac{\partial \boldsymbol{\Phi}}{\partial \mathbf{x}} = \mathbf{F}^{-1},\qquad \det\mathbf{f} = J^{-1}

The inverse motion problem

The spatial domain $\mathcal{B}$ is given. Find the field $\mathbf{X}(\mathbf{x})$ defined on it. Nothing about the physics changes — only which of the two fields, $\mathbf{x}(\mathbf{X})$ or $\mathbf{X}(\mathbf{x})$, is data.

Stresses: the Cauchy stress $\boldsymbol{\sigma}$ acts on the current area, the first Piola–Kirchhoff stress $\mathbf{P}$ on the reference area, and for a hyperelastic material $\mathbf{P}$ is the derivative of the strain energy with respect to $\mathbf{F}$.

\mathbf{P} = \frac{\partial W}{\partial \mathbf{F}} = J\boldsymbol{\sigma}\mathbf{F}^{-T},\qquad \boldsymbol{\sigma} = J^{-1}\mathbf{P}\mathbf{F}^{T} = J^{-1}\frac{\partial W}{\partial \mathbf{F}}\mathbf{F}^{T}

Equilibrium can be written referentially or spatially. The two statements are equivalent when both configurations are known; when $\mathcal{B}_0$ is unknown only the spatial one can be discretised directly, because only it integrates over a domain you actually have. The referential statement is not thereby wrong or unusable — it has to be pulled onto $\mathcal{B}$ by a change of variables first, and once it has been, it is the spatial form.

\operatorname{Div}\mathbf{P} + \rho_0\mathbf{B} = \mathbf{0}\ \ \text{on }\mathcal{B}_0, \qquad \operatorname{div}\boldsymbol{\sigma} + \rho\,\mathbf{b} = \mathbf{0}\ \ \text{on }\mathcal{B}

Key insight — write equilibrium on the domain you know

The spatial form $\operatorname{div}\boldsymbol{\sigma} + \rho\mathbf{b} = \mathbf{0}$ is the one that survives when $\mathcal{B}_0$ is unknown, because every integral in its weak form is taken over $\mathcal{B}$ — the measured domain. That single observation is the whole of Section 5's finite element method. It is also why the load side gets easier, not harder: the surface on which the pressure acts is known before the solve starts.

Nanson's Formula, and Why Tractions Are Not Portable

Areas transform between configurations by Nanson's formula, which relates an oriented reference area element to its spatial image:

\mathrm{d}a\,\mathbf{n} = J\,\mathbf{F}^{-T}\,\mathrm{d}A\,\mathbf{N}

This is bookkeeping until you try to move a load between configurations, at which point it becomes the difference between two different problems. A pressure $\bar{\mathbf{t}} = -p\,\mathbf{n}$ specified on the current surface is not the same boundary condition as a dead traction of magnitude $p$ specified on the reference surface: as the body deforms the first changes both its direction and the area it acts on, and the second does neither. Pulling a traction back or pushing it forward without applying Nanson's transformation silently replaces your boundary-value problem with a neighbouring one whose solution is a different reference configuration.

Trap — careless traction pull-back changes the inverse problem

In the forward problem a sloppy traction transformation costs you accuracy. In the inverse problem it costs you the answer: the recovered $\mathcal{B}_0$ is defined as the configuration that maps to $\mathcal{B}$ under the given loads, so mis-specifying the loads returns a confidently converged, wrong geometry. Decide explicitly, and write down, whether each load is a pressure on the current surface, a dead traction on the reference surface, or a total force; state it in the method section; and check that whatever route you use — direct inverse solve or fixed-point loop — applies the same one. Do not expect Shield's theorem (§4) to sort this out for you: it is a body-force-free statement about which deformations a dual material supports, and it says nothing about how your prescribed loads transform.

The Trick That Does Not Work

The most common intuition about this whole subject is: apply $-p$ to the CT geometry and call the result the unloaded shape. It is exact in linear elasticity, where superposition holds and geometry is not part of the operator. In finite strain it is wrong for two independent reasons.

The tangent is configuration-dependent. The stiffness that resists deflation at the imaged state is not the stiffness that resists inflation at the true reference state. Deflate by $-p$, then re-inflate the result by $+p$, and you do not get the imaged geometry back.
Pressure is a follower load. Its direction is normal to the current surface and the area it acts on changes with deformation — exactly the bookkeeping Nanson's formula does and the naive backward-load trick ignores.

Resist the urge to attach a headline error percentage to this. The size of the deflate-then-reinflate error depends on the material model, the stretch level and the wall geometry, and no single number is defensible across cases. The Unloaded-Shape Laboratory below computes it for an incompressible neo-Hookean thick-walled tube, so you can watch the error grow as you soften the wall or raise the pressure, instead of taking a figure on trust.

3. Four Configurations, Not Two

The literature routinely calls the output of an unloading algorithm “the stress-free reference configuration”. For a rubber part that is fine. For an artery or a ventricle it is false, and the falsehood matters. There are four names in circulation and you should be able to say which one a paper — or your own code — means. Count them carefully, though: they label three distinct configurations, because “residually stressed” is not a fourth body but the load-free body described through its stress field rather than its shape.

Configuration	External load	Internal stress	How you get it
Loaded / imaged	$p \ne 0$	$\boldsymbol{\sigma} \ne 0$	Measured. The data.
Load-free (“unloaded”)	$p = 0$	$\boldsymbol{\sigma} \ne 0$ only if the model carries a residual-stress variable	Inverse elastostatics (§5) or fixed-point unloading (§6).
↳ Residually stressed a description, not a fourth body	$p = 0$	$\boldsymbol{\sigma} \ne 0$, self-equilibrated	The same object as the row above, viewed through its stress field rather than its shape.
Stress-free	—	$\boldsymbol{\sigma} = 0$ pointwise	Cut the body until it relaxes. Generally not a single connected configuration.

The classical demonstration is the arterial ring: cut a load-free arterial ring radially and it springs open into a sector. What the cut releases is the residual stress that was locked into the closed, load-free ring; what you observe is the resulting change in geometry. The opening angle fixes the stress-free sector, and converting it into a residual strain and then a residual stress takes the open and closed radii, the axial stretch and a constitutive law — a model, not a reading.

Nothing in the inverse motion problem recovers any of this, and the reason is worth stating precisely because it is easy to get backwards. An inverse elastostatic solve returns the configuration that the constitutive law you supplied takes as its reference. With an ordinary $W(\mathbf{F})$ and no residual-stress variable, that reference is stress-free by construction and is simultaneously the configuration in equilibrium with zero external load — the method cannot return anything else, and calling its output “residually stressed” is then simply wrong. Give the model an explicit opening-angle, prestrain, growth or initial-stress field — exactly what the $\Theta_0$ slider supplies in the laboratory below — and the two part company: what the solve returns is then the stress-free reference, and the closed load-free body is a separate zero-pressure solve on top of it, which is where the residual stress lives. So the distinction load-free $\ne$ stress-free is real, but it is a statement about your model, not something the solver discovers. Watch it in the widget: at $\Theta_0 = 0$ the recovered reference and the load-free shape are the same numbers to bisection tolerance; set $\Theta_0 = 60^\circ$ and they separate.

Key insight — load-free is not stress-free

Inverse elastostatics answers “what shape does the constitutive law I supplied take as its reference?” With a plain $W(\mathbf{F})$ that shape is simultaneously the one in equilibrium with zero external load and a stress-free one, so the question it never answers is “how much stress does the real tissue carry at zero load?” Residual stress is an additional, separate unknown, parameterised in the classical treatment by an opening angle. If you repeat the phrase “stress-free reference configuration” from a paper's title, flag in your own text that the method cannot recover a residual-stress field it was never given a variable for — that needs an opening-angle, prestrain, growth or initial-stress model on top. The Unloaded-Shape Laboratory lets you switch a nonzero opening angle on and watch the load-free tube acquire a compressive-inner / tensile-outer residual hoop stress while its shape stays perfectly ordinary.

The configuration ladder. Unloading algorithms take you from the left box to the middle box. The third box is a separate object only when the model carries residual stress: with a plain $W(\mathbf{F})$ and no residual-stress variable the middle and right boxes are the same body, and the ladder has two rungs. When residual stress is present in the real organ but absent from the model, getting to the right box requires cutting the specimen — and the stress-free state of a residually stressed organ is generally not one connected configuration at all. When residual stress is instead modelled, by supplying an opening angle, the arrow relabels: the solve's unknown becomes the right box (the stress-free sector), and the middle box is recovered from it by one further zero-load solve.

4. Shield Duality and the Eshelby Stress

The structural insight behind everything that follows is a duality identity of Shield (1967). It says that the inverse motion is not a new kind of object requiring new kinematics — its stored energy is that of an ordinary hyperelastic material with a dual energy $W^{*}$. What Shield's identity does not hand you is the loaded boundary-value problem: as we are about to see, the loads must be imposed in the spatial configuration, and taking that step is Govindjee & Mihalic's contribution rather than a corollary of the energy identity.

Shield's inverse deformation theorem (1967)

If a deformation $\boldsymbol{\varphi}$ is supported without body force by a homogeneous hyperelastic material with strain energy $W(\mathbf{F})$, then the inverse motion $\boldsymbol{\Phi} = \boldsymbol{\varphi}^{-1}$ is supported, again without body force, by the dual material with strain energy $W^{*}(\mathbf{f}) = \det(\mathbf{f})\,W(\mathbf{f}^{-1})$, $\mathbf{f} = \mathbf{F}^{-1}$. Read the hypotheses: homogeneity and no body force are needed, and the theorem is a statement about which deformations are supported. The surface data that accompany the dual problem are the tractions conjugate to $\mathbf{f}$ — the configurational (Eshelby) tractions derived below — not the physical tractions you measured. That is exactly why the identity does not hand you a loaded boundary-value problem you can pose and solve.

The dual energy is not an exotic construction. It is the ordinary strain energy measured per unit current volume instead of per unit reference volume, re-expressed in the inverse gradient. Change variables in the total stored energy:

\int_{\mathcal{B}_0} W\,\mathrm{d}V = \int_{\mathcal{B}} J^{-1}W\,\mathrm{d}v = \int_{\mathcal{B}} \det(\mathbf{f})\,W\!\left(\mathbf{f}^{-1}\right)\mathrm{d}v

So the stored energy has a clean representation on the known domain $\mathcal{B}$. The tempting next move — append the usual external-work terms, written against the spatial positions, and call the result a potential — is a trap worth walking into deliberately, because it fails for an instructive reason:

\Pi[\boldsymbol{\Phi}] \;\overset{?}{=}\; \int_{\mathcal{B}} \det(\mathbf{f})\,W\!\left(\mathbf{f}^{-1}\right)\mathrm{d}v \;-\; \underbrace{\int_{\mathcal{B}}\rho\,\mathbf{b}\cdot\mathbf{x}\,\mathrm{d}v \;-\; \int_{\partial\mathcal{B}_t}\bar{\mathbf{t}}\cdot\mathbf{x}\,\mathrm{d}a}_{\text{contains no }\boldsymbol{\Phi}}

Key insight — why that functional cannot be right

On $\mathcal{B}$ the spatial coordinate $\mathbf{x}$ is the independent variable; the unknown field is $\boldsymbol{\Phi}:\mathbf{x}\mapsto\mathbf{X}$. Both external terms above are therefore constants under $\delta\boldsymbol{\Phi}$ and contribute nothing to the Euler–Lagrange equations: the "loaded" potential is load-independent, and stationarity would return the same reference configuration whether the body carries $200$ mmHg or nothing at all. Writing the loads against $\mathbf{X}$ instead does not rescue it either — that is a different physical problem (dead loading in the reference configuration), not the one you measured.

The formulation that does work imposes equilibrium directly in the spatial weak form and treats $\boldsymbol{\Phi}$ as the unknown. For every admissible virtual displacement $\delta\mathbf{u}$ on the known deformed body:

\mathcal{R}(\boldsymbol{\Phi};\delta\mathbf{u}) \;=\; \int_{\mathcal{B}} \boldsymbol{\sigma}(\boldsymbol{\Phi})\!:\!\nabla_{\!\mathbf{x}}\,\delta\mathbf{u}\;\mathrm{d}v \;-\; \int_{\mathcal{B}} \rho\,\mathbf{b}\cdot\delta\mathbf{u}\,\mathrm{d}v \;-\; \int_{\partial\mathcal{B}_t} \bar{\mathbf{t}}\cdot\delta\mathbf{u}\,\mathrm{d}a \;=\; 0

Here the loads are prescribed on the deformed configuration — where they physically act and where you actually measured them — and the Cauchy stress is evaluated from the inverse deformation gradient, $\mathbf{F} = \mathbf{f}^{-1}$ with $\mathbf{f} = \nabla_{\!\mathbf{x}}\boldsymbol{\Phi}$, through $\boldsymbol{\sigma} = J^{-1}\mathbf{P}\mathbf{F}^{\mathsf T}$. Now $\boldsymbol{\Phi}$ appears in the internal term and the loads survive into the residual. This is the route Govindjee & Mihalic take, and they are explicit that the classical Eshelby / conservation-law machinery is inadequate for exactly this data — prescribed deformed geometry with prescribed Cauchy tractions.

The Stress Conjugate to $\mathbf{f}$

Differentiate $W^{*}$ and a familiar object falls out. Every transpose matters here, so the derivation is worth doing in full rather than quoting.

Full derivation: $\partial W^{*}/\partial\mathbf{f}$ and the Eshelby tensor

Two variational identities do all the work:

\delta(\det\mathbf{f}) = \det(\mathbf{f})\,\mathbf{f}^{-T}\!:\!\delta\mathbf{f} = J^{-1}\mathbf{F}^{T}\!:\!\delta\mathbf{f}, \qquad \delta\!\left(\mathbf{f}^{-1}\right) = -\mathbf{F}\,\delta\mathbf{f}\,\mathbf{F}

The first is Jacobi's formula with $\mathbf{f}^{-T} = \mathbf{F}^{T}$ and $\det\mathbf{f} = J^{-1}$. The second is the variation of a matrix inverse, $\delta(\mathbf{A}^{-1}) = -\mathbf{A}^{-1}\delta\mathbf{A}\,\mathbf{A}^{-1}$, with $\mathbf{A} = \mathbf{f}$ and $\mathbf{f}^{-1} = \mathbf{F}$.

Now vary $W^{*}(\mathbf{f}) = \det(\mathbf{f})W(\mathbf{f}^{-1})$ by the product rule. The first term gives

\delta(\det\mathbf{f})\,W = J^{-1}W\,\mathbf{F}^{T}\!:\!\delta\mathbf{f}

The second term uses $\partial W/\partial \mathbf{F} = \mathbf{P}$:

\det(\mathbf{f})\,\mathbf{P}\!:\!\delta\!\left(\mathbf{f}^{-1}\right) = -J^{-1}\,\mathbf{P}\!:\!\left(\mathbf{F}\,\delta\mathbf{f}\,\mathbf{F}\right) = -J^{-1}\left(\mathbf{F}^{T}\mathbf{P}\mathbf{F}^{T}\right)\!:\!\delta\mathbf{f}

The last step is the index identity $P_{iA}F_{iB}F_{jA} = (\mathbf{F}^{T}\mathbf{P}\mathbf{F}^{T})_{Bj}$ — write both sides out once and you will never doubt it again. Adding the two terms:

\delta W^{*} = J^{-1}\!\left[W\mathbf{F}^{T} - \mathbf{F}^{T}\mathbf{P}\mathbf{F}^{T}\right]\!:\!\delta\mathbf{f} = J^{-1}\!\left[W\mathbf{I} - \mathbf{F}^{T}\mathbf{P}\right]\mathbf{F}^{T}\!:\!\delta\mathbf{f}

The bracket is the Eshelby energy–momentum tensor.

\frac{\partial W^{*}}{\partial \mathbf{f}} \;=\; J^{-1}\left[\,W\mathbf{I} - \mathbf{F}^{T}\mathbf{P}\,\right]\mathbf{F}^{T} \;=\; J^{-1}\,\boldsymbol{\Sigma}\,\mathbf{F}^{T}, \qquad \boldsymbol{\Sigma} \;=\; W\mathbf{I} - \mathbf{F}^{T}\frac{\partial W}{\partial\mathbf{F}}

$\boldsymbol{\Sigma}$ is the Eshelby energy–momentum tensor, also called the configurational stress. It is the natural stress measure of the inverse problem for exactly the reason its name suggests: it is conjugate to changes in material position rather than spatial position, and material position is what you are solving for.

Trap — get the Eshelby tensor exactly right

The referential Eshelby stress is $\boldsymbol{\Sigma} = W\mathbf{I} - \mathbf{F}^{T}\mathbf{P}$. It is not $W\mathbf{I} - \boldsymbol{\sigma}$ and not $W\mathbf{I} - \mathbf{P}\mathbf{F}^{T}$; and the conjugate stress is $J^{-1}\boldsymbol{\Sigma}\mathbf{F}^{T}$, with the $\mathbf{F}^{T}$ on the right. Transpose and push-forward errors here are endemic, because $\mathbf{f}$ is a two-point tensor whose legs live in different configurations, so a wrong transpose still type-checks dimensionally. Do not conflate $\boldsymbol{\Sigma}$ with the Mandel stress either — different object, different conjugate pairing.

Two Consequences, One Good and One Bad

Implementation. An existing nonlinear finite element code gets a long way toward inverse elastostatics by reinterpreting the nodal unknowns as material coordinates and evaluating the ordinary Cauchy stress from $\mathbf{F} = (\mathbf{f}^{h})^{-1}$: the element technology, the Newton loop and the assembly all survive. That — not a $W^{*}$ material routine — is the residual §5 actually assembles, and the distinction matters: $\partial W^{*}/\partial\mathbf{f}$ is a configurational quantity conjugate to $\mathbf{f}$, so a code returning it from its material routine would need configurational rather than Cauchy boundary data. The two routes are not interchangeable at the same boundary conditions. What does not come for free is the loading side. The residual and tangent must be assembled over the known deformed body with the tractions applied there, so a code that integrates loads over the reference configuration — as most do — needs its load and integration-domain handling rewritten too, not just its material model. Section 5 walks through what that costs.

Well-posedness. The Shield transform preserves polyconvexity, rank-one convexity and ellipticity — plain convexity is the notion it does not preserve. Objectivity and isotropy need a caveat: inversion exchanges them. Substituting $\mathbf{F}\mapsto\mathbf{Q}\mathbf{F}$ gives $\mathbf{f}\mapsto\mathbf{f}\mathbf{Q}^{T}$, so objectivity of $W$ buys isotropy of $W^{*}$, and isotropy of $W$ buys objectivity of $W^{*}$; both survive only when the original law is both objective and isotropic. For the anisotropic laws of §5 that means the structure tensors have to be transformed explicitly, not carried across. The dual material is therefore not pathological, and in the homogeneous, load-independent case the standard existence machinery still applies. The exposure lies elsewhere: growth and coercivity conditions transform too and must be rechecked, and in a real inverse problem the transformed loads and any material inhomogeneity depend on the unknown inverse map — that is what costs you the clean variational structure.

Trap — Shield duality is an identity, not a well-posedness theorem

Shield's theorem relates equilibrium states. It relates equilibrium states; it is not by itself an existence or stability theorem for the loaded inverse problem. Note the direction, though: $W$ polyconvex does give $W^{*}(\mathbf{f}) = \det(\mathbf{f})W(\mathbf{f}^{-1})$ polyconvex, so the difficulty is not a lost convexity notion but the load and inhomogeneity dependence on the unknown map. An impeccably well-posed forward problem can have an inverse problem with no minimiser. In practice this does not announce itself with a visible pathology — it shows up as Newton refusing to converge, which you will be tempted to blame on your mesh.

5. Inverse Elastostatic Finite Elements

Govindjee & Mihalic (1996) turned Shield's identity into a finite element method. The recipe is short, and every step is a re-labelling of something you already have.

Mesh the known deformed body $\mathcal{B}$. The node positions $\mathbf{x}_A$ are data. They never move.
The unknowns are the material coordinates $\mathbf{X}_A$ of those fixed spatial nodes. Interpolate them with the same shape functions, now written as functions of $\mathbf{x}$.
Enforce spatial equilibrium weakly on $\mathcal{B}$, with the Cauchy stress evaluated from $\mathbf{f}^h$.
Solve the nonlinear algebraic system by Newton–Raphson. It costs one nonlinear solve, comparable to a single forward analysis.

\mathbf{X}^{h}(\mathbf{x}) = \sum_{A} N_A(\mathbf{x})\,\mathbf{X}_A, \qquad \mathbf{f}^{h} = \sum_{A} \mathbf{X}_A \otimes \nabla N_A(\mathbf{x}), \qquad \mathbf{F} = \left(\mathbf{f}^{h}\right)^{-1}, \quad J = \frac{1}{\det\mathbf{f}^{h}}

\int_{\mathcal{B}} \boldsymbol{\sigma}\!\left(\mathbf{f}^{h}\right) : \operatorname{grad}\boldsymbol{\eta}\;\mathrm{d}v \;=\; \int_{\mathcal{B}} \rho\,\mathbf{b}\cdot\boldsymbol{\eta}\;\mathrm{d}v \;+\; \int_{\partial\mathcal{B}_t} \bar{\mathbf{t}}\cdot\boldsymbol{\eta}\;\mathrm{d}a

\mathbf{r}_A(\mathbf{X}) = \int_{\mathcal{B}} \boldsymbol{\sigma}\,\nabla N_A\,\mathrm{d}v \;-\; \mathbf{f}^{\text{ext}}_A = \mathbf{0}, \qquad \mathbf{K}_{AB} = \frac{\partial \mathbf{r}_A}{\partial \mathbf{X}_B}

Why the inverse formulation is not the forward one run backwards. The mesh is drawn on the body you actually have; the nodal unknowns are relabelled from spatial to material coordinates; and the load that was the awkward part of the forward problem becomes the easy part.

Key insight — the pressure load becomes trivial

In the inverse formulation the loaded surface is the known deformed surface. A pressure load $\bar{\mathbf{t}} = -p\,\mathbf{n}$ is therefore fully specified before the solve begins: it is a dead load, and the follower-load nonlinearity and pressure load-stiffness term simply do not appear in $\mathbf{K}$. This is a genuine structural advantage over every method that works by repeated forward solves — in those, the pressure genuinely follows the deforming surface, every iteration. The converse of the trap is the one people actually fall into: forgetting that pressure is a follower load inside every forward solve of a Sellier loop (§6), and applying it to the fixed target surface instead of the current iterate's deformed surface, which converges neatly to the wrong fixed point.

Boundary Conditions

Essential boundary conditions are imposed on $\mathbf{X}$: “this node's material coordinate equals this value”. This is how you anchor the reference configuration's otherwise free rigid-body and translation modes — the reference geometry is only determined up to a rigid motion unless you pin it. In vascular work the vessel outlets are typically fixed in space; Ramella et al. (2024) fixed all outlets in their aortic models.

Quasi-Incompressibility, and Locking That Does Not Care Which Way You Solve

Nearly all soft materials of interest are close to incompressible, $\nu \to 0.5$. In the inverse formulation the constraint reads $\det\mathbf{f} = 1$, which is exactly as stiff a constraint as $J = 1$ — it is the same constraint. A purely displacement-type inverse element locks exactly as its forward counterpart does. Govindjee & Mihalic wrote a second paper (1998) for precisely this case; use a mixed $\mathbf{X}$–$p$ or $\bar{\mathbf{F}}$-type formulation with the standard volumetric–isochoric split.

W^{*}(\mathbf{f}) = \det(\mathbf{f})\left[\;\bar{W}\!\left(\bar{\mathbf{F}}\right) + U(J)\;\right], \qquad J = \left(\det\mathbf{f}\right)^{-1}, \quad \bar{\mathbf{F}} = J^{-1/3}\mathbf{F}

Trap — “the mesh doesn't move, so it can't lock”

Locking is not about nodes moving; it is about the element's ability to represent an isochoric deformation field with the interpolation it has. $\det\mathbf{f} = 1$ constrains the discrete field the same way $J = 1$ does. If you would not use a plain displacement element for a forward incompressible problem, do not use its inverse twin either.

Anisotropy, and Where the Fibres Live

Fachinotti, Cardona & Jetteur (2008) formulate the inverse design problem for large-deformation anisotropic hyperelasticity. Their formulation expresses equilibrium on the distorted (known, deformed) configuration while keeping a standard constitutive library written in the usual undistorted-configuration measures, so that — in their words — “the modifications to a standard finite elements code are limited to the routines for the computation of the finite element internal forces and tangent matrix.” That is a strong practical statement: the inverse capability is an element-level change, not a solver rewrite. An inverse-elastostatic capability also ships in commercial FE (ANSYS Mechanical), which is how Ramella et al. (2024) used it.

Trap — anisotropy bookkeeping

Fibre directions, collagen orientations and any structure tensor are defined in the reference configuration — the one you are solving for. If you take a fibre field from an image or an atlas, decide explicitly whether it is a material direction (which must be pulled back with $\mathbf{f}$) or a spatial one. Treating an imaged fibre field as though it already lived in the reference configuration silently changes the constitutive response, and the resulting error is invisible in the residual: the solve converges, to the wrong body.

6. Fixed-Point Unloading: Relaxation, Acceleration, Continuation

Suppose you cannot modify the constitutive kernel of your solver — commercial code, contact, complex multiphysics, a validated pipeline you are not allowed to touch. You can still recover the reference configuration using only forward solves. Sellier (2011) proposed the fixed-point scheme now used widely in cardiovascular and cardiac mechanics.

The Sellier iteration

Let $\mathbf{x}^{*}$ be the measured deformed nodal positions. Initialise $\mathbf{X}_0 = \mathbf{x}^{*}$. At iteration $k$: (i) run a forward analysis using $\mathbf{X}_k$ as the reference mesh under the true loads, obtaining deformed positions $\mathbf{x}_k$; (ii) form the residual $\mathbf{R}_k = \mathbf{x}_k - \mathbf{x}^{*}$; (iii) update $\mathbf{X}_{k+1} = \mathbf{X}_k - \alpha_k\,\mathbf{R}_k$. Stop when $\max\lVert\mathbf{R}_k\rVert < \varepsilon$. The sought reference configuration is precisely the fixed point of $\mathbf{G}(\mathbf{X}) = \mathbf{X} - \alpha\left(\boldsymbol{\varphi}_{\mathbf{X}}(\mathbf{X}) - \mathbf{x}^{*}\right)$.

\mathbf{R}_k = \mathbf{x}_k - \mathbf{x}^{*}, \qquad \mathbf{x}_k = \boldsymbol{\varphi}_{\mathbf{X}_k}\!\left(\mathbf{X}_k\right), \qquad \mathbf{X}_{k+1} = \mathbf{X}_k - \alpha_k\,\mathbf{R}_k

Because this is a Picard iteration, its asymptotic behaviour is governed by the contraction factor of the update map:

\left\lVert \mathbf{X}_{k+1} - \mathbf{X}^{\star} \right\rVert \;\le\; L\,\left\lVert \mathbf{X}_{k} - \mathbf{X}^{\star} \right\rVert, \qquad L = \sup_{\mathbf{X}\in\,U}\left\lVert \mathbf{I} - \alpha\,\frac{\partial \mathbf{R}}{\partial \mathbf{X}}(\mathbf{X}) \right\rVert

Read the definition of $L$ carefully: it is a bound holding uniformly over a convex neighbourhood $U$ that contains the iterates, which is what the mean-value argument needs. Everything that follows in this section is an attempt to keep $L$ below one without knowing $\partial\mathbf{R}/\partial\mathbf{X}$. Keep two numbers apart while you do it, because they are routinely conflated and they are not close:

\rho\!\left(\mathbf{I} - \alpha\,\frac{\partial \mathbf{R}}{\partial \mathbf{X}}\right) \;\le\; L, \qquad \rho < 1 \ \Longrightarrow\ \text{local, asymptotic convergence only}

A uniform bound $L < 1$ over such a neighbourhood gives the per-step guarantee written above: every step shrinks the error. The spectral radius $\rho$ gives only the asymptotic rate once the iterate is already near $\mathbf{X}^{\star}$, and it can be comfortably below one while the norm is several times larger — in the laboratory below, at the default settings, $\rho = 0.39$ against $\lVert\mathbf{I}-\partial\mathbf{R}/\partial\mathbf{X}\rVert_2 = 2.11$. Note what that second number is not: it is the norm at the recovered $\mathbf{X}^{\star}$ alone, a single-point evaluation, and the same iteration converges monotonically from $\mathbf{X}_0 = \mathbf{x}^{*}$ in eight forward solves anyway. Evaluated at one point the norm is as local a statement as $\rho$ is; only a bound holding across the iterates certifies anything per step. So a green $\rho$ is not a promise that an iteration started at $\mathbf{X}_0 = \mathbf{x}^{*}$ will get there: at $p = 150$ mmHg and $\mu = 150$ kPa the widget prints $\rho = 0.92$ next to an $\alpha = 1$ trace badged DIVERGED. Both statements are correct. The panel prints $\rho$, and labels it $\rho$.

Relaxation: $\alpha = 1$ Is a Choice, Not a Law

Sellier fixed $\alpha = 1$. Rausch, Genet & Humphrey (2017) showed that $\alpha = 1$ is not necessarily the best choice, that the optimal $\alpha$ is problem dependent and not predictable a priori — in their study the iteration count depends strongly on $\alpha$ and is never better than their augmented method. Their fix is Aitken's $\Delta^2$ dynamic relaxation, which makes convergence almost independent of the initial $\alpha$; they validated it on four cases — a murine venous thrombus, a healthy mouse aortic segment, an ovine anterior mitral valve leaflet, and a murine aortic dissection with intramural thrombus.

\alpha_k \;=\; -\,\alpha_{k-1}\,\frac{\mathbf{R}_{k-1}\cdot\left(\mathbf{R}_k - \mathbf{R}_{k-1}\right)}{\left(\mathbf{R}_k - \mathbf{R}_{k-1}\right)\cdot\left(\mathbf{R}_k - \mathbf{R}_{k-1}\right)}

Implementation note for the laboratory below: the widget clips $\alpha_k$ to $[-2.5, 2.5]$ so a single wild update cannot throw the trace off the plot. That clip fires for roughly one setting in fifty of the sliders' range — softer walls at high pressure with a large opening angle. Unclipped Aitken, which is what the cited papers use, can take much larger steps, so read the orange trace as a lightly projected variant of the formula above rather than the formula itself.

Trap — the iteration is not automatically convergent

Presenting Sellier's iteration with $\alpha = 1$ as unconditionally convergent is a genuine error, not a simplification. Rausch et al. (2017) explicitly showed the optimum is problem-dependent and unpredictable a priori; Marx et al. (2022) documented that even the Aitken-augmented approach still diverged for some of their left-ventricular models. Aitken $\Delta^{2}$, Armijo damping and Anderson acceleration are corrections to a real failure mode, not decorations on a method that works anyway. The laboratory below lets you drive a configuration into divergence and then rescue it.

Line Search and Anderson Acceleration

Marx et al. (2022) added an Armijo line search that damps unfavourable search directions. They report that for some of their LV models the Aitken-augmented approach still diverged, whereas their algorithm was robust for all cases. On a 19-patient left-ventricular cohort (12 aortic-valve-stenosis, 7 aortic-coarctation) it converged within 10 iterations for every case, all 19 run with the same reduced Holzapfel–Ogden law. Keep that separate from their constitutive-law comparison, which was a different experiment: eight laws (Guccione, Usyk, general / reduced / original Holzapfel–Ogden, HO with dispersion, one-fibre HO, Demiray) applied to just two representative cases, the best fit (06-CoA) and the worst (03-AS). The cohort run took 10.85–74.18 min per case for unloading plus parameter fitting against 6.31–24.50 min for the corresponding validation inflation — roughly a small multiple of a single forward solve — on meshes up to 1.5 M elements, with joint identification of passive parameters.

Barnafi, Petras & Gerardo-Giorda (2024) extend the reference-configuration problem to fully nonlinear poroelasticity, and accelerate it with Anderson acceleration, reporting a reduction in the number of iterations of up to 80%. Attach that result to the right object: their fixed point is not Sellier's geometry residual. They rewrite the equations in terms of the reference porosity and define a time-dependent problem whose steady state is that porosity, and it is the pseudo-time map to the steady state that Anderson accelerates. Anderson acceleration itself is also a different kind of correction from everything else in this section: it mixes a history of $m$ previous iterates and residuals, rather than choosing a scalar $\alpha_k$ for the current step, so it does not slot into the “choose $\alpha_k$” box below. It is a general fixed-point accelerator and there is no obstacle to using it on the Sellier map — but that is a transfer, not something this paper reports. That paper also identifies an inconsistency in primal formulations caused by second-order derivatives of the displacement in the nonlinear mass-conservation equations, resolved with mixed formulations.

Load Continuation

There is a third knob that costs nothing to try and is not a research contribution in any of the papers above — it is ordinary nonlinear-solver practice, and it is worth naming because it rescues cases that no relaxation parameter will. Instead of solving the unloading problem at the full load in one go, solve a sequence of problems at $p_1 < p_2 < \dots < p_N = p$, each one starting from the previous answer. Every sub-problem has the same target $\mathbf{x}^{*}$; only the load level changes. With no residual stress and no retained loading, the $p \to 0$ answer is $\mathbf{X} = \mathbf{x}^{*}$, so the first sub-problem starts essentially at its own solution, and each subsequent one starts inside the basin of attraction of the next. Continuation does not change the contraction factor at the solution — it changes whether your iterate ever gets close enough for that contraction factor to matter.

Read that hypothesis, because it is the one people drop. If the model carries prestrain — an opening angle, a growth field, an initial stress — then the reference that maps to $\mathbf{x}^{*}$ at zero pressure is not $\mathbf{x}^{*}$, and initialising at the image starts you outside the very basin continuation was supposed to hand you. In the laboratory below, at $\Theta_0 = 60^\circ$ the $p\to0$ reference is $R_i = 12.15$, $R_o = 13.65$ mm against an image of $10.00 / 11.50$ mm. Continuation must be warm-started from a prestrain-consistent zero-load state, which is what the widget does. Nor does small-step continuation come with a guarantee: it keeps the sequence in the basin along a regular branch for sufficiently small load steps, and near a limit point no finite step count secures that.

The black-box unloading loop and its accelerations. Every box marked FORWARD is a complete nonlinear analysis; the acceleration literature is entirely about reducing how many of them you need and making sure the sequence converges at all.

Why it matters — the trade-off, stated honestly

Inverse elastostatics: one nonlinear solve, exact treatment of pressure, but you need access to the element/constitutive layer. Sellier family: $N$ forward solves, works around any black box, handles contact and plasticity and multiphysics, but may diverge and must handle the follower pressure correctly at every iteration. Do not describe Govindjee–Mihalic as “iterative” because Newton iterates inside it — that is the inner loop that every nonlinear solve has, and comparing it to Sellier's outer loop misrepresents the cost difference by an order of magnitude.

Interactive Tool — Unloaded-Shape Laboratory

An incompressible neo-Hookean thick-walled tube under internal pressure, plane strain, solved analytically — no FE, no mesh, nothing hidden in a discretisation. Read the numbers with the right precision attached: the equilibrium radii and pressures are bisected to machine tolerance, but the conditioning number, the spectral radii and $\alpha^{\star}$ are finite-difference and grid-search estimates (central difference with $h = 10^{-4}R_i$; $\alpha^{\star}$ scanned on a 0.01 grid over the relaxation slider's own range, so that whatever the panel recommends you can actually select), and the solve counts are iterative outcomes. The panel labels which is which. The imaged inner and outer radii $a^{*}, b^{*}$ are the data; the reference radii are the unknowns. Neo-Hookean is a demonstrator, not a model of arterial tissue: real vessel wall stiffens steeply with stretch, so read the stretch levels here as a numerical experiment rather than a physiological prediction. Drive the load up and the stiffness down and watch four separate things break.

Internal pressure $p$ (mmHg) 80

Shear modulus $\mu$ (kPa) 150

Imaged inner radius $a^{*}$ (mm) 10.0

Imaged wall $h^{*} = b^{*}-a^{*}$ (mm) 1.50

Opening angle $\Theta_0$ (°) 0

Relaxation $\alpha$ 1.00

Load-continuation steps $N$ 1

Aitken $\Delta^2$ trace show “image is stress-free” model

Try this: (1) At the defaults, read the “apply $-p$” error in the panel: about $+3.6\%$ — modest, but a whole voxel and more on a 10 mm radius. Now drag $\mu$ down toward 60 kPa, or $p$ up toward 180 mmHg, and watch it pass $20\%$. Both numbers are computed, not quoted: the trick's error has no universal magnitude. Note the second figure beneath it, the CT-as-reference over-inflation — a different and much larger quantity, and the one people accidentally report. (2) Set $\Theta_0 = 60^\circ$. The shape of the load-free tube barely notices, but two things move: the recovered reference and the load-free radii separate (they were identical at $\Theta_0 = 0$, which is the whole content of “load-free vs stress-free”), and the residual hoop-stress curve in the middle panel becomes compressive at the inner wall and tensile at the outer. The two hoop curves now sit on different spans of the same true $r/a^{*}$ axis, because the load-free wall really is at smaller radii. No unloading algorithm recovered $\Theta_0$; you supplied it. (3) Soften the wall and raise the pressure until the $\alpha = 1$ trace in the convergence plot flattens or turns up. Then watch the Aitken trace. Then set $\alpha$ by hand near the computed $\alpha^{\star}$ shown in the panel — when there is one. Soften far enough ($\mu = 10$ kPa at $p = 80$ mmHg, say) and the panel reports instead that no $\alpha$ in the slider's range brings $\rho$ below one, and every trace fails on its first forward solve: at that stiffness the imaged tube used as its own reference cannot support 80 mmHg at any radius, so the iteration has nowhere to start. The direct inverse solve still returns a reference configuration. That contrast — one nonlinear solve succeeding where the whole fixed-point family has no admissible starting point — is worth sitting with. Nothing here is hard-coded; the computation decides what diverges. (4) With a diverging $\alpha = 1$ case on screen, raise the continuation steps $N$. Continuation buys you the basin of attraction, at the price of more forward solves; the trace is plotted against cumulative forward solves so the price is visible. Watch for the BUDGET SPENT badge: at $p = 160$ mmHg, $\mu = 150$ kPa, with $\alpha = 1$ and everything else at its default, $N \ge 4$ uses all 60 of the run's allowed solves and still stops short of the full load, and the panel says which load level it actually reached. The continuation run uses whatever $\alpha$ you have set, so that outcome is a statement about the pair $(\alpha, N)$, not about $N$ alone — the same $N = 4$ at $\alpha = 0.68$ reaches the full load in 20 solves. A run that stops because it ran out of budget has not converged, and the widget will not pretend otherwise. (5) Compare the two cost numbers in the panel: one nonlinear solve for the direct route, versus the forward-solve count the fixed point actually needed.

7. Prestressing Instead of Un-Deforming

Often you do not actually want the unloaded geometry. You want the in vivo stress state on the imaged mesh — because that is what a device-deployment simulation or a rupture-risk calculation starts from. Then you can leave the geometry alone and install a prestress instead. This is a different problem with different failure modes, and conflating the two is a standard error.

MULF — Modified Updated Lagrangian Formulation

Keep the imaged geometry as the FE reference mesh and seek a prestretch field $\mathbf{F}_{\text{pre}}(\mathbf{x})$ such that the imaged geometry is in equilibrium with the in vivo load. The total deformation splits multiplicatively, $\mathbf{F} = \hat{\mathbf{F}}\,\mathbf{F}_{\text{pre}}$. The algorithm iterates: solve for displacement increments under the in vivo load, accumulate the increment into $\mathbf{F}_{\text{pre}}$, reset the geometry to the image, repeat until the displacement increment vanishes. At convergence $\hat{\mathbf{F}} = \mathbf{I}$ at the imaged shape and the stress is nonzero. Gee, Förster & Wall (2010).

\mathbf{F} \;=\; \hat{\mathbf{F}}\,\mathbf{F}_{\text{pre}}, \qquad \hat{\mathbf{F}}\big|_{\text{converged}} = \mathbf{I}\ \text{ on the imaged geometry}

\int_{\mathcal{B}_{\text{image}}} \boldsymbol{\sigma}\!\left(\hat{\mathbf{F}}\,\mathbf{F}_{\text{pre}}\right) : \operatorname{grad}\boldsymbol{\eta}\;\mathrm{d}v \;=\; \int_{\partial\mathcal{B}_t} \bar{\mathbf{t}}\cdot\boldsymbol{\eta}\;\mathrm{d}a

Two routes from the same imaged geometry. They answer different questions and fail in different ways; the choice is made by what you have to hand over at the end, not by which one is more elegant.

The Decision Rule

Geometry is the deliverable (turbine cold shape, mold for a silicone skin, 3D-print compensation) ⇒ you need inverse elastostatics. Prestressing cannot help you: it never changes the shape.
Stress on the imaged geometry is the deliverable (device deployment, rupture risk) ⇒ prestressing is simpler, preserves the mesh quality you paid a segmentation specialist for, and avoids the risk of producing a self-intersecting reference mesh.
Both ⇒ inverse elastostatics; the prestress falls out for free as $\boldsymbol{\sigma}(\mathbf{f})$ evaluated on the image.

Key insight — the two are not numerically equivalent

MULF converges to a prestretch field that need not be compatible with any global reference configuration; the inverse solve returns an actual configuration. For a residually stressed organ that incompatibility is a feature, not a bug — a residually stressed body has no single stress-free configuration to be compatible with (§3). But it also means you cannot “read off” the unloaded geometry from a converged MULF prestretch field, and you should not report one.

8. Uniqueness, Conditioning, and Static Determinacy

When is the inverse motion problem solvable, and how badly conditioned is it? Four separate answers, of which only the last is good news.

Existence of the Inverse Map

$\boldsymbol{\varphi}^{-1}$ exists only if $\boldsymbol{\varphi}$ is a diffeomorphism: $\det\mathbf{F} > 0$ everywhere and $\boldsymbol{\varphi}$ globally injective. In the discrete inverse solve the first requirement is usually screened at the quadrature points, and the second is a global test that screening does not cover.

\underbrace{\det \mathbf{f}^{h}(\mathbf{x}_q) > 0 \ \ \forall q}_{\text{quadrature screen only}}, \qquad \kappa\!\left(\mathbf{K}_{\text{fwd}}\right) = \lVert \mathbf{K}_{\text{fwd}}\rVert\,\lVert \mathbf{K}_{\text{fwd}}^{-1}\rVert \;\longrightarrow\; \infty \ \text{ at a limit point of the forward branch}

Take that screen for the real condition and it will let an inverted element through. Sampling $\det\mathbf{f}^h$ at quadrature points certifies the whole element only when the Jacobian is constant over it — affine simplices. For a bilinear, trilinear or higher-order map the determinant varies inside the element and can be positive at every Gauss point and negative between them. The bilinear map $x = \xi$, $y = (0.7 + 0.8\xi)\eta$ on $[-1,1]^2$ has $\det = 0.7 + 0.8\xi$: at the four $2\times2$ Gauss abscissae it takes the values $0.238$ and $1.162$, comfortably positive, while at the $\xi = -1$ edge it is $-0.1$. What actually rules inversion out is an elementwise certificate — Bézier/Bernstein coefficient bounds on the determinant, or adaptive subdivision until the sign is decided — on top of a genuine global injectivity test. Sampling the determinant at the nodes is not that certificate either, except for element classes where the determinant is provably affine over the element: a trilinear hexahedron can have $\det\mathbf{f}^h > 0$ at all eight corners and $\det\mathbf{f}^h < 0$ inside.

Mind the subscript. $\mathbf{K}_{\text{fwd}}$ is the forward equilibrium tangent, not the inverse tangent $\mathbf{K}_{AB} = \partial\mathbf{r}_A/\partial\mathbf{X}_B$ of §5, and a limit point is a property of the forward load–displacement branch. The matrix the inverse solve is actually conditioned on is $\partial\mathbf{x}/\partial\mathbf{X}$ — a third matrix again, and the one this section later insists you check rather than infer.

Trap — the Newton solve will happily converge to a non-physical body

At large strain the recovered reference mesh can invert ($\det\mathbf{f}^h \le 0$) or self-intersect, and neither shows up as a failure of the nonlinear solve. Check $\det\mathbf{f} > 0$ everywhere — which on non-affine elements means an elementwise bound or subdivision, not a Gauss-point sample — and test the recovered reference surface for self-intersection. A residual of $10^{-10}$ is not evidence that the answer is a body.

Loss of Variational Structure

As established in Section 4, the Shield transform preserves polyconvexity, so $W^{*}$ inherits it; the exposure is in coercivity/growth and in loads that depend on the unknown map. The inverse problem can therefore still lack a minimiser even when the forward problem is impeccably well posed. In practice this shows up as Newton failing to converge, not as a visible pathology — which is why it gets misdiagnosed as a meshing problem.

Non-Uniqueness Through the Forward Map

If the forward problem admits multiple equilibria for the same load — buckling, snap-through, wrinkling of a thin skin — treat it as a warning sign, and be careful about what it actually licenses. Forward multiplicity says one reference can produce several images. That permits, but does not imply, the thing the inverse problem cares about: several references producing one fixed image. The latter is non-injectivity of the reference-to-image map, and it has to be tested, not inferred. Likewise, a forward tangent $\mathbf{K}_{\text{fwd}}$ approaching singularity at a limit point does not by itself establish that $\partial\mathbf{x}/\partial\mathbf{X}$ — the Jacobian the inverse problem is actually conditioned on — is near-singular; the two are related but not the same matrix. What you must check is injectivity of the reference-to-image map and the conditioning of $\partial\mathbf{x}/\partial\mathbf{X}$ at the recovered solution, which is the number the laboratory above prints as $\kappa(\partial\mathbf{x}/\partial\mathbf{X})$. Near a limit point expect both to degrade together; do not assume it.

Geometry and Material Are Not Jointly Identifiable

This is the nonuniqueness that bites in practice, and it is structural rather than numerical. For essentially any assumed stiffness there exists a reference geometry that maps to the imaged shape under the given load: soften the material and the recovered reference shrinks to compensate, exactly. A single loaded state therefore determines a curve of $(\text{material}, \text{reference geometry})$ pairs, all of which reproduce the image perfectly.

Eliminating the degeneracy needs a second, genuinely independent load–response measurement: an image at a second pressure. What is usually available instead is weaker, and worth distinguishing. Marx et al. (2022) take image data at one time point plus a single end-diastolic pressure–volume point, then generate an empirical Klotz EDPVR curve from that one point and fit against it — a population-level prior rather than a second measured load state. Note how much of that one point is itself measured: the LV pressure was obtained by invasive catheterisation for the seven coarctation cases, while for the twelve aortic-stenosis cases it was assigned empirically from a reference pool of 290 patient records. So for most of the cohort the prior enters twice over. That regularises the degeneracy; it does not remove it, and the recovered stiffness inherits whatever the Klotz population relation assumes. The explorer below shows the degenerate curve and what a real second load state does to it.

Interactive Tool — Joint Identifiability Explorer

Same analytic tube. A “true” wall with shear modulus $\mu_{\text{true}}$ produces the imaged geometry at $p_1 = 80$ mmHg. Now pretend you do not know $\mu$. For every candidate $\mu$ the inverse solve returns a reference radius $A(\mu)$ that reproduces the image exactly — zero misfit, every time. Then add a second image at pressure $p_2$ and watch what survives.

True $\mu_{\text{true}}$ (kPa) 150

Second pressure $p_2$ (mmHg) 140

Misfit tolerance / noise floor (mm) 0.05

Load states available

Try this: (1) Switch to “One image”. The misfit is identically zero across the whole $\mu$ axis: every stiffness in the range is a perfect fit, with its own reference geometry. Nothing about the data prefers the truth. (2) Switch back to two images and set $p_2 = 80$ mmHg — the same pressure as the first. The misfit collapses to zero again. A second image at the same load state is not a second load state. (3) Push $p_2$ far from $p_1$: the misfit valley narrows and the interval printed in the panel shrinks. Separation in load, not number of images, is what buys identifiability. (4) Raise the tolerance. The slider moves the acceptance threshold, not the data — no radius is perturbed — but the interval where the misfit is indistinguishable from measurement error widens, which is the same statement: identifiability depends on data quality, not only on experiment design. (5) Drag $p_2$ up from 80 mmHg and watch the hatched zone grow from the soft end. It opens just under 100 mmHg, and at the widget's own default $p_2 = 140$ mmHg it already covers about a third of the $\mu$ axis — you are not waiting for it to appear, you are watching it spread. Those candidate moduli have no equilibrium at all at $p_2$, and not because the tube passes a limit point: for this model $p(a)$ climbs monotonically to a finite ceiling $p_\infty = \mu\ln(R_o/R_i)$, so a soft enough candidate never reaches $p_2$ at any radius. An asymptote, not a fold. A failed solve is a rejection, not a perfect fit, and the green band never crosses the hatching.

Static Determinacy — the One Piece of Good News

For a thin pressurised membrane the equilibrium equations close without any constitutive law. With negligible bending resistance there are three equilibrium component equations for the three in-plane stress components, driven by geometry and pressure alone:

\frac{\sigma_{1}}{\rho_{1}} + \frac{\sigma_{2}}{\rho_{2}} \;=\; \frac{p}{h}

\text{sphere } \rho_1=\rho_2=R \Rightarrow \sigma = \tfrac{pR}{2h}; \qquad \text{cylinder } \rho_2\to\infty \Rightarrow \sigma_\theta = \tfrac{pR}{h}

Zhao, Raghavan & Lu (2011) state that under traction boundaries the wall stress is completely independent of the material properties, and verified it numerically: the forward reference solution used an anisotropic Holzapfel model while the inverse analysis used a neo-Hookean auxiliary model, and the recovered principal stresses differed by less than 1% through most of the aneurysm sac — while a 100-fold increase in the auxiliary parameters moved the first and second principal stresses by at most 0.4% and 0.45% respectively, and then only near the boundary.

The exception is a boundary layer of roughly four element layers near the clamped edges, where the stress does become material-dependent and where the second principal stress deviates by more than 10%; the authors exclude that layer from identification.

Key insight — static determinacy giveth and taketh away

This is why aneurysm wall-stress estimates are trusted despite unknown patient-specific material properties: get the geometry and the pressure right and the stress follows, whatever the constitutive law. It is simultaneously why stress alone can never identify stiffness — the same independence that makes the stress robust makes it uninformative about the material. It is also what makes Zhao et al.'s two-stage pointwise identification possible: compute stress with an arbitrary auxiliary model, then regress stress against strain independently at each material point, so the problem size depends on the constitutive model rather than the mesh and no globally coupled optimisation is needed.

9. Practice: Choosing a Method, and the Pitfalls That Bite

Selection

Method	Cost	Needs	Returns	Main risk
Inverse elastostatics (Govindjee–Mihalic)	One nonlinear solve	Access to the element/constitutive layer, or a code that ships it (e.g. ANSYS Mechanical)	Geometry and stress; exact treatment of pressure loads	Coercivity/growth and map-dependent loading are not inherited, so the inverse variational problem may have no minimiser (§4); recovered mesh may invert
Fixed point (Sellier + Aitken / Armijo / Anderson)	$N$ forward solves, with no bound in general — $\le 10$ with Armijo on the 19 LV cases of Marx et al. (2022)	Nothing but a forward solver	Geometry, and stress from the final forward solve	May diverge; follower pressure must be re-applied every iteration
MULF prestressing (Gee et al.)	Cheapest	Ability to accumulate a prestretch field	Stress only; mesh preserved exactly	Prestretch need not correspond to any configuration
Learned surrogate	Amortised: seconds of inference after an offline training campaign	A large simulation dataset in the target family	Geometry, at inference speed	Only as trustworthy as the training distribution

Convergence Criteria People Actually Use

Rausch et al. (2017) test $\max\lVert\mathbf{R}_k\rVert < \varepsilon$ at $\varepsilon = 0.01$, $0.001$ and $0.0001$ mm. Marx et al. (2022) require four conditions simultaneously: max nodal error $< 0.1$ mm; an end-diastolic volume residual and a Klotz zero-pressure-volume residual, both below $0.5\%$ of the measured end-diastolic volume; and a parameter-update magnitude $< 0.001$. The two volume tests are distinct, and collapsing them into one “volume difference” drops the $V_0$ check — which is the term doing the regularising work in §8.

\begin{aligned}\text{stop when}\quad & \max_A \left\lVert \mathbf{x}^{k}_A - \mathbf{x}^{*}_A \right\rVert < 0.1\ \text{mm}, \qquad \max\left(\lvert 1 - a_{\text{scale}}\rvert,\ \lvert 1 - b_{\text{scale}}\rvert\right) < 10^{-3},\\[3pt] & \left\lvert V^{\text{dat}}_{\text{ed}} - V^{\text{sim}}_{\text{ed}} \right\rvert < 0.005\,V^{\text{dat}}_{\text{ed}}, \qquad \left\lvert V^{\text{Klotz}}_{0} - V^{\text{sim}}_{0} \right\rvert < 0.005\,V^{\text{dat}}_{\text{ed}}\end{aligned}

Why it matters — use a geometric tolerance in physical units

Both groups stop on a distance in millimetres, not on a relative residual drop. A relative criterion says nothing about whether you have converged to something the imaging modality could even resolve: a six-order-of-magnitude residual drop can still leave you sub-voxel nonsense, and a loose relative tolerance can stop you a whole voxel away. Millimetres are comparable to the voxel size; residual ratios are comparable to nothing.

The Pitfall List

Locking. $\det\mathbf{f} = 1$ is as stiff a constraint as $J = 1$. Use mixed or $\bar{\mathbf{F}}$ inverse elements (Govindjee & Mihalic, 1998). (§5)
Follower loads. Exactly dead in the inverse formulation because the surface is known; genuinely following in every forward solve of the Sellier loop. Re-apply the pressure normal to the current iterate's deformed surface or you converge to the wrong fixed point. (§5, §6)
Traction pull-back. A pressure on the current surface and a dead traction on the reference surface are different boundary conditions related by Nanson's formula. Getting this wrong changes which problem you solved. (§2)
Anisotropy bookkeeping. Fibre directions live in the reference configuration you are solving for; decide whether an imaged or atlas fibre field is a material or a spatial direction. Fachinotti et al. (2008) formulate the anisotropic inverse design problem precisely for this setting. (§5)
Joint identifiability. Material parameters and the reference configuration are not jointly identifiable from a single loaded state. Get an image at a second pressure; failing that, say plainly that a population prior such as the Klotz EDPVR curve is regularising the fit rather than identifying it. (§8)
Mesh sanity. Check $\det\mathbf{f} > 0$ — over each element, not only at its quadrature points, unless the elements are affine — and test the recovered reference surface for self-intersection. At large strain this is a real failure mode, not a formality. (§8)
Regularisation of segmentation noise. The imaged surface carries voxel-scale noise and the inverse solve amplifies its high-frequency content into the reference geometry. Smooth the target surface before, not the answer after.
Unfitted / immersed geometry. For heterogeneous cross-sections — calcified plaque, lipid pools — fitting a conforming mesh is often the real bottleneck rather than the inverse solve.
Cost conflation. One nonlinear solve versus $N$ forward solves. Newton's inner iterations are not the same level of iteration as Sellier's outer loop. (§6)

Symbol collision — two different $\bar{\mathbf{F}}$

In the volumetric–isochoric split, $\bar{\mathbf{F}} = J^{-1/3}\mathbf{F}$ is the isochoric part of the deformation gradient (§5). In F-bar element technology, $\bar{\mathbf{F}} = (J_0/J)^{1/3}\mathbf{F}$ is a stabilised gradient built from the element-centroid Jacobian $J_0$. Both appear in this module's subject matter and they are different objects. Check which one a paper means before porting a formula.

\bar{\mathbf{F}} = \left(\frac{J_0}{J}\right)^{1/3}\mathbf{F} \quad \text{(F-bar \textit{stabilisation})} \qquad \ne \qquad \bar{\mathbf{F}} = J^{-1/3}\mathbf{F} \quad \text{(isochoric part, \S5)}

10. Applications and the 2023–2026 Frontier

Vascular

Lu, Zhou & Raghavan (2007) gave the founding demonstration: inverse elastostatic stress analysis of abdominal aortic aneurysms from the CT (pre-deformed) geometry, for membrane-like biological structures, removing the conventional assumption that the imaged geometry is stress-free. Zhao, Raghavan & Lu (2011) extended it to pointwise identification of heterogeneous anisotropic properties in cerebral aneurysms, using static determinacy to decouple the stress computation from the material regression.

Device Deployment

Ramella et al. (2024) quantified what happens if you skip the step. The zero-pressure configuration was obtained by applying 80 mmHg diastolic pressure on the internal lumen of the segmented CT geometry with the inverse elastostatic method implemented in ANSYS Mechanical, all vessel outlets fixed in space, full Newton–Raphson. Against the CT-as-reference result, the pre-stressed model predicted contact pressure higher by 48% (TAVI) and 55% (TEVAR); aortic wall stresses higher by 162% / 157%; aortic wall strains higher by 77% / 21%. Valve-stent stresses increased by 48% on average with the pre-stressed aorta in the TAVI simulations, while TEVAR stent-graft stresses increased by only about 2.5% on average.

Why it matters — a methods argument became a design requirement

Note which quantities moved and which did not. The stent-graft stresses in the TEVAR model barely changed; the vessel-side quantities changed by more than a factor of two. Prestress is not a uniform correction factor you can fold into a safety margin — it redistributes load between the device and the tissue, which is exactly the quantity a device designer is optimising.

Cardiac, and the Learned Shortcut

Marx et al. (2022): a 19-patient LV cohort, reduced Holzapfel–Ogden constitutive law, joint recovery of unloaded reference geometry and passive parameters, within 10 Armijo-damped fixed-point iterations, 10.85–74.18 min per case on meshes up to 1.5 M elements.

Mu, Chan & Yap (2025) then replaced the iteration entirely with HeartUnloadNet, a graph-attention network mapping an end-diastolic mesh plus diastolic pressure, stiffness scaling and fibre orientation directly to the unloaded mesh. Trained and tested on a 20,700-simulation FE dataset — split by anatomy, 42 shapes / 13,890 samples for training against 18 unseen shapes / 6,810 samples held out — and weakly supervised by cycle consistency (loading and unloading learned bidirectionally in one model), it reaches DSC 0.986 and Hausdorff distance 0.083 cm at 0.02 s per case — over $10^{5}\times$ faster than the iterative solver — and still reaches about 97% DSC when trained on only 200 samples.

Key insight — what a surrogate does and does not inherit

A learned unloading map amortises the cost of the iteration, not the identifiability of the problem. The network is trained on FE simulations whose reference configurations were generated, so the training pairs are exact by construction; nothing in that pipeline resolves the geometry–material degeneracy of §8, which is why HeartUnloadNet takes the stiffness scaling as an input. Speed-ups of $10^{5}$ are real and they are not a substitute for a second load state.

Turbomachinery

Fachinotti, Cardona & Jetteur (2008) pose the inverse design problem as “determining the initial shape of the piece, such that it attains the designed shape under the effect of service loads”, targeted at markedly anisotropic parts “like laminated turbine blades”, and demonstrate it on an industrially relevant turbine-blade design example — the cold, manufacturing shape that becomes the aerodynamically correct hot shape. This is the cleanest statement of the case where only inverse elastostatics will do: nobody can machine a prestress field.

Soft Robotics and Tactile Skin

The mold-design problem for a vision-based tactile sensor is the same problem with a different vocabulary: cut the mold so that the cured, gravity-loaded, membrane-pretensioned skin takes the designed geometry.

\text{mold design:}\quad \text{find } \mathcal{B}_0 \ \text{ s.t. } \ \boldsymbol{\varphi}_{\mathcal{B}_0}\!\left(\mathcal{B}_0;\ \rho\mathbf{g},\ p_{\text{pretension}}\right) = \mathcal{B}_{\text{design}}

A sibling problem holds the geometry fixed and varies the material instead. Awasthi, Lamperski & Kowalewski (2023) solve it: given the undeformed shape, the external load, and a desired final shape, find the spatially varying material fields (Mooney–Rivlin $c_1, c_2$) that realise it — a simultaneous-analysis-and-design formulation with analytical gradients for both the objective and the continuum-mechanics constraints, interior-point optimisation with a BFGS Hessian approximation, and a choice of $L^2$ (smooth) or TVD (sharp-contrast) regularisation. It is the design-side twin of a material-identification problem, and it is listed here as a secondary application: the configuration inverse and the material inverse are different problems (§1), and mixing them is how projects lose a year.

Where the Field Is

Two things have happened since 2023. The fixed-point family got robust: Sellier (2011, $\alpha = 1$) → Rausch, Genet & Humphrey (2017, Aitken $\Delta^2$) → Marx et al. (2022, Armijo line search, robust on all 19 LV cases where the Aitken-augmented approach still diverged for some). That lineage ends at Marx: Barnafi, Petras & Gerardo-Giorda (2024) belong to a parallel development rather than the next rung of it, extending the reference-configuration problem to fully nonlinear poroelasticity and applying Anderson acceleration to the pseudo-time map whose steady state is the reference porosity — not to Sellier's geometry residual — with iteration counts reduced by up to 80%. And the iteration started being replaced outright by amortised surrogates (Mu, Chan & Yap, 2025).

The honest summary: for smooth single-phase hyperelasticity the configuration-inverse problem is routinely tractable — one inverse-elastostatic solve, or, on the cohorts reported so far (Marx et al.'s 19 LV cases with one constitutive family), within about ten accelerated fixed-point iterations. Do not upgrade that to “solved”: there is still no convergence bound for the Sellier family in general, which is why this page keeps saying $N$ forward solves and no guarantee. The open problems have moved to harder physics (poroelastic, heterogeneous, unfitted geometry), to amortisation, and to identifiability theory for the joint configuration-plus-material problem, which Section 8 shows is not a detail.

Flashcards

References

Shield, R. T. (1967). Inverse deformation results in finite elasticity. Zeitschrift für angewandte Mathematik und Physik (ZAMP) 18(4), 490–500. doi:10.1007/BF01601719
Govindjee, S.; Mihalic, P. A. (1996). Computational methods for inverse finite elastostatics. Computer Methods in Applied Mechanics and Engineering 136(1–2), 47–57. doi:10.1016/0045-7825(96)01045-6
Govindjee, Sanjay; Mihalic, Paul A. (1998). Computational methods for inverse deformations in quasi-incompressible finite elasticity. International Journal for Numerical Methods in Engineering 43(5), 821–838. doi:10.1002/(SICI)1097-0207(19981115)43:5<821::AID-NME453>3.0.CO;2-C
Barbone, P. E.; Bamber, J. C. (2002). Quantitative elasticity imaging: what can and cannot be inferred from strain images. Physics in Medicine and Biology 47(12), 2147–2164. doi:10.1088/0031-9155/47/12/310
Oberai, Assad A.; Gokhale, Nachiket H.; Feijóo, Gonzalo R. (2003). Solution of inverse problems in elasticity imaging using the adjoint method. Inverse Problems 19(2), 297–313. doi:10.1088/0266-5611/19/2/304
Barbone, P. E.; Gokhale, N. H. (2004). Elastic modulus imaging: on the uniqueness and nonuniqueness of the elastography inverse problem in two dimensions. Inverse Problems 20(1), 283–296. doi:10.1088/0266-5611/20/1/017
Lu, J.; Zhou, X.; Raghavan, M. L. (2007). Inverse elastostatic stress analysis in pre-deformed biological structures: Demonstration using abdominal aortic aneurysms. Journal of Biomechanics 40(3), 693–696. doi:10.1016/j.jbiomech.2006.01.015
Fachinotti, V. D.; Cardona, A.; Jetteur, P. (2008). Finite element modelling of inverse design problems in large deformations anisotropic hyperelasticity. International Journal for Numerical Methods in Engineering 74(6), 894–910. doi:10.1002/nme.2193
Gee, M. W.; Förster, Ch.; Wall, W. A. (2010). A computational strategy for prestressing patient-specific biomechanical problems under finite deformation. International Journal for Numerical Methods in Biomedical Engineering 26(1), 52–72. doi:10.1002/cnm.1236
Sellier, M. (2011). An iterative method for the inverse elasto-static problem. Journal of Fluids and Structures 27(8), 1461–1470. doi:10.1016/j.jfluidstructs.2011.08.002
Zhao, X.; Raghavan, M. L.; Lu, J. (2011). Identifying heterogeneous anisotropic properties in cerebral aneurysms: a pointwise approach. Biomechanics and Modeling in Mechanobiology 10(2), 177–189. doi:10.1007/s10237-010-0225-7
Doyley, M. M. (2012). Model-based elastography: a survey of approaches to the inverse elasticity problem. Physics in Medicine and Biology 57(3), R35–R73. doi:10.1088/0031-9155/57/3/R35
Rausch, M. K.; Genet, M.; Humphrey, J. D. (2017). An augmented iterative method for identifying a stress-free reference configuration in image-based biomechanical modeling. Journal of Biomechanics 58, 227–231. doi:10.1016/j.jbiomech.2017.04.021
Marx, L.; Niestrawska, J. A.; Gsell, M. A. F.; Caforio, F.; Plank, G.; Augustin, C. M. (2022). Robust and efficient fixed-point algorithm for the inverse elastostatic problem to identify myocardial passive material parameters and the unloaded reference configuration. Journal of Computational Physics 463, 111266. doi:10.1016/j.jcp.2022.111266
Awasthi, C.; Lamperski, A.; Kowalewski, T. M. (2023). Multi-material inverse design of soft deformable bodies via functional optimization. Inverse Problems 39(3), 035006. doi:10.1088/1361-6420/acaa31
Barnafi, N. A.; Petras, A.; Gerardo-Giorda, L. (2024). Fully nonlinear inverse poroelasticity: Stress-free configuration recovery. Computer Methods in Applied Mechanics and Engineering 425, 116960. doi:10.1016/j.cma.2024.116960
Ramella, A.; Lissoni, V.; Bridio, S.; Rodriguez Matas, J. F.; Trimarchi, S.; Grossi, B.; Stefanini, G. G.; Migliavacca, F.; Luraghi, G. (2024). On the necessity to include arterial pre-stress in patient-specific simulations of minimally invasive procedures. Biomechanics and Modeling in Mechanobiology 23(2), 525–537. doi:10.1007/s10237-023-01789-0
Mu, S.; Chan, W. X.; Yap, C. H. (2025). HeartUnloadNet: A Weakly-Supervised Cycle-Consistent Graph Network for Predicting Unloaded Cardiac Geometry from Diastolic States. arXiv preprint arXiv:2507.18677.

← 8. Load & Contact Force Reconstruction 10. Adjoints & Differentiable FEM →