Chih-Chen Kao

Wei-Chung Hsu

An Adaptive Heterogeneous Runtime for Irregular Applications in the Case of Ray-Tracing (Extended Abstract)

Heterogeneous architecture has been widely adopted in various computing systems, from mobile devices to servers. However, optimizing the performance for such platforms remains challenging in three aspects: the control flow divergence decreases the utilization of SIMD components, the significant memory copy overhead between computing devices consuming precious memory bandwidth and the load imbalance that degrades the overall performance. In this paper, we proposed three methodologies: Intermediate Feedback, Dynamic Task Partitioning and Heterogeneous Runtime that work collaboratively to overcome the aforesaid problems. We adopted and implemented these methodologies in a heterogeneous runtime library derived from Intel Embree [START_REF] Woop | Embree ray tracing kernels for cpus and the xeon phi architecture[END_REF] and compared the performance results of the two frameworks running Ray-Tracing [START_REF] Purcell | Ray tracing on programmable graphics hardware[END_REF] on various scenes. Experiment results have shown that the performance gain from the proposed methods is significant, especially in complex scenes with a large amount of objects or with large input data sizes the CPU cannot handle efficiently.

Due to the performance and power efficiency potentials of GPGPU and heterogeneous systems, a wide variety of applications, which include molecular simulation, fluid dynamics, biomedical image processing and computer vision, have been developed by leveraging the aforementioned programming models. However, for this type of heterogeneous configuration, many challenges remain before the performance potential can be fully unleashed. Despite significant advantages of GPU programming, writing high-performance heterogeneous programs still require programmers to be familiar with GPU architecture. The performance of a heterogeneous program is significantly influenced by how the computations are mapped into threads and how those threads are scheduled onto distinct cores, the usage of GPU registers and memory hierarchy, the synchronization among all threads and data access, the data transfer between host and GPU memories and the control flow branch divergence issue. In short, a program which benefits from GPGPU must contain explicit parallelism, high regularity and data reuse. The irregularities in an application may throttle the expected performance of the GPU by as much as an order of magnitude.

The terms regular and irregular are often used in compiler literature. For example, in regular code, control flow and data memory references are not data dependent. Dense matrix multiplication operations are good examples of regular code. On the other hand, in irregular code, both control flow and data memory references could be data dependent. For example, graph-based applications are considered irregular, because the connectivity and values of the nodes in a graph are unknown before the input graph is available, and the connectivity and values of the nodes determine which graph elements are accessed [START_REF] Burtscher | A quantitative study of irregular programs on gpus[END_REF]. Regular programs can often be efficiently mapped onto the GPGPU computing unit, and many of them have been ported to heterogeneous systems to benefit from increased performance and power efficiency. However, duplicating the success of heterogeneous computing from regular programs to irregular programs is a challenge. Irregular programs are often operated around pointer-based data structures such as graphs, trees or linked lists that are difficult to map to conventional regular GPU architecture or SIMD components. The execution paths of an irregular program are often unpredictable and could also vary dramatically during runtime [START_REF] Zhang | On-the-fly elimination of dynamic irregularities for gpu computing[END_REF]. Furthermore, in contrast to regular programs, the data dependencies in irregular programs can only be resolved dynamically at runtime, making it difficult to design a proper task scheduler [START_REF] Nasre | Data-driven versus topology-driven irregular computations on gpus[END_REF] [START_REF] Monteiro | A pattern language for parallelizing irregular algorithms[END_REF].

The obstacles of efficiently adopting irregular programs to heterogeneous system can be categorized as follows: Firstly, the irregular memory access (e.g. indirect array references, sparse matrix references) could lead to low effective memory bandwidth, and there is no proper static solution to adaptively determine which type of memory should be used for irregular programs during execution. Secondly, the dynamically varying control flows create thread divergences, which reduces the level of parallelism and SIMD lane utilization in GPU. Thirdly, the chain of input dependencies often causes load imbalance among multiple cores, which degrades the overall throughput. Finally, memory-bound pointer chasing exhibits low data locality and exposes increased data access latency on GPGPU. Therefore, mapping irregular code efficiently onto a heterogeneous system remains difficult [START_REF] Pingali | The tao of parallelism in algorithms[END_REF] [START_REF] Kulkarni | How much parallelism is there in irregular applications?[END_REF].

What type of heterogeneous computing model could handle irregular programs more efficiently is still debatable. For example, Ray-Tracing is an irregular program that is intensively used for global illumination in multimedia applications and has been adopted for implementation on different heterogeneous models. At the early stage of GPGPU computing, Ray-Tracing is designed to run on a GPGPU because a GPGPU is capable of handling a large number of rays in parallel. However, the potential irregularities of Ray-Tracing decrease the utilization of SIMD/SIMT components in GPU, and thereby offer little performance/power advantages [START_REF] Aila | Understanding the efficiency of ray traversal on gpus[END_REF].

Recently, Intel has announced a Ray-Tracing framework called Embree [START_REF] Woop | Embree ray tracing kernels for cpus and the xeon phi architecture[END_REF] with a different design philosophy that is optimized for traditional CPU architectures augmented with medium size (i.e. 512 bits) of SIMD capability. The Embree framework leverages CPU threads with wider SIMD units by using a compiler framework called Intel SPMD Program Compiler. (ISPC) [START_REF] Pharr | ispc: A spmd compiler for high-performance cpu programming[END_REF]. The framework was originally designed for single ray traversal using SSE or AVX-enabled CPUs but has been extended to support Intel Xeon Phi architectures [START_REF] Chrysos | Intel R xeon phi coprocessor (codename knights corner)[END_REF]. Embree features spatial acceleration structures and traversal algorithms and claim to support efficient Ray-Tracing with MIMD architectures and medium size SIMD capability. However, we believe that, with proper runtime and innovative method-ologies, we could make traditional CPU/GPU type of heterogeneous system more competitive on irregular applications. In this work, we use Ray-Tracing as a case study to show the potential of our proposed approaches.

In this work, we address the forenamed problems of irregular programs and proposed a feedback tuning mechanism that can be used to model a specific category of heterogeneous program where the input data is recursively modified and added back to the commonly shared database for the next computation. Previous research applied a statistic approach and heuristic. Our method is based on analyzing and monitoring the communication protocol and behavior among all modules in a program. A new methodology is introduced in this research that encodes the representative feature of a heterogeneous program gathered at runtime and sends them to adjacent modules for adjusting iterative computation to fit the given platform configuration in order to gradually fine tune the system performance. We proposed a dynamic task partitioning mechanism and heterogeneous thread pool in order to resolve issues related to branch divergence and load imbalance. To demonstrate the effectiveness of our proposed schemes, we evaluate our performance gain by implementing the methodologies into a runtime library which is derived from the Intel Embree Ray-Tracing framework and compare the execution time with the original version.

In Conclusion, we explore the performance potential of mapping irregular programs onto heterogeneous systems by taking a Ray-Tracing algorithm as a case study. Three methodologies, Intermediate Feedback, Dynamic Task Partitioning and Heterogeneous Thread Pool, are introduced in this framework. The experiment results shows that our proposed methods could benefit heterogeneous computing resources and thereby increase the system performance, especially for handling complex scenes and for large input data sizes. We believe that the proposed methods could be applied to other heterogeneous frameworks to address the challenges of branch divergences, memory copy overhead and load imbalance when mapping irregular applications to GPGPU.

The benefit of the intermediate feedback is significant. Without this mechanism, the program must be built by pure heuristic approaches, which have limited success and are often effective for only certain specific configurations. The optimized execution setup will be lost if the system setting is changed. For instance, in the case of Embree, the hybrid packet/single-ray tracing algorithm is implemented by utilizing a specific type of BVH tree. However, determining the appropriate tracing method in the algorithm is based on heuristic. The process starts with a 16-wide packet traversal which performs 16-wide box tests. At any point in time, the bit in an active mask will be counted to indicate how many of the packet's rays are still active for a subtree. If this number falls below a given threshold, which is set to 7, the process leaves the packet traversal mode and sequentially traces all active rays in the single-ray mode [START_REF] Benthin | Combining single and packet-ray tracing for arbitrary ray distributions on the intel mic architecture[END_REF]. The drawback of this method is that the threshold may need to change and the program will require recompilation if the system is moved to a machine with a shorter SIMD lane. Also, this method does not prevent any possible divergent execution that lowers the effectiveness of the SIMD engine. Our intermediate feedback differs from the above method by setting the correlation based on the actual activities and data analysis. It is effective for any kind of system configuration since the runtime would automatically adjust itself. If this method were moved to another platform, the program which utilized our runtime library could adapt and gradually arrives at an optimized setting.


Mingming Sun
email: mmsun@mail.ustc.edu.cn

Changlong Li

Xuehai Zhou
email: xhzhou@ustc.edu.cn

Kun Lu

Hang Zhuang
email: zhuangh@mail.ustc.edu.cn

DLBer: A Dynamic Load Balancing Algorithm for the Event-driven Clusters

The event-driven programming model has been proposed to efficiently process iterative applications and incremental applications. In clusters based the event-driven model, applications are structured as a series of triggers, each of which will be invoked when associate events are trigged. And framework assigns a newly submitted trigger to a node where the relevant datasets set. Unfortunately it may lead to load imbalance because associate events occur by chance. Numerous triggers in a node may be simultaneously invoked but other nodes have no triggers running. Jobs composed of short, sub-second triggers present a difficult balancing challenge. To the end, we design DLBer, a new dynamic load balancing algorithm for the event-driven clusters to maximize improve the utilization of node resources.

Introduction

The synchronous data-flow model such as MapReduce, Dryad and their variants is deficient for iterative applications since synchronous computation and lockstep across rounds. It also is not suitable for incremental applications because it processes total dataset for every increment, incurring a significant performance penalty. Therefore, the event-driven programming model is proposed as an asynchronous computation model such as Percolator [START_REF] Peng | Large-scale Incremental Processing Using Distributed Transactions and Notifications[END_REF], Oolongr [START_REF] Mitchell | Oolong: asynchronous distributed applications made easy[END_REF] and Dominor [START_REF] Dai | Domino: an incremental computing framework in cloud with eventual synchronization[END_REF] et al. This model follows the Event-Condition-Action (ECA) rule [START_REF] Mccarthy | The architecture of an active database management system[END_REF]. ECA rules are straightforward: when the event occurs, evaluate the condition; if the condition is fulfilled, execute the action automatically. Applications can be expressed in term of triggers, user-specified code blocks that can be invoked whenever the associated datasets modified. Each trigger completes a computation task. In iterative applications, the result of past iteration can be immediately used to determine the course of current execution. Incremental applications only need to process the updated dataset and some relevant data, which is affected by the updated dataset, instead, recalculates entire dataset.

