Michael Döhler
email: michael.doehler@inria.fr

Laurent Mevel
email: laurent.mevel@inria.fr

Fault Isolation and Quantification from Gaussian Residuals with Application to Structural Damage Quantification

Keywords: Fault isolation, residual evaluation, statistical tests, damage quantification, structural health monitoring

Fault detection for structural health monitoring has been a topic of much research during the last decade. Localization and quantification of damages, which are linked to fault isolation, have proven to be more challenging, and at the same time of higher practical impact. While damage detection can be essentially handled as a data-driven approach, localization and quantification require a strong connection between data analysis and physical models. This paper builds upon a hypothesis test that checks if the mean of a Gaussian residual vector -whose parameterization is linked to possible damage locations -has become non-zero in the faulty state. It is shown how the damage location and extent can be inferred and robust numerical schemes for their estimation are derived based on QR decompositions and minmax approaches. Finally, the relevance of the approach is assessed in numerical simulations of two structures.

INTRODUCTION

Damage in structures impairs their operation if not detected early enough, in which case the operational state can be improved and the remaining life of the structure lengthened. For this reason, the question of vibrationbased damage diagnosis, i.e. damage detection, localization and quantification, has been studied from many different angles in the last decade [START_REF] Fan | Vibration-based damage identification methods: a review and comparative study[END_REF].

For example, purely data-driven empirical approaches based e.g. on mode shapes curvatures or the wavelet transform have been considered for certain classes of structures such as beams. On the other side, model-based approaches rely on an analytical structural model, on the premise that damage causes changes in the parameters of such a finite element model (FEM), involving the mass, damping and stiffness matrices of the structure. Based on a well-fitted reference model, the problem of damage localization and quantification may be resolved through sensitivity analysis and model updating with measurements from the faulty state [Teughels et al., 2002, Viet Hà and[START_REF] Viet Hà | Localization and quantification of damage in beam-like structures using sensitivities of principal component analysis results[END_REF]. However, FEM updating is often an ill-posed problem, since usually the dimension of the FEM is much larger than the number of identified parameters from measurements.

Alternative methods for damage localization with a theoretical background avoid the model updating step, while linking FEM information to data from both safe and damaged structures in statistical tests. For example, changes in the structural flexibility are interrogated [START_REF] Marin | Uncertainty quantification for stochastic damage localization for mechanical systems[END_REF] or the FEM parameters of a Gaussian residual vector are tested for a change based on the local approach [Ben-veniste et al., 1987[START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF][START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios[END_REF][START_REF] Döhler | Subspace-based fault detection robust to changes in the noise covariances[END_REF][START_REF] Döhler | Efficient computation of minmax tests for fault isolation and their application to structural damage localization[END_REF], on which this paper is based. While many fault detection and isolation (FDI) methods concern additive faults, e.g. [Dong andVerhaegen, 2009, Dong et al., 2012], our approach assumes a general description of the fault affecting the state space matrices.

In this work, FDI is based on a Gaussian residual vector with zero mean in the reference state and non-zero mean in the faulty state. The fault isolation problem consists in deciding which components in a related parameterization have changed by hypothesis testing [START_REF] Basseville | Information criteria for residual generation and fault detection and isolation[END_REF]Nikiforov, 1993, Basseville, 1997]. For FEM parameterizations, this is related to damage localization in [START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios[END_REF][START_REF] Balmès | Statistical model-based damage localization: a combined subspace-based and substructuring approach[END_REF]. It was demonstrated in [START_REF] Döhler | Efficient computation of minmax tests for fault isolation and their application to structural damage localization[END_REF] that the minmax test has more appropriate properties and is better suited for the damage localization problem when the parameters are not fully independent and multiple parameters are changing. Minmax approaches have also been applied to damping monitoring [START_REF] Zouari | Subspace-based damping monitoring[END_REF] or damage detection with temperature rejection [START_REF] Balmès | Handling the temperature effect in vibration-based monitoring of civil structures: a combined subspace-based and nuisance rejection approach[END_REF]. From a numerical point of view, QR approaches have been proved to be a tool to enforce stability for such algorithms [START_REF] Zhang | Advanced numerical computation of χ 2 -tests for fault detection and isolation[END_REF]Basseville, 2003, Döhler et al., 2014b]. While this framework allows fault isolation for Gaussian residual vectors, it has never been attempted to quantify the absolute change in the faulty parameter components in this context. In this paper, we extend the local approach to damage quantification. This paper is organized as follows. In Sections 2 and 3, structural modelization and numerical schemes for hypothesis testing are recalled. In Section 4, fault quantification is addressed and applied to structural damage diagnosis in Section 5. Numerical applications of the proposed scheme are reported in Section 6.

MOTIVATION: STRUCTURAL DAMAGE DIAGNOSIS BASED ON LOCAL APPROACH

The behavior of linear time-invariant dynamical structures subject to unknown ambient excitation can be described by system

M Ẍ (t) + C Ẋ (t) + KX (t) = υ(t)
(1) where t denotes continuous time; M, C, K ∈ R m×m are mass, damping, and stiffness matrices, respectively; the (high dimensional) state vector X (t) ∈ R m is the displacement vector of the m degrees of freedom of the structure; and υ(t) is the external unmeasured force (noise).

Observing system (1) at r sensor positions (e.g. displacement, velocity or acceleration sensors) at discrete time instants t = kτ (with sampling rate 1/τ ), it can be transformed to an equivalent discrete-time state space system [START_REF] Juang | Applied system identification[END_REF] x

k+1 = Ax k + v k y k = Cx k + w k (2)
with the states x k ∈ R n , the measured outputs y k ∈ R r , the state transition matrix

A = exp 0 I -M -1 K -M -1 C τ ∈ R n×n
and the observation matrix

C = L d -L a M -1 K L v -L a M -1 C ∈ R r×n , where n = 2m is the model orderand L d , L v , L c ∈ {0, 1}
r×m are selection matrices indicating the positions of displacement, velocity or acceleration sensors. The state noise v k and output noise w k are unmeasured and assumed to be Gaussian, zero-mean, white.

Our system parameter θ of interest is linked to the structural properties of system (1), and can be e.g. mass parameters of the elements in a finite element model (FEM), i.e. M = M(θ), or parameters corresponding to element stiffness such as Young's modulus of the individual structural elements, i.e. K = K(θ). In particular, A = A(θ) in the corresponding system (2). Structural damage leads to changes in θ, and hence the considered faults in system (2) are changes in θ. They are not additive. It is assumed that structural damage is local and thus affects only a small part of the parameters in θ. After damage has been localized, i.e. it has been decided which parameters in θ are faulty ("fault isolation") based on measurements {y k } k=1,...,N of the faulty system, we propose a strategy to determine the change δ in the faulty parameters in order to perform damage quantification.

To achieve this task, our work is based on the asymptotic local approach for change detection [START_REF] Benveniste | The asymptotic local approach to change detection and model validation[END_REF]. It assumes the close hypotheses H 0 : θ = θ 0 (reference system),

H 1 : θ = θ 0 + δ/ √ N (faulty system), (3)
where vector θ 0 denotes the system parameter in the reference state and vector δ is unknown but fixed. With this statistical framework, very small changes in the system parameter θ can be detected if the number of measurements N is large enough. For FDI, a subspace-based residual vector ζ N has been proposed in [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF] and is detailed in Section 5. Using the structural parameterization θ as explained above, it yields the central limit theorem (CLT)

ζ N d -→ N (0, Σ) under H 0 N (J δ, Σ) under H 1 (4)
for N → ∞, where J and Σ are the asymptotic sensitivity and covariance, respectively. In this framework, we recall suitable statistical tests for fault isolation in the following section.

BASIC FAULT ISOLATION TESTS

Definitions

Let a parameter vector θ ∈ R l and a Gaussian residual vector ζ = ζ(θ) ∈ R h be given, yielding (4). Its sensitivity matrix J ∈ R h×l is assumed to have full column rank and its covariance Σ ∈ R h×h is positive definite. Many methods have been proposed in the literature for fault detection (decide if there is a change in θ) and isolation (decide which elements of vector θ have changed) by the means of such a Gaussian residual vector ζ [START_REF] Basseville | Information criteria for residual generation and fault detection and isolation[END_REF]Nikiforov, 1993, Basseville, 1997].

For both problems of fault isolation and quantification, different partitions of the vector δ into two subvectors are considered. Without loss of generality, let this partition be

δ = δ a δ b .
For fault isolation, a decision between δ a = 0 and δ a = 0 is made for each partition. Let F = J T Σ -1 J be the Fisher information matrix of the parameter θ contained in vector ζ, and let J and F be partitioned accordingly as

J = [J a J b ], F = F aa F ab F ba F bb = J T a Σ -1 J a J T a Σ -1 J b J T b Σ -1 J a J T b Σ -1 J b .
In the next sections, two isolation tests are recalled and extended to the quantification of the faults, i.e. the estimation of δ a in Section 4.

Sensitivity Test

The simplest possibility for testing δ a = 0 (no change in parameter subset 'a') against δ a = 0 is to assume δ b = 0 and thus ζ ∼ N (J a δ a , Σ). The corresponding Generalized Likelihood Ratio (GLR) test writes as [START_REF] Basseville | Information criteria for residual generation and fault detection and isolation[END_REF]

t sens = ζ T Σ -1 J a J T a Σ -1 J a -1 J T a Σ -1 ζ, (5)
which is called sensitivity test. The test variable t sens is χ 2 -distributed with dim(θ a ) degrees of freedom and noncentrality parameter δ T a F aa δ a , if δ b = 0 is actually true. For a decision, the test variable is compared to a threshold. Remark 1. If the assumption δ b = 0 does not hold, the non-centrality parameter of t sens follows as

δ T a F aa δ a + 2δ T a F ab δ b + δ T b F ba F -1 aa F ab δ b . (6)

Minmax Test

Instead of assuming δ b = 0, the variable δ b is replaced by its least favorable value for a decision about δ a , as follows.

Define the partial residuals as

ζ a def = J T a Σ -1 ζ, (7a) ζ b def = J T b Σ -1 ζ, (7b
) the robust residual as [START_REF] Basseville | Information criteria for residual generation and fault detection and isolation[END_REF]

ζ * a def = ζ a -F ab F -1 bb ζ b (8) and F * a def = F aa -F ab F -1 bb F ba .
(9) Then, the mean of the robust residual ζ * a is sensitive to changes δ a , but blind to δ b , and it holds

ζ * a ∼ N (F * a δ a , F * a ) .
(10) The corresponding GLR test for δ a = 0 (no change in parameter subset 'a') against δ a = 0 writes as

t mm = ζ * T a F * -1 a ζ * a , (11)
which is called minmax test. The test variable t mm is χ 2 -distributed with dim(θ a ) degrees of freedom and noncentrality parameter δ T a F * a δ a , independently of δ b . For a decision, the test variable is compared to a threshold. Note that the invertibility of all matrices in the computation is guaranteed, since J is assumed to have full column rank and Σ is positive definite.

Numerically Efficient Computation

Both the sensitivity and the minmax test require a number of matrix inversions, which may be numerically critical due to possible ill-conditioning of the covariance matrix Σ. In [START_REF] Zhang | Advanced numerical computation of χ 2 -tests for fault detection and isolation[END_REF]Basseville, 2003, Döhler et al., 2014b] numerically more robust computations of both tests were suggested, which are recalled in the following. They make use of the decomposition Σ -1 = (Σ -1/2 ) T Σ -1/2 , (12) which can be obtained e.g. from an SVD of Σ = U ∆U T as Σ -1/2 = ∆ -1/2 U T , as well as QR decompositions. Note that an efficient computation of Σ -1/2 is detailed in [Döhler andMevel, 2011, Döhler et al., 2014a].

