Heng Wang
email: hwang@nlpr.ia.ac.cn

Muhammad Muneeb Ullah
email: muhammad.muneeb.ullah@inria.fr

Alexander Klaser
email: alexander.klaser@inria.fr

Ivan Laptev
email: ivan.laptev@inria.fr

Cordelia Schmid
email: cordelia.schmid@inria.fr

Alexander Kläser

Evaluation of local spatio-temporal features for action recognition

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Local image and video features have been shown successful for many recognition tasks such as object and scene recognition [START_REF] Fergus | Object class recognition by unsupervised scale-invariant learning[END_REF][START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF] as well as human action recognition [START_REF] Laptev | Learning realistic human actions from movies[END_REF][START_REF] Schüldt | Recognizing human actions: A local SVM approach[END_REF]. Local space-time features capture characteristic shape and motion in video and provide relatively independent representation of events with respect to their spatio-temporal shifts and scales as well as background clutter and multiple motions in the scene. Such features are usually extracted directly from video and therefore avoid possible failures of other pre-processing methods such as motion segmentation and tracking.

Many different space-time feature detectors [START_REF] Dollar | Behavior recognition via sparse spatio-temporal features[END_REF][START_REF] Jhuang | A biologically inspired system for action recognition[END_REF][START_REF] Laptev | Space-time interest points[END_REF][START_REF] Oikonomopoulos | Spatio-temporal salient points for visual recognition of human actions[END_REF][START_REF] Willems | An efficient dense and scale-invariant spatio-temporal interest point detector[END_REF][START_REF] Wong | Extracting spatio-temporal interest points using global information[END_REF] and descriptors [START_REF] Kläser | A spatio-temporal descriptor based on 3Dgradients[END_REF][START_REF] Laptev | Local descriptors for spatio-temporal recognition[END_REF][START_REF] Laptev | Learning realistic human actions from movies[END_REF][START_REF] Scovanner | A 3-dimensional SIFT descriptor and its application to action recognition[END_REF][START_REF] Willems | An efficient dense and scale-invariant spatio-temporal interest point detector[END_REF] have been proposed in the past few years. Feature detectors usually select c 2009. The copyright of this document resides with its authors. It may be distributed unchanged freely in print or electronic forms.

spatio-temporal locations and scales in video by maximizing specific saliency functions. The detectors usually differ in the type and the sparsity of selected points. Feature descriptors capture shape and motion in the neighborhoods of selected points using image measurements such as spatial or spatio-temporal image gradients and optical flow.

While specific properties of detectors and descriptors have been advocated in the literature, their justification is often insufficient due to the limited and non-comparable experimental evaluations used in current papers. For example, results are frequently presented for different datasets such as the KTH dataset [START_REF] Dollar | Behavior recognition via sparse spatio-temporal features[END_REF][START_REF] Jhuang | A biologically inspired system for action recognition[END_REF][START_REF] Kläser | A spatio-temporal descriptor based on 3Dgradients[END_REF][START_REF] Laptev | Learning realistic human actions from movies[END_REF][START_REF] Schüldt | Recognizing human actions: A local SVM approach[END_REF][START_REF] Willems | An efficient dense and scale-invariant spatio-temporal interest point detector[END_REF][START_REF] Wong | Extracting spatio-temporal interest points using global information[END_REF], the Weizmann dataset [START_REF] Blank | Actions as space-time shapes[END_REF][START_REF] Scovanner | A 3-dimensional SIFT descriptor and its application to action recognition[END_REF] or the aerobic actions dataset [START_REF] Oikonomopoulos | Spatio-temporal salient points for visual recognition of human actions[END_REF]. For the common KTH dataset [START_REF] Schüldt | Recognizing human actions: A local SVM approach[END_REF], results are often non-comparable due to the different experimental settings used. Furthermore, most of the previous evaluations were reported for actions in controlled environments such as in KTH and Weizmann datasets. It is therefore unclear how these methods generalize to action recognition in realistic setups [START_REF] Laptev | Learning realistic human actions from movies[END_REF][START_REF] Rodriguez | Action mach: A spatio-temporal maximum average correlation height filter for action recognition[END_REF].

Several evaluations of local space-time features have been reported in the past. Laptev [START_REF] Laptev | Local Spatio-Temporal Image Features for Motion Interpretation[END_REF] evaluated the repeatability of space-time interest points as well as the associated accuracy of action recognition under changes in spatial and temporal video resolution as well as under camera motion. Similarly, Willems et al. [START_REF] Willems | An efficient dense and scale-invariant spatio-temporal interest point detector[END_REF] evaluated repeatability of detected features under scale changes, in-plane rotations, video compression and camera motion. Local spacetime descriptors were evaluated in Laptev et al. [START_REF] Laptev | Local descriptors for spatio-temporal recognition[END_REF], where the comparison included families of higher-order derivatives (local jets), image gradients and optical flow. Dollár et al. [START_REF] Dollar | Behavior recognition via sparse spatio-temporal features[END_REF] compared local descriptors in terms of image brightness, gradient and optical flow. Scovanner et al. [START_REF] Scovanner | A 3-dimensional SIFT descriptor and its application to action recognition[END_REF] evaluated the 3D-SIFT descriptor and its two-dimensional variants. Jhuang et al. [START_REF] Jhuang | A biologically inspired system for action recognition[END_REF] evaluated local descriptors in terms of the magnitude and orientation of space-time gradients as well as optical flow. Kläser et al. [START_REF] Kläser | A spatio-temporal descriptor based on 3Dgradients[END_REF] compared space-time HOG descriptor with HOG and HOF descriptors [START_REF] Laptev | Learning realistic human actions from movies[END_REF]. Willems et al. [START_REF] Willems | An efficient dense and scale-invariant spatio-temporal interest point detector[END_REF] evaluated the extended SURF descriptor. However, evaluations in these works were usually limited to a single detection or description method as well as to a single dataset.

The current paper overcomes above-mentioned limitations and provides a fair comparison for a number of local space-time detectors and descriptors. We evaluate performance of three space-time interest point detectors and six descriptors along with their combinations on three datasets with varying degree of difficulty. Moreover, we compare with dense features obtained by regular sampling of local space-time patches, as recently excellent results were obtained by dense sampling in the context of object recognition [START_REF] Fei-Fei | A Bayesian hierarchical model for learning natural scene categories[END_REF][START_REF] Jurie | Creating efficient codebooks for visual recognition[END_REF]. We, furthermore, investigate the influence of spatial video resolution and shot boundaries on the performance. We also compare methods in terms of their sparsity as well as the speed of available implementations. All experiments are reported for the same bag-of-features SVM recognition framework. Among interesting conclusions, we demonstrate that regular sampling consistently outperforms all tested space-time detectors for human actions in realistic setups. We also demonstrate a consistent ranking for the majority of methods across datasets.

The rest of the paper is organized as follows. In Section 2, we give a detailed presentation of the local spatio-temporal features included in our comparison. Section 3 then presents the experimental setup, i.e., the datasets and the bag-of-features approach used to evaluate the results. Finally, Section 4 compares results obtained for different features and Section 5 concludes the paper with the discussion.

Local spatio-temporal video features

This section describes local feature detectors and descriptors used in the following evaluation. Methods were selected based on their use in the literature as well as the availability of the implementation. In all cases we use the original implementation and parameter settings provided by the authors.

Detectors

The Harris3D detector was proposed by Laptev and Lindeberg in [START_REF] Laptev | Space-time interest points[END_REF], as a space-time extension of the Harris detector [START_REF] Harris | A combined corner and edge detector[END_REF]. The authors compute a spatio-temporal second-moment matrix at each video point µ(•; σ , τ) = g(•; sσ , sτ) * (∇L(•; σ , τ)(∇L(•; σ , τ)) T ) using independent spatial and temporal scale values σ , τ, a separable Gaussian smoothing function g, and space-time gradients ∇L. The final locations of space-time interest points are given by local maxima of H = det(µ)k trace3 (µ), H > 0. The authors proposed an optional mechanism for spatio-temporal scale selection. This is not used in our experiments, but we use points extracted at multiple scales based on a regular sampling of the scale parameters σ , τ. This has shown to give promising results in [START_REF] Laptev | Learning realistic human actions from movies[END_REF]. We use the original implementation available on-line1 and standard parameter settings k = 0.0005, σ2 = 4, 8, 16, 32, 64, 128,

τ 2 = 2, 4.
The Cuboid detector is based on temporal Gabor filters and was proposed by Dollár et al. [START_REF] Dollar | Behavior recognition via sparse spatio-temporal features[END_REF]. The response function has the form: R = (I * g * h ev ) 2 + (I * g * h od ) 2 , where g(x, y; σ ) is the 2D spatial Gaussian smoothing kernel, and h ev and h od are a quadrature pair of 1D Gabor filters which are applied temporally. The Gabor filters are defined by h ev (t; τ, ω) =cos(2πtω)e -t 2 /τ 2 and h od (t; τ, ω) =sin(2πtω)e -t 2 /τ 2 with ω = 4/τ. The two parameters σ and τ of the response function R correspond roughly to the spatial and temporal scale of the detector. Interest points are the local maxima of the response function R. We use the code from the authors' website 2 and detect features using standard scale values σ = 2, τ = 4.

The Hessian detector was proposed by Willems et al. [START_REF] Willems | An efficient dense and scale-invariant spatio-temporal interest point detector[END_REF] as a spatio-temporal extension of the Hessian saliency measure used in [START_REF] Beaudet | Rotationally invariant image operators[END_REF][START_REF] Lindeberg | Feature detection with automatic scale selection[END_REF] for blob detection in images. The detector measures the saliency with the determinant of the 3D Hessian matrix. The position and scale of the interest points are simultaneously localized without any iterative procedure. In order to speed up the detector, the authors used approximative box-filter operations on an integral video structure. Each octave is divided into 5 scales, with a ratio between subsequent scales in the range 1.2 -1.5 for the inner 3 scales. The determinant of the Hessian is computed over several octaves of both the spatial and temporal scales. A non-maximum suppression algorithm selects joint extrema over space, time and scales: (x, y,t, σ , τ). We use the executables from the authors' website 3 and employ the default parameter setting.

Dense sampling extracts video blocks at regular positions and scales in space and time. There are 5 dimensions to sample from: (x, y, t, σ , τ), where σ and τ are the spatial and temporal scale, respectively. In our experiments, the minimum size of a 3D patch is 18 × 18 pixels and 10 frames. (In Section 4.4, we evaluate different spatial patch sizes for dense sampling.) Spatial and temporal sampling are done with 50% overlap. Multi-scale patches are obtained by multiplying σ and τ by a factor of √ 2 for consecutive scales. In total, we use 8 spatial and 2 temporal scales since we consider the spatial scale to be more important than the time scale. We consider all combinations of spatial and temporal scales, i.e., we sample a video 16 times with different σ and τ parameters.

Descriptors

For each given sample point (x, y,t, σ , τ), a feature descriptor is computed for a 3D video patch centered at (x, y,t). Its spatial size ∆ x (σ ), ∆ y (σ ) is a function of σ and its temporal length ∆ t (τ) a function of τ. Dollár et al. [START_REF] Dollar | Behavior recognition via sparse spatio-temporal features[END_REF] proposed the Cuboid descriptor along with the Cuboid detector. The size for the descriptor is given with

