Huu Phuoc Bui

Satyendra Tomar

Stéphane Pierre

Alain Bordas

P A Bordas

Stéphane P A Bordas
email: stephane.bordas@alum.northwestern.edu

Corotational cut finite element method for real-time surgical simulation: Application to needle insertion simulation

Keywords: Cut finite element method, Unfitted FEM, Corotational CutFEM, Needle insertion, Real-time simulation

come   L'archive ouverte pluridisciplinaire

Introduction

Nowadays real-time simulation plays an important role in different fields: graphic animation [START_REF] Musse | Hierarchical model for real time simulation of virtual human crowds[END_REF], fracture of stiff materials [START_REF] Müller | Real-Time Simulation of Deformation and Fracture of Stiff Materials[END_REF], and surgical training and simulation [START_REF] Cotin | Real-time elastic deformations of soft tissues for surgery simulation[END_REF][START_REF] Cotin | A hybrid elastic model for realtime cutting, deformations, and force feedback for surgery training and simulation[END_REF][START_REF] Monserrat | A new approach for the real-time simulation of tissue deformations in surgery simulation[END_REF][START_REF] Courtecuisse | Real-time simulation of contact and cutting of heterogeneous soft-tissues[END_REF], to name a few. In the medical context, surgical simulations are not only useful for training, but also helpful for pre-operative planning, and intra-operative guidance. Surgical simulations have to take into account interactions between a surgeon or an interventional radiologist with a deformable organ via surgical instruments (e.g. a needle), and also interactions between the organ with its neighbouring structures. To be useful in real-life situations, it is required that the computations are performed in real-time. To achieve real-time performance, some advanced solvers, e.g. GPU-based computation [START_REF] Courtecuisse | GPU-based real-time soft tissue deformation with cutting and haptic feedback[END_REF], or asynchronous solver [START_REF] Courtecuisse | Real-time simulation of contact and cutting of heterogeneous soft-tissues[END_REF], can be used. The computational time can also be reduced by using coarse meshes. However, using coarse meshes, some geometric details are often lost. Simulations using coarse meshes are thus only suitable for targeted surgical training which relies more on visual reality than on actual one. Such simulations fall short for surgical planning or guidance where computations must provide accurate results. A model order reduction technique can also be used to solve system equations with significantly less time but it results in lower accuracy [START_REF] Niroomandi | Real-time simulation of surgery by reduced-order modeling and x-fem techniques[END_REF]. To reduce the computational effort, Quesada et al [START_REF] Quesada | Computational vademecums for real-time simulation of surgical cutting in haptic environments[END_REF] propose a computational parametric meta-model which is computed offline, and is only evaluated online.

Medical simulations have to deal with complex anatomical structures, e.g. prostate, blood vessel, liver, brain, brain ventricle, etc. When the patient-specific geometry is considered, the mesh of the organ needs to be reconstructed since the organ geometry is different from one patient to another. Misra and coworkers [START_REF] Misra | The importance of organ geometry and boundary constraints for planning of medical interventions[END_REF] have shown that the geometry of the organ and boundary conditions surrounding the organ are the most important factors influencing the organ deformation, and thus have a direct impact on the accuracy of simulation and planning. And when supercomputers or parallel computation is not considered, in order to response to real-time simulation requirements, computations involving interactions with surgical tools and/or cutting operations are performed on coarse meshes while applying pre-computed deformation from fine meshes, see e.g., [START_REF] Cotin | A hybrid elastic model for realtime cutting, deformations, and force feedback for surgery training and simulation[END_REF]. However, these coarse and fine meshes may not conform to each other, resulting in different geometric descriptions and different boundary conditions, and thus the question of simulation accuracy demands further investigation.

The use of geometric description in a computational method, so that the users only need to provide a background mesh which is not necessarily conforming to the boundary geometry of the simulated object, can dramatically reduce computational cost of preprocessing. This is the idea behind the fictitious domain method [START_REF] Burman | Fictitious domain finite element methods using cut elements: I. A stabilized Lagrange multiplier method[END_REF][START_REF] Burman | Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method[END_REF], or the cut finite element method (Cut-FEM) [START_REF] Burman | CutFEM: Discretizing geometry and partial differential equations[END_REF]. The extended finite element method (XFEM) [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], and the generalised finite element method (GFEM) [START_REF] Strouboulis | The design and analysis of the generalized finite element method[END_REF] have been developed to deal with crack surfaces, or material interfaces evolving during simulation. These approaches allow for computation of coupled physical processes on distinct subregions of the total volume, and do not require an absolute conformity between the meshes from mesh construction as in [START_REF] Tabor | Mesh construction from medical imaging for multiphysics simulation: Heat transfer and fluid flow in complex geometries[END_REF][START_REF] Moumnassi | Finite element analysis on implicitly defined domains: An accurate representation based on arbitrary parametric surfaces[END_REF]. In the context of real-time patient-specific surgical simulations, using a unified geometric description in a finite element code can help in [START_REF] Musse | Hierarchical model for real time simulation of virtual human crowds[END_REF] reducing preprocessing cost of mesh generation, and (2) implicitly improving accuracy of simulation and planning (compared to the situation in which coarse meshes are employed in an inappropriate manner to achieve real-time requirements).

The main contribution of our paper is to propose the corotational CutFEM approach which is suitable for real-time patient specific simulations. The corotational model [START_REF] Moita | A finite element formulation for 3-D continua using the co-rotational technique[END_REF] is widely used for the treatment of large rotations of soft tissues, see e.g., [START_REF] Courtecuisse | GPU-based real-time soft tissue deformation with cutting and haptic feedback[END_REF][START_REF] Suwelack | Quadratic Corotated Finite Elements for Real-Time Soft Tissue Registration[END_REF][START_REF] Haouchine | Image-guided simulation of heterogeneous tissue deformation for augmented reality during hepatic surgery[END_REF][START_REF] Dequidt | Vascular Neurosurgery Simulation with Bimanual Haptic Feedback[END_REF][START_REF] Sase | GPU-accelerated surgery simulation for opening a brain fissure[END_REF][START_REF] Morin | Brain-shift compensation using intraoperative ultrasound and constraint-based biomechanical simulation[END_REF]. By using the proposed methodology, geometries of simulated objects can be captured automatically. We also show that, for a given accuracy, using the proposed method is computationally more advantageous than employing the standard FEM. On the implementation aspect of the method, we propose the so-called multilevel embedding approach to correctly integrate implicit boundaries of the simulated organ, and to accurately capture implicit interfaces of e.g., a tumour. The implementation is verified by convergence studies, which exhibit optimal rates. We use Lagrange multipliers to impose Dirichlet boundary conditions on implicit boundaries, whereas traction or Neumann boundary conditions, which is/are applied on parts of the surface, can be distributed to the background nodes using shape functions. We demonstrate the performance of the corotational CutFEM through various applications: from needle insertion simulations (e.g. for biopsy or brachytherapy) to simulation of electrode lead implantation in Deep Brain Stimulation (DBS) procedure. The algorithm is implemented in the open-source SOFA framework [START_REF] Faure | SOFA: A Multi-Model Framework for Interactive Physical Simulation[END_REF] 1 .

The remaining of the paper is organised as follows. In Section 2, we describe the formulation of the needle insertion problem into soft tissues, together with its discrete form. We then present the geometry discretisation of cut elements, and the algorithm for multilevel embedding approach to correctly integrate implicit boundaries/interfaces. It is then followed by the discussion on implementation aspects. Next, a corotational formulation for CutFEM is briefly introduced. It is followed by the description of how boundary conditions are applied on implicit surfaces. This section is concluded with the discussion on solving the system equations with constraints. Numerical results are presented in Section 3, which demonstrate the capabilities of the approach through various needle insertion problems into liver and brain models. Finally, conclusions are drawn in Section 4.

Methods

Problem setting

