Frédéric Sur

Benoît Blaysat

Michel Grédiac

Determining displacement and strain maps immune from aliasing effect with the grid method

Keywords: aliasing, spectrum folding, grid method, displacement and strain maps

Spatial aliasing may affect methods based on grid image processing to retrieve displacement and strain maps in experimental mechanics. Such methods aim at estimating these maps on the surface of a specimen subjected to a loading test. Aliasing, which is often not noticeable to the naked eye in the grid images, may give spurious fringes in the strain maps. This paper presents an analysis of aliasing in this context and provides the reader with simple guidelines to minimize the effect of aliasing on strain maps extracted from grid images.

Introduction

Non-interferometric full-field measurement techniques are spreading more and more in the experimental mechanics community. Digital image correlation (DIC) [START_REF] Sutton | Image Correlation for Shape, Motion and Deformation Measurements[END_REF] is certainly the most popular one thanks to its versatility. Indeed, this technique does not only work with randomly marked patterns, which are the most common ones because they are quite easy to deposit on the external surface of specimens, but also with regular patterns such as dots or grids. In this latter case, however, another possibility consists in processing the corresponding grid images using spectral methods to obtain in-plane displacement and strain distributions. These spectral methods present several advantages: for instance, metrological performances can be quantified with closed-form expressions [START_REF] Grédiac | Effect of sensor noise on the resolution and spatial resolution of displacement and strain maps estimated with the grid method[END_REF], no optimization algorithm is needed, and the retrieved strain maps are not the result of an interpolation between remote points where the measurement is performed. Nevertheless, a phenomenon often occurs when acquiring images of grids with a digital camera: the presence of parasitic fringes, that is, quasi-periodic patterns superposed to the expected grid image. Although these fringes may be hardly noticeable in the acquired images, they cause spurious patterns in the estimated displacement maps, and mainly in the strain maps. This phenomenon is due to the fact that regular patterns are shot by a camera sensor, which is made of an array of pixels. The signal is thus sampled, and aliasing may occur if the Nyquist condition is not satisfied in the Petersen-Middleton sampling theorem [START_REF] Petersen | Sampling and reconstruction of wave-number-limited functions in N -dimensional Euclidean spaces[END_REF] (a generalization to images of the 1D Nyquist-Shannon sampling theorem [START_REF] Mallat | A wavelet tour of signal processing[END_REF]).

It is worth mentioning that these parasitic fringes are of low amplitude and even often imperceptible to the naked eye in grid images. The impact of these parasitic fringes on displacement and strain maps is, however, potentially significant, these maps being retrieved from actual changes between images which are generally of low amplitude. This is precisely the case in many situations of deformed materials or structures studied with full-field measurement techniques. This phenomenon also occurs when local strain increments are to be calculated, as in plasticity. Indeed, strain increments between consecutive images may be small even though the strain itself is significant. As a consequence, strain increment maps are also potentially affected by spurious fringes. The effect of aliasing on strain and strain increment maps is illustrated in Figure 1.

A first solution to avoid these fringes is to defocus the lens in order to blur the acquired images. This approach limits aliasing by smoothing out the high-frequency components of the imaged grid, thus enforcing the Nyquist sampling condition. Nevertheless, defocusing reduces the sharpness and the contrast of the acquired image, which causes the signal-to-noise ratio in final strain maps to decrease according to the discussion in [START_REF] Grédiac | Effect of sensor noise on the resolution and spatial resolution of displacement and strain maps estimated with the grid method[END_REF][START_REF] Sur | Towards deconvolution to enhance the grid method for in-plane strain measurement[END_REF]. It has been shown, in particular, that the standard deviation of the noise in final displacement and strain maps increases when the contrast decreases. While the anti-aliasing filter (an optical low-pass filter) which equips many digital cameras is supposed to eliminate aliasing, it is often not present in high-end cameras, and it does not completely eliminate aliasing, since a trade-off between sharpness and Nyquist condition must be reached. Another potential solution, effective in some cases, consists in eliminating aliasing through a post-processing step, for instance by applying a suitable notch filter in the frequency domain, as explained in [START_REF] Grédiac | Removing quasi-periodic noise in strain maps by filtering in the Fourier domain[END_REF] in the case of strain maps. Designing such notch filters for eliminating so-called pseudo-periodic noise in images is the aim of several papers in the image processing or microscopy community [START_REF] Aizenberg | Frequency domain median-like filter for periodic and quasiperiodic noise removal[END_REF][START_REF] Aizenberg | A windowed Gaussian notch filter for quasi-periodic noise removal[END_REF][START_REF] Hudhud | Digital removal of power frequency artifacts using a Fourier space median filter[END_REF][START_REF] Ionita | Automatic moiré pattern removal in microscopic images[END_REF][START_REF] Sur | An a-contrario approach to quasi-periodic noise removal[END_REF][START_REF] Wang | 2D FFT for periodic noise removal on strain images[END_REF][START_REF] Wei | A median-Gaussian filtering framework for moiré pattern noise removal from X-ray microscopy image[END_REF], but it seems that the procedures enabling the automatic design of such filters are less efficient in the case of displacement and strain maps, mainly because these maps feature statistical properties which are different from those of natural images [START_REF] Sur | Automated removal of quasiperiodic noise using frequency domain statistics[END_REF]. In addition, some information is irremediably lost in such a procedure. Note finally that parasitic fringes in final 2 strain maps obtained from grid images may be due to other reasons than aliasing, for instance:

• regular grid defects due to printing. Processing the images using the method proposed in [START_REF] Badulescu | Investigation of the grid method for accurate in-plane strain measurement[END_REF][START_REF] Badulescu | A procedure for accurate onedimensional strain measurement using the grid method[END_REF] enables one to get rid of this effect; • unsuitable window and/or settings when using the windowed Fourier transform to process grid images, as explained in [START_REF] Sur | Influence of the analysis window on the metrological performance of the grid method[END_REF].

These last two cases are not addressed here: we only consider the case of fringes impairing grid images caused by aliasing. These fringes propagate up to the final displacement or strain maps when processing the grid images.

In this context, the aim of this paper is to propose and analyze a simple procedure to avoid this undesirable phenomenon: rotating the grid with respect to the camera sensor by a certain angle, whose determination will be detailed. The benefit, compared to the aforementioned solutions, is to keep the highest contrast in the grid images (and thus the lowest noise level in the displacement and strain maps), and to avoid the awkward design of notch filters in the Fourier domain.

The paper is organized as follows: Section 2 briefly recalls the basics of the spectral methods used classically to process grid images. The continuous (in the sense that the space variable is continuous) and discrete (the space variable is sampled) spectra of grid images are then discussed in Section 3, with a special emphasis put on an important consequence of sampling, namely aliasing. We propose in Section 4 a numerical example, which shows that rotating the grid by a certain angle allows us to avoid aliasing, and thus the appearance of spurious fringes in strain maps. The effectiveness of the method is illustrated by real strain maps obtained with two grids deposited on a specimen subjected to a tensile test (one is rotated, the other is not). Section 5 discusses the optimum rotation angle between the imaged grid and the pixel grid. Some demonstrations and reminders are gathered in separate annexes.

2 Grid methods in experimental mechanics 2.1 Retrieving displacement and strain maps from grid images using spectral methods

In experimental mechanics, grids have long been used to measure in-plane displacement and strain components by processing deformed and nondeformed grid images. A detailed review being given in [START_REF] Grédiac | The grid method for in-plane displacement and strain measurement: a review and analysis[END_REF], we only recall here the main available techniques, and we describe in more detail the one used in the present work. Grids being deformed, they can no more be considered as a perfect periodic signal but instead as a phase-modulated signal. The first step is thus to estimate phase modulations. Basically, two spectral techniques have been proposed to deduce displacement and strain maps from grid images, namely geometric phase analysis (GPA), and methods based on windowed Fourier analysis at the nominal frequency of the grid. With GPA, the phase is estimated by taking the Fourier transform of images of a regular pattern, filtering the corresponding spectra with a notch filter centered around the spike of the principal frequency of the pattern, performing the inverse Fourier transform to retrieve the corresponding phase distributions and unwrapping the obtained results. This technique was inspired by procedures employed when processing fringes obtained in interferometry, as explained in the seminal papers [START_REF] Takeda | Spatial-carrier fringe-pattern analysis and its applications to precision interferometry and profilometry: An overview[END_REF][START_REF] Takeda | Fourier transform profilometry for the automatic measurement of 3-d object shapes[END_REF] for instance. It was then adapted to lattice images observed at the microscopic scale in [START_REF] Hÿtch | Quantitative measurement of displacement and strain fields from HREM micrographs[END_REF], and finally extended at the macroscopic scale on grids directly deposited on specimens in [START_REF] Dai | The geometric phase analysis method based on the local high resolution discrete Fourier transform for deformation measurement[END_REF][START_REF] Dai | Geometric phase analysis based on the windowed Fourier transform for the deformation field measurement[END_REF]. The windowed version of the Fourier transform can also be employed to process grid images, as proposed in [START_REF] Surrel | Photomechanics[END_REF]. The region surrounding the spikes of the principal frequencies is considered in [START_REF] Dai | The geometric phase analysis method based on the local high resolution discrete Fourier transform for deformation measurement[END_REF][START_REF] Dai | Geometric phase analysis based on the windowed Fourier transform for the deformation field measurement[END_REF], but only the nominal frequency of the grid is taken into account in [START_REF] Surrel | Photomechanics[END_REF], which dramatically accelerates the calculations.

Recent applications mainly deal with the technique presented in [START_REF] Surrel | Photomechanics[END_REF]. In addition, it has been shown that closed-form expressions for various metrological parameters can be elaborated in this case [START_REF] Grédiac | Effect of sensor noise on the resolution and spatial resolution of displacement and strain maps estimated with the grid method[END_REF][START_REF] Grédiac | The grid method for in-plane displacement and strain measurement: a review and analysis[END_REF][START_REF] Sur | Towards deconvolution to enhance the grid method for in-plane strain measurement[END_REF], which is important to eventually give the measurement uncertainty of the obtained results, as defined in [START_REF]International vocabulary of metrology. basic and general concepts and associated terms (VIM)[END_REF]. For all these reasons, we employ the windowed Fourier transform (WFT) with the nominal frequency of the grid to process the grid images discussed in this paper. The analysis of aliasing proposed here is, however, still valid when the grid images are processed with other spectral methods, for instance GPA [START_REF] Dai | The geometric phase analysis method based on the local high resolution discrete Fourier transform for deformation measurement[END_REF][START_REF] Dai | Geometric phase analysis based on the windowed Fourier transform for the deformation field measurement[END_REF]. Indeed and as we shall see, the procedure presented here permits to mitigate the effect of aliasing on the useful spectrum component. It is independent of the spectral method used to process the grids and retrieve displacement and strain maps. In contrast, methods that process grid images globally (for instance digital image correlation [START_REF] Sutton | Image Correlation for Shape, Motion and Deformation Measurements[END_REF], or the methods described in [START_REF] Allais | Experimental characterization of the local strain in a heterogeneous elastoplastic material[END_REF] or [START_REF] Du | Grid-based matching for full-field large-area deformation measurement[END_REF]), instead of focusing on a limited part of the spectrum, are potentially still affected by aliasing.