In clusters based event-driven model, a newly submitted trigger is assigned to the node where the relevant datasets set. However triggers are not instantly executed after their submission, that is to say, specific events are needed to invoke triggers. This may lead to load imbalance because during execution time the cluster resource usage cannot be determined. Extreme case, a node has lots of triggers invoked in the waiting queue, but other nodes are idle. So framework needs a load balancer to achieve better utilization of node resources and a high system throughout and quick response time of user requests. However, trigger may be short second or sub-second task. This presents a difficult balancing challenge. To address the above problems, we design DLBer, a new centralized dynamic load balancing algorithm for the event-driven clusters. Dynamic load balancing is essential for such systems since unpredictable load estimates. Our load balancer aims at maximize improve the utilization of node resources.

DLBer Design

As shown as Fig. 1, custer based event-driven model runs with a T riggerM aster, which assigns a newly coming trigger to workers where the relevant datasets of the trigger set, and multiple T riggerW orkers, each of which monitors the modification on dataset, and executes trigger tasks. Our DLBer consists of two main modules: a LoadBalancer in T riggerM aster, which coordinates and redistributes triggers invoked among T riggerW orkers according to their workloads, and LocalScheduler in each T riggerW orker. LocalScheduler contains W orkloadEstimator, which automatically predicts the execution time of trigger and evaluates current total workload, and T riggerGrouper, which picks out several suitable triggers to become a group as the basic unit of transfer. As trigger in the event-driven clusters may be short second or sub-second task, we divide triggers invoked into groups to avoid repeatedly network transmission and reduce pressure for T riggerM aster. Once the balancing decisions are made in LoadBalancer, a high-load T riggerW orker will receive transfer instructions and trigger LocalScheduler to transfer triggers. LocalScheduler first selects one or more trigger group, and then transfer trigger group to low-load destination T riggerW orker according to transfer instructions. DLBer employs centralized load balancing policy because the accumulation of workoad information can be achieved by heartbeat message.

Trigger Group Policy

W orkloadEstimator calculates the mean completed time of same triggers, which have already completed, as the execution time of trigger. As triggers continue to complete, the execution time of trigger will be recalculated to get more precise time. Each T riggerW orker workload is the sum of the execution time of all triggers invoked in waiting queue.

The role of T riggerGrouper is to divide triggers invoked into groups as the basic unit of transfer. We define T riGroup, the size of trigger group, which depends on the worker computing capacity and network bandwith. Upon the arrival of a transfer request from T riggerM aster, T riggerGrouper sorts triggers in descending order of predicted execution time into trigger sequence. Then it selectes triggers from head of trigger sequence into a trigger group according to T riGroup. Then trigger group contains minimal triggers to reduce network transmission overhead.

Load Balancing Algorithm

DLBer focuses on T riggerW orkers, whose workload is lower than a threshold called tunder. tunder is greater than the sum of load balancing policy time, the heartbeat interval (default as 5 second) and transmission time of a trigger group. Then the low-load T riggerW orker can be still busy during the transfer triggers.

Algorithm 1 illustrates the process of load balancing algorithm. If workload of worker n is below the tunder, load balancing policy first considers the heaviest load worker a. Supposing worker a will transfer out a trigger group, the workload of worker a needs to exceed the load of worker n, otherwise it is not necessary to transfer tasks. When all nodes in cluster are busy, there is almost no extra scheduling overhead. At the other extreme, our method also avoid frequent task transfering since the defined tunder.

Evaluation

In order to test DLBer using a realistic workload, we ported Domino by writing a Domino load balancing plugin. Domino is an open-source trigger-based programming framework. We use application PageRank on a cluster, which has a master and 8 workers. PageRank in the event-driven model contains two triggers: P ageRankDist and P ageRankSum. Once pagerank value of a page is modified, P ageRankDist will be invoked to change pagerank weights for all its relevant out-degree pages. Then the P ageRankSum for such relevant outdegree pages will be invoked for these changes. The function of P ageRankSum is to simply sum up all the rank values generated from different in-degree pages.

Algorithm 1 : Load Balancer

Input:

Array A: nodes queue by descending workload; Iteration:

while n.load < tunder ( n : the tail of A) do if a.load > tunder & a.load -T riGroup > n.load + T riGroup ( a : the head of A) then task f etch(a, n); else break; end if Reorder array by descending workers' workload; end while P ageRankDist will be invoked again. This continues until a tolerable error defined by user. Our results shows our load balancer outperforms default Domino about 100%.

Fig. 1 .

 1 Fig. 1. The architecture of DLBer

Acknowledgment Our work was supported by the National Science Foundation of China under grants No. 61272131 and No. 61202053, China Postdoctoral Science Foundation grant No. BH0110000014, Fundamental Research Funds for the Central Universities No. WK0110000034, and Jiangsu Provincial Natural Science Foundation grant No. SBK201240198.


Anders Thorin

O Tatin

Nonsmooth modal analysis from elementary impact oscillators to turbomachines: limitations and avenues

Keywords:

In mechanical engineering, nonlinear modes characterise the behaviour of nonlinear vibratory systems. In the current state-of-the-art, they are well defined for smooth nonlinear systems (of moderate size) of Ordinary Differential Equations governing the dynamics. They are defined as continua of periodic orbits forming two-dimensional invariant manifolds in the state space. This framework has lately been extended to nonsmooth mechanical systems involving impact dynamics.

From the theoretical standpoint, strong mathematical results about existence, uniqueness and analytical or approximate equations of nonsmooth modes have recently been established on a linear spring-mass chain undergoing a Newton elastic impact law on one of its masses. From the industrial point of view, their application is not straightforward. This paper investigates the possibilities and the limitations of these tools with a practical end: the issue of blade-casing unilateral contact interactions in turbomachines.

The main differences between a simple spring-mass chain and complex blade-casing models are investigated point by point: non-diagonal mass matrix, geometrical differences, scalability for thousands of degrees of freedom, convergence with respect to the number of dofs, stability and relationships with force and damped mechanical systems.

It is found that the proposed formulation and corresponding solutions apply to very general springmass systems, involving non-diagonal mass matrices and non-tridiagonal stiffness matrices, opening avenues to investigate complex industrial systems. However, the relationship between the forced and damped behaviour and nonsmooth modes is not fully understood yet.

INTRODUCTION

Nonlinear modes of vibration provide useful means of investigating the dynamics of nonlinear systems [START_REF] Rosenberg | On nonlinear vibrations of systems with many degrees of freedom[END_REF]. They can be defined as a two-dimensional continuum of periodic trajectories lying on invariant manifolds in the state space, whose coordinates are the chosen generalised coordinates and generalised velocities. The "invariant"' property refers to the fact that if the state of a mechanical system is on its nonlinear mode at a given time, it will remain on it as time unfolds. Thus, nonlinear modes are particular surfaces in the state space which characterise, locally, the behaviour of a dynamical system.

For smooth nonlinear systems, nonlinear modes can be approximated locally by choosing a pair of master coordinates [START_REF] Shaw | Non-linear normal modes and invariant manifolds[END_REF][START_REF] Kerschen | Nonlinear normal modes, part I: A useful framework for the structural dynamicist[END_REF], writing the generalised coordinates and velocities as unknown functions of this pair and calculating a Taylor expansion of these unknown functions. This results in a parametric equation formulated in the pair of master coordinates, yielding one two-dimensional surface in the state space: a nonlinear mode.

A considerable limitation of this method is that it relies on the smoothness of the differential equations (ODEs) and no longer applies to nonsmooth mechanical systems, i.e. systems undergoing impact or dry friction. Yet, many vibratory systems undergo unilateral contact conditions such as the common blade-casing contact occurrences in turbomachines. The operating clearance between blades and surrounding casing has to be as small as possible for efficiency purposes, potentially initiating contact induced vibratory phenomena [START_REF] Almeida | Experimental analysis of dynamic interaction between a centrifugal compressor and its casing[END_REF]. Because of the velocity discontinuities induced by impacts, such events are no longer captured by ODEs and require more complex mathematical tools such as linear complementary conditions, inclusions in cones, inclusions of differential measures or variational inequalities [START_REF] Acary | Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics[END_REF]. Modal analysis is currently unable to properly handle such systems.

A simple spring-mass oscillator has been recently studied [START_REF] Legrand | Nonsmooth modal analysis of a N degree-of-freedom system undergoing a purely elastic impact law[END_REF][START_REF] Thorin | Nonsmooth modal analysis: Investigation of a 2-dof spring-mass system subject to an elastic impact law[END_REF]. The oscillator is subject to unilateral contact with a Newton impact law: v(t + imp ) = -v(t - imp ) where t imp is any time of impact, v(t - imp ) is the pre-impact velocity and v(t + imp ) is the post-impact velocity. The number of impact(s) per period arises a as a natural criterion for the classification of nonsmooth modes. Periodic solutions with one impact per period (ipp) have been thoroughly investigated in [START_REF] Legrand | Nonsmooth modal analysis of a N degree-of-freedom system undergoing a purely elastic impact law[END_REF], leading to 1 ipp nonsmooth-modes. In [START_REF] Thorin | Nonsmooth modal analysis: Investigation of a 2-dof spring-mass system subject to an elastic impact law[END_REF], a formulation for 2 ipp is proposed. More general results for n dof and k ipp have been found and will be published soon.

This contribution targets the implementation of modal analysis in the framework of rotating machinery where unilateral contact events are known to emerge. It focuses on current challenges and limitation by paying attention to simple qualitative behaviours. Questions of interest are the theoretical and practical applicability of nonsmooth modes to generic mechanical systems featuring potential large arbitrary mass and stiffness matrices.

The main mathematical results, unpublished yet, of nonsmooth modal analysis for a spring-mass chain are provided in Section 1. They serve as a starting point for the present study and a basic example is proposed for illustrative purposes. In Section 2, the main steps required to extend the results to a full blade-casing contact model are identified and investigated.

