J.-C Passieux

J.-N Périé

M Salaün

A dual domain decomposition method for finite element digital image correlation

Keywords: parallel computing, global DIC, high performance computing, substructuring, finite elements, domain decomposition

The computational burden associated to finite element based digital image correlation methods is mostly due to the inversion of finite element systems and to image interpolations. A non-overlapping dual domain decomposition method is here proposed to rationalize the computational cost of high resolution finite element digital image correlation measurements when dealing with large images. It consists in splitting the global mesh into submeshes and the reference and deformed states images into subset images. Classic finite element digital image correlation formulations are first written in each subdomain independently. The displacement continuity at the interfaces is enforced by introducing a set of Lagrange multipliers. The problem is then condensed on the interface and solved by a conjuguate gradient algorithm. Three different preconditionners are proposed to accelerate its convergence. The proposed domain decomposition method is here exemplified with real high resolution images. It is shown to combine the metrological performances of finite element based digital image correlation and the parallelisation ability of subset based methods.

Introduction

The analysis of the mechanical behavior of materials and structures rely more and more on full-field measurements. It is mainly due the large quantity of data they provide, which are particularly interesting for parameter identification purposes for instance. Among them, Digital Image Correlation (DIC) has become one of the most popular because of a favorable ease of use to accuracy ratio [START_REF] Sutton | Determination of displacements using an improved digital correlation method[END_REF][START_REF] Sutton | Image correlation for shape, motion and deformation measurements: Basic Concepts, Theory and Applications[END_REF] and its ability to deal with 3D measurements on the surface (stereo DIC [START_REF] Luo | Accurate measurement of three-dimensional deformations in deformable and rigid bodies using computer vision[END_REF]) and also in the bulk with Digital Volume Correlation (DVC [START_REF] Bay | Digital volume correlation: three-dimensional strain mapping using x-ray tomography[END_REF][START_REF] Rannou | Three dimensional experimental and numerical multiscale analysis of a fatigue crack[END_REF]). In continuum solid mechanics, finite element based DIC (FE-DIC [START_REF] Sun | Finite element formulation for a digital image correlation method[END_REF][START_REF] Besnard | finite-element displacement fields analysis from digital images: Application to Portevin-le Châtelier bands[END_REF][START_REF] Fehrenbach | A fast algorithm for image registration[END_REF][START_REF] Hild | Comparison of local and global approaches to digital image correlation[END_REF]) has proved to be a relevant choice since (a) it allows for interpolation-free communications with finite element simulations, (b) it significantly reduces the measurement uncertainties with respect to classical subset based approaches [START_REF] Hild | Comparison of local and global approaches to digital image correlation[END_REF], since prescribed continuity of the unknown displacement fields acts as a regularisation. The drawback of FE-DIC over subset-DIC [START_REF] Sutton | Image correlation for shape, motion and deformation measurements: Basic Concepts, Theory and Applications[END_REF] is the computational cost when high resolution is sought for [START_REF] Leclerc | Voxel-scale digital volume correlation[END_REF][START_REF] Passieux | High resolution digital image correlation using Proper Generalized Decomposition: PGD-DIC[END_REF]. Indeed, subset based DIC approaches lead to a set of small independent nonlinear systems of equations which are highly parallelisable, whereas FE-DIC method lead to one global non-linear system whose inversion can become prohibitive with a large number of degrees of freedom [START_REF] Passieux | High resolution digital image correlation using Proper Generalized Decomposition: PGD-DIC[END_REF]. In addition, due to the constantly increasing resolution of photographic sensors (standard sensors provide nowadays 29 million pixels, which can be extended up to 260 million pixels using a piezoelectric pixel shift), the manipulation and interpolation of images (large dense matrices) become more and more an issue. This problem is even more acute with tomographic images.

The two main causes responsible for the high computational cost of high resolution FE-DIC are thus system inversions and image interpolations. A method that addresses the first issue consists in reducing the dimensionnality of the problem. A variable separation technique based on the Proper Generalised Decomposition (PGD) has been proposed recently to decrease the numerical complexity of FE-DIC [START_REF] Passieux | High resolution digital image correlation using Proper Generalized Decomposition: PGD-DIC[END_REF][START_REF] Gomes Perini | A multigrid PGD-based algorithm for volumetric displacement fields measurements[END_REF], while preserving the above mentioned advantages. It was shown that it may drastically reduce the computation time associated to the resolution of the FE systems. Another attempt based on GPU implementation proposed recently to adress both sources of computational cost [START_REF] Leclerc | Digital volume correlation: What are the limits to the spatial resolution[END_REF]. However, the limitations remain when the size of the region of interest increases.