Notations and reminder

In this article, the Fourier transform of any integrable 2D function f (x, y) is denoted by F (f ) (u, v) and is defined as:

F (f ) (u, v) = R 2 f (x, y)e -2iπ(xu+yv) dx dy (1)
Similarly, the windowed Fourier transform is defined at any point (ξ,η) by

F w (f )(ξ, η, u, v) = R 2 f (x, y)w(x -ξ, y -η)e -2iπ(xu+yv) dx dy ( 2
)
where w is the analysis window. This is an isotropic and positive function integrating to 1. The Gaussian window defined by

w(x, y) = 1 2πσ 2 e -x 2 +y 2 2σ 2 (3)
is employed here, where σ denotes its standard deviation. We keep the same notations for the Fourier transform of tempered distributions, which allows us to deal with Dirac's delta or constant functions.

Both the 2D Fourier transform and its windowed version are covariant with any rotation of the image domain. More precisely, for any θ ∈ R, if f R (x, y) = f (x cos(θ)+y sin(θ), -x sin(θ)+y cos(θ)), we have

F (f R ) (u, v) = F (f ) (u cos(θ) + v sin(θ), -u sin(θ) + v cos(θ)) (4)
and

F w (f R )(ξ, η, u, v) = F w (f )(ξ cos(θ) + η sin(θ), -ξ sin(θ) + η cos(θ), u cos(θ) + v sin(θ), -u sin(θ) + v cos(θ)) (5)
The proof of this property is reminded in Appendix A.

The modulus of any complex number z is denoted by |z|, and if z = 0, its argument is denoted by arg(z).

The Dirac distribution centered at (a, b) ∈ R 2 is denoted by δ a,b . The convolution of two integrable functions or of two tempered distributions is denoted by f * g. The relation F (f * g) = F (f ) F (g) holds.

Estimating phases through windowed Fourier transform

We consider that we have a set of grid images shot during a mechanical test. The bidimensional signal corresponding to each image is denoted by s(x, y).

Let us assume for a while that grids, that is, orthogonal sets of parallel lines, are aligned with the 2D coordinate axes. Following [START_REF] Badulescu | Investigation of the grid method for accurate in-plane strain measurement[END_REF], a continuous grid can be modeled as follows:

s(x, y) = A(x, y) 2 2 + γ • (2πf x + φ 1 (x, y)) + γ • (2πf y + φ 2 (x, y)) (6)
where

• A(x, y) > 0 is the global field illumination, which may slightly vary across the field of view;

• γ ∈ [0, 1] is the contrast of the oscillatory pattern;

• the line profile is a 2π-periodic real function with a peak amplitude equal to 1 and average value 0;

• f is the nominal frequency of the carrier, the grid pitch being p = 1/f ;

• φ 1 (x, y) and φ 2 (x, y) are the carrier phase modulations due to specimen surface displacements along the x and y-axes respectively. In practice, their spatial derivatives are very small with respect to 2πf .

In real grid images, the amplitude of the signal is not exactly twice that of each of the two line networks considered separately. Although this does not affect the conclusion of the calculation performed below, we discuss this point in Section 4.

As explained above, we consider here the nominal frequency of the grid to process the grid images with the WFT recalled in [START_REF] Aizenberg | A windowed Gaussian notch filter for quasi-periodic noise removal[END_REF]. The value of the WFT is thus calculated twice at any pixel (ξ, η) of the image, once for (f, 0) and once for (0, f ), where f is equal to the nominal frequency of the grid (referred to as the principal harmonic in the remainder of the paper). This frequency is assumed to be known a priori. The standard deviation σ in ( 3) is equal to 1/f as requested in [START_REF] Sur | Towards deconvolution to enhance the grid method for in-plane strain measurement[END_REF]. We obtain two complex numbers defined at each pixel (ξ, η):

F w (s)(ξ, η, f, 0) and F w (s)(ξ, η, 0, f ).
It is generally assumed that the arguments of these complex numbers are merely the two sought phases φ 1 and φ 2 . More rigorously, these quantities are, up to first order approximations, the sought phases convolved by the window w (up to an additive constant, see [START_REF] Sur | Towards deconvolution to enhance the grid method for in-plane strain measurement[END_REF]). The convolution slightly blurs the actual distributions. Note that these phases have to be unwrapped. In the case of a rotated grid, the phases must be brought back along the horizontal and vertical axis through the opposite rotation:

Φ 1 (x, y) = φ 1 (x, y) cos(θ) -φ 2 (x, y) sin(θ) Φ 2 (x, y) = φ 1 (x, y) sin(θ) + φ 2 (x, y) cos(θ) (7)
Both phases Φ 1 (Φ 2 , respectively), taken after and before deformation, are subtracted and finally multiplied by -p 2π to deduce the in-plane displacement at any pixel along the horizontal axis (the vertical axis, respectively). Note that the movement of the physical points shall be taken into account to get rid of the grid defects, as proposed in [START_REF] Badulescu | Investigation of the grid method for accurate in-plane strain measurement[END_REF]. The interested reader is referred to the review paper [START_REF] Grédiac | The grid method for in-plane displacement and strain measurement: a review and analysis[END_REF] for more details. The linearized strain components are finally obtained by calculating the symmetric part of the displacement gradient.

The remainder of the paper examines how seeking the displacement and strain distributions may be disturbed by aliasing, which potentially affects the grid images. Without loss of generality, the discussion of aliasing is led with the Fourier transform. It is indeed still valid when the signal is localized through windowing.

Spectrum of a grid

This section discusses the spectrum of a grid, within the model of [START_REF] Dai | The geometric phase analysis method based on the local high resolution discrete Fourier transform for deformation measurement[END_REF]. Let (x) = k∈Z d k e ikx be the Fourier expansion of the 2π-periodic function. Since is a 0-mean real function, d 0 = 0 and for any k, d k is the complex conjugate of d -k . A pure sine is theoretically the ideal line profile [START_REF] Grédiac | Effect of sensor noise on the resolution and spatial resolution of displacement and strain maps estimated with the grid method[END_REF] (in which case for any k = ±1, d k = 0), but it is not possible to print such grids in practice.

Continuous spectrum

The Fourier transform of s writes (see Appendix B):

F (s) (u, v) = F (A) (u, v) + γ 2 k∈Z d k F (A) * F e ikφ1 (u -kf, v) + γ 2 k∈Z d k F (A) * F e ikφ2 (u, v -kf ) (8)
If A can be considered to be constant, then F (A) = A δ 0,0 and (8) simplifies into the following equation:

F (s) (u, v) = A δ 0,0 + γA 2 k∈Z d k F e ikφ1 (u -kf, v) + γA 2 k∈Z d k F e ikφ2 (u, v -kf ) (9)
Assuming that both partial derivatives of φ l (l = 1 or 2) are bounded, Proposition C.1 (in Appendix C) proves that, for any k, the support of F e ikφ l is included in

- |k| p M x l , |k| p M x l × - |k| p M y l , |k| p M y l ( 10
)
where M x l = (p/2π) max(|∂φ l /∂x|) and M y l = (p/2π) max(|∂φ l /∂y|). Note that with these notations, M x l and M y l can be considered as bounds for the strain components, under the hypothesis that the grid on the specimen is perfectly aligned with the pixel grid in the reference state, the phase modulation being produced only by the deformation of the specimen.

This means that a grid image spectrum is made of a central spike caused by the F (A) term and separated spikes centered at multiples of the frequency f along the horizontal and vertical axes of the spectrum. These spikes are caused by the F e ikφ1 and F e ikφ2 terms, and by the harmonics of the line profile . These spikes are therefore mostly included in

- |k| p M x 1 + kf, |k| p M x 1 + kf × - |k| p M y 1 , |k| p M y 1 ( 11
)
along the horizontal direction, and

- |k| p M x 2 , |k| p M x 2 × - |k| p M y 2 + kf, |k| p M y 2 + kf (12)
along the vertical direction. We write "mostly" because of the blur caused by the convolution of the Fourier transform F (A) of the illumination field in [START_REF] Du | Grid-based matching for full-field large-area deformation measurement[END_REF], in addition to the approximate calculation in the stationary phase method. Moreover, this property also means that the support of the spikes is all the narrower as the upper bounds of the partial derivatives of the phase are small, and as k ≥ 1 is small. Figure 2 illustrates a typical spectrum of a grid. Since a spectrum of a rotated grid is simply the rotated spectrum of a grid aligned with the pixels (see Section 2.2 above), the preceding discussion adapts easily to rotated grids.

u v f 3 2 1 2 1 3 f 1 2 3 1 2 3 k 2|k| p M y 1 2|k| p M x 1 k 2|k| p M y 2 2|k| p M x 2 Figure 2:
An illustration of the continuous spectrum of a grid, supposed here to generate only three pairs of harmonics along each axis.

Discrete spectrum

Digital sensors give sampled information. Assuming square pixels, which is a sound assumption in most real cameras, Nyquist-Shannon sampling theorem states that the continuous signal is well sampled if and only if its spectrum is supported by the box

[-f s x /2, f s x /2] × [-f s y /2, f s y /2],
where f s x , f s y are the sampling frequencies along both discretization directions. In what follows, f s x = f s y = f s = 1 pixel -1 . Note that the case of non-square pixels is still tackled by Peterson-Middleton theorem, the localization condition becoming that the continuous spectrum must be included in the so-called reciprocal cell.

Let us assume that the camera gives a discrete image s of size M × N (M and N are assumed to be even), which samples the continuous s discussed in the preceding section. The question is: how the spectrum of the discrete image (obtained with the discrete Fourier transform of s) is related to the spectrum of the continuous image given by ( 8)?

Let c m,n be the discrete Fourier coefficient (with The support of the continuous image s such that s(k, l) = s(k, l) is supposed to be included in

-M/2 ≤ m ≤ M/2 -1 and -N/2 ≤ n ≤ N/2 -1) of s, that is c m,n = M/2-1 k=-M/2 N/2-1 l=-N/2 s(k, l)e -2iπ(km/M +ln/N ) (13)
[-M/2, M/2 -1] × [-N/2, N/2 -1] (that is, s = 0 outside of this box). From Appendix D, for any m, n such that -M/2 ≤ m ≤ M/2 -1 and -N/2 ≤ n ≤ N/2 -1, the Fourier coefficient c m,n satisfies: c m,n = (k,l)∈Z 2 F (s) m M + k, n N + l (14) = F (s) m M , n N + (k,l)∈Z 2 \{0,0} F (s) m M + k, n N + l (15)
This means that the discrete spectrum (i.e., the set of the M ×N coefficients c m,n ) corresponds to the discretized continuous spectrum F (s), to which are added copies of the discretized F (s) translated by an integer value (the so-called aliases, or folded components). Aliasing thus manifests itself by impairing the spectrum of the sampled image by the components of the continuous spectrum lying outside [-1/2, 1/2] × [-1/2, 1/2], thus breaking the Nyquist condition.

