<figure><img src="image_1.png" /></figure>

HAL
open science

# Haptic Rendering of Hyperelastic Models with Friction

Hadrien Courtecuisse, Yinoussa Adagolodjo, Hervé Delingette, Christian Duriez

► To cite this version:

Hadrien Courtecuisse, Yinoussa Adagolodjo, Hervé Delingette, Christian Duriez. Haptic Rendering of Hyperelastic Models with Friction. 2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Sep 2015, Hamburg, Germany. pp.591-596, 10.1109/IROS.2015.7353432. hal-01184113

HAL Id: hal-01184113

https://hal.science/hal-01184113

Submitted on 14 Oct 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Haptic Rendering of Hyperelastic Models with Friction

Hadrien Courtecuisse<fnref n="1" />, Yinoussa Adagolodjo<fnref n="1" />, Hervé Delingette<fnref n="3" />, Christian Duriez<fnref n="2" />

**Abstract**— This paper presents an original method for inter-
actions' haptic rendering when treating hyperelastic materials.
Such simulations are known to be difficult due to the non-linear
behavior of hyperelastic bodies; furthermore, haptic constraints
enjoin contact forces to be refreshed at least at 1000 updates per
second. To enforce the stability of simulations of generic objects
of any range of stiffness, this method relies on implicit time
integration. Soft tissues dynamics is simulated in real time (20
to 100 Hz) using the Multiplicative Jacobian Energy Decomposi-
tion (MJED) method. An asynchronous preconditioner, updated
at low rates (1 to 10 Hz), is used to obtain a close approximation
of the mechanical coupling of interactions. Finally, the contact
problem is linearized and, using a specific-loop, it is updated
at typical haptic rates (around 1000 Hz) allowing this way new
simulations of prompt stiff-contacts and providing a continuous
haptic feedback as well.

I. INTRODUCTION

During the past decade simulators have become more
and more important in the medical field, especially for
their applications in educational and learning processes for
surgical schools. In the context of minimally invasive pro-
cedures, where using the instrument is the only way for
the surgeon to be in contact with anatomical parts, it is
essential for simulators to provide an extreme realistic haptic
feedback for each single action. The realistic haptic rendering
of interactions between rigid instruments and deformable
organs is particularly an issue due to several reasons: one
above all the strongly non-linear behavior of soft organs.
Furthermore, from the computational point of view, it is
necessary to update contact forces information as rapidly
as possible, since it is common knowledge that providing
a smooth haptic feedback requires at least 500 Hz. Finally,
in the typical simulation environment, rigid interactions on
soft bodies cause stability issues.

In this paper we propose a new solution to simulate a
realistic haptic feedback of interactions with hyperelastic
soft tissues. Implicit time integration is used to enforce the
stability of the simulation in case of arbitrary contacts with
objects of any range of stiffness. Such integration scheme
allows a relatively large time step but it requires to addi-
tionally evaluate a global matrix and the solution of a linear
system of equations for each time step. The Multiplicative
Jacobian Energy Decomposition (MJED) [18] algorithm is
used to provide an implementation of hyperelastic materials
in real time. This method discretizing non-linear hyperelastic

materials on linear tetrahedral meshes leads to a faster
stiffness matrix assembly for a large variety of isotropic
and anisotropic materials (Costa, Veronda Westmann, Boyce
Arruda, StVenant Kirchhoff, NeoHookean, Ogden, Mooney
Rivlin). An asynchronous preconditioner is updated at low
frequency (5-10 Hz) [9] not only to accelerate the conver-
gence of iterative solvers but also to evaluate the compliance
of the system (i.e. the mechanical coupling between the
contacts on deformable bodies). Finally the contact problem
is updated (assuming a linearization) and solved at haptic rate
(1000 Hz) [11] to provide a continuous and realistic haptic
rendering. The method allows to simulate prompt interactions
as well as non-linear friction problems.

II. RELATED WORKS

Finite element methods (FEM) provide high bio-
mechanical realism, mainly because the soft-tissues complex
non-linear behaviour is directly explained through constitu-
tive relations. In [6], real-time computations for linear elastic
materials are achieved by precomputing Greens functions in
order to accelerate online simulations models. Even though
it allows to drastically decreasing computational expenses,
linear models are limited to small strain and small displace-
ments assumption. The *co-rotational* was introduced by [12]
within the field of numerical methods. In this formulation, the
stiffness of each element is assumed linear within the local
frame, described by its rotated state. This approach allows
simulating geometric non-linearities (i.e. large displacements
and rotations) while maintaining the smallest algorithmic
complexity. However, it is restricted to small strain ranges
excluding important physical aspects such as volume preser-
vation during deformations.