STARTING POINT: MATHEMATICAL RESULTS

In this section, the main mathematical results about nonsmooth modes are given. For the sake of conciseness, only the main ingredients are listed in a quite abstract fashion, limiting the presentation to the important results. In substance, the latter prove existence and uniqueness, under some conditions, of periodic solutions to the problem described in subsection 1.1 for k impact(s) per period. They also yield an analytical (k = 1) or approximate (k ≥ 2) expressions of the nonsmooth modes.

Model and dynamical equations

Consider the spring-mass system in Fig. 1. The quantities

m i m 1 k 1 m 2 k 2 m n-1 k n-1 m n k n d Figure 1.
Base model for the study of nonsmooth modes. are masses and k i are stiffnesses. The position of the mass i is stored in x i and by convention x i = 0 at rest, for i ∈ 1, n . The position x n of the last mass is constrained by an obstacle located at position d, which is described mathematically by the following Signorini condition such that ∀t ≥ 0:

d -x n (t) ≥ 0 (1a) λ(t) ≥ 0 (1b) λ(t)(d -x n (t)) = 0 (1c)
This unilateral condition has to be supplemented by an impact law between the n-th mass and the obstacle: to this end, the simple Newton impact law is selected so that the overall dynamics writes:

∀t, d -x n (t) ≥ 0 (2a) ∀t such that x n (t) < d, Mẍ + Kx = 0 (2b) ∀t such that x n (t) = d, ẋn (t + ) = -ẋn (t -) (2c)
with the mass matrix M = diag(m 1 , . . . , m n ), the stiffness matrix

K = tridiag(-k i-1 , k i-1 + k i , -k i ) and the vector of positions x = [x 1 , . . . , x n ] .

Preliminary definitions

Let us define:

• B = - √ M -1 √ K; Q an orthogonal matrix such that BB = Q∆ 2 Q -1 with ∆ a diagonal matrix; ω i > 0 the element (i, i) of ∆.
• k the number of impacts per period and t 1 , . . . , t k the instants of impacts. The duration between two successive impacts denoted by σ i = t i -t i-1 with t 0 = 0. • The period T; as by convention contact is activated at t = 0, the period is equal to the time of the k-th impact:

T = t k . • The function Q T (σ) = n i=1 sin(ω i σ/2) sin(ω i T/2) q n, j 2 (3)
where q n, j is the element (n, j) of Q. • The following vectors:

W n = - 1 √ m n              q n,1 /ω 1 0 . . . q n, n /ω n 0              R = - 1 √ m n              q n,1 0 . . . q n, n 0              (4)
• The matrix exponential

S(t) = exp t 0 I -M -1 K 0 (5)
• The skew-symmetric matrix Π populated by

Π i j (t 1 , . . . , t k ) = Q T (2(t j -t i ) -T ) (6)
for i < j ≤ k. • The symmetric matrix Σ populated by the elements