Impact on phase estimation from grid images

Aliasing causes spurious fringes in the phase maps, which propagate to the displacement and mainly strain maps through the process described in Section 2.3. Since these fringes are not identical after and before deformation, they do not vanish when subtracting the phases between the current and the reference states. In this section, we focus on the impact of aliasing on the phase maps. The method described in [START_REF] Badulescu | Investigation of the grid method for accurate in-plane strain measurement[END_REF][START_REF] Grédiac | The grid method for in-plane displacement and strain measurement: a review and analysis[END_REF] (after [START_REF] Surrel | Photomechanics[END_REF]) for finding the phase is based on the estimation of F (s) (f, 0) and F (s) (0, f ), and the method described in [START_REF] Dai | The geometric phase analysis method based on the local high resolution discrete Fourier transform for deformation measurement[END_REF][START_REF] Dai | Geometric phase analysis based on the windowed Fourier transform for the deformation field measurement[END_REF] (after [START_REF] Takeda | Fourier transform profilometry for the automatic measurement of 3-d object shapes[END_REF][START_REF] Takeda | Spatial-carrier fringe-pattern analysis and its applications to precision interferometry and profilometry: An overview[END_REF]) is based on the estimation of F (s) in neighborhoods of (f, 0) and (0, f ) (in order to obtain estimations of F e ikφ1 and F e ikφ2 in (8)). As mentioned above, these methods are actually based on windowed Fourier analysis, but the discussion of aliasing in the previous section is still valid locally.

Nevertheless, the input datum (the imaged grid) is acquired by a digital sensor and is discrete. The continuous spectrum F (s) is not available, and the calculations are led with the discrete Fourier coefficients c m,n . As a consequence of the discussion in Section 3.2, aliasing impairs the retrieved phase maps as soon as the Fourier coefficients of interest are impaired by the terms F (s) m M + k, n N + l , with (k, l) = (0, 0). Even if the "spike" around (f, 0) or (0, f ) is not aliased (or, in other words, if it is fully included in [-1/2, 1/2] × [-1/2, 1/2]), it may be covered by aliased harmonic spikes outside of this box in the original spectrum of s. In this case, aliasing appears in the retrieved phases and subsequently in the estimated displacement and strain maps.

In the case of grid images aligned with the pixel grid of the sensor, this means that the k-th harmonics of the line profile give aliasing as soon as k ≥ f s /(2f ). For instance, with a typical value of f = 1/5 pixel -1 , meaning that the pattern pitch is encoded on five pixels, harmonics such that k ≥ 3 are aliased and thus potentially cover the spikes of interest. With a 45-degree rotation, harmonics such that k ≥ √ 2f s /(2f ) (thus k ≥ 4) are aliased for the same f /f s value. Note that decreasing f (i.e., increasing the grid pitch p = 1/f , thus the number of pixels imaged over a period) also decreases the number of aliased harmonics. This solution is often not possible in practice because of the trade-off to be satisfied between the pixel density of the sensor and the field of view of the lens. Figure 3 shows, however, that aliasing may be got rid of by rotating the Figure 3: An illustration of the discrete spectrum of a grid image, with an aligned-to-pixel grid (on the left) and a rotated grid (on the right). With the notations of ( 15), the translated copies of F (s) for (k, l) ∈ {(0, 1), (0, -1), (1, 0), (-1, 0)} are shown in a light shade, superposed to the non-translated F (s) in the central, gray square. The discrete Fourier coefficients c m,n are the sum of all the components superposed in the gray square. We can see in the left sketch that the principal harmonic is covered by an aliased spike. On the contrary, the right sketch shows that turning the grid prevents the aliased spikes from covering the principal harmonic. grid instead of increasing the grid pitch. According to the property reminded in Section 2.2 above, a rotation of the grid image is indeed equivalent to a rotation of the spectrum F (s). It can be noted that the rotation angle θ and the indices (m * , n * ) of the Fourier coefficient corresponding to the principal harmonic satisfy tan(θ) = (m * N )/(n * M ). The resulting procedure is investigated in the following section.

Illustrative examples

We propose herein two illustrative examples: a numerical experiment and a real-case study. We demonstrate that rotating the grid image with respect to the pixel grid significantly reduces the effect of aliasing on the estimation of the displacement and strain maps.

Numerical experiment

This section provides the reader with a numerical illustration of aliasing in the context of displacement and strain map estimation. Grid images are simulated according to the model of ( 6), with a constant illumination A = 2 11 , a contrast γ = 0.9, a line profile (x) = sin 7 (x) (chosen to generate realistic sharp lines and several harmonics), and p = 1/f the grid pitch. Such grid images mimic the output of a standard 12-bit camera. Concerning the phase modulations, we choose φ 1 = 0, and φ 2 such that ∂φ 2 /∂y is a sine wave along the x-direction, with a linearly varying period along y. The amplitude of this sine wave is constant, equal to 5 • 10 -3 pixel -1 . While (6) gives a continuous aligned-to-pixel grid, a clockwise rotation of 10 degrees around the origin of the pixel coordinates gives the rotated grid. Both grids are then sampled to produce an N × N image. Here, the size of an image is chosen to linearly increase with the pattern pitch p in order to simulate an increased image resolution of a given field of view. We set N = 200(p -4) so that the minimum pattern pitch (p = 5) gives a 200 × 200 image. This procedure permits to generate images of fixed grids, their resolution increasing with the pitch p. Figure 4 illustrates the generated grids and the corresponding maps of ∂φ 2 /∂y, for p = 7.75 pixels. Note that the phase derivatives are brought back in the same coordinate axis through the procedure described in Section 2.3 above. The function being odd, its Fourier coefficients d k are such that d k = 0 for any odd integer k, and the calculation shows that d k = 0 only for k = ±1, ±3, ±5, ±7. The spectrum of the alignedto-pixel grid image, as shown in Figure 5, complies with the theory of Section 3. Spikes can be seen on the horizontal axis in spite that φ x = 0 because of sampling, which results in a slight phase Frequencies are in cycles per image. The principal harmonic is circled in green. At the bottom: retrieved phase derivative from the aligned-to-pixel grid (on the left) and from the rotated grid (on the right). These figures should be compared to the ground truth presented in Figure 4. Aliasing in the aligned-to-pixel grid dramatically transfers to the phase derivative maps as parasitic spurious fringes.

modulation. Along the vertical axis, four "spikes" can be seen at a distance of |k|/p = 0.129 • k (for k = ±1, ±3) cycles per image from the frequency origin, and four spikes correspond to the folding of the spikes expected for k = ±5, ±7, for which |k/p| > 0.5. Folded spikes can be found respectively at a distance of 1-0.129•5 = 0.355 and 1-0.129•7 = 0.097 of the frequency origin. Since the folded spike corresponding to k = 7 overlaps the principal harmonic at 1/p = 0.129, aliasing is expected to impair the retrieved phase maps. Concerning the rotated grid image, spectrum folding can still be noted. Nevertheless, a visual inspection shows that, in this case, folded spikes do not overlap the principal harmonic. Aliasing is thus expected to be far less noticeable in the retrieved phase maps. Figure 5 also shows the estimated phase derivative maps ∂φ y /∂y. While aliasing indeed produces strong fringes as expected when the phase is estimated from the aligned-to-pixel grid, no such spurious phenomenon can be seen when the estimation is performed on the rotated grid.

A companion video (file video_aliasing_synth.avi) shows the evolution of the retrieved phase derivative with such synthetic grid images. Both the spectrum and the retrieved phase derivative are shown for grid pitch values increasing from 5 to 10 pixels, for the aligned-to-pixel and the rotated grids. Aliased spikes can be seen. They propagate to the retrieved phase derivative map as soon as they cover the principal harmonic, marked with a green circle. In this case, dramatic interferences can be noted. Rotating the grid adequately (here, with a 15-degree angle) permits to shift away the aliased spikes outside the principal harmonic in the spectrum. In this case, no spurious fringe can be seen in the retrieved phase derivative.

It is interesting to note that, when an aliased spikes (corresponding to the k-th harmonic) is perfectly superposed to the principal harmonic, the Fourier coefficient c f M,0 corresponding to the normalized frequency (f, 0) is given by:

γA 2 d 1 F e iφ1 (0, 0) + d k F e ikφ1 (0, 0) (16)
(we use the simplified expression [START_REF] Grédiac | Effect of sensor noise on the resolution and spatial resolution of displacement and strain maps estimated with the grid method[END_REF], knowing that the non-aliased harmonics vanish at (f, 0) in this expression).