Sensitivity Test

Using the thin QR decomposition [START_REF] Golub | Matrix computations[END_REF] Σ -1/2 J a = QR, (13) the sensitivity test in (5) writes as [START_REF] Zhang | Advanced numerical computation of χ 2 -tests for fault detection and isolation[END_REF]]

t sens = α T α, where α = Q T Σ -1/2 ζ.
(14) Here, the number of numerically critical matrix inversions and multiplications is limited to a minimum compared to a direct computation in (5).

Minmax Test Let the thin QR decomposition of

Σ -1/2 [J b J a ] = [Q b Q a ] R bb R ba 0 R aa (15)
be given and partitioned accordingly. Then, the minmax test in (11) writes as [START_REF] Döhler | Efficient computation of minmax tests for fault isolation and their application to structural damage localization[END_REF]]

t mm = β T β, where β = Q T a Σ -1/2 ζ. (16)

FAULT QUANTIFICATION

In this section, estimates of the fault δ a are derived based on the properties of the sensitivity and minmax tests, respectively. Special care is also taken of a numerically sensible computation of these estimates.

Sensitivity Approach

Theorem 2. Define

δ sens a def = (J T a Σ -1 J a ) -1 J T a Σ -1 ζ. (17) Then, under the assumption δ b = 0, δ sens a ∼ N δ a , F -1 aa .
Proof. Since ζ ∼ N (J δ, Σ), it follows under the assumption