Σ i j (t 1 , . . . , t k ) = 2W n S( T 2 ) -S(-T 2 -1 S(t i -t j + T 2 )R (7) for i < j ≤ k. • The constant vector j = [1, . . . , 1] ∈ R k .

Mathematical results

A necessary condition for the existence of a periodic solution with k ipp is that there exists a λ ∈ R k such that:

Π(t 1 , . . . , t k )λ = 0 (8) Σ(t 1 , . . . , t k )λ = dj (9)
If such a condition is satisfied, then let us define:

φ p,T (σ) = - 1 √ m p n j=1 cos(ω j σ/2) ω j sin(ω j T/2) q p, j q n, j (10)
and

ψ p,T (σ) = 2 dφ p,T dt (t) = 1 √ m p n j=1 sin(ω j σ/2) sin(ω j T/2) q p, j q n, j (11)
Also, consider ∀τ ∈ [0, σ i+1 ] the p-th position:

x p (τ + t i ) = i j=1 λ j φ p,σ (2(τ + t i -t j ) -T )+ k j=i+1 λ j φ p,σ (2(τ + t i -t j ) + T ) (12) If ∀i ∈ 1, k , ∀τ ∈ [0, σ i+1 ], x n (τ + σ i+1 ) ≤ d, then x is a periodic solution of problem (2) with k impact times t 1 , . . . , t k . The velocity p is given ∀τ ∈ [0, σ i+1 ] by y p (τ + t i ) = i j=1 λ j ψ p,T (2(τ + t i -t j ) -T ) + k j=i+1 λ j ψ p,σ (2(τ + t i -t j ) + T ) (13)

Example

For illustration purposes, arbitrary choices n = 2 and k = 2 are made with

k 1 = k 2 = 1 and m 1 = m 2 = 1/2.
The definition of matrix Π, see Eq. ( 6), with t 2 = T leads to

Π(t 1 , T ) = 0 Q T (T -2t 1 ) -Q T (T -2t 1 ) 0 (14)
Condition ( 8) yields either the trivial solution λ = 0 or the more interesting equation

Q T (T -2t 1 ) = 0 (15)
From the definition of Q T given in Eq. ( 3), and the eigenfrequencies of the system ω 1 , ω 2 , it is possible to approximate (t 1 , T ) solution of (15). Note that the resonant cases are ex-

cluded (t 1 , T { 2π ω 1 , 2π ω 2 }). The roots of the numerator of Q T (T -2t 1 ) is displayed in Fig. 2.
To correctly interpret Eq. ( 8), is has to be understood that 2 ipp motions do not exist for any (t 1 , T ). Eq. ( 8) gives a necessary (but not sufficient) condition for the existence of such solutions. Note that the line of equation Tt 1 = t 1 in Fig. 2 corresponds to 1 ipp motions of period T/2 but seen as T-periodic motions.

For chosen (t 1 , t 2 ) satisfying Eq. ( 8), Eq. ( 9) yields the values of the vector λ. This forces the impact conditions at t 1 and t 2 to coincide with x n = d. The relation between the positions x p and the components of λ is in Eq. ( 12).

Both conditions ( 8) and ( 9) are not sufficient: among the motions they lead to, they do not exclude those which can violate the constraint x n ≤ d. This condition has to be tested numerically using Eq. ( 12). If it is fulfilled, then the corresponding x is a solution. By continuity, there are other solutions in the neighbourhood of the selected (t 1 , T ), yielding a continuum of periodic trajectories with two impacts per period. Samples of the continuum are plotted in Fig. 3, corresponding to points chosen on the red curve of Fig. 2. This continuum of solutions can also be represented as an Example of periodic trajectories with two impacts per period. They correspond the red curve parametrised by (t 1 , T ) on Fig. 2 and delimited by 1 and 2 . 1 corresponds to a 1 ipp trajectory taken on two impacts. 2 corresponds to the fifth harmonic of the second linear mode. In between are purely 2 ipp motions. invariant manifold, gathering the continuum of periodic orbits in the state state. This two-dimensional object, embedded in a 2n dimensional state space (n positions, n velocities) can be projected into the three-dimensional space (x 1 , x 2 , ẋ2 ). Such a projection is represented in Fig. 4. This example was used to illustrate the mathematical results and show, in practice, how they lead to the construction of nonsmooth modes for the spring-mass chain system of Fig. 1. The next part consists in exploring the applicability of these developments to more challenging models such as a blade-casing model of a turbomachine.

GENERIC MECHANICAL SYSTEMS

Though very fruitful, these results were proven for the springmass system of Fig. 1. FEM models used for turbomachines differ from such a simple model: they involve a priori nondiagonal mass matrices, contain thousands of degrees-offreedom, target much more complicated geometry involving more than two connexion per element, etc. making this method not directly applicable. These points are discussed in the present part. Nevertheless, contact is assumed to occur only on one node.

In all the computations, k = 1, m = 1/N and d = 1.

Non-diagonal mass matrices

For linear spring-mass systems, the equations governing the dynamics of every degree-of-freedom can be gathered in a matrix form with a diagonal mass matrix. When dealing with finite-elements, this property is no longer conserved: basis functions are in general not orthogonal (for the canonical scalar product), yielding off-diagonal terms. However, the mass matrix remains symmetric1. During phases of free flight (x n < d), this symmetry does not invalidate the results of Section 1. Since the proof is not given, the reader will not be able to verify this statement; the key idea is that the mass matrix is a real symmetric matrix, so is diagonalisable in an orthogonal basis. But t no longer holds for the Newton impact law: the transformation v n + = -v n -corresponds to a reflection matrix N = diag(1, . . . , 1, -1) of size 2n which is an isometry for the energy if and only if M n j = 0 for all j n. In other words, the Newton impact does not preserve energy if the mass which impacts the obstacle is not punctual. This is 1The symmetry of the mass matrix can be proven with Maxwell-Betti's theorem.

not a severe restriction. Nevertheless, our first investigations on defining an elastic impact law as an isometry for the energy and a orthogonal symmetry with respect to an hyperplane of the state space are promising and could raise this restriction. This is a complicated topic on itself and is not discussed here. Also, note that non-diagonal mass matrix are commonly lumped resulting in a diagonal mass matrix.

Geometry of the mechanical system

Here, it is shown that the results of Section 1 apply for more complex geometries than that of Fig. 1. Two changes of geometry are investigated:

• the addition of a tangential velocity at the impact;

• the connection of a mass to more than two springs.

A tangential velocity is introduced in a spring-mass system. The goal is to show that the results of Section 1 still apply. The model to capture tangential velocity is as simple as possible, see Fig. component of the velocity is unchanged while the sign of its normal component is changed. This writes:

v + = v - t -e v - n (16)
where + ( -) denotes the post-impact (pre-impact) velocity. Eq. ( 16) can be expanded as

ẋ+ ẏ+ = 1 || t|| 2 ([ ẋ ẏ] • t) t - e || n|| 2 ([ ẋ ẏ] • n) n (17) = 1 1 + a 2 1 -ea 2 a(1 + e) a(1 + e) a 2 -e N ẋ- ẏ- (18)
Note that the normal component of the velocity can be expressed in terms of the gap function g(x, y) between the mass and the obstacle, equal to

g(x, y) = ax + b -y √ 1 + a 2 (19) This relation is v n = ( grad g• v) n/|| n||. With e = 1, N becomes N = 1 1 + a 2 1 -a 2 2a 2a a 2 -1 (20)
The eigenvalues of N are +1 and -1 (with corresponding eigenvectors [1/a, 1] and [-a, 1]). Additionally, N is symmetric and involutory (N 2 = I) hence is an orthogonal symmetry. The matrix of change of basis which diagonalises N does not affect the diagonal nature of the mass matrix, provided its diagonal elements are equal to each other-condition which is fulfilled here. It is therefore possible to write the initial dynamical problem as an equivalent problem corresponding to a 2-dof spring-mass system in chain (as in Fig. 1). More precisely, the equivalent system has the same mass matrix; its stiffness matrix is the conjugation of the initial stiffness matrix with

P = -a 1/a 1 1 (21)
and the new coordinates are given by P • [x, y] . The problem hence falls within the context of the mathematical assumptions and results given in Section 1.

Complex models of aircraft engines are not unidimensional; one node can be connected to multiple neighbouring nodes, breaking the tridiagonal structure of the stiffness matrix. An is illustrated in Figure 6. The stiffness matrix K is x y Figure 6. One mass connected to more than two other masses, breaking the tridiagonal structure of the stiffness matrix. no longer tridiagonal but remains symmetric. Also, since M is diagonal, M -1 K is also symmetric. It can therefore always be transformed in a tridiagonal matrix, leading to a spring-mass system in chain. The linear transformation affects the impact law but, from the previous result, the degrees-of-freedom of the contacting node can be rearranged to match with the generic system illustrated in Fig. 1.

The conclusion of this part is that the results extend to any spring-mass system-as long as its free flight exhibit linear dynamic behaviour. The key reason behind this is that any impact law written as an orthogonal symmetry with respect to a hyperplane can be diagonalised in an orthogonal basis where M -1 K remains symmetric positive definite.

Scalability to large systems

Turbomachine models may include thousands degrees-offreedom. It is natural to investigate how the results of Section 1 scale with the number of degrees of freedom. For motions with one impact per period, this is straightforward: for any T, expression (12) gives a potential solution which will be an actual solution if x n (t) ≤ d for all t. For two impacts per period, solutions only exist for specific times of impact t 1 , t 2 and an additional step is required to find such appropriate t 1 , t 2 . This consists in finding roots of the nonlinear equation ( 8 desktop PC are reported in Fig. 7 using the computer program Mathematica ® . The left diagram shows the computation time for the matrix reduction, i.e. the calculation of the different matrices involved (namely: M, K, B and Q). For a given number of degrees of freedom n, this is to be computed once. The right diagram shows the computation time of one coordinate in the state space, for a given time. Note that these timings would probably be reduced with a lower-level implementation and performance-tuning. Typical timing for 2000 dofs is 2 ms per coordinate. To approximate a nonlinear mode of a model with 2000 dofs with a few thousand states of the system would thus take a few seconds.

The conclusion of this subsection is that the number of degrees-of-freedom is not limiting for the computation of nonsmooth modes.

Sensitivity to number of dofs

It is natural to investigate the sensitivity of nonsmooth modes to the considered number of degrees-of-freedom. This is a challenging point: even for very large numbers of dofs, the existence of nonsmooth modes and their shape is extremely sensitive to n. An illustration is given in Fig. 8 where a mode exist in the neighbourhood of a period T when n = 100, while for the same T, the position calculated with formula (12) violates the constraint x n ≤ d with n = 101, proving no mode locally exists around T. The reason for this is that existence of nonsmooth modes is largely determined by the harmonics (see denominators of Eqs. (3), (10), and (11). The set of the periods of all harmonics {2kπ/ω i , k ∈ N * , i ∈ 1, n } becomes so dense that the existence of nonsmooth modes is ensured only on very tiny domains in the impact times t 1 , t 2 . As opposed to linear modes, nonsmooth modes with an impact law do not converge with the number of degrees-of-freedom. is therefore essential. In linear dynamics, the vibratory response to a periodic excitation is governed by the autonomous (i.e. unforced) response [START_REF] Meirovitch | Methods of analytical dynamics[END_REF]. Vibratory resonance regimes are observed if the forcing frequency lies in the vicinity of the natural frequencies of the system. However, the superposition principle no longer applies when impacts are accounted for, making the response of nonsmooth systems under periodic forcing much more challenging to understand. Additionally, it is not known whether periodically-forced nonsmooth oscillators always exhibit a (periodic) steady-state regime. This subsection aims at comparing the energy of the nonsmooth modal responses to the periodically forced solutions of the same system under a purely sinusoidal force F (t) = F 0 cos(Ωt) acting to every mass and a damping matrix simply chosen as C = αK.

Behaviour of forced and damped system

Using an event-driven algorithm [START_REF] Acary | Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics[END_REF], the time-evolution of the impact oscillator of Fig. 1 was calculated with n = 5 for several values of the excitation period T = 2π/Ω. After 500 such periods, it was observed that the system usually had a periodic motion with 1 impact per period T, 2 impacts per period T or 2T or 3 impacts per period 2T. For such periodic motions, the energy is plotted in Fig. 9. For the other ones, it is not known whether they could eventually reach the periodic regime.

It seems the peaks of energy for the forced systems (for T ≈ 2.4) are almost aligned with the 1 ipp energy curve. Nevertheless, this observation cannot be generalised to all resonances of the nonsmooth modes. It is therefore not possible to conclude that resonance phenomena of forced systems always coincide with that of the nonsmooth modes. The response for two applied forces and two damping coefficients are shown in Fig. 10. Further investigations are required to capture the variety of the observed behaviours in motion and in energy.

Exploratory stability analysis

The question of stability is a common question in the analysis of dynamical systems. Indeed, it characterises the dynamics of a system in the vicinity of its fixed points or periodic orbits. For periodic orbits, Lyapunov stability is by the fact that the response to any perturbation from a given periodic orbit remains in a neighbourhood or the periodic orbit [START_REF] Guckenheimer | Nonlinear oscillations, dynamical systems, and bifurcations of vector fields[END_REF]. Asymptotic stability is stronger and adds the condition that the system converges towards the orbit.

For smooth nonlinear systems, Lyapunov stability can be proven by linearising the dynamics in the neighbourhood of a known periodic solution. Similarly, here, the first points of the Poincaré section x n = d are plotted for a perturbed initial condition on the section. In other words, a point on the section corresponding to a nonsmooth mode (or more exactly, one periodic orbit on the manifold defining the mode) is selected. It is slightly perturbed (except the n-th position which is maintained equal to d) to provide an initial condition inserted in an event-driven algorithm to solve the dynamic equations: the dynamics is time-integrated and when x n = d, the sign of the velocity ẋn is inverted (Newton impact law). The time-integration then leads to a non-periodic trajectory. Every time the system hits the hyperplane x n = d in the state space, the pre-impact velocity ẋn and the position x n-1 are collected. The obtained results are illustrated in Fig. 11 for n = 5. First, it appears that the period T of the nonlinear mode plays an important role in the non-periodic motion with perturbed trajectories as depicted in Fig. 11, top. More than half of the 5000 free flights have a duration in [0.9T, 1.1T]. Second, the 5000 points of the Poincaré map appear to evolve in a fixed neighbourhood of the initial conditions. There is no tendency to move away from the neighbourhood with time: red and orange points are not particularly further from the center of Fig. 11 (bottom) than blue points. This behaviour has been observed for several numerical experiments, with various perturbations and n.

Although not a proof, these observations tend to indicate that orbits of nonsmooth modes are Lyapunov-stable, and not asymptotically stable.

CONCLUSION

The main mathematical results regarding nonsmooth modes, proven in a publication to come, are summarized. They prove existence and uniqueness, under some conditions, of periodic motions of an autonomous spring-mass chain, the last mass of which is subject to a Newton impact law. They also give the analytical expression of the solutions, when they exist and allow for an analytical calculation of nonsmooth modes with 1 impact per period, and approximation of nonsmooth modes for multiple impacts per period.

The goal of the present work is to identify possibilities and limitations to apply these results to the common blade-casing contact problem within turbomachines. The existence of the fundamental skew-symmetric matrix Π, governing nonsmooth modes, was shown to be still ensured for non-diagonal mass matrices, opening doors to FEM applications. Also, the unidimensional geometry of the spring-mass chain proved not to be required for the results to hold, making the study of complex geometries possible, the limitation being that the impact law must have the form of a orthogonal symmetry with respect to a hyperplane. The computation time of a coordinate is about 2 ms for 2000 dofs without performance-tuning and on a standard PC, which is very reasonable and makes the study of complex turbomachines models conceivable. First numerical experiments tend to show that nonsmooth orbits are Lyapunov stable and not asymptotically stable.

Altogether, this investigation shows that the gap between the theoretical results for the specific spring-mass chain and real applications to turbomachines is not as wide as initially expected. However, some limitations remain. The first one is the possibly tremendous difference in shape between nonsmooth modes of two similar systems differing by a few degrees-offreedom. This is explained by the extremely rich frequency range excited by impact events, and the underlying lack of convergence in space discretization. This is not an issue when the number of dof is fixed, as for example in a reduced-order model of a turbomachine where only specific modes of the complex structure are considered. The second current limitation is that the relationship between the energy-frequency graph of the forced system and its nonsmooth modes are not fully understood yet. Nonsmooth modal analysis can therefore not yet be used to predict the forced response of turbomachines and the forcing frequencies leading to high energy motions, when the forcing amplitude is too large.

1 Figure 2 .Figure 3 .

 123 Figure 2. Numerically calculated roots of the numerator of Q T (T -2t 1 ). In red: continuum of solutions corresponding to Fig. 3 and 4. Dashed line: T = 5T 2 where T 2 = 2π/ω 2 is the second eigenfrequency.

Figure 4 .

 4 Figure 4. Example of nonsmooth mode corresponding to (t 1 , T ) on the red curve of Fig. 2. The green curve corresponds to an periodic solution of the nonsmooth dynamics. The red curve highlights the velocity discontinuities occurring when x 2 = d.

Figure 5 .

 5 Figure 5. One mass with two degrees of freedom, subject to an elastic impact law.

Figure 7 .

 7 Figure 7. Computation time for construction of matrices of Section 1.2 and calculation of coordinates in phase space. [ ] 1 ipp. [ ] 2 ipp.

Figure 8 .

 8 Figure 8. Non-convergence with respect to spatial discretization. Curves show the periodic motion of the n-th mass (1 ipp). For n = 100 [ ], the position is an actual solution. For n = 101 [ ], the position is completely different and is not even a solution: it violates the constraint x n ≤ d.

Figure 9 .Figure

 9 Figure 9. Energies of nonsmooth modes for 1 ipp [ ] and 2 ipp [ ]. Energy of periodic solutions for the forced system with F 0 = 1, α = 0.005 [ ] and α = 0.01 [ ].

Figure 11 .

 11 Figure 11. Motions with perturbated initial conditions. Top: Duration of the 5000 first free flights, normalised. Bottom: Poincaré section x n = d.


Hao Hu
email: hao.hu@ensam.eu

Mathias Kleiner
email: mathias.kleiner@ensam.eu

Jean-Philippe Pernot
email: jean-philippe.pernot@ensam.eu

ON THE DETECTION OF OVER-CONSTRAINED SUBPARTS OF CONFIGURATIONS WHEN DEFORMING FREE-FORM CURVES

Keywords: Free-form curves, locally unsatisfiable or redundant configurations, Dulmage-Mendelsohn decomposition, Jacobian matrix analysis

Today, designers use CAD modelers to define and modify NURBS surfaces involved in the design of complex shapes like car bodies or turbine blades. The generated shapes often result from the use of variational modeling techniques where user-specified constraints define the shapes. However, for free-form curve/surfaces, if too much constraints are added to subparts of a configuration, the system will not be solvable even if it is globally well-/underconstrained. When this happens, it is useful to identify locally unsatisfiable subparts of configurations and provide the user feedback for adjustment. Currently, in the domain of geometric constraint solving, techniques are mainly developed for Euler geometries rather than parametric entities like free-form curves/surfaces. In this paper, we apply the Dulmage-Mendelsohn decomposition method to isolate structural over-constrained subparts of configurations. Since structural overconstraints do not necessarily mean unsatisfiable, a Jacobian matrix analysis approach is taken to further detect the inconsistent constraints. Indeed, these numerical methods can be generalized to detect overconstraints on free-form curves. We illustrate our approach on different examples where results show that Gauss elimination, though restricted to linear cases, is more relevant in our context than Dulmage-Mendelsohn decomposition.

INTRODUCTION

With the evolution of manufacturing technologies and advances in the domain of new materials, today's CAD modelers need to develop high quality and complex shapes of products. Free-form shapes, especially in the consumer products industry, not only need to meet aesthetic criteria but also satisfy functional criteria.

The shape of a dashboard, a wing or a turbine blade are examples of complex shapes obtained by modifying NURBS surfaces to satisfy user-specified constraints [START_REF] Pernot | Fully Free Form Deformation Features for Aesthetic and Engineering Designs[END_REF]. To satisfy these requirements, designers are relying on computer aided design (CAD) software.

Free-form curves are often designed using deformation through position, distance tangency, and/or continuity constraints. However, unlike 2D sketching where inconsistent constraints can be identified interactively during the design process, there is no indicator to analyze the constrained status of free-form curves. This is partly due to the absence of suitable techniques for solving algebraic problems that underlie more complicated specifications. In this paper, we seek to extend the shape vocabulary of geometric constraint solvers to free-form class by applying methods to analyze the structure of system of equations.

Therefore, our goal is to provide a diagnosis to the user through the identification of inconsistent constraints, subparts, and corresponding explanations.

Plan of the article

Section 2 summarizes over-constraints detection methods from the geometric constraint solving domain. In section 3, we manually apply a degree of freedom analysis on different B-spline curves in order to introduce the examples used in the following sections. In section 4, we illustrate the use of the Dulmage-Mendelsohn (D-M) decomposition algorithm to identify over-constrained parts of the corresponding equation systems. In section 5, we show how a Gauss elimination can be used in the linear case to further discriminate inconsistent and redundant constraints. Finally, we discuss about future work.

RELATED WORK

Selection of geometric constraint methods for overconstraints detection

An over-constrained system is a system which contains redundant or conflicting constraints. We propose to classify the existing constraint solving methods into two categories with respect to the type of over-constraints detected. Graph based decomposition allow to identify structural overconstraints. However these methods will not detect constraints that are numerically inconsistent with others in the system [START_REF] Michelucci | Geometric constraint solving:the witness configuration method[END_REF]. Algebraic methods, such as Gauss elimination or Grobner bases, thus generalize the notion of over-constraints. However time complexity is exponential in the non-linear case, thus limiting their use to small systems [START_REF] Cox | Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra[END_REF]. For a general constraint solver, common practice is to first exploit geometric domain knowledge to transform a constraint system into a set of small solvable subsystems and then compute and assemble the solutions of its subsystems. Our presentation in the following part is intended to reflect schemes that are more suited to isolate over-constraints within the two phases.

Geometric constraint system decomposition

Jermann and Trombettoni [START_REF] Jerman | Decomposition of geometric constraint systems: A survey[END_REF] have proposed to classify the existing decomposition approaches into 4 categories according to the way they operate.

1) The recursive division methods

These methods, first introduced by Owen [START_REF] Owen | Algebraic solution for geometry from dimensional constraints[END_REF], work by recursively splitting the constraint system into subgraphs. It handles 2D constraint systems formed by points and lines linked by distance and angle constraints. These methods are used to detect whether a system is solvable, but will not provide explanations on which constraints cause the inconsistencies.

2) The recursive assembly methods These methods adopt a bottom-up scheme by firstly finding rigid subsystems, called cluster and then iteratively aggregate rigid components into bigger ones [START_REF] Hoffmann | Decomposition plans for geometric constraint problems, Part II. New Algorithms[END_REF]. This allows to find minimal rigid (dense) subgraphs. Later, Hoffman & al. [START_REF] Hoffmann | Bo Yuan Making constraint solvers more usable: overconstraint problem[END_REF] extends the method to locate 1-overconstrained dense subgraphs by modifying the incremental network maximum flow algorithm. Although this is sufficient when incrementally placing constraints during the design, it cannot be used when combining or editing constrained designs.

3) The single-pass methods

