A Lebée

K Sab

A Bending-Gradient theory for thick laminated plates homogenization

This work presents a new plate theory for out-of-plane loaded thick plates where the static unknowns are those of the Love-Kirchhoff theory, to which six components are added representing the gradient of the bending moment. The Bending-Gradient theory is an extension to arbitrary multilayered plates of the Reissner-Mindlin theory which appears as a special case when the plate is homogeneous. The new theory is applied to multilayered plates and its predictions are compared to full 3D Pagano's exact solutions and other approaches. It gives good predictions of both deflection and shear stress distributions in any material configuration. Moreover, under some symmetry conditions, the Bending-Gradient model coincides with the second-order approximation of the exact solution as the slenderness ratio L/h goes to infinity.

Introduction

Laminated plates are widely used in engineering applications. For instance angleply carbon fiber reinforced laminates are commonly used in aeronautics. However, these materials are strongly anisotropic and the plate overall behavior is difficult to capture. The most common plate theory is the Love-Kirchhoff plate model. However, it is well-known that, when the plate slenderness ratio L/h is not large enough, transverse shear stresses which are not taken into account in the Love-Kirchhoff theory have an increasing influence on the plate deflection.

In recent decades many suggestions have been made to improve the estimation of transverse shear stresses. Two main approaches can be found: asymptotic approaches and axiomatic approaches. The first one is mainly based on asymptotic expansions in the small parameter h/L [START_REF] Caillerie | Thin elastic and periodic plates[END_REF][START_REF] Lewinski | Effective Models Of Composite Periodic Plates .1. Asymptotic Solution[END_REF]. However, no distinction between relevant fields and unknowns was made. The second main approach is based on assuming ad hoc displacement or stress 3D fields. These models can be "Equivalent Single Layer" or "Layerwise". Equivalent single layer models treat the whole laminate as an equivalent homogeneous plate. However, when dealing with laminated plates, these models lead systematically to discontinuous transverse shear stress distributions through the thickness as indicated by Reddy [START_REF] Reddy | On Refined Computational Models Of Composite Laminates[END_REF]. In Layerwise models, all plate degrees of freedom are introduced in each layer of the laminate and continuity conditions are enforced between layers. The reader can refer to Reddy [START_REF] Reddy | On Refined Computational Models Of Composite Laminates[END_REF] and Carrera [START_REF] Carrera | Theories and finite elements for multilayered, anisotropic, composite plates and shells[END_REF] for detailed reviews of kinematic approaches and to [START_REF] Diaz | Model for laminates[END_REF][START_REF] Caron | Multi-particle modelling for the prediction of delamination in multi-layered materials[END_REF][START_REF] Diaz | Analytical determination of the modes I, II and III energy release rates in a delaminated laminate and validation of a delamination criterion[END_REF] for static approaches. Layerwise models lead to correct estimates of local 3D fields. However, their main drawback is that they involve a number of degrees of freedom proportional to the number of layers. The limitation is immediately pointed out with functionally graded materials, where the plate constituents properties vary continuously through the thickness [START_REF] Nguyen | First-order shear deformation plate models for functionally graded materials[END_REF][START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF].

Based on Reissner [START_REF] Reissner | The effect of transverse shear deformation on the bending of elastic plates[END_REF] paper, we suggest an Equivalent Single Layer higherorder plate theory which gives an accurate enough estimate of transverse shear stresses in the linear elasticity framework. For this, we are motivated by two observations. The first one is that Love-Kirchhoff strain fields have clearly been identified as good first-order approximation for slender plates thanks to asymptotic expansion approaches. The second one is that the 3D equilibrium plays a critical role in the estimation of transverse shear stress in all the existing approaches. We show in this work that revisiting the use of 3D equilibrium in order to derive transverse shear stress as Reissner [START_REF] Reissner | The effect of transverse shear deformation on the bending of elastic plates[END_REF] did for homogeneous plates leads to a full bending gradient plate theory. The Reissner-Mindlin theory is as a special case of the new Bending-Gradient theory when the plate is homogeneous.

In Section 2 notations are introduced. In Section 3, we resume Reissner's procedure for deriving transverse shear stress extended to laminated plates. This lead to the Bending-Gradient plate theory detailed in Section 4. Finally, in Section 5 the Bending-Gradient plate theory is applied to fibrous laminates under cylindrical bending and compared to the exact solution and other Single Equivalent Layer approaches.

Notations

Vectors and higher-order tensors are boldfaced and different typefaces are used for each order: vectors are slanted: T T T , u u u. Second order tensors are sans serif: M M M, e e e. Third order tensors are in typewriter style: Φ Φ Φ, Γ Γ Γ. Fourth order tensors are in calligraphic style D D D, c c c. Sixth order tensors are double stroked , Ï Ï Ï.

When dealing with plates, both 2-dimensional (2D) and 3D tensors are used. Thus, T T T denotes a 3D vector and T T T denotes a 2D vector or the in-plane part of T T T . The same notation is used for higher-order tensors: σ σ σ is the 3D second-order stress tensor while σ σ σ is its in-plane part. When dealing with tensor components, the indexes specify the dimension: a i j denotes the 3D tensor a a a with Latin index i, j, k.. = 1, 2, 3 and a αβ denotes the 2D a a a tensor with Greek indexes α, β , γ..

(δ αβ = 1 if α = β , δ αβ = 0 otherwise).
The transpose operation t • is applied to any order tensors as follows:

( t A) αβ ...ψω = A ωψ...β α .
Three contraction products are defined, the usual dot product ( a a a • b b b = a i b i ), the double contraction product ( a a a : b b b = a i j b ji ) and a triple contraction product (A A A ∴ B B B = A αβ γ B γβ α ). In these definitions Einstein's notation on repeated indexes is used. It should be noticed that closest indexes are summed together in contraction products. Thus,

Φ Φ Φ • n n n = Φ αβ γ n γ is different from n n n • Φ Φ Φ = n α Φ αβ γ . The derivation operator ∇ ∇ ∇
is also formally represented as a vector: a a a • ∇ ∇ ∇ = a i j ∇ j = a i j, j is the divergence and a a a ⊗ ∇ ∇ ∇ = a i j ∇ k = a i j,k is the gradient. Here ⊗ is the dyadic product. Finally, the integration through the thickness is noted • :

h 2 -h 2 f (x 3 )dx 3 = f .

Revisiting the Reissner-Mindlin plate theory