∆ x (σ ) = ∆ y (σ ) = 2 • ceil(3σ ) + 1 and ∆ t (τ) = 2 • ceil(3τ) + 1.
We follow the authors' setup and concatenate the gradients computed for each pixel in the patch into a single vector. Then, principal component analysis (PCA) is used to project the feature vector to a lower dimensional space. We download the code from the authors' website 2 and use the default settings (e.g., the size of descriptor after PCA projection is 100). The PCA basis is computed on the training samples.

The HOG/HOF descriptors were introduced by Laptev et al. in [START_REF] Laptev | Learning realistic human actions from movies[END_REF]. To characterize local motion and appearance, the authors compute histograms of spatial gradient and optic flow accumulated in space-time neighborhoods of detected interest points. For the combination of HOG/HOF descriptors with interest point detectors, the descriptor size is defined by ∆ x (σ ) = ∆ y (σ ) = 18σ , ∆ t (τ) = 8τ. Each volume is subdivided into a n x × n y × n t grid of cells; for each cell, 4-bin histograms of gradient orientations (HOG) and 5-bin histograms of optic flow (HOF) are computed. Normalized histograms are concatenated into HOG, HOF as well as HOG/HOF descriptor vectors and are similar in spirit to the well known SIFT descriptor. In our evaluation we used the grid parameters n x , n y = 3, n t = 2 as suggested by the authors. We noticed low dependency of results for different choices of the scale factor for σ , τ in general. We use the original implementation available on-line 1 . When computing HOG/HOF descriptors for Hessian detectors, we optimized the mappings σ = ασ h and τ = β τ h w.r.t. α, β for HOG/HOF scale parameters σ , τ and the scale parameters σ h , τ h returned by the Hessian detector. For the Cuboid detector (computed at low space-time scale values) we fixed the scales of HOG/HOF descriptors to σ 2 = 4 and τ 2 = 2.

The HOG3D descriptor was proposed by Kläser et al. [START_REF] Kläser | A spatio-temporal descriptor based on 3Dgradients[END_REF]. It is based on histograms of 3D gradient orientations and can be seen as an extension of the popular SIFT descriptor [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] to video sequences. Gradients are computed using an integral video representation. Regular polyhedrons are used to uniformly quantize the orientation of spatio-temporal gradients. The descriptor, therefore, combines shape and motion information at the same time. A given 3D patch is divided into n x × n y × n t cells. The corresponding descriptor concatenates gradient histograms of all cells and is then normalized. We use the executable from the authors' website 4 and apply their recommended parametric settings for all feature detectors: descriptor size ∆ x (σ ) = ∆ y (σ ) = 8σ , ∆ t (τ) = 6τ, number of spatial and temporal cells n x = n y = 4, n t = 3, and icosahedron as polyhedron type for quantizing orientations.

Willems et al. [START_REF] Willems | An efficient dense and scale-invariant spatio-temporal interest point detector[END_REF] proposed the extended SURF (ESURF) descriptor which extends the image SURF descriptor [START_REF] Bay | SURF: Speeded up robust features[END_REF] to videos. Like for previous descriptors, the authors divide 3D patches into n x × n y × n t cells, The size of the 3D patch is given by ∆ x (σ ) = ∆ y (σ ) = 3σ , ∆ t (τ) = 3τ. For the feature descriptor, each cell is represented by a vector of weighted

Experimental setup

In this section we describe the datasets used for the evaluation as well as the evaluation protocol. We evaluate the features in a bag-of-features based action classification task and employ the evaluation measures proposed by the authors of the datasets.

Datasets

We carry out our experiments on three different action datasets which we obtained from the authors' websites. The KTH actions dataset [START_REF] Schüldt | Recognizing human actions: A local SVM approach[END_REF] 5 consists of six human action classes: walking, jogging, running, boxing, waving, and clapping (cf. Figure 1, top). Each action class is performed several times by 25 subjects. The sequences were recorded in four different scenarios: outdoors, outdoors with scale variation, outdoors with different clothes and indoors. The background is homogeneous and static in most sequences. In total, the data consists of 2391 video samples. We follow the original experimental setup of the authors, i.e., divide the samples into test set (9 subjects: 2, 3, 5, 6, 7, 8, 9, 10, and 22) and training set (the remaining 16 subjects). As in the initial paper [START_REF] Schüldt | Recognizing human actions: A local SVM approach[END_REF], we train and evaluate a multi-class classifier and report average accuracy over all classes as performance measure.

The UCF sport actions dataset [START_REF] Rodriguez | Action mach: A spatio-temporal maximum average correlation height filter for action recognition[END_REF] 6 contains ten different types of human actions: swinging (on the pommel horse and on the floor), diving, kicking (a ball), weight-lifting, horse-riding, running, skateboarding, swinging (at the high bar), golf swinging and walking (cf. Figure 1, middle). The dataset consists of 150 video samples which show a large intraclass variability. To increase the amount of data samples, we extend the dataset by adding a horizontally flipped version of each sequence to the dataset. Similar to the KTH actions dataset, we train a multi-class classifier and report the average accuracy over all classes. We use a leave-one-out setup and test on each original sequence while training on all other sequences together with their flipped versions (i.e., the flipped version of the tested sequence is removed from the training set).

The Hollywood2 actions dataset [21]7 has been collected from 69 different Hollywood movies. There are 12 action classes: answering the phone, driving car, eating, fighting, getting out of the car, hand shaking, hugging, kissing, running, sitting down, sitting up, and standing up (see Figure 1, bottom). In our experiments, we used the clean training dataset (the authors also provide an automatic, noisy dataset). In total, there are 1707 action samples divided into a training set (823 sequences) and a test set (884 sequences). Train and test sequences are obtained from different movies. The performance is evaluated as suggested in [START_REF] Marszałek | Actions in context[END_REF], i.e., by computing the average precision (AP) for each of the action classes and reporting the mean AP over all classes (mAP).

Evaluation framework

A video sequence is represented as a bag of local spatio-temporal features [START_REF] Schüldt | Recognizing human actions: A local SVM approach[END_REF]. Spatiotemporal features are first quantized into visual words and a video is then represented as the frequency histogram over the visual words. In our experiments, vocabularies are constructed with k-means clustering. We set the number of visual words V to 4000 which has shown to empirically give good results for a wide range of datasets. To limit the complexity, we cluster a subset of 100,000 randomly selected training features. To increase precision, we initialize k-means 8 times and keep the result with the lowest error. Features are assigned to their closest vocabulary word using Euclidean distance. The resulting histograms of visual word occurrences are used as video sequence representations.

For classification, we use a non-linear support vector machine [START_REF] Chang | LIBSVM: a library for support vector machines[END_REF] with a χ 2 -kernel [START_REF] Laptev | Learning realistic human actions from movies[END_REF] K(H i , H j ) = exp(-

1 2A V ∑ n=1 (h in -h jn ) 2 h in + h jn ) , (1)
where H i = {h in } and H j = {h jn } are the frequency histograms of word occurrences and V is the vocabulary size. A is the mean value of distances between all training samples [START_REF] Zhang | Local features and kernels for classification of texture and object categories: A comprehensive study[END_REF]. For multi-class classification, we apply the one-against-rest approach and select the class with the highest score.

Experimental results

This section presents experimental results for various detector/descriptor combinations. Results are presented for the different datasets in Sections 4.1-4.3. Section 4.4 evaluates different parameters for dense sampling. The computational complexity of the tested methods is evaluated in Section 4.5 Due to high memory requirements of some descriptor/detector code, we subsample original UCF and Hollywood2 sequences to half spatial resolution in all our experiments. This enables us to compare all methods on the same data. We evaluate the effect of subsampling for the Hollywood2 data set in Section 4. evaluated for other detectors than those used in original papers. Unfortunately, separate implementations of these descriptors were not available. Note that due to random initialization of k-means used for codebook generation, we observed a standard deviation of approximately 0.5% in our experiments.

KTH actions dataset

KTH actions [START_REF] Schüldt | Recognizing human actions: A local SVM approach[END_REF] 5 is to date the most common dataset in evaluations of action recognition. Among recently reported results, Laptev et al. [START_REF] Laptev | Learning realistic human actions from movies[END_REF] obtain 91.8% using a combination of HOG and HOF descriptors, while Kläser et al. [START_REF] Kläser | A spatio-temporal descriptor based on 3Dgradients[END_REF] get 91.4% with the HOG3D descriptor. Both methods use the Harris3D detector and follow the original experimental setup of [START_REF] Schüldt | Recognizing human actions: A local SVM approach[END_REF]. Adopting the Cuboid detector, Liu and Shah [START_REF] Liu | Learning human actions via information maximization[END_REF] report 94.2%, and Bregonzio et al. [START_REF] Bregonzio | Recognising action as clouds of space-time interest points[END_REF] obtain 93.2% with a 2D Gabor filter based detector. Note, however, that these results were obtained for a simpler Leave-One-Out Cross-Validation setting and are not directly comparable to results in this paper.

Our results for different combinations of detectors and descriptors evaluated on KTH are presented in Table 1. The best results are obtained for Harris3D + HOF (92.1%) and HOG/HOF (91.8%). These results are comparable to 91.8% reported in [START_REF] Laptev | Learning realistic human actions from movies[END_REF] for Harris3D + HOG/HOF. For Harris3D + HOG3D, we only reach 89.0%, about 2.5% lower than the original result in [START_REF] Kläser | A spatio-temporal descriptor based on 3Dgradients[END_REF]. This could be explained by the different strategy of codebook generation (random sampling) used in [START_REF] Kläser | A spatio-temporal descriptor based on 3Dgradients[END_REF]. For the Cuboid detector, the best result 90.0% is obtained with the HOG3D descriptor. The performance of Hessian and Dense detectors are below Harris3D and Cuboids. The low performance of dense sampling on KTH may be explained by the large number of features corresponding to the static background. When comparing performance of different descriptors, we note that HOG/HOF and HOF give best results in combination with Harris3D, Hessian and Dense features.

UCF sports dataset

The results for different combinations of detectors and descriptors evaluated on UCF sport actions are illustrated in Table 2. The best result 85.6% over different detectors is obtained by the dense sampling. This can be explained by the fact that dense features capture different types of motions. Furthermore, they also capture background which may provide useful context information. Scene context indeed may be helpful for sports actions which often involve specific equipment and scene types as illustrated in Figure 1. The second-best result

Hollywood2 dataset

Finally, evaluation results for Hollywood2 actions are presented in Table 3. As for the UCF dataset, the best result 47.4% is obtained for dense sampling while interest point detectors demonstrate similar and slightly lower performance. We assume dense sampling again benefits from a more complete description of motions and the rich context information. Among different descriptors, HOG/HOF performs best. Unlike in results for KTH actions, here the combination of HOF and HOG improves HOF with about 2 percent. The HOG3D descriptor performs similar to HOF. Shot boundary features. Since action samples in Hollywood2 are collected from movies, they contain many shot boundaries, which cause many artificial interest points. To investigate the influence of shot boundaries on recognition results, we compare in Table 4 the performance of the Harris3D detector with and without shot boundary features. Results for HOG/HOF and HOG demonstrate 0.5% and 2% improvement respectively when removing shot boundary features while the change in performance for other descriptors is minor. We conclude that shot boundary features do not influence our evaluation significantly.