In the context of needle insertion into soft tissue, we model both, the needle and the tissue, as dynamic deformable objects. Creating a good quality conforming mesh, when the tissue is modelled as heterogeneous material with complex internal structures, or when the tissue has a complex geometry, is a very complicated and challenging task. To overcome this burden, CutFEM has been proposed as it does not require conforming meshes. Figure 1a schematically shows a problem in which an interface is immersed into a tissue geometry for simulating e.g. a tumour geometry. Figure 1b shows a problem in which the tissue is simulated with an implicit boundary (i.e. the computational mesh is not fitted to the tissue geometry). Let Ω represent the tissue domain, and ∂Ω denote its boundary. The tissue undergoes an imposed displacement ū on the boundary part Γ u , and a traction force t on the boundary part Γ t . The governing equation of the problem reads

ū t λ -λ Γ u Γ t Γ ∂Ω (a) Immersed interface problem. ∂Ω ū Γ u t Γ t λ -λ Ω (b) Implicit boundary problem.
div σ + b + λ = ρ ü in Ω, ( 1a
) σ • n = t on Γ t , (1b)
u = ū on Γ u , ( 1c
)
where u is the displacement field of the object, σ is the Cauchy stress tensor, b is the generalised body force vector, ρ is the mass density, ü is the second partial derivative of u with respect to time, n denotes the outward unit normal vector on Γ t , and λ denotes the interaction force between the needle and the tissue. Within the dynamic behaviour consideration, the generalised body force b is expressed as

b = b -c u, ( 2
)
where b is the body force per unit volume, -c u expresses the resistance opposite to the motion velocity u, and c is the viscosity parameter which can be considered to be numerically given value [START_REF] Chopra | Dynamics of Structures, Prentice-Hall international series in civil engineering and engineering mechanics[END_REF], see Section 2.2. The strain, which is a deformation measure, is expressed through the displacement as

= 1 2 ∇ u + (∇ u) T . (3)
In this paper, we employ a linear elastic material based on corotational formulation, see also Section 2.4, so that the stress-strain relationship is expressed as

σ = λdiv(u)I + 2µ , (4)
where λ and µ are the Lamé's coefficients, and I denotes the identity matrix. The Lamé's coefficients are computed from the Young's modulus E and Poisson's ratio ν as

λ = Eν (1 + ν)(1 -2ν) , µ = E 2(1 + ν) . (5)
The interaction force λ is defined from the interaction law between the needle and the tissue. Three type of constraints between the needle and the tissue are defined during needle insertion simulations: constraint for puncturing at tissue surface, constraint at needle tip, and constraints along the needle shaft, see [START_REF] Bui | Real-time error control for surgical simulation[END_REF][START_REF] Bui | Controlling the error on target motion through real-time mesh adaptation: Applications to deep brain stimulation[END_REF] for details. For the sake of completeness, these constraints are recalled here. When the needle has not been penetrated into the tissue yet, the contact between them verifies the Kuhn-Tucker condition in the direction n (normal to tissue surface)

δ n ≥ 0; λ ts n ≥ 0; δ n • λ ts n = 0, (6)
where λ ts n stands for the normal component of the constraint on the tissue surface, δ n denote the distance between the needle tip and the tissue surface in the direction n. This condition expresses the fact that the contact force only exists when the needle tip is in contact with the tissue surface (δ n = 0). In the tangential direction, a frictional contact is expressed through the Coulomb's friction law

λ ts t < µ f λ ts n (stick); λ ts t = µ f λ ts n (slip), (7)
with µ f denoting the coefficient of friction. When the normal component of the constraint is greater than the tissue puncture strength λ p 0 , which is read

λ ts n > λ p 0 , (8)
the needle tip can penetrate into the tissue. Once inside the tissue, at the needle tip, a constraint is defined between the latter and the tissue. This constraint expresses the phenomenon that the tip is stuck, or can cut and advance inside the tissue. The condition for this constraint reads

λ nt n < µλ nt t + λ c 0 (stick); λ nt n ≥ µλ nt t + λ c 0 (cut and slip), (9)
where n now denotes the normal direction which is the direction of the needle shaft at the tip, the superscript nt indicates the constraint at the needle tip, λ c 0 is the cutting strength of the tissue. Since the needle shaft should follow the trajectory created inside the tissue by the needle tip when the latter advances, constraints along the trajectory are defined between the tissue and the needle shaft which fulfil the Coulomb's friction law λ ns n < µλ ns t (stick); λ ns n = µλ ns t (slide), [START_REF] Misra | The importance of organ geometry and boundary constraints for planning of medical interventions[END_REF] where the superscript ns indicates the constraints along the needle shaft. These three types of constraint are schematically illustrated in Figure 2. It is noted that the total interaction force between the needle and the tissue, at each time step, is nothing different from the normal force f n computed at a cross-section of the needle model which has not been penetrated into the tissue yet, as illustrated in Figure 2.

n t n n n t t t t d d d f n f n Figure 2:
Schematic representation of constraint at the tissue surface λ ts (shown by red colour), constraints defined along needle shaft λ ns (illustrated by green colour), and constraint at the needle tip λ nt (depicted by blue colour). These constraints are defined on a corresponding local frame nt, and they are located at the corresponding constraint points separated by a distance d. The interaction force between the needle and the tissue is simply represented by the normal force f n of the needle model, computed at the non-penetrated (according to the tissue) cross-section of the needle shaft.

Weak form and FEM discretisation

For the numerical studies in the following, we employ homogeneous boundary condition on the Dirichlet boundary Γ u i.e. ū = 0. We denote the Sobolev space in three dimensions as (H 1 (Ω)) 3 , and then define the space of admissible displacements V as

V = {v ∈ (H 1 (Ω)) 3 : v = 0 on Γ u }. ( 11
)
To begin with, for the sake of simplicity, we omit the interaction force λ in Equation (1a), and only consider the tissue/needle without any interaction with the other. The constraints for the interaction force will be dealt later. By multiplying Equation (1a) by a test function v ∈ V , and integrating by parts on the domain Ω, we obtain

Ω divσ • v dΩ + Ω b • v dΩ = Ω ρ ü • v dΩ. ( 12
)
After performing the integration by parts to the integrand of the first term of the left-hand side in Equation [START_REF] Burman | Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method[END_REF], and noting that v vanishes on the boundary part Γ D , we end up with this variational formulation

         Find a displacement u ∈ V such that, ∀v ∈ V , Ω ρ ü • v dΩ + Ω σ(u) : ε(v) dΩ + Ω c u • v dΩ q(u,v) = Ω b • v dΩ + Γt t • v dΓ l(v) (13
) To solve the problem described by Equation ( 13) numerically, we need to transform the continuous variational problem [START_REF] Burman | CutFEM: Discretizing geometry and partial differential equations[END_REF] to a discrete variational one. This is performed by introducing a finite-dimensional space, denoted by V h ⊂ V . It is noted that in the CutFEM framework, for tissue deformation simulations, the space V h is defined on the mesh obtained from elements located inside the tissue surface and from elements cut by the tissue surface, see Section 2.3. The problem now reads

Find a displacement u h ∈ V h such that, ∀v h ∈ V h , q(u h , v h ) = l(v h ). ( 14
)
The space V h is constructed from a mesh Ω h obtained by spatially discretising the domain Ω into n e finite elements Ω e , e = 1, 2, . . . , n e , using n n nodes. For each element Ω e , we define an interpolation function (called shape function), denoted by N e = [N e j ], j = 1, 2, . . . , m e with m e the number of nodes of the element e. Mathematically, these are linear (P 1 ) elements in three-dimensions. By integrating every term of Equation ( 14), and then assembling for the whole domain Ω h , we obtain the following discrete problem, in which the subscript h is dropped for the sake of simplicity (see, e.g., [START_REF] Zienkiewicz | The Finite Element Method: Solid mechanics, Referex collection[END_REF][START_REF] Liu | Chapter 3 -Fundamentals for Finite Element Method[END_REF])

