S Niroomandi

I Alfaro

E Cueto

F Chinesta

Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models

Keywords: Real-time simulation, model reduction, large deformations, haptic devices, Kirchhoff-Saint Venant

Model reduction techniques have shown to constitute a valuable tool for real-time simulation in surgical environments and other fields. However, some limitations, imposed by real-time constraints, have not yet been overcome. One of such limitations is the severe limitation in time (established in 500 Hz of frequency for the resolution) that precludes the employ of Newtonlike schemes for solving non-linear models as the ones usually employed for modeling biological tissues. In this work we present a technique able to deal with geometrically non-linear models, based on the employ of model reduction techniques, together with an efficient non-linear solver. Examples of the performance of the technique over some examples will be given.

Introduction

Minimally invasive surgery (MIS) gets more and more common nowadays. A minimally invasive procedure typically involves use of laparoscopic devices and remote-control manipulation of instruments with indirect observation of the surgical field through an endoscope or similar device, and are carried out through the skin or through a body cavity or anatomical opening. This may result in less pain, less strain of the organism, small injuries (aesthetic reasons), economic gain because of shorter hospital stays, etc. On the other hand, there exist some important difficulties for the surgeon as a result of his restricted vision of the organs, difficult handling of the instruments, very restricted mobility, difficult hand-eye coordination and no tactile perception. Therefore, an important training phase is required before a surgeon acquires the skills necessary to adequately perform minimally invasive surgery. Currently, surgeons are trained to perform minimally invasive surgery by using mechanical simulators -just like plane pilots have been trained for many years-or living animals. There is some consensus on the limited realism of the mentioned simulators due to the complexity of the simulations to be carried [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF]. The latter training method consists in practicing simple or complex surgical procedures on living animals (often pigs for abdominal surgery).

Because of the limitations of current training methods, there is a large interest in developing surgery simulation software, possibly with haptic feedback, for providing efficient and quantitative gesture training systems [START_REF] Delingette | Towards realistic soft tissue modeling in medical simulation[END_REF]. A simulator consists of several parts: input devices, the core of the simulator and output devices [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF]. The input devices consist of a force-feedback (haptic) device with which the user moves the virtual surgery tools, a mouse and a keyboard to have a good view of the operation zone. The core of the simulator does the geometrical and the physical modeling. Geometrical models are obtained using medical imagery like CT-scan and MRI after they have been converted to standard 3D graphics formats. Physical models depend on the specific soft tissue. Many existing, commercial, surgery simulators employ linear elastic models or even mass-spring models. The core of the simulator also detects collision between virtual tools and the soft tissue and then calculates the corresponding reaction forces on the tools. The computed data is sent to the output devices which consist of a screen and the haptic device. The haptic device transforms reaction forces and moments to the user and makes him/her have the perception of virtual contact with the tissue. There are some requirements for visual and haptic feedback that make simulations difficult to perform in real time. A surgery simulator must provide a realistic visual presentation of the surgical procedure. Visual feedback is specially important in video-surgery because it helps the surgeon to acquire a tridimensional perception of the environment. If the positions, orientations and deformations of objects on the screen are updated at a rate less than about 30 times per second, users will no longer perceive the simulation as continuous motion. Haptic feedback provides the sensation of the movement to the user and therefore it significantly enhances his surgical performance. But it should have a frequency between 300 and 1000 Hz that is very difficult to achieve for complex tissue behavior.

In order to comply with these requirements different methods have been used during last years often with trade-off between time and accuracy. Among them are mass-spring methods, surface models and finite element or boundary element methods. Spring-mass models consist of points with mass that are linked by springs and dashpots. They have been used extensively for simulating soft tissues [START_REF] Waters | A physical model of facial tissue and muscle articulation derived from computer tomography data[END_REF][START_REF] Kenedi | Tissue mechanics[END_REF][START_REF] Delingette | A craniofacial surgery simulation testbed[END_REF]. The main advantage of spring models is their ease of implementation and they have been used for statics and dynamics analysis. Another advantage is their ability to model cutting and suturing simply by removing or adding connections between vertices. However, there are some disadvantages with these models. The most important is that, in general, they do not reproduce even the laws of linear Elasticity [START_REF] Meier | Realtime deformable models for surgery simulation: a survey[END_REF].

The models of the last category are based on the laws of continuum mechanics. Among the available numerical methods for solving these models, finite element and boundary element methods can be cited [START_REF] Monserrat | A new approach for the real-time tissue deformations in surgery simulation[END_REF]. A good overview of different soft tissue models based on continuum mechanical is given in [START_REF] Famaey | Soft tissue modeling for applications in virtual surgery and surgical robotics[END_REF]. Very recently, a system based upon the use of neural networks has been presented [START_REF] Deo | Phyness: A physics-driven neural networks-based surgery simulation system with force feedback[END_REF]. In it, the system is trained with a large set of possible load states in order to achieve real-time performance in the execution loop.