The *Total Lagrangian explicit dynamics* (TLED) [20]
allows modeling both geometric and material non-linearities
in real time. The main limitation for this method is to be
based on an explicit time integration scheme, since it entails
very short time steps in order to keep the computation stable,
especially for stiff materials. [18] proposed the *Multiplicative
Jacobian Energy Decomposition* (MJED): a general algo-
rithm to implement hyperelastic materials based on TLED
with implicit time integration schemes.

Haptic rendering enables physical interaction with simu-
lated objects in a virtual environment. Somehow, it represents
a means to evaluate the level of accuracy of physics-based
models in an interactive simulation. However, this link be-
tween the quality of the haptic feedback and the accuracy of
interactive simulation is not direct. Indeed, the challenge is
to maintain the same quality in modelling, given that real-
time constraints of haptic feedbacks are stricter than the ones

This work was in part supported by the project RESET(ANR-14-CE17-0007), grant ANR-2014 of the French Agence Nationale de la Recherche.

<footnote n="1">¹AVR Team-Project, CNRS Strasbourg and Strasbourg University</footnote>

<footnote n="2">²Shacra Team-Project, INRIA Lille - Nord-Europe and Lille University</footnote>

<footnote n="3">³Asclepios Team-Project, INRIA Sophia Antipolis</footnote>
of models used in the visualization of the simulation.

In the context of soft-tissues simulations with haptic feedback, deformable models and haptic renderings share the same fundamentals: in both cases, mechanical forces need to be computed. The manipulated objects interact with each other generating reaction forces which affect the physics-based model, and more importantly the controller of the haptic display. In this last case, the computation of the reaction forces is limited by (i) the need of high frequency refresh rates, (ii) the stability of the device control law, and (iii) the fidelity of the haptic rendering as perceived by the user. Several approaches have been proposed to deal with some or all of these constraints.

In [7], the concept of *virtual coupling* is proposed. Parameters can be tuned to guarantee the stability of haptic interactions but it often introduces ghost-damping forces. A *god-object* approach is introduced in [30], to increase the visual perception of stiffness, enforcing non-penetration. This method is extended by so called *virtual proxy* in [26]. These methods can be gathered as *virtual coupling network* as noted in [3] and [2] to further improve both the stability and the performance of the rendering. Extensions of the god-object to rigid body with 6DoF (constraining position and orientation) is proposed in [22].

In the pioneering works [6], [8], a *displacement-driven* interaction is used instead of a contact modelling. Positions are applied as bilateral constraints (equality conditions) and solved by Lagrange multipliers method. The approach was originally based on a superposition principle, and was further extended [24] using a force extrapolation method. Other methods, such as a finite element model handling geometric non-linearities is proposed in [29]. The model employs mass-lumping and explicit time integration for real-time simulation of dynamic behavior. Due to explicit integration, rendering is possible only on very soft models being computed on coarse meshes.

Similarly, the *small area paradigm* relies on linear modelling [25], for which the equality of boundary conditions affects only a small number of surface nodes. The response forces can be computed at each step of the haptic loop, through a simple update of inverse stiffness matrix. In [17], [28] and [23], precomputations based methods are proposed; stable haptic rendering is efficiently achieved calculating response forces from *precomputed data*, e.g. by interpolations performed directly in the haptic loop. Although visco and hyperelastic models are employed, only point-based interaction is taken into account due to the limitations given by the precomputations.

[5] propose a scheme to simulate forces reflecting deformable objects, based on the simultaneous computation of forces and displacements on the physics-based model. The computation of the FEM model is shortened using spectral Lanczos decomposition method. In [16], a unified approach to the interaction with elastostatic contacts simulation is presented. In this case, contact resolution is based on *calapitance matrix* which relates imposed displacements and response forces. The method is used for single-point as well

as grasping interactions, where haptic rates are achieved using precomputed Green functions. In [4], the contact problem is solved using a *penalty-based method* allowing for multiple contacts and self-collisions. Both the tool and the obstacle are deformable and simulated by optimized finite elements methods thanks to a model reduction. However, penalty-based methods cannot guarantee non-interpenetration: they are very sensitive to the choice of penalty parameters and they do not integrate properly static/dynamic friction. Furthermore, in the context of haptic rendering, it could lead to additional stability problems.

Therefore, based on an example of virtual snap-in simulations (where stiff contacts and deformations are closely linked), [10] introduced Signorini's model for contact handling in the field of haptic rendering. In [11], an extension to friction contact response is proposed. The computation performance is obtained through the use of a precomputed and condensed compliance matrix. Models are limited to small displacements deformations, which is not realistic for soft-tissues models.