The 3D model We consider a linear elastic plate of thickness h occupying the 3D domain Ω = ω×]h/2, h/2[, where ω ⊂ R 2 is the mid-plane of the plate (Figure 1). Cartesian coordinates (x 1 , x 2 , x 3 ) in the reference frame ( e e e 1 , e e e 2 , e e e 3 ) are used. The constitutive material is assumed to be invariant with respect to translations in the (x 1 , x 2 ) plane. Hence, the stiffness tensor C C C is a function of x 3 only. The plate is loaded on its upper and lower faces ω ± = ω × {±h/2} with the distributed force T T T ± . There are no body forces and the plate is clamped on its lateral edge, ∂ ω×]h/2, h/2[ where ∂ ω is the edge of ω. The 3D problem P 3D is summarized as follows:

P 3D                    σ σ σ • ∇ ∇ ∇ = 0 on Ω . (1a) σ σ σ = C C C (x 3 ) : ε ε ε on Ω . (1b) σ σ σ • e e e 3 = T T T ± on ω ± . (1c) ε ε ε = 1 2 ∇ ∇ ∇ ⊗ u u u + u u u ⊗ ∇ ∇ ∇ on Ω . ( 1d
) u u u = 0 on ∂ ω×] -h/2, h/2[. ( 1e
)
where u u u is the 3D displacement vector field, ε ε ε is the strain tensor field and σ σ σ is the stress tensor field.

Fig. 1 The Plate Configuration

Reissner-Mindlin statically compatible fields

We recall here briefly the procedure for the derivation of Reissner-Mindlin equilibrium equations [START_REF] Reissner | The effect of transverse shear deformation on the bending of elastic plates[END_REF][START_REF] Caron | Un nouveau modèle de plaque multicouche épaisse. A new model for thick laminates[END_REF][START_REF] Nguyen | A model for thick laminates and sandwich plates[END_REF]. The generalized Reissner-Mindlin stresses associated to the 3D stress field σ σ σ are:

N αβ (x 1 , x 2 ) = σ αβ (2)
M αβ (x 1 , x 2 ) = x 3 σ αβ (3) Q α (x 1 , x 2 ) = σ α3 (4)
where N N N is the membrane stress, M M M the bending moment, and Q Q Q the shear forces. Reissner-Mindlin equilibrium equations are obtained by integrating equations (1a) and x 3 ×(1a) with respect to x 3 . Taking also into account boundary conditions (1c) yields:

     N αβ ,β + q α = 0 (5a) Q α,α + q 3 = 0 (5b) M αβ ,β -Q α + µ α = 0 (5c)
where q i = T + i + T - i are symmetric loadings per unit surface and µ i = h 2 (T + i -T - i ) are skew-symmetric loadings per unit surface. More precisely, q q q = (q α ) are membrane loadings per unit surface, q 3 is the out-of-plane loading per unit surface, µ µ µ = (µ α ) are couples per unit surface and µ 3 is the transverse bulk loading.

Since in-plane loadings (q q q, µ µ µ) and out-of-plane loadings (q 3 , µ 3 ) are not of the same order in the asymptotic analysis of the plate as h/L goes to 0 (see [START_REF] Lewinski | Effective Models Of Composite Periodic Plates .1. Asymptotic Solution[END_REF]), and for the sake of simplicity, we focus only on the out-of-plane loading q 3 (q α = µ i = 0).

Localization

The second step of Reissner's approach consists in deriving the stress energy per unit surface w * RM (N N N,M M M,Q Q Q) from the 3D model. As in many homogenization procedures, the derivation of w * RM is based on an approximation scheme for the real 3D stress fields in terms of the generalized plate stress fields:

σ σ σ RM (x 1 , x 2 , x 3 ) = σ σ σ (N) (x 1 , x 2 , x 3 ) + σ σ σ (M) (x 1 , x 2 , x 3 ) + σ σ σ (Q) (x 1 , x 2 , x 3 )
where σ σ σ (N) , σ σ σ (M) , and σ σ σ (Q) are 3D stress fields generated by N N N, M M M and Q Q Q as follows:

         σ (N) i j = s (N) i jαβ (x 3 )N αβ (x 1 , x 2 ) (6a) σ (M) i j = s (M) i jαβ (x 3 )M αβ (x 1 , x 2 ) (6b) σ (Q) i j = s (Q) i jα (x 3 )Q α (x 1 , x 2 ) (6c)
where s

(N) i jαβ (x 3 ), s (M) i jαβ (x 3 ) and s (Q) i jα (x 3
) are localization tensors depending only on the x 3 coordinate. This can be rewritten using contraction products as:

σ σ σ RM = s s s (N) : N N N + s s s (M) : M M M + s s s (Q) •Q Q Q
Once this approximation of stress fields is set, the stress potential energy density w * RM (N N N,M M M,Q Q Q) is defined simply as the quadratic form:

w * RM (N N N(x x x),M M M(x x x),Q Q Q(x x x)) = 1 2 σ σ σ RM ( x x x) : S S S (x 3 ) : σ σ σ RM ( x x x) (7)
Hence, a consistent choice for s s s (N) , s s s (M) and s s s (Q) is critical.

Love-Kirchhoff fields

The derivation of s s s (N) and s s s (M) is based on the Love-Kirchhoff plate theory. According to this theory, plane-stress is assumed and the in-plane part of the strain is linear in x 3 :

ε ε ε LK = e e e + x 3 χ χ χ (8)
where e e e is the membrane strain and χ χ χ the curvature. We draw the reader's attention to the fact that strain components ε i3 are not null in the general case. Membrane stress N N N and bending moments M M M are linearly dependent on e e e and χ χ χ:

N N N = A A A : e e e + B B B : χ χ χ (9a) M M M = t B B B : e e e + D D D : χ χ χ (9b)
with:

(A A A,B B B,D D D) = 1, x 3 , x 2 3 c c c(x 3 ) (10)
Using 3D constitutive equation under plane-stress assumption, Love-Kirchhoff constitutive equation ( 9) and in-plane strains definition [START_REF] Diaz | Analytical determination of the modes I, II and III energy release rates in a delaminated laminate and validation of a delamination criterion[END_REF], it is possible to express Love-Kirchhoff stress fields as functions of N N N and M M M:

σ σ σ (N) (x 1 , x 2 , x 3 ) = c c c(x 3 ) : a a a + x 3 t b b b : N N N(x 1 , x 2 ) and σ (N) i3 = 0 (11a) σ σ σ (M) (x 1 , x 2 , x 3 ) = c c c(x 3 ) : (b b b + x 3 d d d ) : M M M(x 1 , x 2 ) and σ (M) i3 = 0 (11b)

Stress field generated by a linear variation of the bending moment

The main idea of Reissner's method is to recall that the shear forces are related to the bending moment through the plate equilibrium [START_REF] Carrera | Theories and finite elements for multilayered, anisotropic, composite plates and shells[END_REF]. With a homogeneous plate, combining both 3D equilibrium and plate equilibrium enables the derivation of a stress field directly depending on shear forces. However, with laminated plates it is not possible to bring out shear force with this procedure. Here, we suggest considering a more general shear variable for laminates, the full bending gradient: R R R = M M M ⊗∇ ∇ ∇. In the following, we resume the procedure from Reissner [START_REF] Reissner | The effect of transverse shear deformation on the bending of elastic plates[END_REF] for deriving shear fields in the case of laminated plates.

We have σ σ σ

(M) • ∇ ∇ ∇ = 0 if M M M is (x 1 , x 2 )-invariant. When M M M is function of x 1 and x 2 ,
we have:

σ σ σ (M) • ∇ ∇ ∇ = s (M) i jβ α (x 3 )M αβ (x 1 , x 2 )∇ j = s (M) i jβ α M αβ ,γ δ jγ = s (M) iγβ α R αβ γ f (R) i = s (M)
iγβ α R αβ γ is the force per unit volume generated by first order variations of the bending moment R R R. R αβ γ is a third-order tensor which respects M αβ symmetries (R αβ γ = R β αγ ). Using σ σ σ (M) definition (Equation 11b) and assuming that each layer follows monoclinic symmetry we identify the force per unit volume as:

f f f (R) = c c c(x 3 ) : (b b b + x 3 d d d ) ∴ R R R and f (R) 3 = 0 (12)
Then, we define σ σ σ (R) the 3D stress generated by a (x 1 , x 2 )-invariant bending gradient R R R associated to the localization tensor s (R)

i jαβ γ such as σ σ σ (R) = s s s (R) ∴ R R R.
This stress field is derived through the auxiliary problem:

σ σ σ (R) • ∇ ∇ ∇ + f f f (R) = 0 0 0 (13a) σ σ σ (R)
• e e e 3 = 0 0 0 for

x 3 = ±h/2 (13b)
The (x 1 , x 2 )-invariant solution of this problem is easily found, leading to the explicit determination of s s s (R) :

s (R) α3ηζ ε (x 3 ) = - x 3 -h 2 c αηγδ (z) b δ γεζ + z d δ γεζ dz , s (R)
αβ ηζ ε = 0 and s (R)

33ηζ ε = 0 (14)
We have derived a localization tensor s s s (R) which depends on all bending gradient components: R αβ γ = M αβ ,γ . Accordingly we define a new approximation of stress fields involving all bending gradient components:

σ σ σ BG = σ σ σ (N) + σ σ σ (M) + σ σ σ (R)
and a new stress energy density identical to Definition 7:

w * BG (N N N,M M M,R R R)
Actually σ σ σ BG approximation for 3D stress fields is a higher-order gradient theory, as described in [START_REF] Boutin | Microstructural effects in elastic composites[END_REF] for 3D continuum and [START_REF] Buannic | Higher-order effective modeling of periodic heterogeneous beams. I. Asymptotic expansion method[END_REF] for periodic beams. However, to be consistent with higher-order theories, we should have taken into account the gradient of other static unknowns such as the membrane stress gradient for instance. It is the choice of the authors to limit the number of static variables only to those which have a contribution to the macroscopic equilibrium of the plate. Thus the number of unknowns remains limited and adapted to engineering applications, contrary to asymptotic expansions and other rigorous approaches in which no distinction is made between significant static unknowns.

The mechanical meaning of the gradient of the bending moment

The full bending gradient R R R has six components (R 111 , R 221 , R 121 , R 112 , R 222 , R 122 )
whereas Q Q Q has two components. Thus, using the full bending gradient as static unknown introduces four static unknowns which a priori are not related to plate equilibrium (5c). Only (N N N,M M M,Q Q Q) appeared while integrating 3D equilibrium equation (1a) through the thickness in Section 3.

Let us derive generalized stresses associated to σ σ σ (R) . Using Equation 14 and integrating by parts when necessary leads to:

σ (R) αβ = 0, x 3 σ (R) αβ = 0, s (R) α3β γδ = i αβ γδ (15)
and we have:

σ (R) α3 = i i i ∴ R R R = Q Q Q. Only Q Q Q remains after integrating σ σ σ (R) through
the thickness and the four other static unknowns are self-equilibrated stress. More precisely we have:

R 111 R 221 R 121 R 112 R 222 R 122 σ 13 s (R) 13111 = 1 s (R) 13122 = 0 s (R) 13121 = 0 s (R) 13211 = 0 s (R) 13222 = 0 s (R) 13221 = 1/2 σ 23 s (R) 23111 = 0 s (R) 23122 = 0 s (R) 23121 = 1/2 s (R) 23211 = 0 s (R) 23222 = 1 s (R) 23221 = 0

The Bending-Gradient plate model

Once stress energy density w * BG (N N N,M M M,R R R) and plate equilibrium equations 5 to which is added R R R = M M M ⊗ ∇ ∇ ∇ are given it is possible to build a complete plate theory using conventional variational tools. The reader is refered to [START_REF] Lebée | A bending gradient model for thick plates, Part I: Theory[END_REF] for details.

Summary of the plate model

Equilibrium equations and boundary conditions involving stress fields are gathered in the set of statically compatible fields:

                   N N N •∇ ∇ ∇ = 0 0 0 on ω (16a) M M M ⊗ ∇ ∇ ∇ -R R R = 0 on ω (16b) (i i i ∴ R R R) •∇ ∇ ∇ = -q 3 on ω (16c) N N N •n n n = V V V d on ∂ ω s (16d) M M M = M M M d on ∂ ω s (16e) (i i i ∴ R R R) •n n n = V d 3 on ∂ ω s (16f)
where ∂ ω s is the portion of edge on which static boundary conditions apply:

V V V d
is the force per unit length and M M M d the full bending moment enforced on the edge. This set of equations is almost identical to Reissner-Mindlin equations where shear forces have been replaced by the bending gradient R R R.

Generalized stresses N N N, M M M, and R R R work respectively with the associated strain variables: e e e, the conventional membrane strain, χ χ χ the conventional curvature and Γ Γ Γ the generalized shear strain. These strain fields must comply with the following compatibility conditions and boundary conditions:

              
e e e = i i i :

(∇ ∇ ∇ ⊗U U U) on ω (17a) χ χ χ = Φ Φ Φ •∇ ∇ ∇ on ω (17b) Γ Γ Γ = Φ Φ Φ + i i i •∇ ∇ ∇U 3 on ω (17c) Φ Φ Φ •n n n = H H H d on ∂ ω k (17d) U U U = U U U d on ∂ ω k (17e)
where U U U is the 3D displacement of the mid-plane of the plate and Φ Φ Φ is the generalized rotation. Γ Γ Γ and Φ Φ Φ are 2D third-order tensors with the following symmetry: Φ αβ γ = Φ β αγ . Moreover, ∂ ω k is the portion of edge on which kinematic boundary conditions apply: U U U d is a given displacement and H H H d is a symmetric second-order tensor related to a forced rotation on the edge. These fields are almost identical to Reissner-Mindlin kinematically compatible fields where the rotation pseudo-vector is replaced by the generalized rotation Φ Φ Φ.

Finally, for constitutive material following local monoclinic symmetry with respect to (x 1 ,

Γ Γ Γ = ∴ R R R, where (Á Á Á -∴ ) ∴ Γ Γ Γ = 0 (18c)
where conventional Love-Kirchhoff stiffnesses are defined as:

(A A A,B B B,D D D) = 1, x 3 , x 2 3 c c c(x 3 )
and is the generalized shear compliance tensor1 defined as:

= h 2 -h 2 x 3 -h 2 t b b b + zd d d : c c c(z)dz •S S S(x 3 ) • x 3 -h 2 c c c(z) : (b b b + zd d d ) dz dx 3 ( 19
)
where S S S = S αβ = 4S α3β 3 is the out-of-plane shear compliance tensor. Since is not always invertible, we introduced Moore-Penrose pseudo inverse for the shear stiffness tensor :

= lim κ→0 ( ∴ + κÁ Á Á) -1 ∴
where Á Á Á is the identity for 2D sixth-order tensors following the generalized shear compliance minor and major symmetries (Á αβ γδ εζ = i αβ εζ δ γδ ). The solution of the plate model must comply with the three sets of equations [START_REF] Lebée | A bending gradient model for thick plates, Part I: Theory[END_REF][START_REF] Lebée | A Cosserat multiparticle model for periodically layered materials[END_REF][START_REF] Lebée | A bending gradient model for thick plates, Part II: Closedform solutions for cylindrical bending[END_REF]. The compliance is positive. However when is not definite, there is a set of solutions, up to a self-stress field.

Projection of the Bending-Gradient plate model

In some cases, the Bending-Gradient is turned into a Reissner-Mindlin plate model. This is the case for homogeneous plates. Thus, we need a means to estimate the difference between both plate models. It is possible to define the exact projection of the Bending-Gradient model on a Reissner-Mindlin model. The Reissner-Mindlin part of is:

RM = 2 3 i i i •i i i ∴ ∴ 2 3 i i i •i i i ( 20
)
RM can be considered as the restriction of when setting warping unknowns to zero. Consequently, we introduce the pure warping part of as the orthogonal complement of RM : W = -RM (21)

Finally we suggest the following relative distance between the Bending-Gradient plate model and the Reissner-Mindlin one:

∆ RM/BG = W , where = αβ γδ εζ ζ εδ γβ α (22)
∆ RM/BG gives an estimate of the pure warping fraction of the shear stress energy and is a criterion for assessing the need of the Bending-Gradient model. When the plate constitutive equation is restricted to a Reissner-Mindlin one, we have exactly ∆ RM/BG = 0.

Application to laminates

Plate configuration

We consider angle-ply laminates. Each ply is made of unidirectional fiber-reinforced material oriented at θ relative to the bending direction x 1 . All plies have the same thickness and are perfectly bounded. A laminate is denoted between brackets by the successive ply-orientations along thickness. For instance [0 • , 90 • ] denotes a 2-ply laminate where the lower ply fibers are oriented in the bending direction. The constitutive behavior of a ply is assumed to be transversely isotropic along the direction of the fibers and engineering constants are chosen similar to those of [START_REF] Pagano | Exact Solutions for Composite Laminates in Cylindrical Bending[END_REF]:

E L = 25 × 10 6 psi, E T = 1 × 10 6 psi, G LT = 0.5 × 10 6 psi, G T T = 0.4 × 10 6 psi, ν LT = ν T T = 0.25
where G T T has been changed to preserve transversely isotropic symmetry. L is the longitudinal direction oriented in the (x 1 , x 2 ) plane at θ with respect to e e e 1 , T is the transverse direction.

Distance between the Reissner-Mindlin and the Bending-Gradient model

In Table 1, are given the values of ∆ RM/BG for the laminates considered in this work. For a single ply, the criterion is zero since the Bending-Gradient model is exactly a Reissner-Mindlin model in this case. However, when there are several plies, the distance is greater than 10%. Thus with these laminates, the shear constitutive equation cannot be reduced to a Reissner-Mindlin behavior.

Stack [0 • ] [0 • , 90 • ] [30 • , -30 • ] s [30 • , -30 • ]
∆ RM/BG 0 16.0% 16.0% 23.9 % Table 1 The criterion ∆ RM/BG for several laminates

Cylindrical bending

Pagano [START_REF] Pagano | Exact Solutions for Composite Laminates in Cylindrical Bending[END_REF] gives an exact solution for cylindrical bending of simply supported composite laminates. We choose the same configuration for the Bending-Gradient model. The plate is invariant and infinite in x 2 direction. It is out-of-plane loaded with q 3 (x 1 ) = -q 0 sin κx 1 where λ = 1/κ is the wavelength of the loading (Figure 2). The plate is simply supported at x 1 = 0 and x 1 = L with traction free edges:

U 3 (0) = 0, U 3 (L) = 0, M M M (0) = 0, M M M (L) = 0, N N N (0)
•e e e 1 = 0 0 0, N N N (L) e e e 1 = 0 0 0 (23)

M 22 (0) = M 22 (L) = 0 is the additional boundary condition compared to the Reissner-Mindlin plate model. This additional boundary condition takes into account free edge effects similar to those described in [START_REF] Lebée | A Cosserat multiparticle model for periodically layered materials[END_REF] for periodically layered laminate. The resolution is provided in details in [START_REF] Lebée | A bending gradient model for thick plates, Part II: Closedform solutions for cylindrical bending[END_REF]. Closed-form solutions using the Reissner-Mindlin model were also derived in order to compare them with the Bending-Gradient. The work of Whitney [START_REF] Whitney | Stress Analysis of Thick Laminated Composite and Sandwich Plates[END_REF] was used for deriving transverse shear stress distributions and shear correction factors were taken into account into the shear constitutive equation of the Reissner-Mindlin plate model.

A comparison with a finite elements solution was also performed on ABAQUS [START_REF] Abaqus | ABAQUS/Standard user's manual[END_REF]. Since the Bending-Gradient is an Equivalent Single Layer theory, conventional shell elements were chosen (3 displacements and 3 rotations). Transverse shear fields with shell elements in ABAQUS are derived using an approach very similar Fig. 2 Pagano's cylindrical bending configuration to [START_REF] Whitney | Stress Analysis of Thick Laminated Composite and Sandwich Plates[END_REF] where it is furthermore assumed that the plate overall constitutive equation is orthotropic with respect to the main bending direction. S4, linear quadrangle with full integration elements, were used. A convergence test was performed. This study enforced the typical size of an element l char = h/5 where h is the plate thickness. For instance when the slenderness is h/L = 1/4 there are 20 elements. Figure 3 shows a typical deformation of this mesh. Periodicity was enforced on lateral edges of the strip in Figure 3 by equating corresponding rotations and displacements. Finally, section integration is performed during the analysis. Fig. 3 Finite Element undeformed and deformed mesh for an anisotropic laminate Two error estimates are introduced: the first one for the transverse shear part of the stresses for which we introduce the following seminorm:

σ σ σ 2 = L 2 -L 2 h 2 -h 2 σ α3 S α3β 3 σ β 3 dx 3 dx 1
and we define the relative error as: ∆ (σ σ σ ) = σ σ σ Ex -σ σ σ σ σ σ Ex , where σ σ σ Ex is the exact shear stress distribution from Pagano [START_REF] Pagano | Exact Solutions for Composite Laminates in Cylindrical Bending[END_REF][START_REF] Pagano | Exact Solutions for Rectangular Bidirectional Composites and Sandwich Plates[END_REF][START_REF] Pagano | Influence of Shear Coupling in Cylindrical. Bending of Anisotropic Laminates[END_REF]. The second one is the mid-span deflection relative error:

∆ (U 3 ) = U Ex 3 (L/2)-U 3 (L/2) U Ex 3 (L/2)
, where U Ex 3 (x 1 ) is the plate deflection taken for the exact solution.

Results

First, we consider a skew-symmetric cross ply [0 • , 90 • ] laminate. In this case, the plate configuration fulfills the assumptions made for the finite elements approximation (orthotropic laminate). In Figure 4, shear stress distribution in both directions are plotted for the exact solution from Pagano [START_REF] Pagano | Exact Solutions for Composite Laminates in Cylindrical Bending[END_REF] σ σ σ Ex , the Bending-Gradient so- lution σ σ σ (R) , Whitney's shear distribution σ σ σ (Q),W and the finite elements solution σ σ σ (Q),FE . The slenderness ratio is set to L/h = 4 as conventionally done when benchmarking plate models. The three approximate solutions yield the same distribution. The discrepancy with the exact solution is well-known and associated to edge effects. In Figure 5 the transverse shear stress distribution error ∆ (σ σ σ ) versus the slen- In this case, Whitney's solution converges with L/h whereas finite elements and Bending-Gradient approximations do not converge and lead to rather small errors (≃ 10 -3 ). In Figure 6 the mid-span deflection error is also plotted versus the slenderness ratio. The three approximate solutions yields almost the same error.

0 0.5 1 1.5 2 2.5 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 hσ13/(λq0) x3/h (a) σ Ex 13 σ (R) 13 σ (Q),FE 13 σ (Q),W 13
We consider now a symmetric and non-orthotropic [30 • , -30 • ] s laminate. This configuration does not comply with the assumptions made for the finite elements approach. In Figure 7 shear distributions are compared to the exact solution. The Bending-Gradient solution remains close to the exact solution. However finite elements and Whitney's solution yield different distributions which are not as accurate as the Bending-Gradient. More precisely, in Direction 2, the FE solution does not capture the change of slope associated to the change of ply orientation. In Direction 1 the macroscopic equilibrium is respected for all approximated solutions ( σ 13 = Q 1 ). However in Direction 2 we can see that σ FE 23 = Q 2 for both finite elements and Whitney's solution. In Figure 8 error versus the slenderness ratio is plotted. Contrary to the finite elements solution and Whitney's solution, the Bending-Gradient solution converges when the plate is slender. More precisely we have: ∆ (σ σ σ BG ) ∝ h L 2 in this case. In Figure 9 the mid-span deflection error is also plotted versus the slenderness ratio. The Bending-Gradient solution is the most accurate one for conventional slenderness.

Finally, in Figure 10 the comparison is made for a non-symmetric and nonorthotropic ply [30 • , -30 • ]. Again, this configuration does not comply with the assumptions made for the finite elements approach. The Bending-Gradient solution remains close to the exact solution and Whitney's solution yields acceptable results (except a mismatch for σ W 23 ). However in this case, finite elements yields inappropriate results: in Direction 1 the stress distribution does not respect macroscopic equilibrium σ FE 13 = Q 1 . We checked nevertheless that FE nodal forces ful- the transverse shear stress distribution error versus the slenderness ratio whereas the Bending-Gradient converges as ∆ (σ σ σ BG ) ∝ h L 2 and both the Whitney and finite elements solutions lead to non-negligible errors. Again, in Figure 12, the deflection error indicates that FE are too compliant and that the Bending-Gradient is more accurate than the Reissner-Mindlin solution.

(Q Q Q),

Discussion

We have compared three approaches for deriving an approximation of the exact solution for cylindrical bending suggested by Pagano [START_REF] Pagano | Exact Solutions for Composite Laminates in Cylindrical Bending[END_REF][START_REF] Pagano | Exact Solutions for Rectangular Bidirectional Composites and Sandwich Plates[END_REF][START_REF] Pagano | Influence of Shear Coupling in Cylindrical. Bending of Anisotropic Laminates[END_REF].

First we derived closed-form solutions for the Reissner-Mindlin model using shear correction factors and shear distributions from [START_REF] Whitney | Stress Analysis of Thick Laminated Composite and Sandwich Plates[END_REF]. This approach yields a fair estimation of the deflection and shear distributions in cylindrical bending but it is not as accurate as the Bending-Gradient approximation in most cases. The main limitation of this approach is the cylindrical bending assumption. It is not sure that shear correction factors and shear distributions will remain valid with more general plate boundary conditions, especially involving torsion, whereas the Bending-Gradient theory is not limited to cylindrical bending.

Second, we implemented a finite elements approximation, using conventional shell elements. This approach assumes both cylindrical bending and orthotropy in the same direction. When these assumptions are not valid, the results might be really affected both for deflection and stress distribution as demonstrated for the [30 • , -30 • ] laminate.

Finally, the Bending-Gradient solution was presented. This approach enables the derivation of stress distributions and gives good enough deflection and stress distribution estimates whatever the plate configuration and the bending direction are. Moreover, it was numerically demonstrated that in some configurations the Bending-Gradient solution converges with the slenderness ratio.

Let us state this convergence condition precisely. We chose to neglect the gradient of membrane stress N N N ⊗ ∇ ∇ ∇ since it is not related to macroscopic stress. In the cylindrical bending configuration, the membrane stress is reduced to N 22 . When N 22 = 0, the membrane stress gradient vanishes. This is the case for [30 • , -30 • ] s and [30 • , -30 • ] since they are balanced laminates (as many θ plies as -θ plies).

In these cases, the Bending-Gradient solution converges because the stress fields related to the membrane stress gradient do not contribute to the final solution. It is possible to generalize this result to any boundary conditions with mirror symmetric laminates (B B B = 0) for which the membrane problem is fully uncoupled from the out-of-plane problem. With these laminates, the Bending-Gradient solution is the Saint-Venant solution for an out-of-plane loaded plate.

Conclusion

In the present paper, we derived a new plate theory, the Bending-Gradient theory, which is the extension of Reissner-Mindlin theory to laminated plates. Comparison with the exact solution for cylindrical bending of cross ply laminates demonstrates that the Bending-Gradient gives good predictions of both deflection and shear stress distributions in any material configuration. It is also the Saint-Venant solution when membrane stresses are fully uncoupled from bending moments and generalized shear stresses. Finally, with usual laminated plates, we pointed out that the Bending-Gradient cannot be reduced to a Reissner-Mindlin plate model.

Several outlooks are under consideration. First, this plate theory can be extended to periodic plates such as sandwich panels [START_REF] Lebée | Reissner-Mindlin Shear Moduli of a Sandwich Panel with Periodic Core Material[END_REF][START_REF] Lebée | Transverse shear stiffness of a chevron folded core used in sandwich construction[END_REF]. Second, the estimation of the influence of the membrane stress gradient on the quality of the shear stress estimation should be studied in detail.

  = 1, 2. C C C = C i jkl is the fourth-order 3D elasticity stiffness tensor. S S S = S i jkl = C C C -1 is the fourth-order 3D elasticity compliance tensor while c c c = c αβ γδ denotes the plane-stress elasticity tensor. c c c is not the in-plane part of C C C but it is the inverse of the in-plane part of S S S : c c c = S S S -1 . The identity for in-plane elasticity is i αβ γδ = 1 2 δ αγ δ β δ + δ αδ δ β γ , where δ αβ is Kronecker symbol

  where a a a, b b b and d d d are the reciprocal compliance tensors of the constitutive equation (9).

R

  111 and R 222 are the cylindrical bending part of shear forces Q 1 and Q 2 , R 121 and R 122 are the torsion part of shear forces and R 112 and R 221 are linked to strictly self-equilibrated stresses (warping).

x 2 )

 2 plane (uncoupling between R R R and (N N N,M M M)) the Bending-Gradient plate constitutive equations are written as:    N N N = A A A : e e e + B B B : χ χ χ (18a) M M M = t B B B : e e e + D D D : χ χ χ (18b)

Fig. 4 Fig. 5 Fig. 6

 456 Fig. 4 Normalized shear distribution σ 13 at x 1 = 0 for a [0, 90 • ] laminate, L/h = 4

Fig. 7

 7 Fig. 7 Normalized shear distribution in both directions at x 1 = 0 for a [30 • , -30 • ] s laminate, L/h = 4, a) σ 13 b) σ 23 .

