\makeatletter
\@ifundefined{HCode}                  
{\documentclass[screen,CRMECA,Unicode,biblatex]{cedram}
\addbibresource{crmeca20240698.bib}
\newenvironment{noXML}{}{}%
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\sfrac#1#2{{#1}/{#2}}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}\hline}
\def\xtbody{\noalign{\relax}}
\def\tabnote#1{\vskip4pt\parbox{.51\linewidth}{#1}}
\def\xxtabnote#1{\vskip4pt\parbox{.65\linewidth}{#1}}
\def\mathbi#1{\boldsymbol{#1}}
\def\xmathop#1{\mathop{\text{#1}}}
\def\rmDelta{\Delta}
\def\myRef{\ref}
\def\inlinefig#1{\includegraphics{#1}}
\usepackage{etoolbox}
\usepackage{upgreek}
%
\usepackage{graphicx}
\usepackage{amsfonts}
\usepackage{amssymb} 
\usepackage{bm} 
%
\csdef{Seqnsplit}{\\}
\def\botline{\\\hline}
\newcounter{runlevel}
\usepackage{multirow} 
\def\nrow#1{\@tempcnta #1\relax%
\advance\@tempcnta by 1\relax%
\xdef\lenrow{\the\@tempcnta}}
\def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}}
\let\MakeYrStrItalic\relax
\def\jobid{crmeca20240698}
%\graphicspath{{/tmp/\jobid_figs/web/}} 
\graphicspath{{./figures/}} 
%\def\0{\phantom{0}}
\def\rmmu{\upmu}
\def\mn{\phantom{$-$}}
\def\refinput#1{}
\def\back#1{}
\newcommand{\bnabla}{\boldsymbol{\nabla}}
\newcommand{\bgamma}{\boldsymbol{\gamma}}
\newcommand{\bdelta}{\boldsymbol{\delta}}
\newenvironment{inftab}{}{}
\def\Custref#1#2{\hyperlink{#2}{#1}}
\def\og{\guillemotleft}
\def\fg{\guillemotright}
\makeatletter
\g@addto@macro{\UrlBreaks}{\UrlOrds}
\gappto{\UrlBreaks}{\UrlOrds}
\makeatother
\input pdf-trans
\newbox\qbox
\def\usecolor#1{\csname\string\color@#1\endcsname\space}
\newcommand\mycolsplit[2]{\colorlet{tmpcolor}{#2}\edef\tmp{\usecolor{tmpcolor}}%
  \def\tmpB{}\expandafter\colsplithelp\tmp\relax%
  \ifnum0=#1\relax\edef\fillcol{\tmpB}\else\edef\bordercol{\tmpC}\fi}
\def\colsplithelp#1#2 #3\relax{%
  \edef\tmpB{\tmpB#1#2 }%
  \ifnum `#1>`9\relax\def\tmpC{#3}\else\colsplithelp#3\relax\fi
}
\newcommand\outline[1]{\leavevmode%
  \def\maltext{$#1$}%
  \setbox\qbox=\hbox{\maltext}%
  \boxgs{Q q 2 Tr \bbthickness\space w \fillcol\space \bordercol\space}{}%
  \copy\qbox%
}
\newcommand\bordercolor[1]{\mycolsplit{1}{#1}}
\newcommand\fillcolor[1]{\mycolsplit{0}{#1}}
\bordercolor{black}
\fillcolor{white}
\newcommand\bbthickness{.45}
\newcommand\mymathbb[1]{%\def\mydelim{$}
\outline{#1}} %$
\newcommand{\bbPi}{\mymathbb{\Pi}}
\DOI{10.5802/crmeca.279}
\datereceived{2024-08-26}
\daterevised{2024-12-05}
\dateaccepted{2024-12-06}
\ItHasTeXPublished
}
{
\documentclass[crmeca]{article}
\usepackage[T1]{fontenc} 
\usepackage[T2A,T1]{fontenc}
\usepackage{upgreek}
\def\CDRdoi{10.5802/crmeca.279}
\def\newline{\unskip\break}
\def\bbPi{\mathbb{\Pi}}
\def\sfrac#1#2{{#1}/{#2}}
\def\xxtabnote#1{\tabnote{#1}}
\def\xmathop#1{\mathop{{#1}}}
\def\xtbody{\tbody}
\let\printbibliography\relax
\def\hypertarget#1#2{\label{#1}}
}

\usepackage{moreverb}

\newcommand{\MEDISIM}{M$\mathsf{\Xi}$DISIM\xspace}
\newcommand{\newefvec}[1]{\underline{\mathbb{#1}}}
\newcommand{\newefmatrix}[1]{\underline{\underline{\mathbb{#1}}}}
\newcommand{\neweftensorthree}[1]{\underline{\underline{\underline{\mathbb{#1}}}}}
\newcommand{\newefmatrixtilde}[1]{\underline{\underline{\tilde{\mathbb{#1}}}}}

\begin{document}

\begin{noXML}

\CDRsetmeta{articletype}{research-article}

\title{Finite strain formulation of the discrete equilibrium gap
principle: application to direct parameter estimation from large
full-fields measurements}

\alttitle{Formulation en transformation finie du principe d'\'ecart
d'\'equilibre discret : application \`a l'estimation directe de
param\`etres \`a partir de mesures de champs}

\author{\firstname{Alice} \lastname{Peyraut}\CDRorcid{0009-0000-9445-8060}\IsCorresp}
\address{Solid Mechanics Laboratory, \'Ecole Polytechnique/IPP/CNRS, France}
\address{MEDISIM Team, INRIA, France}
\email[A. Peyraut]{alice.peyraut@polytechnique.edu}

\author{\firstname{Martin} \lastname{Genet}\CDRorcid{0000-0003-2204-201X}}
\addressSameAs{1}{Solid Mechanics Laboratory, \'Ecole Polytechnique/IPP/CNRS, France}
\addressSameAs{2}{MEDISIM Team, INRIA, France}
                                                 
\begin{abstract}
The Equilibrium Gap Method (EGM) is a direct model parameter
identification method, i.e., that does not require any resolution of
the model. It has been extensively studied in the context of small
strains but not thoroughly investigated for large strains. In this
article, we propose a novel formulation of the EGM, valid in large
strains, and applicable to both boundary and body forces, when
full-field measurements are available. {Our formulation is based on
a recently proposed continuous formulation and consistent
discretization of the equilibrium gap principle.} Additionally, we
developed an estimation pipeline to quantify the robustness of our new
EGM formulation to noise, and we compared its performance to other
classical estimation methods, namely the Finite Element Model Updating
(FEMU) method and the Virtual Fields Method (VFM). Our robustness
quantification pipeline involves generating synthetic data {from a
reference model} through two methods: {by adding noise to the
reference displacement}, or by generating noisy images and
{performing} motion tracking with the Equilibrium Gap
{principle} used as mechanical regularization. While the quality of
estimation using our new EGM formulation is poor with the first data
generation method, it improves drastically with the second method.
Since the second method of synthetic data generation closely mimics
experimental processes, the EGM, when combined with motion tracking
with Equilibrium Gap regularization, demonstrates reasonable noise
robustness. Thus, it is a promising option for direct parameter
estimation from full-field measurements.
\end{abstract}

\begin{altabstract}
La m\'ethode de l'\'ecart d'\'equilibre (EGM) est une m\'ethode directe
d'identification de param\`etres de mod\`eles, c'est-\`a-dire qu'elle
ne n\'ecessite aucune r\'esolution du mod\`ele. Elle a \'et\'e
largement \'etudi\'ee dans le contexte des petites d\'eformations, mais
n'a pas fait l'objet d'un examen approfondi pour les grandes
d\'eformations. Dans cet article, nous proposons une nouvelle
formulation de l'EGM, valable pour les grandes d\'eformations et
applicable \`a la fois aux forces surfaciques et aux forces volumiques,
lorsque des mesures plein champ sont disponibles. Notre approche est
bas\'ee sur une formulation continue et une discr\'etisation
consistante du principe de l'\'ecart d'\'equilibre r\'ecemment
propos\'ees. En outre, nous avons d\'evelopp\'e un pipeline pour
quantifier la robustesse de notre nouvelle formulation EGM au bruit, et
nous avons compar\'e ses performances \`a d'autres m\'ethodes
d'estimation classiques, \`a savoir la m\'ethode \og Finite Element
Method Updating \fg{} (FEMU) et la m\'ethode des champs virtuels (VFM).
Notre pipeline de quantification de la robustesse implique la
g\'en\'eration de donn\'ees synth\'etiques \`a partir d'un mod\`ele de
r\'ef\'erence via deux m\'ethodes potentielles : soit en ajoutant du
bruit au d\'eplacement de r\'ef\'erence, soit en g\'en\'erant des
images bruit\'ees et en effectuant un suivi de mouvement avec le
principe de l'\'ecart d'\'equilibre utilis\'e comme r\'egularisation
m\'ecanique. Alors que la qualit\'e de l'estimation utilisant notre
nouvelle formulation EGM est faible avec la premi\`ere m\'ethode de
g\'en\'eration de donn\'ees, elle s'am\'eliore consid\'erablement avec
la seconde m\'ethode. \'Etant donn\'e que la deuxi\`eme m\'ethode de
g\'en\'eration de donn\'ees synth\'etiques reproduit fid\`element les
processus exp\'erimentaux, l'EGM, lorsqu'elle est associ\'ee au suivi
des mouvements avec r\'egularisation de l'\'ecart d'\'equilibre, fait
preuve d'une robustesse raisonnable face au bruit. Il s'agit donc d'une
option prometteuse pour l'estimation directe de param\`etres \`a
partir de mesures plein champ.
\end{altabstract}

\keywords{\kwd{Equilibrium Gap Method}\kwd{Inverse problems}\kwd{Large
strains}\kwd{Parameter identification}\kwd{Uncertainty quantification}}

\altkeywords{\kwd{M\'ethode de l'\'ecart \`a
l'\'equilibre}\kwd{Probl\`emes inverses}\kwd{Grandes
transformations}\kwd{Identification}\kwd{Quantification d'incertitude}}

\thanks{French National Research Agency (ANR) (Grant no: ANR-19-CE45-0007)}

\maketitle

\end{noXML}

%\tableofcontents

\section{Introduction}

Inverse problems represent a critical issue in many mechanical
problems. This is especially true in biomechanics for example, where
{in-vivo} testing is rarely feasible~\cite{patte2022estimation}.
As a result, measurements derived from clinical imaging techniques,
such as CT-scans or MRI, can be used to identify parameters such as
tissue stiffness~\cite{laville2023comparison}. Thus, inverse problems
are a valuable alternative to experimental methods for parameters
identification. 

Different estimation methods have been developed for inverse problems,
which all present advantages and drawbacks in terms of accuracy,
robustness to noise and model errors, and computation time. In the
scope of this study, we will assume to have access to full-field
measurements as input for the different estimation techniques. A review
of the most common identification methods using full-field measurements
can be found in~\cite{avril2008overview} and~\cite{roux_optimal_2020}.
Note however that methods such as the Finite Element Model Updating
(FEMU) method, do not require full-field measurements to perform well,
and can be applied to more limited data~\cite{collins1974statistical}.

In this article we propose a new formulation of the Equilibrium Gap
Method (EGM), a direct method---i.e., which does not require direct
model resolution---, valid in large strains, with a consistent FE
discretization. We then quantify the estimation quality of the EGM
applied to large motion measurements, for data presenting high noise
levels, and for boundary force and body force tests. If some estimation
methods, such as the Virtual Fields Method (VFM) or FEMU, have been
extensively studied in small strains~\cite{bonnet2005inverse}, as well
as in large strains~\cite{wittek2016finite, deng_fenics_2023}, the EGM
has only been thoroughly studied in small
strains~\cite{claire2004finite,roux2008digital}, but not in large
strains.

If \cite{boddapati2023single} already proposed a generalization of the
VFM in large strains including all discrete virtual fields, which then
represents an EGM formulation, we propose in this article an
alternative formulation of the method, based on the novel Equilibrium
Gap formulation of~\cite{genet_finite_2023}. Our new formulation
presents the advantage to be fully consistent: it converges when the FE
mesh is refined.

Additionally, the EGM has mostly been studied for boundary forces,
while considering that no body force was
applied~\cite{claire2004finite, crouzeix2009orthotropic}. Our new
formulation allows to take into account full body fields. In order to
develop a general formulation of the EGM, we therefore choose to take
into account both boundary forces and body forces.

To quantify the performance of our new formulation, we study its noise
robustness, and compare it to two other classical methods: the FEMU
method~\cite{kavanagh1971finite, pagnacco2013parameter}, which is a
very robust but indirect method, i.e.{,} a method requiring a model
computation at each iteration, and the Virtual Fields Method
(VFM)~\cite{rossi_effect_2015}, which displays a lower noise
robustness~\cite{avril_sensitivity_2004}, but is potentially a direct,
and hence less expensive method. Two families of virtual fields will be
investigated for the VFM. The first family of virtual fields is
generated with plane waves{, which are not optimal for the simple
problem considered here, but allow to study the impact of virtual field
frequency on estimation performance}. In this case the VFM is a direct
method. The second family of virtual fields is generated based on the
method proposed by~\cite{deng_fenics_2023}{, which was designed to
apply the VFM to complex geometries}. In this case, model computations
are required, and the VFM is not a direct method anymore.

In this article, we therefore propose a performance quantification
method, based on synthetic data generation. Two methods of data
generation, based on generating noisy displacement fields and noisy
images, are studied. We then apply the pipeline to our new formulation
of the EGM, to the FEMU method and to the VFM, and study the impact of
noise on the accuracy of the estimation.

\section{Methods}\label{sec:methods}
\subsection{Estimation framework}\label{sec:estimationframework}

The focus of this article is the introduction of a new formulation of
the Equilibrium Gap Method (EGM) for model parameter identification,
and to compare quantitatively its performances to other usual
estimation method. In order to formally introduce the formulation of
the different methods, we start by introducing classical variables and
operators. 

\paragraph*{Kinematics}

The mapping from a material point at a position $\underline{X}$ in the
reference configuration $\Omega$ (with $|\Omega|$ representing the
reference volume) to a position $\underline{x}$ in a deformed
configuration $\omega$ is represented by the function
$\underline{\chi}$, such that
{\begin{equation}
\underline{\chi} : \underline{X} \mapsto \underline{x}.
\end{equation}}\unskip
The displacement field from the reference to the deformed configuration
$\underline{U}$ is then expressed by
{\begin{equation}
\underline{U} := \underline{x} - \underline{X}.
\end{equation}}\unskip
In order to define the deformation gradient tensor, we classically
derive it from the mapping:
{\begin{equation}
\underline{\underline{F}} := \underline{\underline{\nabla}} \,
\underline{\chi} = \underline{\underline{1}} +
\underline{\underline{\nabla}} \, \underline{U}.
\end{equation}}\unskip
The volume change $J$ is in turn defined through its relationship to the
deformation gradient tensor:
{\begin{equation}\label{eq:jacobian}
J := \mathrm{det}(\underline{\underline{F}}).
\end{equation}}\unskip
We can define two additional quantities also related to the deformation
tensor: the right Cauchy-Green tensor, expressed as
{\begin{equation}
\underline{\underline{C}} := 
\underline{\underline{F}}^{\mathrm{T}} \cdot
\underline{\underline{F}},
\end{equation}}\unskip
and the Green--Lagrange strain tensor, expressed as
{\begin{equation}
\underline{\underline{E}} := \tfrac{1}{2} (\underline{\underline{C}} -
\underline{\underline{1}}).
\end{equation}}\unskip
We can also define the linearized strain tensor
$\underline{\underline{\epsilon}}$ as
{\begin{equation}
\underline{\underline{\epsilon}} := \tfrac{1}{2}
(\underline{\underline{\nabla}} \, \underline{U} +
\underline{\underline{\nabla}} \, \underline{U} ^{\mathrm{T}}).
\end{equation}}\unskip

Throughout this article, we distinguish quantities defined in the
deformed or the reference configuration. In particular, $\Gamma$ and
$\gamma_t$ refer to the reference and deformed boundary where a load is
applied. Similarly, we call $\Gamma_0$ and $\gamma_{t0}$ the reference
and deformed boundaries where a displacement is prescribed. 

Additionally, the outer normal in the reference configuration is
denoted $\underline{N}$, while the outer normal in the deformed
configuration is denoted $\underline{n}$, which can be pulled back to
the reference configuration:
{\begin{equation}
\underline{n} = \frac{\underline{\underline{F}}^{-T}\,\underline{N}}
{\Vert \underline{\underline{F}}^{-T} \underline{N}\Vert}.
\end{equation}}\unskip

\paragraph*{Stresses}

Let us also introduce various stress tensors. Both first and second
Piola--Kirchhoff stress tensors, denoted as $\underline{\underline{P}}$
and $\underline{\underline{\Sigma}}$, can be related to the Cauchy
stress tensor $\underline{\underline{\sigma}}$ using the operators
previously introduced, omitting the mapping between the reference and
the deformed configuration to simplify notations:
{\begin{equation}
\underline{\underline{\Sigma}} := J 
\underline{\underline{F}}^{-1} \,
\underline{\underline{\sigma}} \, 
\underline{\underline{F}}^{-T},
\end{equation}}\unskip
and
{\begin{equation}\label{eq:sigma_2nd_principle}
\underline{\underline{P}} := J \, 
\underline{\underline{\sigma}} \, 
\underline{\underline{F}}^{-T}.
\end{equation}}\unskip

\paragraph*{Loading}

Lastly, let us introduce two vector fields, which both represent forces
applied to the studied system. These two vector fields will be
investigated in this article, and represent a particular case, which
can however straightforwardly be generalized. The first vector field
represents a boundary force applied to a given boundary of the system,
and is denoted as $\underline{T}$ in the reference configuration and as
$\underline{t}$ in the deformed configuration. Note that
$\underline{T}=J||\underline{\underline{F}}^{-T} \, \underline{N}||
\underline{t}$. The second vector field, denoted as $\underline{B}$ in
the reference configuration and as $\underline{b}$ in the deformed
configuration, represents a body force applied to the system. Note that
$\underline{B}=J\underline{b}$. We can also define the vector field
$\underline{U}_i$, which represent imposed displacements in the
reference configuration.

\paragraph*{Equilibrium}

The equilibrium, which derives from the principle of conservation of
momentum, can be expressed in terms of any stress {measurement}
through the balance of internal and external virtual works:
{\begin{equation}\label{eq:vritual_work}
\underline{\underline{\sigma}}\,/\,\underline{\underline{\Sigma}}\,/\, 
\underline{\underline{P}}\quad | \quad 
W_{\mathrm{int}}\left(\underline{\underline{\sigma}}\,/\,
\underline{\underline{\Sigma}}\,/\,\underline{\underline{P}}; 
\underline{U}^*\right) =  W_{\mathrm{ext}}\left(\underline{U}^*\right)\,
\forall\,\underline{U}^*,
\end{equation}}\unskip
where $\underline{U}^*$ represents any kinematically admissible to zero
virtual field. Let us recall that when mentioning in this article that
a virtual field $\underline{U}^*$ is kinematically admissible to zero,
it implies that it is smooth enough on $\Omega$, and that is satisfies
homogeneous Dirichlet conditions on the boundaries where a displacement
is prescribed, i.e., $\underline{U}^*=\underline{0}$ on
$\underline{\Gamma}_0$.

Classically, $W_{\mathrm{int}}$ is expressed as: 
{\begin{eqnarray}
W_{\mathrm{int}}\left(\underline{\underline{\sigma}}\,/\, 
\underline{\underline{\Sigma}}\,/\,\underline{\underline{P}}; 
\underline{U}^*\right) &:=& 
\displaystyle \int_{\omega}^{} \underline{\underline{\sigma}}: 
\underline{\underline{\epsilon}} (\underline{U}^*)\,\mathrm{d}\omega \nonumber \\
&=& \displaystyle \int_{\Omega}^{} \underline{\underline{\Sigma}}:
d_{U} \underline{\underline{E}} \cdot \underline{U}^* \,\mathrm{d}\Omega \nonumber \\
&=& \displaystyle \int_{\Omega}^{} \underline{\underline{P}}: 
\underline{\underline{\nabla}}\,\underline{U}^* \,\mathrm{d}\Omega,
\end{eqnarray}}\unskip
where $d_{U} \underline{\underline{E}} \cdot \underline{U}^* =
(\underline{\underline{F}}^{\mathrm{T}} \cdot \underline{\underline{\nabla}} \,
\underline{U}^*)_{\mathrm{sym}}$ is the first differential of the
Green--Lagrange strain tensor.

Additionally, $W_{\mathrm{ext}}$ can be expressed as:
{\begin{eqnarray}\label{eq:W_ext_boundaryforce}
W_{\mathrm{ext}}\left(\underline{U}^*\right)
&:=& \displaystyle \int_{\gamma_t}^{} \underline{t} \cdot 
\underline{U}^* \,\mathrm{d}\gamma + \displaystyle 
\int_{\omega}^{} \underline{b} \cdot \underline{U}^* \,\mathrm{d}\omega 
\nonumber\\
&=& \displaystyle \int_{\Gamma}^{} \underline{T} \cdot 
\underline{U}^* \,\mathrm{d}\Gamma + 
\displaystyle \int_{\Omega}^{} \underline{B} \cdot 
\underline{U}^* \,\mathrm{d}\Omega.
\end{eqnarray}}\unskip

\paragraph*{Behavior}

According to the second principle of thermodynamics, the second
Piola--Kirchhoff stress tensor $\underline{\underline{\Sigma}}$, which
depends on the state variable $\underline{\underline{E}}$, derives from
the material free energy potential $\Psi$:
{\begin{equation}
\underline{\underline{\Sigma}} = 
\frac{\partial \Psi}{\partial \underline{\underline{E}}}.
\end{equation}}\unskip

\paragraph*{Problem formulation}

Based on the equilibrium and behavior defined previously, the problem
consists in solving the following equation:
{\begin{equation}\label{eq:vritual_work-pb}
\underline{U}\quad | \quad 
W_{\mathrm{int}} \left(\underline{\underline{\Sigma}}(\underline{U}) ; \underline{U}^*\right) =
W_{\mathrm{ext}} \left(\underline{U} ; \underline{U}^*\right) 
\, \forall \, \underline{U}^*.
\end{equation}}\unskip

\paragraph*{Linearization}

In this study, we will also perform analyses in the small strains
setting, aiming to simplify the discussion. When considering linearized
quantities, the internal virtual work is expressed as
{\begin{equation}
W_{\mathrm{int}}(\underline{U}, \underline{U}^*) \approx 
\displaystyle \int_{\Omega}^{}  \underline{\underline{\sigma}}
(\underline{U}) : \underline{\underline{\epsilon}}
(\underline{U}^* )\,\mathrm{d}\Omega,
\end{equation}}\unskip
with
$\underline{\underline{\sigma}}(\underline{U}) = 
\underline{\underline{\underline{\underline{K}}}}\,: 
\,\underline{\underline{\epsilon}}(\underline{U})$, 
where $\underline{\underline{\underline{\underline{K}}}}$ is the
material stiffness tensor, i.e., such that
{\begin{equation}
\Psi = \tfrac{1}{2} \underline{\underline{\epsilon}} : 
\underline{\underline{\underline{\underline{K}}}} : 
\underline{\underline{\epsilon}}.
\end{equation}}\unskip
Similarly, the external virtual work can be expressed as:
{\begin{equation}
W_{\mathrm{ext}}(\underline{U}^*) \approx 
\displaystyle \int_{\Omega}^{} \underline{b} \cdot 
\underline{U}^* \,\mathrm{d}\Omega + \displaystyle \int_{\Gamma}^{} 
\underline{t} \cdot \underline{U}^* \,\mathrm{d}\Gamma.
\end{equation}}\unskip

\paragraph*{Finite elements resolution}

Using finite elements discretization, the linearized direct problem
reduces to
{\begin{equation}
\newefvec{U} \, \quad | \quad \,  \newefmatrix{K} \cdot 
\newefvec{U} = \newefvec{F},
\end{equation}}\unskip
where $\newefmatrix{K}$ is the system stiffness matrix, defined as 
{\begin{equation}
\newefmatrix{K} := \displaystyle \int_{\Omega}^{} 
\neweftensorthree{Q} : 
\underline{\underline{\underline{\underline{K}}}} : 
\neweftensorthree{Q}^{\mathrm{T}}, 
\end{equation}}\unskip
with $\neweftensorthree{Q} $ the array of shape functions symmetric
gradients. $\newefvec{F}$ represents the system force vector, which can
include both body and boundary forces, and is expressed as
{\begin{equation}
\newefvec{F}:= \displaystyle \int_{\Omega}^{} \newefmatrix{N} \cdot
\underline{b}\,\mathrm{d} \Omega +\displaystyle \int_{\Gamma}^{}
\newefmatrix{N} \cdot \underline{t} \,\mathrm{d} \Gamma,
\end{equation}}\unskip
with $\newefmatrix{N}$ the array of the shape functions. Additionally,
$\newefvec{U}$ is the nodal displacement vector, defined through the
approximation of the displacement field $\underline{U}$, such that
{\begin{equation}
\underline{U} = {}^{\mathrm t}{}{\newefmatrix{N}} 
\cdot \newefvec{U}.
\end{equation}}\unskip

We can also introduce a diagonal matrix $\newefmatrix{D}$, which
selects the degrees of freedom where the force is imposed. It is
particularly useful for distinguishing cases where only boundary forces
or only body forces are applied. For example, if only boundary forces
are applied, $\newefmatrix{D}$ is equal to 1 on the boundary degrees of
freedom and 0 everywhere else. 

\paragraph*{Estimation problem formulation}

In general, not all model parameters can necessarily be identified
based on the available data. As a result, for inverse problems, some
parameters, which are called $\underline{\theta}$, are identified,
while other parameters, called $\underline{\eta}$, are chosen at
reference values during the estimation process. Note that estimated
parameters can be material or loading parameters.

The estimation process consists in finding the optimized parameters
$\underline{\theta}^{\mathrm{sol}}$, which minimize a cost function
${L}$. The expression of the cost function depends on the chosen
estimation method. In general, the estimation problem can therefore be
written as:
{\begin{equation}
\underline{\theta}^{\mathrm{sol}} := 
\underset{\underline{\theta}}{\mathrm{argmin}}
\{L(\underline{\theta}, \underline{\eta})\}.
\end{equation}}\unskip

This cost function implicitly includes a displacement field
$\underline{U}_{\mathrm{meas}}$, extracted from experimental
{measurements}. Its finite element discretization is denoted as
$\newefvec{U}_{\mathrm{meas}}$.

\subsection{Recall of existing identification methods}

This section recalls the principle and formulation of three classical
identification methods, whose robustness to noise and model errors is
investigated in this article. 

\subsubsection{The Finite Element Model Updating (FEMU) method}
\paragraph*{Principle}

The ``displacement'' identification method, a.k.a. FEMU (Finite Element
Model Updating), first introduced in~\cite{kavanagh1971finite}, has
been widely used in estimation problems, and is the most robust
identification method~\cite{roux_optimal_2020, pierron2023material}.
The FEMU method estimation relies on a very intuitive process, which
consists in finding the parameters minimizing the distance between the
displacement computed with the model, which depends on the different
parameters, and the measured displacement, which is typically extracted
from experimental images~\cite{avril2008overview}. This method
therefore requires a finite element computation at each iteration; as
such, FEMU is also one of the most expensive identification methods.
Note that FEMU is the only method presented in this article, which does
not require full-field measurements for the estimation. It can also
perform well when only partial data is available.

\paragraph*{Cost function}

The identification problem can be formalized as
follows~\cite{roux_optimal_2020}:
{\begin{eqnarray}  
{L}^{\mathrm{FEMU}}(\underline{\theta}, \underline{\eta}) 
&:=& \tfrac{1}{2} 
\Vert \underline{U}(\underline{\theta}, 
\underline{\eta}) - \underline{U}_{\mathrm{meas}}\Vert^2_{L_2(\Omega)} 
\nonumber\\ 
&\approx& \tfrac{1}{2} 
\Vert \newefvec{U}(\underline{\theta},\underline{\eta}) -
\newefvec{U}_{\mathrm{meas}} \Vert^2_{l_2}. 
\end{eqnarray}}\unskip

\subsubsection{The Virtual Fields Method (VFM)}
\paragraph*{Principle}

The FEMU method, {requiring the resolution of many direct
problems}, is {computationally} expensive. Direct methods,
{which do not require any resolution of the considered model}, are
theoretically cheaper~\cite{roux_optimal_2020}. The Virtual Fields
Method (VFM) is a potentially direct identification method, which has
been widely used in various problems~\cite{grediac_virtual_2006}. It
has notably allowed to successfully identify parameters in small
strains for anisotropic elasticity~\cite{sutton_identification_2008},
plasticity~\cite{avril_stress_2008} or elasto-plastic
behavior~\cite{louedec_identification_2013}. The VFM has also been used
for inverse problems in large strains~\cite{rossi_application_2016}.
Additionally, most identification problems use experimental data, which
are frequently noisy. Being theoretically a direct method, VFM should
be less expensive but also more sensitive to noise than FEMU. The
sensitivity of the VFM to noise has been extensively studied
in~\cite{avril_sensitivity_2004, rossi_effect_2015, avril_stress_2008}.

The VFM is based on the principle of virtual work, presented in
Equation~(\ref{eq:vritual_work}). This method relies on defining at
least as many kinematically admissible to zero virtual displacement
fields as parameters to identify. Each virtual displacement field is
then plugged into Equation~(\ref{eq:vritual_work}) as the virtual field
$\underline{U}^*$, which leads to a system of equations to solve. Note
that when as many virtual fields as parameters to identify are defined,
the resolution may be explicit for simple material laws, e.g., when the
stress tensor depends linearly on the parameters.

\paragraph*{Cost function}

The VFM relies on the principle of virtual work: 
{\begin{equation}\label{eq:equilibrium_vfm}
{L}^{\mathrm{VFM}}(\underline{\theta}, \underline{\eta}) := 
\displaystyle \sum_{i=1}^{m} \left(W_{\mathrm{int}}
(\underline{U}_{\mathrm{meas}}; \underline{U}^*_i; 
\underline{\theta}, \underline{\eta}) - 
W_{\mathrm{ext}}(\underline{U}_{\mathrm{meas}};
\underline{U}^*_i ; \underline{\theta}, \underline{\eta})\right)^2,
\end{equation}}\unskip
where $\underline{U}^*_{i=1, \ldots, m}$ represents the chosen
kinematically admissible to zero virtual fields.

Note that at least as many kinematically admissible to zero virtual
fields as parameters should be chosen to define our equation system
(i.e., $m \geq \mathrm{dim}(\underline{\theta})$). If it is possible to
take into account more virtual fields than parameters ($m >
\mathrm{dim}(\underline{\theta})$), this however leads to an over
determined system, and we will limit our study to the standard setting
$m=\mathrm{dim}(\underline{\theta})$.

\subsubsection{The Equilibrium Gap Method (EGM)}
\paragraph*{Principle}

The Equilibrium Gap Method (EGM) is another potentially direct
identification method that has been widely used for inverse problems in
small strains~\cite{claire2004finite}. As a direct method, EGM should
be computationally more efficient, although less robust to noise and to
model errors than methods such as FEMU~\cite{roux_optimal_2020}. The
EGM can be seen as a generalization of the VFM, where all---in the
discrete sense---virtual fields are
considered~\cite{boddapati2023single}.  The accuracy of the VFM depends
on the choice of the virtual field. If the virtual field is
well-adapted, the estimation is hence much more accurate than for less
judicious choices of virtual fields. The EGM takes into account all
possible virtual fields, the best suited as the worst ones; in average,
the estimation is therefore poorer than for a ``good'' virtual field,
but better than for a ``bad'' choice of virtual field. Note that the
EGM has been mainly used for problems that do not include body forces.

\paragraph*{Cost function}

The EGM consists in finding the parameters that generate stresses that
best respect the equilibrium. In the context of the small perturbation
hypothesis, at the discrete level, the estimation has therefore been
written as~\cite{roux_optimal_2020}:
{\begin{equation} \label{eq:EGM-smallstrains}
{L}^{\mathrm{EGM}, I} (\underline{\theta}, \underline{\eta})
:= \tfrac{1}{2} \Vert \newefmatrix{D}\, \newefmatrix{K}
(\underline{\theta}, \underline{\eta}) \,
\newefvec{U}_{\mathrm{meas}} - \newefmatrix{D}\, 
\newefvec{F}(\underline{\theta}, \underline{\eta})\Vert^2_{l_2},
\end{equation}}\unskip
``$I$'' denoting here the cost function for infinitesimal strains.

\subsubsection{Brief review of other existing methods}

For the sake of simplicity, only the methods described previously will
be used for comparison with our new formulation of the EGM. There are
however many other methods, which are used for inverse problems. In
particular, \cite{ladeveze1994error} proposed the Constitutive Relation
Error (CRE) method, from which other methods, such as the modified
Constitutive Relation Error (mCRE) method \cite{azzouna2015robust,
feissel2007modified} or the Constitutive Equation Gap
Method~\cite{avril2008overview} derive. Additionally, since the EGM
usually displays a very high noise
sensitivity~\cite{roux_optimal_2020}, other formulations, based on the
EGM, were developed to address this problem and decrease the noise
sensitivity of the method. In particular, the Reconditioned Gap Method
was developed in~\cite{roux2008digital, roux_optimal_2020},
and~\cite{bui2004numerical} proposed the Reciprocity Gap Method (RGM)
as estimation method. However, these methods require finite element
computations, and are hence not direct. In the scope of this article,
we will hence only take into account a model updating method (FEMU),
and a potentially direct method (VFM) for comparison with our new EGM
formulation, which is detailed hereafter. {Note that other methods,
such as the Reconditioned Gap Method are also very interesting, as they
are faster than FEMU for example, and more robust than raw EGM. For the
sake of simplicity, however, we only focus on the three specific
methods previously mentionned.}

\subsection{General formulation of the Equilibrium Gap Method}
\paragraph*{Principle}

If the EGM has been widely studied for small strains, there has not
been yet to our knowledge a formulation in large strains, for
hyperelastic problems, including both surface and body 
forces. Boddapati et~al.\ 
\cite{boddapati2023single} has proposed a generalization of VFM in
large strains, by including all discrete virtual fields, which then
represents an EGM formulation. We propose in this article an
alternative formulation of the EGM in large strains, based on a novel
Equilibrium Gap formulation introduced in~\cite{genet_finite_2023} for
motion tracking regularization. Our formulation is fully consistent, in
the sense that it generalizes to any number of degrees of freedom, as
it derives from a proper discretization of a continuous term. Note also
that our formulation linearizes into
Equation~(\ref{eq:EGM-smallstrains}), with an additional mass matrix
that makes the term consistent.

As for the classical formulation of the EGM, our proposed new
expression in large strains consists in finding the parameters
$\underline{\theta}$, which generate stresses that best verify the
equilibrium. Notably, in the general case, the parameters
$\underline{\theta}$ should verify the equilibrium of body and boundary
forces.

\paragraph*{Formulation}

In the scope of our problem, we will only consider body and boundary
forces. Note that imposed displacements are not considered; it would
however be trivial to take them into account~\cite{roux_optimal_2020}.
As a reminder, the strong form of the equilibrium, when applied to the
measured displacement field $U_{\mathrm{meas}}$, is:
{\begin{equation} \label{eq:system-strong-formulation}
\begin{cases}
\underline{\nabla} \cdot \underline{\underline{P}} 
(\underline{U}_{\mathrm{meas}}) + \underline{B} = 0 & \quad \mbox{in } \Omega \\
\underline{\underline{P}}(\underline{U}_{\mathrm{meas}}) \cdot 
\underline{N} = \underline{T} & \quad \mbox{on } \Gamma.
\end{cases}
\end{equation}}\unskip

First, we want to quantify the non-verification of the first Equation
of~(\ref{eq:system-strong-formulation}) by integrating  the norm of its
residual over the mesh. However, the divergence of
$\underline{\underline{P}}$ is not defined for usual finite element
discretizations on edges in 2D and on faces in 3D. Hence,
\cite{genet_finite_2023} proposed an approach to consider the
divergence of a tensor in the finite element formulation. Following
this approach, let us introduce a new vector $\underline{\Pi}$, which
corresponds to the projection of $\underline{\nabla} \cdot
\underline{\underline{P}} + \underline{B}$ onto a space of square
integrable functions $V$. As a result, we can write $\underline{\Pi}$ as:
{\begin{equation}
\underline{\Pi} \quad|\quad \displaystyle \int_{\Omega} \underline{\Pi}
\cdot \underline{\Pi}^* \,\mathrm{d}\Omega = {-} \displaystyle
\int_{\Omega} \left(\underline{\nabla} \cdot \underline{\underline{P}} +
\underline{B}\right) \cdot \underline{\Pi}^* \,\mathrm{d}\Omega \quad 
\forall \underline{\Pi}^* \in V,
\end{equation}}\unskip
where $\underline{\Pi}^*$ is any kinematically admissible field to zero
belonging to $V$. Integrating by parts this expression yields:
{\begin{equation}
\underline{\Pi} \quad|\quad \displaystyle \int_{\Omega} \underline{\Pi}
\cdot \underline{\Pi}^* \,\mathrm{d}\Omega =  \displaystyle \int_{\Omega}
\underline{\underline{P}} : \underline{\underline{\nabla}} \,
\underline{\Pi}^* - \displaystyle \int_{\Omega}^{} \underline{B} \cdot
\underline{\Pi}^* \,\mathrm{d}\Omega \quad \forall \underline{\Pi}^* \in V.
\end{equation}}\unskip

Once discretized, this equation simply leads to the following linear
system:
{\begin{equation}
\newefmatrix{M} \cdot \underline{{\bbPi}} = \newefvec{R},
\end{equation}}\unskip
where $\newefmatrix{M} = \int_{\Omega}^{} \newefmatrix{N}
\cdot \newefmatrix{N}^{\mathrm{T}}$, $\underline{\bbPi}$ is such that 
$\underline{\Pi} =  \newefmatrix{N}^{\mathrm{T}} \cdot \underline{\bbPi}$, and
$\newefvec{R} $ is defined as:
{\begin{equation}
\mathbb{R}_i := \begin{cases}
\displaystyle \int_{\Omega} \underline{\underline{P}} :
\underline{\underline{\nabla}} \, \underline{N}_i \,\mathrm{d}\Omega -
\displaystyle \int_{\Omega} \underline{B} \cdot \underline{N}_i
\,\mathrm{d}\Omega &\quad \mbox{if }i \mbox{ is a body d.o.f} \\
0 &\quad \mbox{if }i \mbox{ is a boundary d.o.f},
\end{cases}
\end{equation}}\unskip
with $\underline{N}_i$ the shape function associated to the $i$th degree
of freedom.

We can therefore define the contribution of the volume force to the
cost function ${L}^{\mathrm{EGM}, F}_{\Omega}$---``$F$'' denoting the
cost function for finite strains---as 
{\begin{eqnarray}
{L}^{\mathrm{EGM}, F}_{\Omega} : &=&
\tfrac{1}{2} \underline{\bbPi}^{\mathrm{T}} \cdot \newefmatrix{M} \cdot \underline{\bbPi} \nonumber\\
&=& \tfrac{1}{2} \newefvec{R}^{\mathrm{T}} \cdot \newefmatrix{M}^{-1} \cdot \newefvec{R}.
\end{eqnarray}}\unskip

Secondly, let us now consider the boundary force. In this case, we
simply need to quantify the non-verification of the second Equation
of~(\ref{eq:system-strong-formulation}):
{\begin{equation}\label{eq:boundaryterm_egm}
{L}^{\mathrm{EGM}, F}_{\Gamma} := \tfrac{1}{2} \displaystyle 
\int_{\Gamma}^{} (\underline{\underline{P}} \cdot \underline{N} - 
\underline{T})^2 \,\mathrm{d}\Gamma.
\end{equation}}\unskip
This term corresponds to the contribution of the boundary force to the
cost function.

\paragraph*{Cost function}

Based on the developments previously presented, the cost function can
be written as:
{\begin{equation}  
{L}^{\mathrm{EGM}, F}(\underline{\theta}, \underline{\eta}) := 
{L}^{\mathrm{EGM}, F}_{\Omega}(\newefvec{U}_{\mathrm{meas}}, 
\underline{\theta}, \underline{\eta}) + {L}^{\mathrm{EGM}, F}_{\Gamma}
(\newefvec{U}_{\mathrm{meas}}, \underline{\theta}, \underline{\eta}).
\end{equation}}\unskip
As a reminder, ``$F$'' denotes the cost function for finite strains. Note
that the cost function ${L}^{\mathrm{EGM}, F}$ depends on
$\underline{\theta}$ and on $\underline{\eta}$ through the material
law---included in $\underline{\underline{P}}$---and through the
loading constants. As mentioned previously, note also that our
formulation linearizes into Equation~(\ref{eq:EGM-smallstrains}), with
an additional mass matrix that makes the term consistent.

In the following sections of this article, any reference to our EGM
will pertain to the new formulation. The FEMU method and the VFM will
be employed as benchmarks to compare estimation performances of our new
formulation.

\subsection{Performance quantification method}

All estimation methods presented previously, excepted FEMU, are only
applicable for problems for which full-field displacement measurements
are available. The estimation relies on finding the parameters
$\underline{\theta}$ minimizing the cost functions defined in the
previous section, which all take in argument (directly, e.g., in FEMU,
or indirectly, e.g., in EGM) the displacement field measured from
experimental images and $\underline{\theta}$. Note that for the VFM,
when as many virtual fields as parameters are used, and for particular
cases (e.g., elastic or simple hyperelastic law), the expression of the
identified parameters can be explicit, and the method hence does not
require any optimization process.

In general, all parameters cannot be identified. As reminded in
Section~\ref{sec:estimationframework}, some parameters,
$\underline{\theta}$, are hence estimated, while all other parameters,
$\underline{\eta}$, are fixed at reference values $\underline{\eta}_0$.
The optimized cost function then allows to access the estimated
parameters $\underline{\theta}^{\mathrm{sol}}$. The estimation process is
described Figure~\ref{fig:estimation-pipeline-a}. However, when using
experimental data, this process does not give any information on the
error committed during the estimation, since the ground-truth
parameters are unknown.

\begin{figure}
\includegraphics{fig01}
\caption{\label{fig:estimation-pipeline-a}Method for identifying the
parameters $\underline{\theta}$. ``tol'' is the tolerance for the
minimization. Note that for the VFM, if there are as many virtual
fields chosen as parameters to identify, and if $W_{\mathrm{int}}$ and
$W_{\mathrm{ext}}$ are linear in $\underline{\theta}$, identifying
$\underline{\theta}$ is equivalent to solving a linear system.}
\end{figure}

In order to quantify the reliability of the estimated parameters, we
would like to investigate the impact of noise on the estimation. To do
so, since the ground-truth values of the estimated parameters are
unknown for real experimental images, we generated synthetic---and
potentially noisy and biased data. Two different methods of data
generation are studied here. The first approach involves computing a
displacement field with our finite element model, using parameters set
to benchmark values, to which noise is possibly added. The second
approach consists in generating images with the model, and adding noise
to these images, following the approach of~\cite{genet_finite_2023}.
Both approaches are detailed hereafter.

\paragraph*{Estimation with noise added to the displacement fields}

In order to create our synthetic data, a first---very
classical~\cite{crouzeix2009orthotropic}---option consists in using a
displacement field computed with our model. Noise is then added to the
displacement field in order to be closer to experimental data, which
always contains noise. To do so, we first choose the parameters to
identify $\underline{\theta}$ at reference values
$\underline{\theta}_0$. We also fix the parameters that are not
identified $\underline{\eta}$ at reference values $\underline{\eta}_0$.
The associated displacement field
$\underline{U}_0=\underline{U}(\underline{\theta}_0,
\underline{\eta}_0)$ is then computed with our finite element model. 

The noise added to the displacement field corresponds to a random
Gaussian field $\newefvec{G}$, defined at the discrete level. We can
also define a continuous random Gaussian field $\underline{g}$. For
each degree of freedom, a random value is generated. This random value
is drawn from a Gaussian distribution with a null mean and a standard
deviation depending on the desired noise level, defined through the
Signal-to-Noise-Ratio (SNR). The noisy displacement field corresponds
to our synthetic {measurement} $\newefvec{U}_{\mathrm{synth}} =
\newefvec{U}_0 + \newefvec{G}$.

Additionally, when using real experimental data, the parameters
$\underline{\eta}$ are fixed at reference values during the estimation
process, but these reference values are very likely different from the
ground-truth values $\underline{\eta}_0$. This is especially true in
biomechanics, for example, where the material parameters are strongly
patient-dependent~\cite{patte2020}. Fixing $\underline{\eta}$ at a
reference value will therefore very likely introduce errors into our
model. 

To study the impact of such model errors on the estimation,  we can use
biased values of $\underline{\eta}$ instead of the ground-truth values
$\underline{\eta}_0$ in the estimation. We propose here to alter
$\underline{\eta}_0$ with a bias $\underline{w}$. To do so, the
synthetic {measurement} is generated with
$\underline{\eta}=\underline{\eta}_0$. We then define a new parameter
$\underline{\eta}_w$, by adding a bias $\underline{w}$ to
$\underline{\eta}_0$. $\underline{\eta}_0$---the parameter used for
building the synthetic measure---is then replaced by the parameter
$\underline{\eta}_w = \underline{\eta}_0+\underline{w}$ during the
estimation. This method allows to investigate the impact of model
errors on the estimation, and is illustrated in
Figure~\ref{fig:estimation-pipeline-b}.

\begin{figure}
\includegraphics{fig02}
\caption{\label{fig:estimation-pipeline-b}Pipeline for generating a
synthetic measure, by adding random Gaussian noise to a reference
displacement field (left) and to a reference image (right). When adding
noise to a displacement field, a displacement field is computed with
the model for reference parameter values. A random Gaussian field is
then added to the displacement field. When adding noise to images, a
reference displacement field is computed with the model by fixing the
parameters at reference values. Associated reference and deformed
images---$I_0$ and $I_t$---are then computed. A random Gaussian field
is then generated to alter the pixel of both images, leading to noisy
images $\tilde{I}_0$ and $\tilde{I}_t$. The synthetic displacement
field is then retrieved with tracking.}
\end{figure}

\paragraph*{Estimation with noise added to the images}

Generating a noisy synthetic {measurement} by adding a random
Gaussian field directly to the displacement field is classical,
although not very realistic. This method indeed leads to a completely
unstructured noise, which might be non-physical. As a result, we
subsequently chose to investigate a second method of synthetic data
generation closer to experimental processes. As such, let us first
describe classical methods used to retrieve experimental displacement
fields from experimental data.

The first step of this method consists in, similarly to the first
method, creat{ing} a synthetic displacement field
$\underline{U}_0=\underline{U}(\underline{\theta}_0,
\underline{\eta}_0)$, computed with the model. Images of the reference,
$I_0$, and of the deformed, $I_t$, configurations are then created. To
create images, different textures can be used, such as no texture,
MR-tagging-like texture, or sinusoidal
patterns~\cite{berberouglu2021silico}. $I_0$ will be defined later. The
deformed image is defined as
{\begin{equation}
I_t(\underline{x}) = I_0(\underline{\chi}^{-1}(\underline{x})),
\end{equation}}\unskip
where $\underline{\chi}$ corresponds to the mapping associated to the
synthetic displacement $\underline{U}_0$.

The next step consists in adding noise to the images: the pixels values
of the images are altered with a random Gaussian noise, with null mean
and a standard deviation chosen depending on the desired noise level.
This leads to noisy images, denoted as $\tilde{I}_0$ and $\tilde{I}_t$.
The synthetic ``measured'' displacement field
$\underline{U}_{\mathrm{synth}}$ is then retrieved by applying a
tracking process, described hereafter.

In general, displacement fields are extracted from images, which
contain noise, through motion tracking. The tracking process applied
between two images can be formulated as finding the smooth mapping
$\underline{\chi}$, or the smooth displacement $\underline{U}$---which
is equivalent---, between the material points of the deformed image to
the reference image, through the minimization
problem~\cite{genet_equilibrated_2018, berberouglu2019validation}:
{\begin{equation}
\underline{U}_{\mathrm{meas}} = \underset{\underline{U}}{\mathrm{argmin}}
\left\{(1-\beta) \Xi^{\mathrm{ima}}(\underline{U}) + \beta
\Xi^{\mathrm{reg}}(\underline{U}) \right\},
\end{equation}}\unskip
with $\Xi^{\mathrm{ima}}$ the image correlation energy,
$\Xi^{\mathrm{reg}}$ the regularization energy, and $\beta$ the
regularization strength. The regularization is added to implement
physical constraints in the resolution, and hence provide structure to
the measurements errors. \cite{leclerc2011voxel}
and~\cite{mendoza2019complete} indeed showed that implementing
regularization allowed to decrease the displacement resolution space,
and to find a solution that follows a mechanical law. This method hence
allows to filter high frequencies, as thoroughly explained in Section
2.2 (correlation procedure) of~\cite{tomivcevc2013mechanics}, which
explains that the Equilibrium Gap Regularization acts as a fourth-order
low-pass filter.

Even though more complex correlation energy terms, designed to decrease
the impact of artifacts during tracking, have been proposed
in~\cite{vskardova2019mechanical, berberouglu2022quantification}, we
will consider here the following correlation energy:
{\begin{equation}
\Xi^{\mathrm{ima}}(\underline{U}) := \tfrac{1}{2} \displaystyle 
\int_{\Omega}^{} (\tilde{I}_t(\underline{X} + \underline{U}(\underline{X}))
- \tilde{I}_0(\underline{X}))^2\,\mathrm{d}\Omega.
\end{equation}}\unskip
This approach is commonly used in image intensity-based
approaches~\cite{roux2008digital, vskardova2019mechanical}, and
presents the advantage of being easily differentiable.
{Additionally, it is particularly adapted to our problem, as this
expression represents the maximal likelihood (in the sense of the least
squares) for Gaussian noise~\cite{roux_optimal_2020}.} This tracking
process is therefore applied, with different regularization levels.
Note that in the scope of this article, we will only use the
equilibrium gap regularization, introduced in large deformation
by~\cite{genet_equilibrated_2018, lee2019validation}, and {later}
on improved in~\cite{genet_finite_2023}. {The formulation is
thoroughly described in the latter article.}

This regularization states that the displacement field obtained should
be as close as possible to a problem with arbitrary boundary traction
values. Creating a synthetic {measurement} with our second method
relies on this tracking process. This synthetic {measurement} can
then be used to study the impact of noise, but also the impact of model
errors---by fixing the parameter $\underline{\eta}$ at an arbitrary
value, which is different from the ground truth value
$\underline{\eta}_0$---on the estimation. {Note  that if the
regularization term for the body force introduced
in~\cite{genet_finite_2023} might seem similar to the term implemented
in the cost function of our new formulation of the EGM, the two are in
fact very different, as they serve different purposes (tracking and
estimation respectively) and take different variables as argument.}

For all methods, the noise added to the displacement field and to the
images is a random Gaussian field. For a given noise level, the
estimation might strongly depend on the random field generated. This is
especially true for high noise levels, which can potentially lead to
either overestimate or underestimate the identified parameter for some
estimation methods. The estimation is therefore performed for many
realizations of a given noise level, until convergence of the estimated
parameters distributions.

\subsection{Illustrative framework}

This article aims at comparing the performance of our new formulation
of the EGM with those of two other classical estimation methods---FEMU
and {VFM}. We choose to perform the estimation on a specific
problem, fairly simple for illustration purposes, described hereafter.

\subsubsection{Geometry}

The estimation is performed in plane strains on a square of
characteristic length of 0.6 mm, meshed with triangular elements. The
estimation is performed on a coarse mesh of 144 elements, as well as on
a finer mesh of 556 elements. {This simple geometry was chosen to
perform a proof of concept of our method and validate its performances;
it is indeed not straightforward that the method will work in large
strains, for high noise levels and fine meshes. To illustrate that the
approach can be extended to more complex geometries, we also propose in
Appendix~\ref{app:appendix3d} the application of the method on a 3D
problem.}

\subsubsection{Constitutive law}

Since we would like to study inverse problems in large hyperelastic
deformations, we choose a Neohookean law to describe the behavior of
our square. In this case, the second Piola--Kirchhoff stress tensor
$\underline{\underline{\Sigma}}$ can be written as
{\begin{equation} \label{eq:Sigma}
\underline{\underline{\Sigma}} = \dfrac{1}{2}
\frac{E \nu}{(1+\nu)(1-2\nu)} (J^2-1) \underline{\underline{C}}^{-1}
+ \frac{E}{2(1+\nu)}(\underline{\underline{1}} - 
\underline{\underline{C}}^{-1}),
\end{equation}}\unskip
where $E$ corresponds to the Young modulus, and $\nu$ corresponds to 
the Poisson ratio.

The inverse problem solved in this article consists in trying to
identify the Young modulus, whose ground truth value is fixed at 1~kPa.
The ground truth value of the Poisson ratio is chosen at 0.3. Note that
for such a simple model, the stress tensor is linear in the Young
modulus.

\subsubsection{Loading}

For the EGM, the estimation with a boundary force has been widely
studied, which is not the case for body forces. Since boundary and body
forces might have a different impact on the estimation---notably, the
cost function differs greatly between the two cases---, we choose to
consider two decoupled problems. In the first problem, the only force
considered is a boundary force, applied to the right boundary of the
square. This boundary force is a compression pressure
$\underline{t}=-\tau \underline{n}$, with $\tau =0.3$ kPa. In the
second problem, only a body force $\underline{b} = \mu
\underline{e}_x$, with $\mu=0.3$~mN/mm$^3$ is applied. In both
problems, the left boundary of the cube is clamped.

\subsubsection{Choice of virtual fields for the VFM}

Finding a kinematically admissible virtual field is rarely
straightforward. Any virtual field chosen should indeed be smooth
enough on the whole domain, but also respect the imposed displacement
fields, i.e., be equal to zero anywhere a displacement is imposed. In
our example, our virtual fields should be null on the left boundary of
the square to respect the imposed boundary conditions. Additionally,
$\underline{U}^*$ should not cancel $W_{\mathrm{ext}}$, in order to
avoid an estimation that would necessarily be wrong. In the scope of
our problem, we studied two different virtual fields.

\paragraph*{Plane waves as virtual fields}

The first virtual field family corresponds to plane waves, written as
{\begin{equation} \label{eq:VF-planewave}
\underline{U}^* = \underline{d} \cdot \sin(\underline{k} \cdot
(\underline{X}-\underline{X}_0)),
\end{equation}}\unskip
where $\underline{d}$ is the plane wave direction, chosen here equal to
the wave propagation direction $\underline{e}$ for the sake of
simplicity. This wave corresponds to a tension-compression wave.
$\underline{k}$ represents the wave vector, and is defined as:
{\begin{equation}
\underline{k} := \frac{2 \uppi}{\Delta} \underline{e},
\end{equation}}\unskip
with $\Delta$ the period of the wave.

This choice was possible given the simplicity of our geometry. It was
indeed easy to find a continuous and differentiable function, null on
the left boundary. {Note that this virtual field was chosen as it
is easy to define, and as it includes a parameter (the frequency),
which can be easily varied, and that we study in this article. However,
a more simple virtual field could be defined, whose perfomances are
studied in Appendix~\ref{app:appendixC}}. It is however much more
difficult to find kinematically admissible displacement fields for more
complex geometries.

\paragraph*{Virtual fields obtained by solving the model}

This is the reason, which motivated \cite{deng_fenics_2023} to propose
a method to build virtual displacement fields, applicable to more
complex problems. This method consists in particular in building a
virtual field from an intermediate position, close to the deformed
position, with parameters close to the deformed configuration as well.
We adapted the code provided in~\cite{deng_fenics_2023_code} to our
problem, and used this virtual displacement field in our study for
comparison with our virtual field built as a plane wave. 

Note that if this method can be applied to any geometry, the
construction of the virtual field requires the computation of a finite
element problem (to compute the intermediate configuration). As a
result, with this choice of virtual field, the VFM is not a direct
method anymore. Increasing the complexity of the geometry does
therefore not always allow to build directly a kinematically admissible
virtual field to zero. {Note also that this virtual field is suited
to investigate more complex problems~\cite{deng_fenics_2023} and hence
presents a complexity which is not required for our simple problem.
However, it remains important to characterize its performance with
respect to simpler virtual fields.}

\paragraph*{Identification in the scope of our hyperelastic problem}

In the scope of the problem solved in this article, the second
Piola--Kirchhoff stress tensor depends linearly on the Young modulus
$E$. For the VFM, the estimation is direct, and consists solely in
solving Equation~(\ref{eq:equilibrium_vfm}): the right-hand side of the
equation is divided by the left-hand side of the equation---as a
reminder, the expression of $\underline{\underline{\Sigma}}$ is
detailed in Equation~(\ref{eq:Sigma}). Both sides are then multiplied
by~$E$.  

In the case of a boundary force problem, the Young modulus is denoted
as $E^{\Gamma}$, and $E^{\Omega}$ in the case a body force problem.
Based on the previous observations, $E^{\Gamma}$ and $E^{\Omega}$ can
directly be expressed as:
{\begin{equation} \label{eq:E_surface/volume}
E^{\Gamma/\Omega} = \frac{W^{\Gamma /
\Omega}_{\mathrm{ext}}}{W_{\mathrm{int}}(E=E_0)},
\end{equation}}\unskip
with $E_0$ the ground-truth value of the Young modulus.

The last step for identifying $E$ is to define a kinematically admissible
virtual displacement field to zero.

\subsubsection{Synthetic {measurements} generation}

For the generation of synthetic {measurements}, the displacement
fields are first computed with the model with the parameters' ground
truth values, i.e., a Young modulus of 1 kPa and a Poisson ratio of
0.3. The noise added to the displacement field and to the images
consists in a Gaussian random field with a null mean, and with SNRs of
${+} \infty$---i.e., no noise---, 10, 5 and 3.3 respectively. Examples of
meshes associated to the noisy {measurements}, obtained by adding
noise to the displacement fields are presented
Figure~\ref{fig:meshes_from_disp}. 

\begin{figure}
\includegraphics{fig03}
\caption{\label{fig:meshes_from_disp}Deformed configuration for
different SNRs---synthetic measure---: ${+} \infty$ (red), 10 (blue), 5
(green), 3.3 (purple), with an element size of 0.1~mm for the boundary
force test, at different steps of the deformation.} 
\end{figure}

As for the generation of noisy measurements, obtained by adding noise
on images, the tracking is done with a regularization using a
Neohookean law, with a Young modulus of 1 kPa and a Poisson ratio of
0.3. Additionally, the images are generated using the following
intensity field for the image $I_0$:
{\begin{equation}
I_0(\underline{X}) = \begin{cases}
0 \quad & \mbox{if } \underline{X} \notin \Omega \\
\sqrt{\sin \left(\frac{\uppi\underline{X}_0}{s}\right) \, 
\sin \left(\frac{\uppi\underline{X}_1}{s}\right)} \quad & \mbox{if }\underline{X} \in \Omega,
\end{cases} 
\end{equation}}\unskip
where $s$ corresponds to the tagging period, chosen as $s=0.1$ in the
scope of our illustrative framework. This expression ensures that only
the body is textured. The images created are presented in
Figure~\ref{fig:images_w_noise}. 

\begin{figure}
{\vspace*{-4pt}}
\includegraphics{fig04}
{\vspace*{-4pt}}
\caption{\label{fig:images_w_noise}Images generated for different SNRs,
with superimposed meshes of the ground-truth solution: ${+}\infty$
(red), 10 (blue), 5 (green), 3.3 (purple), at different steps of the
deformation when the load applied is a boundary force.}
{\vspace*{-4pt}}
\end{figure}

For the same noise level, the deformation depends strongly on the
regularization level applied. Examples of the impact of the
regularization levels on the displacement field extracted from images
are presented Figure~\ref{fig:example-impact-regul}.

\begin{figure}[t!]
{\vspace*{-4pt}}
\includegraphics{fig05}
{\vspace*{-4pt}}
\caption{\label{fig:example-impact-regul}Tracking solutions for
different regularization levels (regularization $=$\ 0 on top,
regularization $=$\ 0.2 at the bottom) for SNRs of 10 (blue), 5
(green), and 3.3 (purple), at different steps of the deformation when
the load applied is a boundary force.}
{\vspace*{-4pt}}
\end{figure}

\subsection{Implementation}

The results generated in this article were obtained through the use of
different libraries and open-access code. In particular, the
implementation of our methods relies heavily on the libraries
``dolfin\textunderscore mech''~\cite{genet_2024_10533833} and
``dolfin\textunderscore warp''~\cite{genet_2023_8010786}, which are
based on FEniCS~\cite{LoggWells2010, AlnaesEtal2015} and
VTK~\cite{vtkBook}.

Additionally, we studied the noise robustness of the VFM, using the
virutal fields proposed by~\cite{deng_fenics_2023}. To do so, we
adapted the code available in~\cite{deng_fenics_2023_code}. Note that
this virtual field requires the computation of the problem and is hence
close to FEMU in terms of computation costs.

The code used to generate the figures presented in the Results section
is freely available. In particular, the different estimation methods
are implemented in the library dolfin\_estim~\cite{peyraut_2024}. The
code to generate the figures is reproduced in the appendix of the
paper, and is currently available interactively at
\url{https://apeyraut.gitlabpages.inria.fr/identification-methods-paper-demos}.

\section{Results}

This section presents the estimation results for the different methods
presented in Section~\ref{sec:methods}. The obtained results allow us
to compare the performance of our new formulation of the EGM with
another direct method, the VFM with a plane wave, and two non-direct
method, the VFM using a virtual field created
from~\cite{deng_fenics_2023}, and the FEMU method. 

Two different problems are considered: a boundary force problem and a
body force problem. Within these two problems, two cases are also
studied: identification using the synthetic {measurements}
generated by altering directly the displacement field computed with the
model, and identification using the synthetic {measurements}
generated by adding noise to the images generated with the model.
Additionally, the estimation is performed for two different mesh sizes:
a mesh with characteristic element size of 0.1~mm, and a mesh with
characteristic element size of 0.05~mm.

Note that for the case where the identification is conducted on the
images altered with noise, different regularization levels were
investigated. However, for the sake of simplicity, only two
regularization levels are presented in this section. The first
regularization level is 0, to illustrate the poor identification
results obtained when no regularization is added. The second case
corresponds to a regularization level of 0.2, for which the
identification is the best for all tests performed---except for the
last figure, which will be detailed later in the article. 
{\vspace*{-4pt}}

\subsection{Estimation performance in the case of boundary forces}
{\vspace*{-4pt}}
\subsubsection{Error distribution on the coarser mesh}
{\vspace*{-4pt}}

This section compares the estimation performance of our proposed EGM
with the performances of the FEMU method and the VFM, when the
synthetic {measurements} are created by adding noise to the
displacements, and by adding noise to the images---for a regularization
level of 0.0 and 0.2. The estimation was conducted using a mesh with a
characteristic element {length} of 0.1~mm. The load considered here
is a boundary force. The {obtained} results are presented in
Figure~\ref{fig:compression_test_coarse}.

\begin{figure}
\includegraphics{fig06}
\caption{\label{fig:compression_test_coarse}Estimation error for the
boundary force test, with a characteristic element size of 0.1~mm. This
figure presents the error distribution on the estimated Young modulus
depending on the noise, for the EGM (in green), the VFM
from~\cite{deng_fenics_2023} (in red), the VFM {with} a plane wave
(in blue), and the FEMU method (in yellow). {For the} top row,
noise {is} directly added to the displacement field, {and for}
the bottom row, the synthetic measurements are generated by adding
noise to the images, with a regularization of 0.0 (left) and 0.2
(right). For all methods, for noise applied {to} the displacements,
{for a} SNR going to ${+}\infty$, the Young modulus is perfectly
identified. When noise increases (smaller SNRs), the error on the
estimated parameter increases as well. FEMU is the most robust method,
while the EGM and the VFM are far less robust. For noise added to the
images, a notable error {is} introduced for a SNR of ${+}\infty$ for
the EGM. When regularization increases, the estimation drastically
improves for all methods.}
\end{figure}

This figure illustrates the error committed on the estimated parameter
depending on the noise. The lines and transparent areas (green for the
EGM, blue for the VFM using a plane wave as the virtual field, red for
the VFM using~\cite{deng_fenics_2023}'s method, and yellow for the FEMU
method) correspond to the averages and the standard deviations of the
distributions respectively.

\paragraph*{Estimation using the synthetic {measurements} generated
by adding noise to the displacements}

When noise is added directly to the synthetic displacement field, the
estimation quality depends strongly on the chosen method. The
same---classical---pattern can however be \mbox{observed} for all methods:
when no noise is applied to a synthetic measurement, the parameter can
be perfectly identified. As noise added to the data increases---i.e.,
as the SNR decreases---, the quality of the estimation decreases: the
average of the distribution grows further from the ground-truth value
of the parameter. Additionally, the uncertainty on the parameter
increases---i.e., the standard deviation is higher.

The FEMU method is the most accurate. It displays high robustness to
noise, with a mean error remaining close to 0 for all noise levels, and
a standard deviation remaining under 3\% for the smaller SNR
studied---i.e., for the highest noise. On the other hand, the EGM
displays {the highest} noise sensitivity. If the Young modulus is
perfectly identified when no noise is applied, the estimation becomes
very poor when noise increases, {with for example and error of
about ${-}$40\%~$\pm$~20\%} for a SNR of 10. Additionally, for all SNRs,
the parameter is underestimated; the smaller the SNR, the higher the
negative error. 

{Additionally, the results obtained for the} VFM illustrate well
the impact of the choice of the virtual field on the estimation: the
distribution is very different whether the virtual field was chosen as
a plane wave or based on the work of~\cite{deng_fenics_2023}. For
example, for all SNRs, the Young modulus tends to be overestimated with
the virtual field from~\cite{deng_fenics_2023}, while it is
underestimated for the plane wave. For both studied virtual fields, the
estimation remains good for SNRs higher or equal to 10. However, for
smaller SNRs, the error on the estimation becomes much higher.

\paragraph*{Estimation using synthetic {measurements} generated by
adding noise to the images}

{The previous results were obtained by altering displacements with
noise, which lead to unstructured non-physical displacement fields. In
order to quantify the estimation performance on more realistic data, we
then chose to generate synthetic measurements by adding noise to
images, and to retrieve the displacement fields with a tracking
process, using two regularization levels: no regularization and a
regularization level of 0.2.} For a null regularization, the behavior
of all methods is very close to what was observed {in the top row
of Figure~\ref{fig:compression_test_coarse}}. The main difference
consists in the introduction of an error of ${-}$25\% for the EGM when no
noise is applied {to} the {measurements}---SNR of ${+} \infty$.
This error is due to errors committed during the tracking process,
{as} the estimation improves drastically {when regularization
increases}.  

For all methods, {for a regularization of 0.2}, the mean error of
the estimation remains close to 0 for all SNRs, and the standard
deviation is much smaller {than for a null regularization}. Note
that the improvement of the estimation is marginal for the FEMU method,
for which the estimation {is robust for every cases studied}. The
estimation with the VFM improves for both virtual fields. The EGM
remains the less accurate method. However, the larger error, around
${-}$10\% $\pm$ 8\%, is obtained for very small SNRs, {and the} error
is much smaller for experimental values of SNRs (around 10 for
CT-scans, for example). The EGM could therefore be used for inverse
problems with our new formulation, since it shows a very good
robustness to noise, as long as the Equilibrium Gap regularization is
also used for the tracking.

\subsubsection{Error distribution on the finer mesh}

The results presented in the previous section were obtained on a rather
coarse mesh, {for which} estimation is known to perform well{.
However, it} can be expected to perform more poorly when the mesh is
refined~\cite{genet_finite_2023}. {W}e therefore {also} chose
to conduct {our study} on a finer mesh, in order to be closer to
applications on real problems. The estimation is performed on a mesh
with a characteristic element size of 0.05~mm. The results are
presented in Figure~\ref{fig:estimated_Young_h=0.05-compression}.

\begin{figure}
\includegraphics{fig07}
\caption{\label{fig:estimated_Young_h=0.05-compression}Estimation error
for the boundary force test, with a characteristic element size of
0.05~mm. This figure presents the error distribution on the estimated
Young modulus depending on the noise, for the EGM (in green), the VFM
from~\cite{deng_fenics_2023} (in red), the VFM {with} a plane wave
(in blue), and the FEMU method (in yellow). Three cases are
investigated: (i)~noise added to the synthetic displacement field (top
row), (ii) noise added to images (middle row), (iii) noise added to
images, {refined} using the multiresolution. If the estimation
quality is very poor {for} the first two case{s}, the
estimation improves drastically with the multiresolution.}
\end{figure}

\paragraph*{Estimation using synthetic {measurements} generated by
adding noise to the displacements}

Similarly to previous observations, when noise is added directly to the
synthetic displacement field, the quality of the estimation depends
strongly on the method chosen. The estimation remains fairly stable
compared to the estimation performed on the coarser mesh{: there
is} no notable difference for the distribution obtained with the FEMU
method{, and the} global behavior of all direct methods remains
stable, {even thoug} the{ir} estimation quality is far poorer.
For {example, for} the VFM, the average identified Young modulus
deviates from its ground truth value more rapidly than with the coarser
mesh with both chosen virtual fields{:} for a SNR of 10, the
average error for the VFM with a plane wave as virtual field is ${-}$6\%
with the coarser mesh, while it is around ${-}$20\% with the finer mesh.
The estimation is also far better with the coarser mesh for the EGM.
The average error for a SNR of 10 is around ${-}$40\% in that case, while
it is almost equal to ${-}90$\% for the finer mesh. If refining the mesh
does not seem to impact the estimation quality with the FEMU method, it
seems to deteriorate greatly the estimation with the direct methods.

\paragraph*{Estimation using synthetic {measurements} generated by
adding noise to the images}

{The estimation was performed for measurements obtained with a null
regularization and a regularization of 0.2.} {For a null
regularization, the} FEMU method and the VFM which uses a plane wave as
virtual field {give} good {results} when no noise is added to
the data. {However, the estimation quality deteriorates rapidly as
noise increases, especially for the VFM. For the FEMU method, the
estimation error also increases with noise, although the standard
deviation of the distribution remains fairly small.} For the VFM
from~\cite{deng_fenics_2023}, the estimation did not converge, even for
a SNR of ${+}\infty$, due {to} tracking errors. For the EGM, when no
noise is applied, an initial error of ${-}$100\% can be observed, due to
tracking errors. This error is much higher than for the coarser mesh
{and deteriorates when the SNR decreases.}

When the regularization increases, the estimation improves for all
direct methods. In particular, for the EGM, the error observed for a
SNR of ${+} \infty$ is corrected{.} The same behavior can be observed
for the VFM from~\cite{deng_fenics_2023}, for which the parameter
becomes identifiable when no noise is added to the images. For the VFM
with a plane wave as virtual field, the estimation improves slightly,
with a smaller standard deviation for {small} SNRs, and an average
error closer to 0. However, the quality of the estimation becomes
poorer when noise increases {for all these cases}. Although the
estimation is  better for direct methods, it becomes far less accurate
for the FEMU method. If the Young modulus is still identifiable when no
noise is added to the images, the estimation becomes less accurate when
noise increases: the average error increases, as well as the standard
deviation. Overall refining the mesh leads to a poorer estimation
{and} does therefore not seem to be a good option for all
estimation methods, with both synthetic data generation methods.

\paragraph*{Estimation based on tracking using the multiresolution
method}

The decrease in accuracy of the estimation when using synthetic
{measurements} with a refined mesh could have been anticipated,
{as using a finer mesh for the tracking process necessarily leads
to a increased noise sensitivity of the solution.} While mechanical
regularization is designed to prevent {such phenomenon}, it may not
be sufficient. Consequently, the multiresolution method, {which
consists in applying the tracking process on successively refined
meshes,} helps to avoid this problem. To perform the tracking on a
finer mesh, the solution at a given frame is initialized by the
solution of the tracking on a coarser mesh, at the same frame. The
results for the estimation using synthetic {measurements} generated
by adding noise to images on a finer mesh obtained with the
multiresolution method, are presented in
Figure~\ref{fig:estimated_Young_h=0.05-compression}, on the bottom row.

The multiresolution improves slightly the results for a null
regularization, compared to the results obtained when performing the
tracking directly with the refined mesh. The FEMU method displays the
best results: its average is close to the ground truth value, and the
standard deviation is small for all SNRs. The initial error for a SNR
of ${+}\infty$ is also smaller for the EGM---even though it is still
high. As noise increases, the average error is also high for this
method, even though the standard deviation remains fairly small.
{Similarly, the estimation is of good quality for the VFM, when no
noise is added to the images, but deteriorates as noise increases.} 

Increasing the regularization, together with the multiresolution,
however improves greatly the estimation, leading to error distributions
close to the results found with the coarser mesh. The FEMU method and
the VFM display the same behavior: for all SNRs, the average of the
distribution is close to the ground-truth value and the standard
deviation is very small. A very good robustness to noise is observed
for all methods. As before, the EGM presents the highest sensitivity to
noise. However, the error committed during the estimation is highest
for a SNR of 3.3 (12\% $\pm$ 8\%) and remains reasonable. As a result,
the EGM can also be considered as a robust estimation method for finer
meshes, should the multiresolution been used in addition to the
Equilibrium Gap regularization.

\subsection{Estimation performance in the case of body forces}

The previous section focused on presenting results allowing the
comparison of our new formulation of the EGM with the VFM and the FEMU
method, for a {boundary force test}. {As the form of the
boundary force term and of the body force term differ greatly in the
loss function, we would also like to evaluate the perfomance of our new
formulation of the EGM for a body force test, which is presented in
this section.}

\subsubsection{Error distribution on the coarser mesh}

This section presents the error distribution on the identified Young
modulus on a mesh with a characteristic element {size} of 0.1~mm,
for the four methods, with a synthetic {measurement} created by
adding noise to the displacements, and with a synthetic
{measurement} created by adding noise to the images---for a
regularization level of 0.0 and 0.2. The only load studied in this
section is a body force. The obtained results are presented in
Figure~\ref{fig:estimated_Young_h=0.1-refined=False-volume}. {Note
that in this section, we compare the results obtained with the boundary
force test and with the body force test, even though the two
displacement fields are different. Indeed, their norm are of the same
order (0.06 for the boundary test and 0.02 for the body force test).
Additionally, the noise added to images is proportional to the
displacement norm. As such, the results obtained with the body and the
boundary force tests can be quantitatively compared.}

\begin{figure}
\includegraphics{fig08}
\caption{\label{fig:estimated_Young_h=0.1-refined=False-volume}Estimation 
error for the body force test, with a characteristic element size of
0.1~mm. This figure presents the error distribution on the estimated
Young modulus depending on the noise, for the EGM (in green), the VFM
from~\cite{deng_fenics_2023} (in red), the VFM {with} a plane wave
(in blue), and the FEMU method (in yellow), for synthetic
{measurements} generated by adding noise to the displacements (top
row), and synthetic {measurements} generated by adding noise to
images (bottom row), for a regularization level of 0.0 (left) and 0.2
(right). The estimation with the body force test is globally better
than with the boundary force test, for all methods except the EGM, when
noise is added to the displacement field. When adding noise to the
images, the results with no regularization are very close to the
results generated by adding noise to the displacements. Adding
regularization improves the estimation for all methods, {and in
particular for the EGM}.} 
\end{figure}

\paragraph*{Estimation using synthetic {measurements} generated by
adding noise to the displacements}

The sensitivity to noise is very similar to what was observed for the
boundary force test, for all methods except for the EGM. {However,
for the VFM}, if the global behavior of the distribution remains
similar, the error is however globally smaller for the body force test
than for the boundary force test (average closer to the ground-truth
value and smaller standard deviation). On the contrary, the estimation
of the Young modulus using the EGM method is poorer with the body force
than for the boundary force test. If the parameter can indeed be
identified perfectly when no noise is applied to the displacement
field, the error becomes very high very fast as noise increases
{and} the Young modulus is systematically underestimated.

\paragraph*{Estimation using the synthetic {measurements} generated
by adding noise to the images}

Two cases are investigated: the estimation is first conducted with a
null regularization, and then with a regularization of 0.2. For a null
regularization, the behavior of all methods is very close to what was
observed {for} synthetic {measurements} built by adding noise
to the displacement field. The main difference with the previous case
consists in the introduction of an error of ${-}$90\% for the EGM
{and} for a SNR of ${+} \infty$ due to tracking errors. Note that
this error is much higher than for the boundary force test.  When the
regularization increases, the estimation improves drastically. As
before, the VFM display{s} a behavior very close to the FEMU
method. As for the EGM, the initial error is corrected, and the method
seems far less sensitive to noise.

\subsubsection{Error distribution on the finer mesh}

As for the boundary force test, we chose to conduct the estimation on a
refined mesh. Knowing that the results would be poor if the estimation
was performed directly on a finer mesh, we chose to only present the
results obtained on a finer mesh generated with the multiresolution
method. The results obtained are presented in
Figure~\ref{fig:estimated_Young_h=0.1-refined=True-volume}. Note that
contrarily to previous cases, the regularization added for this case is
0.8. A regularization of 0.2 indeed did not lead to a good estimation
quality. Hence, the regularization strength had to be increased for
better results.

\begin{figure}
\includegraphics{fig09}
\caption{\label{fig:estimated_Young_h=0.1-refined=True-volume}Error
distribution on the estimated Young modulus depending on the noise for
the body force test, with an element size of 0.5~mm obtained with the
multiresolution method, for the EGM (in green), the VFM
from~\cite{deng_fenics_2023} (in red), the VFM with plane waves (in
green), and the FEMU method (in yellow). The synthetic
{measurement} used for the estimation was generated by adding noise
to images, with a regularization of 0 (left) and 0.8 (right). When no
regularization is added to images, the estimation is very good for the
FEMU method but very poor for all direct methods. The estimation
improves for all direct methods when regularization increases.}
\end{figure}

When no regularization is applied during the tracking process, the
estimation is of very poor quality for all methods {except FEMU},
{but improves when regularization increases}. {In particular,
for a null regularization, the} mean and standard deviation for the EGM
do not even appear on the graph because the error committed during the
estimation process is too high. The estimation quality also improves
drastically {with regularization} for the EGM. For SNRs larger than
5, the standard deviation remains small, and the average error close to
0. For a very small SNR---3.3---, the average error is however high,
and the standard deviation large: the parameter is identified with an
error of ${-}$25\% $\pm$ 25\%. A SNR of 3.3 however corresponds to very
noisy data, and is out of the scope of many mechanical  applications.
Additionally, for small SNRs, the estimation quality for the EGM is
slightly poorer for the body force compared to the boundary force test
(the standard deviation is wider).

\subsection{Impact of the choice of virtual field for the VFM}

The presented results illustrate that the estimation quality with the
Virtual Fields Method (VFM) varies significantly depending on the
virtual field selected for the identification. Our results indeed
indicate that the error distribution of the identified Young modulus
differs depending on whether the virtual field was chosen as a plane
wave or constructed based on the method described
in~\cite{deng_fenics_2023}. For most tests and synthetic
{measurements} considered, both virtual fields indeed tend to have
a different impact on the estimation. For example, if the parameter is
underestimated with the plane wave, it is often overestimated with the
virtual field built from~\cite{deng_fenics_2023}, and conversely. 

This variation in estimation quality is also noticeable with the
virtual field chosen as a plane wave, when modifying the value of
$\Delta$, which is a constant chosen empirically, and represents the
period of the wave, which is directly linked to the plane wave
frequency, as presented in Equation~(\ref{eq:VF-planewave}). The error
distribution on the estimated Young modulus with the virtual field
chosen as a plane wave is presented in Figure~\ref{fig:impact-delta}
for different values of $\Delta$, depending on the SNR. The error
distribution was investigated for {measurements} created by adding
noise to the displacements, for a boundary force test. Since the Young
modulus was not identifiable in that case for SNRs higher than 10, as
presented in Figure~\ref{fig:compression_test_coarse}, the study was
conducted for SNRs ranging from 10 to ${+}\infty$. The values of
$\Delta$ studied range from $0.01l$ to $1000.1l$, where $l$ represents
the characteristic length of the square. Note that the higher $\Delta$,
the smaller the plane wave frequency.

\begin{figure}
\includegraphics{fig10}
\caption{\label{fig:impact-delta}Error distribution for the Young
modulus identified with the VFM, with a plane wave as virtual field,
for different values of the period $\Delta$---ranging from 0.1l to
1000.0l, with l the characteristic length of the cube---, for a
boundary force test and with noise directly applied to the displacement
field. The results are presented for SNRs ranging from 10 to
${+}\infty$. When no noise is applied {to} the synthetic measure,
all values of $\Delta$ lead to a perfect estimation of the Young
modulus. As noise increases, smaller values of $\Delta$ lead very large
estimation error.}
\end{figure}

Figure~\ref{fig:impact-delta} presents the error distribution of the
estimated Young modulus for different values of $\Delta$, {to which
different colors are associated}. The highest value of $\Delta$ appears
in yellow, while the smallest value of $\Delta$ is in blue. The error
distribution is represented with error bars: the mean value of the
distribution is represented as a square, while the bars represent the
mean ${\pm}$ the standard deviation of the distribution. For a SNR
of ${+}\infty$, the Young modulus is perfectly identified for all
values of $\Delta$. For a SNR of 100, i.e., a very small level of
noise, the parameter is still identifiable with all $\Delta$, except
for the smallest value ($\Delta=0.1l$). As noise increases, the smaller
values of $\Delta$ lead to large errors on the identified parameter.
For example, for a SNR of 20, which is still fairly small, the error on
the identified parameter is high for all $\Delta$ greater than 1.7$l$,
i.e.{,} the standard deviation is greater than~25\%. 

However, the estimation quality decreases for high values of $\Delta$
when noise increases. For example, the estimation is better for
$\Delta=2.5l$ than for $\Delta=1000l$ for a SNR of 10, as the average
of the distribution is not centered in zero for the larger value of
$\Delta$. As $\Delta$ increases, the plane wave {indeed} becomes
flatter---i.e., the oscillations of the plane wave are barely
noticeable{--- and} the impact of the noise is much higher.
Choosing the value of $\Delta$ hence consists in finding a good
trade-off between precision and robustness.

Additional studies, which are not directly related to the main topic of
this article but present interesting results regarding the limitations
of our method, are included in the appendices. In particular, the
robustness to model errors of the different methods is studied.
Appendix~\ref{app:appendixA}  presents the impact of model errors on
the Poisson coefficient on the estimated Young modulus, while
Appendix~\ref{app:appendixB} presents the impact of the value of the
volume regularization on the estimation.

\section{Discussion}

The results obtained have highlighted that our new formulation of the
EGM leads to very different results depending on the method of
synthetic data generation, the mesh size, the noise level, and the type
of applied load. Generally, the VFM leads to better results, but seems
to be very dependent on the chosen virtual field. Lastly, the FEMU
method performs better than the two other {investigated} methods.
These results can be explained by observing the analytical forms of the
cost functions of the methods.

\subsection{Estimation performance of the EGM}

For the sake of simplicity, we study here the cost function of the EGM
in small elastic strains. The presented work can easily be extended to
more complex problems. Let us first write the finite element
discretization of the EGM cost function presented in
Equation~(\ref{eq:EGM-smallstrains}):
{\begin{equation} \label{eq:discretization}
{L}^{\mathrm{EGM}} = 
\left(\newefmatrix{D} \, \newefmatrix{K} \, \newefvec{U}_{\mathrm{meas}} - \newefmatrix{D} \, \newefvec{F} \right)^{\mathrm{T}} \newefmatrix{M} 
\left(\newefmatrix{D} \, \newefmatrix{K} \, \newefvec{U}_{\mathrm{meas}} - \newefmatrix{D} \, \newefvec{F}\right).
\end{equation}}\unskip
In linear elasticity, the stiffness matrix $\newefmatrix{K}$ depends
linearly on the Young modulus $E$. Let us thus define a matrix
$\newefmatrixtilde{{K}}$, defined as
{\begin{equation} \label{eq:stiffnessmatrix}
\newefmatrixtilde{{K}} = \frac{1}{E} \, \newefmatrix{K}.
\end{equation} }\unskip
Injecting~(\ref{eq:stiffnessmatrix}) into~(\ref{eq:discretization})
hence yields the following expression of the cost function:
{\begin{equation}  
{L}^{\mathrm{EGM}} =
\left(E\, \newefmatrix{D} \, \newefmatrixtilde{{K}}\, \newefvec{U}_{\mathrm{meas}} - \newefmatrix{D} \, \newefvec{F} \right)^{\mathrm{T}} \newefmatrix{M}
\left(E \newefmatrix{D} \, \newefmatrixtilde{{K}} \, \newefvec{U}_{\mathrm{meas}} - \newefmatrix{D}\, \newefvec{F} \right).
\end{equation}}\unskip

The estimation consists in finding the parameter $E$ minimizing the cost
function, i.e.{,} the parameter $E$ verifying
{\begin{equation} \label{eq:derivthe EGM=0}
\frac{\partial {L}^{\mathrm{EGM}}}{\partial E} = 0,
\end{equation}}\unskip
which leads to the following expression of $E$:
{\begin{equation} \label{eq:expression_E}
E = \frac{1}{2} \frac{\newefvec{F}^{\mathrm{T}}  \, \newefmatrix{D}^{\mathrm{T}} \,
\newefmatrix{M} \, \newefmatrix{D} \, \newefmatrixtilde{{K}} \,
\underline{\newefvec{U}}_{\mathrm{meas}} +
\newefvec{U}_{\mathrm{meas}}^{\mathrm{T}} \, \newefmatrixtilde{{K}}^{\mathrm{T}} \,
\newefmatrix{D}^{\mathrm{T}} \, \newefmatrix{M} \,  \newefmatrix{D} \,
\newefvec{F} }{\newefvec{U}_{\mathrm{meas}}^{\mathrm{T}} \, 
\newefmatrixtilde{{K}}^{\mathrm{T}} \, \newefmatrix{D}^{\mathrm{T}} \, \newefmatrix{M} \,
\newefmatrix{D} \, \newefmatrixtilde{{K}} \,
\newefvec{U}_{\mathrm{meas}}}.
\end{equation}}\unskip

Let us first investigate the case where the synthetic {measurement}
is generated by adding noise to the displacements.
$\newefvec{U}_{\mathrm{meas}}$ therefore corresponds to the
approximated measured---and noisy---displacement field and can be
expressed as
{\begin{equation} \label{eq:u_with_noise} 
\newefvec{U}_{\mathrm{meas}} = \newefvec{U}_0 + \newefvec{G},
\end{equation}}\unskip
where $\newefvec{G}$ represents the noise added to the approximated
reference---i.e., associated to $\underline{\theta}_0$ and
$\underline{\eta}_0$---displacement field $\newefvec{U}_0$. Injecting
Equation~(\ref{eq:u_with_noise}) into Equation~(\ref{eq:expression_E})
therefore gives the following expression for the Young modulus:
{\begin{equation}\label{eq:expression_E_with_noise}
E = \frac{1}{2}
\frac{\newefvec{F}^{\mathrm{T}} \,
\newefmatrix{D}^{\mathrm{T}} \,
\newefmatrix{M}\,
\newefmatrix{D}\,
\newefmatrixtilde{{K}}\,
(\newefvec{U}_0 + \newefvec{G}) -
(\newefvec{U}_0 + \newefvec{G})^{\mathrm{T}}\,
\newefmatrixtilde{{K}}^{\mathrm{T}}\, 
\newefmatrix{D}^{\mathrm{T}}\,
\newefmatrix{M}\, 
\newefmatrix{D}\, 
\newefvec{F}}
{(\newefvec{U}_0 + \newefvec{G})^{\mathrm{T}}\,
\newefmatrixtilde{{K}}^{\mathrm{T}}\, 
\newefmatrix{D}^{\mathrm{T}}\, 
\newefmatrix{M}\,
\newefmatrix{D}\,
\newefmatrixtilde{{K}}\, 
(\newefvec{U}_0 + \newefvec{G})}.
\end{equation}}\unskip

According to Equation~(\ref{eq:expression_E_with_noise}), the Young
modulus depends linearly on the noise on the numerator and
quadratically on the denominator. As a result, the more the noise
increases, the more the Young modulus is underestimated, {as the}
quadratic dependency of the denominator amplifies the impact of the
high frequencies of noise on the estimation. When no noise is applied
to the displacements, the parameter should moreover be identified
perfectly. This behavior is consistent with the results {obtained.}

{This} is especially true when the noise is added directly to
the displacements{, resulting in highly unstructured measurements.
As a result, this method for data generation, even though a classical
approach for performance quantification in inverse problems, is
non-physical. Using our second generation method, which consists in
adding noise to images, would be a more consistent alternative.} 

When noise is added directly to the images {without regularization,
the generated measurement} is similar to the synthetic measurement
generated by adding noise to displacements: there is no additional term
enforcing a physical behavior during the tracking process, {and}
the noise is hence highly unstructured. Additionally, an initial error
appears for a SNR of ${+} \infty$, due to tracking errors: the EGM
being very sensitive to any small variation of the displacement field,
small tracking errors can indeed lead to significant errors in the
estimation.  This initial error is moreover greater for the body force
test than for the boundary test, {as} only the boundary where the
pressure is applied is considered {for the boundary force test},
whereas the entire domain is taken into account for the body force
test. As a result, more tracking errors are included in the body force
test, leading to a greater initial error.

Increasing the mechanical regularization however allows to correct this
initial error, smooths the form of the noise, and improves drastically
the estimation for all noise levels. This result is also consistent
with~\cite{genet_finite_2023}, who observed that the regularization
allowed to decrease the error on the displacement field obtained with
the tracking process for all noise levels. This result is consistent
with the form of Equation~(\ref{eq:expression_E}) as well. Noise is
indeed included in $\newefvec{U}_{\mathrm{meas}}$, but corrected during
the tracking. As a result, it is more structured and hence leads to a
better estimation. Note that the estimation for the boundary force test
is slightly better than for the body force test, when regularization is
added during the tracking process, as well as when noise is added
directly to displacements, {for the same reasons as those
previously mentioned.}

Overall, our new formulation of the EGM provides a less accurate
estimation than the other methods studied for comparison. This finding
aligns with classical behaviors documented in the
literature~\cite{roux_optimal_2020}, which suggest that direct methods
typically yield significantly poorer results compared to indirect
methods, such as the FEMU method, {although it is significantly
less expensive.} Furthermore, a carefully selected virtual field for
the VFM is expected to display a better performance than the EGM, as
the latter considers all virtual fields, including the suboptimal ones.
{However, defining a good virtual field  might be challenging for
complex geometries.} {Additionally}, the EGM is less accurate than
{other} methods {only} when noise is added to the displacement
field{:} it is not the case anymore when noise is added to images
{f}or a good regularization choice. The EGM can hence be very well
adapted for the estimation, depending on the desired level of accuracy.
Even when higher accuracy is wanted, the EGM may be very useful,
insofar as it can be used to initialize the estimation and therefore
improve and fasten the identification process.

\subsection{Estimation performance of the VFM}

The results presented in this article showed that the estimation with
the VFM leads a better accuracy and noise robustness than the EGM. This
is consistent with the fact that well-chosen virtual fields should lead
to a better estimation quality than the EGM, which includes all
possible virtual fields, even the ones leading to a poor estimation
quality. The noise sensitivity however differs for both studied virtual
fields---the plane wave and the virtual field constructed
from{~\cite{deng_fenics_2023}}. This difference is noticeable for
both boundary force and body force tests, and is in accordance with the
literature:~\cite{avril_sensitivity_2004} has indeed shown that all
virtual fields do not display the same sensitivity to noise. If new
approaches have been designed to create virtual fields more robust to
noise~\cite{grediac_virtual_2006, marek2017sensitivity}, there does not
exist a general method to create an optimal virtual field---i.e.{,}
the virtual field, which displays the highest robustness to
noise---yet.  

The variability in accuracy was for example illustrated in this work
when the virtual field was chosen as a plane wave, by varying the
frequency. This variation of accuracy depending on the frequency of the
plane wave {can} be explained by studying the cost function of the
method. As for the EGM, since the generalization is
straightforward---the expression of the identified Young modulus is
obtained directly from Equation~(\ref{eq:E_surface/volume})---, let us
conduct the analysis for linearized quantities without loss of
generality. Additionally, for illustrative purposes, let us study the
expression of the Young modulus for a boundary force---which can easily
be replaced by the body force. In that case, we have shown that 
{\begin{equation} \label{eq:E_volume/surface}
E^{\Gamma} = \frac{W_{\mathrm{ext}}^{\Gamma}(\underline{U}^*)}
{W_{\mathrm{int}}(\underline{U}_{\mathrm{meas}}; E=1, \underline{U}^*)} = 
\frac{\displaystyle \int_{\Gamma}^{} \tau \underline{n} \cdot \underline{U}^* \,\mathrm{d}\Gamma}
{\displaystyle \int_{\Omega}^{}\underline{\underline{\sigma}}(\underline{U}_{\mathrm{meas}} ;E=1) :
(\underline{\underline{\nabla}} \, \underline{U}^*)_{\mathrm{sym}}\,\mathrm{d} \Omega}.
\end{equation}}\unskip

Both numerator and denominator depend linearly on $\underline{U}^*$. In
order to find $E^{\Gamma}$, an admissible virtual field
$\underline{U}^*$  needs to be defined. As reminded in
Equation~(\ref{eq:VF-planewave}), the family of virtual fields chosen
as a plane wave can be written $\underline{U}^* = 
\sin(\underline{k}\cdot(\underline{X}-\underline{X}_0)) \underline{e}$.
As a result, the gradient of $\underline{U}^*$ is:
{\begin{equation}
\underline{\underline{\nabla}} \, \underline{U}^* =
|\underline{k}| \cos(\underline{k} \cdot
(\underline{X}-\underline{X}_0)) \underline{e} \otimes \underline{e}.
\end{equation}}\unskip

Note that in order to be well-defined, our virtual field should not
cancel $W_{\mathrm{ext}}$ at the right boundary. As a result, $\Delta$
({appearing through the relation $|k| = {2 \uppi}/{\Delta}$})
should not be of the form ${\alpha}/{l}$, where $l$ is the
characteristic length of the cube, and $\alpha$ an integer. Based on
the previous equation, ${|\underline{k}|}$ can  be factorized at
the denominator, such that:
{\begin{equation}\label{eq:expression-E-planewave}
\displaystyle
E^{\Gamma} = \frac{\displaystyle \int_{\Gamma}^{} \tau \underline{n}
\cdot \sin(\underline{k}\cdot(\underline{X}-\underline{X}_0))\underline{e} \,
\mathrm{d}\Gamma}{|\underline{k}| \displaystyle
\int_{\Omega}^{}\underline{\underline{\sigma}}
(\underline{U}_{\mathrm{meas}} ;E=1) : 
\cos(\underline{k} \cdot (\underline{X}-\underline{X}_0)) 
\underline{e} \otimes \underline{e}\,\mathrm{d} \Omega}.
\end{equation}}\unskip

By observing Equation~(\ref{eq:expression-E-planewave}), we can note
that noise is present at the denominator, where ${1}/{\Delta}$ is a
factor. Noise is indeed included in the term
$\underline{\underline{\sigma}}(\underline{U}_{\mathrm{meas}}; E=1)$,
through $\underline{U}_{\mathrm{meas}}$, which depends on the noise
$\underline{g}$ through the relation:
{\begin{equation}
\underline{g} = \underline{U}_{\mathrm{meas}} - \underline{U}_0,
\end{equation}}\unskip
where $\underline{U}_0$ is the {measurement} without noise. As a
result, Equation~(\ref{eq:expression-E-planewave}) can also be written:
{\begin{equation}\label{eq:expression-E-planewavebis}
\displaystyle
E^{\Gamma} = \frac{\displaystyle \int_{\Gamma}^{}\tau\underline{n}
\cdot \sin(\underline{k}\cdot(\underline{X}-\underline{X}_0))
\underline{e}\,\mathrm{d}\gamma}{|\underline{k}| \displaystyle
\int_{\Omega}^{}\underline{\underline{\sigma}}(\underline{U}_0
+\underline{g} ; E=1) : \cos(\underline{k} \cdot
(\underline{X}-\underline{X}_0)) \underline{e} \otimes \underline{e} \,
\mathrm{d}\Omega}.
\end{equation}}\unskip
For linear elasticity, $\underline{\sigma}(\underline{U}_0
+\underline{g} ; E=1)$ depends in reality linearly on the first
derivative of $\underline{U}_{\mathrm{meas}}$, and hence linearly on
the first derivative of $\underline{U}_0 + \underline{g}$. It is
therefore clear that the impact of noise will increase as the value of
$\Delta$ decreases. This phenomenon will in particular be noticeable
for unstructured noise, i.e.{,} random Gaussian noise added
directly to the displacement field. However, in general cases, and
contrarily to the EGM, it is impossible to predict whether the
parameter will be overestimated or underestimated with the VFM, with a
virtual field chosen as a plane wave. In large strains, the external
work indeed depends on the measured---noisy---displacement field. For
example, in the case considered here, $W_{\mathrm{ext}}^{\Gamma} =
\int_{\Gamma}^{} T \underline{\underline{F}}^{-T}\cdot \underline{N}
\cdot \underline{U}^* J \,\mathrm{d} \Gamma$ (in large strains) depends
on the noisy displacement trough $\underline{\underline{F}}^{-T}$ and
$J$. The expression of the Young modulus therefore depends on the noise
through its numerator and denominator. It is hence not straightforward
to predict the behavior of the estimation.

{T}he estimation quality improves greatly for the VFM when noise is
added to images---for a good regularization---compared to the case
where synthetic data is generated by adding noise to
displacements{,} due to the smoothing of noise with the second data
generation method. Moreover, and contrarily to the EGM, the estimation
seems more robust with the VFM when a body force is applied than for
the boundary force test. We assume it is because less noise is
considered at the numerator for the body force test---only included in
$J$--- than for the boundary test---included in $J$ and in
$\underline{\underline{F}}^{-T}$---, as detailed in the Methods
section.

Overall, the VFM seems to give better results than the EGM. However, it
was possible to define the plane wave as virtual field only given the
simplicity of the considered geometry. Choosing a good kinematically
admissible virtual field for more complex problems is rarely as
straightforward. Methods have been proposed to overcome this issue,
such as the one studied in this article, proposed
by~\cite{deng_fenics_2023}. If this particular method allows to build a
kinematically admissible displacement field, with good estimation
performances, it however transforms the VFM in a {in}direct method,
i.e., the method is not entirely direct anymore as it requires finite
element computations for given steps of the estimation. As a result, if
our new formulation of the EGM leads to a lower estimation quality, it
may still be a better option than the VFM { for computationally
expensive problems}.

\subsection{Estimation performance of the FEMU method}

In general, the noise robustness of indirect methods is higher than for
direct methods. These methods are however very expensive, since they
require finite element computations at every step of the optimization.
The robustness to noise of the FEMU method was studied here for
comparison. The FEMU method is indeed a direct method known to be very
robust to noise~\cite{roux_optimal_2020}, which is what is observed for
all cases of our study. The only case where the estimation is not good
is when noise is added to images, for a fine mesh generated without the
multiresolution method. As explained in the Results section,
\cite{genet_finite_2023} noticed this phenomenon in his study. The
finer the mesh, the bigger the function space for the solution, which
can lead to ill displacement fields solutions. The problem is however
solved when using the multiresolution method.

Using {measurements} generated by adding noise to the displacements
or to the images does not seem to significantly impact the estimation
quality of the FEMU method. We assume this is due to the form of the
cost function: unlike the direct methods studied here, which include
first derivatives of the noise in their cost functions, and are
therefore strongly impacted when the noise is unstructured, the FEMU
method only aims to find the Young modulus that leads to a computed
displacement field, which best matches the noisy measurements.
Consequently, even when noise is added directly to the displacement
field, the {latter} will not be able to fit exactly the
unstructured data, {and} the solution displacement field matching
best the data will automatically ``average'' the noise, which limits
{its} impact on the estimation.

The robustness to noise of the FEMU method is therefore much higher
than for the VFM and the EGM, and should be preferred when high
accuracy is needed. This method might however be too expensive. In this
case, our new formulation of the EGM may be a good alternative. It
indeed displays a good noise robustness and presents the advantage of
being a direct method, and is therefore far less expensive than the
FEMU method.

\subsection{Comparison of the estimation performance of the three
methods}

The choice of an estimation method is rarely simple. It indeed usually
requires a good trade-off between accuracy and computation time. For
example, the FEMU method represents the best option in terms of
accuracy and robustness to noise and model errors. This method is
however very expensive, since it requires a finite element computation
at each iteration. Consequently, for complex problems requiring many
inverse problems, the FEMU method will induce very high computation
times and might not be adapted. In this case, direct methods {such
as the VFM and the EGM} may be better options.

Based on the results presented in this article, the VFM seems at first
sight to be the best option. The estimation with this method indeed
leads to more accurate results: in general, the estimated Young modulus
is closer to the ground-truth value, and is identified with less
uncertainty than with the EGM.  However, the results presented with the
VFM were obtained for two well-chosen virtual fields. The choice of a
plane wave as virtual field was only possible given the simplicity of
the geometry. Defining a kinematically admissible displacement field is
however much more difficult for more complex geometries. In order to
build a good virtual field in such cases, the {method} proposed
by~\cite{deng_fenics_2023} {is particularly appropriate but is (as
mentionned before)} not a fully direct method anymore. Even if less
accurate than the VFM for a well-chosen virtual field, the EGM is
therefore a good option, since building a good virtual field for a
complex problem is not straightforward {or might be expensive}.

Note however that {the} robustness to noise {of the EGM} is
very low when the synthetic {measurements} used for the estimation
are generated by adding noise directly {to} the displacement
field{, but improves greatly when }noise {is added} to images,
and the displacement field {retrieved} with tracking associated to
the Equilibrium Gap regularization. Other methods, such as filtering
high noise frequencies or generating the random Gaussian noise on a
coarser mesh and projecting it on the finer mesh, could be considered
to {create synthetic measurements with smooth and physical }noise.
We however preferred our method of generating noisy images, because of
its closeness to real data processing. 

When using the second method of synthetic data generation, the
estimation with the EGM becomes very good. Until high noise levels, the
average of the estimated Young modulus distribution is close to the
ground truth value, and the standard deviation is relatively small.
However, the EGM remains overall less accurate than the FEMU method,
which may be a better option when high accuracy is needed, even though
it is very expensive.

\section{Conclusion}

When full-field measurements are available for inverse problems in
large hyperelastic strains, different estimation methods can be
considered. All methods present advantages and drawbacks in terms of
accuracy and in terms of noise and/or model errors robustness. Four
different methods were investigated in this article: the FEMU method,
very robust even though quite expensive, the VFM, with two different
virtual fields---a direct method, with virtual fields chosen as plane
waves, and an indirect method, with virtual fields
from~\cite{deng_fenics_2023}---and the EGM, which is less robust but
also less expensive than the indirect methods. Note that no general
formulation of the EGM has been proposed in large hyperelastic strains.
The originality of our work hence lies in a new formulation of the EGM
in large hyperelastic strains ({which is directly based on a
recently proposed continuous formulation and consistent discretization
of the EG principle \cite{genet_finite_2023}}), and its performance
quantification, by comparing {it} with the other three
{aforementioned} methods for {various} noise levels.

The {quantitative} study of the impact of the robustness to noise
of each method relies on the generation of synthetic data. Two methods
for generating such data are investigated here. The first method is
based on a classical approach in uncertainty quantification: a random
Gaussian noise is added to a displacement field computed with the model
using reference parameters. The second method consists also to first
compute a displacement field with the model with parameters chosen at
reference values. Images of the reference and deformed configurations
are then created, to which a random Gaussian noise is added. The
displacement field between the two configurations is finally retrieved
with a tracking process, to which mechanical regularization is added.
The second method therefore allows to generate noisy displacement
fields closer to experimental data, insofar as the noise produced with
the first method is highly unstructured.

The results obtained allow to quantitatively study the estimation
quality of each method, and in particular of our new formulation of the
EGM. As expected, the FEMU method gives very good results: for all
cases studied{,} this method displays a very high robustness to
noise. The VFM is on contrary less accurate, but far less expensive
when using plane waves as virtual fields. Additionally, this method is
in average more accurate than the EGM. 

The VFM however depends strongly on the chosen virtual {field}. It
indeed requires the choice of a ``good'' kinematically admissible
virtual field, which did not represent any difficulty in this article
given the simplicity of the geometry considered, but might be less
straightforward for more complex problems. General methods, such as the
one presented in~\cite{deng_fenics_2023}, can be developed to find
``good'' virtual fields for more complex {geometries}. These
methods however require model computations for building the virtual
fields, which transforms the VFM into a more expensive method.

The EGM therefore represents a good alternative for the estimation
process. As expected, the estimation with the EGM is less accurate than
with the VFM, since it takes into account all kinematically admissible
virtual fields, including the suboptimal ones. Additionally, the EGM is
not robust to noise at all when the synthetic measurements are
generated by adding noise to displacements. However, when the synthetic
data is generated by adding noise to images, with tracking preformed
using the Equilibrium Gap mechanical regularization, the estimation
quality improves drastically. {Therefore, the EGM could be a good
choice for complex estimation problems, even involving high noise
levels.}

{Note that only the EGM, the VFM and the FEMU method have been
investigated in the scope of this article. However, other methods can
also offer a good compromise between robustness and cost. For example,
the Reconditioned EGM~\cite{roux2008digital, roux_optimal_2020} could
be an alternative, with higher noise robustness than the VFM and the
raw EGM, but smaller computational costs than the FEMU method. Note
also that we only identified one parameter in this article. We expect
the uncertainty to increase for our EGM formulation with the number of
identified parameters, which we will investigate in future studies.}

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds
from any organization that could benefit from this article, and have
declared no affiliations other than their research organizations.

\section*{Acknowledgments}

We would like to thank Colin Laville, S\'ebastien Imperiale and
Philippe Moireau for fruitful discussions. We would also like to thank
the French National Research Agency (ANR) for providing financial
support to this project (Grant no: ANR-19-CE45-0007).

\CDRGrant[ANR]{ANR-19-CE45-0007}

\back{}

\appendix{}

\def\appendixlabel{Appendix \Alph{section}}

\section{Quantification of the impact of Poisson ratio errors on the
estimation}\label{app:appendixA}

In this article we have studied the estimation quality of our new
formulation of the EGM in large strains and its robustness to noise,
and have compared it to the performances of the FEMU method and of the
VFM. However, other factor than the noise can influence the performance
quality of an estimation method. In particular, each method presents
different sensitivities to model errors. To investigate the robustness
of the different methods to model errors, we choose to study in this
Appendix the impact of errors of the Poisson ratio value fixed during
the estimation. During the synthetic data generation, the Poisson ratio
is indeed empirically fixed at 0.3. Model errors on this coefficient
can be hence introduced by fixing the Poisson ratio at a value
different from the ground-truth during the estimation, e.g., 0.24 or
0.36. This will automatically induce a bias in the estimation.

Figure~\ref{fig:model-errors-nu} presents the error distribution of the
estimated Young modulus depending on the noise for the four
methods---the EGM, the VFM with a plane wave as virtual field, the VFM
from~\cite{deng_fenics_2023} and the FEMU method---, when committing an
error of ${-}$20\% (in red), ${-}$10\% (in orange), 0\% (in green),
${+}$10\% (in blue), ${+}$20\% (in pink) on the Poisson ratio. The
estimation was performed using synthetic data generated by adding noise
{to} images, for a body force test.

\begin{figure}
\includegraphics{fig11}
\caption{\label{fig:model-errors-nu}Error distribution on the
identified Young modulus depending on the noise, by making an error of
${-}$20\% (red), ${-}$10\% (orange), 0\% (green), ${+}$10\% (blue) and
${+}$20\% (pink) on the fixed Poisson ratio, for the EGM (top left), the
VFM with a plane wave as virtual field (top right), the VFM
from~\cite{deng_fenics_2023} (bottom left) and the FEMU method (bottom
right).}
\end{figure}

For all methods, underestimating the Poisson ratio leads to
underestimating the Young modulus, and conversely. This result can be
explained by the (negative) linear dependency of the Poisson ratio to
the Young modulus for isotropic materials (at fixed bulk modulus): when
the Poisson ratio increases, we can expect the Young modulus to
increase, which is indeed what is observed here. Additionally, adding
an error on $\nu$ during the estimation shifts the average and standard
deviation of the error distribution. The global tendencies of the
distributions however remain similar for all methods.

The results however show that model errors have a great impact on the
estimation accuracy for direct methods. For a Poisson ratio set to
120\% of its ground-truth value during the estimation process, the
error committed for a SNR of ${+}\infty$---i.e., no noise on the
images---is around ${-}$18\% for the EGM and for the VFM using a plane
wave, while it is around 6\% for the FEMU method and the VFM
from~\cite{deng_fenics_2023}.

This result can be explained by the direct dependency of the direct
methods on the Poisson ratio in the cost function. Small variations of
this coefficient therefore directly impact the expression of the Young
modulus. The indirect methods, on the other hand, do not depend
directly on the Poisson ratio, which explains their higher robustness
to model errors.

Our new formulation of the EGM seems therefore very sensitive to model
errors, which is a classical observation for direct methods. This is
however a limitation of our new EGM formulation to keep in mind when
selecting methods for inverse problems.

\section{Quantification of the impact of regularization errors on the
estimation}\label{app:appendixB}

In this appendix, we study the impact of errors on the regularization
value used during the tracking process on the estimation. For the body
test, a volume regularization is indeed added during the tracking
process. This regularization enforces that the divergence of the stress
tensor should be as close as possible to an imposed body force of 0.3
mN/mm$^3$. However, in some experimental tests, the exact value of the
body force might not be known, and the values set for the estimation
might therefore differ from the ground-truth value of 0.3 mN/mm$^3$,
which would induce model errors.

We therefore choose to study the impact of error on the regularization
term on the estimation. Figure~\ref{fig:model-errors-b} presents the
error distribution of the estimated Young modulus depending on the
noise for the four methods---the EGM, the VFM with a plane wave as
virtual field, the VFM from~\cite{deng_fenics_2023} and the FEMU
method---, when committing an error of ${-}$20\% (in red), ${-}$10\% (in
orange), 0\% (in green), ${+}$10\% (in blue), ${+}$20\% (in pink) on
the volume regularization term. The estimation was performed using
synthetic data generated by adding noise to images for a body force
test.

\begin{figure}
\includegraphics{fig12}
\caption{\label{fig:model-errors-b}Error distribution on the identified
Young's modulus, by making an error of ${-}$20\% (red), ${-}$10\% (orange),
0\% (green), ${+}$10\% (blue) and ${+}$20\% (pink) on the volume
regularization value.} 
\end{figure}

The model errors on the volume regularization term seems to mostly
impact the EGM; the impact on other methods is indeed marginal. For a
volume regularization term set to 120\% of the body force ground-truth
value, the average error committed for a SNR of  ${+}\infty$---i.e., no
noise on the images---is around ${-}$10\% for the EGM and close to 0\%
for all other methods. This phenomenon can be explained by the high
sensitivity to the measured displacement field of the EGM, already
observed in the Results section. Any bias introduced on the
regularization term will induce a variation of the synthetic measure.
Given its sensitivity to such variations, a high impact on the
estimation with the EGM can therefore be expected. Since all other
methods are far less sensitive to small accuracy variations of the
measured displacement field, estimation with such methods should be
robust to a biased regularization, which is what is indeed observed.

Our new formulation of the EGM is therefore more sensitive to
regularization errors than non-direct methods. This limitation should
be kept in mind when selecting this method to solve problems for which
the load in not perfectly known.

\section{{Estimation performance for a study in 3D}}
\label{app:appendix3d}

{In this article, we proposed a study on a very simple geometry (a
2D square), for a simple problem (only one parameter identified) as a
proof of concept. Indeed, the EGM had never been used in finite
deformation cases. Since the method is very sensitive to noise, and as
we studied the estimation performances for high noise levels,
non-linear problems and fine meshes, it was not straightforward that
the estimation would perform well. However, the method can easily be
extended to more complex cases. As an illustration, Figure~\ref{fig:3d}
presents the estimation quality, when the problem is in 3D. In this
case, a cube is clamped on the left face, and a homogeneous pressure of
0.3~kPa is applied on the right face.}

\begin{figure}
\includegraphics{fig13}
\caption{\label{fig:3d}Estimation error using the EGM for a boundary
force test for 2D test (dark green) and a 3D test (light green), when
noise is added to displacements (left) and when noise is added to
images with a regularization of 0.2 (right).} 
\end{figure}

The estimation quality is similiar for both cases. In particular,
for a good regularization level, the estimation for the EGM in 3D
displays a good accuracy.

\section{{Estimation performance for a perfectly and a partially
known distribution of boundary forces}}\label{app:appendixC}

{We studied in this article the estimation performance of the
different methods for a known boundary forces distribution. This case
was chosen to evaluate the estimation performance of the boundary term
introduced in our new EGM formulation, for cases where the boundary
force is known, which is the case for some biomechanical
applications~\cite{patte2022estimation}. However, the whole
distribution is not always known for all experiments: sometimes, only
the resultant load is known. As a result, we studied in
Figure~\ref{fig:known-unknown_resultant} the impact on the estimation
performance of the EGM using the exact boundary force distribution or
only the resultant. In the case where only the resultant is known,
Equation~(\ref{eq:boundaryterm_egm}) becomes:
{\begin{equation}\label{eq:boundaryterm_egm_known_resultant}
{L}^{\mathrm{EGM}, F}_{\Gamma} := \frac{1}{2} \left(\displaystyle
\int_{\Gamma}^{} (\underline{\underline{P}} \cdot \underline{N} -
\underline{T})\,\mathrm{d}\Gamma\right)^2.
\end{equation}}}\unskip

\begin{figure}
\includegraphics{fig14}
\caption{\label{fig:known-unknown_resultant}Estimation error using the
EGM for a boundary force test when the whole distribution is known
(blue) and when only the resultant is known (green), when noise is
added to displacements (left) and when noise is added to images with a
regularization of 0.2 (right).}
\end{figure}

{When noise is added to images, the estimation quality is
equivalent whether the exact distribution is known or not. However,
when noise is added to displacements, the estimation is of bad quality
for both cases: the standard deviation for the case where only the load
resultant is known is additionally much higher, which could have been
expected as in that case, the estimation holds less information than
for the case where the whole distribution is known.When using the EGM
with noise added to images and an appropriate regularization, the
estimation quality is however very good. Our method could therefore be
a good option, even for experimental data which only provide the
resultant load and not the whole boundary force distribution.}

\begin{figure}
{\vspace*{-4pt}}
\includegraphics{fig15}
{\vspace*{-4pt}}
\caption{\label{fig15}Estimation error using the VFM for a boundary
force test when noise is added to displacements (left) and when noise
is added to images with a regularization of 0.2 (right), for the
virtual field chosen as a plane wave (blue), the virtual field
from~\cite{deng_fenics_2023} (red) and the virtual field introduced in
Equation~(\ref{eq:simple-expression}) (pink).}
{\vspace*{-4pt}}
\end{figure}

\section{{Estimation performance of an additional virtual field}}
\label{app:appendixD}

{The two virtual fields used in this article were chosen because on
the one hand, the plane wave was easy to define on our geometry, and
one parameter (the frequency) could easily be varied to  study its
impact on the estimation, and on the other hand the virtual field
from~\cite{deng_fenics_2023} can be applied to problems with more
complex geometries. However, as highlighted by one of our reviewers,
these fields are respectively more adapted to non homogeneous problems
and to more complex problems. As a result, in this Appendix, we also
study the estimation quality for the simplest virtual field possible,
expressed as:
{\begin{equation}\label{eq:simple-expression}
U^* = \begin{pmatrix}
\frac{(\underline{X}[0]-\underline{X}_0[0])}{l} \\
0
\end{pmatrix}.
\end{equation}}}\unskip

The performances of the virtual field introduced in
Equation~(\ref{eq:simple-expression}) are presented in Figure~\ref{fig15}, 
and are very similar to those of the
two other virtual fields. Indeed, for noise added to images, there is
no noticeable difference between the three virtual fields.
Additionally, when noise is added to displacements, for SNRs larger
than 10, the performance of the new virtual field is similar to the
performances of the plane wave. As a result, this virtual field
displays a good estimation quality and its simplicity makes it very
adapted to our problem.

{\vspace*{-4pt}}
\section{Code for
Figures~\myRef{fig:compression_test_coarse}, 
\myRef{fig:estimated_Young_h=0.05-compression}, 
\myRef{fig:estimated_Young_h=0.1-refined=False-volume}, 
\myRef{fig:estimated_Young_h=0.1-refined=True-volume}}

{\vspace*{-4pt}}
\subsection{Imports}
\noindent\inlinefig{fx01}

\subsection{Parameters}
\subsubsection{Geometry}
\noindent\inlinefig{fx02}

\subsubsection{Loading and boundary conditions}
\noindent\inlinefig{fx03}

\subsubsection{Material behavior}
\noindent\inlinefig{fx04}

\subsubsection{Images}
\noindent\inlinefig{fx05}

\subsubsection{Varying parameters for the estimation}
\noindent\inlinefig{fx06}

\subsection{Synthetic measurements}
\subsubsection{Synthetic images}
\noindent\inlinefig{fx07}

\subsubsection{Ground-truth motion}
\noindent\inlinefig{fx08}

\subsubsection{Tracking}
\noindent\inlinefig{fx09}
\clearpage

\noindent\inlinefig{fx10}

\clearpage
\subsection{Generating plots}
\subsubsection{Helpers}
\noindent\inlinefig{fx11}

\clearpage
\noindent\inlinefig{fx12}

\noindent\inlinefig{fx13}

\clearpage

\subsubsection{Plots}
%{\vspace*{-.5pc}}
\noindent\inlinefig{fx14}

{\vspace*{-1pc}}
\section{Code for Figure~\myRef{fig:impact-delta}}
{\vspace*{-.5pc}}
\subsection{Imports}
\noindent\inlinefig{fx15}

{\vspace*{-1pc}}
\subsection{Parameters}
{\vspace*{-.5pc}}
\subsubsection{Geometry and varying parameters}
\noindent\inlinefig{fx16}

{\vspace*{-.5pc}}
\subsubsection{Material behavior}
\noindent\inlinefig{fx17}

%{\vspace*{-.5pc}}
%\clearpage
\subsubsection{Loading}
\noindent\inlinefig{fx18}

%%%{\vspace*{-.5pc}}
\subsection{Generating plots}
%%%{\vspace*{-.25pc}}
\subsubsection{Helpers}
\noindent\inlinefig{fx19}

\clearpage
\noindent\inlinefig{fx20} %%%%% merged with 19

%\clearpage

\subsubsection{Plots}
\noindent\inlinefig{fx21} %21 changed to 20

%\clearpage

\section{Code for Figures~\myRef{fig:model-errors-nu}--\myRef{fig:model-errors-b}}
\subsection{Imports}

\noindent\inlinefig{fx22}

\subsection{Parameters}
\subsubsection{Geometry}
\noindent\inlinefig{fx23}

\subsubsection{Loading and boundary conditions}
\noindent\inlinefig{fx24}

\subsubsection{Material behavior}
\noindent\inlinefig{fx25}

\subsubsection{Images}
\noindent\inlinefig{fx26}
\clearpage

\subsubsection{Bias}
\noindent\inlinefig{fx27}

\subsection{Synthetic measurements}
\subsubsection{Synthetic images}
\noindent\inlinefig{fx28}

\subsubsection{Ground-truth motion}
\noindent\inlinefig{fx29}

\subsubsection{Tracking}
\noindent\inlinefig{fx30}

\clearpage

\noindent\inlinefig{fx31}

\subsection{Generating plots}
\subsubsection{Helpers}
\noindent\inlinefig{fx32}

\noindent\inlinefig{fx33}

\noindent\inlinefig{fx34}

{\vspace*{4pt}}

\noindent\inlinefig{fx35}

\subsubsection{Plots}
\noindent\inlinefig{fx36}

{\vspace*{-1pc}}
%\bibliographystyle{crunsrt}
%\bibliography{crmeca20240698}
\printbibliography
\refinput{crmeca20240698-reference.tex}

\end{document}