These methods decompose a system in a single iteration where all the components are produced simultaneously. There are two kinds of methods: maximum matching and weighted maximum matching. Serrano [START_REF] Serrano | D Automatic dimensioning in design for manufacturing[END_REF] describes the detection of over-constraints by finding unmatched constraints after applying maximum matching to the constraint system, where all constraints and entities have one degree of freedom. Also, D-M decomposition algorithm can decompose an equation system into over-,well-and under-constrained subsystem based on maximum matching [START_REF] Dulmage | On the inversion of sparse matrices[END_REF]. Latham and Middleditch extends the work of Serrano by proposing weighted maximum matching to identify overconstraints with an arbitrary number of degrees of freedom [START_REF] Latham R S | Connectively analysis: a tool for processing geometric constraints[END_REF]. It is performed directly on the constraint graph of the geometric system, avoiding the need to translate it into an equation system.

4) The propagation of degree of freedom approaches(PDOF)

According to Trombettoni [START_REF] Trombettoni | GPDOF: a fast algorithm to decompose under-constrained geometric constraint systems: Application to 3D model reconstruction[END_REF], these approaches are more suitable for under-constrained systems. Redundant equations must be detected and removed before launching PDOF. Therefore, they can not be used for overconstraints detection.

Algebraic methods for detecting inconsistent system of equations

In general, almost all of geometric constraints can be translated mechanically into a set of non-linear equations. The equations are usually algebraic, and non-algebraic formulations involving trigonometric functions can be avoided in nearly all cases [START_REF] Hoffmann | Geometric constraint decomposition[END_REF]. Therefore, at the equation level, detecting geometric overconstraints turns into identifying a set of inconsistent or redundant equations among a system of equations.

1) Numerical methods

Equations are said independent if the Jacobian matrix has a full rank at the searched root. For linear systems, the Jacobian matrix is composed of coefficients of equation systems, which is independent of the roots. Light and Gossard use Gauss Elimination(GE) to compute the rank as well as further identify invalid equations [START_REF] Light | Variational geometry: A new method for modifying part geometry for finite element analysis[END_REF]. Serrano extends this work to check existence of overconstraints within strong connected components of equation systems [START_REF] Serrano | Constraint Management in Conceptual Design[END_REF]. However, such roots are hard to find in non-linear systems, therefore computing the Jacobian matrix is a difficult task. In rigidity theory, classical Numerical Probabilistic Method(NPM) studies the structure of the Jacobian at a random configuration [START_REF] Graver | Combinatorial Rigidity[END_REF]. Michelucci extends the NPM to study Jacobian structure at a configuration satisfying the incidence constraints (collinearities and coplanarities) [START_REF] Foufou | Numerical decomposition of geometric constraints[END_REF].Based on the extended NPM, he then develops the Witness Configuration Method(WCM) to identify subtle dependences between geometric constraints [START_REF] Michelucci | Interrogating witness for geometric constraint solving[END_REF].

2) Symbolic methods

Thanks to the consistency algorithm from algebraic geometry [START_REF] Cox | Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra[END_REF],

for polynomials

, , [ , , ]

s n f f x x ⋅⋅⋅ ∈ ⋅⋅⋅ �
,we compute a reduced Grobner basis(GB) of the ideal with respect to any ordering. If this basis is {1}, then the polynomials have no common zero in � ; otherwise, they must have a common zero. In theory, this algorithm can be used to identify all the over-constraints, but the computation cost is exponential and thus only limited to small systems. Other symbolic methods like Wuritt as well as theory of determinants(TD) are also able to characterize over-constrained sets but have the same practical limitation as Grobner base.

Constraint solving in terms of curves D.Podgorelec extends geometric elements to Bezier curves and ellipses by mapping them into auxiliary lines and points and uses local propagation to solve the resulting constraint system [START_REF] Podgorelec | Dealing with redundancy and inconsistency in constructive geometric constraint solving[END_REF]. Y. J. Ahn implements the construction of quadratic Bezier curves subject to tangency and length or energy minimization constraints , which offers the ability to find curves that are to have prescribed clearance from points, circles or straight-line borders [START_REF] Ahn | Constraintbased LN-curves[END_REF]. Iddo Hanniel & al. use symbolic polynomials to represent inequality expressions defined by conditions on curves and apply inequality constraints to coefficients of symbolic polynomials [START_REF] Hanniel | Solving global geometric constraints on free-form curves[END_REF]. These methods are focusing on solving specific issues of CSP(constraint satisfaction problems) on curves, thus lacking the generality. To the best of our knowledge, general methods have not emerged yet.

