Inverse deformation analysis: an experimental and numerical assessment using the FEniCS Project
Arnaud Mazier, Alexandre Bilger, Antonio E. Forte, Igor Peterlik, Jack S. Hale, Stéphane P.A. Bordas, Institute of Computational Engineering, Department of Engineering, University of Luxembourg, Esch-sur-Alzette, Luxembourg., Harvard University, Cambridge, USA., Department of Electronics, Information, Bioengineering, Politecnico di Milano, Milan, Italy., Institute of Computer Science, Masaryk University, Czech Republic., Institute of Research, Development Duy Tan University, Danang, Vietnam
HHighlights
Inverse deformation analysis: an experimental and numerical as-sessment using the FEniCS Project (cid:63)
Arnaud Mazier, Alexandre Bilger, Antonio E. Forte, Igor Peterlik, Jack S.Hale, St´ephane P.A. Bordas • Inverse deformation algorithm can retrieve the undeformed configura-tion of soft objects • Implementation of the inverse deformation algorithm in the FEniCSProject software a r X i v : . [ c s . C E ] F e b nverse deformation analysis: an experimental andnumerical assessment using the FEniCS Project Arnaud Mazier a , Alexandre Bilger a , Antonio E. Forte b,c , Igor Peterlik d ,Jack S. Hale a , St´ephane P.A. Bordas a,e, ∗ a Institute of Computational Engineering, Department of Engineering, University ofLuxembourg, 6, avenue de la Fonte, L-4364 Esch-sur-Alzette, Luxembourg. b Harvard University, 29 Oxford St, Cambridge MA 02138, USA. c Department of Electronics, Information and Bioengineering, Politecnico di Milano,Milan, 20133 Italy. d Institute of Computer Science, Masaryk University, Czech Republic. e Institute of Research and Development Duy Tan University, K7/25 Quang Trung,Danang, Vietnam.
Abstract
In this paper we develop a framework for solving inverse deformation prob-lems using the FEniCS Project finite element software. We validate ourapproach with experimental imaging data acquired from a soft silicone beamunder gravity. In contrast with inverse iterative algorithms that require mul-tiple solutions of a standard elasticity problem, the proposed method cancompute the undeformed configuration by solving only one modified elas-ticity problem. This modified problem has complexity comparable to thestandard one. The framework is implemented within an open-source pipelineenabling the direct and inverse deformation simulation directly from imaging (cid:63)
This study was supported by European Union’s Horizon 2020 research and innovationprogramme under the Marie Sklodowska-Curie grant agreement No. 764644, No. 798244and the financial support of the European Research Council Starting Independent ResearchGrant (ERC StG grant agreement No. 279578). Jack S. Hale is supported by the NationalResearch Fund, Luxembourg, and cofunded under the Marie Curie Actions of the EuropeanCommission (FP7-COFUND) Grant No. 6693582. ∗ Corresponding author
Email addresses: [email protected] (Arnaud Mazier), [email protected] (Alexandre Bilger), [email protected] (Antonio E. Forte), [email protected] (Igor Peterlik), [email protected] (Jack S.Hale), [email protected] (St´ephane P.A. Bordas)
Preprint submitted to Elsevier March 1, 2021 ata. We use the high-level Unified Form Language (UFL) of the FEniCSProject to express the finite element model in variational form and to au-tomatically derive the consistent Jacobian. Consequently, the design of thepipeline is flexible: for example, it allows the modification of the constitutivemodels by changing a single line of code. We include a complete workingexample showing the inverse deformation of a beam deformed by gravity assupplementary material.
Keywords:
Inverse deformation, rest position, undeformed configuration,SOFA, FEniCS Project.
1. IntroductionMotivation.
The organization of a standard biomechanical deforma-tion analysis pipeline typically proceeds as follows. First, by using imagingtechniques such as Magnetic Resonance Imaging (MRI) a segmented imageof the region of interest is obtained. This segmented image is then meshedso that it can be used as input for a finite element simulation. The meshis considered as the initial or undeformed (or reference) configuration of anelastic body. Then, by applying external forces to this elastic body we canfind its deformed (or current) equilibrium configuration.Conversely, an inverse deformation analysis allows us to find the unde-formed configuration of a body knowing its deformed configuration. In thecase of an object subject to gravity, the undeformed configuration can beseen as a theoretical gravity-free configuration. Consequently, determiningthe rest-position of an organ is of interest in many (bio)mechanical prob-lems. For example, in abdominal aortic aneurysms to compute the residualstresses [1, 2], or in breast cancer as an intermedial configuration betweenthe imaging and surgical stance [3]. Besides, this approach can also be usedin problems of industrial interest such as tire or turbine blade design [4, 5].
Problem statement.
The objective of inverse deformation analysis isto determine the undeformed configuration of an object such that it attainsa known deformed configuration under the action of a known loading. It isimportant to note the distinction between inverse deformation analysis andcommon inverse problems. In a typical inverse problem, we might assume weknow the applied forces, the initial and deformed configuration, and the goalis to determine the model parameters that minimize some distance (met-ric) between initial and deformed configurations. In an inverse deformation2nalysis, we assume we know the applied forces, boundary conditions, modelparameters, and the deformed configuration. The objective is to determinethe undeformed configuration that would lead to the deformed configurationif the external forces were to be applied.
Background.
Several authors have tackled the problem of inverse defor-mation analysis using a variety of strategies. To the best of our knowledge [6]was the first to propose exchanging the role of the deformed and undeformedconfigurations, i.e. to express the displacement of the body as a function ofthe deformed state. The study was limited to plane strain deformations anduniform extension. [7] applied the same formalism to a homogeneous elasticmaterial, without body force. He showed the equivalence of the equilibriumequations if the initial and deformed configurations are interchanged as wellas the volumetric strain energies. The results provided by this approach areshown to be commensurate with those of [6] but are based on dual relationsbetween the initial configuration and the deformed configuration. [8] useda variational principle to achieve the same as [7] and showed the validityof the approach for different elastic materials. More recently [9] mathemati-cally analyzed the Schield transformation and the proven inverse deformationtheorem. The theorem states that if a particular deformation is supportedwithout body force for a specific strain energy W , then the inverse deforma-tion is another energy W ∗ , derived from the first: W ∗ ( F ) = det( F ) W (cid:0) F − (cid:1) ,where F is the deformation gradient.[10, 11] introduced the reparameterization of the weak form of the for-ward problem of finite elasticity as a solution method for the inverse prob-lem. This approach only requires C continuity and has a direct physicalconnection to the problem. Additionally, the procedure eliminates boundarycondition difficulties, can be straightforwardly implemented using standardforward numerical methods, and can deal with both compressible or incom-pressible materials.Inspired by [11] (Eulerian model) and [12] (Arbitrary Lagrangian-Eulerian(ALE)), [13] rewrote the constitutive equations in terms of Lagrangian vari-ables. This manipulation makes the inverse analysis code changes limited tothe finite element residual and Jacobian computations, contrary to Eulerianor ALE variables. The formulation is convenient and allows to solve inversedesign problems such as finding the unloaded shape of a turbine blade underknown loading. But few drawbacks arise such as the difficulty of deriving andimplementing the consistent Jacobian of the finite element formulation. De-spite the usefulness of the approach, to our knowledge, this type of analysis3s still not available in any widely used commercial simulation software.Iterative methods identify the undeformed configuration based on severalforward calculations. The algorithm is introduced by [14] with a fixed-pointmethod for elastostatic problems and then generalizes as the backward dis-placement method by [15] for patient-specific blood vessel simulations. Theiterative algorithm of Sellier has been widely applied to many image-basedbiomechanical simulations, mainly thanks to its algorithmic simplicity and itsability to use a standard non-linear elasticity simulation software [3]. How-ever, when applied to strongly non-linear problems resulting from materialor geometric non-linearity, the algorithm lacks robustness. Furthermore, it-erative methods usually require at least one non-linear elasticity problemsolution, resulting in higher costs compared with the approach of [13].In the computer graphics community, [16] used Asymptotic NumericalMethods (ANM) to compute the rest-shape of elastic objects with a neo-Hookean material model. The ANM considers a parametrized version ofthe static equilibrium: f ( x, X ) + λg = 0, where g is gravity, λ a loadingparameter and f are the internal forces with the given deformed configuration x and the unknown rest-configuration X . Then, the algorithm incrementallycomputes the asymptotic expansion of the curve in ( X, λ ) space until λ = 1,which corresponds to the rest-position. In this study, ANM offers superiorperformance, robustness, and convergence speed over traditional Newton-type methods for highly nonlinear material models. But the major drawbackof the method is the complexity of changing the model formulation. Indeed,using a different material model implies to establish a different quadraticrelationship between Cauchy stresses and the rest-position, then derivingthe asymptotic local expansion. More recently, [17] developed an inversionalgorithm applicable to geometrically non-linear thin shells, including theeffects of contact and dry friction with an external body. Contribution.
In this paper we propose to use the Lagrangian formu-lation of [13] coupled with automatic code generation tools provided by theFEniCS Project finite element software [18] to compute the rest or unde-formed configuration of an object knowing the deformed configuration, theexternal loads and the material properties. We show experimental validationthat the methodology is effective at recovering the undeformed configurationfrom imaging data. The formulation requires only a few minor modificationsof the direct simulations, making it easy to implement. The automated dif-ferentiation tools from FEniCS Project provide a great deal of flexibility, forexample, permitting users to quickly and easily modify the material model4o suit their own problem.
Outline.
This paper is organized as follows; first, we give a descriptionof the finite strain elasticity formulation and the constitutive equations used.Next, we explain the inverse deformation analysis method. We test ourformulation on some simple analytical cases described in [19, 20]. Then,we show in some numerical examples how our variational formulation cansurpass the iterative algorithm proposed by [14]. Finally, we demonstratea relevant real-world application by retrieving the undeformed configurationof a Polydimethylsiloxone (PDMS) beam under the action of gravity fromimaging data. 5
ODELEXTERNAL FORCES M O D E L EQUILIBRIUM (1/2)INTERMEDIATERESULTRUN SIMULATIONEQUILIBRIUMSEGMENTATIONMEDICAL IMAGEEQUILIBRIUM IMAGEPROCESS EQUILIBRIUM (1/2)INTERMEDIATERESULTINITIAL CONFIGURATIONNUMERICALMODELINITIAL CONFIGURATIONNUMERICALMODEL ? RUN SIMULATIONMODEL EXTERNAL FORCES SURGERYENVIRONMENT RUN SIMULATIONEQUILIBRIUM (2/2)SIMULATIONRESULTEQUILIBRIUM (2/2)SIMULATIONRESULTRUN SIMULATION Figure 1: Our pipeline starts with the acquisition of a medical image, in which the organcan be segmented. The organ is observed at equilibrium under the effect of external forcesdue to its environment (e.g. gravity). In the usual pipeline (1, in red), the segmentedgeometry is considered as the initial geometry. Then, the external forces are applieduntil equilibrium to obtain the intermediate geometry used for simulating the procedure.We propose an alternative approach (2, in green), where we constrain the intermediategeometry to be identical to the segmented geometry. It involves the computation of anew geometry (represented by the question mark symbol), which is the organ geometrysuch that it would deform to the segmented/imaged configuration if external forces wereapplied. Here, the final result takes into account the undeformed geometry of the organ. . Finite strain elasticity formulation Consider a deformable body B . We denote the undeformed configurationΩ . The location of a particle of B in Ω is denoted X . Conversely, thedeformed configuration is noted Ω, and the location of a particle of B in Ω isnoted x . A one-to-one mapping φ maps the position of a particle X in Ω to the position of the same particle x in Ω, i.e x = φ ( X ). The configurationΩ can be obtained by φ (Ω ) = { φ ( X ) | X ∈ Ω } . These definitions aredepicted in figure 2. Figure 2: In a standard deformation analysis we compute the displacement vector u fromknowledge of the undeformed configuration Ω . In the inverse deformation analysis the goalis to compute the displacement vector u (cid:48) from knowledge of the deformed configurationΩ. Let us introduce the deformation gradient F that maps a line elementd X in Ω to a line element d x in Ω: d x = F · d X . (1)We can write the deformation gradient and the Jacobian as F = ∂ φ ∂ X , (2)7 = det F . (3)As φ is a one-to-one mapping, F is not singular and can be inverted,resulting in J (cid:54) = 0. The Jacobian maps a volume element dΩ in Ω to avolume element dΩ in Ω dΩ = J · dΩ . (4)For each point, we introduce the displacement u as the position differencebetween the deformed and the undeformed configuration u ( X ) = x − X = φ ( X ) − X (5)The deformation gradient can also be written as a function of the displace-ment such as F = ∂ φ ∂ X = ∂ u ∂ X + I = ∇ u + I , (6)where ∇ ( • ) is the gradient in Ω , with respect to the initial spatial position.The gradient in Ω, with respect to the deformed spatial position, is denoted ∇ ( • ). I is the usual second-order identity tensor. Similarly to the strainmeasure F , we introduce the right Cauchy-Green strain tensor C , its conju-gate the left Cauchy-Green strain tensor B and the Green-Lagrange straintensor E C = F T F , (7) B = F F T , (8) E = 12 ( C − I ) . (9)Hyperelastic material laws commonly use invariants of C and B to definetheir elastic energy I C = tr( C ) , (10)II C = 12 (cid:0) (tr( C )) − tr (cid:0) C (cid:1)(cid:1) , (11)III C = det C . (12)8 .2. Strong form At equilibrium in the deformed configuration, the balance of momentumcan be written as follows ∇ · σ + ρ b = , (13)where σ is the Cauchy stress tensor, ρ is the density of the material in thedeformed configuration and b are the external forces in the deformed config-uration. Equation 13 is called the strong form and is written in the deformedconfiguration Ω. To write the strong form in the initial configuration Ω , weintroduce the density of the material ρ in the undeformed configuration andthe first Piola-Kirchhoff stress tensor P ∇ · P + ρ b = , (14)where σ and P are related by the Piola transform σ = 1 J P F T . (15) The weak form is obtained by multiplying the strong form by test func-tions η and integrating over the whole domain. Equation 14 is written in theinitial configuration and leads to − (cid:90) Ω ( ∇ · P ) · η dΩ = (cid:90) Ω ( ρ b ) · η dΩ . (16)By using the divergence theorem we obtain − (cid:90) Ω ( ∇ · P ) · η dΩ = (cid:90) Ω P : ∇ η dΩ − (cid:90) ∂ Ω ( P · n ) · η d ∂ Ω , (17)where the colon operator : is the inner product between tensors, n is theoutward unit normal at the boundary and ∂ Ω the surface boundary of Ω .The quantity P · n is the traction boundary condition. We here assume thatit is prescribed on a part Γ of the boundary as P · n = t . On the remainingpart of the boundary, we assume that the value of the displacement is given,i.e. a Dirichlet condition. We then obtain the equilibrium in the referenceconfiguration (cid:90) Ω P : ∇ η dΩ = (cid:90) Ω ρ b · η dΩ + (cid:90) Γ t · η dΓ . (18)9ote that the boundary integral on the remaining part ∂ Ω \ Γ vanishesdue to the Dirichlet condition. By injecting equation 15 in the last equation18, we obtain the weak form in the deformed configuration (cid:90) Ω σ : ∇ η dΩ = (cid:90) Ω ρ b · η dΩ + (cid:90) Γ t · η dΓ . (19) For many materials, simple elastic models such as the St. Venant Kirch-hoff model are not sufficient to describe the observed behavior. More complexhyperelastic models provide a mechanism of modeling the stress-strain be-havior of complex materials such as elastomers or biological tissues.
Neo-Hookean
A neo-Hookean solid is a hyperelastic material model thatcan be used for predicting the nonlinear stress-strain behavior of materialsundergoing large deformations. Its strain energy density is defined as: ψ NH = µ B − − µ ln( J ) + λ J ) , (20)where λ and µ are material constants called the Lam´e parameters. Mooney-Rivlin
A Mooney–Rivlin solid is a hyperelastic material modelwhere the strain energy density function ψ MR is a linear combination of twomodified invariants of the left Cauchy–Green deformation tensor B . Rubber-like materials are often modeled using the Mooney–Rivlin model with strainenergy density ψ MR = C (cid:0) I B − (cid:1) + C (cid:0) II B − (cid:1) + D ( J − , (21)with the modified invariants I B = J − I B , II B = J − II B and where C , C , D are material constants. All material models previously introduced were intended for compressiblematerials, i.e. materials where the volume may change during deformation.Conversely, some materials such as living tissues or rubbers can be assumedto be nearly-incompressible or even completely incompressible, i.e. volumeis preserved during deformation J ∼ locking . Simply put, locking oc-curs when too many constraints are imposed on the discrete formulation andits overall approximation power is destroyed.To overcome these difficulties, mixed formulations have been developed.In these formulations, the variational principle is modified by writing thepotential energy functional similar to equation 21, except that the strainenergy is expressed in terms of the deviatoric component only and the in-compressibility constraint is explicitly enforced using a Lagrange multiplierwith physical meaning akin to pressure ( p ). It turns out that the Lagrangemultipliers can be expressed as a function of the hydrostatic pressure values f ( p ). It can be shown that ψ ( u , p ) = ψ ( u ) − f ( p ) , (22) f ( p ) := − σ hydro = −
13 tr( σ ) = − ∂ψ∂J . (23) Neo-Hookean
By calculating f ( p ), we can deduce the mixed displacement-pressure formulation of a nearly-incompressible Neo-Hookean material ψ NH ( u , p ) = µ B − − µ ln( J ) + p ln( J ) − λ p . (24) Mooney-Rivlin
By calculating f ( p ), we can deduce the mixed displacement-pressure formulation of a nearly-incompressible Mooney-Rivlin material ψ MR ( u, p ) = C (cid:0) I B − (cid:1) + C (cid:0) II B − (cid:1) + p ( J − − D p . (25)These nearly-incompressible energy densities are used to generate the FEn-iCS Project results in this paper. We use the FEniCS Project finite element software [18] to discretise boththe standard finite strain elasticity problem and the inverse finite strain elas-ticity problem that we will outline in the next section. We use a mixed11isplacement-pressure finite element formulation with second-order continu-ous Lagrangian finite elements for displacement u and first-order continuousLagrangian finite elements for pressure p . This pairing is well-known to beinf-sup stable and relatively robust with respect to numerical locking.The variational forms of the residual equations 18 and 29 are defined inthe Unified Form Language (UFL) [21] and symbolically differentiated toderive an expression for consistent Jacobian. The FEniCS Form Compiler(FFC) [22] is used to automatically generate low-level C ++ code from thehigh-level UFL description that can calculate the Jacobian and residual celltensors. The overall solution process is driven by the DOLFIN finite elementlibrary [23]. We use a standard Newton-Raphson algorithm with continuationin the loading parameter. The linear system within the Newton-Raphsonalgorithm is solved using the direct solver MUMPS via PETSc [24]. Thecomplete implementation of the standard or inverse problem is around 100lines of Python code that closely follows the mathematical structure of theproblem. We refer the reader to the supplementary material [25] for furtherdetails.
3. Inverse finite strain elasticity formulation
This section presents two methods to compute the undeformed configu-ration knowing the deformed configuration under known loading. We firstintroduce our methodology derived from [13], then we briefly outline a simpleiterative geometric algorithm described in [14].
In section 2, we introduced how to compute the deformed configurationof a body undergoing external forces. The inputs were the undeformed ge-ometry and the external forces, which means X , the rest-position was knownand x , the deformed position, was unknown. In this section, we introduceour method to compute the undeformed configuration of a body undergoingexternal forces. The inputs of the inverse deformation formulation are thedeformed geometry and the forces applied to the body. The most intuitiveapproach is to solve equation 18 or equation 19 for the unknown X . Thisapproach has the advantage of being based on classical mechanical princi-ples. However, mechanical quantities such as strains or stresses are defineddepending on X . This approach requires few straightforward modificationsto the equations in order to solve equation 18 or 19. In the inverse approach,12he initial geometry is replaced by the deformed geometry ( x ). We redefinethe displacement of equation 5 as u (cid:48) ( x ) = X − x . (26)Note that trivially u (cid:48) + u = . (27)This redefinition does not modify the classical finite element pipeline:the unknown position is still the first term in which the known position issubtracted. Notice the selection of the gradient compared to equation 6: wenow compute gradients in the deformed configuration and this necessitatesthe redefinition of the deformation gradient F = ∂ φ ( X ) ∂ X = (cid:18) ∂ X ∂ x (cid:19) − = (cid:18) ∂ u (cid:48) ( x ) ∂ x + I (cid:19) − = ( ∇ u (cid:48) + I ) − . (28)Henceforth, when performing an inverse deformation analysis, F and all de-rived quantities (strain measures, invariants, energy densities, stress measuresetc.) are always computed using the above redefinition in terms of u (cid:48) .The goal then is to solve equation 18 with x known and X unknownfor u (cid:48) . The weak equilibrium in the inverse deformation is expressed in thedeformed configuration (cid:90) Ω σ : ∇ η dΩ = (cid:90) Ω ρ b · η d Ω + (cid:90) Γ t · η dΓ . (29)We can notice three main differences compared to equation 16: (1) Theintegration domain is no longer the undeformed domain but the deformeddomain. (2) The gradient of the trial function is in the deformed configu-ration. (3) The external forces are written in the deformed configuration.This new formulation requires us to make one change compared to the directpipeline; rewrite F in terms of u (cid:48) . The computation of the gradient and theintegration domain in the deformed configuration in the inverse analysis isequivalent to the computation of the gradient and the integration domainin the undeformed configuration in the direct analysis. That is why thesechanges in the formulation do not require a significant modification of a codeto perform the inverse analysis. This formulation can find the undeformedconfiguration of an object, knowing only the deformed configuration and theapplied forces. The process is ”one-shot” based on the equation of continuummechanics. 13 .2. Iterative geometric algorithm [14] proposed an Iterative Geometric Algorithm (IGA, not to be confusedwith Isogeometric Analysis). The algorithm is simple to implement and onlyrequires an existing (standard) forward deformation solver. The algorithmstarts with an initial guess for the undeformed configuration (usually chosen,for lack of a better choice, the deformed one) and applies successive displace-ment fields to it until a convergence criterion is reached. The sequence ofdisplacement fields is obtained from the direct simulations of the current rest-configuration undergoing external forces. The shape of the object after thedirect simulation provides an error compared to the exact rest-configurationby measuring the distance to the initial configuration. An updated estimateof the undeformed configuration is calculated by correcting the previous guesswith the difference between the computed and deformed configuration. Thealgorithm stops when the error (computed using the l -norm) is below adefined threshold (cid:15) or a maximum number of iterations N B max has beenreached. The process is outlined in algorithm 1.
Algorithm 1:
Iterative geometric algorithm from [14]. X ← X ini run direct simulation 0 with X the initial configuration u ← x − X err ← error between x and X ini j ← while err > (cid:15) and j < N B max do X j ← X ( j − − u ( j − run direct simulation j with X ( j − the initial configuration u j ← x j − X j err ← error between x j and X ini j ← j + 1 end4. Numerical results The inverse deformation framework is very similar to the traditional directframework. To assess the numerical precision of the inverse method, we firstapply a series of tests to verify the soundness of the direct approach in whichan analytic solution is known. 14 .1.1. Shear deformation
Simple shear : Simple shear deformation is a popular benchmark test [19].The initial geometry is a unit cube with prescribed Dirichlet boundary condi-tions u = ( y · k, , T with y the y -coordinate and k a constant, as illustratedin figure 3. Figure 3: 2D plane cut of a simple shear deformation of a unit cube. An x -displacementof y · k is applied on the boundary. For simple shear deformation, the deformation gradient is equal to F = ∇ u + I = k
00 1 00 0 1 . (30)Now, let us consider a cube made of a Mooney-Rivlin material. By replacingthe deformation gradient F in the equation 21, we obtain the value of thestrain energy density function in the cube. Following [19] we can obtain theenergy density and the components of the Cauchy stress tensor σ ψ = k ( C + C ) , (31)15 = k (2 C + 4 C )3 ,σ = − k (4 C + 2 C )3 ,σ = k (2 C − C )3 ,σ = k (2 C + 2 C ) ,σ = σ = 0 . (32)The values of ψ and σ have been evaluated in our framework with several val-ues of k , degrees of discretization, and constitutive parameters. The relativeerror (by using the L -norm) in strain energy and Cauchy stress tensor, com-pared to the analytical values, shows the exactness of the direct deformationframework to machine precision (10 − magnitude error). Generalized shear
The generalized shear deformation test is similar tothe simple shear deformation [19]. The initial geometry is a unit cube withprescribed Dirichlet boundary conditions u = ( y · k, , T with y the y -coordinate and k a constant, as illustrated in figure 4. Figure 4: 2D plane cut of a generalized shear deformation of a unit square. An x -displacement of y · k is applied to the boundary of the cube. For generalized shear deformation, the deformation gradient is equal to F = ky
00 1 00 0 1 . (33)16n the same manner as in the simple shear deformation, we consider acube made of a Mooney-Rivlin material and can apply the same methods tofind the analytical strain energy density function ψ and the Cauchy stresstensor components σ ψ = (cid:90) [ C ( I −
3) + C ( II − (cid:90) k y ( C + C ) dy= 4 k ( C + C )3 , (34) σ = k (8 C + 16 C )9 ,σ = − k (16 C + 8 C )9 ,σ = k (8 C − C )9 ,σ = k (2 C + 2 C ) ,σ = σ = 0 . (35)We realize the same tests as the simple shear (different k values, mesh pre-cision, and mechanical parameters) and evaluate the identical quantities, ψ and σ values. We observed an impact of the mesh on the strain energy andthe Cauchy stress. The error quickly decreases on mesh refinement to reachrelative errors under 2%. This section presents a series of tests to verify the consistency of ourinverse method with the direct approach. More precisely, we show that theundeformed configuration corresponds to the initial configuration used to de-form it. During these tests, we also compare our method to the IGA methodpresented in section 3.2 and evaluate their performance and convergencerates.
This test is based on the direct shear deformation verification performedin section 4.1.1. We verify that the inverse deformation of the simple shear17nd the generalized shear is consistent with the direct finite element analysis.The idea is to start the test with the deformed configuration and apply theinverse deformation to verify that the rest-configuration corresponds to theinitial geometry of the direct deformation. Since both shear deformations areentirely determined by a displacement field, the inverse deformation consistsof applying the opposite displacement field. It is then trivial to claim that thegeometry will be recovered, i.e. a unit cube. However, this test also verifiesthe deformation gradient, the strain energy, and stress tensors are sound.As explained previously, those measures should be equal in both inverse anddirect deformation. We verify these statements numerically in these tests.
Inverse simple shear : For the inverse simple shear deformation, the ma-terial points are now shifted by − k · y on the x -axis while the bottom is fixed( y = 0). As illustrated in figure 5. Figure 5: 2D plane cut of an inverse simple shear deformation of a unit cube. An x -displacement of − y · k is applied on the boundary. We calculate the deformation gradient F which is equal to the deforma-tion gradient in equation 30, as expected ∇ u (cid:48) + I = − k
00 1 00 0 1 , (36)18 = ( ∇ u (cid:48) + I ) − = k
00 1 00 0 1 . (37)Therefore, the strain energy, which is usually defined depending on F , isequal to the strain energy in equation 31, and the stress tensor of equation32 remains valid. Since the deformation is homogeneous (constant deforma-tion gradient), our quadratic finite element method is able to reproduce theanalytical solution down to machine precision. Inverse generalized shear : Similarly, the inverse version of the generalizedshear deformation leads to the same deformation gradient tensor (equation33), then to the same strain energy density function (equation 34). Therelative error is evaluated with different discretizations of the initial meshbut the same parameters set and we obtain with high precision the initialgeometry.
Part I : Let us consider a mesh with a single unit tetrahedron with a linearLagrangian finite element space. Its domain is denoted Ω T . The nodal coor-dinates are [0 , , T , [1 , , T , [0 , , T and [0 , , T . The nodes with y = 0are fixed, leaving only one free node. A uniform force f is applied along the y -axis. The tetrahedron is deformed so that the free node moves along the y -axis.In a first step, we compute the deformation φ with the direct method.A displacement u is computed for the free node. The deformed domain isΩ T = φ (Ω T ). In a second step, the initial geometry is the deformed geometryΩ T , i.e. a unit tetrahedron with the nodes y = 0 fixed, and the remainingnode displaced from u . The same uniform force f (cid:48) = f is applied. Aninverse simulation is computed so that the displacement of the free node is u (cid:48) . This example is depicted in figure 6a. Part II : We consider the same unit tetrahedron, with the same boundaryconditions. A uniform force f (cid:48) is applied along the y -axis.In the first step, an inverse simulation is computed, leading to a displace-ment of u (cid:48) . In the second step, the resulting geometry is deformed with adirect simulation leading to a displacement of u . This part of the exampleis depicted in 6b. 19 igure 6: Single tetrahedron test, plane view. The difference with the first part of the test is the order of the successivesimulations. In part I , the inverse simulation is performed after the directsimulation. In part II , it is the opposite. In both parts of the test, the goalis to verify that the following relationship: u (cid:48) = − u .Furthermore, the inverse simulation is computed with IGA to comparethe results and performance with our method. In this test, the error measureis defined as: (cid:107) u (cid:48) + u (cid:107) l . We measured this error with different constitutiveequations and varying their associated mechanical parameters. In total, weperformed 153 tests and provided a statistical analysis in table 1.We observe that the accuracy of the iterative algorithm depends on thenumber of iterations, but it also increases the computational cost becauseeach iteration calls a direct simulation. Our method provides high accuracywhile requiring only the solution of a problem with similar complexity to asingle iteration of IGA. Beyond the numerical results, one point is that in7 tests over the 153 of the part II , the iterative algorithm was not able toreach the accuracy of our method within 50 iterations.20art I Part IIPB IGA (1) IGA (2) PB IGA (1) IGA (2)average error 4.49E-12 2.12E-6 2.28E-12 5.22E-12 1.98E-6 2.28E-12SD 1.07E-11 1.09E-6 6.69E-12 1.25E-11 1.05E-6 2.41E-12minimum 9.26E-22 5.44E-8 4.15E-35 1.04E-21 5.47E-8 6.76E-12maximum 5.52E-11 4.48E-6 5.11E-11 7.26E-11 3.99E-6 5.04E-11avg Table 1: Benchmark results on a single tetrahedron simulation. We compare our physics-based method (PB), with the iterative geometric algorithm (IGA (1)) with the arbitraryconvergence criterion 10 − , and with the IGA at the same accuracy than PB (IGA (2)).
5. Experimental results
In this section, we will demonstrate that our inverse simulation methodcan match the outcome of a real experiment and therefore has value as apredictive modelling tool.We fixed one extremity of a beam made from Polydimethylsiloxane (PDMS)to a vertical support and allowed it to deform under gravity as shown in figure7). To extract the mesh of the deformed configuration from the image, weused the software Blender and contoured the beam on 2D images by hand,as shown in figure 8. This mesh will be called the ”reference” and used asground-truth for this section.To run the inverse deformation algorithm, we need three input parame-ters: the applied force field, the deformed configuration, and the mechanicalproperties. In this section, the force field is gravity and the deformed con-figuration was obtained by manual processing. A separate experiment wasperformed to obtain the mechanical properties and will be detailed in thefollowing section. igure 7: Experimental set-up: Initially straight PDMS beam clamped on the left sideand deformed by gravity. igure 8: Manual process in Blender to delineate the contours and extract the 3D meshof the deformed configuration. We used a PDMS (Sylgard 184, Ellsworth Adhesives) cylinder of density965 kg / m of undeformed dimensions 182 mm and 8 . ± ™ mechanicaltesting system (Biomomentum, Canada) as a testing rig for the unconfinedcompression tests. We used the following protocol:23 A 1 . µ N was used tomeasure the vertical force. • The vertical displacement was measured by the moving stage of the rigwith a resolution of 0 . µ N. • To minimize friction, paraffin oil was used between the sample and thecompression platens. • One loading cycle was executed on each specimen. To detect the re-sponse of the material at large strains, the samples were compressedat a constant speed of 0 .
083 mm / s until a displacement correspondingto 30% of the measured height was achieved. Particular attention wasused to monitor the samples that had uniformly expanded in the radialdirection and that their upper and lower faces remained adhered to themoving platen and the fixed platform for the entire duration of the test. • The Abaqus evaluation routine was used to fit the true stress - truestrain experimental curves with a Mooney-Rivlin model. Abaqus em-ploys a linear least-squares fit for the Mooney-Rivlin form to find theoptimal model parameters.In our case the optimal parameters are: D = 7 . × − Pa, C = 101709 .
668 Pa, C = 151065 .
460 Pa.24 igure 9: Biomomentum Mach-1 ™ mechanical testing system used for characterizing me-chanical properties of the PDMS. To verify the mechanical properties we compare the output of three dif-ferent simulation softwares all using an incompressible Mooney-Rivlin modeland boundary conditions imitating the setup shown in figure 11.
FEniCS:
We used the same model as described in section 2.4.2.
Abaqus:
We use a static step with a gravity load to solve the beam defor-mation in Abaqus. Abaqus/Standard uses Newton’s method as a numeri-cal technique for solving the nonlinear equilibrium equations. We employedC3D8RH elements, an 8-node linear brick, hybrid/mixed, constant pressure,reduced integration with hourglass control. The hybrid/mixed formulationis needed because of the material’s near-incompressibility.
SOFA:
We employed the
Multiplicative Jacobian Energy Decomposition method(MJED) which is an optimized algorithm for building the stiffness and tan-gent stiffness matrices of non-linear hyperelastic materials [27]. An MJED25mplementation is available in SOFA [28] for finite element formulation us-ing linear tetrahedral elements. The linear system of equations was solvedin every step of quasi-static simulation using a fast in-house linear equationsolver based on the Cholesky decomposition.For each model, we perform a mesh convergence analysis shown in figure10 where we plot the maximum deformation of the beam (located at the tip)for different mesh resolutions.
Number of points M a x i m u m d e f o r m a t i o n o f t h e t i p [ mm ] ExperimentAbaqusFEniCSSOFA
Figure 10: Mesh convergence analysis of the forward simulation. We calculate the maxi-mum deformation of the tip of beam for several level of refinement of the mesh.
We observe that the tip displacement for the three software converge tosimilar solutions (FEniCS: 132 .
52 mm, Abaqus: 132 .
71 mm, SOFA: 130 . .
68 mm. We observe a small differ-ence between the numerical solutions and the experiment.FEniCS and Abaqus give similar results while SOFA is 2 mm off. Weobserve in figure 10 that FEniCS and Abaqus converged with 60,000 pointswhile SOFA is still not converged with 160 ,
000 points. One reason is thatSOFA is usually designed for real-time simulation and only uses dynamicsolvers which can lead to inaccuracy compared with static solvers from FEn-iCS and Abaqus. Furthermore, the differences between numerical solutionscan be explained by the use of three slightly different formulations of theMooney-Rivlin law as well as different solvers for solving the equation.26 igure 11: 3D plot of the forward simulation. From top to down, in magenta: the experi-mental data, in blue: the SOFA simulation, in wire-frame green: the FEniCS simulationand in red: the Abaqus simulation.
Some factors can explain the difference between the numerical solutionsand the experimental value. For instance, the variation may be explained byinadequate constitutive equations or boundary conditions. Then, uncertain-ties in the mechanical properties measures may also be a factor, especiallybecause the PDMS might exhibit slightly asymmetric behaviour under com-pression and tension. Finally, we obtained the reference mesh of the unde-formed configuration manually based on 2D imaging data where inaccuraciescan be introduced.
In the previous section, we compared the forward simulations of threedifferent software with our experimental solution. In this section, we wantto verify the possibility of retrieving the undeformed configuration of our ex-perimental solution knowing only the surface of the deformed configuration,the known applied loads and the material properties.For this, we converted our experimental surface mesh of the deformedconfiguration into a volumetric mesh and applied our inverse deformation27lgorithm implemented using FEniCS. We previously showed a deformationdifference of 4 .
84 mm for the forward simulation in FEniCS. Of course, donot expect to obtain a perfectly straight beam (the ideal undeformed config-uration), but rather an error on the same order as in the forward simulation.
Figure 12: 3D plot of the inverse simulation. From down to top, in magenta: the exper-imental data we wish to retrieve the undeformed configuration, in wire-frame black: thetheoretical straight beam , in yellow: the result of the FEniCS inverse simulation.
We show in figure 12 the result of the inverse deformation algorithm. Asexpected, the inverse simulation (in yellow) applied to the experimental data(deformed configuration in magenta) is slightly different from the theoreti-cal straight beam that we should obtain (in black). To be more precise, weachieve an error of 5 .
36 mm compared with the idealised straight beam. Asmentioned previously, we expect an error on the order of that for the standarddeformation problem (4 .
84 mm) due to the inherent parametric and model-ing uncertainties (material model, material properties, boundary conditions,geometry) already discussed. We therefore judge that the proposed method-ology has strong potential for prediction of the undeformed configuration ofa soft body. 28 . Conclusions
In the present paper we performed a numerical and experimental studyof the inverse deformation problem .Our study used the Lagrangian formulation of [13] as a basis for imple-menting the inverse algorithm in the FEniCS Project finite element software.We took advantage of the automatic differentiation and code generation capa-bilities to bypass the difficulties of deriving and implementing the consistentJacobian. The user must then supply the deformed configuration, the me-chanical properties, boundary conditions and the applied forces. The usercan easily modify the input mesh, run the code efficiently in parallel, changethe constitutive model or change the boundary conditions according to theirneeds. We have made the code and data available in the supplementarymaterial.We applied the approach to simple academic examples where we con-sidered two different incompressible hyperelastic models (neo-Hookean andMooney-Rivlin) and different boundary conditions. We demonstrated on asimple test case that our method is more efficient in terms of robustness andaccuracy than the IGA method of [14]. We have only compared with the clas-sical IGA method of Sellier but other works like [29] have improved on thisalgorithm. However, we can say that unless an iterative approach requiresonly one forward model solution, in most circumstances the mechanics-basedapproach detailed here is likely to be faster and more robust.Finally we applied the method to an experiment with a PDMS beamdeformed under gravity. We verified and quantified the performance of thedirect simulations of three different widely-used software (Abaqus, FEniCS,SOFA). Using the inverse deformation algorithm we achieve an error of 5 . Supplementary material
The reference [25] (doi:10.6084/m9.figshare.14035793) contains a full im-plementations of the forward and inverse deformation problems using theFEniCS Project finite element software. The latest version is also availableon GitHub at https://github.com/Ziemnono/fenics-inverseFEM
References [1] M. L. Raghavan, B. Ma, M. F. Fillinger, Non-invasive determinationof zero-pressure geometry of arterial aneurysms, Annals of BiomedicalEngineering 34 (2006) 1414–1419. doi:10.1007/s10439-006-9115-7 .[2] J. Lu, X. Zhou, M. L. Raghavan, Inverse elastostatic stress analysisin pre-deformed biological structures : Demonstration using abdominalaortic aneurysms, Journal of Biomechanics 40 (2007) 693–696. doi:10.1016/j.jbiomech.2006.01.015 .[3] A. Mˆıra, A. K. Carton, S. Muller, Y. Payan, A biomechanical breastmodel evaluated with respect to MRI data collected in three differentpositions, Clinical Biomechanics 60 (2018) 191–199. arXiv:1811.10221 , doi:10.1016/j.clinbiomech.2018.10.020 .[4] M. Koishi, S. Govindjee, Inverse design methodology of a tire, TireScience and Technology 29 (2001) 155–170. doi:10.2346/1.2135236 .[5] V. D. Fachinotti, A. Cardona, P. Jetteur, Finite element modelling of in-verse design problems in large deformations anisotropic hyperelasticity,International Journal For Numerical Methods In Engineering 74 (2008)894–910. doi:10.1002/nme.2193 .[6] J. Adkins, A reciprocal plane property of the finite plan strain equations,Journal of the Mechanics Physics of Solid 6 (1958) 267–275. doi:10.1016/0022-5096(58)90002-4 .[7] R. T. Schield, Inverse deformation results in finite elasticity, Zeitschriftf¨ur angewandte Mathematik und Physik ZAMP 18 (1967) 490–500. doi:10.1007/BF01601719 . 308] D. E. Carlson, T. Shield, Inverse deformation results for elastic materi-als, Zeitschrift f¨ur angewandte Mathematik und Physik ZAMP 20 (1969)261–263. doi:10.1007/BF01595564 .[9] M. M. Carroll, F. J. Rooney, Implications of Shield ’ s inverse defor-mation theorem for compressible finite elasticity, Zeitschrift f¨ur ange-wandte Mathematik und Physik ZAMP 56 (2005) 1048–1060. doi:10.1007/s00033-005-2023-0 .[10] S. Govindjee, P. A. Mihalic, Computational methods for inverse finiteelastostatics, Computer Methods in Applied Mechanics and Engineering136 (1996) 47–57. doi:10.1016/0045-7825(96)01045-6 .[11] S. Govindjee, P. A. Mihalic, Computational methods for in-verse deformations in quasi-incompressible finite elasticity, Interna-tional Journal For Numerical Methods In Engineering 43 (1998)821–838. doi:10.1002/(SICI)1097-0207(19981115)43:5<821::AID-NME453>3.0.CO;2-C .[12] T. Yamada, Finite element procedure of initial shape determination forhyperelasticity, Structural Engineering and Mechanics 6 (1998) 173–183. doi:10.12989/sem.1998.6.2.173 .[13] A. Albanesi, V. Fachinotti, A. Cardona, Design of compliant mech-anisms that exactly fit a desired shape, Mec´anica Computacional 28(2009) 3191–3205.[14] M. Sellier, An iterative method for the inverse elasto-static problem,Journal of Fluids and Structures 27 (2011) 1461–1470. doi:10.1016/j.jfluidstructs.2011.08.002 .[15] J. Bols, J. Degroote, B. Trachet, B. Verhegghe, P. Segers, J. Vieren-deels, A computational method to assess the in vivo stresses and un-loaded configuration of patient-specific blood vessels, Journal of Com-putational and Applied Mathematics 246 (2013) 10–17. doi:10.1016/j.cam.2012.10.034 .[16] X. Chen, C. Zheng, W. Xu, K. Zhou, An asymptotic numerical methodfor inverse elastic shape design, ACM Transactions on Graphics 33(2014). doi:10.1145/2601097.2601189 .3117] M. Ly, R. Casati, F. Bertails-Descoubes, M. Skouras, L. Boissieux, In-verse elastic shell design with contact and friction, SIGGRAPH Asia2018 Technical Papers, SIGGRAPH Asia 2018 37 (2018). doi:10.1145/3272127.3275036 .[18] M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,C. Richardson, J. Ring, M. E. Rognes, G. N. Wells, The FEniCSproject version 1.5, Archive of Numerical Software 3 (2015). doi:10.11588/ans.2015.100.20553 .[19] L. A. Mihai, A. Goriely, Numerical simulation of shear and the Poyntingeffects by the finite element method: An application of the generalisedempirical inequalities in non-linear elasticity, International Journal ofNon-Linear Mechanics 49 (2013) 1–14. doi:10.1016/j.ijnonlinmec.2012.09.001 .[20] C. K. Lee, L. Angela Mihai, J. S. Hale, P. Kerfriden, S. P. Bor-das, Strain smoothing for compressible and nearly-incompressible fi-nite elasticity, Computers and Structures 182 (2017) 540–555. doi:10.1016/j.compstruc.2016.05.004 .[21] M. S. Alnæs, A. Logg, K. B. Ølgaard, M. E. Rognes, G. N. Wells, Unifiedform language: A domain-specific language for weak formulations ofpartial differential equations, ACM Trans. Math. Softw. 40 (2014) 9:1–9:37. doi:10.1145/2566630 .[22] A. Logg, K. B. Ølgaard, M. E. Rognes, G. N. Wells, FFC: the FEniCSform compiler, in: A. Logg, K.-A. Mardal, G. Wells (Eds.), AutomatedSolution of Differential Equations by the Finite Element Method, Lec-ture Notes in Computational Science and Engineering, Springer BerlinHeidelberg, 2012, pp. 227–238.[23] A. Logg, G. N. Wells, DOLFIN: Automated finite element comput-ing, ACM Trans. Math. Softw. 37 (2010) 20:1–20:28. doi:10.1145/1731022.1731030 .[24] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschel-man, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, D. Karpeyev,D. Kaushik, M. G. Knepley, D. A. May, L. C. McInnes, R. T. Mills,32. Munson, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang,H. Zhang, PETSc Web page, (2019).[25] A. Mazier, A. Bilger, A. E. Forte, I. Peterlik, J. S. Hale, S. P. Bordas,Supplementary material for inverse deformation analysis: an experimen-tal and numerical assessment using the FEniCS project, (2021).[26] A. E. Forte, S. Galvan, F. Manieri, F. Rodriguez y Baena, D. Dini, Acomposite hydrogel for brain tissue phantoms, Materials and Design 112(2016) 227–238. doi:10.1016/j.matdes.2016.09.063 .[27] S. Marchesseau, T. Heimann, S. Chatelin, R. Willinger, H. Delingette,Fast porous visco-hyperelastic soft tissue model for surgery simulation:Application to liver surgery, Progress in Biophysics and Molecular Bi-ology 103 (2010) 185–196, special Issue on Biomechanical Modelling ofSoft Tissue Motion. doi:10.1016/j.pbiomolbio.2010.09.005 .[28] F. Faure, C. Duriez, H. Delingette, J. Allard, B. Gilles, S. Marchesseau,H. Talbot, H. Courtecuisse, G. Bousquet, I. Peterlik, et al., Sofa: Amulti-model framework for interactive physical simulation, in: Soft tis-sue biomechanical modeling for computer assisted surgery, Springer,2012, pp. 283–321. doi:10.1007/8415\_2012\_125 .[29] M. K. Rausch, M. Genet, J. D. Humphrey, An augmented iterativemethod for identifying a stress-free reference configuration in image-based biomechanical modeling, Journal of Biomechanics 58 (2017) 227–231. doi:10.1016/j.jbiomech.2017.04.021doi:10.1016/j.jbiomech.2017.04.021