Very recently, General Purpose Graphic Processing Units (GPGPUs) have been incorporated into the field of surgical simulation. This graphic hardware is able to perform very fast operations, and a complete nonlinear, explicit, finite element code has been implemented in that work with a gain on speed on the order of 20 times [START_REF] Taylor | High-speed nonlinear finite element analysis for surgical simulation using graphics processing units[END_REF]. This approach has allowed to implement medium-sized models (16000 tetrahedra) including neo-hookean behavior in an explicit total Lagrangian approach. The main drawback of this elegant approach to the problem is the conditional stability of explicit integration procedures and also the explicit evaluation of the constitutive laws.

Other methods like mesh free methods have been recently used for realtime simulation of soft tissues. Lim and De applied the point collocationbased method of finite spheres (PCMFS) technique to simulate tissue deformations that are geometrically nonlinear [START_REF] Lim | Real time simulation of nonlinear tissue response in virtual surgery using the point collocation-based method of finite spheres[END_REF]. The technique is based on a novel combination of a multi-resolution approach coupled with a fast analysis scheme in the nonlinear model.

Capel et al. [START_REF] Capell | Interactive skeleton-driven dynamic deformations[END_REF] computed the equations of motion for a linearly elastic material using FEM and then by using coarse volumetric meshes to do simulations at interactive rates. Krysl et al. applied model order reduc-tion approach for dimensional reduction of nonlinear finite element models of solid dynamics [START_REF] Krysl | Dimensional model reduction in non-linear finite element dynamics of solids and structures[END_REF]. Although they proposed a general approach the computations cannot be performed in real time. Barbic and James, however, exploited dimensional model reduction techniques to build reduced-coordinate deformable models for material objects with complex geometry using St. Venant-Kirchhoff constitutive laws. In this way they have performed simulations at real-time update rates of geometrically-nonlinear materials [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF]. Unfortunately the method seems not to be generalized to more complex material behaviors like neo-Hookean constitutive laws, etc.

In a previous paper the authors introduced a new technique for the realtime simulation of non-linear tissue behavior based on a model reduction technique known as Proper Orthogonal (POD) or Karhunen-Loève decompositions. The technique is based upon the construction of a complete model (using Finite Element modeling or other numerical techniques, but data could also be extracted from experimental results, if available) and the extraction and storage of the most relevant information in order to construct a model with very few degrees of freedom, but that takes into account the highly non-linear response of most living tissues. It was applied to the simulation of palpation a human cornea composed of a nonlinear, hyperelastic matrix reinforced with collagen fibers [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF]. But the reduced model is actually linear although the best set of basis vectors, in statistical sense, is employed, if no updating of the tangent stiffness matrix is performed, as in a standard Newton-Raphson procedure. So the results had some error, more important the higher the strain is, in comparison with a standard FEM solution.

In this paper we apply an advanced non-linear solver, the asymptotic numerical method (ANM) to simulate the nonlinear behavior of soft tissues in a neighborhood of a known equilibrium state. In this method complex equilibrium paths are sought in the form of asymptotic expansions, and they are determined by solving several linear problems with a single tangent stiffness matrix [START_REF] Cochelin | The asymptotic-numerical method: an efficient perturbation technique for non-linear structural mechanics[END_REF]. Then we apply model reduction techniques (in this work we have used again POD) to the obtained system of equations. These reduced models are obtained for some "most probable" load states and then are properly interpolated for any other load state. In this way the system can be solved in very short times. The organization of the paper is as following: the fundamentals of model reduction techniques are explained in section 2. In section 3 first the asymptotic numerical method is introduced and then the formulation for geometrically-nonlinear problems is obtained. In section 4 it is shown how ANM formulation in reduced coordinates of POD can be obtained and in section 5 several examples are given to show the performance of the method.

Model reduction techniques

The essential idea of model reduction techniques in the context of Continuum Mechanics is to extract the most relevant information of a given simulation, that is performed off-line and stored for this particular application, in order to construct a good-quality Ritz-like basis to perform reduced-model (on-line) simulations with very few degrees of freedom [START_REF] Ryckelynck | On the a priori Model Reduction: Overview and recent developments[END_REF] [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF]. These basis functions are global and, in statistical terms, of very good quality. This is in sharp contrast with the finite element method, for instance, that employs general-purpose, locally supported, piecewise polynomial shape functions.

Fundamentals: Karhunen-Loève or Proper Orthogonal decomposition

The technique we employed is known by a wide variety of names, since it has been employed and re-discovered in many branches of science and engineering. Maybe the most common names are Karhunen-Loève decomposition [START_REF] Karhunen | Uber lineare methoden in der wahrscheinlichkeitsrechnung[END_REF] [19], proper orthogonal decomposition (POD) [START_REF] Sirovich | Turbulence and the dynamics of coherent structures part I: coherent structures[END_REF] or empirical orthogonal functions [START_REF] Lorenz | Empirical Orthogonal Functions and Statistical Weather Prediction[END_REF].