Influence of subsampling. We also investigate the influence of reduced spatial resolution adopted in our Hollywood2 experiments. In Table 4 recognition results are reported for videos with full and half spatial resolution using the Harris3D detector. The performance is consistently and significantly increased for all tested descriptors for the case of full spatial resolution. Note that for full resolution, we obtain approximately 3 times more features per sequence than for half resolution.

Dense sampling parameters

Given the best results obtained with dense sampling, we further investigate the performance as a function of different minimal spatial sizes of dense descriptors (cf. present results for Hollywood2 and UCF videos with half spatial resolution. We observed no significant improvements for different temporal lengths, therefore we fixed the temporal length to 10 frames. The overlapping rate for dense patches is set to 50%. We can see that the performance increases with smaller spatial size, i.e., when we sample denser. However, the performance saturates in general at a spatial size of 24 × 24 for Hollywood2 and 18 × 18 for UCF.

Computational complexity

Here we compare the tested detectors by their speed and the number of detected interest points. The comparison was performed on a set of videos from Hollywood2 with spatial resolution of 360 × 288 pixels (half resolution) and about 8000 frames length in total. The run-time estimates were obtained on a Dell Precision T3400 Dual core PC with 2.66 GHz processors and 4GB RAM. Table 6 presents results for the three detectors and dense sampling in terms of average number of frames per second and average number of features per frame. Note that feature computation is included in the run time. Among the detectors, Cuboid extracts the densest features (44 features/frame) and it is the slowest one (0.9 frames/second). Hessian extracts the sparsest features (19 features/frame) and is consequently the most efficient (4.6 frames/second). As for the dense sampling, since there was no feature detection as such, the overall computational time was mainly spent on the feature description. Obviously, dense sampling extracts many more features than interest point detectors. Note that the time of descriptor quantization was not taken into account in this evaluation.

Conclusion

Among the main conclusions, we note that dense sampling consistently outperforms all tested interest point detectors in realistic video settings, but performs worse on the simple KTH dataset. This indicates both (a) the importance of using realistic experimental video data as well as (b) the limitations of current interest point detectors. Note, however, that dense sampling also produces a very large number of features (usually 15-20 times more than feature detectors). This is more difficult to handle than the relatively sparse number of interest points. We also note a rather similar performance of interest point detectors for each dataset. Across datasets, Harris 3D performs better on KTH dataset, while the Cuboid detector gives better results for UCF and Hollywood2 datasets. Among the tested descriptors, the combination of gradient based and optical flow based descriptors seems to be a good choice. The combination of dense sampling with the HOG/HOF descriptor provides best results for the most challenging Hollywood2 dataset. On the UCF dataset, the HOG3D descriptor performs best in combination with dense sampling. This motivates further investigations of optical flow based descriptors.

Figure 1 :

 1 Figure 1: Sample frames from video sequences of KTH (top), UCF Sports (middle), and Hollywood2 (bottom) human action datasets. sums v = (∑ d x , ∑ d y , ∑ d t ) of uniformly sampled responses of the Haar-wavelets d x , d y , d t along the three axes. We use the executables from the authors' website 3 with the default parameters defined in the executable.

Table 1 :

 1 [START_REF] Blank | Actions as space-time shapes[END_REF]. The ESURF and Cuboid descriptors are not Average accuracy for various detector/descriptor combinations on the KTH dataset.

		HOG3D HOG/HOF	HOG	HOF	Cuboids ESURF
	Harris3D	89.0%	91.8%	80.9% 92.1%	-	-
	Cuboids	90.0%	88.7%	82.3% 88.2%	89.1%	-
	Hessian	84.6%	88.7%	77.7% 88.6%	-	81.4%
	Dense	85.3%	86.1%	79.0% 88.0%	-	-
		HOG3D HOG/HOF	HOG	HOF	Cuboids ESURF
	Harris3D	79.7%	78.1%	71.4% 75.4%	-	-
	Cuboids	82.9%	77.7%	72.7% 76.7%	76.6%	-
	Hessian	79.0%	79.3%	66.0% 75.3%	-	77.3%
	Dense	85.6%	81.6%	77.4% 82.6%	-	-

Table 2 :

 2 Average accuracy for various detector/descriptor combinations on the UCF dataset.

Table 3 :

 3 Mean AP for various detector/descriptor combinations on the Hollywood2 dataset.

		HOG3D HOG/HOF	HOG	HOF	Cuboids ESURF
	Harris3D	43.7%	45.2%		32.8% 43.3%	-	-
	Cuboids	45.7%	46.2%		39.4% 42.9%	45.0%	-
	Hessian	41.3%	46.0%		36.2% 43.0%	-	38.2%
	Dense	45.3%	47.4%		39.4% 45.5%	-	-
			HOG3D HOG/HOF	HOG	HOF
	reference			43.7%	45.2%	32.8% 43.3%
	w/o shot boundary features	43.6%	45.7%	35.3% 43.4%
	full resolution videos		45.8%	47.6%	39.7% 43.9%

Table 4 :

 4 Comparison of the Harris3D dector on (top) videos with half spatial resolution, (middle) with removed shot boundary features, and (bottom) on the full resolution videos. 82.9% is obtained for the Cuboid detector. Also above 80% are dense points in combination with HOG/HOF and HOF. Harris3D and Hessian detectors perform similar at the level of 80%. Among different descriptors, HOG3D provides best results for all detectors except Hessian. HOG/HOF gives second-best result for UCF. The authors of the original paper[START_REF] Rodriguez | Action mach: A spatio-temporal maximum average correlation height filter for action recognition[END_REF] report 69.2% for UCF. Their result, however, does not correspond to the version of UCF dataset available on-line 6 used in our evaluation.

Table 5

 5

	Spatial		Hollywood2		UCF
	Size	HOG3D HOG/HOF HOG HOF	HOG3D HOG/HOF HOG HOF
	18 × 18 24 × 24 36 × 36 48 × 48 72 × 72	45.3% 45.1% 44.8% 42.8% 39.7%	47.4% 47.7% 39.4% 45.8% 39.4% 45.5% 47.3% 36.8% 45.6% 46.5% 35.8% 45.5% 45.2% 32.2% 43.0%	85.6% 82.0% 78.6% 78.8% 77.8%	81.6% 77.4% 82.6% 81.4% 76.8% 84.0% 79.1% 76.5% 82.4% 78.6% 73.9% 79.0% 78.8% 69.6% 78.4%

). As before, further spatial scales are sampled with a scale factor of √ 2. As in Sections 4.2 and 4.3, we

Table 5 :

 5 Average accuracy for dense sampling with varying minimal spatial sizes on the Hollywood2 and UCF sports dataset.

		Harris3D + Hessian + Cuboid Detector + Dense +	Dense +
		HOG/HOF	ESURF	Descriptor	HOG3D HOG/HOF
	Frames/second	1.6	4.6	0.9	0.8	1.2
	Features/frame	31	19	44	643	643

Table 6 :

 6 Average speed and average number of generated features for different methods.

http://www.irisa.fr/vista/Equipe/People/Laptev/download.html#stip

http://vision.ucsd.edu/~pdollar/toolbox/doc/index.html

http://homes.psat.kuleuven.be/~gwillems/research/Hes-STIP/

http://lear.inrialpes.fr/software

Available at http://www.nada.kth.se/cvap/actions/

Available at http://www.cs.ucf.edu/vision/public_html/

Available at http://www.irisa.fr/vista/Equipe/People/Laptev/download.html

Acknowledgements. This work was partially funded by the European research project CLASS, the MSR/INRIA joint project, and the QUAERO project.


D Doyen

A Ern

S Piperno

Time-integration schemes for the finite element dynamic Signorini problem

Keywords: elastodynamics, frictionless unilateral contact, time-integration schemes, finite elements, modified mass method

The discretization of the dynamic Signorini problem with finite elements in space and a time-stepping scheme is not straightforward. Consequently a large variety of methods for this problem have been designed over the last two decades. Up to date, no systematic comparison of such methods has been performed. The aim of the present work is to classify and compare them. For each method, we discuss the presence of spurious oscillations and the energy conservation. For explicit approaches, the stability condition on the time step is also discussed. Numerical simulations on two 1D benchmark problems with analytical solutions illustrate the properties of the different methods. Most of the discretizations considered herein can be found in the literature, but the semi-explicit modified mass method is new and features, in our opinion, several attractive properties.

INTRODUCTION

The design of robust and efficient numerical methods for dynamic contact problems has motivated a large amount of work over the last two decades and remains a challenging issue. Here, we focus on the dynamic Signorini problem, which models the infinitesimal deformations of a solid body that can come into contact with a rigid obstacle. This problem is the simplest dynamic contact problem, but also the first step toward more complex situations, such as multi-body problems, large deformation problems, contact with friction, etc... For an overview of the different contact problems, we refer to [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF][START_REF] Laursen | Computational contact and impact mechanics[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF]. In structural dynamics, the usual space-time discretization combines finite elements in space and a time-stepping scheme. In this framework, the discretization of the dynamic Signorini problem involves mainly three choices: (i) the finite element space; (ii) the enforcement of the contact condition; (iii) the time-stepping scheme. The combination of these three ingredients is not straightforward. For instance, it is well-known that the combination of an exact enforcement of the contact condition and an implicit Newmark scheme yields spurious oscillations as well as poor energy conservation and long-time behavior. Moreover, the combination of an exact enforcement and an explicit scheme is not straightforward, whereas the use of a penalty contact condition tightens the stability condition of explicit schemes. Consequently, various alternative discretizations have been designed for the dynamic Signorini problem. However, up to date, no systematic comparison of such methods has been performed. The aim of the present work is to classify and compare the main discretizations for the dynamic Signorini problem. Most of the discretizations considered herein can be found in the literature. The semi-explicit modified mass method (Discretization 7.2 below) is new and features, in our opinion, several attractive properties. We classify the different discretizations into four groups. The first three groups correspond to different ways of enforcing the contact condition: exact enforcement [START_REF] Carpenter | Lagrange constraints for transcient finite element surface contact[END_REF][START_REF] Deuflhard | A contact-stabilized Newmark method for dynamical contact problems[END_REF][START_REF] Kane | Finite element analysis of nonsmooth contact[END_REF][START_REF] Laursen | Improved implicit integrators for transient impact problems-geometric admissibility within the conserving framework[END_REF][START_REF] Paoli | A numerical scheme for impact problems. I. The one-dimensional case[END_REF][START_REF] Paoli | A numerical scheme for impact problems. II. The multidimensional case[END_REF][START_REF] Taylor | On a finite element method for dynamic contact/impact problems[END_REF][START_REF] Vola | Consistent time discretization for dynamical contact problems and complementarity techniques[END_REF], enforcement with penalty [START_REF] Armero | Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems[END_REF][START_REF] Belytschko | Contact-impact by the pinball algorithm with penalty and lagrangian methods[END_REF][START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF] and enforcement with contact condition in velocity [START_REF] Ayyad | A frictionless viscoelastodynamic contact problem with energy consistent properties: Numerical analysis and computational aspects[END_REF][START_REF] Belytschko | Contact-impact by the pinball algorithm with penalty and lagrangian methods[END_REF][START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF]. The fourth approach is based on a modification of the mass matrix [START_REF] Hager | A stable energy-conserving approach for frictional contact problems based on quadrature formulas[END_REF][START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF]; it can be seen as an alternative choice of the finite element space. These four classes yield different semidiscrete problems in space which in turn can be discretized in time using various time-stepping schemes, either implicit or (semi-)explicit. In our opinion, a reliable discretization should not trigger spurious oscillations during the contact phases and achieve a tight energy conservation. Energy conservation is important to guarantee a good behavior in long time. By energy conservation we mean that the variation of the energy is equal to the work of the external forces (the contact forces should not work). For each scheme, the presence of spurious oscillations and the energy conservation are examined. In addition, for explicit approaches, the stability condition on the time step is also discussed. Furthermore dynamic contact problems yield shock waves, and spurious oscillations appear near the shock in the numerical solutions, owing to the so-called Gibbs phenomenon. This issue being important but not specific to dynamic contact problems, it will be briefly addressed in the appendix. To illustrate the behavior of the different discretizations, numerical simulations on two 1D benchmark problems with analytical solutions have been performed. The first benchmark, the impact of an elastic bar, is well-known in the literature and allows to detect spurious oscillations. The second one, the bounces of an elastic bar, is geared towards the energy conservation and long-time behavior of the discretizations. It is new to our knowledge. Although both benchmarks are 1D, we believe they are representative of the various difficulties encountered in higher space dimension. Furthermore, the mathematical analysis of the different methods is beyond the scope of this article, but we mention, whenever they exist, the theoretical results (well-posedness of the discrete problems and convergence of the discrete solutions). The article is organized as follows. We formulate the dynamic Signorini problem in the continuous setting (Section 2.1) and we introduce the main ingredients for its approximation (Sections 2.2 and 2.3). We present the two benchmark problems with their analytical solutions (Section 3). We describe the four classes of discretizations together with numerical results: exact enforcement of the contact condition (Section 4), enforcement with penalty contact condition (Section 5), enforcement with contact condition in velocity (Section 6) and modification of the mass matrix (Section 7). The Appendix is devoted to the Gibbs phenomenon.