Fig. 8 Fig. 9 Fig. 10

 8910 Fig. 8 Shear stress distribution error versus slenderness ratio for a [30 • , -30 • ] s laminate

Fig. 11 Fig. 12

 1112 Fig. 11 Shear stress distribution error versus slenderness ratio for a [30 • , -30 • ] laminate

αβ γδ εζ follows major symmetry: αβ γδ εζ = ζ εδ γβ α and minor symmetry αβ γδ εζ = β αγδ εζ . Thus there are only

independent components


Fanny Buyens
email: fanny.buyens@cea.fr

Michele A Quinto

Dominique Houzet
email: dominique.houzet@inpg.fr

Adaptive Mesh Reconstruction in X-Ray Tomography

Keywords: X-ray tomography, mesh representation, level-set

This paper presents an X-ray tomographic reconstruction method based on an adaptive mesh in order to directly obtain the typical gray level reconstructed image simultaneously with its segmentation. It also leads to reduce the number of unknows throughout the iterations of reconstruction and accelerates the process of algebraic algorithms. The process of reconstruction is no more based on a regular grid of voxels but on a mesh composed of non regular tetraedra that are progressively adapted to the content of the image. Each iteration is composed by two main steps that successively estimate the values of the mesh elements and segment the sample in order to make the grid adapted to the content of the image. The method was applied on numerical and experimental data. The results show that the method provides reliable reconstructions and leads to drastically reduce the memory storage compared to usual reconstructions based on pixel representation.