Dulmage-Mendelsohn Decomposition

Dulmage-Mendelsohn decomposition algorithm allows to decompose a system (based on bipartite graph) into over-constrained, well-constrained and under-constrained part [START_REF] Pothen | Computing the Block Triangular Form of a Sparse Matrix[END_REF]. The well-constrained part can be further induced into irreducible subparts.

Here we extend this method to analysis system of constraints represented in equation level.

The decomposition has the following properties:

1) G3 is structural over-constrained part, G1 is structural well-constrained part, G2 is structural under-constrained part

2) The unmatched equations(UE) are either conflicting or redundant 3) G1 , G2 and G3 are independent of the choice of the maximum matching Since here the numerical information of equation system is ignored, structural over-constrained part is subset of over-constrained parts of a system, that is to say, except G3, it is also possible for G1 to have overconstraints. Algebraic methods are proposed to further analysis the irreducible subparts contained in G3. For a discussion on this, the reader is referred to the work of Serrano [START_REF] Serrano | Constraint Management in Conceptual Design[END_REF]. In section 4.4, we will give such examples on B-splines.

Jacobian structure of linear constraints

The Jacobian matrix

The Jacobian matrix is a matrix (n*m) containing the partial derivatives of each constraint equation with respect to each degree of freedom. It is usually defined and arranged as following:

D-M decomposition Figure 1 1 1 1 1 2 1 2 n m m m n f f f x x x J f f f x x x ∂ ∂ ∂     ∂ ∂ ∂     =       ∂ ∂ ∂   ∂ ∂ ∂   ��� � � � � ��� (1)
Where th i constraint is in the form:

( ) 1 2
, 0

i j f x x x = ��� (2)
Identify inconsistent equations Equations ( 2) are independent if the Jacobian has full rank at the searched root. For a linear constraint system, Jacobian matrix is the coefficient matrix of the equation system and by using Gauss Elimination methods, we can check the rank as well as further analyze its structure.

Constructing the Jacobian of (2), the following system is obtained:

1 x x J f - ′ = - (3)
Reformulate equation ( 3), we have:

J x f ∆ =- (4)
Where x x x ′ ∆ = -.

Use the transformation R and C, equation ( 4) becomes:

T RJCC x Rf ∆ =- (5)
The structure of equation ( 5) is as following:

(6) Equation ( 7) represents the th i row of equation ( 6), where subscript i = r+1…n, indicating a specific row of the matrix.

( ) ( )

T i i i RJC C x R f ∆ = - (7)
Since these rows of the product RJC are all zeros, in order to have a solution for ( )

T i C x
∆ , the following must be satisfied:

0 i R f -= i=r+1…n (8)
Hence, the non-zero terms of the th i row in equation ( 8) indicate an inconsistent set of constraint equations.

TEST CASES

In this section, we consider constraints satisfaction problems only on curves since it allows us to illustrate our problems clearly.

Simple Curves

In this section, we set up three constraint satisfaction problems on B-splines with control points as variables. For the three configurations, constraints of curve 1 are the same while for curve 2, constraints are added incrementally. The objective is to minimize the deformation of the two curves, while satisfying the constraints. However, none of the [START_REF] Piegl | The NURBS book[END_REF], these examples will be used to illustrate several key concepts.

Complex curves

Here is a B-spline curve with 13 control points(Fig 7). 10 of them (green ones) are free to move in order to satisfy 9 position constraints (p1-p9). Again, we want to minimize the deformation of the curve while satisfying the 9 constraints. Here again, the constraints can not be satisfied.

Degree of freedom analysis based on local support property

We apply degree of freedom(DOF) analysis on these examples manually based on local support property.

DOF analysis of configuration two-configuration four

According to the local support property of B-splines, constrained points on curve 1 and curve 2 are influenced by certain area of control points (Table 1). Each free control point has a DOF of 2 and position constraints constrains 1 degree of freedom. Table 2 summarizes the DOF analysis results. The curve is divided into 10 segments (Seg1-Seg10). End points of each segment are calculated by P(Ui) and P(Ui+1), where Ui and Ui+1 are knots. Analysis result is shown in Table 3.

DOF analysis of configuration six

From Table 2 to Table 3, we know that a constraint system can be locally over-constrained but globally be under-, well-and over-constrained. Since all these constraint systems have no solution, an over-, underand well-constrained system described from degree d d - Affected by

1 3 1 1 d d - 1 3 1 1 d d - 0 2 2 2 d d - 1 3 2 2 d d -
0 0 0 3 d d - 0 0 1 4 d d - 0 0 2 5 d d - 0 0 3 6 d d - 0 0 4 7 d d - 0 0 5 8 d d - 0 0 6 9 d d - 0 0 7 10 d d - 0 0 8 11 d d - 0 0 9 12 d d - Dof 2*

IDENTIFY LOCAL UNSATISFIABLE PARTS OF CONFIGURATIONS

In this section, we try to automatically identify additional overconstraints in the previous configurations. As our initial trial, identification process is conducted at the level of equations.

Equation-based representation of the constraint systems

Representation of constraints for configuration two 0.4(0. + 0.6x1) + 0.6(0.7x1 + 0.3x2) = 1 (e1) 0.4(0. + 0.6y1) + 0.6(0.7y1 + 0.3y2) = 1 (e2)

(0.5x1 + 0.5x2) = 1.5 (e3)
••• (e18)
Where e1 and e2 correspond to p1 position constraint, e3 and e4 correspond to p2 position constraint. For the other constraints, we use ei to represent them orderly.

For brevity, only configuration two and configuration six are taken as examples for this method of identification.

Bipartite graph representation of the constraint systems

Representation of constraints for configuration two

D-M decomposition of the systems

Decomposition of configuration two Fig. 11 shows that the curve 1 is over-constrained which is consistent with the result shown in Table 2.

For configuration six, Fig. 13 shows that the constraint system is divided into four parts: two under-constrained parts, one over-constrained part and the rest is well constrained. The result matches well with Table 3.

Limitation of applying D-M decomposition algorithm

The over-/under-/well-constrained sub-systems decomposed by D-M algorithm is structurally over-/under-/well-constrained. As stated in section 2.2, it is also possible for well-constrained sub-systems to have overconstraints. Here we show such examples on B-splines.

On the segment 3 of configuration seven, the constrained point is restricted to two different positions at the same time(conflicting constraints:p3 and p4); On the segment 9 of configuration eight, the constrained point is restricted to the same positions twice(redundant constraints:p8 and p9) From Fig. 16, we can know that even if we have modified several constraints in configuration seven and eight, results after D-M decomposition are the same: Seg 3 in configuration seven and Seg 9 in configuration eight are still well-constrained! This is contradictory to our modification in both configurations. In order to accurately find all the overconstraints, we resort to the algebraic methods.

DETECTION USING JACOBIAN MATRIX ANALYSIS

To further detect over-constraints in these two configurations (including both conflicting and redundant constraints), we apply Gauss elimination on their corresponding equation systems. As illustrated in Fig. 17 and Fig. 18, conflicting

CONCLUSION AND FUTURE WORK

In this paper we have tried to evaluate different methods for detecting over-constrained subparts of free-form based geometric configurations. First, considering the local support property of free-form curves, we have applied a DOF analysis and shown that it is not sufficient to detect all inconsistencies. Then, we have applied D-M decomposition to the corresponding equations systems in order to automatically detect structural over-constraints. Since these methods ignore the numerical information, inconsistent constraints may still be undetected. Therefore, we have shown that a Jacobian Matrix analysis can be used to fully discriminate conflicting and redundant constraint. However, the time complexity of the numerical analysis restricts its use to small cases when facing non-linear equations.

In the geometric constraint solving domain, a 2-step approach is often used to circumvent this, where a graph-based method is first realized to decompose the system into subparts on which algebraic solving methods are applied. In future work, we plan to adapt these approaches to the detection of inconsistent and redundant constraints.

Graph methods such as MFA are an efficient way of decomposing into irreducible subsystems [START_REF] Hoffmann | Finding solvable subsets of constraint graphs[END_REF]. However, in terms of constraint solving on free-form curves, as far as we know, there is no graphreduction methods that are able to structurally decompose a system while considering its geometric information. A perspective could be to develop a graph-based decomposition that takes into account the local support property. Combined with an algebraic method such as witness configuration [START_REF] Thierry | [END_REF][START_REF] Mireille | Defining tools to address overconstrained geometric problems in Computer Aided Design[END_REF], this could provide a general and reliable approach for the identification of overconstraints.

Figure 2 :Figure 3 :

 23 Figure 2: Configuration One

Figure 6 :

 6 Figure 6: Configuration five-a B-spline curve

Figure 4 :Figure 5 :

 45 Figure 4: Configuration Three

Figure 7 :

 7 Figure 7: Configuration six-The curve with 9 position constraints

  (0.5y1 + 0.5y2) = 1.25 (e4) 0.4 (1.2 + 0.6 x2) + 0.6 (0.3 x1 + 0.7 x2) = 2.0 (e5) 0.4 (0. + 0.6 y2) + 0.6 (0.3 y1 + 0.7 y2) = 1.0 we use equations to represent the 9 position constraints. They are: 0.125 x1 + 0.59 x2 + 0.26 x3 + 0.02 x4 = 1.0 (e1) 0.125 y1 + 0.59 y2 + 0.26 y3 + 0.02 y4 = 1.5 (e2) 0.085 x3 + 0.63 x4 + 0.28 x5 = 3.5 (e3) 0.085 y3 + 0.63 y4 + 0.28 y5 = -0.5 (e4)

Figure 8 :

 8 Figure 8: Bipartite graph of configuration two Representation of constraints for configuration six

Figure 9 :

 9 Figure 9: Bipartite graph representation of configuration sixBased on the equations described in section 4.1, bipartite graph representation of both systems is shown in Fig.8and Fig.9, where vertices are used to represent variables and equations while edges indicate their relationship.