As explained in the introduction, the estimation of the phase φ 2 is performed by taking the argument of the coefficient given in [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. Moreover, F e ikφ2 (0, 0) is simply the weighted average value of e ikφ2 through the analysis window centered at (ξ, η), the argument being approximately equal to w * φ 2 (ξ, η), see [START_REF] Sur | Towards deconvolution to enhance the grid method for in-plane strain measurement[END_REF].

The argument of c f M,0 can eventually be written as:

arg(d 1 ) + w * φ 2 (ξ, η) + arg 1 + d 1 F e ikφ2 (0, 0) d k F (e iφ2 ) (0, 0) (17)
this equation holding modulo 2π. Now, the ratio F e ikφ2 (0, 0)/F e iφ2 (0, 0) is a complex number of modulus approximately equal to 1, and of argument (1-k)w * φ 2 (ξ, η), which is small in the case of small displacements. Since, with the synthetic grids, d 1 /d k is always a real number ( (x) = sin 7 (x) being an odd function, its Fourier coefficients are all imaginary numbers), the last argument in [START_REF] Petersen | Sampling and reconstruction of wave-number-limited functions in N -dimensional Euclidean spaces[END_REF] can be considered to be 0. This discussion shows that, when an aliased spike is perfectly superposed to the principal harmonic, the resulting interference is negligible. This can be seen in the companion video when the grid pitch is set exactly to 6 pixels. Of course, this is only valid in the present case of synthetic grids: real grid profiles are likely to show asymmetries which prevent from being modeled as a strictly odd function.

Real-case study

A tensile test on a specimen equipped with two grids (one is aligned with the pixel grid, the other is rotated) was carried out to illustrate the benefit of rotated grids with respect to aliasing. An open-hole specimen in polymethylmethacrylate (PMMA) was considered here. The dimensions are 220 × 50 × 4 mm 3 . The diameter of the hole is 4 mm. This hole was drilled to cause a heterogeneous strain distribution to locally appear in the specimen. The two grids were transferred on both sides of the hole following the procedure proposed in [START_REF] Piro | Producing and transferring low-spatial-frequency grids for measuring displacement fields with moiré and grid methods[END_REF] (see also [START_REF] Grédiac | The grid method for in-plane displacement and strain measurement: a review and analysis[END_REF]), as illustrated in Figure 6. The nominal pitch of these grids was 0.2 mm. On the right-hand side of the hole, the grid was rotated by 19.5 degrees while the lines of the grid were parallel to the borders of the specimen on the left-hand side. The specimen was then subjected to a tensile test with a tensile machine. It can be checked that the material did not plastically deform with the corresponding stress level. A Sensicam QE 12-bit CCD camera shot the images of both grids while testing. The distance between specimen and camera was such that p = 5 pixels were used to encode one grid period. These grid images were then processed separately to retrieve the phase distribution in each case, and then to deduce the ε xx strain distribution by using the procedure described in Section 2.3 above (x represents the horizontal direction in Figure 6). One can notice that, with a small number of pixels per period (typically 5 pixels per period), a rotation of the imaged grid makes the contrast to slightly decrease, as peaks of the line profile are not aligned any more with the pixel grid but instead cross several pixels. As a consequence, the noise floor increases in displacement and strain maps, as explained in [START_REF] Grédiac | Effect of sensor noise on the resolution and spatial resolution of displacement and strain maps estimated with the grid method[END_REF]. A larger number of pixels per period should circumvent this detrimental effect, as the sharpness of the image is likely to increase in this case.

As mentioned in Section 3, pixel intensity at the crossing of the lines is not exactly twice as high as the intensity of the line. This is not crucial since this model is accurate enough to build predictive formula which have been experimentally assessed in [START_REF] Grédiac | Effect of sensor noise on the resolution and spatial resolution of displacement and strain maps estimated with the grid method[END_REF]. This explains, however, that the spectrum of a real grid image is not exactly as the one of the numerical simulation shown in Figure 5. A more accurate model of the continuous grid would be to enhance [START_REF] Dai | The geometric phase analysis method based on the local high resolution discrete Fourier transform for deformation measurement[END_REF] with an additive term B * k,l δ kp,lp , where B is the corrective intensity at the crossing of any two grid lines, approximately lying at pixels (kp, lp). The Fourier transform of s becomes (8), to which the following term is added:

F   B * k,l δ kp,lp   = f 2 F (B) k,l δ kf,lf (18)
The Fourier transform of the 2D Dirac comb k,l δ kp,lp is indeed f 2 k,l δ kf,lf . Figure 6 shows the two grids on the sample. The spectra of the aligned-to-pixel grid and of the rotated grid are shown in Figure 7. The spikes observed around the frequencies (kf, lf ) are consistent with the previous discussion. The frequencies around (f, 0) and (0, f ) being covered by aliased spikes in the aligned-to-pixel grid, aliasing is expected to impair the retrieved strain map in this case. On the contrary, the aliased spikes do not seem to cover these frequencies in the spectrum of the rotated grid. The respective strain maps are shown in Figure 8. As expected from the aspect of the spectrum, aliasing affects the strain maps obtained from the aligned-to-pixel grids, and not those obtained from the rotated grids. The reader is invited to watch the second and third companion videos (available in supplementary files video_aliasing_real_epsxx.avi and video_aliasing_real_epsxy.avi). It can be seen that the strain distributions provided by the aligned-to-pixel grid are impaired by parasitic fringes, while those provided by the inclined one are not. The effect of aliasing is even stronger in the strain field increment. These fringes scroll while testing (thus as the load increases), which is due to the fact that the actual grid pitch slightly changes from one image to another because of the deformation of the specimen. In particular, the fringes disturb the actual strain distribution around the hole. The "flashing dots" correspond to the spatially correlated noise caused by sensor noise, which propagates from the images to the final strain maps whatever the orientation of the grid.

Towards an optimum rotation angle

The discussion of Section 3.2 tells us that if the line profile shows many harmonics of significant amplitude, then aliasing is mostly unavoidable. Nevertheless, if we assume that the line profile is made of only a few harmonics, then the corresponding spikes sparsely populate the discrete spectrum. It is possible to prevent aliased spikes to cover the principal harmonic by rotating the grid, as illustrated in Section 4 above. Intuitively speaking, we would like to set the rotation angle in such a way that the aliased spikes are as much as possible separated from the principal harmonic.

Let θ be the rotation angle. Thanks to the covariance property reminded in Section 2.2, the principal frequencies are thus (f cos(θ), -f sin(θ)) (to retrieve φ 1 in ( 6)) and (f sin(θ), f cos(θ)) (to retrieve φ 2 ). Without loss of generality, we focus on (f cos(θ), -f sin(θ)), and assume that 0 ≤ θ ≤ π/2. The aliased harmonics that are the most likely to interfere with this frequency come from the nearest copies of F (s) to the principal frequencies, that is the ones translated from (1, 0) and (0, -1) in [START_REF]International vocabulary of metrology. basic and general concepts and associated terms (VIM)[END_REF]. The discussion of Section 4.2 shows that spikes from these translated copies are expected to be centered at (kf cos(θ) + lf sin(θ) + 1, -kf sin(θ) + lf cos(θ)) [START_REF] Sur | An a-contrario approach to quasi-periodic noise removal[END_REF] and

(k f cos(θ) + l f sin(θ), -k f sin(θ) + l f cos(θ) -1) (20)
respectively, for any integer values k, l.

If k, l had real values, the center of an aliased spike would be superposed to the principal harmonic provided that:

f cos(θ) = kf cos(θ) + lf sin(θ) + 1 -f sin(θ) = -kf sin(θ) + lf cos(θ) (21)
which gives after solving the linear system:

k = 1 -cos(θ)/f l = -sin(θ)/f (22)
The nearest aliased spike to the principal frequency, among the copy of F (s) translated from (1, 0), is therefore such that:

k = 1 -[cos(θ)/f ] l = -[sin(θ)/f ] (23)
where [x] rounds x to the nearest integer.

The same calculation shows that the nearest aliased spike to the principal frequency, among the copy of F (s) translated from (0, -1), is such that:

k = 1 -[sin(θ)/f ] l = [cos(θ)/f ] (24)
In the first case, the squared distance between the principal harmonic and an aliased spike is

D 2 f (θ) = (kf cos(θ) + lf sin(θ) + 1 -f cos(θ)) 2 + (-kf sin(θ) + lf cos(θ) + f sin(θ)) 2 (25) = (k -1) 2 f 2 cos 2 (θ) + l 2 f 2 sin 2 (θ) + 1 + 2(k -1)lf 2 cos(θ) sin(θ) + 2(k -1)f cos(θ) +2lf sin(θ) + (k -1) 2 f 2 sin 2 (θ) + l 2 f 2 cos 2 (θ) -2(k -1)lf 2 cos(θ) sin(θ) (26) = (k -1) 2 f 2 + 2(k -1)f cos(θ) + 1 + l 2 f 2 + 2lf sin(θ) (27)
This expression simplifies into ((k -1)f + cos(θ)) 2 + (lf + sin(θ)) 2 .

In the second case (to retrieve φ 2 ), the same calculation gives ((k -1)f +sin(θ)) 2 +(l f -cos(θ)) 2 . With (k, l) as in [START_REF] Surrel | Photomechanics[END_REF] and (k , l ) as in [START_REF] Sutton | Image Correlation for Shape, Motion and Deformation Measurements[END_REF], these two squared distances are identical and we can conclude that the squared distance between the nearest aliased spike and the principal frequency is simply: where {x} = x -[x] denotes the distance between any real number x and the nearest integer to x. Figure 9 shows three plots of D f (θ) for three typical values of the grid pitch p = 1/f (5, 8, and 10 pixels). For example, for p = 1/f = 5 pixels, an angle of 45 o is optimum, but an angle of 30 o is still appropriate. It is interesting to note that, with p = 10 pixels, an angle of 45 o is not appropriate and an angle of 40 o should be preferred. As we can see, the value of the maxima of D f decreases as p increases. This is logical since the spikes are more densely spread over the spectrum as p increases, the distance between the principal frequency and the nearest aliased spike thus decreases. Plotting [START_REF] Wei | A median-Gaussian filtering framework for moiré pattern noise removal from X-ray microscopy image[END_REF] gives an indication on how to set the rotation angle, by choosing an angle maximizing D f . However, it may be practically more convenient to use a local maximum corresponding to a wide peak instead of the global maximum, which may correspond to a sharp peak. With a sharp peak, a small error on the angle is indeed likely to give aliasing. Moreover, the user shall not forget that only "first order" aliased spikes are considered in the present calculation, spikes coming from translation of the continuous spectrum larger than 2 in (15) are not considered. This simplification may be over-optimistic if the grid profile is not smooth enough, which would prevent the Fourier coefficients d k from quickly vanishing and would produce many harmonics with a significant amplitude. In addition, it should be mentioned that the spreading of the spectrum spikes is not taken into account in the preceding discussion: moving the center of the aliased spikes away from the principal frequency may not be sufficient because of the width of the aliased spikes which may be too large.

D 2 f (θ) = f 2 {cos(θ)/f } 2 + {sin(θ)/f } 2 (28)

Conclusion

In this paper, we analyzed spatial aliasing in the context of grid images acquired through a digital camera, and the resulting effect on the displacement and strain maps estimated from these images. The discussion is independent from the spectral technique used, even though the WFT set on the nominal frequency of the grid was employed here. We have proved in a theoretical model that rotating the grid permits to reduce aliasing. This approach was assessed both in a numerical study and in a real experiment. We have also proposed a way to optimally set the rotation angle so that the aliased spikes are as much as possible separated from the principal harmonic.

A Spectrum of a rotated image

Let θ ∈ R and R θ = cos(θ) sin(θ) -sin(θ) cos(θ) be the corresponding rotation matrix. For any 2D function f representing a continuous image, the Fourier transform of the rotated image writes:

F (f • R θ ) (u) = R 2 f (R θ x)e -2iπ<x,u> dx (29)
where u ∈ R 2 , < •, • > is the standard scalar product, and • denotes the composition between functions. A change of variable gives

F (f • R θ ) (u) = R 2 f (x)e -2iπ<R -1 θ x,u> dx (30)
Since any rotation is a unitary transform, we deduce < R -1 θ x, u >=< x, R θ u >, therefore (4) holds. Concerning the windowed Fourier transform, the same argument gives

F w (f • R θ )(ξ, u) = R 2 f (R θ x)w(x -ξ)e -2iπ<x,u> dx (31) = R 2 f (x)w(R -1 θ x -ξ)e -2iπ<x,R θ u> dx (32)
The analysis window being isotropic, w(R

-1 θ x -ξ) = w(R -1 θ (x -R θ ξ)) = w(x -R θ ξ), therefore (5) holds.

B Spectrum of a continuous grid

Let us recall that, with the definition of the Fourier transform given in (1), the following relations hold for tempered distributions:

F (1) = δ 0,0 (33)
F e -2πi(ax+by) = δ a,b (34)
F (f • g) = F (f ) * F (g) (35)
where δ a,b is the Dirac delta function concentrated at (a, b) ∈ R 2 . With this formula, the Fourier transform of (6) writes:

F (s) = F (A) + γ 2 k∈Z d k F A(x, y)e 2iπkf x+ikφ1(x,y) + F A(x, y)e 2iπkf y+ikφ2(x,y) (36) = F (A) + γ 2 F (A) * k∈Z d k F e 2πif kx * F e ikφ1 (37)
+ γ 2 F (A) * k∈Z d k F e 2πif ky * F e ikφ2
Consequently,

F (s) (ξ, η) = F (A) (ξ, η)+ γ 2 k∈Z d k F (A) * F e ikφ1 (ξ -kf, η)+ γ 2 k∈Z d k F (A) * F e ikφ2 (ξ, η-kf ) (38) Indeed, F e 2πif kx = δ -kf,0 .
In practice, the illumination field A varies slowly across the imaged grid, giving a concentrated spectrum around the mean value F (A) (0, 0). In the limit-case of a constant illumination A(x, y) = A, F (A) = A δ 0,0 , and the preceding formula simplifies into This proves ( 8) and ( 9).

C Stationary phase method

A simplified evaluation of F e iψ (with the notations of Section 3, ψ = kφ 1 or = kφ 2 ) is available through the stationary phase method, which we briefly recall here. The basic idea of this method is that sine functions with a quickly varying phase over a given region integrates to 0 on this region. Since F e iψ (u, v) = R 2 e iψ(x,y)-2iπ(xu+yv) dx dy (40) this principle justifies that the only significant contribution is caused by regions in R 2 where the phase ϕ(x, y) = ψ(x, y) -2π(xu + yv) does not oscillate (hence "stationary phase"), that is, where the partial derivatives of the phase ϕ are zero. Consequently, the Fourier transform of e iψ can be approximated by where the (x i , y i ) are such that ∂ψ/∂x(x i , y i ) = 2πu and ∂ψ/∂y(x i , y i ) = 2πv, and ε > 0 defines the half-width of the (small) region giving the contribution to the integral over R 2 . Under the same principle, F e iψ (u, v) = 0 if 2πu (resp. 2πv) is not within the range of the values taken by ∂ψ/∂x (resp. ∂ψ/∂y), which means that the phase derivatives cannot vanish.

Approximating ψ(x, y) -2π(xu + yv) as a function of (x, y) by its second-order Taylor expansion around (x i , y i ) gives, the gradient being null at these points: ψ(x, y) -2π(xu + yv) ψ(x i , y i ) -2π(x i u + y i v) + 1 2 (x -x i , y -y i ) H ψ (x i , y i ) (x -x i , y -y i ) T (42) where H ψ (x i , y i ) is the 2×2 Hessian matrix of ψ at (x i , y i ). Plugging (42) into (41) and substituting in the integral x -x i by x and y -y i by y gives: which is [START_REF] Ionita | Automatic moiré pattern removal in microscopic images[END_REF].

Figure 1 :

 1 Figure 1: Strain and strain increments measured on an open-hole specimen during a tensile test. Spurious fringes can be seen. They are caused by aliasing, which affects the grid image. This experiment is discussed in Section 4.2.

  Since the complex exponential function is 2iπ-periodic, c m,n can be considered as an (M, N )-periodic two-dimensional series along both indices M and N . With these notations, for any m, n such that -M/2 ≤ m ≤ M/2 -1 and -N/2 ≤ n ≤ N/2 -1, c m,n is the Fourier coefficient associated with a sinusoidal component of respective frequencies |m| and |n| cycles-per-image along horizontal and vertical axes. In this unity, frequencies are bounded by max(M, N )/2. Normalized frequencies |m|/M and |n|/N are in cycles-per-pixel, and are bounded by 1/2.

Figure 4 :

 4 Figure 4: Synthetic grid. Upper row: aligned-to-pixel grid. Bottom row: rotated grid. On the left: grid image. On the right: derivative of the imposed phase modulation.

3 Figure 5 :

 35 Figure5: Synthetic grid. On the top: spectrum of the aligned-to-pixel grid image (on the left) and of the rotated grid image (on the right). Frequencies are in cycles per image. The principal harmonic is circled in green. At the bottom: retrieved phase derivative from the aligned-to-pixel grid (on the left) and from the rotated grid (on the right). These figures should be compared to the ground truth presented in Figure4. Aliasing in the aligned-to-pixel grid dramatically transfers to the phase derivative maps as parasitic spurious fringes.

Figure 6 :

 6 Figure 6: Real-case study. The two grids on a sample.

Figure 7 :

 7 Figure 7: Real-case study. Spectra of the aligned-to-pixel grid (top-left) and of the rotated grid (top-right). Frequencies are in cycles per pixel. As can be seen in the bottom figure, aliased spikes are actually present in the spectrum of the aligned-to-pixel grid. These aliased spikes are shifted outside the principal harmonic when rotating the grid: in the annotated bottom spectrum, the aliased spikes are distributed in the straight lines on either sides of the middle line. In the aligned-to-pixel grid, the aliased spikes are almost superposed to the non-aliased ones.

5 # 10 - 3 Figure 8 :

 51038 Figure 8: Real-case study. Retrieved strain map ε xx from the aligned-to-pixel grid (on the left) and from the rotated grid (on the right). The large parasitic fringes that can be seen on the left are not present on the right. These fringes are caused by aliasing. Note the effect of the noise, inherent to any digital acquisition device.

Figure 9 :

 9 Figure 9: Plot of D f (θ) against θ ∈ [0, 90 o ] for three typical values of the grid pitch p = 1/f .

F

  (s) (ξ, η) = A δ 0,0 (ξ, η) + γA 2 k∈Z d k F e ikφ1 (ξ -kf, η) + γA 2 k∈Z d k F e ikφ2 (ξ, η -kf ) (39)

Fe

  e iψ (u, v) iψ(x,y)-2iπ(xu+yv) dx dy (41)

F 2 e i( 1 2 2 s 2 F 2 F

 22222 e iψ (u, v)(xi,yi)e iπψ(xi,yi)-2iπ(xiu+yiv) y) H ψ (xi,yi) (x,y) T ) dx dy (43)The domain of the preceding integrals can be extended to R 2 , if H ψ is non-degenerated, through the stationary phase principle (the contribution of the terms outside [-ε, ε]×[-ε, ε] being negligible because of the quickly varying phase), which gives:F e iψ (u, v) (xi,yi) e iπψ(xi,yi)-2iπ(xiu+yiv) R (x,y) H ψ (xi,yi) (x,y) T ) dx dy (44) = (xi,yi) 2πe iσ(H ψ (xi,yi))π/4 | det H ψ (x i , y i )| 1/2 e iπψ(xi,yi)-2iπ(xiu+yiv)(45)where det(A) denotes the determinant of any square matrix A and σ(A) is the sign of det A, since the integral over R 2 in (44) is a multidimensional Fresnel integral.If no stationary point of the phase in (41) exists, which means that the system ∂ψ ∂x (x, y) -2πu = 0 ∂ψ ∂y (x, y) -2πv = 0 (46) has no solution in (x, y), we obtain, for any u, v, F e iψ (u, v) 0 (47) In other words, we have proved the following proposition. Proposition C.1 If the absolute value of the derivatives ∂Ψ/∂x and ∂Ψ/∂y are bounded from above by B x and B y respectively, then for any (u, v) such that |u| > B x /(2π) and |v| > B y /(2π), F e iΨ (u, v) = 0.This means that the support of F e iΨ is included in-B x 2π , B x 2π × -B y 2π , B y 2π .D Link between discrete spectrum and Fourier transform of a continuous 2D signal Since s(k, l) = s(k, l) and s(x, y) = 0 as soon as (x, y) / ∈ [-M/2, M/2 -1] × [-N/2, N/2 -(k, l)e -2iπ(km/M +ln/N ) (49) We obtain the following expression with the Poisson summation formula (under mild assumptions on f , (k,l)∈Z 2 f (k, l) = (k,l)∈Z 2 F (f ) (k, l) holds) and the property of the Fourier transform of a translated function reminded in Appendix B: c m,n = (k,l)∈Z se -2iπ(•m/M +•n/N ) (k, l) (s) (k + m/M, l + n/N ) (51)

Acknowledgments. This work is partially funded by GdR CNRS ISIS (Timex project). The authors acknowledge useful discussions with Johan Hoefnagels. They also thank Mohamed Cheikh Teguedi for his assistance in the experimental part of the study.


Roland Becker

Erik Burman

Peter Hansbo

A hierarchical NXFEM for fictitious domain simulations

We suggest a fictitious domain method, based on the Nitsche XFEM method of (Comput. Meth. Appl. Mech. Engrg 2002; 191:5537-5552), that employs a band of elements adjacent to the boundary. In contrast, the classical fictitious domain method uses Lagrange multipliers on a line (surface) where the boundary condition is to be enforced. The idea can be seen as an extension of the Chimera method of (ESAIM: Math. Model Numer. Anal. 2003; 37:495-514), but with a hierarchical representation of the discontinuous solution field. The hierarchical formulation is better suited for moving fictitious boundaries since the stiffness matrix of the underlying structured mesh can be retained during the computations.

Our technique allows for optimal convergence properties irrespective of the order of the underlying finite element method.

KEY WORDS: Nitsche's method; fictitious domain; extended finite element method 1. the (structured) mesh on which the problem is set up and 2. a narrow band of elements that overlays the first mesh.

INTRODUCTION

Fictitious domain methods were introduced in order to be able to use Cartesian meshes also for solving problems on domains with complex boundaries. The idea is to enforce Dirichlet boundary conditions on a given curve (surface) that is discretized independently of the mesh, cf. Glowinski et al. [START_REF] Glowinski | A fictitious domain method for Dirichlet problem and applications[END_REF]. On this curve, the boundary condition is enforced, typically by the use of Lagrange multipliers. The system of equations can then be set up on a Cartesian mesh and the degrees of freedom falling outside of the boundary are discarded. The problem with this approach is that the derivatives of the finite element solution normal to the curve cannot accommodate the jump necessary to achieve optimal order convergence, cf. [START_REF] Girault | Error analysis of a fictitious domain method applied to a Dirichlet problem[END_REF]. Another problem is how to choose the discretization of the curve relative to the elements it crosses in order for the problem to be well posed. Guidelines for this are given in [START_REF] Girault | Error analysis of a fictitious domain method applied to a Dirichlet problem[END_REF] but they are by necessity rather vague.

In this paper we introduce an alternative method based on the use of Nitsche's method in the vein of Hansbo et al. [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF] (building on [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF][START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF]), where overlapping meshes were considered. We shall also employ overlapping meshes in the form of (Figure 1) This allows for the direct use of the method proposed in [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF], where the elements on the underlying mesh were cut by the overlying and the solution pasted together by the use of Nitsche's method [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF]. The outer boundary of the band can then be used as the Dirichlet boundary. We remark at this point that another strength of the approach in [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF] is that any boundary condition can be applied at the outer boundary of the band. This is not so straightforward in a classical fictitious domain method.

However, in [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF] this was achieved by modifying the elements of the underlying mesh, which does not allow for the system matrix on the underlying mesh to remain unchanged. If the boundary was to move continuously, the entries in the system matrix would also have to be changed continuously. It is desirable to have a fixed matrix for the underlying problem and to see the imposition of the boundary condition as a set of additional degrees of freedom (as in the original fictitious domain approach where the additional degrees of freedom are the Lagrange multipliers).

As has been noted by Areias and Belytschko [START_REF] Areias | A comment on the article: a finite element method for simulation of strong and weak discontinuities in solid mechanics[END_REF], the method of [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF] can be interpreted as an extended finite element method (XFEM) method (which instead adds degrees of freedom in a hierarchical fashion) by a reordering of the degrees of freedom. In this paper, we seize on the hierarchical concept to construct a fictitious domain method which is formally identical to that of [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF] (thus benefitting from the optimal order error analysis therein) while still keeping the underlying system matrix unchanged by the location of the boundary.

Our method shares some characteristics with the fat boundary method proposed by Maury [START_REF] Maury | A fat boundary method for the Poisson problem in a domain with holes[END_REF], which also makes use of an auxiliary band-like domain. However, in [START_REF] Maury | A fat boundary method for the Poisson problem in a domain with holes[END_REF], the original Poisson problem is recast on the continuous level and requires a fixed point iteration scheme to solve, like in classical domain decomposition methods. In contrast, our method is defined on the discrete level and can be solved monolithically.

MODEL PROBLEM AND FINITE ELEMENT SPACES

The continuous problem

As a model problem, we consider the Poisson equation

-u = f in , u = 0 on * , ( 1
)
where is a domain in R 2 with polygonal boundary * and f is a given forcing term. We embed in a larger rectangular domain ˆ so that is completely contained in the interior on ˆ . Finally, we introduce a third domain 1 consisting of a band (the width of which may be mesh dependent) whose outer boundary coincides with * and whose inner boundary forms a line, denoted , in the interior of . The remainder of is denoted by 2 := \ 1 , cf. Figure 2). Clearly, the extension to three dimensions is straightforward.