M ü + C u + Ku = f ext , ( 15
)
where M is the mass matrix, K is the stiffness matrix, C is the damping matrix, and f ext is the external force vector of the whole system. These quantities are assembled from the corresponding element quantities: M e = Ωe N eT ρN e dΩ, C e = Ωe N eT cN e dΩ, K e = Ωe B eT EB e dΩ, in which B e is the strain-displacement matrix computed from the derivative of the shape function N e , and E is the fourth order tensor describing the stress-strain relationship σ = E : ε. Moreover, f eext is computed from the body force b and the traction t acting on the element e. In this paper, a lumped mass matrix is employed, where a diagonal mass matrix (from the mass density ρ) is integrated over the volume of each element. The stiffness matrix K is computed using the corotational FEM for both, the needle and the tissue, see Section 2.4. In practice, it is difficult to determine the damping matrix C as knowledge of viscous parameter c is lacking. In our simulation, we employ Rayleigh damping which is a linear combination of stiffness and mass matrices, see [START_REF] Zienkiewicz | The finite element method: Its basis and fundamentals[END_REF] C = αM + βK.

We set α = 0.01, and β = 0.1 in our simulations.

The relation [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] can be rewritten as

M a = f (x, v), (17)
where x, v = u, and a = ü, represent the position, velocity, and acceleration, respectively, and f (x, v) = f ext -Ku -Cv represents the net force (the difference of the external and internal forces) applied to the object. It is noted that we will use corotational formulation, see [START_REF] Felippa | A unified formulation of small-strain corotational finite elements: I. theory[END_REF], to better capture the tissue/needle deformation under large rotational deformation, see also [START_REF] Müller | Stable real-time deformations[END_REF], which often occurs in surgical simulations. Therefore, the system term Ku is computed from the contribution of each element as will be described in Section 2.4. For temporal discretisation of Equation [START_REF] Tabor | Mesh construction from medical imaging for multiphysics simulation: Heat transfer and fluid flow in complex geometries[END_REF], i.e. to numerically solve the problem in time, an implicit backward Euler scheme [START_REF] Baraff | Large steps in cloth simulation[END_REF] is used to discretize the velocity and position as follows ut+τ = ut + τ üt+τ ;

u t+τ = u t + τ ut+τ , (18)
where τ denotes the time step. Inserting Equation [START_REF] Moumnassi | Finite element analysis on implicitly defined domains: An accurate representation based on arbitrary parametric surfaces[END_REF] into Equation [START_REF] Tabor | Mesh construction from medical imaging for multiphysics simulation: Heat transfer and fluid flow in complex geometries[END_REF] yields the final discrete system

(M + τ C + τ 2 K) A dv = τ f (x t , v t ) -τ 2 Kv t b ( 19
)
or simply Adv = b, where dv = v t+τ -v t . Note that here we use b to denote the right-hand side of Equation [START_REF] Moita | A finite element formulation for 3-D continua using the co-rotational technique[END_REF], and it does not mean the generalised body force described in Equation (1a).

When taking into account the interaction between the needle and the tissue, Equation ( 19) is rewritten as

Adv = b + H T λ, (20)
in which H T provides the direction of the constraints. The interaction constraint λ, between the needle and the tissue, is computed using Lagrange multipliers, see Section 2.6. To understand the matrix H, we suppose that the needle is going to have contact with the tissue surface at the point P , as illustrated in Figure 3. When they have been in contact, the relative displacement between the needle tip Q and the contact point P in the direction of the contact n vanishes

δ n = (u P -u Q ) • n = 0. (21)
Furthermore, the displacement at the contact point P on the tissue surface can be interpolated from the nodal displacements of the corresponding triangle. Also, the displacement of the needle tip can be interpolated from the corresponding needle element (1D element embedded in 3D space, described later)

u P = i N 1 i (P )u i ; u Q = i N 2 i (Q)u i , (22)
where N 1 i , and N 2 i are the shape functions defined for the elements of the tissue surface and the needle, respectively. By introducing Equation [START_REF] Dequidt | Vascular Neurosurgery Simulation with Bimanual Haptic Feedback[END_REF] into Equation ( 21), we can write

H 1 u 1 -H 2 u 2 = 0, ( 23
)
in which u 1 , and u 2 are the nodes involving due to contact of the tissue surface and of the needle, respectively. By a similar approach, we can write the interaction force due to contact applied to the involving nodes on the tissue and the needle as

H 1 T λ = -H 2 T λ, (24)
with λ the interaction force in the contact space. It is noted that, when using CutFEM for simulation of soft tissue, the tissue surface is not conforming to the computational mesh, the force H 1 T λ is again interpolated into the involving nodes of the computational mesh, see also Section 2.5. After solving [START_REF] Suwelack | Quadratic Corotated Finite Elements for Real-Time Soft Tissue Registration[END_REF] for dv, the position and velocity for needle and tissue are updated as

v t+τ = dv + v t ; x t+τ = x t + τ v t+τ . (25)
The tissue domain is discretised by a tetrahedron mesh. For the tetrahedra that intersect with the immersed interface/implicit boundary, an embedded element set is employed to facilitate the integration, see Section 2.3.

Since the length of the needle is much greater than the dimensions of its cross section, it can be safely assumed that the cross section, which is normal to the mid-line (i.e. the neutral line) of the needle, remain planar and normal to the mid-line during deformation. This makes it possible to use the Euler-Bernoulli beam theory [START_REF] Timoshenko | History of Strength of Materials: With a Brief Account of the History of Theory of Elasticity and Theory of Structures[END_REF] to describe the behaviour of the needle. The needle is then discretised by one dimensional elements, and suitable shape functions (C 1 continuous, e.g. Hermite) are employed for the interpolation of the displacement field so that the bending behaviour of the needle is well taken into account. By this, each node of the elements used for the needle has 6 degrees of freedom (3 translations and 3 rotations).

Geometry discretisation for immersed/implicit boundaries

The discretisation technique presented here is applicable for both, fictitious domain problems [START_REF] Burman | Fictitious domain finite element methods using cut elements: I. A stabilized Lagrange multiplier method[END_REF][START_REF] Sotiropoulos | Immersed boundary methods for simulating fluidstructure interaction[END_REF] and nonconforming interface problems, see, e.g. [START_REF] Qin | A Cartesian grid nonconforming immersed finite element method for planar elasticity interface problems[END_REF].

Geometry discretisation

In fictitious domain or nonconforming interface approaches, the boundary of a given domain is embedded into a computational mesh. The latter is used to approximate the solution of the governing PDEs. In general, the boundary surface of the given domain is represented by a very fine mesh, and the boundary surface mesh and the computational mesh are not assumed to be conforming. In fictitious domain method, the governing equations are integrated only on the inside volume (i.e. Ω 2 , see Figure 4a) bounded by the surface boundary, whereas in interface problems, during integration, different mechanical properties are assigned to the outside volume Ω 1 and inside volume Ω 2 , which are separated by the surface, see Figure 4a. The domain is discretised by a mesh which is nonconforming with respect to the surface Γ, as schematically shown in Figure 4b.

In order to accurately integrate on a given domain, the elements are firstly classified into three categories: cut elements, inside elements (i.e. inside of Γ) and outside elements (i.e. outside of Γ), see Figure 4c. To facilitate integration, each cut element is then embedded with a sub-element set consisting of elements which are conforming with the surface Γ. It is important to note that the sub-element set is only used for integration purpose, and since the degrees of freedom of the cut elements are still defined only on their nodes, the sub-element embedding approach does not affect the approximation properties of the discretisation at all. To identify the cut elements, we use the level-set method [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF]. Figure 5 schematically shows an element which is potentially cut by the interface. To know if the element is cut by the interface, we define the level-set function as a signed distance function from the nodes of the element to the interface, and check the sign of that function. An edge P i P j of the element is cut by the interface if and only if

Ω 1 Ω 2 Γ (a) (b) (c)
δ(P i ) • δ(P j ) < 0, (26)
where P i and P j denote the position of the two ends of the edge, and δ is the signed distance defined from the points to the interface, see Figure 5. The intersection between the interface and an edge is approximated by the zero level-set. When Equation ( 26) is fulfilled, the barycentric coordinate of the intersection reads

ξ = |δ(P i )| |δ(P i )| + |δP j | , (27)
and the position of the intersection is computed as In general applications, the immersed interface has a complex geometry, and it is not necessarily a convex surface, see Figure 9a for a liver surface as an example. Therefore, using the level set function to identify the inside and outside elements may not be efficient, especially when the elements are far from the interface. In Section 2.3.4, an efficient algorithm is presented to identify inside and outside elements when the cut elements are already marked using the level-set method. Moreover, in Section 2.3.4, the details of embedding each cut element by a conforming subtriangulation for integration purpose is presented.