THE DYNAMIC SIGNORINI PROBLEM

Governing equations

We consider the infinitesimal deformations of a body occupying a reference domain Ω ⊂ R d (d = 1, d = 2 or d = 3) during a time interval [0, T ]. The tensor of elasticity is denoted by A and the mass density is denoted by ρ. An external load f is applied to the body. Let u : (0, T ) × Ω → R d , ǫ(u) : (0, T ) × Ω → R d,d and σ(u) : (0, T ) × Ω → R d,d be the displacement field, the linearized strain tensor and the stress tensor, respectively. Denoting time-derivatives by dots, the equilibrium equation is

ρü -div σ = f, σ = A : ǫ, ǫ = 1 2 (∇u + T ∇u) in Ω × (0, T ). ( 1
)
The boundary ∂Ω is partitioned into three disjoint open subsets Γ D , Γ N and Γ c . Dirichlet and Neumann conditions are prescribed on Γ D and Γ N , respectively,

u = u D on Γ D × (0, T ), σ • n = f N on Γ N × (0, T ), (2)
where n denote the outward unit normal to Ω. We set u n := u |∂Ω • n and σ n := n • σ |∂Ω • n, the normal displacement and the normal stress on ∂Ω, respectively. On Γ c , a unilateral contact condition, also called Signorini condition, is imposed,

u n ≤ 0, σ n (u) ≤ 0, σ n (u)u n = 0 on Γ c × (0, T ). (3)
At the initial time, the displacement and velocity fields are prescribed,

u(0) = u 0 , u(0) = v 0 in Ω. (4)
Problem ( 1)-( 4) is an evolution partial differential equation under unilateral constraints. Here, the equation is second-order in time and the constraint holds on the displacement; this is not the most favorable case. The existence and uniqueness of a solution has only been proven in 1D, when the contact boundary is reduced to a point [START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF][START_REF] Eck | Unilateral contact problems[END_REF]. In 1D, it has also been proven that the variation of energy is equal to the work of the external forces; the contact force does not work [START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF][START_REF] Eck | Unilateral contact problems[END_REF]. In higher dimension, the existence of a solution is proven in the case of a viscoelastic material [START_REF] Eck | Unilateral contact problems[END_REF].

Basic time-integration schemes in linear elastodynamics

In this section, we briefly recall some basic facts about time-integration schemes in linear elastodynamics; most of this material can be found in [START_REF] Hughes | The finite element method[END_REF]. Firstly, we discretize the problem in space with a finite element method. The number of degrees of freedom is denoted by N d . Let K, M , and F (t) be the stiffness matrix, the mass matrix, and the column vector of the external forces, respectively. The space semi-discrete problem consists in seeking u :

[0, T ] → R N d such that, for all t ∈ [0, T ], M ü(t) + Ku(t) = F (t), (5)
with the initial conditions u(0) = u 0 and u(0) = v 0 . For solving such a system of ODEs (ordinary differential equations), linear one-step schemes are the most frequently used. The interval [0, T ] is divided into equal subintervals of length ∆t. We set t n = n∆t and denote by u n , un , and ün the approximations of u(t n ), u(t n ), and ü(t n ), respectively. We define the convex combination n+α := (1-α) n +α n+1 , where stands for u, u, ü or t, and α ∈ [0, 1]. We also set F n+α := F (t n+α ). At time t n , the energy of the system is

E n := 1 2 T un M un + 1 2 T u n Ku n .
It is also convenient to define the quadratic form

E n β,γ := E n + ∆t 2 2
β -1 2 γ T ün M ün . Now we can formulate some of the most common time-stepping schemes in linear elastodynamics.

Discretization 2.1 (HHT-Newmark) Seek u n+1 , un+1 , ün+1 ∈ R N d such that M ün+1 + Ku n+1+α = F n+1+α , (6)
u n+1 = u n + ∆t un + ∆t 2 2 ün+2β , (7)
un+1 = un + ∆t ün+γ , (8)
where α, β, γ are parameters.

Discretization 2.2 (Midpoint) Seek u n+1 , un+1 , ün+ 1 2 ∈ R N d such that M ün+ 1 2 + Ku n+ 1 2 = F n+ 1 2 , (9)
u n+1 = u n + ∆t un+ 1 2 , ( 10
) un+1 = un + ∆t ün+ 1 2 . ( 11
)
Discretization 2.3 (Central differences) Seek u n+1 ∈ R N d such that M u n+1 -2u n + u n-1 ∆t 2 + Ku n = F n . (12)
When α = 0, the schemes ( 6)-( 8) are called Newmark schemes. For these schemes, the energy balance reads [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF],

E n+1 β,γ -E n β,γ = T 1 2 (F n+1 + F n ) + γ - 1 2 (F n+1 -F n ) (u n+1 -u n ) -γ - 1 2 T (u n+1 -u n )K(u n+1 -u n ) + β - 1 2 γ T (ü n+1 -ün )M (ü n+1 -ün ) . ( 13
)
The particular choice β = 1/4, γ = 1/2 yields an implicit, unconditionally stable, and secondorder scheme. It is energy-conserving in the sense that

E n+1 -E n = T F n+1 (u n+1 -u n ). ( 14
)
When α ∈ [-1/3, 0], β = 1/4(1-α) 2 , γ = 1/2-α, the schemes ( 6)-( 8) are called HHT schemes or α-methods. Such schemes are implicit, unconditionally stable, second-order and dissipative in the high frequencies. The amount of dissipation is controlled by the parameter α. The midpoint scheme is implicit, unconditionally stable, and second-order. It is energy-conserving, in the sense that

E n+1 -E n = T F n+ 1 2 (u n+1 -u n ). ( 15
)
The central difference scheme is explicit (provided the mass matrix is lumped), conditionally stable and second-order. Here it is written as a two-step linear scheme involving only the displacement but it can be formulated as a one-step scheme. Actually, it is a Newmark scheme with parameters β = 0, γ = 1/2; the associated velocity and acceleration are then un = 1 2∆t (u n+1 -u n-1 ) and ün = 1 ∆t 2 (u n+1 -2u n + u n-1 ). The central difference scheme conserves a quadratic form which is nearly the energy,

E n+1 0, 1 2 -E n 0, 1 2 = T F n+1 + F n 2 (u n+1 -u n ). ( 16
)
There exist also second-order explicit schemes with high-frequency dissipation, such as the Chung-Hulbert schemes [START_REF] Hulbert | Explicit time integration algorithms for structural dynamics with optimal numerical dissipation[END_REF].

Enforcing the contact condition

The enforcement of a contact condition in a finite element setting has been widely studied in the case of the static Signorini problem [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]. We assume that the mesh is compatible with the partition of the boundary. Let N c be the number of nodes lying on the contact boundary. We define the linear normal trace operator on Γ c , g : v -→ -v |Γ • n and the associated matrix G.

Note that the dimension of G is N c × N d . We denote by {G i } 1≤i≤Nc the rows of the matrix G. Thus, G i u yields the value of the normal displacement at the ith node of the contact boundary. With an exact enforcement, the static Signorini problem consists in seeking a displacement u ∈ R N d and a contact pressure r ∈ R Nc such that

Ku = F + T Gr, (17)
Gu ≥ 0, r ≥ 0, T rGu = 0. ( 18
)
Here the problem is formulated as a complementarity problem. Other formalisms can be found in the literature, e.g., variational inequality, Lagrangian formulation and formulation with subderivatives. If the matrix K is positive definite, problem ( 17)-( 18) has a unique solution. For solving this problem, a large variety of methods has been developed [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF]: Uzawa algorithms, active set methods, semi-smooth Newton methods, Lemke algorithm, etc... Penalty formulations are another classical way of dealing with constrained problems. We have to define a penalty function R ǫ : R Nc → R Nc . For instance we can choose R ǫ (v) = 1 ǫ (v) -, where (v) -denotes the negative part of v. The penalized static Signorini problem consists now in seeking u ∈ R N d such that

Ku = F + T GR ǫ (Gu). (19)
A third way of enforcing the contact condition, specific to the dynamic problem, is to replace the Signorini condition by an approximation involving the velocity instead of the displacement [START_REF] Eck | Unilateral contact problems[END_REF]. Assume that u n = 0 at a certain time t c . Then, on a short time interval afterwards, u n ≈ (t -t c ) un . This motivates the following contact condition in velocity,

un ≤ 0, σ n (u) ≤ 0, σ n (u) un = 0, on Γ c . ( 20
)
It must be stated that this condition is not rigorously equivalent to the Signorini condition.

BENCHMARK PROBLEMS

To compare the different methods, we test them on two 1D problems. Both problems can be formulated in the same setting. We consider an elastic bar dropped against a rigid ground. The bar is dropped, undeformed, from a height h 0 , with an initial velocity -v 0 , under a gravity g 0 ≥ 0. The length of the bar is denoted by L, the Young modulus by E and the density by ρ.

Let c 0 := E ρ denote the wave speed. The reference domain is Ω = [0, L]. The displacement u of the bar and the contact pressure r satisfy the following equations,

ρü -E ∂ 2 u ∂x 2 = -ρg 0 + r, in Ω × (0, T ), (21)
u(0, t) ≥ 0, r(t) ≥ 0, r(t)u(0, t) = 0 on (0, T ), ( 22) ∂u ∂x (L, t) = 0 on (0, T ), ( 23)

u(•, 0) = h 0 , u(•, 0) = -v 0 . (24)
Problem ( 21)-( 24) has a unique solution and the variation of the energy is equal to the work of the gravity force [START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF],