δ b = 0 (see Section 3.2) that ζ ∼ N (J a δ a , Σ) . It follows easily J T a Σ -1 ζ ∼ N (J T a Σ -1 J a δ a , J T a Σ -1 J a )
and, since F aa = J T a Σ -1 J a , the assertion follows.

Since the sensitivity approach requires the assumption δ b = 0, which is in reality never guaranteed, the potential error can be analyzed in the following corollary. Corollary 3. For δ b arbitrary, it holds

δ sens a ∼ N δ a + F -1 aa F ab δ b , F -1 aa . ( 18
)
Proof. Since ζ ∼ N (J a δ a + J b δ b , Σ), the assertion follows from the definition of δ sens a in (17).

A numerically sensible computation of δ sens a can be achieved as follows. Corollary 4. With the QR decomposition (13), namely Σ -1/2 J a = QR, it holds

δ sens a = R -1 Q T Σ -1/2 ζ, coinciding with R -1 α in the test computation in (14).
Proof. With decomposition (12), δ sens a in (17) can be written as

δ sens a = (Σ -1/2 J a ) T Σ -1/2 J a -1 (Σ -1/2 J a ) T Σ -1/2 ζ.
Plugging in the QR decomposition Σ -1/2 J a = QR, where

Q T Q = I and R is invertible, yields δ sens a = R T Q T QR -1 R T Q T Σ -1/2 ζ,
leading to the assertion.

Minmax Approach

Theorem 5. Let the variables for the minmax test be given in (7a)-( 9). Define

δ mm a def = (F * a ) -1 ζ * a . (19) Then, δ mm a ∼ N (δ a , (F * a ) -1 ).
Proof. The assertion follows immediately from property (10).

Similar to the minmax test (11), the computation of δ mm a in ( 19) is a numerical challenge since the computation of ζ * a and F * a requires a number of matrix operations (inversions, subtractions) that may be critical especially if Σ is badly conditioned. In ( 15)-( 16) a solution for the minmax test computation was recalled from [START_REF] Döhler | Efficient computation of minmax tests for fault isolation and their application to structural damage localization[END_REF] that is based on the QR decomposition, requiring a minimum of the former numerically critical operations. Based on these results, a simple computation for δ mm a is based on the same QR decomposition as follows. Corollary 6. With the QR decomposition (15), it holds

δ mm a = R -1 aa Q T a Σ -1/2 ζ, coinciding with R -1
aa β in the test computation in (16).

Proof. In the proof of [Döhler et al., 2014b, Thm. 1] it was shown that QR decomposition (15) yields

ζ * a = R T aa Q T a Σ -1/2 ζ, F * a = R T aa R aa .
Substituting these results in (19) leads immediately to the assertion.

STRUCTURAL DAMAGE LOCALIZATION AND QUANTIFICATION

We now return to our application on structural damage localization and quantification, as laid out in Section 2. Based on the asymptotic local approach for change detection [START_REF] Benveniste | The asymptotic local approach to change detection and model validation[END_REF], we assume the close hypotheses (3), where θ is a parameter vector containing local structural parameters, as for example the stiffness corresponding to each structural element.

In this framework, a subspace-based residual vector has been defined in [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF] based on outputs {y k } k=1,...,N of system (2) as

ζ N = √ N vec(S T H p+1,q ), (20)
where H p+1,q is an estimate of the Hankel matrix

H p+1,q def =     R 1 R 2 . . . R q R 2 R 3 . . . R q+1 . . . . . . . . . . . . R p+1 R p+2 . . . R p+q    
containing the output correlations R i = E(y k y T k-i ), and S is the left null space of H p+1,q from the reference system.