P = (1 -ξ)P i + ξP j . ( 28
)
P 3 P 1 P 2 P 1 P 2 n δ(x) δ(P 1 ) < 0 δ(P 3 ) < 0 δ(P 2 ) > 0

Refinement for invalid cut elements

In Section 2.3.1, we assume that elements of the background mesh are cut by the surface interface with valid cases: i.e. the intersection between a tetrahedron element and the surface is either a triangle or a quadrilateral, see Figure 9b. The number of intersection between the surface and the element edges is either 3 or 4. However, invalid cutting cases arise when, for example, an edge of an element is cut by the surface with more than one intersection, or there are only two, or more than 4 edges of a tetrahedron, which are cut by the surface. These invalid cases can arrive when a coarse background mesh is used together with a curved interface surface. Figures 6a,6b, and 6c schematically show, in two dimensions, some invalid cut cases between a triangle and an interface curve.

To overcome this issue, one solution, as in [START_REF] Fries | Higher-order accurate integration of implicit geometries[END_REF], is to recursively embed the invalid cut element with a set of sub-elements until all (sub) cut elements are valid, see Figures 6d, 6e, and6f. In what follows, this procedure is called refinement. However, it is important to note that the invalid cut element is not actually refined since we do not introduce any new degree of freedom into the background mesh. The number of level of refinement needed to get all valid cut elements depends not only on the coarseness of the background mesh and curvature of the interface surface, but also on the relative location between the cut element and the surface. In Figures 6d,6e, only one refinement level is needed in order to get valid cut elements, whereas in Figure 6f, two refinement levels are necessary. Once we get all valid cut elements, they are embedded by a conforming subtriangulation as usual for integration reason, described above. Figures 6g,6h, and 6i show the conforming embedded elements for valid cut (sub) elements.

As an example to demonstrate the implemented algorithm which works on tetrahedra, Figure 7a shows a spherical surface which is immersed into a background mesh, and valid and invalid cut elements with the spherical surface are embedded with sub-tetrahedra shown in Figure 7c. Figure 7b shows four levels of refinement needed to capture the intersections between an invalid cut tetrahedron and the spherical surface.

Stabilisation by moving background nodes

There may be situations in which an element is cut by the surface with a very small intersection, as schematically shown in Figure 8a. These situations can strongly affect the stability of the system matrix. To remedy this stability issue, we propose to move the concerned node by a random distance, which is proportional to the element size, defined by the vector with 0 < β < 1 a random variable, l e the element size, and n the outward unit normal of the surface. In practice, β is set between 0.1 and 0.9 which can efficiently overcome the issue. This approach is schematically shown in Figure 8b. We define a criterion to detect when the intersection between an element and a surface is small, and thus the moving node approach is applied. This criterion reads

(a) (b) (c) (d) (e) (f) (g) (h) (i)
d = - 1 2 β • l e • n, (29)
V i V e < α, (30)
where V i is the volume of the smaller part of the element cut by the surface, and V e is the total volume of the element. In this study, we set α to 0.05, i.e., when the volume of the smaller part of the element is less than 5% of the total volume of the element, we move the node.

It is important to note that this moving node approach is also applicable when refinement method described in Section 2.3.2 fails to determine valid cut elements. In fact, if after five levels of refinement, invalid cut elements are still present, we simply move, by a distance vector expressed in Equation ( 29), the element node which is closest to the intersection between the element and the surface.

Implementation aspects

To discuss the implementation aspects, as a starting point, we assume that an arbitrary surface is immersed into a computational background mesh. An example of a liver surface being immersed in a computational tetrahedron mesh is shown in Figure 9a.

Since the interface surface is arbitrary, and may be concave, due to the change of orientation of the outward normal of the surface, it may not be efficient to use only level-set function to distinguish between three types of elements described above. Moreover, to mark those elements which are far from the surface, using the distance function as level-set function, in combination with outward normal of the surface, raises the ineffectiveness of the algorithm. To remedy this issue, we mark different types of elements using the following approach, see also Figure 10a.

Step 1 Using level-set function (as described in Section 2.3.1), mark the cut elements by checking the potential intersections between each triangle of the surface mesh with immediate neighboring elements (tetrahedra) of the triangle. Label these elements by 1,

Step 2 Mark the outside elements by propagating the marking procedure, starting from the elements located at the boundary ∂Ω of the domain until a cut element is reached. Label these elements by 0,

Step 3 Mark the remaining elements as inside. Label these elements by 2.

(a)

Triangle Quadrilateral

(b)

Figure 9: A liver surface is immersed in a computational tetrahedron mesh (a), and intersection cases between a surface and a tetrahedron (b).

To identify which tetrahedra are located around each triangle of the surface, to check for potential intersections in Step 1, and to classify the tetrahedra located at the boundary domain ∂Ω (for marking procedure to propagate from) in Step 2, we use the following approach. We first compute the bounding box of the domain. The bounding box is then subdivided into subcubes, as shown in Figure 10a. All subcubes which are incident to each triangle of the surface, see Figure 10b, are then computed. All tetrahedra which are incident to each subcube are also figured out, see Figure 10c. From these two data structures, all cut elements can be marked because one can easily access the tetrahedra around each triangle on the surface to check for intersections. On the other hand, to propagate marking procedure for outside elements, all tetrahedra located at the boundary ∂Ω can also be easily accessed from the boundary subcubes.

Figure 11 shows three type of tetrahedra marked when a liver surface is immersed in a tetrahedron mesh. inside Once all cut elements are identified, to facilitate the integration, a set of sub-tetrahedra with conforming nodes (regarding the interface surface) is embedded in each cut tetrahedron. There are only two kinds of intersection between a surface and a tetrahedron: a triangular intersection or a quadrilateral intersection, see Figure 9b. Therefore, it is sufficient to employ a set of eight tetrahedra to embed for each cut tetrahedron, see Figure 12a. This set is called the template set. Depending on each real case, where the tetrahedron is cut by a surface with a triangular intersection or a quadrilateral one, the template set is rotated and then mapped into the cut tetrahedron geometry using the mapped mesh method [START_REF] Grosland | Automated hexahedral meshing of anatomic structures using deformable registration[END_REF][START_REF] Bui | Real-time error control for surgical simulation[END_REF]. It is noted that the template nodes 4, 5, 6, 7, 8, 9, see Figure 12a, are located at the middle of their corresponding edges. Having determined the real intersections between the cut tetrahedron and the interface, as described in Section 2.3.1, if any edge is intersected by the interface surface, the corresponding midpoint is moved to match the real intersection. Once a template set is embedded for each cut tetrahedron, the last step consists in computing the relative location (inside or outside) of the subtetrahedra of the template with respect to the interface surface. For this purpose, we again use the same level-set method which was used for computing the relative location (inside or outside) of the tetrahedra. Figure 13 shows the sub-tetrahedra which are marked inside and outside with respect to the liver surface.

As described in Section 2.3.2, when an invalid cut element arises, it is recursively embedded with a tetrahedron set from a predefined template. It results in a tree data structure as schematically shown in Figure 14. To efficiently handle this type of data structure in implementation, we use a STL-like C++ tree class http://tree.phi-sci.com/. Using this container, it allows to recursively add embedded elements as children of the cut element (regarded as parent) very easily. The container also provides different kinds of iterators to access desired elements efficiently.

Numerical integration

For those elements that are fully contained in the domain Ω 1 or Ω 2 , see Figure 4a, integration of element stiffness and mass matrices are performed normally, as in classical FEM. Only for cut elements, the integration is split into two parts which are related to the inside and outside sub-tetrahedra with respect to the interface. A two dimensional schematic representation of the integration using natural coordinates on reference element is shown in Figure 15.

The stiffness matrix of the cut element reads

K e = K Ω 1 e + K Ω 2 e , (31)
where K Ω 1 e and K Ω 2 e denote the stiffness contributions of the part belonging to Ω 1 and Ω 2 , respectively, to the element e. The stiffness matrix on each part is computed by summing the contributions from their sub-elements. For Figure 14: Schematic representation of a tree data structure obtained when cut elements are embedded with a subtetrahedron set or when an invalid cut element is recursively refined to get all valid cut elements before embedding with a subtetrahedron set. The elements of the background mesh are indexed at level 0, while embedded elements are indexed at subsequent levels.