In the community of numerical simulation, a family of solvers, referred to as domain decomposition methods (DDM) has been developed for high performance computing on parallel computer architectures. The domain decomposition methods, such as Schwarz methods [START_REF] Glowinski | First International Symposium on Domain Decomposition Methods for Partial Differential Equations[END_REF], were originally based on overlapping partitions of the studied region. Then non-overlapping methods were preferred because of their ease of implementation. The dual domain decomposition (such as Finite Element Tearing and Interconnecting, FETI [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF]) seeks iteratively the displacement continuity at the interface assuming the equilibrium of the subdomains. On the contrary, primal domain decomposition methods (such as Balancing Domain Decomposition, BDD [START_REF] Mandel | Balancing domain decomposition[END_REF]) prescribe displacement continuity and get iteratively equilibrium at the interface. Mixed [START_REF] Ladevèze | Nonlinear computationnal structural mechanics-New approaches and non-incremental methods of calculation[END_REF][START_REF] Roux | Digital image correlation and fracture: an advanced technique for estimating stress intensity factors of 2D and 3D cracks[END_REF][START_REF] Passieux | A scalable time-space multiscale domain decomposition method: adaptive time scale separation[END_REF] and hybrid [START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF] alternatives have also been proposed that mixes primal and dual approaches. Most often, the global problem is condensed onto the interface degrees of freedom and solved by a Krylov iterative algorithm [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Mandel | Balancing domain decomposition[END_REF]. The latter can be computed in parallel since it involves independent problems on each subdomain.

In the same way, a dual non-overlapping domain decomposition algorithm is proposed herein, in the context of FE-DIC, to alleviate both finite element system solving and image manipulation issues at the same time. The method consists in splitting the domain into a set of subdomains (with submeshes and subset images). The algorithm based on the preconditionned conjugate gradient thus involves a set of smaller independent problems (based on small meshes and small subset images) well suited for parallel processing. Three different preconditionners are proposed with different extra-costs and efficiencies. After describing its principle, the method is used to analyse high resolution images of a real experiment.

The proposed domain decomposition approach 2.1 Digital image correlation

Let us consider two grayscale images f (x) and after g(x) corresponding to the reference and deformed states of the specimen respectively. The displacement field between those two states is denoted u(x), where x ∈ Ω is a point in the region of interest Ω of the image. Given f and g, the DIC problem consists in finding the displacement field u ∈ L 2 (Ω) (where L 2 (Ω) defines Hilbert space [L 2 (Ω)] d with d = 2 for DIC and 3 for DVC) that best fullfils the graylevel conservation equation [START_REF] Horn | Determining optical flow[END_REF]:

u = arg min u(x)∈L 2 (Ω) Ω f (x) -g(x + u(x)) 2 dx (1)
Problem ( 1) is nonlinear. Assuming that g (and thus f ) is differentiable, it can be solved by a Newton-like algorithm. At iteration k the unknown displacement correction δu k = u ku k-1 is assumed to be small enough to allow for a first order Taylor expansion of the deformed state image

g(x + u k ) ≈ g(x + u k-1 ) + δu k T ∇g(x + u k-1 ) (2)
In equation [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF], ∇g(x + u k-1 ) denotes the gradient of image g at non integer pixels positions x + u k-1 . Since g(x + u k-1 ) is supposed to converge to f (x), the gradient ∇g(x + u k-1 ) is approximated by ∇f (x), which does not depend on the unknown and which can then be computed once and for all. The stationnarity conditions associated to the minimisation of the linearised version of Problem (1) yields the following variationnal formulation:

a(δu k , v) = l k (v) ∀v ∈ L 2 (Ω) (3)
where

∀(u, v) ∈ L 2 (Ω) × L 2 (Ω), a(u, v) = Ω v T ∇f ∇f T u dx and l k (v) = Ω v T ∇f f (x) -g(x + u k-1 ) dx
As such, this problem is ill-posed in Hadamar's sense, i.e. the displacement cannot be found pixel-wise. This problem needs to be regularised. For that the unknown displacement u(x) is searched for in an approximation subspace V ⊂ L 2 (Ω) whose dimension is much smaller than the number of pixels in Ω. Generally, V is defined by a set of basis functions N i (x). The approximation of the displacement thus reads δu k (x) = i N i (x) q k i where q k i are the corresponding coefficients of the linear combination, referred to as degrees of freedom (dof) of the displacement correction at iteration k.

By applying the Galerkin method, Problem (3) produces, at iteration k, the following linear system of equations:

M q k = b k (4)
where q k is a vector collecting the dof q k i and

M ij = a(N i (x), N j (x)) and b k j = l k (N j (x)) (5)
The correlation operator M is symmetric positive definite [START_REF] Fedele | Global 2D digital image correlation for motion estimation in a finite element framework: a variational formulation and a regularized, pyramidal, multi-grid implementation[END_REF].

The choice of the approximation subspace V or more precisely of its basis functions N i (x) may be varied: piecewise constant [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] or polynomial functions [START_REF] Sutton | Determination of displacements using an improved digital correlation method[END_REF], Fourier series [START_REF] Roux | Correlation image velocimetry : a spectral approach[END_REF][START_REF] Mortazavi | Image-based continuous displacement measurements using an improved spectral approach[END_REF], rigid body translations and rotations [START_REF] Passieux | Multiscale displacement field measurement using digital image correlation: Application to the identification of elastic properties[END_REF][START_REF] Passieux | A digital image correlation method for tracking planar motions of rigid spheres: application to medium velocity impacts[END_REF], B-splines [START_REF] Cheng | Full-field speckle pattern image correlation with b-spline deformation function[END_REF][START_REF] Elguedj | Isogeometric analysis for strain field measurements[END_REF], NURBS [START_REF] Dufour | Integrated digital image correlation for the evaluation and correction of optical distortions[END_REF], finite elements [START_REF] Sun | Finite element formulation for a digital image correlation method[END_REF][START_REF] Besnard | finite-element displacement fields analysis from digital images: Application to Portevin-le Châtelier bands[END_REF][START_REF] Fehrenbach | A fast algorithm for image registration[END_REF], eXtended finite elements [START_REF] Réthoré | From pictures to extended finite elements: extended digital image correlation (X-DIC)[END_REF][START_REF] Rannou | Three dimensional experimental and numerical multiscale analysis of a fatigue crack[END_REF], separation of variable [START_REF] Passieux | High resolution digital image correlation using Proper Generalized Decomposition: PGD-DIC[END_REF][START_REF] Gomes Perini | A multigrid PGD-based algorithm for volumetric displacement fields measurements[END_REF], mechanical [START_REF] Roux | Stress intensity factor measurements from digital image correlation: post-processing and integrated approaches[END_REF][START_REF] Hild | Calibration of constitutive models of steel beams subject to local buckling by using digital image correlation[END_REF][START_REF] Blaysat | Interface debonding characterization by image correlation integrated with double cantilever beam kinematics[END_REF] and optical analytical function [START_REF] Réthoré | 3D displacement measurements using a single camera[END_REF] or numerical precomputed functions [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF].

As regards computational costs, the so-called subset methods, that rely on piecewise polynomial functions, involve a bloc diagonal correlation operator M. Its resolution is therefore highly parallelisable. In opposition, with almost all the other interpolation types, and in particular with finite elements, the operator M is global. Its size, which is equal to the dimension of the approximation subspace (or the number of dofs), can become prohibitive when high resolution is required [START_REF] Leclerc | Voxel-scale digital volume correlation[END_REF][START_REF] Passieux | High resolution digital image correlation using Proper Generalized Decomposition: PGD-DIC[END_REF][START_REF] Gomes Perini | A multigrid PGD-based algorithm for volumetric displacement fields measurements[END_REF].

As stated in the introduction, the second main computational limitations, which is not peculiar to FE-DIC, is the computation of the right hand side b k . The latter needs to evaluate the deformed state image g at the non-integer pixels values x + u k-1 (x). This implies using interpolation algorithms whose computational complexity scales in O(n log(n)) where n is the number of points.

The computational cost of FE-DIC may become significant when high resolution images and meshes are jointly used. In the following, a domain decomposition algorithm is proposed to alleviate both computational issues associated to the resolution of large systems and the manipulation of large images.

Remark. The subset based approaches [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF][START_REF] Sutton | Determination of displacements using an improved digital correlation method[END_REF] could be seen as a kind of overlapping partition of the region of interest since the size of the step between two subsets can be set smaller than the size of the subset itself. However, all the problems on the subsets are solved independently and the displacement continuity across interfaces is not searched for.

Domain decomposition

Let us now consider a non-overlapping partition of the region of interest Ω into a set of subdomains Ω s . In the same way any quantity with subscript s defines the restriction to Ω s of quantity . For the sake of simplicity, the method will be first presented in the case of two subdomains Ω 1 and Ω 2 without loss of generality. In each subdomain, the linearised correlation problem at iteration k reads: find

δu k s ∈ L 2 (Ω s ), such that ∀s, ∀v s ∈ L 2 (Ω s ), a s (δu k s , v s ) = l k s (v s ) (6)
Problem ( 6) must be solved under the constraint of displacement continuity at the interface Γ = Ω 1 ∩ Ω 2 , see Figure 1. For that a Lagrange multiplier λ ∈ L (L being the ad-hoc space) is introduced to ensure this constraint in a weak sense:

u 1 -u 2 , λ = λ T (u 1 | Γ -u 2 | Γ ) = 0 (7)
The three fields formulation thus reads,

∀(v 1 , v 2 , µ) ∈ L 2 (Ω 1 ) × L 2 (Ω 2 ) × L, a 1 (δu k 1 , v 1 ) + a 2 (δu k 2 , v 2 ) + δu k 1 -δu k 2 , µ + v 1 -v 2 , λ = l k 1 (v 1 ) + l k 2 (v 2 ) + u k-1 2 -u k-1 1 , µ (8)
When introduced in formulation ( 8), the finite element interpolation leads to the following coupled linear system:

  M 1 0 C T 1 0 M 2 C T 2 C 1 C 2 0     q k 1 q k 2 Λ   =   b k 1 b k 2 -C 1 p 1 -C 2 p 2   (9)
where Λ, p 1 and p 2 are the dof vectors corresponding to the FE interpolation of λ, u k-1 1 and u k-1 2 respectively. C s denote rectangular signed boolean operators such that if a dof belong to Γ, his value is set to 1 and the corresponding dof on the other side of Γ is set to -1.

As such, Problem ( 9) is coupled and involves all the unknown and the data of the image in the entire region of interest Ω. It is thus obviously not solved directly. In order to rationalize the computational cost of FE-DIC, the problem is split into the following coupled equations:

q k s = M -1 s (b k s -C T s Λ) ∀s ( 10
) s C s q k s + p s = 0 (11)
The dual interface problem is then obtained by the condensation of Problem ( 9) on Interface Γ. It consists in substituting equations [START_REF] Fedele | Global 2D digital image correlation for motion estimation in a finite element framework: a variational formulation and a regularized, pyramidal, multi-grid implementation[END_REF] into equation [START_REF] Fehrenbach | A fast algorithm for image registration[END_REF]:

S d Λ = t k (12)
where the condensed right hand side t k and the so-called dual Shur complement operator S d read

S d = s C s M -1 s C T s ( 13
)
t k = s C s p s + M -1 s b k s (14)
Remark. It is easy to see that Problem [START_REF] Glowinski | First International Symposium on Domain Decomposition Methods for Partial Differential Equations[END_REF] remains in the same form when more than two subdomains are considered. Nevertheless, attention must be paid to the interface dof belonging to more than two subdomains (multiplicity greater than 2). If an interface node is linked to n subdomains, n -1 Lagrange multipliers are introduced in each dimension.

Iterative solution algorithm

Introducing a decomposition of the domain and writing the condensed problem at the interface is not sufficient to make this method parallelisable. More, operator S d is formally based on the inversion of the local correlation operators M s whose computation cost would be prohibitive. Obviously, it is actually neither computed nor assembled. Indeed, the condensed symmetric linear system [START_REF] Glowinski | First International Symposium on Domain Decomposition Methods for Partial Differential Equations[END_REF] is solved thanks to a preconditionned conjugate gradient (CG) solver. After initialisation r 0 = t k -S d Λ 0 and d 0 = P -1 r 0 , one iteration of the preconditionned CG reads:

α = r T i P -1 r i d T i S d d i (15)
Λ i+1 = Λ i + α d i (16) r i+1 = r i -α S d d i (17)
β = r T i+1 P -1 r i+1 r T i P -1 r i (18)
d i+1 = P -1 r i+1 + β d i ( 19
)
where P is a preconditionner whose choice is discussed later. As written in the above algorithm, the conjugate gradient only requires matrix-vector products between operator S d and vector d i . As S d is expressed as a sum on subdomains, this product corresponds to the following computation:

S d d i = s C s M -1 s C T s d i = s C s y s
where y s are the solutions of the set of local systems M s y s = C T s d i that correspond to local correlation problems with different right hand sides. These system inversions are thus computationally affordable and, above all, independent from one subdomain to the other. That means that parallel processing can be used when computing these local products. Since the right hand side t k is the sum of local contributions t k s , it can also be computed in parallel. Note that operators M s remains unchanged during correlation iterations k. They can be factorised once at the beginning. The stopping criterion of the CG is classically based on the norm of the residual r i+1 .

The resulting domain decomposition DIC method yields exacly the solution one would get with one single mesh and a classical FE-DIC. The metrological performances of this method are thus the same than a classical FE-DIC method with a unique subdomain [START_REF] Hild | Comparison of local and global approaches to digital image correlation[END_REF].

Image decomposition

However, as stated earlier, the computation of the local contributions to the left and right hand side of [START_REF] Glowinski | First International Symposium on Domain Decomposition Methods for Partial Differential Equations[END_REF] requires to manipulate and interpolate images f and g at non-interger pixels positions within subdomain Ω s . In order to reduce the computational burden associated to whole image interpolation, we propose also to subdivide images f and g in a set of rectangular subset images f s and g s surrounding domain Ω s as it is usually done in subset-based approaches [START_REF] Sutton | Image correlation for shape, motion and deformation measurements: Basic Concepts, Theory and Applications[END_REF]. More precisely, g s is interpolated at the position x + u(x). Since u is neither constant within subdomain Ω s nor zero mean, g s is defined up to the average rigid body translation estimated by a Fast Fourier Transform, and it is a bit larger than f s . As in subset based methods, it may appear that the displacement u s is such that x + u s is larger than initially estimated subset image g s . A simple test is done before each interpolation step in order to avoid such a problem by enlarging or shifting the corresponding deformed subset image g s .

Ω Ω Ω 1 2 1 f (x) 2 f (x) f(x) initial FE-DIC pb independant FE-DIC sub-problems iterative solver
Therefore the computational effort associated to both finite elements (system inversions) and image interpolations are local, and thus parallelisable. A scheme of the image and mesh decomposition is provided in Figure 1. Not only these computations can be made on parallel computers, but the complexity and memory requirements of both system inversion and image interpolation are such that this method may be of interest in the context of high resolution analyses even in a sequential implementation.

Initialisation

Initialisation of the DIC iterations. Since a Newton-like solver is used for solving the correlation problem, a good initial value u 0 (x) of the unknwon displacement is required. In practice, the algorithm is here initialised by the solution to DIC problems in each subdomain independently. These problems correspond to classical FE-DIC analyses performed independently on each subdomain. Note that these problems may themselves be also initialised by multigrid and coarse graining techniques [START_REF] Hild | Measuring stress intensity factors with a camera: Integrated digital image correlation (I-DIC)[END_REF][START_REF] Gomes Perini | A multigrid PGD-based algorithm for volumetric displacement fields measurements[END_REF] in order to avoid local minimas. This provides a good initial estimate of the unknown displacement as shown in the examples below. The displacement of the inner nodes in the subdomains are almost converged. The correction computed by the proposed domain decomposition method mostly concerns nodes near the interface. Therefore, only a few extra iterations on the correlation problem are required using the proposed domain decomposition method, to bond the subdomains solutions.

Initialisation of the CG iterations. One conjugate gradient solver is used at each correlation iteration k. Each time, an initial value of Λ is also required. In practice, the first one (iteration k = 0) is set to zero, and the next ones (k > 0) are set to the last value of the previous CG resolution. This is a very naive but effective way to reduce the number of CG iterations, as shown in the examples below.

Convergence acceleration

In addition to a good initialisation, there are obviously many more efficient ways to speed-up a CG resolution. For instance, one can use preconditionning [START_REF] Glowinski | First International Symposium on Domain Decomposition Methods for Partial Differential Equations[END_REF][START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF]. The idea is to improve the condition number of the initial operator S d on which depend the convergence rate of the CG. In practice, a classical choice consists in constructing an operator P that is an approximation of S d but easier to invert. The preconditionned conjuguate gradient solver [START_REF] Gosselet | Total and selective reuse of Krylov subspaces for the resolution of sequences of nonlinear structural problems[END_REF][START_REF] Hild | Measuring stress intensity factors with a camera: Integrated digital image correlation (I-DIC)[END_REF][START_REF] Hild | Comparison of local and global approaches to digital image correlation[END_REF][START_REF] Hild | Calibration of constitutive models of steel beams subject to local buckling by using digital image correlation[END_REF][START_REF] Horn | Determining optical flow[END_REF] consists in solving the problem P -1 S d Λ = P -1 t k s . In the following, we propose to use three basic preconditionners, listed in an increasing order of efficiency:

1. the double diagonal preconditionner which is both cheap to build and to use in the preconditionned conjugate gradient (PCG), since the resulting operator is diagonal:

P dd = s diag(C s diag(M s ) -1 C T s )
2. the diagonal preconditionner which is a bit more expensive to assemble since local operators M s need to be inverted. Its cost is however affordable since it has to be done once and comparable to the cost of the so-called Dirichlet preconditionner of FETI [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF]:

P d = s diag(C s M -1 s C T s )
3. the quasi-diagonal preconditionner. Its assembly is almost as cheap as the double diagonal, but it is not exactly diagonal.

P qd = s C s diag(M s ) -1 C T s
In practice, C s C T s are band matrices and their bandwidth, a priori, depends on the numbering of the Lagrange multipliers and the dimension of the problem (2D in this case). However, it is possible to build a tridiagonal operator by using an appropriate numbering of the Lagrange multipliers. Two rules have to be considered: [START_REF] Bay | Digital volume correlation: three-dimensional strain mapping using x-ray tomography[END_REF] an interface dof in one subdomain must not be associated to more than 2 Lagrange multipliers (which is possible even for multiplicity higher than 2) and ( 2) the Lagrange multipliers must be sorted by physical dimension (first x dofs, then y dofs). It can be noticed that such a method will also work in 3D DVC. Furthermore, the number of non-zero extra diagonal values is equal to the sum of the multiplicity of each interface dof minus the number of interface dof, which, at least in 2D DIC, is very low. However, it is necessary to perform a Cholesky factorisation, for instance, once and for all the correlation iterations. Then, at each CG iteration, a forward and a backward substitution is performed. The associated extra cost is therefore very reasonable.

Remarks Another possibility would consists in computing a lumped correlation operator as what is done in transient dynamics [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF]. But according to our tests, and thanks to the short varying evolution of the image gradients, our implementation of such a technique was not efficient for preconditionning a domain-decomposition DIC problem. Beyond preconditionning, another way to accelerate the CG resolution would be to reuse Krylov subspaces [START_REF] Gosselet | Total and selective reuse of Krylov subspaces for the resolution of sequences of nonlinear structural problems[END_REF] which would be particularly efficient in this case since the operator is constant and the right hand side changes little.

Finally, a scheme of the domain decomposition DIC algorithm is given in Figure 2. 3 Application to a four point bending test

Find the max. on of

Fourier Initialisation

Compute

Experimental setup

The domain decomposition method previously described is now applied to the analysis of a four point bending test performed on an open hole polymethyl methacrylate (PMMA) specimen. The latter is a parallelepipedic coupon of dimension 200 × 20 × 5mm with a 8mm hole drilled in the center. A black and white speckle pattern is sprayed onto the surface of the specimen. A high resolution CCD camera (29 megapixels) has been used to capture images (4384 × 6576 pixels). An unstructured finite element mesh (including both linear triangles and quadrangles of 30 pixels width in average) is adjusted on the image, see Figure 3.

Initialisation and DD solution

A first decomposition into 179 subdomains1 is considered (Fig. 3). The average subdomain size is 160 pixels. The analysis is performed with a stopping tolerance for correlation and CG iterations of 10 -3 . First, the initialisation is performed on each subdomain indepandantly. This stage is highly parallelisable since no communications are required between subdomains. This initial displacement solution is plotted in Figure 4 (top) with an amplification factor of 50. The colorscale corresponds to the horizontal component of the displacement field. It can be noticed that this low-cost solution is a good initialisation for the proposed domain decomposition solver. A zoom in the top right is provided to show that the displacement continuity is actually not verified between subdomains. It seems clear that the extra effort, provided by the domain decomposition solver, to make the displacement continuous between subdomains, is rather slight. In this case, it took 8 extra iterations to converge. The resulting displacement is plotted in Figure 3 (bottom). The zoom on the bottom right shows that the displacement is then continuous.

Convergence acceleration

The variation of the norm of the residual r i as a function of the CG iteration i for the third correlation iteration (k = 3) is plotted in Figure 5 for the five following solvers: • iCG. A standard conjuguate gradient without preconditionning (P set to identity) where the solution Λ of the previous correlation iteration (k -1) is used to initialize the CG at correlation iteration k.

• iCG Pdd. The same as iCG but with the double diagonal preconditionner P dd .

• iCG Pd. The same as iCG but with the diagonal preconditionner P d .

• iCG Pqd. The same as iCG but with the quasi diagonal preconditionner P qd .

It took around 65 iterations for the classic CG to reach convergence. By merely initializing the CG with the previous solution vector Λ, it can be seen that the number of iterations is almost divided by two since the initial estimate of Λ is better by one decade in comparison to Λ = 0. However, this technique is not actually an acceleration technique since the convergence rate is the same than without initialisation. When using the proposed preconditionners, it can be noticed that the convergence rate itself is also improved.

The number of iterations required by the CG to reach the convergence criterion is reported in Figure 6 as a function of the chosen solver for the eight correlation iterations. The total number of CG iteration used by the domain decomposition solver can be drastically reduced when using both initialisation and preconditionning. According to this example, all of the three proposed preconditionners are efficient in reducing the number of iterations. The quasi-diagonal preconditionner seems the more efficient in this case.

In order to illustrate the sparsity of the quasi-diagonal preconditionner, the non-zero values of P qd are potted in black in Figure 7 (left). The matrix is actually tridiagonal, but the number of non-zero extra diagonal elements is very low. It is equal to 2 × 568 (by symmetry) whereas the number interface dof is 5238. It seems that P qd is a cheap but rather good approximation of the dual Schur complement S d even if the latter is not banded at all. Indeed, its non-zero values are plotted in Figure 7 (right) for comparison purposes.

Metrological and computational performances

The measurement uncertainty is finally analysed in Figure 8. As mentionned above, the domain decomposition method actually converges toward the monolithic solution (one single domain). Therefore, the metrological performances are the same. However, an a priori performance analysis is carried out in order to illustrate the effect of bonding the subdomains on measurement uncertainties. The standard displacement uncertainty is estimated by subjecting the reference image to a sub-pixel rigid body shift with a Fast Fourier Transform [START_REF] Bornert | Assessment of digital image correlation measurement errors: methodology and results[END_REF]. Then the measurement is performed between the reference and shifted images. The standard deviation of the error between measured and prescribed displacement is plotted as a function of the average subdomain size in pixels. This standard uncertainty is estimated for the initial displacement estimate (uncoupled subdomains), the solution to the proposed domain decomposition solver (coupled subdomains) and a standard single domain FE-DIC method used as a reference. As expected, the uncertainty of the solution to the domain decomposition solver is equal to the one with a single subdomain solution, regardless of the subdomains size. On the other hand, the solution computed independently on each subdomain (which is used as initialisation), is subject to larger uncertainties, as subset-based DIC methods do [START_REF] Hild | Comparison of local and global approaches to digital image correlation[END_REF], since continuity is not prescribed between subdomains.

At this stage, only a sequential implementation of the method was performed in Matlab. It is subsequently not possible to provide a relevant estimate of the speed-up. However, the average normalised CPU time required to compute one iteration of local FE-DIC (including image interpolation and FE system resolution) is plotted Figure 8: A priori performance analysis: Variation of the standard displacement uncertainty with respect to graylevel interpolation for the initial displacement estimate (uncoupled subdomains), the solution to the domain decomposition solver (coupled subdomains) and a standard single domain FE-DIC method (red dashed line). The gray zone corresponds to subdomain size smaller than the average element size of the mesh.

in Figure 9 versus average sudomain size. It appears that the numerical complexity seems to be a bit less that O(N 2 ). In other words, the CPU time of one iteration is divided by almost 100 when the size of the subdomains is divided by 10. With a sequential implementation of the method, 65% of the CPU time is devoted to compute the initialisation (first bloc of Figure 2), whereas the remaining 35% are used to perform the 8 extra iterations with the DD solver (second bloc of Figure 2). Since the first 65% correspond to computations that are independent by subdomain (it requires no data exchange between subdomains), this first initialisation is thus highly parallelizable and its computational cost, with a parallel implementation, is expected to be divided by the number of subdomains. The second part (35%) is also highly parallelizable, but requires some data exchange between subdomains. The results presented in Figures 8 and9 exemplify that the proposed approach combines the advantages of FE-based DIC and subset-based DIC methods.

Summary and perspectives

A non-overlapping dual domain decomposition method is proposed for finite element based digital image correlation. It consists in splitting the resolution of a global FE-DIC problem into a set of local FE-DIC sub-problems, based on smaller independent submeshes and subset images. The initial problem is condensed on the interface. A well initialised conjuguate gradient solver combined with three possible cheap preconditionners is presented. The method is used to analyse a 2D DIC problem with real high resolution images. The fact that the method has the same metrological performances as a classic FE-DIC method is illustrated. The first CPU estimations gives a glimpse of the potential of a parallel implementation of the method. It should offer the scalability of subset-based methods to global approaches to DIC. The use of such a domain decomposition method should be of great interest for the correlation of high definition digital volume images that are out of reach with standard FE DVC approaches [START_REF] Leclerc | Voxel-scale digital volume correlation[END_REF][START_REF] Leclerc | Digital volume correlation: What are the limits to the spatial resolution[END_REF].

Essentially developed for high performance computing considerations, such a method may have also other applications. Indeed, this partitionning method is also a good candidate for measuring continuous displacements in the case of a multiple-camera instrumentation of a single specimen [START_REF] Orteu | Multiple-camera instrumentation of a single point incremental forming process pilot for shape and 3D displacement measurements: Methodology and results[END_REF][START_REF] Passieux | Multiscale displacement field measurement using digital image correlation: Application to the identification of elastic properties[END_REF]. This kind of instrumentation yields, in essence, a set of independent images, that the proposed domain decomposition method would manage conveniently and efficiently. The extention to stereo digital image correlation could subsequently be a good candidate for measurements on large structures. This will also be a good tool for the high resolution analysis of complex microstructures present, for instance, in composite materials, polycrystaline metals or biomaterials [START_REF] Mortazavi | High resolution measurement of internal fullfield displacements and strains using global spectral digital volume correlation[END_REF][START_REF] Lecomte-Grosbras | Three-dimensional investigation of free-edge effects in laminate composites using x-ray tomography and digital volume correlation[END_REF][START_REF] Gomes Perini | A multigrid PGD-based algorithm for volumetric displacement fields measurements[END_REF]. Such a coupling method could also be used to couple finite elements with other interpolations types, like, for instance, X-FEM [START_REF] Réthoré | Shear-band capturing using a multiscale extended digital image correlation technique[END_REF] and Williams' series expansions [START_REF] Roux | Stress intensity factor measurements from digital image correlation: post-processing and integrated approaches[END_REF], for the estimation of stress intensity factors of curvilinear cracks, following what is done in simulation [START_REF] Passieux | Local/global non-intrusive crack propagation simulation using a multigrid X-FEM solver[END_REF].

Figure 1 :

 1 Figure 1: Global image f and domain Ω decomposition in two subdomains Ω s and two subset images f s . A set of subset of image g are also built. Note that subdomains and submeshes are not overlapping, but subset images are.

Figure 2 :

 2 Figure2: A scheme highlighting the link between the proposed approach, the conventional FE-DIC strategy and piecewise constant Fourier based approach used for initialisation.

Figure 3 :

 3 Figure 3: Example of a four point bending test. Image and domain decomposition in the case of 179 subdomains. The overlapping subset images f s are in yellow dashed line. The finite elements belonging to a same subdomain are plotted with the same color. A zoom is provided on the right.

Figure 4 :

 4 Figure 4: Example of a four point bending test: displacement solution in pixels (amplification factor 50): (top) FE-DIC solution computed independently on each subdomain (bloc Initialis. in Figure 2) and (bottom) continuous solution computed by the domain decomposition method (bloc Correlation iterations in Figure 2)

Figure 5 :

 5 Figure 5: Variation of the norm of the residual r i as a function of the iteration number i for the third correlation iteration k = 3 as a function of the choosen solver

Figure 6 :

 6 Figure 6: Number of iteration required to satisfy the convergence indicator as a function of the correlation iteration k for 5 different solvers.

Figure 7 :

 7 Figure 7: Non-zero elements of the quasi-diagonal preconditionner P qd (left) and comparison to those of the dual Schur complement S d (right).

Figure 9 :

 9 Figure 9: Insight in CPU performance: average CPU time required to compute one iteration of local FE-DIC normalised by the cost of one iteration of the global FE-DIC problem.

Ω is devided into 180 (6 × 30) regions. But no subdomain is associated to the region of the hole since it is not meshed.


Johann Lavaud
email: johann.lavaud@uni-konstanz.de

Robert F Strzepek

Peter G Kroth

Photoprotection capacity differs among diatoms: Possible consequences on the spatial distribution of diatoms related to fluctuations in the underwater light climate

In the present report, we show a fundamental difference between diatom species from different marine habitats in their ability to cope with changes in irradiance.

Estuarine species show a higher and more flexible capacity for photoprotection than oceanic and coastal species, and when exposed to excess light, the impairment of their photosynthetic capacity due to photoinhibition was reduced. This resulted in maintenance of growth in a fluctuating light regime, conferring the estuarine species an adaptive advantage. The ability of diatoms, and to a larger extent other phytoplankton, to occupy a wide range of ecological niches depends critically on their capacity to exploit the differences in underwater light climate. These results may explain how diatoms adapt to the challenge of maintaining optimal photosynthetic production in turbulent waters where the rate of light changes is high.

Introduction

Phytoplankton are passively transported within the water column by turbulent mixing and must cope with highly variable light. Unlike nutrient availability and temperature, fluctuations in light intensity can show very high frequency (few seconds) coupled with high amplitude (from darkness to full sunlight) [START_REF] Lewis | Turbulent motions may control phytoplankton photosynthesis in the upper ocean[END_REF][START_REF] Macintyre | The effect of water motion on short-term rates of photosynthesis by marine phytoplankton[END_REF]. Turbulence and light fluctuations have been identified as key parameters in the structure of phytoplankton communities and have been integrated in recent ecological modeling for species dynamics and bloom development [START_REF] Huisman | Species dynamics in phytoplankton blooms: incomplete mixing and competition for light[END_REF][START_REF] Lichtman | Competition of phytoplankton under fluctuating light[END_REF][START_REF] Huisman | Changes in turbulent mixing shift competition for light between phytoplankton species[END_REF].

Fluctuation in the underwater light climate strongly influences the photosynthetic productivity of phytoplankton [START_REF] Lewis | Turbulent motions may control phytoplankton photosynthesis in the upper ocean[END_REF][START_REF] Cullen | The kinetics of photoadaptation in the context of vertical mixing[END_REF]. Fluctuating light in combination with excessive light is potentially harmful for photosynthesis as it can cause a decrease in productivity and fitness [START_REF] Long | Photoinhibition of photosynthesis in nature[END_REF][START_REF] Külheim | Rapid regulation of light harvesting and plant fitness in the field[END_REF].

Photosynthetic organisms attenuate these undesirable effects by safely dissipating excess energy [START_REF] Niyogi | Safety valves for photosynthesis[END_REF]. The extent and regulation of this photoprotection vary amongst phytoplankton groups (Wagner et al. 2006 and citations within), diatoms being generally considered to have a better ability to cope with fluctuating light [START_REF] Mitrovic | Unusual allometry between in situ growth of freshwater phytoplankton under static and fluctuating light environments: possible implications for dominance[END_REF][START_REF] Wagner | Balancing the energy flow from captured light to biomass under fluctuating light conditions[END_REF].

Diatoms contribute highly to Earth's primary production (25 %) [START_REF] Sarthou | Growth physiology and fate of diatoms in the ocean: a review[END_REF]. Their dominance in turbulent waters, where vertical mixing exposes them to large fluctuations in irradiance [START_REF] Huisman | Changes in turbulent mixing shift competition for light between phytoplankton species[END_REF][START_REF] Tozzi | Historical climate change and ocean turbulence as selective agents for two key phytoplankton functional groups[END_REF], suggests an unusual photosynthetic flexibility. Here we report a habitat-related difference amongst diatom species in their ability to cope with light. Estuarine species show a higher and more flexible capacity for photoprotection than oceanic and coastal species, conferring them an adaptive advantage in a turbulent environment. These results provide an ecophysiological explanation for the spatial distribution of planktonic diatoms as a function of the underwater light climate.

Materials and methods

Culturing

The diatom strains used were Phaeodactylum tricornutum (clone CCAP1052, isolated near Plymouth, England), Skeletonema costatum (clone 'Banyuls'-Laboratoire Arago Collection, Banyuls-sur-Mer, Mediterranean Sea, France), Thalassiosira weissflogii (CCMP1336, estuary Long Island, New York, USA), T. pseudonana (CCMP1335, Moriches Bay, Long Island), and T. oceanica (CCMP1003, Sargasso Sea).

P. tricornutum and S. costatum were grown photoautotrophically at 18°C in continuously flushed airlifts of 300 mL sterile natural seawater F/2 medium. The Thalassiosira species were grown in Aquil medium. The light conditions were 45 µmol photons m -2 s -1 on a 16:8 h light (L):dark (D) cycle. For the growth rate experiments on a 1:1 h L:D 24 h cycle (45 µmol photons m -2 s -1 ), the cells were maintained in this regime for at least two weeks and then transferred to 450 µmol photons m -2 s -1 (same L:D cycle). All measurements were made during the exponential growth phase during the first week of growth in order to catch the immediate response of the cells to this switch in light regime. Specific growth rates, µ (d -1 ), were calculated from regression of natural logarithm of culture chlorophyll a (Chl a) concentration or fluorescence during the exponential growth phase of acclimated cultures.

Photosynthetic architecture and pigment analyses

Steady-state O 2 yield per flash (Y SS ) and the concentration of P 700 per Chl a were used as proxies for photosystems (PS) II and I, respectively and measured as described before (Lavaud et al. 2002a). P 700 was determined with an Aminco DW-2 spectrophotometer (American Instrument Co., USA). The O 2 yield per flash was measured with a home-made flash O 2 electrode (Lavaud et al. 2002b). PSII antenna size (1/I 1/2 of Y SS , the reciprocal of the half-saturating flash intensity of flash O 2 evolution saturation curves) was measured the same way while applying a set of neutral filters to vary the intensity of the flashes (Lavaud et al. 2002a). Xanthophyll analyses were performed by HPLC as described [START_REF] Arsalane | Influence of the pool size of the xanthophyll cycle on the effects of light stress in a diatom: Competition between photoprotection and photoinhibition[END_REF]). The xanthophyll de-epoxidation state (%) was calculated as (DT/DD+DT) x 100 where DD is diadinoxanthin, the epoxidized form and DT is diatoxanthin, the de-epoxidized form.