We assume that the evolution of a certain field u(x, t), governed by a PDE, is known. In practical applications, this field is expressed in a discrete form, that is, it is known typically at the nodes of a spatial mesh and for some time steps of existing simulations u(x i , t n ) ≡ u n i . The same can be wrote introducing a time discretization u n (x) ≡ u(x, t = nΔt); ∀n ∈ [1, . . . , P ]. The main idea of the Karhunen-Loève (K-L) decomposition is how to obtain the most typical or characteristic structure φ(x) among these u n (x) ∀n. This is equivalent to obtaining a function φ(x) that maximizes the functional λ defined by

λ = n=P n=1 i=N i=1 φ(x i )u n (x i ) 2 i=N i=1 (φ(x i )) 2 (1)
The maximization (δλ = 0) leads to:

n=P n=1 i=N i=1 φ(x i )u n (x i ) j=N j=1 φ(x j )u n (x j ) = λ i=N i=1 φ(x i )φ(x i ); ∀ φ (2)
which can be rewritten in the form

i=N i=1 j=N j=1 n=P n=1 u n (x i )u n (x j )φ(x j ) φ(x i ) = λ i=N i=1 φ(x i )φ(x i ); ∀ φ (3)
Defining the vector φ such that its i-th component is φ(x i ), Eq. ( 3) takes the following matrix form

φT c φ = λ φT φ; ∀ φ ⇒ c φ = λφ (4)
where the two-point correlation matrix is given by

c ij = n=P n=1 u n (x i )u n (x j ) ⇔ c = n=P n=1 u n (u n ) T (5)
which is symmetric and positive definite. If we define the matrix Q containing the discrete field history:

Q = ⎛ ⎜ ⎜ ⎜ ⎝ u 1 1 u 2 1 . . . u P 1 u 1 2 u 2 2 . . . u P 2 . . . . . . . . . . . . u 1 N u 2 N . . . u P N ⎞ ⎟ ⎟ ⎟ ⎠ (6)
it is straightforward to verify that the matrix c in Eq. ( 4) results

c = Q Q T (7)
where the diagonal components are given by

c ii = (Q Q T ) ii = j=P j=1 (u j i ) 2 . ( 8
)
Thus, the functions defining the most characteristic structure of u n (x) are the eigenfunctions φ k (x) ≡ φ k associated with the highest eigenvalues.

A posteriori reduced modeling

If some direct, off-line, simulations have been carried out and their results stored in memory, we can determine u(

x i , t n ) ≡ u n i , ∀i ∈ [1, . . . , N],
∀n ∈ [1, . . . , P ], and from these the r eigenvectors related to the r-highest

eigenvalues φ k = φ k (x i ), ∀i ∈ [1, . . . , N], ∀k ∈ [1, . . . , r] (with r N
). Now, we can try to use these r eigenfunctions for approximating the solution of a problem slightly different to the one that has served to define u(x i , t n ). For this purpose we need to define the matrix

A A = ⎛ ⎜ ⎜ ⎜ ⎝ φ 1 (x 1 ) φ 2 (x 1 ) . . . φ r (x 1 ) φ 1 (x 2 ) φ 2 (x 2 ) . . . φ r (x 2 ) . . . . . . . . . . . . φ 1 (x N ) φ 2 (x N ) . . . φ r (x N ) ⎞ ⎟ ⎟ ⎟ ⎠ . (9)
Now, we consider the linear system of equations resulting from the discretization of a partial differential equation (PDE) in the form

K U = F . ( 10
)
Obviously, in the case of evolution problems F contains the contribution of the solution at the previous time step. Then, assuming that the unknown vector contains the nodal degrees of freedom, it can be expressed as

U = i=r i=1 ζ i φ i = A ζ, (11)
the global system of equations now reads

K U = F ⇒ K A ζ = F , (12)
and multiplying both terms by A T it results

A T K A ζ = A T F , (13)
which proves that the final system of equations is of low order, i.e. the dimensions of A T K A are r × r, with r N, and the dimensions of both ζ and A T F are r × 1.

Note, on the other hand, that non-homogeneous boundary conditions (like rigid-body motions, for instance) do not impose any special difficulties to the method. A simple change of variables can deal with this, as proved in [START_REF] González | Recent advances on the use of separated representations[END_REF].

Asymptotic numerical methods

The evolution of a non-linear solid mechanics problem is usually done with the help of iterative methods based upon the Newton-Raphson technique. This iterative structure of the method needs for an update of the tangent stiffness matrix at each iteration. Modified Newton methods, which employ a constant tangent matrix in order to save computation time are also very popular, but fail notably in the neighborhood of buckling points, for instance.

Nowadays different alternatives exist. Potier-Ferry and coworkers have developed in the last years a family of methods based upon asymptotic expansions that are able to accurately trace the non-linear solution of the problem with a minimum number of matrix updates, see for instance [START_REF] Cochelin | The asymptotic-numerical method: an efficient perturbation technique for non-linear structural mechanics[END_REF][START_REF] Zahrouni | Asymptotic numerical method for nonlinear constitutive laws[END_REF]. In this technique, the nodal displacement vector U and the load parameter λ (here it plays the role of a pseudo-time, ranging from 0 to 1) are represented by power series expansions with respect to a control parameter a. Introducing these expansions into the weak form of the problem a sequence of linear problems, that share a unique tangent stiffness matrix, is obtained. Thus, it seems natural to combine this asymptotic technique with a Ritz-like technique to obtain a reduced model that can accurately describe non-linear constitutive equations without the need for iterations.