example, the stiffness matrix K Ω 2 e reads

K Ω 2 e = s∈Ω k 2 Ω k 2 B T s E 2 B s dΩ k 2 = s∈Ω k 2 Np i=1 B s (ξ i ) T E 2 B s (ξ i )ω i det(J(ξ i )), ( 32
)
in which B s is the strain displacement matrix of the sub-element s, E 2 is the material stiffness tensor of the domain Ω 2 , ξ i and ω i are the quadrature coordinates and the corresponding weight parameters, N p is the number of quadrature points used, and J is the Jacobian matrix of the coordinate transformation. Since the stiffness matrix should be expressed on the cut (parent) element where the degrees of freedom are defined, we must compute

P 3 P 1 P 2 interface n Ω 2 Ω 1 ξ η 0 1 1 ξ η 0 1 1 ξ i x i ξi x i = χ(ξ i ) ξi = χ -1 p (x i )
Natural coordinates used for the cut (parent) tetrahedron

Natural coordinates used for subtetrahedron

Figure 15: Integration on cut elements. Note that the degrees of freedom are only defined on the element nodes P 1 , P 2 and P 3 . χ and χ p are the coordinate mapping from the reference element to the physical sub-element, and to the physical cut (parent) element, respectively.

the strain displacement matrix B p of the parent element at the quadrature point ξi corresponding to the physical coordinates x i of the sub-element, see Figure 15. So, Equation (32) reads

K Ω 2 e = p∈Ω k 2 Np i=1 B p ( ξi ) T E 2 B p ( ξi )ω i det(J( ξi )), (33)
where B p is the strain displacement matrix defined on the parent element P 1 P 2 P 3 .

In this paper, we use linear tetrahedra. Therefore, the strain displacement matrix is constant across the element volume, and ω i = 1/6, det(J) = 6V k with V k the sub-element volume, see e.g. [START_REF] Dhatt | Finite Element Method[END_REF] for more details. Thereby, we get

K Ω 2 e = Ω k 2 B T p E 2 B p V k . (34)
The computation for K Ω 1 e can be done using the same concept. Also, integration over cut elements for the mass matrices is performed by the same procedure.

Corotational formulation for CutFEM

In many surgical simulations, tissues undergo large displacements and rotations, see e.g. [START_REF] Plantefève | Patient-specific biomechanical modeling for guidance during minimally-invasive hep-atic surgery[END_REF][START_REF] Bui | Controlling the error on target motion through real-time mesh adaptation: Applications to deep brain stimulation[END_REF]. Using linear elasticity for modelling of soft tissues results in artifacts for large rotational deformation [START_REF] Müller | Stable real-time deformations[END_REF]. To overcome this issue, we compute the stiffness matrix using the corotational formulation as in [START_REF] Felippa | A unified formulation of small-strain corotational finite elements: I. theory[END_REF], in which the rigid motion can be extracted from the total finite element displacements. Indeed, the system term Ku of Equation ( 17) is computed by assembling the element internal force, denoted by f i e . In the following, the superscript i is dropped for simplicity. The element displacement vector u e is computed from the rotated current configuration and the initial one as

u e = R T e x e -x 0e , (35)
where R e stands for the element rotation matrix of the element local frame with respect to its initial orientation, being updated at each time step. The corotated element internal force is then computed as

f e = R e K e u e . (36)
Using polar decomposition, the element rotation matrix R e is computed from the element deformation gradient F e as follows

R e • U e = F e , (37)
where U e is the right stretch tensor that is responsible for deformation. The element deformation gradient F e is computed as

F e = P n • N , (38)
where P n is the element nodal coordinates, and N is the derivative of shape functions. For tetrahedron element, P n and N are 3 × 4 and 4 × 3 matrices, respectively

P n =    x 1 x 2 x 3 x 4 y 1 y 2 y 3 y 4 z 1 z 2 z 3 z 4    , (39)
N =       ∂N 1 ∂x ∂N 1 ∂y ∂N 1 ∂z ∂N 2 ∂x ∂N 2 ∂y ∂N 2 ∂z ∂N 3 ∂x ∂N 3 ∂y ∂N 3 ∂z ∂N 4 ∂x ∂N 4 ∂y ∂N 4 ∂z       . ( 40
)
The corotational formulation is summarised in Figure 16. For elements cut by the surface, since degrees of freedom are only defined on their nodes (not on nodes of subelements), see also Section 2.3.5, the element rotation matrix R e is still computed by the procedure described above.

Boundary conditions on immersed surfaces

When using immersed surface, one needs to impose Neumann and Dirichlet boundary conditions on the immersed surface, which does not conform with the computational mesh. For the Neumann boundary condition, a force applied on the surface is barycentrically mapped onto the nodes of the element from the computational mesh which contains the applied point of that force, see Figure 17.

The concept of this approach is based on the master-slave scheme [START_REF] Rabczuk | Immersed particle method for fluidstructure interaction[END_REF]. In this scheme, the displacement of a point on the surface can be seen as a mapping of the displacement of the computational mesh, which is given by

u S = Ju M , (41)
where u S is the displacement of the surface (considered as slave), u M is the displacement of the computational mesh (regarded as master ), and J denotes

R T e x e -x 0e f e = K e (R T e x e -x 0e ) R e f e = R e K e (R T e x e -x 0e )

Undeformed configuration Deformed configuration

Figure 16: By using corotational formulation, the rigid body motion is removed from the element deformation.

P 3 P 1 P 2 Q P 3 P 1 P 2 Q 3 Q 2 Q 1 (a) (b)
x i Figure 17: A force Q applied on the surface at point x i (a), is barycentrically mapped into the nodes of the element from the computational mesh containing x i (b).

the matrix containing barycentric coordinates (the relative coordinates) of the considered point on the surface with respect to the computational mesh.

For example, if the slave point is located at the centroid of the triangle, the relationship described by Equation ( 41) is explicitly expressed as

u S v S = 1/3 0 1/3 0 1/3 0 0 1/3 0 1/3 0 1/3            u 1 v 1 u 2 v 2 u 3 v 3            , ( 42
)
where u and v are two components of the displacement in two dimensional space, the subscripts 1, 2, and 3 denote three nodes of the triangle. By applying the principle of virtual work

u T M Q M = u T S Q S , (43)
where Q S and Q M are the forces applied on the surface, and the equivalent forces applied at the computational mesh, respectively. By substituting Equation [START_REF] Dhatt | Finite Element Method[END_REF] into Equation ( 43), one obtains

Q M = J T Q S . (44)
Dirichlet boundary conditions at some nodes i on the surface, can be imposed as u i S = ū, where ū is the prescribed displacement. Taking into account Equation ( 41), the Dirichlet boundary condition can be expressed as

Ju M = ū. ( 45
)
Using Lagrange multipliers, one can easily impose the prescribed displacement on the computational mesh by solving the system equation set

A J T J 0 u M λ = b ū , ( 46
)
where λ stands for Lagrange multipliers used for Dirichlet boundary conditions.

Solving system equations with constraints