Introduction

Tomography reconstruction from projections data is an inverse problem widely used in the medical imaging field. With sufficiently large number of projections over the required angle, the FBP (Filtered Back-Projection) algorithms allow fast and accurate reconstructions. However in the cases of limited views (low dose imaging) and/or limited angle (specific constrains of the setup), the data available for inversion are not complete, then the problem becomes more illconditioned and the results show significant artifacts. In these situations, an alternative approach of reconstruction, based on a discrete model of the problem, consists in using an iterative algorithm or a statistical modelisation of the problem to compute an estimate of the unknown object. These methods are classicaly based on a volume discretization into a set of voxels and provide 3D maps of densities. High computation time and memory storage are their main disadvantages. Moreover, whatever the application is, the volumes are segmented for a quantitative analysis. Numerous methods of segmentation with different interpretations of the contours and various minimized energy functions are offered with results that can depend on their use by the application.

During the last decade, several studies focused on new approaches simultaneously performing reconstruction and segmentation. These methods are based on a modelling of the sample and solve an optimization problem for which the unknown variables are the location of the boundaries and the intensity distribution. Some authors proposed the reconstruction of 3D binary images composed of uniform compact objects totally included in uniform backgrounds [START_REF] Soussen | Polygonal and polyhedral contour reconstruction in computed tomography[END_REF]. The scene is reconstructed using parametric surface models (splines). The parameters are estimated from the data but the final results strongly depend on the initialization. Senasli et al. [START_REF] Senasli | 3D reconstruction of vessel lumen from very few angiograms by dynamic contours using a stochastic approach[END_REF] also proposed to use active contours based on a stochastic process to reconstruct binary images in the context of few data. Using these object-based methods for tomographic reconstruction, Charbonnier [START_REF] Charbonnier | Reconstruction d'images: régularisation avec prise en compte des discontinuités[END_REF] proposed to introduce a priori on the shape, to constrain the length of the curve to be appropriated to the estimated object. Yoon et al [START_REF] Yoon | Simultaneous segmentation and reconstruction: A level set method approach fo limited view computed tomography[END_REF] proposed a simultaneous segmentation and reconstruction method in the case of limited angle of view in CT. The segmentation and reconstruction are achieved by alternating a segmentation step using the a two-phase level set function and reconstruction step using a conjugate gradient to estimate the intensity value. An initial image obtained after few iterations of the reconstructed algorithm is used to initialize the level set function. The method is performed using usual pixelised grids.