d dt 1 2 Ω ρ u2 + 1 2 Ω E ∂u ∂x 2 = Ω -g 0 u, ∀t ∈ [0, T ]. (25)
In the first problem, v 0 > 0 and g 0 = 0. This benchmark has been widely used in the literature. It enables us to examine the possible spurious oscillations triggered by the numerical schemes.

In the second problem, v 0 = 0 and g 0 > 0, so that the bar can make several bounces. With a suitable choice of parameters, the motion of the bar is periodic in time and we can calculate the exact solution. It enables us to examine the energy conservation and the long-time behavior of the different numerical methods. This benchmark is, to our knowledge, new.

Impact of an elastic bar

Let us describe the solution of the benchmark problem (Figure 1). Before the impact, the bar is undeformed and has a uniform velocity -v 0 . The bar reaches the ground at time t 1 = h0 v0 . After the impact, the bar stays in contact with the ground. A shock wave travels from the bottom of the bar to the top. Above the shock wave, the velocity is -v 0 ; below, the velocity is zero. Then, the shock wave travels from the top to the bottom. Above the shock wave, the velocity is v 0 ; below, the velocity is still zero. As soon as the wave reaches the bottom, the bar takes off, undeformed, with a uniform velocity v 0 . The speed of the shock wave is c 0 . Thus, the wave takes a time τ w := L c0 to cross the bar, and the bar takes off at time t 2 = t 1 + 2τ w . The analytical solution can be easily expressed using travelling functions. Defining the auxiliary function

H v (x, t) = -v min(x/c 0 , τ w -|t -τ w |), the exact solution is u(x, t) =      h 0 -v 0 t if t ≤ t 1 H v0 (x, t -t 1 ) if t 1 < t ≤ t 2 v 0 (t -t 2 ) if t 2 < t (26)
In particular the displacement at the bottom of the bar and the contact pressure are:

u(0, t) =      h 0 -v 0 t if t ≤ t 1 0 if t 1 < t ≤ t 2 v 0 (t -t 2 ) if t 2 < t r(t) =      0 if t ≤ t 1 Ev0 c0 if t 1 < t ≤ t 2 0 if t 2 < t (27)
These two quantities are represented in Figure 2 (with the parameters chosen in Section 3.3

t = 0 t 1 -v 0 -v 0 v 0 t 2 v 0 -v 0 0 0 0 Figure 1.
Impact of an elastic bar.

Bounces of an elastic bar

In this benchmark problem (Figure 3), the bar is dropped, undeformed, with a zero velocity.

It takes a time τ f := 2h0 g0 to reach the ground. At the impact, at time t 1 := τ f , the bar is undeformed and has uniform velocity -v f , where v f := √ 2h 0 g 0 . After the impact, as in the previous benchmark, the bar stays in contact with the ground during a time 2τ w . During this contact phase, the response of the bar is the superposition of a shock wave due to velocity at the impact and a vibration due to the gravity. When the bar takes off, at time t

2 := t 1 + 2τ w , it has a uniform velocity v f but it is compressed (u(x, t 2 ) = ũ(x) := g0 c0 (x 2 -2Lx)).
Consequently, during the flight phase, the response of the bar, is the superposition of a rigid parabolic motion (due to the gravity and the velocity) and a vibration (due to the initial deformation). If we choose proper parameters (for instance, τ f = 3τ w ), we can ensure that the bar reaches the ground with uniform velocity -v f and with displacement field ũ. By doing so, the second impact occurs at time t 3 := t 2 + 2τ f = t 2 + 6τ w . When the bar takes off again, at time t 4 := t 3 + 2τ w , it is undeformed and has a uniform velocity v f . The next flight phase is a rigid parabolic movement. Then, this sequence of two contact phases and two flight phases repeats periodically. To compute the analytical solution, we use a decomposition on the eigenmodes in addition to the travelling functions. We set t 4k+1 = 3τ w + 16kτ w , t 4k+2 = t 4k+1 + 2τ t 4k+2 = t 4k+1 + 8τ w and t 4k+4 = t 4k+1 + 10τ w . We define also the auxiliary functions

t = 0 t 1 t 2 t 3 t 4 -v f v f v f -v f v 0
P (x, t) = h 0 - 1 2 g 0 (t -τ f ) 2 , (28)
S 1 (x, t) = ∞ n=1 a n (1 -cos(c 0 ν n t))sin(ν n x), (29)
S 2 (x, t) = - 2g 0 L 2 3c 2 0 + ∞ n=1 b n cos(c 0 λ n t)cos(λ n x), (30)
where

a n = -2g0 c 2 0 Lν 3 n , ν n = (n -1 2 ) π L , b n = 4g0 g0λ 2 n , λ n = n π L .
The function S 1 corresponds to the vibration of a bar, clamped at its bottom, initially at rest, under a gravity g 0 . The function S 2 corresponds to the vibration of a bar, free at its two extremities, with the initial displacement ũ, a zero initial velocity and without external force. The computation of series S 1 and S 2 is standard; see [START_REF] Courant | Methods of mathematical physics[END_REF] for instance. The exact solution is

u(x, t) =                P (x, t + τ f ) if t ≤ t 1 H v f (x, t -t 4k+1 ) + S 1 (x, t -t 4k+1 ) if t 4k+1 < t ≤ t 4k+2 P (x, t -t 4k+2 ) + S 2 (x, t -t 4k+2 ) if t 4k+2 < t ≤ t 4k+3 H v f (x, t 4k+4 -t) + S 1 (x, t 4k+4 -t) if t 4k+3 < t ≤ t 4k+4 P (x, t -t 4k+4 ) if t 4k+4 < t ≤ t 4(k+1)+1 (31)
The displacement at the bottom of the bar is represented in Figure 4 (with the parameters chosen in Section 3.3).

Numerical simulations

The parameters used in the numerical simulations are E = 900, ρ = 1, L = 10, h 0 = 5.

In the first benchmark v 0 = 10, g 0 = 0; in the second benchmark v 0 = 0, g 0 = 10. The bar is discretized with a uniform meshsize ∆x and linear finite elements are used. We define ν c := c 0 ∆t ∆x the Courant number, which is the relevant ratio to link the mesh size and the time step. In particular, the central difference scheme with the lumped mass matrix is stable in the linear case under the condition ν c ≤ 1. Let us set some notation which will be used in the description of numerical results. If v is a displacement vector, we denote by v 0 and (Kv) 0 the displacement and the stress at the bottom node of the bar, respectively. In the benchmark problems, the load vector is time-independent; we denote it by F .

DISCRETIZATIONS WITH EXACT ENFORCEMENT OF THE CONTACT CONDITION

In this section we combine standard finite elements in space and an exact enforcement of the contact condition at each node of the contact boundary. This leads to the semi-discrete problem,

M ü(t) + Ku(t) = F (t) + T Gr(t), (32)
Gu(t) ≥ 0, r(t) ≥ 0, T r(t)Gu(t) = 0. (33)
Problem (32)-(33) is a system of differential equations under unilateral constraints. The same kind of formulation arises in rigid-body dynamics with impact [START_REF] Brogliato | Nonsmooth impact mechanics[END_REF][START_REF] Stewart | Rigid-body dynamics with friction and impact[END_REF], so the mathematical results and the numerical methods developed in this framework can in general be applied to our problem. Mathematically, problem (32)-(33) turns out to be delicate. Firstly, the functional framework is not obvious. Due to the unilateral constraints, the velocity can be discontinuous and there is in general no strong solution (i.e. twice differentiable in time) to this problem. We then have to define a concept of weak solution: the displacement u is continuous in time, the velocity u is a function with bounded variation in time, the acceleration ü and the contact pressure r are measures (they contain impulses). Secondly, this weak solution is in general not unique. Consider the simple example of a point mass impacting a rigid foundation. Before the impact, the motion of the point mass is uniquely determined. After the impact, an infinity of velocities and trajectories are admissible (Figure 5). To recover uniqueness, an additional condition, specifying the velocity after an impact, is needed. If we denote by v -the normal velocity before the impact and by v + the normal velocity after the impact, we can prescribe

v + = -ev -, ( 34
)
where e is a non-negative parameter. This impact law is commonly used in rigid-body dynamics. In our case, it seems reasonable to take e = 0. Indeed, in the continuous dynamic Signorini problem, the unilateral constraint holds on the boundary and the boundary does not bounce after an impact. In the discrete setting, the unilateral constraint holds on the nodes 000000 000000 000000 111111 111111 111111 0000000 0000000 0000000 1111111 1111111 1111111 0000000 0000000 0000000 1111111 1111111 1111111 of the contact boundary. If we want these nodes to have the same behavior as the boundary in the continuous setting, we are led to choose e = 0. We can now formulate the semi-discrete problem in a rigorous way.

v + ? v -
Problem 4.1. Seek a displacement u : [0, T ] → R N d and a contact pressure r : [0, T ] → R Nc such that M ü + Ku = F + T Gr, ( 35
)
Gu ≥ 0, r ≥ 0, T rGu = 0, ( 36
) T r i (t)G i u(t + ) = 0 if G i u(t) = 0, ( 37
)
with the initial conditions u(0) = u 0 and u(0) = v 0 .

Most of the mathematical terms in equations ( 35)-(37) must be understood in the sense of measures. In particular T rGu and T r i (t)G i u(t + ) should be defined with suitable duality products. For more details, we refer to [START_REF] Brogliato | Nonsmooth impact mechanics[END_REF][START_REF] Stewart | Rigid-body dynamics with friction and impact[END_REF].

Remark 4.1. The impact law is a consequence of the discretization in space. Indeed, the continuous problem does not need an impact law to have a unique solution. This fact is proven in 1D [START_REF] Eck | Unilateral contact problems[END_REF][START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF]; in higher dimension, the uniqueness is still an open problem, but the difficulty does not seem to come from the absence of an impact law.

Remark 4.2. The semi-discrete solution does not conserve the energy (there is a loss of energy at each impact of a node). This is another difference with the continuous solution.

Remark 4.3. The impact law is different from the concept of persistency condition sometimes encountered in the literature [START_REF] Armero | Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems[END_REF][START_REF] Laursen | Computational contact and impact mechanics[END_REF][START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF][START_REF] Laursen | Improved implicit integrators for transient impact problems-geometric admissibility within the conserving framework[END_REF]. The persistency condition is defined in the continuous setting and in the fully discrete setting. It requires that the contact force does not work. Note that in the continuous problem, the persistency condition seems to be a consequence of the Signorini condition (it is at least proven in 1D).

Implicit schemes

To motivate the discussion, let us begin with an ill-founded discretization. We choose a Newmark scheme (trapezoidal rule) for the elastodynamic part and we enforce the contact condition (36) at a certain time, say t n+1 . We pay no attention to the impact law (37). This choice corresponds to Discretization 4.1 with α = 0, β = 1/4, γ = 1/2.

M ün+1 + Ku n+1+α = F n+1+α + T Gr n+1 , (38)
Gu n+1 ≥ 0, r n+1 ≥ 0, T r n+1 Gu n+1 = 0, (39)
u n+1 = u n + ∆t un + ∆t 2 2 ün+2β , (40)
un+1 = un + ∆t ün+γ .

At each time step, Problem (38)-( 41) is equivalent to a linear complementarity problem and is well-posed. In contrast to the static case, the matrix K does not need to be definite for the problem to be well-posed (Dirichlet boundary conditions are not needed). When this scheme is tested on the first benchmark, we observe large spurious oscillations on the contact pressure and small spurious oscillations on the displacement during the contact phase (Figure 6). On the second benchmark, we observe a poor behavior in long-time and a poor energy conservation (Figure 7). Let us try to explain what happens exactly. Suppose there is contact at the ith node of the contact boundary at time t n+1 (i.e. G i u n+1 = 0), then

G i un+1 = - 1 2∆t G i u n + 1 - β γ G i un + ∆t 2β -γ 2β G i ün , (42)
G i ün+1 = - 1 β∆t 2 G i u n - 1 β∆t G i un - 1 -2β 2β G i ün . (43)
As we can see, the impact law is not satisfied since we would expect that after an impact, G i un = G i ün = 0. During a contact phase following an impact, the velocity and the acceleration oscillate. For the acceleration, the magnitude of the oscillations is vi ∆t , where v i is the velocity before the impact. These oscillations of the acceleration trigger oscillations of magnitude mivi ∆t on the contact pressure, where m i is the mass associated with the node i (Figure 6). Moreover, the energy balance takes the form,

E n+1 -E n = T r n+1 + r n 2 (Gu n+1 -Gu n ) + T F n+1 + F n 2 (u n+1 -u n ), (44)
so that the contact force works when a node changes status. When a node comes into contact (G i u n > 0, r n = 0, G i u n+1 = 0, r n+1 > 0), the work is negative; when a node is released (G i u n = 0, r n > 0, G i u n+1 > 0, r n+1 = 0), the work is positive. As the contact pressure is polluted by large oscillations, this strongly perturbs the rest of the structure (Figure 7). The poor behavior of the Newmark scheme can be summarized by the following diagram: be avoided. Several options can be considered to design better algorithms. The first option consists in using dissipative schemes, such as HHT schemes (Discretization 4.1). The spurious oscillations are then damped (Figure 8), but at the expense of poor energy conservation (Figure 9). First-order schemes like θ-schemes, which are implicit, unconditionally stable, dissipative schemes, yield the same kind of results (Discretization 4.2). Discretization 4.2 (θ-schemes [START_REF] Vola | Consistent time discretization for dynamical contact problems and complementarity techniques[END_REF]) Seek u n+1 , un+1 ∈ R N d , and r n+1 ∈ R Nc such that preferred to second-order schemes for dynamic contact problems, due to the non-smoothness of the solution. We must distinguish two issues: the treatment of the contact condition and the treatment of the shock waves induced by the contact. As discussed previously, a proper treatment of the contact condition is not related to the order of the scheme. As for the shock waves, they require a scheme with dissipation and there exist second-order schemes with dissipation, such as the HHT or Chung-Hulbert schemes. We refer to the Appendix for further discussion.

M ün+ 1 2 + Ku n+θ = F n+θ + T Gr n+1 , (45)
Gu n+1 ≥ 0, r n+1 ≥ 0, T r n+1 Gu n+1 = 0, ( 46
)
u n+1 = u n + ∆t un+θ , (47)
A second option consists in finding a scheme which satisfies the impact law or, more precisely, a scheme which enforces the acceleration to be zero during the contact phases. No implicit Newmark scheme achieves this. An extra-step is needed to enforce the impact law (Discretization 4.3). As illustrated in Figure 10, the large spurious oscillations have disappeared. However, this stabilization takes effect only one step after the impact, which explains the peak in the contact pressure just after the impact. Stabilizations procedures which enforce the impact law as soon as the impact occurs have been proposed in [START_REF] Taylor | On a finite element method for dynamic contact/impact problems[END_REF] and [START_REF] Deuflhard | A contact-stabilized Newmark method for dynamical contact problems[END_REF]. This second problem will be presented at the end of the section (Discretization 4.6). As a third option, it is possible to prevent the oscillations of the acceleration from transferring to the contact pressure by removing the mass at the contact boundary. This approach will be developed in Section 7. The fourth option consists in finding a time discretization where the contact force does not work or is at least dissipative. For instance, the midpoint scheme with an enforcement of the contact condition at time t n+1 achieves the following energy balance,

Discretization 4.3 (naive stabilized Newmark

) 1. Seek u n+1 ∈ R N d , un+1 ∈ R N d , ün+1 ∈ R N d ,
E n+1 -E n = T r n+1 (Gu n+1 -Gu n ) + T F n+ 1 2 (u n+1 -u n ). ( 50
)
It is easy to check that the work of the contact force is always non-positive. As illustrated in Figure 11, the contact pressure still oscillates but the stress is practically free of oscillations. However, the losses of energy are quite important and deteriorate the long-time behavior (Figure 12).

Discretization 4.4 (Midpoint-implicit contact) Seek u n+1 , un+1 ∈ R N d , and r n+1 ∈ R Nc such that M ün+ 1 2 + Ku n+ 1 2 = F n+ 1 2 + T Gr n+1 , ( 51
)
Gu n+1 ≥ 0, r n+1 ≥ 0, T r n+1 Gu n+1 = 0, ( 52
)
u n+1 = u n + ∆t un+ 1 2 , ( 53
) un+1 = un + ∆t ün+ 1 2 . ( 54
)
Another scheme with dissipative contact has been proposed in [START_REF] Kane | Finite element analysis of nonsmooth contact[END_REF]. The Newmark scheme with parameters β = 1/2 and γ = 1 and with an enforcement of the contact condition at time t n+1 yields the following energy balance

E n+1 -E n = T r n+1 (Gu n+1 -Gu n ) - 1 2 T (u n+1 -u n )K(u n+1 -u n ) + T F n+1 (u n+1 -u n ). ( 55
)
The work of the contact force is always non-positive, but there is a strong bulk dissipation. To remove this dissipation, one can, as proposed in [START_REF] Kane | Finite element analysis of nonsmooth contact[END_REF], discretize the acceleration coming from the contact forces with the dissipative parameters (β = 1/2 and γ = 1) and the acceleration coming from the elastic forces with a trapezoidal rule (β = 1/4 and γ = 1/2). This yields Discretization 4.5. With such a discretization, the energy balance is

E n+1 -E n = T r n+1 (Gu n+1 -Gu n ) + T F n+1 (u n+1 -u n ). ( 56
)
The numerical results are similar to those obtained with Discretization 4.4.

Discretization 4.5 (Newmark with dissipative contact [START_REF] Kane | Finite element analysis of nonsmooth contact[END_REF])

Seek u n+1 , un+1 , ün+1 int , ün+1 con ∈ R N d , and r n+1 ∈ R Nc such that M ün+1 + Ku n+1 = F n+1 + T Gr n+1 , ( 57
)
Gu n+1 ≥ 0, r n+1 ≥ 0, T r n+1 Gu n+1 = 0, ( 58
)
u n+1 = u n + ∆t un + ∆t 2 2 ün+2β int + ∆t 2 2 ün+1 con , (59) un+1
= un + ∆tü n+γ int + ∆tü n+1 con , (60)
where ün+1 = ün+1 int + ün+1 con and M ün+1 con = T Gr n+1 . To compensate the losses of energy in schemes with dissipative contact, the so-called velocityupdate method has been introduced in [START_REF] Laursen | Improved implicit integrators for transient impact problems-geometric admissibility within the conserving framework[END_REF]. Applied to Discretization 4.4, this procedure does not improve significantly the long-time behavior on our benchmark. In [START_REF] Deuflhard | A contact-stabilized Newmark method for dynamical contact problems[END_REF], the authors add to Discretization 4.5 a stabilization procedure (Discretization 4.6). Discretization 4.6 (stabilized Newmark [START_REF] Deuflhard | A contact-stabilized Newmark method for dynamical contact problems[END_REF])

1. Seek u n+1 pred ∈ R N d and λ n+1 ∈ R Nc such that M u n+1 pred = M u n + ∆tM un , (61)
Gu n+1 pred ≥ 0, λ n+1 ≥ 0, T λ n+1 Gu n+1 pred = 0. (62) 2. Seek u n+1 , un+1 , ün+1 int , ün+1 con ∈ R N d , and r n+1 ∈ R Nc such that M ün+1 + Ku n+1 = F n+1 + T Gr n+1 , ( 63
)
Gu n+1 ≥ 0, r n+1 ≥ 0, T r n+1 Gu n+1 = 0, ( 64
)
u n+1 = u n+1 pred + ∆t 2 2 ün+2β int + ∆t 2 2 ün+1 con , ( 65
) un+1 = un + ∆tü n+γ int + ∆tü n+1 con , ( 66
)
where ün+1 = ün+1 int + ün+1 con and M ün+1 con = T Gr n+1 .

With this scheme, the contact pressure is now almost free of oscillations (Figure 13). Indeed, if

G i u n+1 = G i u n+1 pred = 0, then G i ün+2β int + G i ün+1 con = 0. ( 67
)
The long-time behavior is still poor (Figure 14). The additional step required by this scheme is not expensive compared with the main step, especially if the mass matrix is lumped.

Semi-explicit schemes

Now, we try to discretize the elastodynamic part of the problem with an explicit scheme, such as the central difference scheme. It is not possible to enforce an explicit contact condition.

The contact condition can be enforced implicitly (Discretization 4.7) or by a projection step (Discretization 4.8).

-0. Discretization 4.7 (Central differences-implicit contact [START_REF] Paoli | A numerical scheme for impact problems. I. The one-dimensional case[END_REF][START_REF] Paoli | A numerical scheme for impact problems. II. The multidimensional case[END_REF])

Seek u n+1 ∈ R N d and r n+1 ∈ R Nc such that M u n+1 -2u n + u n-1 ∆t 2 + Ku n = F n + T Gr n+1 , (68)
Gu n+1 ≥ 0, r n+1 ≥ 0, T r n+1 Gu n+1 = 0. (69) (70)
Discretization 4.8 (Central differences-projected contact [START_REF] Carpenter | Lagrange constraints for transcient finite element surface contact[END_REF])

1. Seek u n+1 pred ∈ R N d such that M u n+1 pred -2u n + u n-1 ∆t 2 + Ku n = F n . ( 71
)
Copyright c 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 00:1-6 Prepared using nmeauth.cls 2. Seek u n+1 ∈ R N d and r n+1 ∈ R Nc such that

u n+1 = u n+1 pred + T Gr n+1 , (72)
Gu n+1 ≥ 0, r n+1 ≥ 0, T r n+1 Gu n+1 = 0. ( 73
)
The second approach seems more explicit, but, actually, the two approaches are equivalent when the mass matrix M is diagonal. Discretization 4.7 has been proposed in [START_REF] Schatzman | Numerical approximation of a wave equation with unilateral constraints[END_REF] in 1D and [START_REF] Paoli | A numerical scheme for impact problems. I. The one-dimensional case[END_REF][START_REF] Paoli | A numerical scheme for impact problems. II. The multidimensional case[END_REF] in the framework of rigid-body dynamics. Discretization 4.7 has been proposed in [START_REF] Carpenter | Lagrange constraints for transcient finite element surface contact[END_REF].

It is easy to check that the acceleration at the contact boundary vanishes during a contact phase. Indeed, if

G i u n+1 = G i u n = G i u n-1 = 0, then G i ün = G i u n+1 -2u n + u n-1 ∆t 2 = 0. (74)
Consequently, there are no spurious oscillations (Figure 15). The energy balance reads

E n+1 0, 1 2 -E n 0, 1 2 = T r n+2 + r n+1 2 (Gu n+1 -Gu n ) + T F n+1 + F n 2 (u n+1 -u n ). ( 75
)
The losses of energy are quite important (Figure 16). In 1D, the convergence of the discrete solutions to the continuous solution, provided a stability condition is met (the same as in the linear case), has been established in [START_REF] Schatzman | Numerical approximation of a wave equation with unilateral constraints[END_REF]. The convergence of the discrete solutions to a semi-discrete solution has been proven in [START_REF] Paoli | A numerical scheme for impact problems. I. The one-dimensional case[END_REF][START_REF] Paoli | A numerical scheme for impact problems. II. The multidimensional case[END_REF].

DISCRETIZATIONS WITH PENALTY CONTACT CONDITION

In this part we combine standard finite elements in space and a penalty approximation of the contact condition. Then, the semi-discrete problem is a mere system of ODEs.

Problem 5.1. Seek a displacement u : [0, T ] → R N d such that, ∀t ∈ [0, T ], M ü(t) + Ku(t) = f (t) + T GR ǫ (Gu(t)), ( 76
)
with the initial conditions u(0) = u 0 and u(0) = v 0 . Problem 5.1, being a system of ODEs, has one and only one solution, which is furthermore twice differentiable in time.

Implicit schemes

To begin with, we discretize Problem 5.1 with an implicit Newmark scheme.

Discretization 5.1 (Newmark) Seek u n+1 , un+1 , ün+1 ∈ R N d , such that M ün+1 + Ku n+1 = F n+1 + T GR ǫ (Gu n+1 ), (77)
u n+1 = u n + ∆t un + ∆t 2 2 ün+2β , (78) un+1
= un + ∆tü n+γ . ( 79
)
We observe that the penalty formulation tends to reduce spurious oscillations (Figure 17). Nevertheless, the oscillations grow with the penalty parameter 1/ǫ (Figure 18). This is not surprising since the penalty contact condition tends to the exact contact condition when 1/ǫ tends to infinity. If the oscillations are too large, stabilization procedures can be used (see for instance the procedure described in [START_REF] Armero | Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems[END_REF]). With the addition of a penalty term, the Newmark scheme (trapezoidal rule) no longer conserves the energy (Figure 19). In [START_REF] Armero | Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems[END_REF][START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF], the authors proposed a discretization of the penalty term which enables to recover energy conservation (Discretization 5.2). It is based on a midpoint scheme. On our benchmark problems, it does not yield significantly better results.

Discretization 5.2 (Energy-conserving midpoint [START_REF] Armero | Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems[END_REF][START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF]) where

Seek u n+1 , un+1 ∈ R N d , such that M ün+ 1 2 + Ku n+ 1 2 = F n+ 1 2 + T G Rǫ (Gu n+1 , Gu n ), (80)
u n+1 = u n + ∆t un+ 1 2 , ( 81
) un+1 = un + ∆tü n+ 1 2 , (82)
( Rǫ (Gu n+1 , Gu n )) i =      1 2ǫ ((Giu n+1 ) -) 2 -((Giu n ) -) 2 Giu n+1 -Giu n if G i u n = G i u n+1 0 if G i u n = G i u n+1 ≥ 0 1 2ǫ (Gu n+1 + Gu n ) if G i u n = G i u n+1 < 0 (83) Setting E n pen := E n + 1 2ǫ ((Gu n ) -) 2
, the energy balance reads

E n+1 pen -E n pen = T F n+ 1 2 (u n+1 -u n ). ( 84
)

Explicit schemes

We can also envisage an explicit scheme for Problem 5.1.

Discretization 5.3 (Central differences)

Seek u n+1 ∈ R N d such that M u n+1 -2u n + u n-1 ∆t 2 + Ku n = F n + T GR ǫ (Gu n ). (85)
The penalty term stiffens the system of ODEs, which limits the stability domain of the schemes. This is not surprising since the penalty contact condition tends to the exact contact condition when 1/ǫ tends to infinity and an explicit exact enforcement of the contact condition is not possible. A more detailed discussion on the choice of the penalty parameter can be found in [START_REF] Belytschko | Contact-impact by the pinball algorithm with penalty and lagrangian methods[END_REF].

DISCRETIZATIONS WITH CONTACT CONDITION IN VELOCITY

In this part, standard finite elements in space are combined with an approximation of the contact condition involving the velocity.

M ü(t) + Ku(t) = f (t) + T Gr(t) (86) G u(t) ≥ 0, r(t) ≥ 0, T r(t)G u(t) = 0 ( 87
)
with the initial conditions u(0) = u 0 and u(0) = v 0 .

With this contact condition in velocity, the semi-discrete problem is much simpler than (32)-(33). Problem 6.1 is still a system of differential equations under unilateral constraints, but the constraint involves now the velocity instead of the displacement. The general theory developed in [START_REF] Duvaut | Inequalities in mechanics and physics[END_REF][START_REF] Glowinski | Numerical analysis of variational inequalities[END_REF] applies to Problem 6.1. The solution u is unique. Furthermore, u is continuous and u is differentiable in time almost everywhere, so that the equations have a sense at almost every time. The mathematical analysis of this problem can be found in [START_REF] Duvaut | Inequalities in mechanics and physics[END_REF].

DISCRETIZATIONS WITH MODIFIED MASS

In the previous three parts, we have considered various ways of enforcing the contact condition.

Here we describe methods based on a modification of the mass matrix. Such methods are thus compatible with any enforcement of the contact condition. For brevity, we restrict ourselves to an exact enforcement of the contact condition. In the modified mass matrix, the entries associated with the normal displacements at the contact boundary are set to zero. The motivation for this modification is very simple: if the mass is removed, the inertial forces and the oscillations are eliminated. This approach has been introduced in [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF] modified mass matrix is defined as

M * = M * * 0 0 0 .
Many choices are possible to build the block M * * . In [START_REF] Hager | A stable energy-conserving approach for frictional contact problems based on quadrature formulas[END_REF][START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF], the authors devise various methods to preserve some features of the standard mass matrix (the total mass, the center of gravity and the moments of inertia). We can also simply take the corresponding block in the standard mass matrix (and this is what will be done in our numerical simulations below). The modified problem reads

M * ü(t) + Ku(t) = F (t) + T Gr(t), (97)
Gu(t) ≥ 0, r(t) ≥ 0, T r(t)Gu(t) = 0. ( 98
) If we set u(t) = u * (t) u c (t) , K = K * * K * c K c * K cc , F (t) = F * (t) F c (t)
, and G = (G * G c ), then equations ( 97) and (98) can be recast as

M * * ü * (t) + K * * u * (t) + K * c u c (t) = F * (t), (99)
K c * u * (t) + K cc u c (t) = F c (t) + T G c r(t), (100)
G c u c (t) ≥ 0, r(t) ≥ 0, T r(t)G c u c (t) = 0. ( 101
)
For a given t and a given u * (t), there exists one and only one u c (t) satisfying ( 100) and (101). Denote by Q :

[0, T ] × R N * d → R Nc the non-linear application such that u c (t) = Q(t, u * (t)). Problem 7.1. Seek a displacement u : [0, T ] → R N d such that, for all t ∈ [0, T ], M * * ü * (t) + K * * u * (t) + K * c Q(t, u * (t)) = F * (t), (102) u c (t) = Q(t, u * (t)), ( 103
)
with the initial conditions u(0) = u 0 and u(0) = v 0 .

The operator Q(t, •) is Lipschitz continuous at each time t, so that equation ( 102) is a Lipschitz system of ODEs. Therefore, it has a unique solution u * , twice differentiable in time. Owing to (103), u c is differentiable in time almost everywhere. The detailed mathematical analysis of the semi-discrete modified mass formulation can be found in [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF][START_REF] Doyen | Convergence of a space semi-discrete modified mass method for the dynamic signorini problem[END_REF]. A result of convergence of the semi-discrete solutions to a continuous solution is proven for viscoelastic materials in [START_REF] Doyen | Convergence of a space semi-discrete modified mass method for the dynamic signorini problem[END_REF].

Remark 7.1. In contrast to the semi-discrete problem with standard mass matrix, the semidiscrete problem with modified mass matrix does not require an impact law and conserves the energy.

Implicit schemes

Let us discretize Problem 7.1 with a HHT-Newmark scheme.

The equations can be recast as a linear complementarity problem,

M * ün+1 + Ku n+1+α = F n+1+α + T Gr n+1 , (108)
Gu n+1 ≥ 0, r n+1 ≥ 0, T r n+1 Gu n+1 = 0, (109) In spite of the modification of the mass matrix, the problem is well posed. In practice, we use this set of equations to compute the solution. As expected, the large oscillations have disappeared during the contact phase (Figure 24). The energy conservation is also very satisfactory (Figure 25), since there holds

u n+1 = u n + ∆t un + ∆t 2 2 ün+2β , (110)
E n+1 * -E n * = T r n+1 + r n 2 (Gu n+1 -Gu n ) + T F n+1 + F n 2 (u n+1 -u n ), (112)
where E n * has the same expression than E n , except that M is replaced by M * .

Semi-explicit schemes

The modified mass method was introduced with implicit schemes. But we can also discretize Problem 7.1 with an explicit scheme, such as the central difference scheme. This yields a new semi-explicit scheme.

M * * u n+1 * -2u n * + u n-1 * ∆t 2 + K * * u n * + K * c Q(t n , u n * ) = F n * , (113)
u n+1 c = Q(t n+1 , u n+1 * ). ( 114
)
In practice, the equations are solved the following way.

Figure 2 .

 2 Figure 2. Impact of an elastic bar. Displacement at the bottom of the bar (left) and contact pressure (right).

Figure 3 .

 3 Figure 3. Bounces of an elastic bar.

Copyright c 2000

 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 00:1-6 Prepared using nmeauth.cls

Figure 4 .

 4 Figure 4. Bounces of an elastic bar. Displacement at the bottom of the bar.

Figure 5 .

 5 Figure 5. Impact of a point mass.

Discretization 4 . 1 (

 41 HHT-Newmark) Seek u n+1 , un+1 , ün+1 ∈ R N d , and r n+1 ∈ R Nc such that

Figure 6 .

 6 Figure 6. Impact of an elastic bar. Displacement u n 0 (left) and contact pressure r n (right). Discretization 4.1 with α = 0, β = 1/4 and γ = 1/2. ∆x = 0.1, ∆t = 0.005, νc = 1.5.

Figure 7 .

 7 Figure 7. Bounces of an elastic bar. Displacement u n 0 (left) and energy E n -T F u n (right). Discretization 4.1 with α = 0, β = 1/4 and γ = 1/2. ∆x = 0.1, ∆t = 0.005, νc = 1.5.

Figure 8 .

 8 Figure 8. Impact of an elastic bar. Displacement u n 0 (left) and contact pressure r n (right). Discretization 4.1 with α = -0.2, β = 1/4(1 -α) 2 , γ = 1/2 -α. ∆x = 0.1, ∆t = 0.005, νc = 1.5.

Figure 9 .

 9 Figure 9. Bounces of an elastic bar. Displacement u n 0 (left) and energy E n -T F u n (right). Discretization 4.1 with α = -0.2, β = 1/4(1 -α) 2 , γ = 1/2 -α. ∆x = 0.1, ∆t = 0.005, νc = 1.5.

Figure 10 .

 10 Figure 10. Impact of an elastic bar. Displacement u n 0 (left) and contact pressure r n (right). Discretization 4.3 with α = 0, β = 1/4 and γ = 1/2. ∆x = 0.1, ∆t = 0.005, νc = 1.5.

Figure 11 . 5 .

 115 Figure 11. Impact of an elastic bar. Contact pressure r n (left) and stress (Ku n )0 (right). Discretization 4.4. ∆x = 0.1, ∆t = 0.005, νc = 1.5.

Figure 12 .

 12 Figure 12. Bounces of an elastic bar. Displacement u n 0 (left) and energy E n -T F u n (right). Discretization 4.4. ∆x = 0.1, ∆t = 0.005, νc = 1.5.

Figure 13 .

 13 Figure 13. Impact of an elastic bar. Displacement u n 0 (left) and contact pressure r n (right). Discretization 4.6 with β = 1/4 and γ = 1/2 (lumped mass matrix). ∆x = 0.1, ∆t = 0.005, νc = 1.5.

Figure 14 .

 14 Figure 14. Bounces of an elastic bar. Displacement u n 0 (left) and energy E n -T F u n (right). Discretization 4.6 with β = 1/4 and γ = 1/2 (lumped mass matrix). ∆x = 0.1, ∆t = 0.005, νc = 1.5.

Figure 15 .

 15 Figure 15. Impact of an elastic bar. Displacement u n 0 (left) and contact pressure r n (right). Discretization 4.7 (lumped mass matrix). ∆x = 0.1, ∆t = 0.0025, νc = 0.75.

Figure 16 .

 16 Figure 16. Bounces of an elastic bar. Displacement u n 0 (left) and energy E n -T F u n (right). Discretization 4.7 (lumped mass matrix). ∆x = 0.1, ∆t = 0.0025, νc = 0.75.

Figure 17 .Figure 18 .

 1718 Figure 17. Impact of an elastic bar. Displacement u n 0 (left) and contact pressure T GRǫ(Gu n ) (right). Discretization 5.1 with α = 0, β = 1/4, γ = 1/2, ǫ = 10 -4 . ∆x = 0.1, ∆t = 0.005, νc = 1.5.

Figure 19 .

 19 Figure 19. Bounces of an elastic bar. Displacement u n 0 (left) and energy E n -T F u n (right). Discretization 5.1 with α = 0, β = 1/4 and γ = 1/2, ǫ = 10 -4 . ∆x = 0.1, ∆t = 0.005, νc = 1.5.

Problem 6 . 1 .

 61 Seek a displacement u : [0, T ] → R N d and a contact pressure r : [0, T ] → R Nc such that, for almost every t ∈ [0, T ],

Figure 20 .

 20 Figure 20. Impact of an elastic bar. Displacement u n 0 (left) and stress (right). Discretization 6.1 with α = 0, β = 1/4 and γ = 1/2. ∆x = 0.1, ∆t = 0.005, νc = 1.5.

Figure 21 .

 21 Figure 21. Bounces of an elastic bar. Displacement u n 0 (left) and energy E n -T F u n (right). Discretization 6.1 with α = 0, β = 1/4 and γ = 1/2. ∆x = 0.1, ∆t = 0.005, νc = 1.5.

Figure 22 .

 22 Figure 22. Impact of an elastic bar. Displacement u n 0 (left) and contact pressure r n (right). Discretization 6.2 (lumped mass matrix). ∆x = 0.1, ∆t = 0.0025, νc = 0.75.

Figure 23 .

 23 Figure 23. Bounces of an elastic bar. Displacement u n 0 (left) and energy E n -T F u n (right). Discretization 6.2 (lumped mass matrix). ∆x = 0.1, ∆t = 0.0025, νc = 0.75.

Discretization 7 . 1 (

 71 HHT-Newmark[START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF])Seek u n+1 ∈ R N d , un+1 * ∈ R N d and ün+1 * ∈ R N d such that M * * ün+1 * + K * * u n+1+α * + K * c Q(t n+1+α , u n+1+α * ) = F n+1+α * ,

  un+1 = un + ∆tü n+γ .(111) Copyright c 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 00:1-6 Prepared using nmeauth.cls

Figure 24 .

 24 Figure 24. Impact of an elastic bar. Displacement u n 0 and contact pressure r n . Discretization 7.1 with α = 0, β = 1/4 and γ = 1/2. ∆x = 0.1, ∆t = 0.005, νc = 1.5.

Figure 25 .

 25 Figure 25. Bounces of an elastic bar. Displacement u n 0 (left) and energy E n * -T F u n (right). Discretization 7.1 with α = 0, β = 1/4 and γ = 1/2. ∆x = 0.1, ∆t = 0.005, νc = 1.5.

Figure 26 .

 26 Figure 26. Bounces of an elastic bar. Displacement u n 0 (left) and contact pressure r n (right). Discretization 7.2 (lumped mass matrix). ∆x = 0.1, ∆t = 0.0025, νc = 0.75.

Figure 27 .

 27 Figure 27. Bounces of an elastic bar. Displacement u n 0 (left) and energy E n * -T F u n (right). Discretization 7.2 (lumped mass matrix). ∆x = 0.1, ∆t = 0.0025, νc = 0.75.

1 .Figure 28 .

 128 Figure 28. Impact of an elastic bar. Stress in the bar at time t=1.0. Discretization 7.1. ∆x = 0.1. ∆t = 0.005, νc = 1.5 (left). ∆t = 0.00125, νc = 0.375 (right).

Figure 29 .

 29 Figure 29. Impact of an elastic bar. Stress in the bar at time t=1.0. Discretization 7.1 with α = -0.06, β = 1/4(1 -α) 2 , γ = 1/2 -α. ∆x = 0.1. ∆x = 0.1. ∆t = 0.005, νc = 1.5 (left). ∆t = 0.00125, νc = 0.375 (right).

  ).

	0000000 0000000 1111111 1111111 0000000 0000000 1111111 1111111 000000 000000 111111 111111 0000000 0000000 1111111 1111111 000000 000000 111111 111111 0000000 0000000 1111111 1111111

  and r n+1 ∈ R Nc satisfying (38)-(41). Prepared using nmeauth.cls 2. If G i u n < 0 and G i u n+1 = 0, then un+1 and ün+1 are modified so that

							G i	un+1 = 0	and	G i	ün+1 = 0.							(49)
		0.4										500
		0.3										400
	displacement	0 0.1 0.2									contact pressure	100 200 300
												0
		-0.1
						ex./naive stabilized Newmark exact solution			-100				ex./naive stabilized Newmark exact solution
		-0.2
		0.4	0.5	0.6	0.7	0.8	0.9	1	1.1	1.2	1.3	0.4	0.5	0.6	0.7	0.8	0.9	1	1.1	1.2	1.3
						time										time
	Copyright c 2000 John Wiley & Sons, Ltd.				Int. J. Numer. Meth. Engng 2000; 00:1-6

  The time discretization

		0.4										500
		0.3										400
	displacement	0 0.1 0.2									contact pressure	100 200 300
												0
		-0.1
								vel./midpoint exact solution			-100						vel./midpoint exact solution
		-0.2
		0.4	0.5	0.6	0.7	0.8	0.9	1	1.1	1.2	1.3	0.4	0.5	0.6	0.7	0.8	0.9	1	1.1	1.2	1.3
						time										time

  ]. Set N * d := N d -N c .For the sake of simplicity, suppose that the degrees of freedom associated with normal displacements at the contact boundary are numbered from N * d + 1 to N d . The

		0.4										500
		0.3										400
	displacement	0 0.1 0.2									contact pressure	100 200 300
												0
		-0.1
						vel./central differences exact solution			-100				vel./central differences exact solution
		-0.2
		0.4	0.5	0.6	0.7	0.8	0.9	1	1.1	1.2	1.3	0.4	0.5	0.6	0.7	0.8	0.9	1	1.1	1.2	1.3
						time										time