We now rephrase the problem (1) as an interface problem. For u in 1 ∪ 2 we define the jump of u on by [u] := u 1 | -u 2 | , where u i = u| i is the restriction of u to i . Conversely, for u i defined in i we identify the pair {u 1 , u 2 } with the function u which equals u i on i . Then we may formulate (1) as:

-u = f in 1 ∪ 2 , u = 0 on * , [u] = 0 on , [* n u] = 0 on . ( 2
)
Here n is the outward pointing unit normal to 1 and * n v = n•∇v.

For a bounded open connected domain D we shall use standard Sobolev spaces H r (D) with norm ||•|| r,D . The inner products in H 0 (D) = L 2 (D) are denoted as (•, •) D .

Finite element spaces

Denote by T h

1 the triangulation of 1 and by T h 2 the triangulation of ˆ . We shall make a discretization on the whole of ˆ even though the solution has no physical significance outside * . This is in line with classical fictitious domain methods and means that the stiffness matrix assembled from T h 2 remains fixed even if the domain should change, as it must do in many applications. We will use the following notation for mesh-related quantities. Let h K be the diameter of an element K ∈ T h i and h = max K ∈T h i ,i=1,2 h K . To distinguish elements from the two meshes, we will sometimes use indexed element notation K i ∈ T h i for clarity. Furthermore, we introduce G h as the set of elements in T h 2 intersected by ,