Figure 10 :Figure 11 :

 1011 Figure 10: Maximum matching of configuration two

Figure 13 :

 13 Figure 13: Associated directed graph of overconstrained/under-constrained sub-graph

Figure 12 :Figure 14 :

 1214 Figure 12: Maximum matching of configuration six

Figure 15 :Figure 16 :

 1516 Figure 15: Configuration eight-a b-spline with 8 constraints

Figure 17 :

 17 Figure 17: Analysis result of configuration seven

Table 1 :

 1 Area of control points that influence the constrained points (configure twoconfigure four)

		Configuration two	Configuration three	Configuration four
		Curve 1	Curve 2	Curve 1	Curve 2	Curve 1	Curve 2
	DOF	4	4	4	4	4	4
	DOC	6	0	6	2	6	4
	Status(local)	Over	Under	Over	Under	Over	Well
	Status(global)	Under		Well		Over

Table 2 :

 2 Constrained status of 10 segments (over/under/well is short for over-/under-/well-constrained)

	Seg1	Seg2	Seg3	Seg4	Seg5	Seg6	Seg7	Seg8	Seg9	Seg10

Table 3 :

 3 Constrained status of 10 segments (over/under/well is short for over-/under-/well-constrained) of freedom point of view does not reflect its solvability.

		4	2*4	2*3	2*3	2*3	2*3	2*4	2*4	2*3	2*2
	Doc	2*1	0	2*3	0	0	0	2*1	0	2*4	0
	Status	Under	Under	Well	Under	Under	Under	Under	Under	Over	Under


Stanislas Rigal

Gilles Cuny

On the rarity of anterior teeth of Asteracanthus magnus (Euselachii; Hybodontiformes)

Keywords: Asteracanthus, dentition, Jurassic, Bathonian, France

A partial dentition found in the Middle Jurassic (Bathonian) of Normandy (France) associates lateral teeth of Asteracanthus magnus with anterior teeth usually attributed to A. tenuis or A. longidens. According to this French specimen, teeth usually identified as anterior of A. magnus represent in fact teeth of the first lateral file. This discovery explains the apparent rarity of isolated anterior teeth of A. magnus in the fossil record, for which no associated dentition were hitherto known. It also demonstrates that the general arrangement of the dentition of A. magnus is closer to the one of A. smithwoodwardi than to the one of A. medius.

Introduction

The genus Asteracanthus was erected by Agassiz in 1837 based on dorsal fin-spines found in the Upper Jurassic of England that he attributed to A. ornatissimus (AGASSIZ 1833-44, vol. 3, p. 31). In the same work, this author also erected the genus Strophodus, which he thought was closely related to Acrodus, based on a partial dentition found in the Bathonian of France, which he attributed to S. longidens (AGASSIZ 1833-44, vol. 3, p. 117). The latter was, however, destroyed during the bombing of the Caen Museum by the Allied in July 1944 [START_REF] Bigot | La destruction des collections et des bibliothèques scientifiques de Caen[END_REF]. The discovery of teeth of Strophodus associated with fin-spines of the type species of Asteracanthus, A. ornatissimus, represented the first piece of evidence that the teeth and finspines could represent in fact a single animal [START_REF] Woodward | On some remains of the extinct Selachian Asteracanthus from the Oxford Clay of Peterborough, preserved in the collection of Alfed N. Leeds, Esq., of Eyebury[END_REF]. The apparent synonymy between the two genera is, however, not so simple as recently demonstrated by [START_REF] Rees | Hybodont sharks of the English Bathonian and Callovian (Middle Jurassic)[END_REF], who erected a new genus, Frangerodus, for the teeth of Strophodus lingualis. No dorsal fin-spines have so far been associated with this new genus.

Frangerodus appears to be less heterodont than Asteracanthus, but the heterodonty pattern of the latter genus is not always fully understood. It has been suggested, for example, that one species of Asteracanthus, A. magnus, a species known from the Bathonian of England and France, might have possessed a single file of anterior teeth although most Jurassic species possess two files [START_REF] Rees | Hybodont sharks of the English Bathonian and Callovian (Middle Jurassic)[END_REF]. We report here a partial dentition of A. magnus from the Bathonian of Normandy, which allows a better understanding of the heterodonty pattern in this species.

Material and methods

The described sample (MPV.2010.3.44) shows 54 teeth included into a piece of Bathonian Calcaire de Caen measuring 19 x 12 x 11 cm (Fig. 1). MPV.2010.3.44 was purchased by the Paléospace, palaeontological museum of Villers-sur-Mer (MPV) during the auction of the Folet collection, a local avocational palaeontologist, at Saint Valery en Caux (Normandy) on January 3 rd , 2010. It was extracted from a quarry in the vicinity of the city of Caen, although the exact location was not provided and remains unknown as many quarries in the area have been closed down and are no longer accessible. Insert Fig. 1 near here Although the dentition is not fully articulated, several groups of teeth are in connection. Together with the fact that the porous and fragile roots of the teeth are perfectly preserved, it indicates that these teeth belonged to a single individual.

All measurements were taken with callipers having an accuracy of 0.1 mm. Pictures of the specimen were taken with a digital camera (Canon G11) in macro mode. The systematic scheme used here follows REES & UNDERWOOD (2008) and descriptive terminology follows CAPPETTA (2012), more particularly concerning the definition of "rows" and "files". The use of the term "root" also follows CAPPETTA ( 2012), but does not imply direct homology with the root of the teeth of other gnathostomes [START_REF] Underwood | Development and evolution of dentition pattern and tooth order in the skates and rays (Batoidea; Chondrichthyes)[END_REF].

Geological settings

The Calcaire de Caen Formation is dated from the beginning of the Middle Bathonian (DUGUÉ et al. 1998). It is made of an homogeneous, fine-grained biomicrite, and its environment of deposition is considered as a low-energy, subtidal one (DUGUÉ et al. 1998). This formation is well-known for its fossils of crocodiles and dinosaurs that have been described during the 19 th Century by Cuvier, Geoffroy-Saint-Hilaire and Eudes-Deslongchamps (BUFFETAUT 2008).

Systematic palaeontology

Class Chondrichthyes HUXLEY, 1880

Cohort Euselachii HAY, 1902 Order Hybodontiformes [START_REF] Patterson | British Wealden sharks[END_REF] Family Hybodontidae OWEN, 1846 Subfamily Acrodontinae CASIER, 1959 sensu [START_REF] Maisey | Hamiltonichthys mapesi, g. and sp. nov. (Chondrichthyes; Elasmobranchii), from the Upper Pennsylvanian of Kansas[END_REF] Genus Asteracanthus AGASSIZ, 1837 Asteracanthus magnus (AGASSIZ, 1838) Description: The dentition MPV.2010.3.44 includes 54 teeth, which have been subdivided into six files: symphyseal, first anterior, second anterior, first lateral, second lateral and posterior ones (Fig. 2). Anterior and symphyseal teeth are arched in labial or lingual view and show enlarged foramina concentrating in the centre of the lower part of the root in lingual view, whereas lateral teeth are flattened with roots less deep than in the anterior teeth. Insert Fig. 2 near here Four teeth identified as potential symphyseal ones are strongly arched in lingual view (Fig. 2.1). The crown is symmetric mesio-distally and shows a reticulated ornamentation as well as a longitudinal crest, which is worn out at the apex of the crown. The mesial, distal, labial and lingual faces are ornamented with irregular vertical ridges, more or less developed. The crown is separated from the root by a shallow groove and overhangs it. The root is at least twice the height of the crown. It shows randomly distributed foramina, with larger ones concentrating in the centre of its lower part in lingual view.

Nine teeth identified as belonging to the first anterior file are also strongly arched in labial or lingual view, but their crowns are asymmetric, the distal part being more elongated and tapered than the mesial one (Fig. 2.2). The distal part is also curved lingually. The ornamentation is reticulated and a longitudinal crest is present, although it tends to disappear on the mesial part of the crown. The ornamentation of the mesial, distal, labial and lingual faces is similar to that of the symphyseal teeth. The crown is separated from the root by a shallow groove and overhangs it. The root has a structure similar to that of the symphyseal tooth.

Twelve teeth identified as belonging to the second anterior file are less arched in labial or lingual views than those of the first anterior file (Fig. 2.3). Apart from that, their asymmetric crown is quite similar to the one of the teeth of the first anterior file. Two teeth in connection show that this shape allowed a tight interlocking between teeth of the same family. The ornamentation of the crown is reticulated, but the longitudinal crest is either absent or restricted to the distal part of the crown. The ornamentation of the mesial, distal, labial and lingual faces is similar to that of the preceding teeth. The crown is separated from the root by a shallow groove and overhangs it. In labial view, the root is divided in two halves by a vertical crest. Large foramina tend to concentrate on its distal half.

Nine teeth identified as belonging to the first lateral file are twice as large as the anterior ones (Table 1, fig. 2.4). In apical view, the crown shows a labial projection so that its outline is nearly triangular. Its ornamentation is reticulated, without a longitudinal crest. The base of the crown, as well as that of the root is nearly flat. The crown is separated from the root by a groove and overhangs it. The root is higher than the crown, but contrary to what is observed in anterior teeth, it never exceeds twice the height of the crown. The root, like the crown, shows a labial projection, better developed in its upper than in its lower part. Large foramina are scattered in the lower part of the root.

Insert Table 1 near here

Eleven teeth identified as belonging to the second lateral file are twice as large as those of the first lateral one (Fig. 2.5). The crown has a rhomboid, almost rectangular, outline in apical view. The ornamentation is reticulated, and often quite worn out on the apical face, whereas the mesial, labial, distal and lingual ones are smooth. Teeth identified as belonging to the lower dentition, see below, are narrower labio-lingually than the ones belonging to the upper dentition (Table 1). The crown is separated from the root by a groove. The root has a rather quadrangular shape and is projected distally. Large foramina are roughly aligned in the middle part of the root.

Nine teeth identified as belonging to the posterior file show a flat, oval-shaped crown, which is slightly compressed mesio-distally (Fig. 2.6). Its ornamentation is coarser than on the anterior and lateral teeth and is made of numerous small anastomosing ridges. The mesial, labial, distal and lingual sides of the crown are smooth. The root is quadrangular in shape with a flat base and is projected lingually. It is at least twice as high as the crown and twice as long labio-lingually. The root is also projected distally. In apical view, the distal projection is separated from the crown by a groove. The labial and lingual faces of the root are oblique so that there is an efficient interlocking between teeth of the same family as demonstrated by teeth still in connection. The root is devoid of enlarged foramina.