In the framework of SPECT tomography, Brankov et al. proposed a method for 2D tomographic image reconstruction based on an irregular image representation [START_REF] Brankov | Tomographic image reconstruction based on a content-adaptive mesh model[END_REF]. The mesh is generated from a FDK reconstruction using a feature map extraction performed on a regular pixelized grid. The classical algebraic ML-EM algorithm is then used to proceed the reconstruction on the irregular mesh. Sitek and al. also proposed a method of reconstruction based on a different representation of the volume: the domain is here defined by a cloud of points and represented by non overlapping tetraedra defined by these points [START_REF] Sitek | Tomographic reconstruction using an adaptive tetrahedral mesh defined by a point cloud[END_REF][START_REF] Pereira | An svd based analysis of the noise properties of a point cloud mesh reconstruction method[END_REF]. From a coarse regular grid, nodes are added to the the point cloud at the locations of large intensity variations during the reconstruction. The method is more complex but strives for a more efficient image representation.

In this context, we propose a method of reconstruction in X-ray tomography that uses a classical iterative algorithm modified to run on an adaptive triangular mesh jointly with a level set method to achieve simultaneous reconstruction and segmentation. This approach provides both the density map of the sample and the segmentation of all the materials that compose the object. The refinement of the mesh, relying on the segmentation, significantly reduces the number of unknowns since the homogeneous regions, containing poor information, are represented by a coarse mesh. No additional post processing steps are required leading to a reliable reproducibility of the process. Besides, the representation of the image by a mesh adapted to its content allows to reduce considerably the memory size of the reconstructed image (reduction of the problem of memory storage) and to deal with highly pixelated detectors.