To overcome the limited refresh rate when dealing with deformable bodies, one of the strategies is to implement an intermediate representation of the constraints provided by the virtual environment, see e.g. [1], [19] or [14]. This representation allows to separate the haptic rendering from the physics engine, accepting the fact of using a simplified model for the simulation [13]. In [15] the simplified model relies on the linearization of non-linear deformable objects that runs at low rates.

When trying to combine constraint-based approaches on deformable models with intermediate representations to handle non-linearities, the main issue is to solve constraints on a reduced model. This problem brings to the condensation of the mechanical behaviour in the constraint space. For instance, in [14], a mixed-LCP formulation is used to solve both unilateral and bilateral constraints. The formulation relies on a simplified inverse matrix similar to the compliance matrix, since only diagonal blocks are taken into account. Nevertheless, indirect contacts cannot be rendered and computation of direct contact forces is not always accurate, especially for light or stiff objects.

On the contrary, the method presented in [27], which is extended in this paper, uses an asynchronous strategy between simulation and haptic rendering, sharing the full compliance matrix in the contact space. However, a precomputation of the compliance is used, preventing from any use of hyperelastic models. Thanks to GPU, it has been shown in [9], that a very good approximation of this compliance matrix can be obtained, even on non-linear models, computing a preconditioner at low rates. This paper shows that we can combine these two last methods with MJED method.

### III. TIME INTEGRATION AND BIOMECHANICAL MODEL

The dynamic equation of simulated bodies can be written using the synthetic formulation, given by Newton's second law:

$$M(\mathbf{q}) \ddot{\mathbf{q}} = \mathbf{P}(t) - \mathbf{F}(\mathbf{q}, \dot{\mathbf{q}}) + \mathbf{R}(\mathbf{q}, \lambda) \quad (1)$$
where **q** is the vector of generalized degrees of freedom, M(**q**) is the inertia matrix (assumed constant and noted M in the rest of the paper). F gives the internal forces according to the position **q** and velocities **q̇** of the deformable objects. P represents the external forces (such as gravity) at time t. R is a non-linear function that gathers constraint forces. These parameters depend on both geometrical aspects (related to the position **q**), and mechanical properties such as intensity of contact forces (given by the Lagrange multipliers λ in the constraint space).

### A. Implicit time integration

Using a backward Euler implicit integration and considering the time interval $[t_a, t_b]$ of length $h = t_b - t_a$, the equation (1) can be rewritten:

$$ 
\begin{aligned}
M(\dot{\mathbf{q}}_b - \dot{\mathbf{q}}_a) &= h(\mathbb{P}(t_b) - \mathbb{F}(\mathbf{q}_b, \dot{\mathbf{q}}_b)) + h\mathbb{R}(\mathbf{q}_a, \boldsymbol{\lambda}_b) \\
\mathbf{q}_b &= \mathbf{q}_a + h\dot{\mathbf{q}}_b
\end{aligned}
\quad (2) $$

Since F is a non-linear function, we perform a Taylor series expansion making a first order approximation:

$$ \mathbb{F}(\mathbf{q}_a + d\mathbf{q}, \dot{\mathbf{q}}_a + d\dot{\mathbf{q}}) = \mathbf{f}_a + \frac{\partial \mathbb{F}}{\partial \mathbf{q}} d\mathbf{q} + \frac{\partial \mathbb{F}}{\partial \dot{\mathbf{q}}} d\dot{\mathbf{q}} \quad (3) $$

where **B** = $\frac{\partial \mathbb{F}}{\partial \dot{\mathbf{q}}}$ and **K** = $\frac{\partial \mathbb{F}}{\partial \mathbf{q}}$ are known respectively as the damping and the stiffness matrices. Replacing (3) in (2) and using $d\mathbf{q} = \mathbf{q}_b - \mathbf{q}_a = h\dot{\mathbf{q}}_b$ and $d\dot{\mathbf{q}} = \dot{\mathbf{q}}_b - \dot{\mathbf{q}}_a$, we obtain:

$$ \left( \frac{1}{h} \mathbf{M} + \mathbf{B} + h\mathbf{K} \right) d\dot{\mathbf{q}} = -h\mathbf{K}\dot{\mathbf{q}}_a - \mathbf{f}_a - \mathbf{p}_b + \mathbb{R}(\mathbf{q}_a, \boldsymbol{\lambda}_b) \quad (4) $$

where $\mathbf{p}_b$ is the value of function P at time $t_b$. R is evaluated with a proximity-based collision detection (see [11] for details) using $\mathbf{q}_a$. It provides a set of potential contacts with the respective contact directions **n**. To simplify the solution process, we assume that contact directions **n** do not change during the contact response. These normal are collected in a matrix **H** (known as the Jacobian of the Contacts) that defines the transformation from the motion space to the constraint space:

$$ \underbrace{\left(\frac{1}{h}\mathbf{M} + \mathbf{B} + h\mathbf{K}\right)}_{\mathbf{A}} \underbrace{\frac{d\dot{\mathbf{q}}}{\mathbf{x}}}_{\mathbf{b}} = \underbrace{-h\mathbf{K}\dot{\mathbf{q}}_a - \mathbf{f}_a - \mathbf{p}_b}_{\mathbf{b}} + \mathbf{H}^T \boldsymbol{\lambda}_b \quad (5) $$

### B. MJED method for Hyperelastic materials

Any hyperelastic material is fully determined by its strain energy function $W_h$ which describes the amount of energy necessary to deform the material. This strain energy function often involves the invariants of the right Cauchy-deformation tensor defined as $C = \nabla\Phi^T\nabla\Phi$, where $\Phi$ is the deformation function from the configuration at rest towards the deformed one. Most of the strain energy of hyperelastic isotropic materials is defined from 2 invariants of the deformation tensor $I_1 = \text{tr}C$ and $I_2 = ((\text{tr}C)^2 - \text{tr}C^2)/2$ and the Jacobian $J = \text{det}\nabla\Phi$.

The MJED method proposes to decouple the invariants of C from the expression of J to avoid assembling the

stiffness matrix and complex derivative expressions. The rationale behind the MJED method is that the derivative of J with respect to the nodal position is trivial for most of the finite element (especially for linear tetrahedra) whereas its derivative with respect to C is not trivial. Since the converse is true for the 2 invariants, the MJED method optimizes the force and the stiffness computation by separating both terms.

To accurately model the viscoelasticity of biological tissues, it is not sufficient to resort to Rayleigh damping. Instead, the MJED method relies on Prony series which consists in adding to the hyperelastic 2nd Piola Kirchoff stress tensor $S_h$ some time-dependent stresses (see [18] for details). Adding the viscous properties through the Prony series does not have a significant impact on the total computation times.

### C. Constraint-based simulation

According to (4), the dynamic behavior of two interacting objects is given by the following system:

$$ \mathbf{A}_1\mathbf{x}_1 - \mathbf{H}_1^T\boldsymbol{\lambda} = \mathbf{b}_1 \quad (6) $$

$$ \mathbf{A}_2\mathbf{x}_2 - \mathbf{H}_2^T\boldsymbol{\lambda} = \mathbf{b}_2 \quad (7) $$

$$ \mathbf{H}_1\mathbf{x}_1 + \mathbf{H}_2\mathbf{x}_2 \geq \boldsymbol{\delta} \quad (8) $$

$$ 0 \geq \boldsymbol{\delta}, \text{ and } \boldsymbol{\lambda} \geq 0, \text{ and } \boldsymbol{\delta}\boldsymbol{\lambda} = 0 \quad (9) $$

where suffix 1 and 2 denotes the matrices of respectively objects 1 and 2. Equation (9) enforces the Signorini law that describes a complementarity between contact force $\boldsymbol{\lambda}$ and penetration $\boldsymbol{\delta}$ (i.e. either objects are distant $\boldsymbol{\delta} > 0$ then the contact force vanish $\boldsymbol{\lambda} = 0$, or a positive contact force $\boldsymbol{\lambda} > 0$ must be applied to cancel the penetration $\boldsymbol{\delta} = 0$)

The previous system of equation is solved using the Schür complement method. Defining $\mathbf{x}_1^{\text{free}} = \mathbf{A}_1^{-1}\mathbf{b}_1$ and $\mathbf{x}_2^{\text{free}} = \mathbf{A}_2^{-1}\mathbf{b}_2$, and replacing (6) and (7) in (8) gives:

$$ \underbrace{(\mathbf{H}_1\mathbf{A}_1^{-1}\mathbf{H}_1^T + \mathbf{H}_2\mathbf{A}_2^{-1}\mathbf{H}_2^T)}_{\mathbf{W}} \boldsymbol{\lambda} = \underbrace{\boldsymbol{\delta} - \mathbf{H}_1\mathbf{x}_1^{\text{free}} - \mathbf{H}_2\mathbf{x}_2^{\text{free}}}_{\boldsymbol{\delta}^{\text{free}}} \quad (10) $$

where **W** is the compliance matrix that models the coupling between contact points in the constraint space. In this equation both $\boldsymbol{\lambda}$ and $\boldsymbol{\delta}$ are unknown.

Adding Coulomb friction the reaction force is enclosed into a cone whose height and direction are given by the normal force $\boldsymbol{\lambda}_n$.