Photosynthetic measurements

Flash O 2 measurements (Y SS , number of active PSII), photosynthesis versus irradiance (P/E) curves and fluorescence yield measurements were performed on darkacclimated samples (15 min). Cells were brought to a final concentration of 10 µg Chl a mL -1 following gentle centrifugation (5 min, 5000 rpm). Cells were resuspended in F/2 medium and allowed to recover for 1 h in ambient light. P/E curves were determined at 18°C with a calibrated Clark O 2 electrode (model DW1, Hansatech Ltd., England) by measuring O 2 evolution normalized to Chl a as a function of light intensity. White light of adjustable intensity was provided by a KL-1500 quartz iodine lamp (Schott, Mainz, Germany). A new sample was used for each irradiance. P max (maximum photosynthetic rate) and E k (light intensity for saturation of photosynthesis) were derived from P/E curves. Chl a fluorescence yield measurements were performed with a PAM-101/103 fluorometer (for P. tricornutum and S. costatum) and a PhytoPAM with an emitterdetector-cuvette assembly (ED101) (for the Thalassiosira species) (Heinz Walz, Effeltrich, Germany). The maximum photosynthetic efficiency of PSII (F v /F m ) was determined by inducing a transient closure of the PSII reaction centers with a saturating light pulse (600 ms, 4000 µmol photons m -2 s -1 ). The capacity for energy dissipation was calculated as NPQ = F m /F m ´-1. For a comparison with the potential for photochemistry qP [START_REF] Ruban | The super-excess energy dissipation in diatom algae: comparative analysis with higher plants[END_REF].