Discussion

Reconstruction of the dentition

Two sets of two second lateral teeth in connection, one with a slightly concave surface (Fig. 3, LUSL), the other with a slightly convex surface (Fig. 3, LLSL), appear to have been in occlusion, indicating that MPV.2010.3.44 encompasses teeth belonging to the lower and upper jaws. If we compare with Agassiz"s description of Asteracanthus (Strophodus) longidens, the convex teeth are likely to belong to the lower jaw, whereas the slightly concave ones belong to the upper jaw (AGASSIZ1833-44, vol. 3, p. 122). The two upper second lateral teeth appear to be in connection with two first lateral (Fig. 3, LUFL), two second anterior (Fig. 3, LUSA) and one first anterior ones (Fig. 3, LUFA). These seven teeth belong to the left part of the upper dentition. Its right part is probably represented by four second lateral teeth, one devoid of crown and one with the crown partly lacking, located on the other side of MPV.2010.3.44 (Fig. 3, RUSL). All their crowns are slightly concave and they appear to belong to the same file. At least four teeth in a file are congruent with the seven teeth per file observed in A. medius by OWEN (1869). These four teeth are still in connection with a first lateral tooth (Fig. 3, RUFL). Five posterior teeth, three of which still in connection, located in the middle part of MPV.2010.3.44 between the left and right upper part of the dentition are therefore likely to belong to the lower jaws (Fig. 3, LP), the lower dentition having slid on the upper one. Apparently, the lower dentition was quite scattered around during the sliding, so that no connection between teeth of different files was preserved. As there is a mixing between lower and upper teeth in the central part of MPV.2010.3.44, it is quite difficult to ascertain whether the anterior and lateral teeth present in this area belong to the lower or upper jaws. The same is true of the symphyseal teeth recovered (Fig. 3,S), which are also located in this area.

Insert Fig. 3 near here When describing the lower jaw of Asteracantus medius, OWEN (1869) noticed that it did not possess a symphyseal file, like in the modern Heterodontus. However, comparisons with this genus seems to be quite hazardous, as different species might have various dental patterns. Indeed, several authors attributed a symphyseal file either on only one jaw (lower jaw of Heterodontus portusjacksoni in PFEIL ( 2010)) or on both jaws [START_REF] Reif | Morphogenesis, pattern formation and function of the dentition of Heterodontus[END_REF]. Woodward, on the other hand, associated the symphyseal file of Asteracanthus ornatissimus with the lower jaw (WOODWARD 1888), as did PFEIL (2010). MPV.2010.3.44 does not allow deciding on which jaw it is the most likely that symphyseal teeth were present.

The distal projection of the root of the posterior teeth seems to indicate that MPV.2010.3.44 possessed a single file of posterior teeth, at least on the lower jaw, like the Liassic species A. smithwoodwardi from Switzerland [START_REF] Peyer | Die schweizerischen Funde von Asteracanthus (Strophodus)[END_REF] although most Middle Jurassic species possessed two files [START_REF] Rees | Hybodont sharks of the English Bathonian and Callovian (Middle Jurassic)[END_REF]. The general dental formula of MPV.2010.3.44 is therefore likely to be S-A1-A2-L1-L2-P for the lower jaw and A1-A2-L1-L2-P for the upper jaw.

Comparisons

Compared with A. tenuis and A. longidens, MPV.2010.3.44 possess less elongated lateral teeth, without a dome and its ornamentation pattern is very different from the one of A. ornatissimus [START_REF] Rees | Hybodont sharks of the English Bathonian and Callovian (Middle Jurassic)[END_REF]. On the other hand, the shape and ornamentation pattern of the teeth of MPV.2010. 3.44 belonging to the second lateral file are very similar to those of A. magnus (AGASSIZ 1833-44, REES & UNDERWOOD 2008). However, the anterior teeth of MPV.2010.3.44 are more arched than the ones usually attributed to the latter species and are more reminiscent of anterior teeth of A. tenuis or A. longidens [START_REF] Rees | Hybodont sharks of the English Bathonian and Callovian (Middle Jurassic)[END_REF]. The tooth P.66481 from the collection of the Natural History Museum, London, identified by REES & UNDERWOOD (2008, pl. 5, figs. 1-2) as an anterior tooth of A. magnus is strikingly similar to the teeth of MPV.2010.3.44 interpreted here as belonging to the first lateral file. They share the same ornamentation pattern and a flat crown with a triangular outline. First lateral teeth with a similar outline were also illustrated for A. subreticulatus by AGASSIZ (1833-44) and [START_REF] Oppenheimer | Der Malm der Schwedenschanze bei Brunn. -Beiträge für Paleontologie und Geologie des Österrreich und des Ungarns[END_REF]in A. smithwoodwardi by PEYER (1946). A. subreticulatus was considered a junior synonym of A. ornatissimus by WOODWARD (1889). We therefore conclude that P.6648 was misidentified as an anterior tooth and represents in fact a first lateral one, leaving the species A. magnus without identified anterior teeth. As there is a single file of first lateral teeth, and not two as for anterior teeth, this explains the apparent rarity of this kind of teeth as noted by [START_REF] Rees | Hybodont sharks of the English Bathonian and Callovian (Middle Jurassic)[END_REF] in the collection of the Natural History Museum in London.

REES & UNDERWOOD (2008) considered A. magnus and A. medius as belonging to the same group based on the following characters: "Teeth with a weak reticulate ornamentation and a dentition with wide, rectangular lateral teeth and weakly arched anterior teeth". MPV.2010.3.44 shows that in fact A. magnus is quite different from A. medius based on the following characters: presence of strongly arched anterior teeth, first lateral teeth much shorter mesio-distally than the second lateral ones (They have approximately the same length in A. medius) and the presence of a single file of posterior teeth. All in all, the general pattern of the dentition of A. magnus is much closer to A. smithwoodwardi than to A. medius. PFEIL ( 2010) described an Asteracanthus dentition of BSPG 2010 I 91 from the Tithonian of South-West Germany, where the symphyseal, anterior and first lateral teeth (symphyseal, anterior and anterio-lateral in PFEIL ( 2010)) are very similar to those of MPV.2010.3.44. However, the shape and the size of both second lateral (lateral in PFEIL ( 2010)) and posterior teeth are fairly different from those of MPV.2010.3.44. The distal edge of the second laterals of (BSPG 2010 I 91) has a blunt tip and are smaller and less elongated (mesio-distal length: 23 mm, labio-lingual width: 17 mm). The posterior teeth have a more elongated (mesio-distal length: 7 mm, labio-lingual width: 20 mm) and more angular shape.

Conclusion

MPV.2010.3.44 represents the first discovery of a dentition of A. magnus partly in connection and allows a complete reappraisal of the dentition pattern of this species. Its first lateral teeth were misidentified as anterior ones and the latter were probably not recognized as such and attributed to A. tenuis or A. longidens because anterior teeth of these three species are quite similar. This explains the apparent rarity of A. magnus anterior teeth in the collection of the Natural History Museum, London noted by REES & UNDERWOOD (2008). The arrangement of the dentition of A. magnus appears quite different from the one of A. medius and is more reminiscent of what can be observed in the older A. smithwoodwardi.

Figure captions Fig. 1 .Fig. 2 .

 captions12 Figure captions

Fig. 3 .

 3 Fig. 3. Interpretation of the dentition of MPV.2010.3.44. LLSL: Left lower second lateral teeth, LP: Lower posterior teeth, LUFA: Left upper first anterior tooth, LUFL: Left upper first lateral teeth, LUSA: Left upper second anterior teeth, LUSL: Left upper second lateral teeth, RUFL: Right upper first lateral tooth, RUSL: Right upper second lateral teeth, S: Symphyseal tooth.

Table File

 File

		Symphyseal First anterior Second	First lateral Second	Posterior
				Anterior		lateral
	Number 4	9	12	9	11	9
			2 of them	3 of them	4 of them	3 of them
			worn	worn	worn	worn
	Maximal	7,2	9,2	9	8	10,2 (average	8,4
	basal		(average on	(average on		on 3 teeth	(average on 6
	height		8 teeth (8,7	5 teeth (8,8		(9,7 to 11))	teeth (8,2 to
			to 9,7))	to 9,3))			8,7))
	Maximal	5,1	4	3,6	6,7	6	2,7
	coronary		(average on	(average on		(average on 3	(average on 6
	height		8 teeth (3,8	5 teeth (3,5		teeth (5,7 to	teeth (2,5 to
			to 4,2))	to 3,8))		6,,2))	3,1))
	Maximal	15,4	15,3	18,6	23,5	45,4 (average	12,2 (average
	mesio-		(average on	(average on	(average on	on 4 teeth	on 3 teeth
	distal		8 teeth (14	7 teeth (17	5 teeth	(44,8 to	(10,8 to
	length		to 16))	to 19,6))	(22,8 to	45,9))	14,9))
					24))
	Maximal	9,1	8 (average	8 (average	14,5	LT: 18,9	Base: 13,5
	labio-		on 8 teeth	on 7 teeth	(average on	(average on 3	(average on 6
	lingual		(7,8 to 8,1))	(7,5 to 8,5))	5 teeth	teeth (18,7 to	teeth (12,1 to
	length				(14,2 to	19,2))	16))
					14,8))	UT: 20,8	Crown: 9,2
						(average on 4	(average on 6
						teeth (20,4 to	teeth (8 to
						21,3))	11,5))

Table 1 .

 1 Numbers and measurements of the different kind of teeth identified in MPV.2010.3.44. When possible, measurements were averaged and the smallest and the largest ones are given between brackets. LT: lower teeth, UT: upper teeth. All measurements are in mm.

Acknowledgements

We are thankful to Laurent Picot (Paléospace, Villers-sur-Mer) who put MPV.2010.3.44 at our disposal for study and to Charlie Underwood for useful discussion. Sarah Rigal provided us with hard-to-find references. We also thank to our two reviewers, Friedrich Pfeil and Jürgen Kriwet, who help us to improve the manuscript. This work was supported by UMR 5276.