$$ 
\begin{aligned}
\dot{\boldsymbol{\delta}}_{\text{T}} &= 0 \Rightarrow \|\boldsymbol{\lambda}_{\text{T}}\| < \mu \|\boldsymbol{\lambda}_{\text{n}}\| \quad (\text{stick}) \\
\dot{\boldsymbol{\delta}}_{\text{T}} &\neq 0 \Rightarrow \boldsymbol{\lambda}_{\text{T}} = -\mu \|\boldsymbol{\lambda}_{\text{n}}\| \frac{\dot{\boldsymbol{\delta}}_{\text{T}}}{\|\dot{\boldsymbol{\delta}}_{\text{T}}\|} = -\mu \|\boldsymbol{\lambda}_{\text{n}}\| \mathbf{T} \quad (\text{slip})
\end{aligned}
\quad (11) $$

where $\mu$ is the friction parameter and $\mathbf{T}$ is the tangential plane to the contact normal $\mathbf{n}$. If the reaction force is strictly included inside the cone, then objects stick together, but if the reaction force is on the cone's border then objects are slipping along the tangential direction. In this last case, the friction force must be directed along motion direction.

Equation (10) describes a NLCP (Non-linear complementarity problem) that is solved using a modified Gauss-Seidel algorithm (see [11] for details). Once $\boldsymbol{\lambda}$ is obtained it is replaced in equations (6) and (7) to compute a corrective
motion:

$$ \begin{aligned} \mathbf{x}_1 &= \mathbf{x}_1^{\text{free}} - \mathbf{A}_1^{-1} \mathbf{H}_1^T \lambda \\ \mathbf{x}_2 &= \mathbf{x}_2^{\text{free}} - \mathbf{A}_2^{-1} \mathbf{H}_2^T \lambda \end{aligned} \tag{12} $$

Using equations (2) and (5), $\mathbf{x}_1$ and $\mathbf{x}_2$ are finally integrated obtaining $\mathbf{q}_{1,a}$ and $\mathbf{q}_{2,b}$, the positions of objects 1 and 2 at the end of the time step that fulfil the laws of contact and friction.

## IV. ASYNCHRONOUS COMPUTATION OF CONSTRAINT

In order to compute the compliance matrix $\mathbf{W}$ in equation (10) the inversion of the large matrix $\mathbf{A}$ must be performed. Its value depends on the stiffness matrix $\mathbf{K}$ and its value may significantly vary due to large deformations of hyperelastic materials. This implies that the inversion of this matrix at haptic rates is not possible. We then introduce our method where the compliance matrix $\mathbf{W}$ is asynchronously evaluated at low rates, while a local contact problem is computed at high rates to update haptic forces.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Multithreaded approach for the haptic rendering of hyperelastic models. A preconditioner representing an approximation of the compliance matrix, is updated at low frequency. A local contact problem is updated at haptic rate. The contact directions and the compliance are assumed constant and only the violation $\delta$ is updated according to the motion of the device.</figcaption></figure>

### A. Approximation of the mechanical coupling

We extended the method proposed in [9] where a preconditioner $\mathbf{P} \simeq \mathbf{A}$ is updated at low frequency. This method relies on the assumption that $\mathbf{A}$ undergoes small perturbations between two consecutive time steps. Therefore if at a specific time $t$, $\mathbf{P}_t = \mathbf{A}_t^{-1}$ is calculated it may be considered a “good” approximation for the following time steps. A sparse factorization $\mathbf{P} = \mathbf{A}_{\text{old}} = \mathbf{L}\mathbf{D}\mathbf{L}^T$ is updated at low frequency on a dedicated CPU thread, and the last preconditioner available is used to let the simulation advance (see fig 1). Moreover, [9] has shown that using GPU parallelization technique, allows to evaluate the compliance matrix $\mathbf{W}$ in real time, by computing columns independently $\mathbf{L}^{-1}\mathbf{H}$, i.e. the product of the inverse of the lower triangular system $\mathbf{L}$ with the Jacobian of the contacts.

With hyperelastic materials, the assumption of small perturbation remains valid if the preconditioner is updated sufficiently fast. This is the case for meshes of reasonable size (i.e. with a number of nodes still compatible with real-time computations). In order to limit the divergence of the preconditioner, we estimate the nodal rotations $\mathbf{R}$ that were introduced from the last update of the preconditioner. Since no computation of the rotation is necessary in MJED

formulation, we use the shape matching method proposed in [21] to evaluate the average nodal rotation matrix. Therefore, the most recent preconditioner $\mathbf{P}_{\text{old}}$ is rotated using the current rotation matrix $\mathbf{R}_{\text{cur}}$ as follows:

$$ \mathbf{P} = \mathbf{R}_{\text{cur}}^{\text{T}} (\mathbf{L}_{\text{old}}\mathbf{D}_{\text{old}}\mathbf{L}_{\text{old}}^{\text{T}}) \mathbf{R}_{\text{cur}} \tag{13} $$

Lastly, this preconditioner is less sensitive to geometrical non-linearities and the resulting compliance matrix provides a close approximation of the exact mechanical coupling between contacts.

### B. Haptic contact force computation

In order to update haptic forces, our approach consists in sharing constraints equations and compliance between a haptic control loop (at high rates $\simeq 1$ kHz) and the simulation (at low rates $\simeq 30$ Hz). After performing the collision detection and the computation of deformable models, the NLCP is defined and solved in the simulation. The compliance matrix $\mathbf{W}$ and the Jacobian of the contacts $\mathbf{H}$ are shared with a separate haptic loop. In the haptic loop, the position of the device is refreshed to update the motion of the virtual instrument driven by the device motions. The displacement provides new violation in the constraint space, then a new value of the contact response is computed and sent to the force feedback device.

Contact constraints are set in the simulation before being actually processed and rendered in the haptic rendering. Therefore, the method relies on the computation of proximity queries. Consequently, if the virtual tool is approaching an obstacle in the simulation, the contact is set and it can be activated in the haptic loop (for instance if the user continues its motion towards the obstacle) before collision in the simulation.

## V. RESULTS

In order to validate the non-linear behavior of MJED models we simulated compression and elongation tests (see fig 2(b)). A cube is fixed onto one of its sides and a displacement is imposed on the opposite side as to extend and compress the material. In figure 3 the intensity of the force is measured depending on the nature of the hyperelastic material. For the co-rotational model the force varies as a linear function of displacement as any rotation is applied, whereas the force other materials is non-linear. As expected, during the extension, the St Venant Kirchhoff material becomes harder while Neo-Hookean and Mooney Rivlin become softer. A converse relationship is obtained when a compression test is performed. The accuracy of the method has been evaluated on the basis of the delay necessary to compute the preconditioner. We produced a simulation with a cube composed of 2058 linear tetrahedra and modeled with Mooney Rivlin material. We recorded the motion of an instrument controlled through the Phantom Omni haptic interface (see fig 2(c)) and we simulated frictionless contacts with the cube. In figure 4, we report the contact force $\lambda$ obtained with 3 different methods: i) the exact inverse of $\mathbf{A}$ computed at each time step ii) the asynchronous preconditioner proposed in this paper
<figure><img src="image_3.png" /><figcaption>Fig. 2. Real-time simulation of MJED hyperelastic models with constraints. 2(b) example of elongation. 2(c) frictionless contact with haptic force feedback. 2(d) grasping simulation using non-linear friction. 2(e) multicontact simulation with mechanical coupling. 2(f) Retina surgery simulation.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Intensity of the force (in Newton) versus displacement for the elongation test in elongation 2(a) and compression 2(b). Parameters have been chosen to have equivalent mechanical behaviours. Corotational: $E = 1$ and $\nu = 0.45$. St. Venant Kirchhoff: $\mu = 0.344$ and $\lambda = 3.10$. Neo Hookean: $\mu = 0.344$ and $k0 = 3.33$. Mooney Rivlin: $c1 = 0.086$ and $c2 = 0.086$ and $k0 = 3.33$.</figcaption></figure>

iii) a precomputed version of the inverse evaluated on the rest position. The precomputed version is equivalent to a contact response with linear materials which overestimates the force for large deformations introducing stability issues (the simulation crashed during the first contact, after two seconds of simulation). Results show a perfect match between contact forces of the exact solution and the asynchronous ones, even for large deformations, but our method is from 3× up to 5× faster than the exact solution. This allows the contact problem to be refreshed in the haptic loop much faster and to produce a realistic haptic rendering of the interactions. Table 5 reports the average computation

<figure><img src="image_5.png" /><figcaption>Fig. 4. Contact force (in Newton) when using a compliance matrix computed with: the exact inverse, a precomputed version from the rest configuration and our asynchronous preconditioner.</figcaption></figure>

time of the preconditioner for various meshes resolutions. Although the LDL<sup>T</sup> factorization becomes more expensive for detailed meshes, the delay to update the preconditioner remains limited, since the simulation is slowing down as well. Indeed, the computation of the factorization is 41× slower for a mesh composed of 2744 nodes than for a mesh composed of 512 nodes, but the delay necessary to update the preconditioner entails only 3.3 further simulation steps.