(= (F m -F s )/(F m ´-F o ´)) it was calculated as qN = (F m -F m ´)/(F m -F o ´) (see
Photoinhibition kinetics were performed by exposing the cells to different irradiances (5 min illumination duration, a new sample was used for each irradiance step) or illumination times (at 2000 µmol photons m -2 s -1 ) directly in the DW1-O 2 electrode chamber and measuring the decrease in the number of active PSII reaction centers (Y SS ) by regular sampling. Using the same protocol, the extent of the PSII electron cycle was determined by measuring the O 2 deficit at the beginning of an O 2 flash sequence (see Lavaud et al. 2002b). The O 2 deficit per PSII reaction center in a 20-flash series was quantified by comparison of the O 2 deficit in illuminated ( L ) and non-illuminated ( D ) samples as follows:

O 2 deficit = [{(20xY SS L )-(Σ (Y 1…20 ) L )} -{(20xY SS D )-( Σ (Y 1…20 ) D )}] / Y SS D .

Results and discussion

The key to our findings came from the measurement of non-photochemical chlorophyll a fluorescence quenching (NPQ), a parameter to evaluate photoprotection.

The NPQ process allows the photosystem II (PSII) light-harvesting complex (LHC) antenna to safely dissipate light excitation energy absorbed in excess of photosynthetic capacity during a sudden increase in irradiance [START_REF] Holt | Toward an understanding of the mechanism of the nonphotochemical quenching in green plants[END_REF]). NPQ is one of the most important mechanisms for rapid (seconds to minutes) regulation of photochemistry.