G h = {K ∈ T h 2 : K ∩ ⊂ ∅},
and the corresponding mesh-dependent domain

G h = G h K .
We shall also need the mesh-dependent boundary * G h , which consists of the edges on elements in G h that form the boundary of G h . This boundary is also split in the part exterior to , * ext

G h and interior to , * int G h .
The nodes on of the elements in T h 1 , together with the points of intersection between elements in T h 2 and , define a partition of , = j∈J h j . Note that each part j belongs to two elements, one from each mesh. We denote these elements by K j 1 and K j 2 , respectively. A local meshsize on is defined by:

h(x) = h K j 1 , x ∈ j .
(3)

For any element K ∈ T h i , let P K = K ∩ i denote the part of K in i . We make the following assumptions regarding the meshes: A1: The triangulations are non-degenerate, i.e.

h K / K C ∀K ∈ T h i , i = 1, 2
, where h K is the diameter of K and K is the diameter of the largest ball contained in K . A2: The meshes have locally compatible meshsize over . More precisely, let

K j 1 ∈ T h 1 and K j 2 ∈ T h
2 be the elements which contain a specific part j of . We assume that ch K j

1 h K j 2 Ch K j 1 ∀ j ∈ J h .
Here and below, C and c denote generic constants. We shall seek a discrete solution

U =(U 1 ,U 2 ,U 3 ,U 4 ) in the space V h =V h 1 ×V h 2 ×V h 3 ×V h 4 ,
where

V h 1 = { ∈ H 1 ( 1 ) : | K is a polynomial of degree p ∀K ∈ T h 1 , | * = 0}, V h 2 = { ∈ H 1 ( ˆ ) : | K is a polynomial of degree p ∀K ∈ T h 2 , | * ˆ = 0}, V h 3 = { ∈ H 1 ( G h ) : | K is a polynomial of degree p ∀K ∈ G h , | * int G h = 0}, V h 4 = { ∈ H 1 ( G h ) : | K is a polynomial of degree p ∀K ∈ G h , | * ext G h = 0}.
Note that functions in V h are, in general, discontinuous across ; the discontinuity is represented by the hierarchical space V h 3 × V h 4 . In Figure 3 we illustrate the concept: a discontinuous function on a one-dimensional element occupying the domain (0, 1) (solid line) can be written as the sum of a continuous function (dashed line) from V 1 h and piecewise continuous functions which are zero in the nodes of the element (dash-dotted line) from V h 3 and V h

The finite element method (FEM)

The method is defined by the variational problem: find U ∈ V h such that

a h (U, ) =l() ∀ ∈ V h , (4)
where, if we denote int

G h := G h ∩ 2 and ext G h := G h \ int G h , a h (U, ) = (∇U 1 , ∇ 1 ) 1 +(∇U 2 , ∇ 2 ) ˆ +(∇U 2 , ∇ 3 ) int G h +(∇U 2 , ∇ 4 ) ext G h +(∇U 3 , ∇ 3 ) int G h +(∇U 4 , ∇ 4 ) ext G h +(∇U 3 , ∇ 2 ) int G h +(∇U 4 , ∇ 2 ) ext G h -(* n U 1 , []) -([U ], * n 1 ) +( h -1 [U ], []) , l() = ( f, 1 ) 1 +( f, 2 ) ˆ +( f, 3 ) int G h +( f, 4 ) ext G h
, with f extended, e.g. by zero outside , and where h is the local meshsize (3). Here, the jump

[U ] is interpreted as U 1 -(U 2 +U 3 )
. The continuity conditions of u and * n u at are satisfied weakly by means of a variant of Nitsche's method [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF] for consistent weak enforcement of Dirichlet boundary conditions. To ensure stability, the parameter has to be taken sufficiently large, cf. [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF].

In order to show consistency, we first have to define in what sense the method is consistent since the bilinear form a h (•, •) is defined also outside of . Thus, we introduce a second bilinear form

a h * (U, ) = (∇U, ∇) 1 ∪ 2 -(* n U 1 , []) -([U ], * n 1 ) +( h -1 [U ], []) ,
and right-hand side

l * () = ( f, ) 1 ∪ 2 .
It is straightforward to show that the method analyzed in [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF] would then read: find U * ∈ V h such that

a h * (U * , ) =l * () ∀ ∈ V h , ( 5
)
where