To solve the system equations with constraints, there are several prominent approaches, each with their own advantages and disadvantages. The penalty method, where the respective matrix entries are multiplied by some large value, suffers from the uncertainty of the penalty parameter (weight). By using a very small value, the constraints may not be accurately imposed, and using a high value may lead to an ill-conditioned system (and thus directly affecting the computational time). The Nitsche's method is problem-dependent since it requires a penalization parameter for stability. Compared to these two approaches, the Lagrange multiplier method offers some advantages. It offers excellent generality and accuracy (no dependence on the problem), and it is almost insensitive to the user's choice (no parameter to choose). However, this method is no panacea either. Using Lagrange multiplier method, we have to solve for additional unknowns. Nevertheless, in our implementation for Lagrange multiplier method, we do not solve the system matrix with constraints by inverting the whole system (augmented with Lagrange multipliers) since it renders the stiffness matrix indefinite. Instead, we solve the constraints in three steps, as described in the following: i) solve the unconstrained system, ii) solve for constraints, and iii) update the solution when Lagrange multipliers are available. These steps are detailed in Algorithm 1. Note that, solving the system equations by using this approach is also relevant for parallelization, which is crucial for real-time simulations. We also remark that we do not employ any stabilization technique in the proposed Lagrange multiplier's approach, such as that used in [START_REF] Burman | Fictitious domain finite element methods using cut elements: I. A stabilized Lagrange multiplier method[END_REF]. As stated in Section 2.2, for the tissue we employ C 0 -nonconforming three-dimensional linear elements, and for the needle we utilise C 1 -conforming one-dimensional linear beam elements (see the end of Section 2.2). The numerical results presented in Section 3, and those presented in [START_REF] Bui | Real-time error control for surgical simulation[END_REF][START_REF] Bui | Controlling the error on target motion through real-time mesh adaptation: Applications to deep brain stimulation[END_REF], show that the proposed Lagrange multiplier approach, together with this choice of spaces, delivers converging results, and that it does not suffer from any locking or instability issues. We intend to conduct a rigorous mathematical analysis of the proposed method in our future work.

The interaction between the needle (denoted by subscript 1) and the tissue (denoted by subscript 2) can be expressed by the following equation set

    A 1 0 H T 1 0 A 2 H T 2 H 1 H 2 0            dv 1 dv 2 λ i        =        b 1 b 2 0        , ( 47
)
where λ i is the Lagrange multiplier representing the interaction between the needle and the tissue. We can see that, Equation [START_REF] Anitescu | Recovery-based error estimation and adaptivity using high-order splines over hierarchical T-meshes[END_REF] describing the interaction between the needle and the tissue, and Equation ( 46) expressing the constraints used for Dirichlet condition on implicit boundaries, have the same general form

A J T J 0 x λ = b c . ( 48
)
Equation ( 48) can be reformulated as

x = A -1 b x f ree -A -1 J T λ xcor , (49a)
JA -1 J T W λ = J A -1 b x f ree -c, (49b)
in which, x f ree can be seen as the solution of the unconstrained system Ax = b, x cor is the corrective solution due to constraints, which can be seen as the solution of the system Ax cor = J T λ. Here W denotes the compliance matrix. As described in [START_REF] Bui | Controlling the error on target motion through real-time mesh adaptation: Applications to deep brain stimulation[END_REF], to be able to solve the Equation (48) in real time, we use the advanced method described in Algorithm 1.

Algorithm 1 Algorithm for solving the system equation with constraints (Equation ( 48)). L denotes a lower unit triangular matrix, and D is a diagonal matrix.

1: LDL T ← decompose(A) (Cholesky decomposition) 2: x f ree ← solveCholesky(Ax f ree = b) (Solve unconstrained problem) 3: λ ← solve(Wλ = Jx f ree -c) (Solve for constraints, see (49b)) 4: x cor ← solveCholesky(Ax cor = J T λ) (Compute corrective solution) 5: x ← x f ree -x cor (Update solution, see (49a))
Therefore, three main steps for solving Equation (48) can be summarised as

Step 1: Decompose the matrix A, and solve for x f ree , described by Lines 1 and 2 of Algorithm 1,

Step 2: Use x f ree , solve for Lagrange multipliers λ described in Line 3,

Step 3: Solve for the corrective solution x cor , and update the solution x, described in Lines 4 and 5.

It is noted that the compliance matrix W = JA -1 J T is computed from the decomposition LDL T of the matrix A, and this computation is parallelized on GPU, as proposed in [START_REF] Courtecuisse | Real-time simulation of contact and cutting of heterogeneous soft-tissues[END_REF][START_REF] Bui | Controlling the error on target motion through real-time mesh adaptation: Applications to deep brain stimulation[END_REF].

Results

In this section, we present numerical results with the following goals: 1. To measure the convergence of the proposed implicit boundary method as compared to its conforming counterpart (classical FEM) in terms of the ability to reproduce the geometry of the organ/domain. This is the focus of Section 3.1, and partly of Section 3.2.

2. To measure the ability of the implicit boundary method in dealing with boundary conditions which are not imposed upon the mesh, but upon an implicitly defined boundary. This is the focus of Section 3.2-Section 3.4.

Convergence study

Before proceeding to study the convergence of the numerical solution obtained using CutFEM, we first study how accurately the non-polygonal geometry is approximated in CutFEM. In order to measure the ability of the method to represent the exact geometry, we compute the volume of the geometry discretised by CutFEM. Since the input of the geometry of the organ, e.g., liver or brain, is given in a discretised form, we compare the volume obtained using CutFEM against the volume obtained from the discretised surface. We consider the discretised surface of a spherical geometry with the diameter d = 1.4. As shown in Figure 18, the volume error, between the volume obtained from the discretised surface (V e = 1.40005) and the volume V of the sphere discretised by CutFEM (based on the embedded discretised surface), converges with an optimal rate, which is approximatively the same rate as the L 2 norm of the displacement error.

In order to verify the reliability of the solution obtained using CutFEM, now we study the convergence of the solution under mesh refinement. We consider two cases, a tensile test and a bending test. As shown in Figure 19, a beam is studied with a spherical surface being embedded inside. As mentioned above, the spherical surface, with the diameter d = 1.4, is given in the discretised form. At one end, the beam is subjected to a uniform horizontal pressure (tensile) or a uniform vertical pressure (bending), while the other end is clamped. For the convergence study, the linear elastic constitutive law is used. The mechanical properties, namely, the Young's modulus and the Poisson's ratio, are denoted by E and ν, respectively. For the material outside the sphere we denote them by E 1 , ν 1 , and for the material inside the sphere we denote them by E 2 and ν 2 . The dimension of the beam is 6 × 2 × 2 mm. In order to compare the convergences rate of the tests with the theoretical ones, we use the same mechanical properties for the material inside and outside the sphere surface. We set thus E 1 = E 2 = 1000 MPa, and

ν 1 = ν 2 = 0.1.
The convergence is studied by computing the solution of the tensile and bending tests employing the tetrahedral meshes consisting of 7 × 3 × 3, 13 × 5 × 5, 25 × 9 × 9, and 49 × 17 × 17 nodes. We propose to use the solution from the classical FEM when employing a very fine mesh (97 × 33 × 33 nodes) as the reference solution. We then study the convergence of the error between the CutFEM solution and the reference FEM solution. This error is measured by using both, the L 2 norm and the energy norm. The L 2 norm of the displacement error reads

η L 2 = Ω (u h -u r ) 2 dΩ Ω u 2 r dΩ , (50)
where u h denotes the displacement solution of the CutFEM, and u r denotes the reference solution of the classical FEM. The energy norm is defined as

η Energy = Ω (σ h -σ r ) • ( h -r ) dΩ Ω σ r • r dΩ , (51)
where σ h and h denote the stress and strain of the CutFEM, respectively, and σ r and r denote their reference values obtained from the classical FEM on a very fine mesh. For tensile and bending tests, Figure 20 shows the convergence rates in the L 2 norm and the energy norm versus the number of DOFs. It is observed that for both tests, the rates of convergence of the L 2 norm and the energy norm agree well with the theoretical rates. Indeed, for 3D problems using linear elements, the L 2 norm of the displacement error is of order O(N -2/3 ), while the energy norm converges with an order of O(N -1/3 ), where N denotes the total number of the degrees of freedom.

A comparison with FEM

The aim is to compare displacement simulation results obtained from the CutFEM using non-conforming meshes with those obtained from the classical FEM using conforming meshes. We carry out the study on two different geometries: the simple beam geometry with an immersed sphere as shown in Figure 19, and the more complex liver geometry as shown in Figure 9a.

For the beam geometry, the dimensions and the mechanical properties are same as those used in Section 3.1. The displacement measured at the centre of the left end of the beam is employed to compare between the CutFEM (where the sphere is modelled implicitly) and the classical FEM. The beam is subjected to a uniformly distributed pressure q = -5 N/mm 2 in the vertical direction at the left end, whereas it is clamped at the right end. Under various mesh refinements, Figure 21 shows the displacement of the point during simulations using the CutFEM and the classical FEM. It is observed that, with the same number of degrees of freedom used, the results obtained from the CutFEM perfectly agree with that of the classical FEM. Also, under mesh refinement, the displacement asymptotically convergences to the solution of the fine mesh.