Problem formulation

We consider here a Lagrangian description of a solid, whose points are described at the initial configuration by the vector X. The solid is denoted by Ω 0 at its initial, undeformed configuration, while its boundary is denoted, as usual, by Γ. The prescribed displacements and tractions are applied to Γ u and Γ t . We assume that Γ=Γ u ∪ Γ t and Γ u ∩ Γ t = ∅. After the body is loaded each material point is described by its position vector, x, in the current deformed configuration.

Following classical approaches to the problem, the deformed configuration of the solid will be given by x = X + u [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF] while the deformation gradient tensor can then be obtained as F = ∂x ∂X . By recalling that the Green-Lagrange strain tensor, E, is given by

E = 1 2 (F T F -1) = γ l (u) + γ nl (u, u) (15)
where we have used the same notation as in [START_REF] Cochelin | The asymptotic-numerical method: an efficient perturbation technique for non-linear structural mechanics[END_REF]:

γ nl (u, v) = 1 2 ∇(u T )∇(v), γ l (u) = 1 2 (∇(u T ) + ∇(u)) (16)
Hyperelastic materials are based upon the definition of a strain-energy function, Ψ. The second Piola-Kirchhoff stress tensor S can be obtained after this function by

S = ∂Ψ ∂E (17)
that is a symmetric tensor and is related to the first Piola-Kirchhoff stress tensor, P , by P = F S.

Using quantities related to the reference configuration the equilibrium equation is as following

∇ • P + B = 0 in Ω 0 ( 18
)
in which B is the body force. Note that, for simplicity, quasi-static processes are assumed. Application of reduced-order modeling to dynamical problems has been successfully treated in [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF].

The boundary conditions of the body are defined by

u(X) = ū on Γ u , P N = λ t on Γ t ( 19
)
where N is the unit vector normal to Γ, t is an applied traction and λ is a loading parameter. The work associated to the equilibrium equation is given by

Ω 0 S : δE dΩ = λ Γ t t • δudΓ ∀δu ∈ H 1 (Ω) (20)
where in the above equation δE is expressed by

δE = 1 2 [F T ∇(δu) + ∇(δu) T F ] = γ l (δu) + γ nl S (u, δu) (21)
where γ nl S (u, δu) is defined by

γ nl S (u, δu) = γ nl (u, δu) + γ nl (δu, u) (22)
The Saint-Venant Kirchhoff model (linear elastic behaviour under large strains) is usually considered as a suitable approximation to most living tissues for real-time purposes, although it seems to be quite limited if a more detailed descriptin is needed. Nevertheless, we can say that the state-of-theart in the field assumes this model as suitable, see for instance, [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF]. It is characterized by the energy function

Ψ = λ 2 (tr(E)) 2 + μE : E (23)
where λ and μ are Lame's constants. The second Piola-Kirchhoff stress tensor can be obtained by

S = ∂Ψ(E) ∂E = C : E (24)
in which C is the classical fourth-order elastic tensor.

An asymptotic numerical method for geometrically nonlinear problems

Asymptotic numerical methos, as created by Potier-Ferry and coworkers [START_REF] Cochelin | The asymptotic-numerical method: an efficient perturbation technique for non-linear structural mechanics[END_REF][24] [START_REF] Yvonnet | A model reduction method for the post-buckling analysis of cellular microstructures[END_REF] (whose notation is also being used here) seek to obtain an analytic expression of the load-displacement curve of a solid or structure in terms of an expansion parameter "a". This expansion of the terms of interest is developed in the neighborhood of a known solution (u n ; S n ; λ n ) at step n. The series is truncated at order N, being λ the loading parameter introduced in the weak form. For Kirchhoff-Saint Venant models it is sufficient to take

⎧ ⎨ ⎩ u n+1 (a) S n+1 (a) λ n+1 (a) ⎫ ⎬ ⎭ = ⎧ ⎨ ⎩ u n (a) S n (a) λ n (a) ⎫ ⎬ ⎭ + N p=1 a p ⎧ ⎨ ⎩ u p S p λ p ⎫ ⎬ ⎭ , (25)
where (u p , S p , λ p ) are unknowns. Above, (u n+1 (a), S n+1 (a), λ n+1 (a)) represents the solution along a portion of the loading curve described continuously with respect to "a". The introduction of Eq. ( 25) into Eq. ( 20) and Eq. ( 24) leads to a series of linear problems with the same tangent operator, a circumstance that seems to be optimal for the purposes of this work.

The series expansion of δE(u) gives

δE n+1 (a) = γ l (δu) + γ nl S (δu, u n ) + n p=1 a p γ nl S (δu, u p ) (26)
The series expansions of S gives in turn

S n+1 (a) = C : E n+1 (a) = C : [γ nl (u n , u n ) + γ l (u n )+ ( 27
) n p=1 a p γ l (u p ) + γ nl S (u n , u p ) + p-1 i=1 γ nl S (u i , u p-i ) ] (28)
at order p we obtain