<table><caption>Fig. 5. Time to update the LDL<sup>T</sup> factorization and corresponding number of simulation step according to the number of nodes in the FE mesh.</caption><tr><td>Number of Nodes</td><td>512</td><td>1000</td><td>1728</td><td>2744</td></tr><tr><td>Time (sec)</td><td>0.03</td><td>0.11</td><td>0.41</td><td>1.23</td></tr><tr><td>Simulation step</td><td>3</td><td>5</td><td>7</td><td>10</td></tr></table>

Finally, we used our method in some more complex scenarios involving haptic feedback with multicontacts and friction. In fig. 2(d) the instrument is controlled through the Phantom Omni haptic interface, and we simulate multiple friction contacts allowing to grasp deformable objects. Such simulation is particularly difficult since the non-linearity of the friction (stick/slip transition) must be established according to the non-linear behaviour of the hyperelastic material. With our method the user can feel the effect of non-linearities thanks to the haptic interface and one can grasp the object in real time. We also performed a simulation involving 5 deformable organs (liver, diaphragm, intestine,
colon, and the stomach) in contact with each others (see fig. 2(e)). The liver is composed of 2000 tetrahedra elements and modeled with Mooney Rivlin material, while other organs are modeled with the co-rotational approach and computed on GPU. Through the haptic interface, the user can feel the effect of the mechanical coupling, for instance when the liver is in compression due to contacts with the diaphragm (that is also deformable), it provides a non-linear contact response force. The overall scene runs at 20 FPS. Finally we produced a simulation of retina surgery 2(f) where two instruments are inserted within the eye through trocars. Although there is no direct contact between them, the instruments are coupled through the non linear model of the eye and users can feel efforts of one instrument when using the other.

## VI. CONCLUSION

We proposed a method to provide a realistic haptic feedback of deformable models based on the MJED hyperelastic formulation. Our method relies on the computation of an asynchronous preconditioner able to provide a close approximation of the mechanical coupling of the contacts (that is subject to strong non linearites due to complex interactions and hyperelastic materials). Another asynchronous loop is used to update the local contact problem at haptic rate (assuming a linearization of contacts). In the future we plan to validate our method by comparing the force provided by the simulation with a force sensor on phantoms.

## REFERENCES