grant sponsor: Publishing Arts Research Council; contract/grant number: 98-1846389

has been extensively studied in [START_REF] Glowinski | Numerical analysis of variational inequalities[END_REF]. Unfortunately, the contact condition in velocity is not equivalent to the Signorini condition. This approximation is valid only on a short time interval after an impact. The strategy adopted is the following: if a node satisfies the non-interpenetration condition, then at the next iteration no constraint is enforced on this node; if a node breaks the non-interpenetration condition, then at the next iteration the contact condition in velocity will be applied to this node. This approach allows for slight interpenetration. At each time step, we define the matrix G n whose rows G n i are

This approach based on a contact condition in velocity has also been widely used in rigid-body dynamics with impacts.

Implicit schemes

A midpoint scheme with a contact condition velocity has been proposed in [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF]. Discretization 6.1 (Midpoint [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF]) Seek u n+1 , un+1 ∈ R N d , and

An interesting feature of this scheme is to be energy-conserving,

The contact pressure does not perturb the structure despite its oscillations (Figure 20). Owing to the energy conservation, the long-time behavior is quite satisfactory (Figure 21).

Semi-explicit schemes

In [START_REF] Belytschko | Contact-impact by the pinball algorithm with penalty and lagrangian methods[END_REF], an explicit scheme based on the contact condition in velocity has been proposed. Discretization 6.2 (Central differences [START_REF] Belytschko | Contact-impact by the pinball algorithm with penalty and lagrangian methods[END_REF])