S p = C : γ l (u p ) + γ nl S (u n , u p ) + p-1 i=1 γ nl S (u i , u p-i ) (29)
Introducing this asymptotic expansion into Eq. [START_REF] Damil | A new method to compute perturbed bifurcation: application to the buckling of imperfect elastic structures[END_REF] gives

Ω 0 S n + N p=1 a p S p : γ l (δu) + γ nl S (u n , δu) + N p=1 a p γ nl S (u p , δu dΩ = (30)
λ n + N p=1 a p λ p Ψ ext (δu) (31)
with Ψ ext (δu) = Γ t t • δudΓ. Introducing Eq. ( 29) into Eq. ( 30) gives a series of problems, one for each level of p (p = 1, . . . , N), that take the form

L(δu, u n ) = λ p Ψ ext (δu) + F nl p (δu) (32)
with

L(δu, u n ) = Ω {S n : γ nl S (u n , u p )+[γ l (δu)+γ nl S (u p , δu)] : C : [γ l (u p )+γ nl S (u n , u p )]}dΩ (33)
and where F nl p (δu) is equal to zero at order one and at order p it reads

F nl p (δu) = - Ω { p-1 i=1 S i : γ nl S (u p-i , δu)+ p-1 i=1 [γ nl S (u i , u p-i )] : C : [γ l (δu)+γ nl S (u n , δu)]}dΩ (34)
Replacing Eq. ( 29) into Eq. ( 30), and discretization of Eq. ( 32) by standard finite element procedure, it leads to a sequence of linear problems in the form order 1

K t u 1 = λ 1 f u T 1 u 1 + λ 2 1 = 1 (35) order p K t u p = λ p f + f nl p (u i ) i < p u T p u 1 + λ p λ 1 = 0 (36)
where K t denotes the tangent stiffness matrix associated with Eq. (33), common to the problems at different orders p. It is the same as the one applied in a classical iterative algorithm like Newton-Raphson, but without updating. Above, u p is the discretized form of the displacement field at order p, f is the loading vector and f nl p is the discretized form associated with F nl p (δu) in Eq. ( 34), which at order p only depends on the values of u i , i < p. The solution of the problem can be obtained as a sequence of p linear problems:

order 1 ⎧ ⎪ ⎨ ⎪ ⎩ û = {K t } -1 f λ 1 = 1 √ ûT û+1 u 1 = λ 1 û (37) order p ⎧ ⎪ ⎨ ⎪ ⎩ u nl p = {K t } -1 f nl p λ p = -λ 1 {u nl p } T u 1 u p = λ p λ 1 u 1 + u nl p (38)

A geometrically nonlinear reduced-order model

As explained before the solution of reduced-order models leads to a very fast calculation of the system of equations but the price to pay is the impossibility of non-linear simulations. On the other hand, by using ANM we can obtain the solution of highly nonlinear problems in the neighbourhood of an equilibrium point very accurately. In this section we combine these two methods to solve geometrically nonlinear problems. This approach was previously established in a very different context by Yvonnet and coworkers [START_REF] Yvonnet | A model reduction method for the post-buckling analysis of cellular microstructures[END_REF] and later generalized in [START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF].

Here we assume that the POD basis has been calculated as explained in section 2. The terms of the asymptotic expansion associated with the displacements are expressed as functions of POD basis as

u p = M m=1 φ m ζ m p = Aζ p , (39)
where ζ p are the degrees of freedom. So the new asymptotic expansion of u is expressed by

u n+1 (a) = A ζ n + N p=1 a p ζ p (40)
where (ζ n , λ n ) represents the previous known solution. Introducing Eq. (40) into Eq. (32) results in

L(Aζ p , Aδζ) = λ p δΨ ext (Aδζ) + F nl p (Aδζ) (41)
after finite element discretization and taking into account that the POD basis is orthonormal, we arrive at the following sequence of linear systems of equations:

at order 1 A T K t A T ζ 1 = λ 1 A T f ζ T 1 ζ 1 + λ 2 1 = 1 (42) at order p A T K t Aζ p = λ p A T f + f nl p ζ T p ζ 1 + λ p λ 1 = 0 (43)
The size of the above equations depends on the number of the POD basis and the chosen level of expansion. In general, very low orders are necessary, as will be demonstrated below. It is noteworthy that these problems have a common tangent stiffness matrix that can be computed off-line with important time savings.

The solution of these equations can be obtained as order 1

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ζ = {A T K t A} -1 A T f λ 1 = 1 √ ζT ζ+1 ζ 1 = λ 1 ζ û = A ζ, u 1 = λ 1 ζ 1 (44) order p ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ζ nl p = {A T K t A} -1 A T f nl p λ p = -λ 1 {ζ nl p } T ζ 1 ζ p = λ p λ 1 ζ 1 + ζ nl p u p = Aζ p (45)

Interpolation of reduced models

From Eq. ( 13) we can notice that the reduced model employs the tangent stiffness matrix, K, linearized from the non-linear problem formulation at a given time instant. Instead of inverting the full stiffness matrix of size N ×N, we employ model reduction techniques to invert the matrix A T K A, of size r × r, that is much lower than the original size, as mentioned before.

However, this tangent stiffness matrix K corresponds to a given state of the model (i.e., a given load position and load value, for instance, in mechanical problems). Different load values would lead to different matrices K along the loading path (this is due to the non-linear character of the problem). This problem is now translated into the sequence of problems (44)-(45). This means that the non-linearity of the problem has somehow been translated to the f nl p term of Eqs. ( 44)-(45), whose evaluation constitutes now the most time-consuming part of the algorithm.

The philosophy of the method proposed in [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF] is to perform some direct simulations for a set of different positions of the surgical tool. For each tool position a reduced-order basis set can thus be computed. But the position of the tool at the on-line simulation does not need to be coincident to the precomputed ones, and therefore some interpolation scheme is needed if the load is placed at positions far from the originally computed ones. This complex interpolation procedure has now two different parts: on one side we need to interpolate the reduced basis sets, from pre-computed complete models, to an arbitrary position of the load, for instance. On the other side, we need to also interpolate the terms f nl p , again from the pre-computed complete models to the arbitrary, reduced ones.

S 0 S 2 S 1 G(n,N) T S0 G(n,N) X = S 1 0 log S ( 1) X S X =Exp ( ) S0

Interpolation of the reduced basis: a geometrical approach

Farhat and coworkers have pointed out that the set of empirical eigenfunctions given by Eq. ( 4) for a given model forms the so-called Grassman manifold G(n, N) [START_REF] Amsallem | An Interpolation Method for Adapting Reduced-Order Models and Application to Aeroelasticity[END_REF]. Therefore, in order to interpolate the set of basis functions, that consequently do not form a vectorial space, we must move to the tangent plane at a point of the manifold, which is a "flat" space, interpolate there, and project back to the manifold, as schematically explained in Fig. 1.

In this way, the columns of A constitute a basis of the subspace S 0 of dimension n of the space R N . At each point S of the manifold G(n, N) one can define a tangent plane of the same dimension, T S , with its points defined by a matrix Γ ∈ R N ×n . The exponential mapping Exp S transforms χ in an n-dimensional subspace S given by a matrix A ∈ R N ×n , such that Γ = U ΣV T (Singular value decomposition)

A = AV cos Σ + U sin Σ
Conversely, the logarithmic mapping Log S , defines a map between a point in the neighborhood of S ∈ G(n, N) and the tangent plane at the origin. Thus, the image of S , in a neighborhood of S, given by the logarithmic mapping, χ = Log S S ∈ T S will be

(I -AA T )A (A T A ) -1 = U ΣV T (Singular value decomposition) Γ = U tan -1 (Σ)V T
So, consider, for instance, A 0 ∈ R N ×n and A 1 ∈ R N ×n , two matrices representing two subspaces, obtained for different parameters of the model (for instance, load positions, but the theory is completely general for other parameters of the model), s 0 and s 1 . Let S 0 , S 1 be the two subspaces originated by considering parameters s 0 and s 1 . Let, in turn, Y(t) be the geodesic line that joins both subspaces (points in the Grassmann manifold), having S 0 as origin. In that case, the initial derivative of the geodesic line, that belongs to the tangent plane at S 0 , will be Ẏ0 = Log S 0 S 1 such that the matrix representing this initial derivative of the geodesic will be (see [START_REF] Amsallem | An Interpolation Method for Adapting Reduced-Order Models and Application to Aeroelasticity[END_REF] for a complete proof of this)

(I -A 0 A T 0 )A 1 (A T 0 A 1 ) -1 = U ΣV T Γ = U tan -1 (Σ)V T
Let S denote the point of the Grassmann manifold representing the reduced-order basis for the new value of the parameter. s. The initial derivative of the new geodesic line, joining S 0 and the sought interpolated subspace S, will be Ẏ0 = r Ẏ0 with r = ss 0 s 1s 0 The computation of the singular value decomposition is not a very timeconsuming task. For instance, on a PC at 2GHz, the computation of the Matlab instruction svd(A), with A a random matrix of 5 × 5000 elements takes on average much less than 1ms, still compatible with real-time constraints.

POD with interpolation

A less rigorous, but much simpler, method to interpolate among previously computed reduced models was established in [START_REF] Ly | Modeling and control of physical processes using proper orthogonal decomposition[END_REF]. This procedure is known in the literature as Proper Orthogonal Decomposition with Interpolation (PODI) [START_REF] Ly | Modeling and control of physical processes using proper orthogonal decomposition[END_REF]. Let [U m (X i )] k i=1 be the response of the system for k different parameter values (in this case, loads at positions X i ), at time step m. The basic algorithm is described as 1. Perform the complete model simulation for each parameter value. 2. Although in the standard PODI technique the POD procedure is applied to the complete set of snapshots of the system to obtain an orthonor-

mal basis A = [φ 1 • • • φ n ],
we proceed by just applying the POD to each complete model (i.e., to each load position). Thus, we obtain an orthonormal basis for each system's parameter value. 3. Interpolate from neighboring basis sets to obtain a reduced basis for the new state of the system. To this end we employ finite element piece-wise linear interpolation.

Results obtained with this technique showed to be much more efficient and accurate than those obtained with the more rigorous technique base upon interpolation on the Grassman manifold. This surprising result is still not well understood by the authors and is currently one of our research topics. In any case, results presented in section 6 were obtained with the PODI technique.

Finally, in order to establish a good set of basis A = [φ 1 . . . φ n ], whatever the chosen interpolation technique, a good sampling strategy should be chosen for the position of the loads in the complete models. In this case, a lattice was established over the domain. Loads were applied at every lattice position. Obviously, other strategies can also be implemented. If the lattice is not dense enough, it is clear that the mentioned interpolation described in this section would lead to less accurate results.

Interpolation of the non-linear force terms

The other big ingredient of the resulting method, namely the non-linear force term f nl p of Eqs. ( 44)-( 45) needs also to be interpolated between precomputed models, since it is too computationally-demanding to be evaluated in real-time. To this end, we establish a piece-wise polynomial (FE-like) interpolation among complete models. This simple procedure has rendered excellent results, as will be shown in the subsequent section.

Numerical examples

In order to show the performance of this method we have applied it to some academic examples devoted to the simulation of soft tissues. The first one is the pinching of a model of the human cornea, while the second is devoted to liver palpation.

Pinching the human cornea

In this example we applied 0.014N to nine neighboring nodes located at the central region of the cornea, that was meshed using trilinear hexahedral elements. It consisted of 8514 nodes and 7182 elements. The mesh is shown in Figure 2 in two views. The material properties of the cornea were assumed to be E = 2MP a and ν = 0.48. A more detailed constitutive model can be found at [START_REF] Alastrue | Biomechanical modeling of refractive corneal surgery[END_REF].

In this example nine modes were applied that provide decent approximation. The solution has been obtained using ANM-POD for p = 1, . . . , 6. The tangent stiffness matrix at the origin, that is independent of the load state of the organ,is employed for all the simulations. In order to verify the results we have computed the solution by full FE Newton-Raphson method. The loading factor (λ) has been plotted versus the maximum displacement (in modulus) in Figure 3.

As it can be seen the results have good accuracy with Newton-Raphson solution. The deformed cornea obtained using ANM-POD for λ = 1 is shown in Figure 4 and the one obtained using full FE Newton-Raphson is depicted in Figure 5.

If we consider a case in which the load is placed at a position for which a reduced model has not been stored, the interpolation procedure, as described in section 5, is applied. As mentioned before, the tangent stiffness matrix at origin is identical for all load states. The procedure follows by (linearly) interpolating the reduced basis set from the four nearest neighbor pre-computed states and also by interpolating the non-linear force term f nl p for each order of the expansion.

The resulting load-displacement curve is depicted in Fig. 6. The excellent agreement with the complete model results (denoted by p FEM in the legend) is noticeable.

Palpation of the liver

In this example we consider the palpation of a liver with a surgical tool, assumed perfectly rigid. The liver is the biggest gland in the human body, after the skin. Liver geometry has been obtained from the SOFA project The resulting behavior is interpolated among four previously stored reduced models.

[30] and post-processed in order to obtain a mesh composed by 2853 nodes and 10519 tetrahedra, whose geometry is shown in Fig. 7. The anterior surface of the liver is considered free, while the posterior one was assumed to be supported over different organs (it is connected to the diaphragm by the coronary ligament, for instance). The inferior vena cava travels along the posterior surface, and the liver is frequently assumed clamped a that location. Although the assumed boundary conditions are not strictly correct from a physiological point of view, our main interest is to show that the model can be solved under real-time constraints with reasonable accuracy.

Although the literature on the mechanical properties of the liver is not very detailed, we have assumed a Young's modulus of 160 kPa, and a Poisson coefficient of 0.48, thus nearly incompressible [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF]. As in the previous example, a load of 1.2N has been applied at an arbitrary point and the reduced model has been constructed, composed in this case by the 9 modes depicted in Fig 8, that capture the 99.9% of the energy of the system.

The FE complete model was left to run until lack of convergence of the Newton-Raphson algorithm. A standard Newton-Raphson algorithm was employed, without any modification such as arc-length methods or anything similar. The load (in terms of load factor λ)-displacement curve obtained by POD-ANM techniques is depicted in Fig. 9. It is worth noting the high accuracy of the results for a wide range of displacements (up to ≈ 5mm under the tool tip) for expansions of order 5 and 6, less than 5%. This accuracy can also be noticed from Figs. 10 and 11, representing the displacements for FE complete mode and POD-ANM (sixth-order expansion) techniques, respectively.

Timings

The examples below have been performed on a laptop running Matlab on a Windows-operated computer at 2 GHz., with 2 Gb RAM. The code prototype is obviously thought for the design and testing of new algorithms. Under these circumstances, the examples ran at around 20 Hz, which is enough for visual perception, but not for haptic environments. This can be improved by considering more sophisticated programming and not a generalpurpose language such as Matlab. The summary of the computational cost for a typical simulation is as follows:

• Interpolation of the basis and the non-linear force term: 3.9 • 10 -3 s.

• Construction of the reduced stiffness matrix A T K t A (note that the complete tangent stiffness matrix at the origin, K t , is identical for all simulations and therefore could be previously stored in memory): 4.7 • 10 -2 s. • Solving the resulting system of equations (6 × 6 in this particular example): 3.1 • 10 -4 s.

The use of Graphics Processing Units [START_REF] Taylor | High-speed nonlinear finite element analysis for surgical simulation using graphics processing units[END_REF] is obviously another possibility to achieve true real-time performance, that seems to be at hand.

Conclusions

In a previous work by the authors [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF], model reduction techniques were applied for real-time simulation of the palpation of the human cornea. It was demonstrated that even for very sophisticated constitutive laws, model reduction techniques provided reasonable results at real-time frequencies for haptic environments. However, it has been noticed that, although reduced models obtained by POD techniques are optimal in some statistical sense, no tangent stiffness matrix updating was possible under such severe constraints imposed by real-time. Thus, the models obtained by employing this technique are still linear.

In this paper we have pursued a different strategy. Through the combination of model reduction techniques and an Asymptotic Numerical Method (ANM) we have shown that it is possible to obtain an expansion of the solution in a neighborhood of a previous equilibrium state that provides a very simple algorithm that fits reasonably well to the true solution within some radius of convergence. The main feature of the method is that it possesses only one tangent stiffness matrix that does not need to be updated within a large interval of strains (much larger that traditional Newton-Raphson schemes). It has been show how this radius of convergence is enough for many applications, involving strains on the order of 10%.

Other approaches based upon reduced models exist, see [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF], for instance. The main difference between the approach followed in this work and that of Barbic and James is that in the latter case a general-purpose reduced basis is employed. These basis functions are obtained from mass-scaled principal component analysis. These basis functions are thought to be optimal for any load state of the solid. In the present approach different sets of basis functions are obtained for different load (and possibly boundary conditions) states that can be changed according with the user experience. Basis sets for any state different from the precomputed ones are obtained after interpolation, as explained in the previous sections. In section 5.1 we have explained how the set of reduced models for different loading states of the organ forms a manifold. Therefore, it is crucial for the method to work well that the set of complete models to be solved to feed the algorithm is chosen adequately. For the moment, the brute-force approach seems to be out of reach. In our opinion, this process should be guided by surgeons, by indicating the most probable loading states for each surgical procedure, distinguishing, for instance, between experienced surgeons and medicine students.

The method has been developed for Saint Venant-Kirchhoff constitutive laws. This kind of laws is among the best state-of-the-art existing models [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF], but is judged to be poor for some applications. It is well-known that it suffers from instabilities when subjected to compression (although no instabilities have been found in the examples tested by the authors), so it would be interesting to extend the proposed technique to some more sophisticated constitutive laws, such as neo-hookean laws, for instance, or other much better suited for specific applications [START_REF] Alastrue | Biomechanical modeling of refractive corneal surgery[END_REF].

Despite the facts commented above, the proposed technique seems to be an appealing method for the simulation of linear elastic materials undergoing large strains at real time. The authors continue to work in the development of the technique. For instance, the issue of the simulation of cutting in the context of reduced models is also a bottleneck of the method. This constitutes

Figure 1 :

 1 Figure 1: Schematic description of the Grassman manifold formed by the set of basis functions of all the reduced modes of a given model.

Figure 2 :

 2 Figure 2: Geometry of the finite element model for the human cornea.

Figure 3 :

 3 Figure 3: The loading factor vs maximum displacement [mm] for the Pinched cornea.

Figure 4 :

 4 Figure 4: u y , [mm], contour of the pinched cornea obtained by ANM-POD.

Figure 5 :

 5 Figure 5: u y , [mm], contour of the pinched cornea obtained by Newton-Raphson.

Figure 6 :

 6 Figure 6: The loading factor vs minimum displacement, [mm], for the Pinched cornea. The resulting behavior is interpolated among four previously stored reduced models.

Figure 7 :

 7 Figure 7: Geometry of the finite element model for the liver.

Figure 8 :

 8 Figure 8: Nine most important eigenmodes for the simulation of the liver palpation. The corresponding eigenvalues are: 2.69E -03 (a), 8.28E -05 (b), 8.82E -07 (c), 1.07E -08 (d), 6.73E -11 (e), 5.21E -13 (f), 2.52E -15 (g).

Figure 9 :

 9 Figure 9: Load-displacement, [mm], curve for the liver palpation for different orders of approximation. The continuous blue line represents the solution for the complete model employing Newton-Raphson methods.

Figure 10 :

 10 Figure 10: Result for the FEM model, displacement u y , [mm].

Figure 11 :

 11 Figure 11: Result for the reduced-ANM model, displacement u y , [mm].

our current effort of research and will be published elsewhere.