For the liver, due to its complex geometry, in order to apply the same boundary conditions acting on different conforming meshes (used for FEM), and on different non-conforming ones (used for the CutFEM), a homogeneous Dirichlet boundary condition is implicitly applied to the model through the points located on an imaginary cutting section, shown by the points in blue colour in Figure 22, whereas a uniformly distributed pressure q = -1 N/mm is implicitly applied to the model in the vertical direction through the mesh shown by green colour in Figure 22. The displacement is measured at the point inside the liver shown by the grey colour in Figure 22. For these simulations, Young's modulus and Poisson's ratio are 1000, and 0.4 respectively. Figure 23 shows the displacement of the point shown in Figure 22. It is observed that, at the same mean size of the elements used, the results obtained from the CutFEM agrees well with those of the classical FEM. It is noted that the number of degrees of freedom is not used as the same input during the comparison between the CutFEM and the classical FEM since, with respect to the liver geometry, it does not characterise the same mesh resolution between conforming mesh (used for FEM) and the non-conforming mesh (used for CutFEM). Instead, the mean element size l is employed (see Figure 23).

We now briefly discuss the computational cost of CutFEM vs FEM. We consider the liver model for this comparison. For the simulation performed by FEM using a mesh with the mean mesh size l = 0.00635283, compared to that performed by CutFEM using a mesh with l = 0.0059953, we observe a very good agreement on the displacement, see Figure 23. However, the use of CutFEM is more advantageous than the use of FEM since the computational times needed to solve the system equation by using FEM and CutFEM are 983.41 ms and 429.85 ms, respectively. Thus, for this mesh size, CutFEM is faster than FEM by a factor of approximately 2.3.

Needle insertion simulations

Having established the accuracy and reliability of CutFEM simulations as compared to the classical FEM, we now employ the CutFEM approach for needle insertion problems.

Immersed interface

The needle, which is initially inclining at an angle of 3.5 degrees, is inserted into a phantom tissue with a simple geometry, as shown in Figure 24. We simulate a spherical inclusion (can be considered as a tumor) which is implicitly immersed in the tissue phantom. The dimension of the tissue phantom is 6 × 2 × 2, while the radius of the inclusion is 0.7. The length of the needle is 2.8, and its cross section radius is 0.05. Young's modulus and Poisson's ratio of the needle is set to 20 000 and 0.2, respectively. These parameters for the phantom tissue are E 1 = 1 000 and ν = 0.4. The Poisson's ratio of the inclusion is also set to 0.4, however, in order to investigate the effect of the inclusion stiffness on the needle-tissue interaction force profile, the Young's modulus of the inclusion E 2 is varied with respect to that of the phantom tissue by a factor of 1, 2, 4 and 8. The penetration strength at the tissue surface is set to 1, and the frictional coefficient between the tissue and the needle shaft is set to 0.5.

Figure 25 shows the needle-tissue interaction force with respect to the displacement of the needle tip, with different ratios E 2 /E 1 . It is observed that when the needle tip reaches the tissue surface, the interaction force between the needle and the tissue occurs. This interaction force continuously increases and Force Displacement

E 2 = E 1 E 2 = 2E 1 E 2 = 4E 1 E 2 = 8E 1
Figure 25: Force displacement curves, with varying Young's modulus ratio E 2 /E 1 , when the needle is inserted into the soft tissue model. When the needle tip displacement reaches about 3, the needle is retracted, and thus generating negative needle-tissue interaction force; at the stage when the needle is completely retracted from the tissue, the interaction force vanishes.

when it reaches the tissue surface penetration strength, the needle penetrates into the tissue. It also reveals that the closer to the inclusion the needle tip is, the more different the interaction force profiles are obtained when the ratio E 2 /E 1 is varied. This is indeed logical due to the higher stiffness of the inclusion as compared to that of the tissue. When the displacement of the needle tip reaches 3, the needle is continuously retracted until completely out of the tissue. During retraction process, the interaction force changes its sign, and is negative, as observed in Figure 25. The same conclusion about the effect of inclusion-tissue stiffness ratio on the interaction force profile can be drawn during the retraction phase as during the insertion stage.

In Figure 26, we present the results of mesh refinement study during needle-tissue interaction, while keeping the ratio E 2 /E 1 = 8 fixed. The curve profiles are qualitatively the same as of those curves presented in Figure 25. Moreover, it shows that when the mesh is refined, the proposed Lagrange multiplier method delivers convergent results (the difference between the red and green curves is smaller than that between the blue and red curves).

Fictitious boundary

We now consider a comparison between a needle insertion simulation into a liver using implicit/fictitious boundaries (CutFEM), and with those using explicit/conforming boundaries (classical FEM). The goal is to show the performance of the CutFEM approach in dealing with Dirichlet boundary conditions defined implicitly on the surface, and in dealing with the interaction between the needle tip and the implicit surface (a kind of Neumann boundary conditions), as compared to those dealt explicitly by FEM. For the liver, Young's modulus and Poisson's ratio are set as 1 200 kPa and 0.4, respectively. For the needle, the same parameters as above are used. The background mesh used for simulation with CutFEM, and the boundary conditions applied to the liver surface, are shown in Figure 27. The point, where the displacement is measured during the insertion and the retraction of the needle, is also shown in Figure 27. As can be seen in Figure 28, the displacement results obtained by CutFEM agree well with those obtained by classical FEM.

Electrode implantation simulation in Deep Brain Stimulation

The CutFEM is now employed to simulate an electrode lead implantation, using in Deep Brain Stimulation (DBS) procedure. We also take into account the brain shift phenomenon, which occurs due to the leak of cerebro-spinal fluid when a craniotomy is performed. The goal of the simulation is to insert an electrode inside the brain until it reaches the target, e.g., subthalamic nucleus (STN) area for treatment of Parkinson's disease. To do that, a cannula is inserted together with the electrode lead through a hole drilled in the skull. When they reach the STN area, the cannula is retracted while keeping the electrode lead inside. As in [START_REF] Bui | Controlling the error on target motion through real-time mesh adaptation: Applications to deep brain stimulation[END_REF], frictional interactions between the brain tissue with the cannula and electrode lead are simulated. Young's modulus of 6 kPa and Poisson's ratio of 0.45 are set to the brain tissue. The cannula and electrode lead are set with Young's modulus of 10 GPa, and with Poisson's ratio of 0.3. The input background mesh used for the CutFEM simulation of brain behaviour is shown in Figure 29a. We consider simple boundary conditions for the brain tissue. Indeed, brain tissue around the optic nerves and the brainstem are considered to be clamped. Moreover, bilateral interaction constraints are considered between the brain surface and the skull. It is noted that these constraints are implicitly applied on the brain surface, as described in Section 2.5. During simulation, the displacement is measured at the point in the STN area as shown in Figure 29b. Figure 30 shows the brain deformation at different stages of the simulation. The horizontal lines in Figure 30 help to show the differences of the brain deformation due to brain shift, and due to cannula insertion. The displacement of the STN target due to brain shift, and due to cannula insertion and full retraction, is presented in Figure 31. It reveals that the brain shift is the origin of the STN displacement which is stabilised about 1.1 cm. When the cannula is inserted inside the brain tissue, due to frictional interaction between them, the displacement of the STN target increases. The cannula is undergoing retraction immediately after the cannula tip has reached the STN target. This induces the decreasing of the STN displacement before the STN target is stabilised around the location, where the STN was found after the brain shift stage. These results show that the CutFEM, even with a nonconforming mesh, is able to simulate the behaviour of the brain while it is in interactions with the skull, the cannula, and the electrode lead. This makes the discretisation as independent as possible from the geometric description. We believe that with such a tool, the patient-specific simulations (in terms of geometries, for instance) can be performed, particularly in the context of real-time simulations.

Conclusions