Consequently, it is one of the first lines of defense that diatoms [START_REF] Arsalane | Influence of the pool size of the xanthophyll cycle on the effects of light stress in a diatom: Competition between photoprotection and photoinhibition[END_REF][START_REF] Lavaud | General features of photoprotection by energy dissipation in planktonic diatoms (Bacillariophyceae)[END_REF]) use to attenuate photoinhibitory oxidative damage generated by light stress [START_REF] Niyogi | Safety valves for photosynthesis[END_REF]. Five plankton diatom strains were chosen as representative of marine habitats (estuary, bay, coast and open ocean; see the Materials and methods section for details) characterized by different rates of water mixing and hence underwater light climates. A light gradient exists in the marine environment that ranges from dynamic and turbid estuarine waters, characterized by light fluctuations with high amplitude and frequency (partly due to an active tidal hydrodynamism), to calm and clear oceanic waters characterized by a more stable average irradiance (illustrated by a seasonal stratification); coastal/bay waters being intermediary [START_REF] Macintyre | The effect of water motion on short-term rates of photosynthesis by marine phytoplankton[END_REF]. To test if the photoprotective capacity of the diatoms reflected their original habitat, NPQ was measured in laboratory isolates exposed to irradiances ranging from darkness to full sunlight (Fig. 1). The estuarine species showed a higher NPQ (up to 3 to 5 fold) than coastal and oceanic species. Within the same genus (Thalassiosira), photoprotective potential clearly reflected the geographical origin of the isolates: NPQ was highest in the estuarine species and lowest in the oceanic species, while T. pseudonana, originating from a bay, showed intermediate NPQ induction. Measurements of NPQ in the field suggest that this process is essential for maintaining the photosynthetic activity of diatoms in a variable light climate [START_REF] Kashino | Strategies of phytoplankton to perform effective photosynthesis in the North Water[END_REF][START_REF] Fujiki | Photoprotective response of xanthophyll pigments during phytoplankton blooms in Sagami Bay, Japan[END_REF], and that it may explain the spatial and temporal distribution of diatom species [START_REF] Meyer | Xanthophyll cycling in Phaeocystis globosa and Thalassiosira sp.: a possible mechanism for species succession[END_REF][START_REF] Serodio | Non-photochemical quenching of chlorophyll fluorescence and operation of the xanthophyll cycle in estuarine microphytobenthos[END_REF]. We therefore examined why an estuarine species is better suited to fluctuating light and what consequences this ability may have on the maintenance of photosynthetic activity and growth.