The residual vector ζ N in (20) satisfies the Central Limit Theorem (4) [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF]. The estimation of its asymptotic sensitivity J is described in detail in [START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios[END_REF][START_REF] Balmès | Statistical model-based damage localization: a combined subspace-based and substructuring approach[END_REF] for a chosen structural parameterization θ, where the residual is first derived wrt. the modal parameters (frequencies, mode shapes) obtained from the eigenvalues and eigenvectors of system (2), which are then derived wrt. the structural parameters using FEM (1). The estimation of Σ is described in [Döhler andMevel, 2011, Döhler et al., 2014a].

Based on these considerations, both the sensitivity and the minmax tests from Section 3 can be applied to residual ζ N (asymptotically, N large enough) for each of the structural parameters in θ = θ 1 . . . θ l to perform damage localization. A parameter whose test value exceeds a threshold is regarded as faulty and corresponds to a damage location in the structure. Finally, the quantification of the damage extent can be performed with the respective estimators derived in Section 4. Note that estimates δ are obtained. The parameter change that we are interested in for damage quantification, however, is the absolute change θ -θ 0 of the structural parameters, which can be obtained from θ -θ 0 = δ/ √ N due to the local hypothesis in (4).

APPLICATION

The damage localization and quantification techniques have been applied to two simulated structures. The first example is a simple mass-spring chain, the second a truss structure. In both cases, output-only datasets with displacement samples were generated at the sensor coordinates in reference and damaged states from white noise excitation. Measurement noise was added with a magnitude of 5% of each generated output signal. The damping was always defined such that all modes have a damping ratio of 2%. All parameters of the tests (S, J , Σ) were estimated based on the dataset in the reference state and the information from the respective FEM.

In the damaged state, both the sensitivity and the minmax tests from Section 3 were applied to the datasets for each structural element, and the damage extent was estimated for the damaged element(s) based on the respective tests as shown in Sections 4-5. In each quantification example, 100 realizations of the simulated time series were used. The quantification results are shown for several damage cases and extents, each as the mean from the 100 estimates together with their standard deviation.

Mass-spring chain

A mass-spring chain with eight elements is considered (Fig. 1) with masses

m 1 = m 3 = m 5 = m 7 = 1, m 2 = m 4 = m 6 = m 8 = 2 and stiffnesses k 1 = k 3 = k 5 = k 7 = 200, k 2 = k 4 = k 6 = k 8 = 100.
Output time series of length N = 100,000 are simulated at time step τ = 0.05 s at four elements. J is computed on all eight modes of the system. We consider three damage cases, where damage is simulated by a stiffness reduction in one or two springs.

Damage localization

Damage in element 4, Fig. 2 While the minmax tests behave perfectly, the sensitivity test for element 3 also reacts due to the violation of δ b = 0. The strong reaction only for element 3 among the undamaged elements becomes clear when evaluating the non-centrality parameter of the test in (6), in particular the value of F ba F -1 aa F ab where 'a' corresponds to any undamaged element and 'b' corresponds to the damaged element 4: this value is at 22.3 for element 3 while being lower than 0.4 for all other undamaged elements.

Damage in elements 2 and 4, Fig. 3 The minmax tests behave perfectly again, while the damage in spring 2 cannot be detected with the sensitivity tests as the test value of element 3 is significant. This is analogous to the previous case, where it was explained why the test for element 3 reacts when element 4 is damaged (see Fig. 2). The minmax test is robust to changes in the non-tested parameters θ b by design.

Damage in elements 3 and 4, Fig. 4 The sensitivity and minmax tests behave very well in this case. Note that the test values for element 3 is in both tests larger than for element 4, while element 4 is more damaged. Indeed, the test values serve only for the decision if the respective element is damaged or not by comparing it to a threshold. They are not directly linked to the damage extent, which is estimated separately in the following as shown in Section 4. isfying, where again only large damage extents are slightly overestimated. The damage quantification from the sensitivity approach is biased (here underestimated), since the assumption δ b = 0 is violated. As shown in Corollary 3, the theoretic bias for δ a is F -1 aa F ab δ b in (18), and F aa > 0 by definition. Indeed, F ab = -0.4 < 0 between the parameters corresponding to damaged elements 2 and 4 in this case, which explains the underestimation.

Damage quantification

Damage in elements 3 and 4, Fig. 5 (right) The results for the minmax approach are similar as in the previous two-damage case, whereas the damage extent is strongly overestimated now when using the sensitivity approach. This is in line with F ab = 2.09 > 0 for the parameters corresponding to elements 3 and 4 in this case.

Truss

A truss structure with 25 elements of equal stiffness properties, having six sensors, has been considered as a more complex example (Fig. 6). Output time series of length N = 50,000 are simulated at time step τ = 0.05 s. At this sampling frequency, we are in the more realistic case of modal truncation, since only ten modes are present in the data that can be taken into account for the computation of J , compared to altogether 25 modes of the analytical model. In general, only a small number of modes can be obtained from measurements of system (2) compared to the number of modes present in the FEM (1). Furthermore, S is estimated on the simulated data in this example and not on the analytical modes of the model. Damage is simulated by decreasing the element stiffness. Damage in element 16, Fig. 7 (left) The estimates from both the sensitivity and the minmax approach are close to the true values. The estimation errors are however larger for large damage extents, maybe due to the nature of the local approach: the relation between ζ N and δ is non-linear and the first-order approximation involving J (computed at δ = 0) may thus be less accurate for large changes δ.

Damage in elements 16 and 23, Fig. 7 (right) The minmax approach yields quite accurate estimates, while the sensitivity approach overestimates the damage extent significantly, similar to the example in Fig. 5 (right). This overestimation can be explained analogously by F ab = 0.2 > 0 between damaged elements 16 and 23. While this example does not satisfy the theoretical assumptions perfectly anymore (modal truncation), it still gives reasonable results. This shows that the developed damage quantification approach is promising, and further investigation is required for application on real structures.

CONCLUSIONS

In this paper, the quantification of parameter changes has been investigated based on a Gaussian residual vector and hypothesis testing. It has been shown that the local hypothesis modeling for a parameter change in θ as δ/ √ N is not only a convenient mathematical tool for establishing the CLT for the Gaussian residual, but also has a sensible and meaningful interpretation leading to an effective quantification of the parameter change. Based on sensitivity and minmax tests, estimators for the parameter change have been derived together with QR based robust numerical schemes for their computation. An application was shown for vibration-based structural damage localization and quantification in a simulation study, where changes in local parameters of a structure indicate damage. It has been noticed that the local approach is more precise if parameter changes are small, as expected. With more than one parameter changing, i.e. more than one damage, the minmax test proved to be more effective than the sensitivity test for both damage localization and quantification, where it avoids false alarms at undamaged elements for localization and it avoids bias in the damage quantification. Further work should link damage quantification and compressed sensing to handle FEM dimensions that largely exceed the modal parameter dimension, which is related to the number of sensors and identified modes.

Fig. 1 .

 1 Fig. 1. Mass-spring chain with four sensors.

Fig. 3 .Fig. 4 .

 34 Fig. 2. Sensitivity tests (left) and minmax tests (right) for 10% damage in element 4.

Fig. 5 .

 5 Fig. 5. Quantification of different damage extents in elements 4 only (top), 2 and 4 (left), 3 and 4 (right).

Fig. 7 .

 7 Fig. 7. Quantification of different damage extents in element 16 (left) and both elements 16 and 23 (right).

This work was partially supported by the European project FP7-PEOPLE-2009-IAPP 251515 ISMS.


Karolina Golec

Matthieu Coquet

Florence Zara

Guillaume Damiand

Improvement of a Topological-Physical Model to manage different physical simulations

Keywords: Physical Simulation, Mass-Spring System, Tensor-Mass Model, Combinatorial Maps

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

In medical applications, the actual challenge is to propose physical models which provide interactive simulations with topological modifications (such as cutting). The goal is to simulate surgical operations for training. In computer graphics, the field of soft bodies simulation covers many methods [NMK + 06]. The two main approaches, namely Finite Element Method (FEM) and Mass-Spring System (MSS), focus either on accuracy or on the performance of the system with visually satisfactory results. MSS is faster and easier to implement than FEM, but provides less accurate results. However, some papers focus to improve the MSS model by integrating some mechanical properties (such as Young's modulus and Poisson's ratio) in the computation of the springs properties [BSSH04, LSH07, BBJ + 09].

Furthermore, as an alternative to the classic FEM, the Tensor-Mass (TM) approach has been introduced by Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee.

Delingette, Cotin and Picinbono [START_REF] Delingette | A Hybrid Elastic Model Allowing Real-Time Cutting Deformations and Force Feedback for Surgery Training and Simulation[END_REF] and extended by Schwartz [SDR + 05]. This approach is based on another way to solve the mechanical equations of the objects, with a direct computation of the forces applied on each node considering the evaluation of the strain energy density. This approach allows to compute elements' information regardless of the neighborhood of this element.

In this context, this paper proposes an extension of the work of Flechon [START_REF] Flechon | A generic topological framework for physical simulation[END_REF][START_REF] Flechon | A unified topological-physical model for adaptive refinement[END_REF] who proposed a unified topological-physical model (called LCC+MSS) adapted for a mass-spring simulation that handles cutting or piercing simulated objects. Only one model is used to define the topology, the geometry and the physical properties of the simulated object, enabling topological modifications during its simulation by updating on the fly. In this paper, our contributions are:

• The improvement of the efficiency of the topological-physical model for cutting by changing the architecture of the model.

• The improvement of the topological-physical model to manage several physical models. We illustrate this issue by integrating the mass-spring system and tensor-mass model which both naturally allow dynamical topological modifications thanks to a local formulation of the strain.

• The use of the topological-physical model to simulate realistic and accurate experiments on soft tissues.

The paper is organized as follows. Section 2 presents the two physical models used. Section 3 presents the unified topological-physical model proposed by Flechon. Section 4 presents our extension of this unified model with discussion in Section 4.3. Section 5 presents some experiments showing that the performances of the LCC+MSS model have been improved by the new architecture of the topological-physical model. Moreover, Section 6 shows that our model can be used to conduce real experiments made on soft tissues and to compare physical models.

PHYSICAL MODELS

Mass-Spring System. In computer graphics, massspring systems have largely been used to simulate the behavior of deformable objects [START_REF] Terzopoulos | Heating and melting deformable models[END_REF][START_REF] Bridson | Robust treatment of collisions, contact and friction for cloth animation[END_REF]. This physical model is based on the discretization of the object into a set of masses (also called particles) interconnected by springs. Then at each time step, the following main steps are performed for each particle:

(1) Computation of the forces applied on the particle due to springs and external forces. Remember that the force involved by a spring connecting particles i and j, with stiffness constant k i j and initial length l i j , is defined at time t by F i j (t) = F e i j (t) + F v i j (t).

• F e i j (t) is the elasticity force of this spring with F e i j (t) = k i j P j (t) -P i (t)l i j u i j (t), F e ji (t) = -F e i j (t),

where P i (t), P j (t) are the positions of i, j and u i j (t) = (P j (t) -P i (t))/ P j (t) -P i (t) .

• F v i j (t) is the viscosity force of this spring with

F v i j (t) = γ i j (V j (t) -V i (t)) • u i j (t) u i j (t),
where V i (t),V j (t) are the velocities of i, j and γ i j = 2 k i j (m i + m j )/2 the spring's viscosity coefficient where m i , m j are the masses of i, j.

(2) Computation of the acceleration according to Newton's second law with d 2 dt 2 P i (t) = F i (t)/m i where F i (t) represents the forces applied on particle i.

(3) Computation of the velocity and position according to the acceleration using a numerical integration scheme (for instance the semi-implicit Euler one).

Tensor-Mass model. The TM approach is based on the discretization of the object into several elements as for the FEM, but then the equations are solved locally, making this approach more suitable for topological modifications of the object. To account for various mechanical behaviours, several formulations have been presented: linear Hooke's model [START_REF] Cotin | A Hybrid Elastic Model allowing Real-Time Cutting, Deformations and Force-Feedback for Surgery Training and Simulation[END_REF], non-linear geometrical model based on Saint Venant-Kirchhoff's elasticity model [START_REF] Picinbono | Real-Time Large Displacement Elasticity for Surgery Simulation: Nonlinear Tensor-Mass Model[END_REF], anisotropic material [START_REF] Picinbono | Non-linear anisotropic elasticity for real-time surgery simulation[END_REF] and non-linear visco-elastic deformations [SDR + 05] with some pre-computations to accelerate the process.

The simulation's loop of the MT approach involves the following main steps for each element E of the object:

(1) Computation of the displacement of a point X inside E with

U E (X) n-1 ∑ j=0 Λ j (X) U j ,
where n is the number of 3D nodes P j defining the element, U j the displacement of each node P j from its initial position, and Λ j (X) some interpolation functions defined according to the type of element used for the discretization.

(2) Computation of the deformation energy according to the displacement of the element's nodes. In our work, we consider a non-linear elasticity behavior using the Green-Saint Venant strain-tensor

ε(X) = 1 2 ∇U T (X) + ∇U(X) + ∇U T (X) ∇U(X) .
The associated strain energy density is defined by

W (X) = λ 2 (tr ε(X)) 2 + µ tr ε(X) 2 ,
where λ and µ are the Lamé coefficients characterizing material stiffness.

(3) Computation of the elasticity force on any node

P j ∈ E for j ∈ [0, n -1] with F E (P j ) = - ∂W E (P j ) ∂U j ,
where W E (P j ) is the energy density of deformation of the considered element evaluated at node P j .

(4) Computation of the acceleration according to Newton's second law (in the same way that for MSS).

(5) Computation of the deformation and displacement of the object using a numerical integration scheme (in the same way that for MSS).

In this paper, we used the approach proposed by Faure [START_REF] Faure | An Implicit Tensor-Mass solver on the GPU for soft bodies simulation[END_REF]. The formulation of the forces applied on each node of the elements was generated thanks to symbolic computation. Thus, the steps (1), (2) and

(3) of the simulation loop were replaced by the direct formulation of the force applied on the nodes.

3D Linear Cell Complex. The unified topologicalphysical model (called LCC+MSS) proposed by Flechon [START_REF] Flechon | A generic topological framework for physical simulation[END_REF][START_REF] Flechon | A unified topological-physical model for adaptive refinement[END_REF] used as underlying data structure the 3D Linear Cell Complex (LCC) [START_REF] Damiand | CGAL User and Reference Manual[END_REF] from the CGAL Open Source geometric algorithms library [START_REF]CGAL User and Reference Manual[END_REF]. With this topological structure, an orientable object is represented using 3D combinatorial maps (called 3-maps) [START_REF] Lienhardt | N-dimensional generalized combinatorial maps and cellular quasimanifolds[END_REF][START_REF] Damiand | Combinatorial Maps: Efficient Data Structures for Computer Graphics and Image Processing[END_REF] as a subdivision of 0-cells (vertices), 1-cells (edges), 2-cells (faces) and 3-cells (volumes). The cells are described by darts, a generalization of half-edges [START_REF] Mäntylä | An Introduction to Solid Modeling[END_REF] to higher dimension, plus pointers between these darts called β i with i ∈ {0, 1, 2, 3}. The 3D linear cell complex data structure is then obtained by associating to 0-cells the coordinates in R 3 of corresponding points. This representation using darts and pointers describes the cells of an object as well as all the incidence and adjacency relationships. If two i-cells share a common (i -1)-cell, they are adjacent (for example two volumes sharing a common face). If one cell belongs to the boundary of another cell, these two cells are incident (for example a vertex belonging to the boundary of a face). This information is very useful to allow topological modifications during the simulation. Lastly, any information can be associated with cells thanks to the attribute notion (e.g. a length to edges or a color to faces). We denote by i-attribute(d) the attribute associated with the i-cell containing dart d.

3D LCC for MSS. Flechon added information to a 3D LCC to construct the LCC+MSS model which is suitable for a MSS based on the Baudet formulation for hexahedral mesh [BBJ + 09]. This formulation is available for basic hexahedra including four inner springs.

In the LCC+MSS model, a data structure Particle is associated to each 0-cell thanks to a 0-attribute and a data structure Spring is associated to each 1-cell thanks to a 1-attribute. Four additional springs inside each hexahedron are also described by the Spring data structure. But contrary to particles and springs previously described, these inner springs are not associated with any cells of the 3-map, but directly associated to the object.

• The Particle structure stores information related to a particle: its mass, velocity, acceleration and the sum of the forces applied to it. Moreover, it stores the list of the inner springs attached to it. • The Spring structure stores the physical properties of a spring: its initial length, its stiffness and its two extremities (thanks to two pointers to the two Particle connected by the considered spring). Interest of using topological model for physical simulation. The classical data-structure used in MSS simulations is a graph where vertices correspond to particles and edges correspond to springs. This simple model allows to directly and easily compute the forces resulting from the springs and to accumulate these forces on the particles. However, the treatment becomes more difficult for topological modifications. Indeed, it is not possible with such a graph, to directly retrieve the elements of the object, nor all the neighborhood relations between the elements, the springs and the particles when these relations are necessary to implement topological modifications.

For example, to cut an object between two adjacent hexahedra (i.e. two hexahedra sharing a face), particles of the face could need to be duplicated (as we will see later with Fig. 4). To detect if a particle needs to be duplicated or not, we need to search for existence of a path of adjacent elements joining the two cut hexahedra. This is not possible directly in the graph. Thus, a data-structure has to be added to represent elements by storing indices of its particles and indices of its adjacent elements. Then, we also need to add to each particle the indices of its incident elements. These additional data-structure and links need to be initialized and updated coherently after each operation. Moreover, other data-structures and links will be required to implement another topological modification.

The main interest of using a topological model is that we are sure that all the topological relations are described and can be retrieved directly. Existing topological constraints ensure the topological validity of the described objects and can be used to test the consistency of implemented operations. Moreover, we can use many existing and proved topological operations that can serve as basic tool to implement high level operation in our physical simulation.

OUR NEW SOLUTION

In this paper, we propose an improvement of the LCC+MSS model. At this time, we consider only hexahedral elements due to the physical models used, but the structure is suitable for any topology. Our model follows two new principles: (1) to store the darts of each element in 3-cells of the 3-map;

(2) to associate the physical information directly in the corresponding cells of the 3-map.

Storing of the darts in their element

For the first principle, we add in each element an array called dart containing all the darts of the element grouped by faces i.e. an array dart[6][4] is associated to each hexahedron. The initialization of these arrays is performed at the construction of the 3-map using the convention given in Fig. 3. With this structure, all points of a given hexahedron can be retrieved thanks to the stored darts. For example, point

X 0 is obtained from dart[0][3] and point X 1 is obtained from dart[4][2]
. Moreover, each point of an hexahedron can be retrieved by three different darts. For example X 0 can also be retrieved thanks to

dart[1][0] and dart[4][1]. f 5 f 3 f 2 f 1 f 4 0 X 5 X 7 X 2 X 1 X 3 X f 0 6 X 4 X
Figure 3: Convention used to store darts of a hexahedron composed of points X i with i ∈ {0, . . . , 7}. f j = dart[j][0] with j ∈ {0, . . . , 5} correspond to the first darts of each face (the bolded arrows).

One major interest of our proposal is to provide a direct access to all the incidence relations between cells, allowing a direct access in O(1) to information required by the different physical models. This direct access is done either by using the relations given by the 3-map or by the relations given by the dart array.

Storing of physical information in the corresponding cells

The second principle is based on the fact that the physical information associated to the i-cells (with i ∈ {0, 1, 2, 3}) of the 3-map is depending on the physical model chosen. Consequently, the direct association of the physical information to the corresponding cells allows to consider easily different physical models. In this subsection, we present the use of our proposal for a mass-spring system and a tensor-mass model, but any physical model can be envisaged. We terminate this subsection by showing the generic computation of the forces accumulation thanks to this second principle. Naturally, the computation of the acceleration of the object and its integration in order to obtain its velocity and its position is the same for all physical objects who are based on the newtonian dynamics.

Mass-Spring System. For a mass-spring system based on the Baudet formulation [BBJ + 09], particles are associated to 0-cells thanks to 0-attributes, springs are associated to 1-cells thanks to 1-attributes, and inner springs are associated to 3-cells thanks to 3-attributes. Indeed, we remember that in the considered MSS, each hexahedron includes 4 inner springs s i with i ∈ {0, 1, 2, 3}. By convention, we fix spring s 0 (resp. s 1 , s 2 , s 3 ) between particles (P 0 , P 7 ) (resp. (P 0 , P 7 ), (P 1 , P 6 ), (P 2 , P 5 ) and (P 3 , P 4 )) where particles P i correspond to point X i with i ∈ {0, . . . , 7} presented in Fig. 3.

To compute the forces involved by springs (supported both by 1-cells and 3-cells), each spring needs to know its two extremity particles which are the parameters of the addForce() function used for this purpose. This can be retrieved thanks to the 3-map and to the dart array, without storing explicitly the link between springs and particles.

• For each spring s associated to a 1-cell, one dart d s of the edge is directly known thanks to the 3-map (see Fig. 2). The first particle extremity of the spring is given by 0-attribute(d s ). The second particle extremity of the spring is given by 0-attribute(d s ) with d s =other_extremity(d s ) the dart associated with the second extremity of the edge. • For any inner spring s i associated to a 3-cell, its two particles extremity are directly retrieved by our convention and the dart array. For example for s 0 , its first extremity is 0-attribute(dart[1][0]) and its second extremity is 0-attribute(dart[3][1]).

Tensor-Mass model. For the tensor-mass model, nodes are associated to 0-cells thanks to 0-attributes, and the formulation of the forces applied on each node of an element are associated to 3-cells thanks to 3-attributes.

In this paper, we only consider hexahedral mesh with a non-linear behavior using the Green-Saint Venant strain-tensor, but it could easily be extended to any topology or physical behavior thanks to the formulation of the forces proposed by Faure [START_REF] Faure | An Implicit Tensor-Mass solver on the GPU for soft bodies simulation[END_REF].

To compute the forces applied on each element inside the addForce() function, each element starts by computing the displacement of their nodes using their current and initial position. Thanks to the dart array, each node of a 3-cell is retrieved in constant time.

Then, the force applied on each node of an element is directly computed using the formulation of forces generated by symbolic computation.

Generic forces accumulation. Thanks to this new way to associate physical information to a LCC, the forces applied on the object are easily computed by performing a loop on all the relevant i-cells, that is the i-cells affected by the force computation (e.g. the 1-cells and 3-cells for MSS; the 3-cells for MT) and calling the appropriate addForce() function on the associated physical information. In practice, it consists to iterate through all the enabled i-attributes (i.e. the attributes which are associated to cells of the 3-map).

Comparison with previous solution

As recalled in Section 3, the initial definition of the LCC+MSS topological-physical model [START_REF] Flechon | A generic topological framework for physical simulation[END_REF] associates physical entities to cells (similarly than in our new solution), except for inner springs which are not associated with elements. Moreover, each spring stores its two extremity particles to be able to accumulate the force of this spring to the particles.

Improvement for cutting. With this previous solution, extremity particles of springs may have to be updated after a cutting. For example in Fig. 4, the cut of the two top hexahedra leads to the split of particles v 1 and v 2 into four particles v 1 , v 1 , v 2 , v 2 . But after this cutting, springs of the top-right hexahedra are still associated with original particles v 1 and v 2 instead of particles v 1 and v 2 (see Fig. 4(a)). To solve this problem after the cutting, a post-processing is applied to update the incorrect extremities of springs (see Fig. 4(b)).

The solution proposed [START_REF] Flechon | A generic topological framework for physical simulation[END_REF] is to update all the springs having as extremity the particles involved in the cutting. This requires to store in each particle the list of the springs touching it. Then, when a particle is duplicated due to a cutting, an iteration is performed through this list to update the extremities of the springs to point to the new particles. Moreover, the list of springs is split in two parts: the first part for springs incident to v 1 and the second one for springs incident to v 1 (the new particle). This treatment requires a complex and time consuming processing.

v 3 v 4 v 1 v' 1 v 2 v' 2 (a) v 3 v 4 v 1 v' 2 v 2 v' 1 (b)
Figure 4:

In the LCC+MSS model (pictures from [START_REF] Flechon | A generic topological framework for physical simulation[END_REF]), the cutting involves the updating of the spring extremities. (a) The bold segments represent the inner springs wrongly attached to particle v 1 after the cutting. (b) The same springs updated after a postprocessing step.

With our new solution, as the darts of an element are constant whatever the topological modification applied on its neighbor elements, no change is required for cutting. For example, the darts of the two hexahedra split in Fig. 4 are still valid after the cutting. Moreover, as springs and inner springs extremities are only defined by using the darts of the hexahedra, they are still valid after the cutting without any updating. Thus our algorithms with topological modifications are simpler than Flechon and our method is faster (as shown in our experiments in Section 5) than Flechon.

Improvement for any physical models. Another advantage of our proposal concerns the ease to consider several physical models as presented in previous subsection. Indeed, as physical information are associated with cells thanks to the 3-map, the same algorithm is used to accumulate forces regardless of the physical model used (in our actual implementation, for both mass-spring system and tensor mass model). Moreover, the 3-map and the dart array are up to date after cutting without any need of post-processing. Thus, only physical properties concerning by the cutting (like the mass) have to be updated. This simplifies the integration of new physical models in our system.

RESULTS

In this section, we present results obtained in a Intel R i5 4690 CPU, 4 cores @3.50 GHz with 16 Go RAM. For the following tests, we consider a gravity put to g = -9.8 m.s -2 and beams with a density ρ = 1000 Kg.m -3 , a Young's modulus E = 400 KPa and a Poisson's ratio ν = 0.3. Time for cutting. We start by a comparison in time for cutting between our approach and that of Flechon [START_REF] Flechon | A generic topological framework for physical simulation[END_REF]. In Fig. 5, the curves represent the time in secondes involved to cut a beam compared to the number of faces disconnected. We consider beams discretized in 200 × 200 × 2, 400 × 400 × 2, 600 × 600 × 2 and 800 × 800 × 2 elements of size 1 × 1 × 1 m with a cutting according to the Z axis. In We observe that our method is nearby twice as fast than the LCC+MSS model as less information has to be updated in our approach to take into account a topological modification of the object. Precision. Fig. 6 presents a comparison for the traction test between the results obtained using the MSS and the MT approach compared to the analytic solution with beams discretized in 5 × 5 × 20 and 10 × 10 × 40 elements. We consider the analytic solution given by

y = ρglh 24 E I (4Lx 3 -6L 2 x 2 -x 4 )
with L = 1 m the length, h = 0.25 m the height, l = 0.25 m the thickness of the beams and I = lh 3 /12 the inertia moment. Fig. 7 presents a visual comparison of the equilibrium state obtained with the 10 × 10 × 40 MSS beam and the 10 × 10 × 40 MT beam.

As expected, we observe that the MT simulation is closest to the analytic solution than MSS simulation. This experiment shows the interest of having several physical models in a same software to simply comparisons between different simulations. Time of execution. Fig. 8 presents a comparison in time for the traction test performed on beams with several discretization. The aim was to vary the number of elements of size 0.025 × 0.025 × 0.025 m. We present sequential and parallel time with a parallelization of the simulation's loop performed using the Intel R TBB (Threading Building Blocks) C++ template library for task parallelism. Four threads were involved for the parallel time. We obtained an average speedup of 1.45 for the MSS and 2.73 for the MT.

The MT simulation is slower than the MSS simulation due to the more complex physical equations. However its parallelization has a better speedup than the MSS one. This is due to the fact that MT uses only once type of attribute for the physical information (associated with 3-cells), while MSS uses two different types of attributes (1-cells for springs and 3-cells for inner springs) requiring two different loops.

REPRODUCE REAL EXPERIMENTS

In this section, we present the use of our model to simulate experiments performed on real tissues. Firstly, we explain the original experiments made on porcine tissues [START_REF] Nicolle | Shear mechanical properties of the spleen: experiment and analytical modelling[END_REF]. Secondly we present our simulation and the preliminary results obtained with the MSS and the MT model. Note that the comparison between two physical models in this simulation is possible thanks to our improvement of the physical-topological model.

Experiment. The original experiments were performed on a porcine liver disks of 2 -3.7 mm thickness and diameter of approximately 15 mm. Fig. 12 presents a schematic picture of the experiment. The disks of porcine tissues were glued between the plates of a rotational rheometer. Then, the top plate was rotating at a frequency varied from 0.1 Hz to 2 Hz, while the bottom one was constrained. This experiment is equivalent to a local deformation of the porcine disks. The authors [START_REF] Nicolle | Shear mechanical properties of the spleen: experiment and analytical modelling[END_REF] of the original experiments derived a non-linear model from the results obtained with the experiments. However, since we only consider the linear case, the following equations are modified accordingly.

The model is represented by the relaxation modulus G(t) and the stress-strain relationship σ [ε(t)]. We have:

G(t) = K Γ(1 -n) t -n
with K the consistency, n the linear constitutive index and t the time [START_REF] Kelly | Fractal ladder models and power law wave equations[END_REF][START_REF] Nicolle | A strain-hardening bi-power law for the nonlinear behaviour of biological soft tissues[END_REF]. The Gamma function is defined by

Γ(x) = ∞ 0 t x-1 e -t dt.
Considering the experiments made on liver tissues, we used as values K = 861 Pa.s n and n = 0.1138 [START_REF] Nicolle | Shear mechanical properties of the spleen: experiment and analytical modelling[END_REF].

The stress-strain relationship is defined for the linear case by

σ [ε(t)] = t 0 G(t -t )ε(t )dt ,
where σ represents the stress, ε(t) the strain and ε a constant strain rate. We considered the same three strain rates as the authors [START_REF] Nicolle | Shear mechanical properties of the spleen: experiment and analytical modelling[END_REF] with ε ∈ {0.0151 s -1 , 0.133 s -1 , 0.67 s -1 }.

Simulation. We used our model to reproduce the experiments by simulation. We assumed that locally the tissue is represented by a rectangular block of infinitely long plates. The bottom of the tissue did not move and the top of the tissue only moved on the abscissa X. For this, we applied a direct displacement to the top particles and studied the response of the tissue. The simulated tissue consists of 3,751 particles with 10×10×30 volumes (see Fig. 13). The unit of a cube is 1 mm. Results. We considered porcine liver which is a noncompressible tissue with E = 1000 Pa, ν = 0.499 and ρ = 1000 Kg/m 3 . The gravity force is set to zero as it is negligible in the real experiment. As we performed a shear experiment in the linear case, the strain is small enough to be negligible. For this reason, we consider the stress-time relationship to analyse our simulations.

Fig. 14 shows the stress according the time obtained for the MT model and the MSS model. For the simulations, the stress is computed by dividing the force applied on the object by the area of the top plate of the tested tissue (in our case 300 mm 2 ). As expected, we observe that MT and MSS have a linear behavior for small deformations.

In the future work, we will improve the MSS solutions by implementing volume preservation and applying some correction forces. We also want to implement a new damping solution to ensure the correctness of the model for linear case.

In this paper, we presented an improvement of the LCC+MSS model proposed by Flechon [START_REF] Flechon | A generic topological framework for physical simulation[END_REF]. Our improvement is based on two new principles: (1) to store the darts of each element in the 3-cells of the 3map;

(2) to associate the physical information directly in the corresponding cells of the 3-map. Thanks to these two new principles, our model allows simpler algorithms since no more modifications are required for topological modification operations. Secondly, our model gives better results since less operations are required. Lastly, our model is more generic since it is now easier to add a new physical model while keeping all the existing operations. We have illustrated these improvements in our experiments with better results obtained for cutting than the previous model. Moreover, two different physical models (namely the mass-spring system and the tensor-mass model) have been presented who are both based on the same unified topological model. With these two physical models, we have illustrated the generic computation of the accumulation of the forces thanks to our two new principles.

In future work, we are considering adding others physical models to easily compare the results of simulations.

The most interesting would be to implement the Finite Element Method and give thought to the Position Based Dynamics approach, like the one of Müller [START_REF] Müller | Position based dynamics[END_REF].

According to the MSS, we plan to include the volume preservation constraint in order to improve its precision. We plan to study the work of Aubert [START_REF] Aubert | Volumepreserving Space Deformation[END_REF] which created a space deformation model (called DOGME), to deal with the volume preservation, and the more recent article of Duan [DHC + 14] who treats about preserving volume by position corrections. Furthermore, to improve the precision of our simulations, it is also important to take a deeper look into constraining the particles to distribute correctly the acting forces. Lastly, we want to extend our simulation of soft tissues to the case of non-linear mechanical behavior to consider bigger deformations. All these future work will benefit of the LCC+MSS model allowing to define simpler and more efficient algorithms to perform topological modifications during simulation for any kind of physical objects.

Fig. 1

 1 represents two adjacent cubes described by a 3map. Darts are displayed as arrows. β 1 (d) is the next dart following dart d in the same face and the same volume. β 0 (d) is the previous dart in the same face and volume. β 2 (d) gives the other dart from d belonging to the same edge, the same volume but not the same face. β 3 (d) gives the other dart from d belonging to the same edge, the same face but not the same volume.

Figure 1 :

 1 Figure 1: Two adjacent cubes described by a 3-map. Each cube has 6 × 4 darts. Purple disks represent vertices associated with points to obtain a 3D LCC.

Fig. 2 Figure 2 :

 22 Fig.2illustrates an example of a LCC+MSS model for two adjacent cubes.

Figure 5 :

 5 Figure 5: Comparison for the time involved for cutting.

Figure 6 :

 6 Figure 6: Comparison between the MSS, the MT and the analytic solution for a traction test.

Figure 7 :

 7 Figure 7: Comparison between the MSS (left) and the MT (right) for a traction test performed on beams discretized in 10 × 10 × 40 elements.

Figure 8 :

 8 Figure 8: Comparison in time for the MSS and MT simulation performed in sequential or in parallel with TBB.

Figure 9 :

 9 Figure 9: Cutting made during a MT simulation of a beam (subject to gravity) fixed to its two extremities.

Figure 10 :

 10 Figure 10: Cutting made during a MSS simulation. Some particles were fixed and the gravity was applied. The density of the letters was put to ρ = 200 Kg.m -3 . Fig. 11 shows several cuts applied on a giraffe mesh (with no simulation) that contains 85,405 elements. The different parts were separated by hand to obtain a proper visualization.

Figure 11 :

 11 Figure 11: Ten successive cuts on a giraffe mesh.

Figure 12 :

 12 Figure 12: Schematic representation of the experiments made on porcine tissues [NNP12].

Figure 13 :

 13 Figure 13: MSS composed of 10 × 10 × 30 hexahedra to reproduce the experiments made on tissues.

Figure 14 :

 14 Figure 14: Evolution of stress in time for respectively the MT model (top) and the MSS model (bottom) for three different strain rates.

Table 1

 1

	Mesh	# faces cut Flechon Our model
	Frog	112	4.0	2.0
	Elephant	592	14.5	7.6
	Mushroom	987	32.2	17.7
	Giraffe	5210	176.3	91.6

, more complex meshes are used for the comparison in time (in ms) for cutting: a frog (2,354 elements), an elephant (125,754 elements), a mushroom (27,000 elements) and a giraffe (85,405 elements).

Table 1 :

 1 Comparison of time (in ms) for cutting.

Acknowledgments

This work was supported by the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). The authors would like to thank the following people for their support: E. Flechon for providing the original version of TopoSim (software implementing the LCC+MSS model), X. Faure and F. Jaillet for providing the code to compute the forces generated in the mass-tensor model.