U * | 1 =U 1 , U * | 2 =U 2 +U 3 , | 1 = 1 and | 2 = 2 + 3
, since we may ignore the solution outside we set 2 | ˆ \( 2 ∪ G h = 0. Note that the integrals here are taken over only the domain , unlike the proposed method. However, since the solution is completely decoupled by the cut, it does not matter what we do outside of the domain (e.g. how we extend f and how the boundary conditions on * ˆ are specified). Explaining the method of Areias and Belytschko [START_REF] Areias | A comment on the article: a finite element method for simulation of strong and weak discontinuities in solid mechanics[END_REF] without the use of function spaces could be done as follows: only the standard FEM is used on the mesh T h 1 on the band, and on the part of the mesh T h 2 on ˆ which is inside . Cut the elements of T h 2 intersecting and discard all elements completely outside . The solution on these two meshes is then joined by Nitsche's method.

In the next section we prove that indeed U * =U | and hence the analysis of [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF] carries over to our formulation.

A PRIORI ERROR ESTIMATES

Consider the following mesh-dependent norms:

v 2 1/2,h, := h(x) -1/2 v 2 0, = j∈J h h -1 K j 1 v 2 0, j , v 2 -1/2,h, := h(x) 1/2 v 2 0, = j∈J h h K j 1 v 2 0, j ,
and

|v | 2 := ∇v 2 0, 1 ∪ 2 + * n v 1 2 -1/2,h, + [v] 2 1/2,h, .

Proposition 3.1

Let U denote the solution of (4) and U * the solution of (5). Then U | =U * .

Proof

The proof proceeds in three steps 1. Show the existence of a subspace Ṽh ⊂ V h such that a h (U, ) = a h * (U, ) and l( ) =l * ( ) ∀ ∈ Ṽh .

Show that

∃ U ∈ Ṽh such that U | =U and that ∃ U * ∈ Ṽh such that U * | =U * .
3. Apply coercivity and Galerkin orthogonality to the discrete error using the results of [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF].

First, let

Ṽh := { ∈ V h : 2 | ˆ \( 2 ∪ ext G h ) = 0 and 4 | ext G h = -2 | ext G h }.
Since the integrals on are the same in a h (•, •) and a h * (•, •) we only need to prove the equivalence of the volume integrals. By using the decomposition [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF] where B (U, ) denotes the integrals over . It then follows that for all ∈ Ṽh

ˆ = ˆ \( 2 ∪ ext G h )⊕( 2 \ int G h )⊕ int G h ⊕ ext G h , we may write for all ∈ V h a h (U, ) = (∇U 1 , ∇ 1 ) 1 +(∇U 2 , ∇ 2 ) ˆ \( 2 ∪ ext G h ) +(∇U 2 , ∇ 2 ) 2 \ int G h +(∇(U 2 +U 4 ), ∇( 2 + 4 )) ext G h +(∇(U 2 +U 3 ), ∇( 2 + 3 )) int G h + B (U, ),
a h (U, ) = (∇U 1 , ∇ 1 ) 1 +(∇U 2 , ∇ 2 ) 2 \ int G h +(∇(U 2 +U 3 ), ∇( 2 + 3 )) int G h + B (U, ) = (∇U 1 , ∇ 1 ) 1 +(∇(U 2 +U 3 ), ∇( 2 + 3 )) 2 + B (U, ) = a h * (U, ). (7)
The equality l( ) =l * ( ) is shown in a similar fashion. Second, observe that Ṽh only imposes constraints on components of the function outside , since the constraint on 2 in int G h is compensated for by the freedom of 3 . Hence, the existence of the sought U and U * .

Finally we recall the following coercivity result from [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF], for some c>0 and for all v ∈ V h there holds

c |v | 2 a h * (v, v). ( 8
)
In particular this holds for v =U -U * and hence, since

U -U * ∈ Ṽh c |U -U * | 2 a h * (U -U * ,U -U * ) = a h * (U -U * , U -U * ) = a h (U, U -U * )-a h * (U * , U -U * ) =l( U -U * )-l * ( U -U * ) = 0.
We have the following consistency relation.

Proposition 3.2

The discrete problem ( 4) is consistent in the sense that, for u solving (2) there holds

a h * (u, ) =l * () ∀ ∈ V h ,
or, equivalently,

a h * (u -U, ) = 0 ∀ ∈ V h . ( 9
)
The proof is given in [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF].

We can now directly take advantage of the theory developed in [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF] which shows that we have optimal error estimates for any polynomial degree of the underlying FEM: Theorem 3.3 For U solving (4) and u solving (2), the following a priori error estimates hold

|u -U | Ch p |u| p+1, , (10)
and

u -U 0, Ch p+1 |u| p+1, . (11)
We refer to [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF] for the proof. Here we shall only point out one of the crucial points in the analysis. Accuracy of the method is expressed by the orthogonality relation ( 9), but to show convergence we also need stability of the discrete problem as expressed by [START_REF] Maury | A fat boundary method for the Poisson problem in a domain with holes[END_REF]. In order to show that our system matrix is positive definite, we rely on the following inverse inequality, see, e.g. Warburton and Hesthaven [START_REF] Warburton | On the constants in hp-finite element trace inverse inequalities[END_REF].

Lemma 3.4

For ∈ V h , the following inverse inequality holds:

* n 1 -1/2,h, C I p ∇ 0, 1 .
The size of the constant C I can be found by solving a small local eigenvalue problem; explicit bounds are discussed in [START_REF] Warburton | On the constants in hp-finite element trace inverse inequalities[END_REF]. The constant C I determines the size of , we are obliged to take >C 2 I p 2 in order to ensure coercivity. Here is the point of the band: if we just cut the mesh with * and apply Nitsche's method on the cut elements, sliver elements would be generated that would require extremely large C I in Lemma 3.4, leading to severe ill conditioning of the discrete system. We illustrate the problem in the case p = 1. Thus, consider using only V h 1 and cutting the mesh with * . Then C I can be found, for a given element K as the largest eigenvalue max in the eigenproblem of finding U ∈ V h 1 and ∈ R such that

(h 1/2 K * n U, * n v) K ∩* = (∇U, ∇v) K ∩ ∀v ∈ V h 1 .
For p = 1, ∇v is constant on each element and thus we have

* n v 2 L 2 (K ∩* ) = meas(K ∩* )|* n v| 2 , ( 12
)
where meas(•) denotes the length, area, or volume of the object in question, and

∇v 2 L 2 (K ∩ ) * n v 2 L 2 (K ∩ ) = meas(K ∩ )|* n v| 2 (13)
and it follows that

h 1/2 K * n v 2 L 2 (K ∩* ) h K meas(K ∩* ) meas(K ∩ ) ∇v 2 L 2 (K ∩ ) , (14)
and thus we must choose The situation is illustrated in Figure 4 where * is represented by the dashed line. We note that as * moves to the left in Figure 4, the area meas(K ∩ ) will approach zero, whereas the length meas(K ∩* ) remains bounded from below, which means that must grow without bound. This situation is remedied by inserting the band of elements between the cut and the boundary.

>C 2 I h K meas(K ∩* ) meas(K ∩ ) .
In Burman and Hansbo [START_REF] Burman | Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method[END_REF], where no band was used, this conditioning problem was instead eliminated by the use of additional stabilizing terms, limiting the analysis to linear polynomials. In the present formulation, we avoid the use of additional terms since the mesh on the band 1 is always shape regular. We also remark in this context that the use of higher-order polynomials does not require a curved interface since the interface is artificial. Thus, the interface can always be assumed to be piecewise affine and no further complexity is introduced by higher-order approximations, cf. also [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF].

IMPLEMENTATION ASPECTS

The difference between the suggested approach and that of [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF] is in the new hierarchical formulation. Conceptually, this means that in the method of [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF], there are only two finite element spaces: the one on the band and the one inside the band consisting of elements up to the interface . This gives two sets of unknowns, u 1 and u 2 say, and the system of equations becomes

S 1 B B T S 2 u 1 u 2 = f 1 f 2
where S 1 is the stiffness matrix from the discretization on the band, S 2 from the cut mesh, B represents the coupling terms, and f 1 , f 2 represents load terms. The degrees of freedom in u 2 outside the cut can then be discarded already at the outset. One problem with this approach is that all matrices B, S 1 , and S 2 change if we want to move the interface. This could be in a time-dependent problem, or if we want to use the scheme for the purpose of shape optimization.

In the present approach we instead obtain the system

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ S 1 B 12 B 13 B 14 B T 12 S 2 B 23 B 24 B T 13 B T 23 S 3 B 34 B T 14 B T 24 B T 34 S 4 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ u 1 u 2 u 3 u 4 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ f 1 f 2 f 3 f 4 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦
where S 2 now denotes the stiffness matrix on the underlying mesh which does not change. Since this matrix is by far the largest of the involved matrices, this means that we only have to recompute small matrices (corresponding in a sense to the multiplier matrices of the original fictitious domain method).

In our implementation, two different approximations of the interface have been used: the first defining the approximation of given by the piecewise affine (internal) boundary of T h 1 , the second defined by the approximation of the intersection of with the elements on T h 2 (which defines a different affine approximation of ). The first has been used as a master surface for the integration of jump and consistency terms and should also, for exactness, be used for the computation of the bilinear form a h (•, •). We have however, for simplicity, used the second affine approximation of the interface for computing a h (•, •) which removes some technical difficulties in the implementation but introduces a conformity error. This error does not, however, affect the convergence rates in our example below. For higher-order methods, the intersection of the meshes must be handled more carefully, although, as noted above, the interface can remain piecewise affine since it is given by the user-defined internal boundary of T h

NUMERICAL EXAMPLE

We consider a problem posed on a disc of radius r 0 = 0.95. With r the length of the radius vector, we use f =r to obtain the exact solution u = u = (r 3 0 -r 3 )/9. The stabilization parameter was set to = 10. In Figure 5 we show the obtained convergence rates in L 2 ( )-and H 1 ( )-norms, which are optimal. An elevation of the solution is given in Figure 6, and an elevation of the solution on the whole of ˆ is shown in Figure 7. Note that we have extended f =r to hold on the whole of ˆ and imposed zero boundary conditions on * ˆ . This is of no consequence since the solution is decoupled at .

CONCLUDING REMARKS

In this contribution we have shown that the NXFEM method is well suited for fictitious domaintype simulations. It has optimal convergence for arbitrary polynomial order and does not require Lagrange multipliers to enforce Dirichlet boundary conditions. Indeed, since the boundary conditions are prescribed on a regular mesh, we can handle all types of boundary conditions in the usual way.

Applications for the proposed method include shape optimization, where the boundary of the domain has to be moved in order to calculate sensitivities, and for computations involving objects moving across a background mesh in general.

Figure 1 .

 1 Figure 1. Underlying and overlying meshes.

Figure 2 .

 2 Figure 2. Boundary and domain definitions.

Figure 3 .

 3 Figure 3. A discontinuous trial function and its split into a continuous and a discontinuous part.

Figure 4 .

 4 Figure 4. An element cut by the boundary.

Figure 5 .

 5 Figure 5. Convergence in broken H 1 and L 2 norms.

Figure 6 .

 6 Figure 6. Elevation of the computed solution on .

Figure 7 .

 7 Figure 7. Elevation of the solution on ˆ .


Silvina Carrizo
email: scarrizo@conicet.gov.ar

Sébastien Velut
email: sebastien.velut@univ-paris3.fr

des sciences du territoire -Frontiers and boundaries of territorial sciences Energy Networks and Regional Integration in South America AUTEURS

Keywords: Énergie, réseaux, intégrations Energy, networks, integration

Dans la dernière décennie, les pays d'Amérique du Sud ont connu une croissance économique robuste fondée sur l'exportation de matières premières, ce qui implique une consommation croissante d'énergie provenant de différentes sources. Les États ont dû aussi faire face à deux défis : les inégalités sociales et territoriales ainsi que les grandes distances et les faibles densités de population. Le premier défi implique que les États doivent mettre en place des politiques publiques pour apporter des réponses spécifiques aux groupes sociaux vulnérables. Le second entraîne des problèmes particuliers d'aménagement du territoire et de construction d'infrastructures. Les relations entre États, dans le cadre politique de l'UNASUR, font jouer un rôle central à l'énergie pour soutenir la croissance et réaliser l'intégration régionale. Cet article montre que, si ces programmes ont apporté des solutions pour fournir de l'énergie bon marché à des grands programmes productifs, ils n'ont pas donné de solutions adéquates à une autre série de problèmes sociaux et géopolitiques.

INTRODUCTION

Regional integration of South American countries has already a long history. Political instability, rivalries between countries as well as the ups and downs of economic cycles led to phases of integration, such as in the beginning of the 1990s with the creation of the Mercosur (Common Market of the South), as opposed to phases of dismantlement of existing mechanisms, like the present terminal crisis of the CAN (Community of Andean Nations). In contrast to these political fluctuations, supra national integration of some precise sectors is sometimes seen as a pragmatic way to build stronger ties between countries, based on common interests. Since the 1970s, energy is one of those sectors, probably the most important for two main reasons. The first is that South America as a whole has a rich energy endowment, ranging from fossil fuels to renewable energies, whose geographic distribution differs from the geography of consumption, making necessary connections between regions and countries. The second reason is that such connections and energy facilities represent long term investments. They are built to stay and operate for several decades, shaping the relationships between spaces and between economic agents.

According to a document published in 2012 by UNASUR (Union of South American Nations) and OLADE (Latin American Energy organization), energy integration should become a tool for social and economic development and contribute to eradicate poverty. The universal access to energy is seen as a citizen's right and it is one of the main challenges for South America, where social and spatial inequality is high and persistent. Most governments implement public policies in order to give specific energy solutions to vulnerable social groups and to connect isolated regions to energy networks. Economic subsidies to companies, social tariffs for the poor and building of large infrastructures across low density regions are the common actions tackled by South American governments. Within the political framework of UNASUR, created in 2010, governments have put energy regional integration in pivotal role to sustain economic growth and social welfare. This paper exposes a set of elements about economic and social challenges in South American countries and energy regional integration. It discusses the results of energy integration policies in relation to two main goals: economic development and social integration.

ENERGY AND DEVELOPMENT

South American exports have risen steadily during the last decade due to the opening of trade barriers, companies outsourcing and the momentum of Asian economies, especially China. Mining activities have expanded in order to satisfy external demands. Only 1/4 of regional exports stays in the subcontinent, which means that exchanges between South American countries remain at a low level. Exports are going principally to North America and Europe, but Asia represents a substantial and growing destination for many countries particularly Brazil and Argentina. Commodities such as ores, metals and agricultural products, notably soy beans, represent the main exported products. Only Brazil is also significantly exporting industrial products, like motor vehicles and electronic devices. These exports mean that South American countries are exporting energy not only because some of them, notably Venezuela and Ecuador, are oil exporters and OPEP members, but also because mining, agriculture and long distance transportation require large quantities of energy. Although no estimations have been made, one can say that exported South American products incorporate a significant energy content and, as a matter of fact, green house gases emissions. Energy availability in good conditions of prices, quality and quantity has become crucial to sustain the current model of economic development.

Agricultural production depends more and more on hydrocarbons, since the so called "green revolution" is based on the use of chemical and machines. Mining and ore processing are important electricity consumer with increasing demands. For instance, in Peru, where mining is growing fast, mines absorb 10% of national power capacity (850 MW of 8,600 MW) and new investments are planned (Latinomineria, 2011). In Chile, copper production takes 1/3 of electricity sold in the country (passing from 27% in 1995 to 34% in 2008). It represents more than 80% of power sold in the Northern regional system and 20% in the Central system, where population is concentrated. Moreover electricity represents half of direct energy consumption for Chilean copper production [START_REF] Cochilco | Demanda de energía eléctrica y seguridad de abastecimiento para la minería del cobre de 08/2007[END_REF][START_REF] Cochilco | Demanda de energía eléctrica en la minería del cobre y perspectivas de seguridad en su abastecimiento de 01/2010[END_REF]. These examples show the enormous quantities of energy needed by mines. Nevertheless, transportation is the main energy user in South America. In most countries, it relies on road transportation for goods and people, since railways networks do not cover the countries and river shipping is limited to some regions. Even in largest South American countries, railway or fluvial transportation is underdeveloped, except on some of the large rivers, like the Amazon and the Parana. As a consequence, most of the energy use by the transportation sector is demanded by trucks carrying products to consumers and exports facilities. Regional integration of transport networks is still to be done.

In spite of the fact that South America has large energy resources used to sustain economic development, energy distribution is still socially unfair. The region shows a high degree of urbanization, above 80% in some countries. It connects to modern energy networks, mainly electricity and gas in an easier way than in countries where the population is disseminated. Although urban energy services are extended, universal energy access has not been accomplished. Energy access is still a problem for two types of populations. Poor city dwellers, and particularly those living in slums, represent the highest figure. People living in low density areas far from large urban centers form the other group. 40 million people lack of access to modern electrical services in Latin American and Caribbean countries (BID, 2012) and many urban inhabitants do not obtain energy in a legal and safe way, using instead illegal and precarious connections to power lines passing by their homes. Bolivia and Peru exhibit the lowest electrification levels : 1/5 of their populations has no connection to power networks. In rural regions, far away from large urban centers, millions of households depend on wood and coal for cooking and heating, and, what's more, through the use of unhealthy and unsecure equipment. Bottled gas is also an important energy source for people that do not have access to modern networks, but its prices increased significantly. Poor people who have access to electricity employ it to heat their house and sanitary water, in an expensive, inefficient and dangerous way.

All South American governments have devised policies for universal access to energy, putting rural communities in first place. Most countries put increasing emphasis on diversifying their energy supply and using renewable energies. They also look for distributed energy production and support investment in order to bring power to isolated or lowincome population. For instance, Argentina operates the PERMER (Programa de energía renovable para mercados rurales/Programme for renewable energy in rural markets) a successful initiative, developed in most Argentinian provinces, that enables to bring energy to isolated population and public buildings (schools, dispensaries, police stations). Thousands of photovoltaic systems installed in localities have changed many people life. Similarly, the Brazilian programme Luz para Todos (meaning "light/power for all"), launched in 2002, has brought energy to millions of rural households in a decade.

In order to contribute to universal energy access, but also to reduce social inequalities,; subsidies are given either to consumers or to companies. Public policies include subsidies for fossil fuels and electricity as a way to control tariffs for final users or help energy companies install equipment. In some cases, public subsidies are universal, aiming to improve social welfare and boost the economy. In other situations, they are restricted to certain sectors of population, productive activities or services like transport. Subsidies are controversial because their impacts on energy efficiency, environment, foreign trade and public finance are negative. Fossil fuel subsidies are costly and socially unfair because rich consumers and large private companies are receiving higher benefits than poor households and small firms. Conversely, charging all consumers with the real cost of energy services might entail social exclusion and trigger political protests. Therefore, energy integration, within the countries and between them, must find its way through many contradictions.

THE REGIONAL INTEGRATION PATCHWORK

Building large scale infrastructure to strengthen energy systems has been a big challenge for South American countries with restricted financial capacity, but also a necessary endeavour to harness huge rivers with dams, exploit fossil fuels in harsh conditions and bring power to far away consumers. In the second half of the 20th century high levels of demand in South American metropolis had justified large investments in energy facilities. Consuming centers import and transform energy to satisfy their needs but also to secure their energy supply. Meanwhile, in many cases, territories that export energy do not have enough resources or proper services. Often the territories crossed by large infrastructures do not benefit from them [START_REF] Carrizo | Aprovisionamiento eléctrico de Buenos Aires y desigualdades regionales entre la metrópolis y el Noreste argentino[END_REF].

Binational hydraulic dams -like Itaipu (Brasil-Paraguay), Yacyreta (Argentina-Paraguay) and Salto Grande (Argentina-Uruguay)-are examples of those projects responding to metropolitan demands, planned in the 1970s. They have considerable impacts on the environment and have been built with little consideration for local people. Moreover, high voltage power lines that cross large distances to take energy to the main consuming region do not provide energy to local consumers adjacent to the lines. For instance, Yacyreta (3,200 MW), built upon the Parana river, is connected to Buenos Aires (700 km) by several power lines; meanwhile the Argentina Northeast -where the dam is locatedexhibits the lowest level of electric service [START_REF] Carrizo | Aprovisionamiento eléctrico de Buenos Aires y desigualdades regionales entre la metrópolis y el Noreste argentino[END_REF]. International dams were historically the first elements for energy integration and they still embody regional problems, not only of territorial inequalities, but also of national asymmetries. For example, for many years Paraguay has been claiming better prices for the electricity from Itaipu and Yacyreta, sold to Brazil and Argentina, and the right to sell it to third countries. In 2011, after a long negotiation, Brazil eventually agreed to triplicate prices but did not allow Paraguay to use its transmission line to sell power to distant markets. Electricity exports represent the main export for Paraguay and the main fiscal income, like gas for Bolivia. For those two poor countries, it is important that energy integration with their richer neighbours improves the way in which they sell their resources and get a fair price helping them to finance needed investments and relieve poverty. If energy integration only serves the interest of richer countries, it can become a source of disagreement between partners.

While hydroelectric binational connections have been planned several decades ago, international gas connections were only conceived at the end of 20th century. In the 1990s, many governments implemented liberal reforms that included the privatization of energy companies, deregulation of activities and decentralization. Multinational companies developed international projects to build larger energy markets. At that time, only Argentina had an important natural gas network that the State started to extend in 1946, building the first national gasline (1,700 km long) to carry gas from Comodoro Rivadavia (Province of Chubut) to Buenos Aires, feeding patagonian and pampean cities located on its way. The State also laid out gas pipelines from the five national hydrocarbon fields (situated at the extreme South of the continent, in central Patagonia on the Atlantic coast, in the Andes in North Patagonia and Cuyo, and in the North of the country) and from Bolivian gas fields. In the 21st century, almost all South American countries extended their gas networks, building pipelines from Bolivia and Argentina to Brazil, and Chile, developing their gas fields and installing liquefied natural gas terminals. Unlike power lines, gas pipelines, between fields and consumption centers, are easy to use by medium-sized towns. Equipment for gas derivation is not as expensive as for electricity. Nevertheless, low population density, scarce economic activity and low incomes of inhabitants often makes it impossible to achieve even small and medium scale projects.

The share of gas is rising in most countries. Significant reserves exist, both conventional and unconventional, scattered all over the continent, but with the exception of Argentina, countries did not use this energy vector until the last two decades. Transnational gas connections have triggered gas use, but they also nurtured the interest for national gas sources. Most countries have built new facilities to satisfy national demands but also for imports and exports. In this context, regional integration stands somewhere between national integration of gas and power networks and participation in the world gas market. The former is consistent with national energy policies seeking to take advantage of all national resources, the latter helps exporting countries to get the best prices, and importing ones to guarantee supply through long term contracts with different providers. South American integration of gas and power networks, in spite of political support and the obvious advantage of cutting distance, seems difficult to achieve due to local resistance to large infrastructure, mistrust of private companies towards governments and the lack of common guidelines for energy policies. Even if regional interconnections have been built in South America, regulation frameworks are still different from one country to another, and between energy subsectors. For instance, differences appear in public and private degrees of participation in energy activities, technologies, standards and in policies aims. Multinational companies' ubiquity, external demands dependence have helped to associate interests and unify some rules. Regional institutions, like OLADE, have tried to promote uniform norms, to negotiate changes in national policies to make them compatible at the South American scale. Nevertheless energy tariffs, prices and taxes differ widely between countries. Differences in tariffs are responsible for local conflicts or tensions in border regions, especially when there are shortages of fuels and smuggling becomes problematic. Prices and taxes are key elements used to stimulate energy production, either for national demands or for exports: cheap energy (e.g., Argentina) improves industry competitiveness and facilitates energy access for low income population but does not stimulate exploration of new resources and energy efficiency. Conversely high prices (e.g., Chile) trigger investment in energy production and make consumers more responsible but they can also become a burden for consumers and firms. Taxes on oil and gas must be carefully crafted in order to provide fiscal revenues without killing the golden egg hen. The lack of dialog on those crucial topics between South American countries, eager to defend their own sovereignty, is a major obstacle to integration. On the political side, UNASUR and multilateral banks are focused on transnational infrastructures: these are certainly useful to make the whole energy system more efficient, but the discussion on energy policies and norms will also have to be tackled.

CONCLUSIONS

Public policies are still mainly national, although they share many common goals and tools: universal access to energy, cheap energy for mines and industries, development of new energy sources. Environmental aspects, such as green house gases emissions, are not on the top of the agenda. Energy efficiency has not been seriously taken into consideration in the most consuming sectors like transportation. On the contrary, many South American countries rely on their large energy endowments from different sources to keep pace with rising consumption. Energy integration between countries is conceived mainly as a matter of infrastructure, devised to bring cheap energy to consumers in large urban regions. Social and spatial inequalities are given specific attention in most countries. Regional integration is still limited to physical exchanges without a common legal framework that might soon prove indispensible. Small countries with high energy endowments perceive physical integration as redolent of colonial intentions from their energy-hungry neighbours. Regional integration in South America must be compatible with bold nationalism exhibited by most countries eager to preserve their own sovereignty, and at the same time take advantage of complementarities between them: a challenge that energy issues exemplify thoroughly.