Materials and method

After a step of initialization, the proposed method is composed of three main stages : reconstruction, segmentation and adaptation of the mesh to the object that successively alternate until convergence.

Initialization

This step consists in generating a first mesh according to the Delaunay's criteria i.e. a triangulation whose circumscribing circle of any facet of the triangulation contains no vertex in its interior (Delaunay triangulation). All the elements of this initial mesh are initialized with the same non-zero constant value. The number of elements N t of this first mesh is chosen as 1 10 N dim < N t < N dim where dim = {2, 3} is the dimension of reconstruction and N = max{N x , N y } with N x and N y the number of pixels through the dimensions of the detector.

Reconstruction

The physics of the radiation determines the nature of the information collected by the detectors of the acquisition system. We consider here the case of X-ray transmission tomography for which every crossed material attenuates the flux of photons according to its nature and its thickness. A difference of intensity is measured in the log scale between the incoming and the outcoming photon flux in the materials which is expressed by the Beer-Lambert law. Discretized, the model is written as following:

p = Hf + ϵ ( 1
)
where f is the vector containing the N discretized values of the object, p is the vector containing the M measurements, ϵ is a vector representing modeling and measurements errors, and H is the M ×N projection matrix. M is equal to n proj × n det where n proj is the number of projections and n det the number of detectors. The elements of matrix H represent the contribution of each element of the mesh to the projection rays. As shown in figure 1, the element h i,j corresponds to the length of the projection ray i intercepted by the element of the mesh j.

Here we investigate a Conjuguate Gradient (CG) optimization algorithm to solve the problem as a least-square estimation. We denote the unknown parameter as f and the estimation is done by minimizing the cost function J(f ):

f = arg min f {J(f )} (2) with J(f ) = ∥p -Hf ∥ 2 (3)
For quadratic cost functions, the iterative CG algorihtm with an optimal descent step is given bellow, with n the indice of the iterations, N max the maximum number of iterations, d the CG descent direction and δ the gradient norm tolerance. g n denotes the gradient of the cost function at the n th iteration g n = ∂J(f ) ∂f .

Fig. 1: The contribution of the element j of matrix H to the projection ray i corresponds to the intercepted length l. n = 0, d 0 = -g 0 while n < N max and norm(g n ) < δ do

α n = -g T n d
⟨dn,H T Hdn⟩

f n+1 = f n + α n d n β = g T n+1 g n+1 g T n g n d n+1 = -g n+1 + βd n n + +
In the above algorithm, the descent step α n is calculated with the goal of minimizing the cost function for the next iteration in the direction of the descent. This means :

α n = arg min α {J(f n + αd n )} (4)
The expression given in the CG algorithm above is the exact solution for eq. 4 while the cost function is quadratic. Even with nonquadratic cost functions, the same strategy can be used to estimate the descent step.

The gradient-based iterative methods of reconstruction require two main steps (a projection and a backprojection) over one iteration [START_REF] Kawata | Constrained iterative reconstruction by the conjugate gradient method[END_REF]. These two main operators are very time consuming. To speed up computation, they were thus implemented on Graphic Process Units (GPU) using a ray-tracing method.

Segmentation

The step of segmentation is performed after the step of reconstruction. It consists in delineating the boundaries of the object(s) which are used for the mesh refinement. The segmentation is performed using the multiple phase level set method. The level set function is defined by an implicit function ϕ(x, t) that moves in R n in relation to time t. The curve described by the interface Γ (t) = {x | ϕ(x, t) = 0} between two regions of the image domain, is determined at time t by zero-level curve of the function ϕ(x, t).

The movement of the curve is defined by the Euler-Lagrange equation :

ϕ t + V x • ∇ϕ = 0 ( 5
)
where ϕ t = ∂ϕ ∂t and V x = ∂x ∂t is the velocity field. We investigate here the work of [START_REF] Chan | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF] and [START_REF] Goldstein | Geometric application of the split bregman method: Segmentation and surface reconstruction[END_REF] whose segmentation approach is based on convex energies. This Globally Convex Segmentation or GCS formulation is based on the observation that the steady state solution of the gradient flow coincides with the steady state of the simplified flow:

∂ϕ ∂t = ( ∇. ∇ϕ |∇ϕ| -µ ( (c 1 -f ) 2 -(c 2 -f ) 2 ) ) (6)
that represents the gradient descent for minimizing the energy:

E(ϕ) = |∇ϕ| + µ < ϕ, r > (7) with r = (c 1 -f ) 2 -(c 2 -f ) 2 .
The global minima is defined by constraining the solution in the interval [0,1] resulting in the optimization problem :

min 0≤ϕ≤1 |∇ϕ| + µ < ϕ, r > (8)
After optimization, the segmented region is found by thresholding the level set function to get :

Γ = x : ϕ(x) > α (9)
In the case of an image composed by only one material, 0 ≤ ϕ(x) ≤ 1 and α is set to 0.5. More generally, if 0 ≤ ϕ(x) ≤ c, with c ≤ 1 and c = max( ϕ(x)), α is set to c 2 . In the case of an image composed of several materials, the regions are then segmented by setting α to the values corresponding to the valleys of the histogram of the image.

Refinement

Given the set of points S defining the above segmented boundaries : S = {x i,j , . . . , x Nj ,j ; . . . ; x i,M , . . . , x NM ,M }, x ∈ R 2 [START_REF] Goldstein | Geometric application of the split bregman method: Segmentation and surface reconstruction[END_REF] where N j is the number of points of the contour j, and M the number of segmented contours. Since the lengths of the contours are different, the optimal number of points describing each boundary is given by: [START_REF] Cgal | Computational Geometry Algorithms Library[END_REF] where, l j is the length of contour j, and l min the length of the smallest contour. k is empirically set to 100. The mesh is then generated using the functions of CGAL library [START_REF] Cgal | Computational Geometry Algorithms Library[END_REF] according to the constrained Delaunay's criteria.

K = k • l j l min

Results

The method was implemented and applied in 2D for the reconstruction of a slice of a statistical model of a human knee which is composed of a single material (bone). Secondly, the method was applied on the common Shepp-Logan phantom which is composed of several materials.

2D Reconstruction of the knee

The first studied phantom is a statistical model of the knee only composed by one material (bone). The projections have been computed using the RT-CT module of the CIVA software [START_REF]CIVA[END_REF] with a typical configuration of a common radiographic exam of the knee: the distance knee -detector is 80 mm and the distance knee -X-ray source 1080 mm. The dimensions of the 2D imaging detector are 151.5mm×151.5mm (with 1024×1024 pixels). Figure 3 shows the results of 2D reconstructions of a slice of the 3D model of the knee using 36 projections with the FDK algorithm (Fig. 3a) and with the ATM method (Fig. 3b) using α = max(ϕ(x)) 2 (one material to be segmented). Figure 3c shows the final mesh Figures 3b and3c correspond respectively The final image and its corresponding final mesh. This latter is composed of 1751 triangles which can be compared to the 1024×1024 pixels (=1048576 pixels) of the FDK reconstruction 3a (about ×600 more elements). As Table 1 sums up, if coded on 64 bits (binary), the pixelized image sizes in memory 8MB while the final irregular meshed image sizes 84KB (ASCII) that is more than hundred times less with in addition the information of segmentation (the mesh images are stored in OFF file format). Figure 3d compares the real shape of the knee (black line) and the level of the central slice with those obtained with FDK reconstruction coupled with an additional step of level set segmentation (green line), and directly with the ATM method (red line)

2D Reconstruction of the Shepp-Logan

The second phantom is a 2D slice of the Shepp-Logan phantom which is shown in figure 6c. It is composed of several structures whose values range from 0 to 2.0. 256 projections were computed with an analytical projector in the case of a fan-beam configuration. The detector is composed by 256 pixels. Figure 4 shows the reconstructions obtained with 4b FDK algorithm and 4c-4d with ATM method (respectively the finale image and its corresponding triangular mesh) using a set of 360 projections. The mesh was initialized with 65536 triangles and reached 9522 triangles at the end on the process. In this case of several materials, parameter α was set to 0.05, 0.18 and 0.4 in order to segment all the sturctures. Figure 4e shows Fig. 4c superimposed with the real contours of the phantom (in red) and Figure 4f shows the profiles obtained from the FDK (red) and ATM (blue) reconstructed images compared with the real one (pink). The ATM profile is quite noiser but

Parameters study

Figure 5 shows the results of 2D reconstructions of the Shepp-Logan phantom when varying the value of parameter k. It shows that this parameter has influence on the definition of small regions. Indeed, when zooming on the top zone of the phantom (images on the right column), we can observe the mesh density increasing with the value of k and therefore a better definition of these three regions. However, the other regions of the phantom (bigger) are not impacted by the variation of k. The value of parameter k has thus to be high when the studied object is composed by small structures. Figure 6 shows the influence of the number of iterations of the first step of reconstruction before segmentation. It shows that a good segmentation of the object is obtained rather fast, between five and ten iterations. The drastic simplication of the mesh generated after segmentation makes the computation faster during the next iterations what speed up the global proposed method.

Conclusions and Perspective

The method described in this paper proposes an alternative image representation in the field of image reconstruction in X-ray tomography. The useful pixel image representation (regular grid) can mismatch with the content of the images and is here replaced by an irregular and adaptive mesh. This kind of representation decreases the number of unknowns since the mesh is fine to well describe the contours of the object and its details, and becomes progressively coarse in the homogeneous regions where there is no information.

The method alternatively runs a current algebraic or statistical algorithm adapted to triangular mesh (in particular the main operators of projection and back-projection), a level set method of segmentation and a step of refinement. In this study, the method has been evaluated on 2D numerical data using an algorithm of conjugate gradient for the reconstruction and a level set approach for the segmentation. The results show that the method provides the common gray level images in addition of a mesh adapted to the content of the images.

The mesh is a means to obtain additional information without further image processing tools application. Each segmented region can be easily extracted and studied regardless of the other ones if necessary.

This kind of algorithm could be used for example in the follow up of pathologies such as the osteoarthritis of the knee. The mesh could thus be initialized with a previous exam to drastically speed up the reconstruction, and the results could be stored in the informatic file of the patient generating less problems of memory size. Indeed Table 1 shows how this representation is efficient for large reconstruction : it decreases the space memory by a factor 1000 for the storage of a 1024×1024 pixels image. Since this representation decreases the number of unknowns, the data sets are light and the parallelization of the two main operators does not significantly accelerate the computation time in 2D. The reconstruction of 3D volumes and the use of 2D detector whose number of pixels increases could thus be efficiently treated with this approach.

At present a 3D projector has been implemented on GPU using a 3D rendering method used in computer graphics adapted to our iterative process and a quad tree hierarchical structure on each face of the volume of reconstruction to speed up the stage of initilization of the rays.

  Fig 2 shows a 3D representation of the phantom and a snapshot of the simulation scene in CIVA.In the follow, let's denote ATM the proposed method for Adaptive Triangular Mesh reconstruction.

Fig. 2 :

 2 Fig. 2: Snapshot of the scene in CIVA software for the simulation of the projections of the phantom showing the 3D phantom of the knee, the detector in green and a part of the trajectory of the X-ray source.

Fig. 3 :

 3 Fig. 3: Reconstructions of a 2D slice of the phantom of the knee using 36 projections with (a) the analytical FDK algorithm, (b) the ATM method -(c) the corresponding irregular mesh of the final ATM image showed in (b) and (d) Comparison of the real contour of the knee (black contour) with the contours obtained from FDK reconstruction followed by a level set segmentation (green contour) and directly with the ATM method (red contour).