[1] Y. Adachi, T. Kumano, and K. Ogino. Intermediate representation for stiff virtual objects. In *Proceedings of the Virtual Reality Annual International Symposium (VRAIS'95)*, VRAIS '95, pages 203–Washington, DC, USA, 1995. IEEE Computer Society.

[2] Richard J. Adams, Blake Hannaford, Member Ieee, and Member Ieee. Control law design for haptic interfaces to virtual reality. *IEEE Transactions on Control Systems Technology*, 10:3–13, 2002.

[3] Richard J. Adams, Manuel R. Moreyra, and Blake Hannaford. Stability and performance of haptic displays: Theory and experiments. In *Proceedings ASME International Mechanical Engineering Congress and Exhibition*, pages 227–234, 1998.

[4] Jernej Barbič and Doug L. James. Six-dof haptic rendering of contact between geometrically complex reduced deformable models. *IEEE Trans. Haptics*, 1(1):39–52, 2008.

[5] C. Basdogan and A.A. Srinivasan. Haptic rendering in virtual environments. In *In Stanney, K.(Ed.)*, *Handbook of Virtual Environments, Lawrence Erlbaum, Inc*, 2002.

[6] Morten Bro-Nielsen and Stéphane Cotin. Real-time volumetric deformable models for surgery simulation using finite elements and condensation. *Comput. Graph. Forum*, 15:57–66, 1996.

[7] J.E. Colgate, M.C. Stanley, and J.M. Brown. Issues in the haptic display of tool use. *Intelligent Robots and Systems, IEEE/RSJ International Conference on*, 3:3140, 1995.

[8] Stéphane Cotin, Hervé Delingette, and Nicholas Ayache. Real-time elastic deformations of soft tissues for surgery simulation. *IEEE Transactions on Visualization and Computer Graphics*, 5(1):62–73, 1999.

[9] Hadrien Courtecuisse, Jérémie Allard, Pierre Kerfriden, Stéphane P a Bordas, Stéphane Cotin, and Christian Duriez. Real-time simulation of contact and cutting of heterogeneous soft-tissues. *Medical Image Analysis*, 18:394–410, 2014.

[10] C. Duriez, C. Andriot, and A. Kheddar. Interactive haptic for virtual prototyping of deformable objects: Snap-in tasks case. In *EUROHAPTICS*. Citeseer, 2003.

[11] Christian Duriez, F. Dubois, a. Kheddar, and C. Andriot. Realistic haptic rendering of interacting deformable objects in virtual environments. *IEEE Transactions on Visualization and Computer Graphics*, 12(1):36–47, 2006.

[12] C a Felippa and B Haugen. Unified Formulation of Small-Strain Corotational Finite Elements : I . Theory Unified Formulation of Small-Strain Corotational Finite Elements : I . Theory. (January), 2005.

[13] Clément Forest, Hervé Delingette, and Nicholas Ayache. Surface contact and reaction force models for laparoscopic simulation. In *International Symposium on Medical Simulation*, June 2004.

[14] Carlos Garre and Miguel A. Otaduy. Haptic rendering of complex deformations through handle-space force linearization. *World Haptics Conference*, 0:422–427, 2009.

[15] P. Jacobs, M.J. Fu, and M.C. Çavuşoğlu. High fidelity haptic rendering of frictional contact with deformable objects in virtual environments using multi-rate simulation. *The International Journal of Robotics Research*, 29(14):1778–1792, 2010.

[16] Doug L. James and Dinesh K. Pai. A unified treatment of elastostatic contact simulation for real time haptics. In *SIGGRAPH '05: ACM SIGGRAPH 2005 Courses*, page 141, 2005.

[17] Mohsen Mahvash and Vincent Hayward. High-fidelity haptic synthesis of contact with deformable bodies. *IEEE Comput. Graph. Appl.*, 24(2):48–55, 2004.

[18] Stéphanie Marchesseau, Tobias Heimann, Simon Chatelin, Rémy Willinger, and Hervé Delingette. Multiplicative Jacobian energy decomposition method for fast porous visco-hyperelastic soft tissue model. *Lecture Notes in Computer Science*, 6361 LNCS:235–242, 2010.

[19] William R. Mark, Scott C. Randolph, Mark Finch, James M. Van Verth, and Russell M. Taylor, II. Adding force feedback to graphics systems: issues and solutions. In *Conference on Computer graphics and interactive techniques, SIGGRAPH '96*, pages 447–452, New York, NY, USA, 1996. ACM.

[20] K Miller, G Joldes, D Lance, and A Wittek. Total Lagrangian explicit dynamics finite element algorithm for computing soft tissue deformation. *Communications in numerical methods in engineering*, 23(2):121–134, 2007.

[21] Matthias Müller, Bruno Heidelberger, Matthias Teschner, and Markus Gross. Meshless deformations based on shape matching. *ACM Transactions on Graphics*, 24:471, 2005.

[22] Michael Ortega, Stephane Redon, and Sabine Coquillart. A six degree-of-freedom god-object method for haptic display of rigid bodies with surface properties. *IEEE Transactions on Visualization and Computer Graphics*, 13(3):458–469, 2007.

[23] Igor Peterlik, Mert Sedef, Cagatay Basdogan, and Luděk Matyska. Real-time visio-haptic interaction with static soft tissue models having geometric and material nonlinearity. *Computers & Graphics*, 34(1):43–54, 2010.

[24] G. Picinbono and J.-C. Lombardo. Extrapolation: a solution for force feedback? In *International Scientific Workshop on Virtual Reality and Prototyping*, pages 117–125, Laval France, June 3-4 1999.

[25] D. C. Popescu and M. Compton. A model for efficient and accurate interaction with elastic objects in haptic virtual environments. In *GRAPHITE '03: Proceedings of the 1st international conference on Computer graphics and interactive techniques in Australasia and South East Asia*, pages 245–250, 2003.

[26] Diego C. Ruspini, Krasimir Kolarov, and Oussama Khatib. The haptic display of complex graphical environments. In *SIGGRAPH '97: Proceedings of the 24th annual conference on Computer graphics and interactive techniques*, pages 345–352, New York, NY, USA, 1997. ACM Press/Addison-Wesley Publishing Co.

[27] Guillaume Saupin, Christian Duriez, and Stephane Cotin. Contact model for haptic medical simulations. In *ISBMS '08: Proceedings of the 4th international symposium on Biomedical Simulation*, pages 157–165, 2008.

[28] Mert Sedef, Evren Samur, and Cagatay Basdogan. Real-time finite-element simulation of linear viscoelastic tissue behavior based on experimental data. *IEEE Comput. Graph. Appl.*, 26:58–68, 2006.

[29] Y. Zhuang and J. Canny. Real-time simulation of physically realistic global deformation. In "Hot-topic", in *IEEE Visualization Conference*, San Francisco, 1999.

[30] C. B. Zilles and J. K. Salisbury. A constraint-based god-object method for haptic display. In *IROS '95: Proceedings of the International Conference on Intelligent Robots and Systems-Volume 3*, page 3146, Washington, DC, USA, 1995. IEEE Computer Society.