Phaeodactylum tricornutum and Skeletonema costatum were chosen as model species adapted to high and low rates of light fluctuation, respectively. P. tricornutum is primarily both tychopelagic and estuarine while S. costatum is often pelagic and coastal.

The rate of light change is higher in estuary than in coastal waters because of the higher vertical diffusivity and the higher light attenuation [START_REF] Macintyre | The effect of water motion on short-term rates of photosynthesis by marine phytoplankton[END_REF]. This choice was justified by the fact that the species had comparable photosynthetic capacity (Table 1) in our growth conditions. This similarity allowed us to rule out an adaptive response of the architecture of the photosynthetic apparatus to iron availability, as has been observed in oceanic species [START_REF] Strzepek | Photosynthetic architecture differs in coastal and oceanic diatoms[END_REF]. Even though we did not perform indepth analysis for the three Thalassiosira species, we observed that in the light conditions we used to grow the cells (45 µmol photons m -2 s -1 , 16:8 h light (L):dark (D)), the three Thalassiosira species showed a similar F v /F m , rate for linear electron transport (ETR = 11  3 µmol m -2 s -1 for all the species) and growth rate as P. tricornutum and S. costatum.

Additionally, a compilation of several growth and photosynthetic parameters from four of the species used here (P. tricornutum, S. costatum, T. weissflogii and T. pseudonana)

showed that there is no clear evidence for systematic differences, especially due to size/shape and pigment absorption, when they were grown under similar light and nutrient conditions [START_REF] Macintyre | Photoacclimation of photosynthesis irradiance response curves and photosynthetic pigments in microalgae and cyanobacteria[END_REF]. Hence, the fact that the global photosynthetic capacity was similar for P. tricornutum and S. costatum on one side, and that some essential photosynthetic features were also similar for the three Thalassiosira species on the other side is of first importance for the comparison of NPQ capacity (Fig. 1) and for deeper characterization of photoprotection ability under different light conditions. It is at least reasonable to argue that the physiological state vs growth irradiance was similar for all the species under the conditions used here and that this allows a reliable comparison of their ability to cope with changes in light intensity. The differences in the response of the species are then likely to be more adaptive (i. e. genetic) than acclimative (i. e. phenotypic), even if the acclimative aspect can not be totally ruled out.

The reaction of P. tricornutum and S. costatum to a gradual increase in incident light was dramatically different. NPQ was induced in S. costatum at low light intensities (typically just over the growth intensity) while in P. tricornutum NPQ was induced only when light intensities were 5-to 6-fold higher than the growth light intensity (inset Fig. 1). The most striking difference between the species was the maximal extent of NPQ which was 2.7 times higher in P. tricornutum (Fig. 1). Both the delay in induction and the higher extent of NPQ in P. tricornutum are characteristic of high light acclimated diatoms (J. Lavaud, unpubl. results). NPQ is dependent on the xanthophyll cycle [START_REF] Holt | Toward an understanding of the mechanism of the nonphotochemical quenching in green plants[END_REF]), an enzymatic interconversion of epoxidized to de-epoxidized xanthophylls regulated in response to changes in irradiance [START_REF] Lavaud | General features of photoprotection by energy dissipation in planktonic diatoms (Bacillariophyceae)[END_REF]). However, our results suggest that the difference in NPQ was not due to a difference in the xanthophyll cycle operation since de-epoxidation levels were even higher in S. costatum (Table 1). The lower NPQ in S. costatum may be due to a different LHC antenna polypeptide composition and organization, and/or to fewer xanthophylls involved in the NPQ process [START_REF] Lavaud | General features of photoprotection by energy dissipation in planktonic diatoms (Bacillariophyceae)[END_REF][START_REF] Ruban | The super-excess energy dissipation in diatom algae: comparative analysis with higher plants[END_REF].

In comparison to S. costatum, P. tricornutum was able to fine tune the dissipation of light excitation energy (qN, similar to NPQ) relative to light-dependent photochemical activity (qP). The relationship between the two parameters was sigmoidal in P. tricornutum with a rapid phase of qN development/dissipation within a very narrow range of qP (Fig. 2, open circles) that may be due to the fast activation of the xanthophyll cycle enzymes [START_REF] Jakob | Unusual pH-dependence of diadinoxanthin deepoxidase activation causes chlororespiratory induced accumulation of diatoxanthin in the diatom Phaeodactylum tricornutum[END_REF][START_REF] Goss | The importance of a highly active and delta-pH-regulated diatoxanthin epoxidase for the regulation of the PS II antenna function in diadinoxanthin cycle containing algae[END_REF]). In S. costatum, the relationship between qN and qP was linear (Fig. 2, open triangles) and lacked the rapid change around E k observed in P. tricornutum. This rapid change in the distribution of photosynthetic (qP) vs photoprotective (qN) allocation is of importance for an organism that is rapidly mixed through a light gradient like in estuarine ecosystems. P. tricornutum out-performed S. costatum not only in having a higher ability for NPQ in high light but also by being able to keep NPQ low up to a high qP value of 0.6 meaning that 40 % of the PSII reaction centers are in a reduced state (also highlighted in [START_REF] Ruban | The super-excess energy dissipation in diatom algae: comparative analysis with higher plants[END_REF]).

Indeed, if NPQ induced near the surface is not dissipated rapidly and nearly completely when the cells are mixed back to the bottom of the water column (i. e., to low irradiances), their light-harvesting ability and consequently their quantum yield for photosynthesis could be reduced. Hence, in P. tricornutum, the sigmoidal relationship between photosynthetic vs photoprotective energy allocation is consistent with adaptation to low mean irradiance punctuated by episodic exposure to very high irradiance that characterizes the thychopelagic/estuarine niche. On the other hand the linear relationship observed in S. costatum is more consistent to adaptation to a more stable average irradiance characterized by slower light fluctuations, as is encountered in pelagic/coastal ecosystems. It is noteworthy that the relationship qN/qP is also linear in Arabidopsis thaliana [START_REF] Ruban | The super-excess energy dissipation in diatom algae: comparative analysis with higher plants[END_REF]).

The divergence in photoprotective versus photochemical regulation observed in P. tricornutum was dependent on other rapid processes. Specifically, energy dissipation via the PSII electron cycle is an important alternative photoprotective mechanism in diatoms (Lavaud et al. 2002b). Measurements of the PSII electron cycle provided additional evidence of the ability of the estuarine species to rapidly and flexibly cope with changes in light excitation energy. PSII electron cycling varied in the same sigmoidal way with photochemistry as NPQ in P. tricornutum but was not detectable in S. costatum (Fig. 2, closed circles and triangles, respectively). These results illustrate how rapid regulatory processes work cooperatively as an adaptive response to changing light intensity. At irradiances up to 3 fold the growth irradiance, P. tricornutum only dissipates excess excitation energy via the PSII electron cycle. This feature probably explains the ability of P. tricornutum to keep NPQ low even if a high proportion of PSII are in a reduced state [START_REF] Ruban | The super-excess energy dissipation in diatom algae: comparative analysis with higher plants[END_REF]). At higher irradiances, NPQ operates in concert with the PSII electron cycling. In S. costatum such a relationship between these mechanisms appears not to exist.