A corotational formulation of the CutFEM has been proposed. Using Lagrange multipliers, we have shown the methods to implicitly apply Dirichlet boundary conditions on the immersed surface, which is not conformed to the background mesh used for the simulations. We verified the implementation by studying convergences through a tensile test and a bending one. We also demonstrated the performance of the CutFEM, compared to the classical FEM, through needle insertion problems in medical simulations (e.g. for biopsy or brachytherapy). We have shown that the interaction between the tissue surface being implicitly defined, and the needle tip, to account for the penetration phenomenon of the latter to the former (a kind of Neumann boundary conditions) during simulations, is working properly, and refers to the real-life situation. Also, by employing the CutFEM, it makes the discretisation as independent as possible from the geometric description, and also minimises the complexity of mesh generation, especially for complex geometries. Because of the constraint of computational time in real-time simulation context, while still preserving the geometric details, the use of coarse meshes becomes possible by using CutFEM. Two kinds of applications using the CutFEM have been studied in this paper: i) immersed interfaces, and ii) fictitious boundaries. The immersed interfaces are useful for simulations of heterogeneity of tissues, e.g., when tumors or internal structures of tissues are considered. Simulations using fictitious boundaries are suitable for applications in which the tissue geometries are complex, and by using the CutFEM, it is possible to integrate only the material inside the tissue surface while using a nonconforming mesh to its boundaries.

In the presented work, for the discretisation of simulated domains, we limited ourselves to linear elements. Another limitation is the embedding of linear subelements for the cut elements to facilitate numerical integration. By doing this, we still cannot precisely capture the surface geometries. For more accurate integration of implicit geometries, one can use higher order elements, as proposed in [START_REF] Moumnassi | Analysis using higher-order xfem: implicit representation of geometrical features from a given parametric representation[END_REF][START_REF] Fries | Higher-order accurate integration of implicit geometries[END_REF]. We note that the isogeometric analysis, in the spirit of integrating finite element analysis and CAD design tools by using higher-order spline/NURBS basis functions from CAD geometry, could also be used for surgical simulations. We intend to study this approach in our future work. The interested readers are referred to [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF] for pioneering work on isogeometric analysis, and a very recent work [START_REF] Anitescu | Recovery-based error estimation and adaptivity using high-order splines over hierarchical T-meshes[END_REF] on recovery-based error estimation and adaptivity over hierarchical T-meshes.

Chouly, Prof. Alexei Lozinski, and Prof. Stéphane Cotin for many helpful discussions.

Figure 1 :

 1 Figure 1: Two dimensional representation of the problems studied in the paper.

Figure 3 :

 3 Figure 3: Schematic representation of contact between the needle (shown by red colour) and the tissue surface (shown by black colour).

Figure 4 :

 4 Figure 4: Two dimensional representation of the problem domains (a), the domains are nonconformingly discretised with respect to the surface Γ (b), the cut elements are highlighted by yellow colour whereas all inside elements of the volume Ω 2 are highlighted by red colour, the remaining elements are marked as outside (i.e. outside of Γ) elements (c).

Figure 5 :

 5 Figure 5: Schematic representation of intersection computation between a triangle and a discretised surface. An embedded subtriangulation is performed on the element for integration purpose.

Figure 6 :

 6 Figure 6: Schematic representation of invalid cut elements in two dimensions: (a), (b), (c); refinements lead to valid cut elements: (d), (e), (f). Valid cut elements are embedded with a subtetrahedron set (shown in green colour) for integration purpose: (g), (h), (i).

Figure 7 :

 7 Figure7: A sphere is immersed into a background mesh (a), an invalid cut element is embedded with sub-tetrahedra after refinement (b), and finally all, invalid and valid, cut elements, intersected with the spherical surface, are embedded with sub-tetrahedra, with or without refinement, respectively, for integration reason (c).

Figure 8 :

 8 Figure 8: When an element is cut by the surface with a small intersection (a), the concerned node is moved by a distance proportional to the element size (b).

Figure 10 :

 10 Figure 10: Subdivision of the bounding box (a), subcubes incident to each triangle of the interface are computed (b), elements incident to a subcube are marked (c).

Figure 11 :

 11 Figure 11: Three type of tetrahedra are marked when a liver surface is immersed in a computational tetrahedron mesh: cut element (a), inside elements (b) and outside elements (c).

Figure 12 :

 12 Figure 12: A set of eight template tetrahedra (a). Depending on whether the intersection between a surface and a cut tetrahedron is a triangular section or a quadrilateral one, the template is rotated and then mapped in the cut tetrahedron geometry (b).

Figure 13 :

 13 Figure 13: Once a template tetrahedron set is embedded for each cut tetrahedron, using the level-set method, sub-tetrahedra are marked as inside (a) or outside (b) with respect to the interface surface.

Figure 18 :

 18 Figure 18: Convergence of the volume of the sphere, discretised by CutFEM from the embedded discretised surface, as compared to the exact volume obtained from the discretised surface (V e = 1.40005), under background mesh refinement.

Figure 19 :

 19 Figure19: A beam is clamped at the right end (marked by the red points), and is subjected to, at the left end (marked by the green triangles), a uniform horizontal pressure (the tensile case) or a uniform vertical pressure (for bending case). A sphere surface is immersed inside the beam geometry.

Figure 20 :

 20 Figure 20: Convergence rates under mesh refinement for the tensile test (a), and for the bending test (b).

Figure 21 :

 21 Figure 21: Displacement (in mm) measured at the centre of the left end of the beam during simulations by the CutFEM and the classical FEM.

Figure 22 :

 22 Figure 22: Liver is implicitly clamped at the points shown in blue colour, and is implicitly subjected to a uniformly distributed vertical pressure acting on the mesh shown in green colour. During the simulations, the displacement is measured at the point shown by the grey colour.

Figure 23 :

 23 Figure 23: Displacement measured (in cm) at the point shown in Figure 22 during simulations by the CutFEM, and by FEM. Different mesh resolutions are denoted by different values of l.

Figure 24 :

 24 Figure 24: Schematic illustration of a needle insertion simulation into a simple tissue geometry. The needle is initially inclining at an angle of 3.5 degrees.

Figure 26 :

 26 Figure26: Force displacement curves, with varying mesh size (coarse mesh: 498 DOFs, finer mesh: 1575 DOFs, even finer mesh: 3183 DOFs) but for fixed Young's modulus ratio E 2 /E 1 = 8, when the needle is inserted into the soft tissue model.

Figure 27 :

 27 Figure 27: Background mesh used for the simulation of the liver behaviour during needle insertion and retraction. Some constraints, implicitly defined on the liver surface using Lagrange multipliers, are shown by the blue points. During the simulation, the displacement is measured at a point located near the needle shaft.

Figure 28 :

 28 Figure 28: Displacement measured at the point during insertion and retraction, using CutFEM and standard FEM. The vertical dash line indicates the moment when the needle is starting to be retracted.

Figure 29 :

 29 Figure 29: Input background mesh used for simulation of brain behaviour (a); during simulation, the brain has interaction with the skull, and the displacement is measured at the point shown in green colour (b).

Figure 30 :

 30 Figure 30: Brain at initial state (a), brain deformation due to brain shift (b), brain deformation due to cannula insertion (c). The blue horizontal line helps to visualise the brain deformation due to brain shift occurring after a craniotomy, whereas the green line helps to visualise the brain deformation due to the cannula insertion.

Figure 31 :

 31 Figure 31: Displacement of the STN target due to brain shift, and due to cannula insertion and retraction.

https://www.sofa-framework.org

Acknowledgements

Stéphane Bordas, Satyendra Tomar and Huu Phuoc Bui thank partial funding for their time provided by the European Research Council Starting Independent Research Grant (ERC Stg grant agreement No. 279578) RealTCut "Towards real time multiscale simulation of cutting in non-linear materials with applications to surgical simulation and computer guided surgery". We are also grateful for the funding from the Luxembourg National Research Fund (INTER/FWO/15/10318764).

The funding from the University of Strasbourg Institute for Advanced Study (BPC 14/Arc 10138) for the first author is gratefully acknowledged.

Huu Phuoc Bui would particularly like to thank Dr. Davide Baroli, Dr. Lionel Untereiner, Dr. Paul Hauseux, Prof. Thomas-Peter Fries, Dr. Franz