Céline Dubois

Steven Le Corre
email: steven.le-corre@ec-nantes.fr

Malek Zarroug

Patrick Rozycki

Nicolas Moës

Impact on highly compressible media in explicit dynamics using the X-FEM

Keywords: X-FEM, Explicit dynamics, Nonlinear elasticity, Remeshing technique, Numerical methods

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Many problems involving very large deformation cannot be efficiently computed by the same Lagrangian mesh throughout the simulation. As the material undergoes severe strains, elements aspect ratios deteriorate [START_REF] Benson | Computational methods in lagrangian and eulerian hydrocodes[END_REF]. This leads to a significant loss of accuracy in the approximation. The mesh is no longer able to represent a good approximation: it fails to satisfy the completeness condition. In worst cases, several elements may invert, the Jacobian determinants become negative, which may lead to the loss of positive definitiveness of the stiffness matrix and thus important local errors that may abort the calculation. Another drawback of the distorsion in explicit dynamics is that the stable time step is considerably reduced [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF].

A classical solution to avoid these problems is to use arbitrary Lagrangian Eulerian methods (ALE). Contrary to purely Lagrangian methods, the mesh motion is independent of the material motion and can be chosen in a way to eliminate distortion. These methods seem very attractive since they take advantage of both Lagrangian and Eulerian descriptions. Nevertheless their practical use is not straightforward. The choice of a judicious mesh motion can be a quite difficult task, particularly due to the fact that boundary nodes or material interfaces nodes in the case of multi material problems can only move in tangential direction [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF]. The extended finite element method (X-FEM) is known to be an efficient solution to deal with arbitrary discontinuities like cracks, voids and inhomogeneities [9, [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]. Within this method, the geometry of the discontinuity is independent of the mesh and is represented by a level-set function [START_REF] Sethian | Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and material science[END_REF]. Several authors proposed to couple X-FEM with ALE to solve bi-material problems [START_REF] Khoei | An extended arbitrary Lagrangian-Eulerian finite element method for large deformation of solids mechanics[END_REF][START_REF] Vitali | An extended finite element formulation for contact in multi-material arbitrary Lagrangian-Eulerian calculations[END_REF]. The results are satisfactory. In ALE methods, in addition to the mesh motion two other essential points require special care. The first one is that the topology of the mesh does not change during the simulation, the number of nodes remains the same so it has to be chosen carefully from the beginning of the calculation. The second one is that ALE equations involve convective terms to take into account transport of material through the mesh. The computation of these terms need the use of advection algorithms which can be computationally expensive.

Another class of methods, referred by Benson [4] as the "rezoning method", can be seen as a good alternative to ALE methods. The basic idea is to run a calculation with a Lagrangian mesh until it becomes excessively distorted and then generate a new undistorted mesh. The solution fields are mapped from the old distorted mesh to this new mesh, which is used as a Lagrangian mesh to continue the calculation. We can notice that the remeshing work has to be done on the deformed configuration in order to ensure a good quality for the new mesh [START_REF] Lee | Errors indicators and adaptive remeshing in large deformation finite element analysis[END_REF]. The new mesh can be built by adapting the old one. Some places can be refined by dividing elements but then some transition elements have to be introduced for the mesh to be conform. This process can really become heavy in three dimensions [START_REF] Cugnon | Automatisation des calculs éléments finis dans le cadre de la méthode -p[END_REF] so that it is often decided to replace the old mesh by a completely new one. This last method is also a better way to get a new mesh with good quality shape.

However, the difficulty arises from the fact that the new mesh has to approximate the boundary of the deformed body as good as possible. Within this context, X-FEM can be seen as a relevant solution since the mesh does not have to conform to the external interfaces. The new mesh can simply be a regular and structured mesh independent of the shape of the deformed body. The only condition to satisfy is that the deformed body lies inside the structured new mesh. This idea has been recently used within the Eulerian Finite Cover Method (FCM) to treat quasi-static equilibrium problems with large deformation [START_REF] Terada | Eulerian finite cover method for quasic-static equilibium of hyperelastic bodies[END_REF]. The authors propose to remesh with a regular and structured mesh after each Newton-Raphson solving step. The regular mesh is the same during the entire simulation, that is why it is called the "Eulerian mesh".

In this work, we use the possibility offered by X-FEM to remesh on structured meshes in the context of dynamics simulations. The remeshing process is not done at each time step but only when the current mesh is no longer able to handle the deformation. The new structured mesh does not have to remain the same during the whole simulation, it can be changed for a finer one for example. We present a specific algorithm to map the quantities from the old mesh to the new one. The main point of our work is to take advantage of the regularity of the new mesh to save computational time. Another great advantage of this method is that the remeshing operation can be automatized. Contrary to most rezoning methods available now, it does not require any user intervention.

X-FEM in non-linear explicit dynamics

Governing equations

Let us consider a body in its initial configuration 0 . The vector X denotes the position of a material point in this reference configuration. The motion of the body is described by x = φ(X, t) where x is the position of the material point X at time t. The displacement of a material point can be written as u(X, t) = x -X = φ(X, t) -X. The deformation gradient is defined by

F = ∇ 0 φ(X, t) = ∂ x/∂ X, ( 1
)
where ∇ 0 indicates the gradient operator with respect to material coordinate X. The Jacobian of the transformation is J = det F. The body is subjected to a force per unit mass b. The strong form of the balance of momentum can be expressed as follows:

∇ 0 • P + ρ 0 b = ρ 0 ü in 0 , ( 2
)
where P is the nominal stress, ρ 0 the density in the reference configuration and ü the acceleration of a material point. The displacement is prescribed as u = u on 0,u and the nominal traction vector as n 0 • P = t 0 on 0,t . Initial conditions are also imposed: u(X, 0) = u 0 (X) and u(X, 0) = u0 (X).

The weak form of the momentum equations including natural (Neumann type) boundary conditions is therefore given by:

0 δu • ρ 0 üd 0 + 0 δ F T : Pd 0 = 0 ρ 0 δu • bd 0 + 0,t δu • t 0 d 0 , (3)
where test functions δu belong to a space of C 0 functions which vanishes on the boundary of imposed displacements 0 u . The first term in the left-hand side of the above equation represents the virtual inertial work, the second term is the virtual internal work. The whole right-hand side term is the virtual external work.

The mechanical behavior of the material is supposed to follow an hyperelastic law. As this kind of constitutive model is path-independent, the second Piola-Kirchhoff stress S can be expressed as the derivative of a stored energy function ψ:

S = 2∂ψ/∂ C where C = F T • F, ( 4
)
is the right Cauchy-green strain tensor. For numerical examples, we use the Neo-Hookean constitutive model in which the stored energy is defined by

ψ(C) = λ(ln J ) 2 /2 -μ ln J + μ(trace(C) -3)/2, ( 5
)
where λ and μ are the Lamé constants. The second Piola-Kirchhoff stress is then given by:

S = μ(I -C -1 ) + λ(ln J )C -1 , ( 6
)
and P can be obtained by P = S • F T .

Spatial discretization with non-matching meshes on interface

The X-FEM relies on the idea of enriching the finite-element approximation by additional functions through the concept of partition of unity [START_REF] Melenk | The partition of unity finite element method: basic theory and applications[END_REF]. The enrichment functions can be used to model discontinuities inside the elements like cracks, voids or inhomogeneities and also to improve accuracy in problems where the behavior of the solution field is known analytically, like stress around crack tip for example. We consider a domain 0 divided into finite elements. Let N = {n 1 , n 2 , . . . , n m } be a set of m nodes in the mesh and let 0 g ⊂ 0 be the region where some enrichment is needed. The general form of the X-FEM displacement approximation can be expressed as follow:

u h (X) = n I ∈N φ I (X)u I + n J ∈N g φ J (X)ψ(X)a J , (7)
where the nodal set N g is defined as

N g = {n J /n J ∈ N , ω J ∩ 0 g = ∅}
where ω J is the support of the nodal shape function φ J (X), i.e. the union of all elements containing n J as a node. The additional degrees of freedom are denoted as a J , they usually do not have a physical significance but they allow one to deal with discontinuities. The function ψ(X) is called the enrichment function.

Modeling holes is a somewhat special since there is no need for additional degrees of freedom. The classical shape functions for the displacement are simply multiplied by an enrichment function V (X) whose value is 1 if the node lies outside the void and 0 if it is inside [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite element method[END_REF] The displacement field can be written as

u h (X) = n I ∈N φ I (X)V (X)u I = n I ∈N φ I (X)u I . ( 8
)
Nodes whose support is entirely inside the void are removed from the calculation (see Fig. 1). The interface between the material and the void is implicitly described as the zero level of a level-set function ϕ(X). This function gives the signed distance of any point X to the interface: ϕ(X) is positive if X lies in the void, negative in the material. As the other fields, ϕ is interpolated on the mesh by:

ϕ(X) = n I ∈N φ I (X)ϕ I (9)
We can notice for example that for linear tetrahedral elements, the level-set function is linear within each element. Substituting the trial and test functions from their approximation [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF] in the weak form (3) and using the arbitrariness of test functions, the following semi-discrete system of equations is obtained:

M • ü + f int = f ext (10)
Mass matrix terms M I J are given by:

M I J = 0 ρ φ I φ J d 0 (11)
The internal nodal force vector f int I is expressed as:

f int I = 0 ∂ φ I ∂ X • P d 0 . ( 12
)
The external nodal force vector f ext I is expressed as:

f ext I = 0 φ I ρ 0 b d 0 + 0 φ I t 0 d 0 (13)
In practice, instead of using φ I (X), elements that intersect the voids are cut along the iso-zero of the level-set and then partitioned in sub-elements. The weak form is only integrated in the sub-elements and elements that belong to the matter part [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite element method[END_REF] (see Fig. 2).

Time discretization

The numerical time integration is performed through the well known central differences explicit scheme. The mass matrix is usually diagonalized by a lumping technique so that the discrete momentum equations do not require the solution of any system of equations. The method is stable if the time step Fig. 2 For elements crossed by the iso-zero of the level-set, integration is done on sub-elements inside matter remains lower than a critical value t crit . In a linear finite element analysis without damping, one can demonstrate that

t crit = 2 ω max ≤ min l e c e , ( 14
)
where ω max is the maximum frequency of the system, l e is a characteristic length of element e in deformed configuration and c e is the current wavespeed in element e. For nonlinear problems, no stability theorem can be found. The critical time step is calculated for the linearized system and then a reduction factor is applied. Rozycki et al. [START_REF] Rozycki | X-fem explicit dynamics for constant strain elements to alleviate mesh constraints on internal or external bounadraies[END_REF] showed that the critical time step in the X-FEM could be the same as the one in classical FEM provided the use of a specific lumping technique. For a complete review of the mass lumping strategies for X-FEM, one can also refer to Menouillard et al. [START_REF] Menouillard | Mass lumping strategies for x-fem explicit dynamics: application to crack propagation[END_REF].

We use the central difference algorithm as detailed in Belytschko et al. [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF] where the time step can be variable. As pointed out by the authors, this is more suitable for finite strains since the stable time step changes during the simulation as the mesh deforms and as the deformation of the material makes the wave speed change. The time increment between the iteration n and n + 1 is denoted by t = t n+1t n . Each time increment can be decomposed as follows:

1. Update of nodal displacements:

u n+1 = u n + t un + t 2 2 ün , (15)
2. First partial update of nodal velocities:

un+ 1 2 = un + t 2 ün , ( 16
)
3. Computation of nodal accelerations:

ün+1 = M -1 • (-f int,n+1 + f ext,n+1 ), ( 17
)
4. Second partial update of nodal velocities:

un+1 = un+ 1 2 + t 2 ün+1 . ( 18
)
3 Remeshing on a regular grid

The main issue in remeshing is to handle data transfer between two meshes. We have to define a way to transmit the solution fields. Two kinds of data are distinguished in the literature: variables stored at nodes and variables stored at integration points. For nodal variables, collocation methods are widely used. They can be decomposed in three steps:

1. Identification: for each new node, find in which old element it is located; 2. Inversion: find the local coordinate of the new node in the old element; 3. Interpolation: use the interpolation on the old mesh to find the respective nodal value of the new node.

Other methods based on averaging procedure can be used to transfer nodal fields between incompatible meshes. For example, the least-square method allows to minimize the square of the difference between two interpolations of the same field. This method is also known as the L2 projection and is widely used in dynamics to transfer initial conditions on a mesh [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF]. Another method based on Lagrange multipliers and called the Mortar method is often used to transfer fields between interfaces in domain decomposition problems or contact problems. In the case of incompatible meshes, this method is very close to the least-square method. Opposite to the collocation method, the least-square and the Mortar methods do not imply any loss of information when transferring from a fine mesh to a coarse one. Moreover they allow to keep the duality between two fields on the new mesh. However, these averaging methods require the inversion of a global mass matrix and the construction of integral cross products of shape functions defined on different meshes. Despite their advantages in terms of accuracy, these methods can make the data transfer very time consuming in an explicit context. Furthermore, in our case, the boundaries of the new and the old mesh will not match exactly so that such projection methods do not seem easily applicable. That is why we choose to use the collocation method. For internal variables, usually known at Gauss points, the task is much more difficult because they have to stay compatible with other variables after the projection. This has to Fig. 3 X-FEM allows to remesh the deformed configuration with a non-conforming regular mesh be taken into account in elasto-plastic models for example. Several authors [START_REF] Boroomand | Recovery procedure in error estimation and adaptivity. Part II: adaptivity in nonlinear problems of elasto-plasticity behaviour[END_REF][START_REF] Mediavilla | A robust and consistent remeshing-transfer operator for ductile fracture simulations[END_REF] proposed methods to ensure the equilibrium during the data transfer. We will not further detail these methods since our study focuses on hyperelastic bodies.

For nodal fields, a similar problem can arise. In explicit dynamics if we transfer all kinematics nodal fields: displacements, velocities and accelerations, we cannot ensure that the conservation of linear momentum (3) will be satisfied on the new mesh. To avoid these problems, we choose to transfer only displacements and velocities. The acceleration field is then deduced from the above by means of the linear momentum equation [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF]. The procedure can therefore be summarized as:

1. Projection of u n+1 () and un+ 1 2 () from mesh 1 to mesh 2; 2. Computation of ün+1 on mesh 2 by (17).

Scan conversion technique

The use of the X-FEM allows to remesh the deformed configuration with regular and non-conforming grids (cf. Fig. 3). We can take advantage of the regularity of the new mesh to deal efficiently with the identification step required for the collocation.

Scan conversion (or rasterization) is a standard technique in computer graphics for displaying filled polygons on a grid of pixels [START_REF] Breen | 3D scan conversion of csg models into distance, closest-point and color volumes[END_REF]. The goal of the method is to determine efficiently pixels which lie inside polygons. Let us consider a convex polygon on a regular 2D grid in the plane (x, y). We begin by classifying the edges of the polygon according to the highest y coordinate of their extremities. For each grid row that intersects the polygon, we determine which edges are intersected and find the left and right intersection points x L and x R (Fig. 4). We can then conclude that each grid point between x L and x R is inside the polygon. As pointed in [START_REF] Breen | 3D scan conversion of csg models into distance, closest-point and color volumes[END_REF] the computational complexity of the scan conversion algorithm is O(e + r + n) where e the number of edges of the polygon, r the number of rows that intersects the polygon and n the number of grid points inside. If the sides of the polygon are not smaller than the grid spacing, r and e become negligible compared to n and the complexity reduces to O(n).

function ScanPolygonConvex(Polygone P) sort the edges of P by their maximum y coordinate; for each scanline y (from top to bottom) do select the two edges that intersect the scanline; compute the intersecting points x L et x R ; add all nodes with a x coordinate between x L and x R to the set of nodes scanned by P; end Algorithm 1: Scan conversion algorithm for a twodimensional convex polygon.

For the 3D Case, the former 2D scan conversion technique is adapted and limited to the case of tetrahedra in a uniform regular grid. We first define the intersection of the considered tetrahedron with a horizontal plane at height z, which therefore forms a 2D polygon (triangle or quadrangle) in the considered plane. This polygon is then scanned by the 2D technique described in the previous section.

Figure 5 represents the slicing of a tetrahedron to form polygons. The complexity of the algorithm is also linear in the number of grid points inside the polyhedron.

function ScanTetrahedron(Tetrahedron T) sort the faces of T by their maximum z coordinate; for each scanplan z (from top to bottom) do select the faces that intersect the scanplan; build the polygon P resulting from the intersection between the scanplan and the selected faces; apply the 2D function ScanPolygonConvex to P; add all nodes scanned by P to the set of nodes scanned by T; end Algorithm 2: Scan conversion algorithm for a tetrahedron.

In our problems, we scan every deformed old elements and sub elements to find which nodes of the new regular mesh are inside. The complexity of the algorithm is linear in the number of nodes of the new mesh recovered by the old mesh.

Building of the level-set on the new mesh

In order to continue the computation after a remeshing step, we also need to locate the position of the level-set on the new mesh. Opposite to the displacement, the level-set is not a Lagrangian field but an Eulerian one. The level-set is known for a given position and not for a given material point. When the mesh moves, the level-set stored at a node still represents the distance between this node and the interface in the initial configuration but not in the deformed one. When the initial mesh is deformed, the level-set field is then no longer a distance function. If this field is transferred directly by collocation, it will lead to a wrong positioning of the interface for the new mesh. This can be shown by a simple example in one dimension in Fig. 6.

One could imagine to reinitialize the level-set to a distance function on the first mesh in its deformed configuration before doing the transfer to the second mesh. The reinitialisation or redistanciation of a level-set is usually done by finding the level-set φ satisfying the above transport equation:

∂φ ∂t + sign(φ o )(|∇φ| -1) = 0 ( 19
)
where φ(x, t = 0) = φ o is the initial value of the level-set. In order to alleviate oscillations phenomenon, this equation is solved by means of specific algorithm (ENO, WENO) which can be computationally expensive. More efficient scheme like the Fast Marching Methods can also be used to accelerate the process, but anyway it will require additional iterations which are not suitable for explicit dynamics.

Another way to obtain a distance function is to discretize the interface and to simply evaluate the distance using its definition:

D(x) = min y |x -y|sign(φ o ) ( 20
)
This method can be directly applied on the new mesh; there is no projection to do. That is the reason why we choose this solution. Moreover, we propose here a more efficient strategy using the information obtained during the scan to rebuild the level-set on the new mesh. First of all, we know if a node has been scanned by a element (or sub-element) in the matter or in the void. We can conclude on the sign of the level-set for this node: negative in the matter and positive in the void. We then have to evaluate the distance. During the integration, the initial mesh has been cut along the iso-zero of the level-set, so we have access to a mesh of the deformed contour. This mesh is called the boundary mesh and is made of triangular facets. For each node of the new mesh, we can loop over boundary facets, calculate the distance to each facet and then select the minimum distance and the closest facet. The cost of this simple algorithm is O(N * F) where N is the number of new nodes and F the number of facets of the boundary mesh. In 3D, this can be expensive and useless because we only need to know the level-set near its iso-zero. We propose here a fast algorithm which takes advantage of the informations obtained during the scan conversion phase.

The first step of this algorithm is to select a group of nodes on the new mesh where the level-set function has to be known. This group will be called the "new narrow band nodes". A simple way to define these nodes is to consider that they are the nodes scanned by old elements close to the level-set. The first task is therefore to select a group of elements close to the level-set in the old mesh which will be called the "old narrow band elements", as illustrated in Fig. 7. This set is made of old elements cut by the iso-zero of the level-set function that represents the boundary of the structure. In order to be sure to have a sufficiently large band, we also include the neighbors of the cut old elements which lie inside the matter.

Then the "new narrow band nodes" can be constructed by taking the nodes scanned in this band of old elements (Fig. 8a). Because of robustness issue that will be described later, the neighboring nodes of these nodes are also included to the new narrow band. This is illustrated on Fig. 8b.

Doing so, we have limited the number of nodes where the level-set will be evaluated. To further improve this procedure, one can also reduce the number of facets to be taken into account for each node. Each facet has been created by the cutting operation of an element of the old mesh so a link Fig. 7 "Old narrow band" elements can be established between the facets and the elements of the "old narrow band elements". A facet is set to be potentially close to the element that it crosses and to the neighbor of this element (see Fig. 9). Thus, an element of the "old narrow band elements" is associated to a group of eventually close facets. Defining those local groups, any node scanned in the narrow band can be related to a small amount of facets which makes the cost of the calculation of the level-set proportional to the number of nodes in the narrow band.

Algorithm robustness

An important issue is to ensure the reliability of the proposed algorithm. First of all, we can easily check that the algorithm proposed allows to obtain exactly the same results when choosing a second mesh identical to the first one.

A second important point is to ensure that no node is forgotten during scan operation. Such an error can come from the fact that our algorithm employs floating point values. Rounding operations during the scan conversion can lead to 9 A facet is potentially close to its "parents" elements and to his neighbors miss a node very close to the edge of an element. We have to introduce a tolerance when we test whether a node lies into an old element. In practice, it can be seen as a small increase in the element size (Fig. 10). Thus, for example when a new node is exactly on one edge common to two old elements, it is scanned by both. In the future, a good improvement would be to design an algorithm using only integers as it is usually done in computer graphics.

In the distance calculation, we have to take care to get all facets possibly close to a node. This is the reason why we take the neighbors of the old elements to construct the narrow band of a facet. This can avoid the critical case illustrated in Fig. 11, when the closest facet to a node does not cross the element by which it was scanned.

Another critical situation that could happen in the X-FEM context is that new nodes having to be considered on the new mesh can be unrecovered by the old mesh. If they are not scanned, no displacement or velocity fields can be transferred to them (see Fig. 12). This case is avoided by the fact that we include the neighbors of the nodes scanned in the "old narrow band elements" to construct the "new narrow band nodes". If new nodes are not scanned, we use their closest Fig. 10 Introduction of a tolerance in order to ensure that no nodes are forgotten Fig. 11 The closest facet to a node does not cross the element by which it was scanned facet to find which is their closest deformed element and then we extrapolate the solution fields (see Fig. 13).

Numerical examples

In this section, we study several examples in order to test the validity of the algorithm described above. The goal of the

Example 1: free fall of a sphere

In order to test the efficiency of our algorithm, we begin with a very basic problem: the free fall of a sphere. The value of the sphere radius is 2.4 m. The problem is first discretized with a regular tetrahedral mesh whose characteristic length is 0.4 m. This mesh was built from a hexahedral one and each hexahedron was cut into six tetrahedra. Two different cutting directions were used for two consecutive hexahedra in order to obtain an isotropic mesh. During the explicit computation, the initial mesh is replaced by a new one built in a similar way but with a characteristic length of 0.2 m (see Fig. 14).

The sphere undergoes a rigid body motion. Displacement and velocity fields are uniform. Data transfer is therefore straightforward. We see in Fig. 15 that those fields perfectly match the analytical solution before and after remeshing.

The first step of the level-set computation is to select the "old narrow band elements" on the initial mesh. This band is made of the elements crossed by the level-set iso-zero and their neighbors (see Fig. 16). This group of elements cover nodes on the new mesh. These nodes and their neighbors make up the "new narrow band nodes" (see Fig. 17). The new iso-zero of the level-set is computed only on this small amount of nodes. In this example, there are 12,980 nodes in the "new narrow band" whereas the initial mesh contained 24,389 nodes. The computational cost to get the level-set is divided by two compared to a straight forward approach.

The second simplification of our algorithm then consists in reducing the number of facets tested to obtain the closest one. The iso-zero of the level-set initially contains 3,456 facets. The number of facets tested for each node of the "new narrow band" varies from one node to another between 20 and 200 facets and depends on the curvature of the iso-zero. An average number of about 100 facets are tested per node.

Thanks to these two simplification processes, the total number of nodes times facets investigated for the levelset computation is about 1,298,000. With the straight forward approach, the distance computation would have been done at 24,389 nodes with 3,456 facets, which leads to 84,288,384 computations. Hence, the presented method enables to reduce the computational effort by a factor 65. This first example is very promising. For industrial problems, we can expect that the number of nodes near the level-set is very small compared to the total number of nodes on the new mesh so that the benefit of our algorithms will be much more significant.

Example 2: dynamic traction test

In order to show the relevance of the data transfer approach, we propose in this section a simulation on a nearly cubic part subjected to dynamic loading. The dimensions of the part are 1 m × 1 m × 1.05 m. The boundary conditions are a complete clamping of the lower face and a sudden traction of 45,000 N on the upper face. The lateral faces can only move in the traction direction. This case can be compared to the solution of a rod under a state of uniaxial strain. For the material, we choose the Neo-Hookean model described in Sect. 2.1, whose parameters are typical of a polymeric foam: ρ = 22 kg m -3 , λ = 66,428 Pa and μ = 16,607 Pa.

The part is discretized with a mesh of characteristic length 0.1 m. The loaded face is not meshed, but implicitly represented by the iso-zero of a level-set. We compare the results obtained without remeshing and with five successive remeshing steps. At each remeshing step, the new mesh used to discretize the deformed configuration is chosen equal to the initial mesh. Displacements of a point on the loaded face obtained with and without remeshing are compared in Fig. 18. Results for both simulations are nearly identical. A small gap appears between the two displacements from We also compare global variables obtained from both simulations. Figure 19 shows the evolution of kinetic and internal energies. The results are also quite similar. Small gaps appear for both energies at different instants. The projection sometimes increases the gap and sometimes decreases it but always remains in an acceptable range.

In order to further analyze the influence of the number of projections, we compute the error of kinetic energy, internal energy (see Fig. 20) and the error on the part volume (see Fig. 21) for different numbers of remeshing steps. We can see that the effect of the remeshing is rather random. The remeshing does not lead to an automatic decrease or increase of energy and volume. Errors induced from a remeshing can be erased by the next one. For some particular numbers of projections (especially 5 and 7), one can notice errors of about 10% which could be worrying. Nevertheless, this is due to the particular design of this test where the part upper boundary is straight and raises progressively on a regular mesh. The projection error is therefore amplified because several elements are cut in the same way. For a more complex geometry, this phenomena is expected to be perfectly random so that the error will not be so exaggerated. The approximation induced by the projection technique is therefore acceptable.

Example 3: dynamic compression test

We consider here the same part than in the previous example but we change the load to a compression of 90,000 N on the upper face. The mesh is chosen identical to the previous case. Displacements obtained without remeshing can be seen in Fig. 22 where a compression ratio of 60% was reached.

123 In order to estimate the distortion of the mesh, we use a simple aspect ratio quality factor defined by:

Q = α h max ρ = α h max S 3V ( 21
)
where α = √ 6

12 is a normalization coefficient allowing to have a quality equal to 1 for an equilateral tetrahedron, h max is the length of the longest edge, ρ is the radius of the inscribed sphere, S is the sum of face surfaces and V is the volume. Our mesh is made of tetrahedra with an initial aspect ratio of 1.71. During the compression stage, from the beginning to t = 0.016 s, the maximum aspect ratio increases and gets close to 5 so that the accuracy of results for these elements is no longer guaranteed. Furthermore, the critical time step, which is directly linked to the elements shape, decreases dramatically, starting from 3.5 × 10 -4 to 3.87 × 10 -5 s. This can be excessively time consuming in an industrial context.

To highlight the way our remeshing strategy enables to overcome those two problems, different situations of remeshing were tested: 2, 3 and 4 successive remeshing steps during the compression stage. For each test, the projection mesh used to remesh the deformed part is identical to the initial one. As shown in Fig. 23, the distortion goes back to its initial value at each remeshing, which is expected because of the perfect regularity of the new mesh. Furthermore, the value of the critical time step is maintained at higher values thanks to the remeshing (Fig. 24). Even if the time step does not go back to its initial value (the wave velocity is not the same in This last test finally validates the principle of the proposed technique and demonstrates its efficiency in the context of dynamic simulations implying large deformations. Nevertheless, the considered situations are still geometrically rather simple and the next section will focus on more complex deformations induced by impacts on compressible bodies.

Impact problems

Impact situations are present in many industrial problems, from manufacturing processes to the crash of automobiles. It seems important to include the modeling of these phenomena in our method.

During an impact, one or several bodies stay in contact until the vanishing of waves results in a release of contact forces. We can distinguish two main approaches to treat contact constraints: the Lagrange multiplier method and the penalty method. These two methods have been used as the basis of most contact algorithms [START_REF] Wriggers | Computational contact mechanics[END_REF].

The main principle of Lagrange methods is to ensure the equality of the displacements at the contact interface between two bodies by means of Lagrange multipliers, which can be seen as normal tensions. The advantage of this approach is the exact respect of the non-penetration condition but its drawback is that it requires the introduction of new degrees of freedom. Furthermore, in explicit dynamics, the Lagrange multiplier method implies the solution of an implicit system at each time step [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF], which can be time consuming if contact interfaces are large. With the usual penalty method, the normal tensions are approximated by forces which are proportional to the penetration. The major drawback of this method is that some interpenetration of the contacting bodies is allowed. Nevertheless, in the case of large deformations we consider here, this approximation can be considered as precise enough. Despite this little drawback, this method is easy to implement and computationally efficient because it does not involve any additional degrees of freedom. The major problem induced by this method is that in contrast to Lagrange multipliers method, the penalty method always decreases the critical time step [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF]. The choice between the two methods depends on the accuracy needed and on the time of computation wanted. In our case, the use of the X-FEM has also to be considered to make a choice.

In the X-FEM context, contact interfaces are not meshed. This particularity has to be taken into account in the contact algorithms. The use of the X-FEM with the Lagrange multipliers method requires special care to satisfy the inf-sup condition needed for mixed formulations. Geniaut et al. [START_REF] Geniaut | A stable 3D contact formulation for cracks using x-fem[END_REF] and Béchet et al. [START_REF] Béchet | A stable lagrange multiplier space for the stiff interface conditions within the extende finite element method[END_REF][START_REF] Moës | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF] have proposed specific algorithms for statics problems to define the appropriate Lagrange multipliers discretization space. For large deformation dynamics problems, the search of this discretization space has to be Fig. 28 One dimensional penalty test case: evolution of the critical time step according to the fraction of matter with added mass for different values of the penalty factor k Fig. 29 One dimensional impact test done at each time step since the contact interfaces are moving so that the computational cost of this operation can become very important in an explicit framework. Furthermore, when remeshing is needed, we cannot ensure that the transfer of Lagrange multipliers leads to the conformance with the infsup condition on the new mesh. In this context, the penalty method seems to be a more suitable choice.

In the following section, we will describe the use of the penalty method in explicit dynamics in classical finite elements. We will then emphasize the special treatments due to the introduction of the X-FEM. Finally, we will discuss on its reliability through some numerical examples.

Penalty method

The penalty method originates in optimization theory. Its application in contact problems consists in adding the above term to the weak form of the linear momentum:

δG P = c k 2 δ(g 2 N )H(g N )d ( 22
)
where g N is the normal gap function between two bodies, c is the contact interface, k is the penalty parameter and H(g N ) is the Heaviside function:

H(g N ) = 1 if g n > 0 0 if g n < 0
The non-penetration condition is supposed to be better satisfied as the penalty parameter k is increased. The penalty method (Fig. 25) can be seen as the addition of a spring between two bodies in order to keep them separated. It increases the stiffness of the elements at the boundary in contact. In the context of explicit methods, this leads to a decrease of the critical time step [START_REF] Belytschko | Contact-impact by the pinball algorithm with penalty and lagrangian methods[END_REF]. This feature can be easily demonstrated in 1D considering the impact of a truss on a rigid wall. The stability is governed by the element in contact associated with the penalty spring.

The mass and the linearized stiffness matrices are respectively given by:

M = ρ Sl 2 1 0 0 1 and K = E S l 1 -1 -1 1 + k
where k = l E S k . By applying the eigenvalue element inequality, we obtain the FEM critical time step t FEM c :

t FEM c = l c 2 2 + k + √ 4 + k 2 (23)
The critical time step decreases as the inverse of penalty parameter k. For example, if k is chosen equal to 10 (i.e. the penalty spring is ten time more stiff than the material), the time step is decreased by 58%.

4.2 Use of the penalty method in explicit dynamics with X-FEM

Stability

Within the X-FEM, the penalty traction has to be applied on a surface which is not meshed. For a one dimensional case this can be modeled by attaching the penalty spring at the boundary of the material part of the X-FEM element.

As clear from Fig. 26, the boundary element is crossed by the interface at a distance l of its length l. For such a case,

M = 1 2 ρ Sl 0 0 K = E S l + (1 -) 2 k -+ (1 -)k -+ (1 -)k + 2 k
where represents the fraction of matter in the X-FEM element.

Applying again the eigenvalues criterion, the critical time step is obtained by:

t c = l c 2 √ 2 + k -2k + 2 2 k + √ 4 2 -8 2 k + 8 3 k + k 2 -4k 2 + 8k 2 2 -8k 2 3 + 4 4 k 2
As illustrated on Fig. 27, if the material fraction tend to zero, then the critical time step tend to zero. The meshing procedure with X-FEM is independent of the position of the boundary, so that the parameter is not controlled and the values of the critical time step can be extremely low. This is not acceptable for industrial applications. To overcome this drawback, we propose to add mass on the X-FEM element. The added mass corresponds to the complementary mass of the "void part" of the X-FEM element. The mass matrix is therefore identical to the one obtained in finite element whereas the stiffness remains the one of the X-FEM element:

M = 1 2 ρ Sl 1 0 0 1 K = E S l + (1 -) 2 k -+ (1 -)k -+ (1 -)k + 2 k
The critical time step is then given by:

t c = l c 2 2 + k -2k + 2 2 k + √ 4 2 -8 2 k + 8 3 k + k 2 -4k 2 + 8k 2 2 -8k 2 3 + 4 4 k 2
When tend to zero, the critical time step tend to a bounded limit equal to l c 2 k . In Fig. 28, we plot the X-FEM critical time step normalized by the one obtained in finite element with penalty t FEM c . We can observe that this ratio is always greater than one. Thanks to mass addition, the X-FEM time step stays greater than the one in finite elements.

This technique is a mass scaling method which has the particularity to add a naturally bounded quantity of mass. This additive mass is low compared to the mass of the complete structure so it is expected not to affect significantly the results.

One dimensional impact example

In order to check the reliability of the proposed method, we consider the 1D case of two linear elastic rods moving with equal speed in opposite directions (see Fig. 29). We call g the half distance initially separating the two rods.

The numerical values are given in Table 1.

We compare the results obtained by classical FEM and by X-FEM. For both cases, each rod is discretized with 20 elements. In the X-FEM case, the impacting boundaries of the rods are not meshed: the fraction of matter for the impacting element of the rod on the left-hand side is denoted as 1 , whereas it is denoted as 2 for the rod at the right-hand side (see Fig. 30).

The stability of the explicit scheme is governed by the two impacting elements properties associated with the penalty spring. The time steps will differ a bit from the previous study since here the two impacted bodies are deformables.

For the FEM case (meshed boundaries), the computation of the critical time step is based on the eigenvalues of the above matrices:

M = ρ Sl 2 ⎡ ⎢ ⎢ ⎣ 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 ⎤ ⎥ ⎥ ⎦ K = E S l ⎡ ⎢ ⎢ ⎣ 1 -1 0 0 -1 1+ k -k 0 0 -k 1 + k -1 0 0 -1 1 ⎤ ⎥ ⎥ ⎦
which lead to the following FEM critical time step:

t FEM c = l c 2 1 + k + √ 1 + k 2 (24)
For the X-FEM case without mass scaling, the element matrices depend on 1 and 2 :

M = ρ Sl 2 ⎡ ⎢ ⎢ ⎣ 1 0 0 0 0 1 0 0 0 0 2 0 0 0 0 2 ⎤ ⎥ ⎥ ⎦ K = E S l ⎡ ⎢ ⎢ ⎣ 1 + k(1 -1 ) 2 -1 + k 1 (1 -1 ) -k 2 (1 -1 ) -k(1 -1 )(1 -2 ) -1 + k 1 (1 -1 ) 1 + k 2 1 -k 1 2 -k 1 (1 -2 ) -k 2 (1 -1 ) -k 1 2 2 + k 2 2 -2 + k 2 (1 -2 ) -k(1 -1 )(1 -2 ) -k 1 (1 -2 ) -2 + k 2 (1 -2 ) 2 + k(1 -2 ) 2 ⎤ ⎥ ⎥ ⎦
An analytical formula of the critical time step for the X-FEM would be difficult to analyse. We just represent in Fig. 31 the evolution of the time step for different values of 1 and 2 with and without mass scaling. This time step is normalized by the time step obtained with classical finite element t FEM c (see Eq. 24). As in the previous study, the X-FEM time step tends to zero when 1 or 2 tend to zero. This drawback disappears when using the mass scaling technique: the X-FEM time step even becomes greater than the FEM one.

For a case where 1 = 0.1 and 2 = 0.9, we compare the results obtained by FEM and X-FEM. The mass scaling for X-FEM leads to a mass addition of 2.61% of the total mass. The evolution of displacements, velocities, accelerations and contact forces at the two contacting boundaries for FEM and X-FEM are represented in Figs. 32 and33.

We can see that displacements, velocities and contact forces obtained with X-FEM are very close to the classical FEM solutions. However, some differences between the results of the two methods can be observed for the acceleration. Analytically, the acceleration is always zero except at the impact and at the release instants where it tends to infinity. Obviously, such an evolution cannot be represented numerically: one only gets maximum values at impact and release times as well as numerical oscillations between these two states. These oscillations are minimized in the X-FEM thanks to the mass addition.

In Fig. 34, we can also observe the differences between X-FEM and FEM in the evolution of kinetic energy. As expected, the addition of mass with the X-FEM leads to a small increase of kinetic energy at the impact time. However, this increase does not involve significant differences between X-FEM and FEM since the mass addition is limited to 2.61% of the total mass. This example shows that the X-FEM with mass addition allows to obtain accurate results while keeping a reasonable time step.

Remeshing during impact

We have demonstrated in the previous paragraph the possibility to treat impact problem within X-FEM. We know focus on the effects of the remeshing process on the impact treatment, therefore combining the remeshing technique of Sect. 3 and the present penalty technique.

Procedure

As mentioned in Sect. 3, for a remeshing occurring at time t n+1 , we only transfer the predictors of displacements and velocities: u n+1 and un+ 1 2 . The acceleration ün+1 is evaluated by resolution of the equilibrium equation on the second mesh. The velocity un+1 is deduced from this acceleration.

The contact condition will be taken into account in the computation of ün+1 . In explicit methods, the acceleration at time t n+1 does not affect the displacement at the same time. That is why Carpenter et al. [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF] pointed that the contact condition has to be treated for the forward increment. This means that the computation of ün+1 will use contact forces depending on the penetration predicted for time t n+2 . We can conclude that there is no contact force to transfer between the old and the new meshes in explicit dynamics.

Due to the change of discretization between both meshes, the penetration on the second mesh at time t n+1 will differ from the penetration at the same time on the first mesh. In a Lagrange multiplier method, the impenetrability condition is ensured on the first mesh whereas interpenetration appears in the second one due to projection. The computation of Lagrange multipliers on the second mesh would be necessary to overcome this issue. In the penalty method, some interpenetration is allowed, so there is no reason to try to eliminate it on the second mesh. Thus, the accuracy of the results will not be affected by the remeshing procedure with the penalty method.

Numerical example

We now consider a 3D example involving remeshing during impact. We aim to simulate the impact of a sphere against a rigid plane. The sphere, of initial radius 0.3 m has an initial Fig. 33 One dimensional impact test: evolution of acceleration (a) and contact force (b) at contacting boundaries, comparison between the FEM and X-FEM solutions for 1 = 0.1, 2 = 0.9 and k = 0.1 velocity of 10 m s -1 . We use a regular isotropic mesh, of characteristic length 0.05 m, made of tetrahedra built according to the process described in Sect. 3.4.1. The geometry of the sphere is described by a level-set function and the plane is not discretized but described by an analytical function. The initial gap between the sphere and the plane is set to 0.005 m so that contact occurs at the beginning of the computation. During the impact, the deformed sphere is remeshed by our projection method on a mesh identical to the initial one as described in Fig. 35.

The wave propagation due to impact is shown in Figs. 36 and 37, where one can see the distribution of the displacements in a cutting plane of the sphere at different time steps.

In Fig. 38, we compare the evolution of kinetic and internal energies with and without remeshing. These energies are also compared to the one obtained with ABAQUS code. We can observe a fairly good agreement between the three simulations. The remeshing leads to a small increase of the kinetic energy and a small decrease of the internal energy. These gaps seem reasonable compared to the amount of energy involved in this problem and to the other discretization errors.

The penetration is almost the same from one mesh to another after the remeshing. This is illustrated in Fig. 39 where the displacement of an impacting point was plotted.

Finally, we can observe the benefit coming from the remeshing on the distortion rate (Fig. 40) and on the critical time step (Fig. 41). Without remeshing, the aspect ratio increases to values close to 4 whereas it is kept under 2.8 by the remeshing operation. Furthermore, as visible from Fig. 41, the remeshing also allows to keep a bigger time step while the material undergoes severe deformations. On this example, at the most deformed state (t = 0.019 s), the time step is multiplied by 4.

Conclusion

In this paper, we propose a new strategy to deal with impact problems on highly compressible media in the framework of explicit dynamics. The basic idea is to use the X-FEM to simplify the remeshing work needed for this kind of problems. This simplification relies on the fact that the X-FEM allows to remesh with meshes that do not match the shape of the deformed structure. A simple structured mesh can be used each time whenever a remeshing is needed.

The data transfer between successive meshes is done by a collocation method. An important step in that process is to locate the new nodes inside the old mesh. The fact of Efficiency and accuracy of the proposed algorithms are tested on several examples. We particularly demonstrate that the remeshing strategy allows to limit the distortion and to keep a reasonable time step when the material is submitted to severe strains. The method therefore brings all the satisfying improvements one expects from a remeshing process. Naturally, compared to other remeshing techniques, our use of regular grids can be considered as non-optimal and yet somewhat academic. The extension to a more adapted mesh, which will be the subject of further works, leads to heavier calculations.

In order to be able to deal with more realistic physical situations, the treatment of an impacting rigid body was also implemented, the contact problem being treated through a penalty method. As the use of such a method with X-FEM in explicit dynamics leaded to a decrease of the critical time step, we proposed a specific mass lumping strategy for the X-FEM elements to overcome this issue. Excellent agreement with reference solution results are obtained for 1D and 3D examples.

Fig. 1

 1 Fig. 1 Nodes of the grid used in the X-FEM approach

Fig. 4

 4 Fig. 4 Principle of the scan conversion of a polygon in 2D. a Case where edge1 and edge2 are activated. b Case where edge2 and edge3 are activated

Fig. 6

 6 Fig.6The level-set field is no longer a distance function in the deformed configuration. If no reinitialisation step is done, collocation leads to a wrong positioning of the iso-zero

Fig. 8

 8 Fig. 8 Construction of the "new narrow band" nodes. a Phase 1. b Phase 2

Fig. 12

 12 Fig. 12 Some non-scanned nodes have to be taken into account for the computation on the new mesh

Fig. 13 Fig. 14 Fig. 15

 131415 Fig. 13 Management of data transfer for unscanned nodes: (a) finding the closest, (b) extrapolation of the fields values

Fig. 16 Fig. 17 2 Fig. 18

 1617218 Fig. 16 Old narrow band elements on mesh 1. a 3D view. b Cut view

Fig. 19

 19 Fig. 19 Dynamic traction test: evolution in time of kinetic (a) and internal energies (b)

Fig. 20 Fig. 21

 2021 Fig. 20 Dynamic traction test: evolution of the relative error of the kinetic (a) and the internal energy (b) with the number of remeshing steps

Fig. 22 Fig. 23

 2223 Fig. 22 Dynamic compression test: displacement of one point of the loaded face

Fig. 24

 24 Fig. 24 Dynamic compression test: evolution in time of the critical time step for different remeshing strategies

Fig. 25

 25 Fig. 25 One dimensional approach of the penalty method

Table 1 1 Fig. 30

 1130 Fig. 30 One dimensional impact test: representation of the X-FEM mesh

Fig. 31

 31 Fig. 31 One dimensional impact test: evolution of the X-FEM critical time step with (a) and without (b) mass scaling the mass and the linearized stiffness matrices are then given by:

Fig. 32

 32 Fig. 32 One dimensional impact test: evolution of displacements (a) and velocities (b) at contacting boundaries, comparison between the FEM and X-FEM solutions for 1 = 0.1, 2 = 0.9 and k = 0.1.

Fig. 34

 34 Fig.34One dimensional impact test: evolution of kinetic energy, comparison between the FEM and X-FEM solutions for 1 = 0.1, 2 = 0.9 and k = 0.1

Fig. 35

 35 Fig. 35 Impact of a sphere: state of the sphere at the beginning of the simulation on mesh 1 (a) and at the beginning of the simulation on mesh 2 (b)

Fig. 39

 39 Fig. 39 Impact of a sphere: evolution of the displacement at a contacting point

Fig. 40

 40 Fig. 40 Evolution of maximum aspect ratio

Fig. 41

 41 Fig. 41 Evolution of the critical time step

Acknowledgments

We would like to thank PSA Peugeot-Citron for their financial support and more particularly M. Laurent Rota for his scientific advice.


Diana Lupan
email: diana.lupan@cti.pub.ro

Mihai Dasc Lu
email: mihai.dascalu@cs.pub.ro

Tr U An-Matu

Philippe Dessus
email: philippe.dessus@upmf-grenoble.fr

Analyzing Emotional States Induced by News Articles with Latent Semantic Analysis

Keywords: emotional state, Latent Semantic Analysis, automatic evaluation of news articles

Emotions are reflected both in verbal and written communication. If in the first case they can be easier to trace due to some specific features (body language, voice tone or inflections), in the second it can be quite tricky to grasp the underlying emotions carried by a written text. Therefore we propose a novel automatic method for analyzing emotions induced by texts, more specifically a reader's most likely emotional state after reading a news article. In other words, our goal is to determine how reading a piece of news affects a person's emotional state and to adjust these values based on his/her current state. From a more technical perspective, our system (Emo2 -Emotions Monitor) combines a context independent approach (actual evaluation of the news employing specific natural language processing techniques and Latent Semantic Analysis) with the influences of user's present emotional state estimated through his/her specific feedback for building a more accurate image of a person's emotional state.

Introduction

Emotions are one of the defining elements of human nature, as they are present in daily life and almost in any context. They add value to our interactions in a unique way and their induced effect can change the meaning of an entire message.

Although emotions are observably present in face-to-face or verbal communication, they can also be traced in written messages. If in the first case emotions are mostly expressed using body language and different voice features (for example, tone, rhythm, frequency), in the second one all these features are eliminated, as the communication channel cannot support them. Therefore, it can be difficult to detect emotions in written communication due to the lack of characteristics and our sole elements of analysis are words, grouped into sentences. Despite this fact, analyzing emotions in written texts can provide a better understanding of the underlying message by predicting its most likely impact on the person reading it.

Our approach covers the automatic evaluation of a person's mood and sentiment change when he/she reads an article. Our focus is to analyze short news consisting of a headline and a brief description that are written with the intention of "provoking" emotions and attracting the reader's attention.

This paper covers two approaches implemented in the Emo2 (Emotions Monitor) system: context independent and environment-influenced evaluation of emotional states. The context independent approach implies an evaluation of the article's content with different natural language processing techniques from the perspective of the intended emotional state induced by the author of the article, while the environmentinfluenced approach takes into account the user's current emotional state for adjusting the results, based on the influence of other similar articles to one's current emotional state.

The remainder of the paper first presents similar related works and a general overview of the architecture describing each module. The fourth section focuses on the context independent evaluation process including details about the word valence and article similarity approaches. The fifth section presents the environmentinfluenced evaluation method, while the last sections comprise the evaluations, results and conclusions of our survey.

Related Work

The UA-ZBSA system [START_REF] Kozareva | UA-ZBSA: A Headline Emotion Classification through Web Information[END_REF] aimed at classifying a set of news headlines into six types of emotions: "anger", "disgust", "fear", "joy", "sadness", "surprise". Our approach is similarly comprised of six types of emotions, but we modified them to best fit our dimensions of analysis. Therefore, the included categories include "in control of the situation", "fear", "joy", "sadness", "excitement" and "boredom" described in detail in the next sections). UA-ZBSA's emotion classification uses a bag of word approach comprising of frequency and co-occurrence counts of concepts collected from the World Wide Web. It applies the idea that "words which tend to co-occur across many documents with a given emotion are highly probable to express that emotion". Moreover, it considers that adjectives with the same polarity tend to appear together. On the other hand, UPAR7's design [START_REF] Chaumartin | UPAR7: A knowledge-based system for headline sentiment tagging[END_REF] is based on the assumption that all the words from an article's title potentially carry emotions. The objective of the system was to identify the expression which carries the main topic of the title as this expression has major importance. Moreover, rules for detecting specific emotions were also developed (ex. a negation may be an indicator of a possible surprise). The system consists of the SS-Tagger (a Part-of-Speech tagger) and the Stanford Parser and used WordNet, SentiWordNet and WordNet-Affect as lexical resources. It also counts lexical elements that are later on classified as good indicators of surprise: negations, modal auxiliaries and question marks.

3

General Overview of the Functional Architecture

The designed system consists of several modules, grouped into processing layers that are presented in Fig. 1:

Fig. 1. System's Architecture
Firstly, regarding the overall processing algorithm, we must determine a person's initial or current emotional state. Therefore, our application tries to determine as accurately as possible the user's current emotional state by analyzing the feedback given by him on a set of 5 news articles, with regards to his/her induced current state. Although these initial news might alter a person's current emotional state, the results are later on used for adapting and personalizing future automatic results provided by the system through regressions in order to best reflect the actual impact of each new article on the user's current state.

From a higher level perspective, our analysis is divided into two separate steps: the initial, objective assessment, performed in the "Core evaluation" module, and the subjective, adaptive component of personalizing the results by considering the current emotional state of the user (integrated within the "Estimate Induced Emotional State" component).

The proposed evaluation method uses multiple inputs: the enriched ANEW database, a trained LSA semantic vector space, a corpus of articles expressing pure states and the news currently analyzed. The enriched ANEW database [START_REF] Bradley | Affective norms for English words (ANEW): Stimuli, instruction manual and affective ratings[END_REF] consists of concepts (words) that have three values associated expressing its emotional state (Happy/ Unhappy, Calm/ Excited, Controlled/ In-Control). The Happy/ Unhappy dimension is quite straightforward as it expresses a state of happiness or sadness. Calm/Excited refers to the intention of the user to keep reading about the information presented in the current news article. The last dimension underlines the feeling of being in control of the situation just read versus being controlled by the events (or persons) described in the present news. The values for each dimension are in the [START_REF] Bradley | Affective norms for English words (ANEW): Stimuli, instruction manual and affective ratings[END_REF][START_REF] Chaumartin | UPAR7: A knowledge-based system for headline sentiment tagging[END_REF] range and the initial database of approximately 1,000 words was enriched using an algorithm that consists of 2 stages that are presented in detail in the following section.

On the other hand, the evaluation process can be divided in two approaches: semantic similarity evaluation between texts by means of LSA and mining key concepts. The later consists of extracting key concepts from both the body and the title of an article and combining their "emotional values" in order to determine the overall induced emotions. Therefore, a set of reference documents (initially the corpus consisted of 10 "pure" documents per emotional state) is used to obtain the semantic dimension of the analysis, while key-concepts mining offers the lexical dimension of our analysis. Key-concept mining is applied separately on the title and on the body of news because we empirically considered that the title contains more valuable information than the rest. The title is usually designed to be short and to attract the reader's attention; moreover, we can consider that the title expresses from a human perspective the essence and the synthetic message of the body. Thus, due to their intrinsic value, the key-concepts belonging to the title have a greater contribution in the formula than the ones extracted from the body.

In the process of evaluating the emotional state by means of Latent Semantic Analysis, the result is determined by using the average value of all documents found similar with the content of the currently assessed article, weighted by the similarity between them. By using the similarity of concepts that induce a certain state, we can imply that documents with similar structure in terms of coverage of key concepts inflict a similar emotional state. The previous statement is not generally valid, but under the consideration of comparing news articles from a specific domain, within a short timeframe, our results from this approach proved to be quite accurate and to level the spikes of different dimensions of certain news articles within a predefined timeframe. In essence, news articles covering the same topic in a short timeframe rarely change their messages in terms of induced emotions.

The values obtained from the two perspectives are combined so the final result approximates, in the best way possible, the emotional state that is most likely to be induced when a person reads the article. The "Current emotional state" module integrates the user's input in our analysis (a questionnaire and 5 news articles) by estimating his present state and the fluctuations observed after reading each article. This set of articles contains texts extracted from the latest news feed into our application so their influence on the assessment of other similar articles can be considered relevant and significant. In our experiments, articles consisting of the title and the content (typically 2 or 3 phrases) were extracted from CNN.com's RSS (Really Simple Syndication) feeds.

Context-Independent Evaluation

The context-independent approach is based on two different stages: the first compares the concepts found in the analyzed article with terms stored in our enriched ANEW database (Word valence approach), while the second determines the similarities between the article and a set of predefined documents that express a dominant emotional state (Article similarity). In the end, the two approaches are combined during the article evaluation step for obtaining the corresponding valences on all 3 dimensions as a final score for each analyzed news article.

Word Valence Approach

Each term is associated with three connotations that express what kind of feelings it arouses when someone reads it. These valences expressed on the 3 dimensions (Happy/ Sad, Calm/ Excited, Controlled/ In-Control) are into a local database, formed through an enrichment process that is described in detail in the following paragraphs.

Initially, the number of stored terms was around 1,000 and the associated values were obtained from the ANEW database (Affective Norms for English Words, [START_REF] Bradley | Affective norms for English words (ANEW): Stimuli, instruction manual and affective ratings[END_REF]) that was developed to provide a set of normative emotional ratings for a large number of words in English. This initial database was afterwards increased by using WordNet [START_REF] Fellbaum | WordNet: An Electronic Lexical Database[END_REF] synsets and Latent Semantic Analysis -LSA [START_REF] Landauer | An introduction to Latent Semantic Analysis[END_REF][START_REF] Wiemer-Hastings | How latent is Latent Semantic Analysis?[END_REF] that enabled the determination of similarities between a "known" concept and a "new" one. By combining metrics of similarity from a lexicalized ontology with cosine similarity between concepts of the semantic vector space [START_REF] Manning | Foundations of statistical Natural Language Processing[END_REF], our results were improved and we were able to better weigh the valence of words not present in the initial ANEW database.

LSA was trained on the Reuters corpus of news that provided adequate materials, relevant from the perspective of topics with regards to currently assessed articles and general enough to grasp underlying connections between concepts. The main reasons of using LSA are that its behavior has been proven to be close to that of human beings [START_REF] Foltz | Latent semantic analysis for text-based research[END_REF] and that it can be successfully used to enrich data obtained from experiments with human subjects and to approximate how reading a news article affects someone's emotional state [START_REF] Wang | Sentiment Classification of Documents Based on Latent Semantic Analysis[END_REF].

The actual enrichment algorithm consists of two stages. During the first stage we select a word and determine its synonyms using synsets from WordNet, whereas in the second stage similar concepts, from similar contexts, are determined through LSA. By considering the cosine similarity as the semantic metric for comparing two concepts, we can introduce new terms into our enriched database based on "known" concepts (already existent in the database), their corresponding values and the similarity with the term to be added. Moreover, only "known" concepts that are similar over a certain threshold are considered synonyms and are taken into account in computations. Therefore, based on the synonyms identified through WordNet of the concepts already present in the ANEW database and on the semantic similarity with k most nearest concepts, we determine a new word's valence and introduce it in our enriched database. After running multiple iterations with increasing k values, in our final experiments we chose k to be the three most dominant and similar concepts by means of cosine similarity to the currently assessed one.

Article Similarity

As previously mentioned, this approach consists of determining similarities between the content of a news article with a corpus of documents that express a pure emotional state. Although the concept of "pureness" is subject to personal evaluation, these selected documents express a dominant dimension, identified and agreed upon by several human evaluators. The initial corpus used within our experiments comprises of 10 general interest news articles per state and LSA was used to compute the similarity between the currently evaluated article and each of the documents within the corpus.

Article Evaluation

The final result for each article is obtained by combing the values of each of the two previous approaches. Firstly, the concepts from the title are considered to be more valuable, in terms of expressed emotions, than the ones from the body.

The evaluation formula from the word valence approach is the following:

(

where and are a linear combination of valences of keyconcepts contained within the title and the body of the article.

The value of p was determined experimentally using multiple iterations within the range [0.55, 0.9], with increments of 0.05. The experiments revealed that 0.6 is the value that best combines the message expressed by the title with the one implied by the content and, in this case, the best overall correlation was obtained between the computer-based results and the human evaluators. On the other hand, the content of an article is considered to be similar to a text from the set if the similarity exceeds a threshold of 0.2 (the similarity can take values in the [-1, 1] range). For each emotional state, the arithmetic mean between the values that meet this requirement is computed. Next, the results are linearly normalized in order to be projected on a [START_REF] Bradley | Affective norms for English words (ANEW): Stimuli, instruction manual and affective ratings[END_REF][START_REF] Chaumartin | UPAR7: A knowledge-based system for headline sentiment tagging[END_REF] value range. The final result is computed using the following weighted formula which is designed to give more credit to the value obtained from analyzing the key concepts than to the value determined through article similarity, intrinsically limited by the bag of words approach:

(2)

The best fit for p was experimentally determined by testing decreasing values within the range [0.55, 0.9] and the experiments revealed that p = 0.65 leads to the best combination of the values associated with the word valence approach and with the semantic similarity between news and pure emotional state documents.

Environment-Influenced Evaluation

The environment-influenced approach comes as an adjustment to the results obtained in the first stage (context-independent evaluation) and can be seen as a feedback loop used for further personalizing the results provided by our system. Our assumption is that the emotional state of a person is differently influenced by a news article depending on his/hers current state and a set of preconceptions (mentality). For example, a person who is initially happy and reads sad news will be less affected and therefore less sad than a person who was initially unhappy. Moreover, the effects of reading the news also depend on specific features of the person reading the article, such as psychological profile, how geographically close he/she is to the events presented in the news. Taking into account these observations, some additional features which rate these aspects were introduced. Firstly, the current emotional state is assessed by asking the reader to rate 5 news articles which are selected in order to be as similar as possible to the news that will be afterwards automatically analyzed and displayed (Fig. 2); later on, the user must answer a question referring to his current geographical location. The five presented articles are chosen by sorting all the news in the database after their similarities with the main news and picking up the top five most similar ones. Although these initial articles might change the user's current emotional state, their automatic assessment compared to the user's feedback provides valuable information on the current emotional state and offers the basis for applying regressions in order to fine-tune the results for the newly assessed articles.

Fig. 2. Gra

The current geographica the assumption that events greater proportion. Concept depicted events are located valences are augmented rel achieve this functionality, G used for determining the loc

Results

As mentioned in the previo from CNN.com. To test th displays an article and thr induced. Each chart is posi state so the user can have sample of the application's depicting the 3 dimensions ous sections, the designed system was tested on RSS fe he application, the user is prompted with a window wh ree charts that represent the most likely emotional s itioned between two images that suggest a pure emotio an accurate representation of the result. Fig. 3 present graphical user interface, including a generated radar gr of the analysis:

providing feedback on a specific article and an example o Our validation survey included 10 participants (5 women and 5 men) of age between 20 and 25 with diverse backgrounds (computer science and medical). The survey consisted of 2 stages, each having: 5 initial evaluation news (as previously mentioned, these articles are selected to be as similar as possible to the news to be automatically evaluated by the application) in order to determine the current emotional state, a question referring to the user's geographical position and 10 additional news articles as assessment for obtaining feedback from users regarding the application's results. Therefore, the tests consisted of 20 RSS feeds from all the categories mentioned above. Synthetic results, covering the three dimensions of the analysis (Pleasure Valence, Arousal Valence, Dominance Valence), are presented in Fig. 4.

The correlations between the approaches and the survey, for all dimensions, are presented in Table 1. As it can be observed from the previous table, Pleasure gets equal correlation between the two approaches and the average correlations on all three dimensions are very similar. But, by also taking into consideration the nature of each dimension, we can extrapolate from our experiment that Arousal is a more environment-dependent value in the sense that arousal may vary more than the others depending on the previous mood states, whereas Dominance is more influenced by the actual context and the topics within the article. Additionally, we obtained very high correlations between the two approaches, which was expected due to the fact that P2 refines the initial results.

On the other hand, there are some situations where the second approach obtains significantly inaccurate results in comparison to the survey (that expresses the actual emotions of the user) and also to the other approach. Besides the subjectivity of the humans' evaluation of their own emotional state, a possible explanation for this situation could be the fact that the news from the questionnaire were insufficiently similar to the ones displayed in the results; so, instead of adjusting the values from the first stage to be closer to reality, the algorithm is misguided and modifies them incorrectly, making the final results even worse. A possible solution to this problem could be enriching the news database with more texts, in order for new articles to be evaluated in comparison to other older, but more similar ones.

Additionally, although location provides some valuable input to our analysis, further factors will be addressed in future versions of our system: topics and personal interests highly impacting one's emotional state, gender and similar factors that can be used to differentiate assessment for different categories of users.

Conclusions

Our goal was to develop an application that estimates how reading a piece of news affects a person's emotional state and our approach consisted of a contextindependent and an environment-influenced analysis of emotions. As technical background, we used several natural language processing techniques, mostly centered on Latent Semantic Analysis. The results were validated through a survey and, through the provided feedback, we can consider that our method is promising and relevant, in an area where subjectivity is the central player. On the other hand, there are some aspects that can be improved; therefore we plan to perform a deeper analysis and understanding of the text and to model the user's emotional state evolution.

As future research, we envision conducting surveys with an increased number of RSS feeds and audience so the obtained feedback would be more accurate. Also, increasing our corpus of "pure" emotional texts and including more psychological factors in our adaptive model could lead to improved results.

Fig. 3 .

 3 Fig. 3. Graphical interface p generated graph

Fig. 4 .

 4 Fig. 4. Pleasure, Arousal and Dominance valences. P1 denotes the context-independent evaluation; P2 denotes the environment-influenced approach.

Table 1 .

 1 Correlations between the automatic results and the survey data

		Pleasure	Arousal	Dominance	Average on all 3
		Valence	Valence	Valence	dimensions
	Survey -P1	0.55	0.48	0.61	0.55
	Survey -P2	0.54	0.58	0.56	0.56
	P1 -P2	0.86	0.91	0.96	0.91

Acknowledgment. This work was partially supported by project 264207, ERRIC-Empowering Romanian Research on Intelligent Information Technologies/FP7-REGPOT-2010-1.