An enhanced ability to cope with excess light exposure protects the photosynthetic apparatus from photoinhibitory damage [START_REF] Neale | Algal photoinhibition and photosynthesis in the aquatic environment[END_REF][START_REF] Demmig-Adams | Photoprotection and other responses of plants to high light stress[END_REF][START_REF] Long | Photoinhibition of photosynthesis in nature[END_REF] as shown previously specifically in diatoms [START_REF] Arsalane | Influence of the pool size of the xanthophyll cycle on the effects of light stress in a diatom: Competition between photoprotection and photoinhibition[END_REF]Lavaud et al. 2002a;[START_REF] Lavaud | General features of photoprotection by energy dissipation in planktonic diatoms (Bacillariophyceae)[END_REF]). The extent of photoinhibition was measured during short light exposures (5 min) of increasing irradiance (Fig. 3A). In both species, photoinhibition began at E k x 3-4 (see Table 1).

Above 800 µmol photons m -2 s -1 a difference between the two species was apparent: photoinhibition reached a maximum (of about 20 %) for P. tricornutum (closed circles) while in S. costatum (open triangles) it increased gradually to nearly 40 % at a light intensity comparable to full sunlight. After 15 min at this light intensity, S. costatum was more photoinhibited than P. tricornutum (Fig. 3B) but the difference started to decrease after 45 min light (only 10 %). Most of the photoinhibitory damage took place during the first 5-10 min of light stress and its extent was higher in S. costatum (by about 20 %) [START_REF] Lavaud | General features of photoprotection by energy dissipation in planktonic diatoms (Bacillariophyceae)[END_REF], thus highlighting the importance of rapid regulatory mechanisms in attenuating the deleterious effects of light stress on photosynthetic activity. Interestingly, the first-order reaction time constants of energy dissipation were shorter than the estimated time (up to a few seconds) needed for the irradiance to double in estuarine waters [START_REF] Macintyre | The effect of water motion on short-term rates of photosynthesis by marine phytoplankton[END_REF] suggesting that the rapidity and flexibility of these mechanisms is well suited for regulating the photosynthetic activity of diatoms living in a turbulent aquatic environment. To distinguish between the energy dissipated via NPQ and the PSII electron cycle, photoinhibition was measured in cells of P. tricornutum exhibiting a higher capacity for NPQ (up to 9-10 obtained by specific culture conditions) but a comparable ability for the PSII electron cycle (Lavaud et al. 2002a;Lavaud et al. 2002b). These cells were less photoinhibited during both a short light exposure (Fig. 3A) and under prolonged light stress (Fig. 3B) than cells with lower NPQ. For a short light exposure, light intensity ≥ 1500 µmol photons m -2 s -1 were needed to observe photoinhibition. When the light stress was prolonged (45 min) under a maximal irradiance (Fig. 3B), the difference between low and high NPQ cells started to be significant (20 % instead of less than 10 %). These results provide evidence that dissipatory pathways protect the photosynthetic apparatus not only during short variations in light intensity (few minutes) but also during a more prolonged exposure to excess light (about 1 hour), and that NPQ plays an important role in the global photoprotection in these light conditions.

In order to probe the response of cells to both a change in the light regime and intensity over longer time scales (generations), we grew P. tricornutum and S. costatum under a frequently changing light regime (1:1 h L:D) designed to mimic the light environment of cells spending 1 h at the surface of the water column followed by 1 h below the euphotic zone and so on. When both species were grown at low irradiance, they had similar F v /F m (0.70 ± 0.02) and growth rates (Table 1). At a 10-fold higher growth irradiance (450 µmol photons m -2 s -1 , just before the start of photoinhibition in both species in the light conditions of Fig. 3A), their growth rates drastically decreased (Table 1), but while P. tricornutum was able to maintain growth as well as F v /F m (0. 70 ± 0.01), S. costatum failed to grow and its photosynthetic potential decreased (F v /F m = 0.56 ± 0.01). Similar results were obtained with the Thalassiosira species (decrease of 20-25 % in the growth rate and F v /F m in the high light 1:1 h L:D regime for T. oceanica compared to T. weissflogii). Even though the light regime used was unlike the underwater light climates typically reported for estuarine ecosystems [START_REF] Macintyre | The effect of water motion on short-term rates of photosynthesis by marine phytoplankton[END_REF], the difference between species in their photoadaptative response was striking and corroborated the measurements performed under shorter light exposures, particularly the fact that at the high irradiance used for this experiment, energy dissipation was higher in P. tricornutum than in S. costatum (see Fig. 2). It is then likely that the higher photosynthetic flexibility and photoprotection ability confer P. tricornutum with an advantage over generations when grown under a frequently changing light regime with periodic exposure to high light.

Our results suggest that the geographical and spatial distribution of diatoms may be influenced by their ability to cope with rapid fluctuations in light intensity that are dictated by the rate of mixing and light attenuation of waters. It has been previously

shown that the availability of iron drives differential adaptive responses and ensures the distribution of diatoms amongst coastal/estuarine and oceanic ecosystems [START_REF] Strzepek | Photosynthetic architecture differs in coastal and oceanic diatoms[END_REF]. It is also well known how iron as well as other nutrients (silica, nitrogen, phosphate) strongly influence the ecophysiology of diatoms, especially their photosynthetic physiology [START_REF] Geider | Response of the photosynthetic apparatus of Phaeodactylum tricornutum (Bacillariophyceae) to nitrate, phosphate or iron starvation[END_REF][START_REF] Sunda | Interrelated influence of iron, light and cell size on marine phytoplankton growth[END_REF][START_REF] Claquin | Uncoupling of silicon compared with carbon and nitrogen metabolism and the role of the cell cycle in continuous cultures of Thalassiosira pseudonana (Bacillariophyceae) under light, nitrogen, and phosphorus control[END_REF][START_REF] Staehr | Photoacclimation of four marine phytoplankton species to irradiance and nutrient availability[END_REF][START_REF] Sarthou | Growth physiology and fate of diatoms in the ocean: a review[END_REF]. All of our experiments were conducted with nutrient replete media. Still, it is possible that for the 1:1 h L:D experiment, the switch in the light regime could have changed the enzymatic uptake of nutrients that could feedback on photosynthetic capacity and ultimately growth. Nevertheless, it is rather difficult to evaluate the real impact of such a process on our measurements as most of the previous studies on the diatom species used here were conducted on nutrientlimited cultures [START_REF] Kolber | Effects of growth irradiance and nitrogen limitation on photosynthetic energy conversion in photosystem II[END_REF][START_REF] Geider | Response of the photosynthetic apparatus of Phaeodactylum tricornutum (Bacillariophyceae) to nitrate, phosphate or iron starvation[END_REF][START_REF] Mckay | Physiological and biochemical response of the photosynthetic apparatus of two marine diatoms to Fe stress[END_REF][START_REF] Strzepek | Photosynthetic architecture differs in coastal and oceanic diatoms[END_REF]. In this report, we further demonstrate that fluctuations in the light underwater climate may also explain the differential spatial distribution of diatom species. It confirms and provides an explanation for the previous results on the ecological niche occupancy as a function of light regime between species [START_REF] Sakshaug | Thalassiosira oceanica and T. pseudonana: two different photoadaptational responses[END_REF], between ecotypes [START_REF] Gallagher | Ecotypic differentiation in the marine diatom Skeletonema costatum: influence of the light intensity on the photosynthetic apparatus[END_REF], and can be extended to other important phytoplankton groups such as the cyanobacteria and prochlorophytes which also use NPQ to acclimate to the underwater light climate [START_REF] Bailey | The occurrence of rapidly reversible non-photochemical quenching of chlorophyll a fluorescence in cyanobacteria[END_REF]. Accordingly, our estuarine model species shows a higher and more rapid capacity for safely dissipating energy when irradiance is suddenly increased, thus ensuring a higher flexibility in photosynthetic regulation, and ultimately the maintenance of photosynthesis. Even though the fluctuating light regime we used (1:1 h L:D, low/high light) to probe the response of the species over generations is not very realistic, it allowed us to conclude that the estuarine species was better able to maintain its growth in these light conditions, and we suggest that this was due to its higher photosynthetic flexibility and photoprotection ability. This attribute may well explain both the ecological success of diatoms in turbulent fluctuating light environment where they often out-compete other phytoplankton groups [START_REF] Mitrovic | Unusual allometry between in situ growth of freshwater phytoplankton under static and fluctuating light environments: possible implications for dominance[END_REF][START_REF] Huisman | Changes in turbulent mixing shift competition for light between phytoplankton species[END_REF], and patterns of species succession related to turbulence [START_REF] Lichtman | Population and community responses of phytoplankton to fluctuating light[END_REF][START_REF] Huisman | Changes in turbulent mixing shift competition for light between phytoplankton species[END_REF]). Together with other recent studies [START_REF] Fietz | Acclimation of the diatom Stephanodiscus neoastraea and the cyanobacterium Planktothrix agardhii to simulated natural light fluctuations[END_REF][START_REF] Wagner | Balancing the energy flow from captured light to biomass under fluctuating light conditions[END_REF], our work also provides an ecophysiological basis for the improvement of predictive models for light competition between phytoplankton species in aquatic ecosystems [START_REF] Huisman | Species dynamics in phytoplankton blooms: incomplete mixing and competition for light[END_REF][START_REF] Lichtman | Competition of phytoplankton under fluctuating light[END_REF]. Our results, as well as previous ones [START_REF] Strzepek | Photosynthetic architecture differs in coastal and oceanic diatoms[END_REF], thus provide an explanation at a physiological level for the early 'Margalef's theory' [START_REF] Margalef | Life-forms of phytoplankton as survival alternatives in an unstable environment[END_REF] proposing that the competition amongst phytoplankton species could be related to the upper-ocean turbulence and the supply of nutrients. Additionally, protection from specific grazers can also be important in determining the spatial as well as temporal distribution of diatom species (see [START_REF] Smetacek | Diatoms and the ocean carbon cycle[END_REF]. It has been recently debated [START_REF] Falkowski | The evolution of modern phytoplankton[END_REF][START_REF] Tozzi | Historical climate change and ocean turbulence as selective agents for two key phytoplankton functional groups[END_REF]) how modifications in the upper-ocean turbulence, through its effect on the supply of nutrients, has shaped the functional and evolutionary ecology of diatoms and currently influences their geographical distribution. We propose that the present report strengthens and expands this theory by showing that the spatial distribution of diatoms might be influenced, via water turbulence, not only by the nutrient availability but additionally by the underwater light climate; the ability for nutrient uptake and the photoadaptative strategy probably influencing each other [START_REF] Kolber | Effects of growth irradiance and nitrogen limitation on photosynthetic energy conversion in photosystem II[END_REF][START_REF] Raven | Predictions of Mn and Fe use efficiencies of phototrophic growth as a function of light availability for growth and of C assimilation pathway[END_REF][START_REF] Sunda | Interrelated influence of iron, light and cell size on marine phytoplankton growth[END_REF][START_REF] Staehr | Photoacclimation of four marine phytoplankton species to irradiance and nutrient availability[END_REF] and ultimately defining the success of a given species/group of organisms in a given ecosystem/niche [START_REF] Strzepek | Photosynthetic architecture differs in coastal and oceanic diatoms[END_REF]. d -1 . Data (± SD) are the average of three to four independent measurements. All measurements were performed on cultures growing at 45 µmol photons m -2 s -1 , 16:8 h light (L):dark (D) except for µ. 1/I 1/2 of Y SS is a measurement of the PSII light-harvesting antenna size (the pigment, especially xanthophyll, content of the species is very similar [START_REF] Lavaud | General features of photoprotection by energy dissipation in planktonic diatoms (Bacillariophyceae)[END_REF]). F v /F m is the maximum photosynthetic efficiency of PSII, P max is the maximum photosynthetic rate, E k is the light intensity for saturation of photosynthesis, XDep is the level of xanthophyll de-epoxidation, µ is the growth rate. 2000 µmol photons m -2 s -1 is equivalent to full sunlight [START_REF] Long | Photoinhibition of photosynthesis in nature[END_REF].