Fig. 4 :

 4 Fig. 4: (a) 2D Shepp-Logan (256x256 pixels) (b) FDK reconstruction, (c) ATM reconstruction and (d) its corresponding final mesh (9522 triangles). (e) ATM reconstruction superimposed with the real contour of the numerical phantom and (f) the centered horizontal profiles obtained from FDK and ATM compared with the real one

Fig. 5 :

 5 Fig. 5: Influence of parameter k on the mesh definition and the representation of the image. On the left column (a), (d) and (g) represent the gray level images corresponding to the meshed images (b), (e) and (h) respectively obtained with k=50, k=100 and k=200. On the right column, (c), (f) and (i) images show a zoom of the three little regions at the top of the phantom.

Fig. 6 :

 6 Fig. 6: Influence of the number of iterations of reconstruction on segmentation and mesh generation : results (a) after 5 iterations, (b) 10 iterations and (c) 15 iterations.

Table 1 :

 1 Memory storage size of the reconstructed images: pixel representation vs ATM

	representation
	Phantom	Image representation	Data type	File size
	Shepp-Logan	256×256 pixels	binary (64 bits) 524 MB
		9522 triangles	ASCII	472 KB
	Knee	1024×1024 pixels binary (64 bits) 8 MB
		1751 triangles	ASCII	84 KB


Building the Assembly De Bruijn Graph from an Implicit Suffix Tree

Suffix trees belong to the most studied indexing data structures for strings. Generalised suffix trees can index the substrings of a set of input words. In its construction algorithm Ukkonen used implicit suffix trees, which relax the assumption that a suffix is not a prefix of another suffix, and thus represent some suffixes by internal nodes rather than by leaves. In computational biology, for the sake of genome assembly, one wishes to assemble a target sequence from a multitude of input strings, usually called the reads. Numerous algorithms build assembly related graphs that represent the overlaps of the input reads: for instance, the overlap graph, or a special version of the De Bruijn graph, in which only k-mers occurring within the reads are represented by nodes. We address the question of algorithms for building the overlap graph and the De Bruijn graph directly from indexing data structures of the read set, and investigate it when the index in an implicit generalised suffix tree. We will present the algorithms for the De Bruijn graph and for its contracted version, and show they run in a time that is linear in the size of the index. Generally, this work establishes algorithmic paths to convert classical indexing data structures in graph representations useful for assembly.

Acknowledgments: This work is supported by ANR blanc Projet ANR-12-BS02-0008 Colib'read, Defi MASTODONS SePhHaDe from CNRS.


M H Jabbar

A M'zah

O Hammami

D Houzet

3D MPSoC Design Using 2D EDA tools: Analysis of Parameters

Design space exploration of 3D MPSoC architecture is reported in this paper analyzing the impact of 2D EDA tools to the 3D architecture performance. In particular, we study how 3D performance is varied when changing the EDA tools options intending to highlight design issues of 3D design. Results show that 3D timing performance is affected greatly compared with power consumption and total wirelength.

I. Introduction and Related Works

The number of processing elements will continue to increase in the coming years as projected by ITRS and interconnect is one of the important design problems that reduces the benefits of traditional CMOS scaling method at each generation. 3D integration has become the alternative technology to continue producing higher performance and capacity electronic devices through its short and high density of vertical interconnect using through silicon via (TSV). Design space exploration is one of the important things to be concerned helping designers to evaluate different possible design implementations before it is implemented in real hardware and is particularly important for 3D architecture to be able to choose the architectural candidate with the most performance gain.

A limited number of works have been reported with regards to the design space exploration of 3D architecture. System level design space exploration for 3D architecture is proposed by [START_REF] Priyadarshi | Pathfinder 3D: A flow for system-level design space exploration[END_REF] enabling exploration of different stacking and partitioning schemes and their affect on the performance, power and temperature. Another design space exploration for 3D stacked architecture is presented in [START_REF] Milojevic | Automated Pathfinding tool chain for 3D-stacked integrated circuits: Practical case study[END_REF] [START_REF] Milojevic | Pathfinding: A design methodology for fast exploration and optimisation of 3D-stacked integrated circuits[END_REF] focusing different 3D packaging solutions with logic and memory integration. Design space exploration of 3D architecture focusing on microprocessor and memory architecture is presented in [START_REF] Xie | Design Space Exploration for 3D Architectures[END_REF]. Our previous work of 3D design space exploration is limited to only single tier and thus does not sufficient to have a complete analysis [START_REF] Hammami | The Impact of EDA Tools in 3D IC Design Space Exploration: A Case Study[END_REF].

Differs from the previous reported works, this study perform design space exploration of 2D EDA tools parameters on the 3D architecture. We have chosen different placement and routing options in SoC Encounter and evaluate their impact on the timing slack, power consumption and wirelength of the 3D architecture. Due to the unavailability of 3D design tools capable of doing 3D synthesis, 3D placement, 3D CTS and 3D routing, designing using 2D EDA tools is the only solution. The aim of this study is to analyze how 2D EDA tools are affecting the overall 3D architecture performance since it will not be an issue when using a true 3D design tool. We have extended our previous work by integrating a complete 3D design exploration flow to get more accurate results and analysis.

II. Design Space Exploration

Our 3D design is based on Tezzaron 3D technology which is a two-tier face-to-face bonding using microbump for the inter-tier connection (Figure 2) [START_REF] Patti | Three-Dimensional Integrated Circuits and the Future of Systemon-Chip Designs[END_REF]. Figure 1 shows the design flow used in this work to explore placement and routing options in the place and route tool. Synopsys Design Compiler was used for the logic synthesis while Cadence SoC Encounter was used for place and route of both tiers that run in parallel during the exploration. 3D timing analysis and power analysis was done using Synopsys PrimeTime and PrimePower tool. The design space exploration is conducted using a script that automatically modifies the EDA tool options at each iteration. The selected options of EDA tool as shown in Table 1 have been chosen based on its properties that strongly effect the quality of final netlist for 2D design and thus make it more important to evaluate it in the 3D design.

RTL Design and partitioning

Logic synthesis with timing budgeting

Placement

III. Resuts and Discussion

The results of the exploration are shown from Figure 3 to Figure 6 for timing slack, power and wirelength of the 3D architecture. It can be seen from the figures that the timing slack for both NoC and processor is varied substantially between each exploration while for power consumption as well as total wirelength does not varied very much. The reason is because we have chosen timing driven option for both placement and routing and power driven is only explored at placement stage.

IV. Conclusion

We have presented a design space exploration of 3D MPSoC architecture to evaluate the impact of 2D EDA tool on its performance such as timing slack, power and wirelength. The results suggest that 3D timing performance is varied widely than power consumption and total wirelength.

Figure 1 :Figure 2 :

 12 Figure 1: Design flow for 3D design space exploration

Figure 3 :Figure 4 :Figure 5 :Figure 6 :

 3456 Figure 3: Timing slack (WNS) of processor clock

Table 1 :

 1 EDA tool options for the exploration

		Placement options	Routing options
					Route
	Design	Timing	Power	Timing	Timing
	ID	Driven	Driven	Driven	Driven
					Effort
	1	False	False	False	5
	2	False	False	False	10
	3	False	False	True	5
	4	False	False	True	10
	5	False	True	False	5
	6	False	True	False	10
	7	False	True	True	5
	8	False	True	True	10
	9	True	False	False	5
	10	True	False	False	10
	11	True	False	True	5
	12	True	False	True	10
	13	True	True	False	5
	14	True	True	False	10
	15	True	True	True	5
	16	True	True	True	10