Numerical simulations suggest that the stability condition of the central difference scheme is not tightened by the contact condition. The numerical results are similar to those obtained with Discretization 4.7 (Figures 22 and23).

Seek u

The first step is explicit, provided the mass matrix M * * is lumped. The second step is a constrained problem which concerns only the variable u c . Discretization 7.2 satisfies the system

and the energy balance,

where E n 0, 1 2 * has the same expression than E n 0, 1 2 , except that M is replaced by M * . We observe numerically that the stability condition on the time step is the same as in the linear case. Compared with Discretizations 4.7 and 6.2, the semi-explicit modified mass method shows a better energy conservation and a better long-time behavior (Figure 27).

APPENDIX

When approximating a solution presenting a shock or a sharp wave front with finite elements in space and a time-stepping scheme, spurious oscillations are observed. This is the wellknown Gibbs phenomenon which is due to the poor approximation of eigenmodes associated with the high frequencies. In the elastodynamic Signorini problem, we must deal with shock waves after the impacts. If we plot the stress distribution computed with the modified mass method and a Newmark scheme, we see this phenomenon (Figure 28). Moreover, the magnitude of the oscillations does not decrease with the time step. For instance, as observed in [START_REF] Grosu | Stability of semidiscrete formulations for elastodynamics at small time steps[END_REF], with small time steps, the oscillations appears in front of the shock instead of behind. These oscillations can be eliminated by using dissipative schemes, such as the HHT or Chung-Hulbert schemes (or, which is equivalent, by filtering). For instance, the damping obtained with a HHT scheme is presented in Figure 29. It may seem a little bit contradictory to discard the use of dissipative schemes for the contact oscillations and to reintroduce them for the Gibbs phenomenon. Actually, the oscillations on the stress caused by the Gibbs phenomenon are far smaller than those caused by the contact condition (compare Figure 6 with Figure 28). The former are proportional to the discontinuity jump, the latter are proportional to mivi ∆t , where v i is the velocity before the impact and m i is the mass associated with the impacting node. Consequently, the amount of dissipation needed is much smaller.