Figure legends

Figure 1. Capacity for the dissipation of the light excitation energy (NPQ) versus irradiance in diatoms isolated from different ecosystems. Two thousand µmol photons m - 2 s -1 corresponds to full sunlight in nature [START_REF] Long | Photoinhibition of photosynthesis in nature[END_REF]. Data (± SD) are the average of three to four independent measurements. T im e of illum ination (m in ) B

Figure 2 .

 2 Figure 2. Relationship between the potential for photochemistry (0-min < qP < 1-max)

Figure 3 .

 3 Figure 3. (A) Extent of photoinhibition (estimated by the loss of PSII active centers as a

Figure

  Figure 1. Lavaud, Strzepek & Kroth

Table 1 .

 1 Photosynthetic properties of the two model diatom species (Phaeodactylum tricornutum and Skeletonema costatum) adapted to high and low rates of light fluctuation,

	respectively.
	Measure		P. tricornutum	S. costatum
	PSII:PSI		1.7:1	2.1:1
	1/I 1/2 of Y SS		0.32	0.33
	F v /F m		0.74  0.02	0.73  0.04
	P max		122.4  13.0	155.6  21.1
	E k		132  10	125  10
	XDep	E k	2.4  2.4	20  3.0
		2000 µmol photons m -2 s -1	27.4  0.2	40.2  2.8
	µ, 1:1 h L:D	45 µmol photons m -2 s -1	+ 1.03  0.08	+ 0.96  0.08
		450 µmol photons m -2 s -1	+ 0.48  0.04	-0.06  0.02
	Units are as follows: P

max , µmol O 2 mg Chl a -1 h -1 ; E k , µmol photons m -2 s -1 ; XDep, %; µ,

Acknowledgements

We thank G. Vétion and Dr. M.-J. Chrétiennot-Dinet for providing the strain of S. costatum and information on it, and B. Rousseau and M. Desnos for technical help. This work was supported by grants from the University of Konstanz (to J.L. and P.G.K.) and by the Foundation for Research Science and Technology (FRST) of New Zealand (to R.F.S.). We thank the two anonymous reviewers for their help in improving the manuscript.


domain: stat math info

text:
The figures show the laser doppler data (blue) and the appoximation (red) given by the learned representation. During learning with Algorithm 1 we included a constant waveform , see green baseline. Note that relation between spike occurence (red asterisk, bottom figures) and doppler response is only clearly visible during low spiking rates.

Results on LFP

Algorithm 1 produces a compact representation of the data, providing additional insight. four spikes were learned on the data, note that their amplitudes are normalized. Their shapes differ mainly in duration. Right: amplitudes of detected spikes (top) are decreasing with time. They appear strongly negatively correlated with their local spiking rate (defined as 1 / distance to previous spike [s]).

Processing LFP Data: Two Approaches

Hemodynamic response is learned using Algorithm 1, where we take as the previously detected spike latencies.

Time-frequency analysis reveals synchronization of spikes with 0.9-1 Hz rhythm in hemodynamics, probably related to respiration.

B) Continuous model: LFP data is composed of spike shapes repeated at different latencies with varying amplitudes plus some noise :

Difficulties:

Large dataset with several hundred epileptic spikes.

Spikes change over time (shape, amplitude, spiking rate).

Overlaps in spikes and hemodynamic responses.

3a Correlation based detection of waveforms, can be solved by sparse coding (e.g. LARS [2]) with all possible shifts of the .

3b Simple averaging if detected waveforms are not overlapping.

Otherwise overlaps are accounted for in linear system (discrete deconvolution), which may be ill-conditioned for strong overlaps. , which is initialized with e.g. an arbitrary spike in the data (2). Then all spikes are detected and their latencies and amplitudes are calculated (3a, details below). The are now used to update waveforms . Steps 3a and 3b are repeated until convergence. If the residual still contains much structure, the previous steps may be repeated with an additional waveform.

Local spiking rates (1/distance [s] to previous spike) are plotted on time-frequency map of laser doppler. During the three spike outbursts around 750, 1250, and 1900 seconds the spikes appear to synchronize with 0.9-1 Hz activity in the doppler. Around 750 seconds local spiking rates even clearly correspond to the subharmonics of this rhythm.

Hemodynamic response learned on laser doppler using Algorithm 1 with and the learned spike latencies. The presence of the oscillations with 7s period probably result from strong overlapping effects due to high spiking rates and long doppler response.

Summary

Proposed algorithm

Produces informative data representations suited to analyze and compare multimodal data.

Stepwise structure makes it flexible and allows integration of prior information.

Electro-Metabolic Coupling

Hemodynamic response was learned that well explains the data. Synchronization of spiking rates with periodic activity of 0.9-1 Hz in laser doppler becomes apparent in timefrequency map.

  Maureen Clerc 1 , Alexandre Gramfort 2 , Sandrine Saillet 3 , Christian Bénar 4 , Théodore Papadopoulo 1 1 Project-Team Athena, Inria Sophia Antipolis, France 2 Telecom ParisTech -CEA Neurospin, Paris, France 3 Institut de Neurosciences des Systèmes, Marseille, France sebastian.hitziger@inria.fr https://team.inria.fr/athena/ suited to learn meaningful representations of the multimodal datasets. As a result, we are able to learn a hemodynamic response and discover spike synchronization with hemodynamic activity. This work was supported by the ANR grants CoAdapt (09-EMER-002-01) and MultiModel (2010-BLAN-0309-04) and the doctoral grant of the region Provence-Alpes-Côte d'Azur.

  : LFP data is segmented into time windows centered around spikes [1] Drawbacks: (a) spikes need to be detected previously, (b) overlaps are not properly treated.

  Outline for learning waveforms (here spikes) and their latencies and amplitudes . Algorithm starts learning a single waveform