David Vanderhaeghe

Romain Vergne

Pascal Barla

Bill Baxter

Keywords: Perception, cognition, interaction Interaction et visualisation Équipe-Projet IPARLA Shading, Expressive Rendering, NPR

We define a new shading tool called a Bidirectional Appearance Distribution Function (BADF) tailored to the direct control of stylized appearance. A BADF can be thought of as defining the appearance of a sphere from all possible illumination directions. Our BADF formulation generalizes and improves upon previous stylized shading techniques by enabling the direct control of shading profiles in screen space, exaggerating surface features in a flexible manner, and letting users control stylized appearance from multiple lighting or viewing directions. This allows users to start from a simple shading behavior, and refine from there towards greater stylization. Our GPU implementation works in real-time, which benefits both editing, and rendering in interactive systems. These features make BADFs an efficient tool for many applications in artistic and scientific illustration domains.

Fonction de distribution bidirectionnelle de l'apparence pour l'éclairement stylisé

Résumé : Nous définissons un nouvel outil d'éclairement appelé Fonction de distribution bidirectionnelle de l'apparence (BADF). Cet outil est conçu pour le contrôle direct de la stylisation de l'apparence. Une BADF peut être vue comme la définition de l'apparence visuelle d'une sphère pour toutes les directions de lumières possibles. Notre formulation de la BADF généralise et améliore les méthodes de stylisation précédentes car elle permet à l'utilisateur de contrôler directement le profil des primitives d'éclairement en espace écran, d'exagérer les détails de surface de manière flexible et de lui laisser la possibilité de styliser plusieurs directions de vue et/ou d'éclairage. Ceci permet à l'utilisateur de partir d'un rendu classique et de le modifier pour obtenir un éclairement riche et stylisé. Notre implémentation sur carte graphique permet un rendu en temps réel, autorisant l'édition interactive et l'utilisation de notre méthode dans des contextes de rendu temps réel. Les BADFs sont donc un outil efficace pour les applications artistiques et de visualisation scientifique.

Mots-clés : Eclairement, Rendu expressif, NPR

Introduction

Shading in Computer Graphics is generally designed to simulate reality by using a combination of complex materials and lighting environments. Although this approach has many applications in architecture or special effects for movies, other applications rather use shading in a stylized manner. In scientific illustrations (e.g., archeology, medicine, botanics, anatomy, mechanics), shading depicts materials and illumination in simplified ways [START_REF] Wood | Scientific Illustration: A Guide to Biological, Zoological[END_REF][START_REF] Hodges | The Guild Handbook of Scientific Illustration[END_REF] to communicate clearly object shape or function. Whereas in artistic illustrations (e.g., comics, anime, graphic novels, illustrated books), conventional shading is radically altered, for instance by creating sharp color transitions that only give a hint of material or illumination. Such highly stylized shading is chosen to set different moods or atmospheres, or exaggerate object structure [START_REF] Hogarth | Dynamic Light and Shade[END_REF][START_REF] Mccloud | Understanding Comics: The Invisible Art[END_REF].

With hand-made stylized shading, material and lighting properties are only implicitly given, and do not need to match physically realistic constraints. These additional degrees of freedom regarding physical accuracy are traditionally used to enhance object shape through shading. Indeed, skilled artists are able to reveal smooth or sharp surface features by creating subtle or strong color variations around them. Such flexibility might even be used to design sophisticated dynamic behaviors, since shading does not need to be consistent across varying viewing or lighting conditions. In scientific illustrations [START_REF] Wood | Scientific Illustration: A Guide to Biological, Zoological[END_REF], variations in shading may provide a means to better depict object parts that are revealed under specific viewpoints. In artistic illustrations [START_REF] Hogarth | Dynamic Light and Shade[END_REF], variations in shading may change the perception of a scene depending on the overall illumination direction.

The main goal of this paper is to provide an approach to stylized shading for both scientific and artistic illustrations. Ideally, users should be able to design stylized shading appearance without having to deal with complex materials or illumination at all. This is the main idea behind the LitSphere of Sloan et al. [START_REF] Sloan | The lit sphere: A model for capturing NPR shading from art[END_REF], which is the previous method most similar to our work. It consists of a single shaded sphere stored in a 2D texture and looked up using screenspace normals. An artist simply draws the desired appearance using painting software or captures it from existing paintings. Variations on this basic idea have been successfully used in different contexts such as volumetric rendering [START_REF] Bruckner | Style transfer functions for illustrative volume rendering[END_REF] or digital sculpture software such as ZBrush ® .

However, because of its reduced dimensionality and simple look-up function, a LitSphere is limited with regards to our target applications: 1) Diffuse or glossy materials cannot properly be imitated since the illumination orientation is "baked" into the model; 2) Disturbing sliding artifacts may occur when high-frequency variations are found in the texture; 3) Shape features are only enhanced by using special (e.g., darker) colors along contours.

Our main contribution to stylized shading is the introduction of Bidirectional Appearance Distribution Functions. A BADF explicitly models the stylized appearance of a shaded sphere under every possible illumination orientation. We show in Section 3.1 that it can be defined as a 4D function that takes a normal and a view direction as input and outputs a color. Increasing dimensionality is the key to address the main limitation of LitSpheres: since multiple illumination orientations are taken into account, a BADF can properly represent the characteristic effects of opaque materials such as diffuse gradients or highlights, and the way they move. Defining such a 4D function with an image-based approach would not be tractable because it would require too much user interaction. Instead, we introduce a generalized Phong model which allows artists to easily achieve a wide range of stylized shading appearances. As explained in Section 3.2, it is able to reproduce Phong, Half-Lambert and cartoon shading as special instances. Our model also permits to depict salient surface details by tweaking shading as detailed in Section 3.3, which would not be possible with an image-based technique. Such an approach enhances shape through shading in ways that more closely resemble hand-made illustrations compared to previous work.

With these features, it takes very few steps to create interesting BADFs. In addition, for more advanced shading behaviors, users may design predefined stylized appearances, and let our model interpolate appearances between key viewing or lighting directions, as explained in Section 3.4. With directional control, the more time spent adding key directions, the richer the resulting stylized appearances can be. We show examples of artistic and scientific illustrations obtained with our technique in Sections 4 and discuss its limitations in Section 5.

Previous work

The earliest form of direct appearance control we are aware of is inverse rendering. Many methods have been proposed, some enabling material editing [BAOR06, CPK06, EBJ * 06], others for the control of illumination [KPC93, SDS * 93, RH01]. The inverse rendering method closest to our approach is the Illumination Brush of Okabe et al. [START_REF] Okabe | Illumination Brush: Interactive Design of All-Frequency Lighting[END_REF]. Starting from a known BRDF, they infer the environment lighting by directly painting on the object the desired result of either diffuse or specular reflections. Even if the paint strokes may be moved around for stylization purpose, the adaptation of the method to the creation of stylized shading is not straightforward, since it relies on BRDF models and realistic illumination and materials. We rather seek a method that is not bound by such physical constraints.

Many stylization techniques have based their approach on a modification of an existing shading model. The simplest one is Toon Shading, which simply segments the result of a diffuse shading function into two or three color bands. The Technical Illustration Shader of Gooch et al. [START_REF] Gooch | A nonphotorealistic lighting model for automatic technical illustration[END_REF] uses a similar approach with an unclamped diffuse term (also called Half-Lambertian) and cool-to-warm color gradients. It has also been found to offer better depiction of shape in video game applications [START_REF] Mitchell | Illustrative rendering in Team Fortress 2[END_REF]. The X-Toon method of Barla et al. [START_REF] Barla | X-Toon: An extended toon shader[END_REF] extends the 1D toon function to a 2D function stored in a texture. The additional dimension is used for creating back-lighting or depth-offield effects, for instance. However, even with local control [START_REF] Todo | Locally controllable stylized shading[END_REF], the range of styles provided by remapping a pseudo-diffuse function to a color gradient is limited.

These shading methods have been used as building blocks in more recent work that focus on shape depiction through shading. The Exaggerated Shading of Rusinkiewicz et al. [START_REF] Rusinkiewicz | Exaggerated Shading for Depicting Shape and Detail[END_REF] employs a Half-Lambertian and a hierarchy of simplified normals to align the light direction with grazing angles. It creates an accurate depiction of surface details, but is tied to the use of a Half-Lambertian shading. In contrast, the Apparent Relief technique of Vergne et al. [START_REF] Vergne | Apparent relief: A shape descriptor for stylized shading[END_REF] INRIA is not correlated to a specific shading model, but it only uses shape features as a mask to blend between a pair of shading appearances.

Other techniques modify Phong shading in order to manipulate highlights specifically. The system of Anjyo and Hiramitsu [START_REF] Anjyo | Stylized Highlights for Cartoon Rendering and Animation[END_REF] provides control over the shape of cartoon highlights: they interpolate parameters of highlight transformations across time. Existing operations include translation, rotation, scaling, splitting, squaring, and boolean operators. The method has been extended [START_REF] Anjyo | Tweakable light and shade for cartoon animation[END_REF] to control the location of cartoon diffuse shading (modifying light direction), and to manipulate directly the shading on an object's surface. However, the manipulations are object-specific, and thus do not easily transfer between different 3D objects. A similar approach for highlights has been subsequently proposed by Pacanowski et al. [START_REF] Pacanowski | Sketch and Paint-based Interface for Highlight Modeling[END_REF]. Their method is not limited to cartoon highlights, but works for any star-shaped highlight with its appearance defined by a color gradient. They also propose a simple technique for varying appearance depending on a few predefined illumination directions.

Our BADF model encompasses techniques based on Lambertian (including half-Lambertian) or Phong shading models and incorporates surface shape features by applying complex shading variations without any apriori on shading. The resulting stylized appearance is controlled from arbitrary key lighting or viewing directions, and works in real-time on modern graphics hardware. It is also easily transferred between 3D objects, as opposed to methods that require establishing a specific surface parametrization.

BADF model

One of the most important considerations to an artist when designing appearance is the location of the dominant light source [START_REF] Wood | Scientific Illustration: A Guide to Biological, Zoological[END_REF]. Given this information, artists are able to shade surface regions depending on their orientation relative to that light, and place secondary lights to create more complex effects. Our approach is to define an illumination reference frame L, and to define all light sources in this frame. End-users then manipulate L in the same way they would manipulate an environment illumination map.

Definitions

The main idea of a BADF is to encapsulate the definition of L directional light sources l i , their interaction with the object material, and stylistic effects, into a single function. Since all light sources are expressed in L, we also express the normal n and view v directions in L and define a BADF as a 4D function of n and v:

ρ A : S 2 × S 2 → [0, 1] 4 n, v → (r, g, b, a) (1)
where S 2 is the sphere of directions.

In this paper, we introduce a layered BADF model that takes inspiration from image and video editing software such as Photoshop ® or Shake ® (see Figure 1 and the supplemental video). Each layer is assigned a light source direction l i , a shading function ρ i and a color c i . The BADF is obtained by combining layers together with user-specified blending functions ⊙ i , one per layer. We use additive, multiplicative, subtractive and alpha blending functions in our system. More formally, our BADF model is given by:

ρ A (n, v) = C L (n, v) C i (n, v) = (ρ i (n, v) c i ) ⊙ i C i-1 (n, v) (2)
where

C 0 (n, v) is a constant background color and i ∈ [1, L].
If, for all i, we choose ⊙ i to be the additive blending function and

ρ i (n, v) = ρ(l i , v)(n•l i )
with ρ a Bidirectional Reflectance Distribution Function (BRDF), then we fall back to the "conventional" shading model. In stylized shading though, it is more interesting for artists to try out different combinations of layers, and to make use of shading functions that depart from realistic BRDFs. Our BADF model provides such flexibility via three key technical contributions: a shading function (Section 3.2) directly controllable in image-space via isocurve widgets; a warping method (Section 3.3) that enhances surface features by manipulating layer color and shading function; and a directional interpolation technique (Section 3.4) that varies layer parameters based on key lighting or viewing directions. When used in combination, these techniques provide for a wide range of stylized appearances, as shown in Figures 6 and7.

Shading function

In our model, each of the L directional light sources not only conveys material effects such as diffuse gradients or highlights, but also stylistic effects such as sharp, cartoon-like transitions, or Half-Lambert fall-offs. This is obtained via the shading function ρ i , defined as the composition of an angular parametrization u i : S 2 × S 2 → [0, π] and a profile function α i : [0, π] → [0, 1]:

ρ i (n, v) = α i • u i (n, v).
(3)

The choice of parametrization u i controls the way shading moves onto the surface and is defined relative to the lighting direction l i : where d i may take on different values depending on desired motion behavior: d i = n in the diffuse case and d i = r in the specular case, with r the mirrored view direction. We also found useful to introduce a screen-aligned parametrization, as with Lit Spheres, to create specific effects. In this case, the parametrization ignores the illumination reference frame orientation. This simple behavior is obtained by setting d i = n with n the normal expressed in screen-space. Figure 2 shows the three parametrizations. The choice of profile function α i controls the way shading intensity falls off away from l i (or u i = 0). We use a profile function that is a generalization of Phong shading:

u i (n, v) = acos(d i • l i ) (4) INRIA 0 π d i = n d i = r d i = n
α i (u) = |β i + (1 -β i ) cos(ζ(u -τ i ))| + γi (5)
where ζ clamps values to [0, π] to ensure a monotonically decreasing function, and τ i , β i and γ i are offset, bias and exponent parameters. Offset controls the angular portion of the parametrization where the shading function is equal to 1. Bias provides a control similar to the Half-Lambertian term in [START_REF] Gooch | A nonphotorealistic lighting model for automatic technical illustration[END_REF][START_REF] Mitchell | Illustrative rendering in Team Fortress 2[END_REF]. Exponent determines the fall-off rate of shading, and is traditionally used in Phong's model to control material shininess.

In practice, the β i and γ i parameters are not easy to manipulate. We thus introduce alternative fall-off and cut-off controls (see Figure 3-left). The falloff f i ≥ 0 controls the parametric location of the profile function inflection point, while the cut-off c i ≥ f i controls the parametric location at which the profile reaches 0. To evaluate u i using these alternative controls, we set β i = β(c i ) and γ i = γ(f i , c i ), as described in Equations 10 and 11 in Appendix A. With these parameters, artists directly control shading appearance with isocurve widgets (see Figure 3-right), instead of having to tune unintuitive parameters. Conventional shading profiles are easily reproduced with this function: Phong's diffuse term corresponds to τ i = 0, f i = c i = π/2; Half-lambertian is same as diffuse, except that c i = π; Phong's specular term is same as diffuse, but leaves f i free to vary in the [0, π/2] range; cartoon shading occurs when τ i = c i , and ambient shading when τ i = π (used as a fill-in layer in our results).

Surface features

In traditional artistic and scientific illustrations, it is common to reveal surface features through the careful use of biased shading variations. Mean Curvature Shading [KWTM03] conveys shape details by assigning a pair of dark and bright colors to surface concavities and convexities respectively. Similarly, a pair of LitSpheres may be interpolated based on curvature information, as in the Z-Brush ® software.

In our approach, we make use of surface curvature not only to modify shading color, but also the extent of shading shape. To this end, we first extract a measure of surface curvature κ from the depicted object. In our system, we use the view-centered mean curvature as defined in the work of Vergne et al. [VPB * 09], because it provides automatic levels-of-detail and is computed dynamically from screen-space normals. However, our approach is independent of the choice of curvature measure. Taking into account surface shape turns a BADF into a 5D function:

ρ A : S 2 × S 2 × R → [0, 1] 4 n, v, κ → (r, g, b, a). (6)
First, the color becomes a function of κ:

c i (κ) = λ i (κ)ĉ i + (1 -λ i (κ)) či λ i (κ) = 0.5 + 0.5 tanh(κs i ) (7)
where ĉi and či are the colors assigned to convex and concave features respectively, λ i : R → [0, 1] is a concave-to-convex transition function, and s i is a user-controlled slope for this transition. Figure 4-left shows the effect obtained by using a pair of slightly different colors. Second, the shading function is warped so that the light direction l i is locally attracted to or repelled from d i . This warping is similar in idea to the Light Warping technique of Vergne et al. [VPB * 09], except that we have an explicit knowledge of where the light direction comes from. This is done by locally altering the parametrization: where T i : R → R is a curvature-dependent translation in parametric space, and t i is a user-controlled translation offset. Note that both c i and ρ i make use of the same concave-to-convex transition function; it ensures that shading color and shape variations are properly correlated, and makes manipulation easier for artists. Figure 4-middle shows the result of using shape features for manipulating shading only, while Figure 4-right shows the two methods in combination.

ρ i (n, v, κ) = α i • (u i (n, v) -T i (κ)) T i (κ) = (2λ i (κ) -1) t i (8) INRIA
Augmented with these controls, our BADF model is able to depict shape features in a variety of ways, ranging from simple changes of color, to applying shading in concavities or convexities only, brightening convexities in otherwise dark regions (or darkening concavities in bright regions), or any other kind of combination.

Directional control

The last two sections described the behavior of each layer, and described a list of shading and shape parameters that we summarize in table 1 Apart from the first two parameters in italics, all other parameters may be considered free variables, and could for instance be key-framed in an off-line animation system. In this paper, we rather focus on real-time stylized illustrations, and therefore we provide another type of control, based on direction. For artistic illustrations, we choose to control appearance based on lighting direction, expressed in screen-space. This choice makes it easy to finely control back-lighting effects in particular. In scientific illustrations, we rather use the direction of the view vector in world space, so that control may be based on the part of the object revealed from a given viewpoint. Other types of directional control might be used within our system, such as the orientation of the sun in outdoor 3D scenes, or objects proximity relations, but we leave their exploration to future work.

Our approach is to let an artist tune a set of BADF parameters for K key (lighting or viewing) directions ω k ; and then interpolate in-between to get a direction-dependent BADF. Figure 5 illustrates this functionality on a simple example. Consider an arbitrary parameter p i for light l i . Its value is now given as a function of ω:

p i (ω) = K k=1 w k (ω)p k i (9)
where p k i stores the parameter value in direction ω k , and the w k are weights ensuring a smooth interpolation and satisfying p i (ω k ) = p k i and k w k = 1. We choose to use a biharmonic interpolation because it satisfies these requirements with C 2 continuity. Each key direction ω k introduces a constraint on the sphere of directions, in the form of a weight vector (w 1 , w 2 , • • • , w n ), with all weights equal to 0 except for w k = 1. We then solve for the biharmonic solution given these constraints using a least-squares minimization approach, as explained in Appendix B. Finally, we constrain the resulting weights to the [0, 1] range with a smooth-step function, and a renormalization step. The resulting values interpolate key direction weights smoothly, and are well-adapted to our needs. However, we note that our shading technique is not dependent upon this choice of interpolation schemes, and other schemes like spherical radial basis functions [START_REF] Tsai | All-frequency precomputed radiance transfer using spherical radial basis functions and clustered tensor approximation[END_REF] would work as well.

Results

We have implemented the layered BADF model described in the previous section mostly on the GPU, with the exception of the biharmonic weight interpolation. Our system runs in real-time (≥ 30fps) on all our examples on a G80 graphics card in 1600 × 1200, which is adequate for BADF editing and creation. With our prototype system, we spent approximately 10 minutes on average for the creation of the BADFs seen in this paper. We also conducted a pilot user study, presented in supplementary material. It points out that despite the drawbacks of our UI, the system is easy to use and fairly intuitive, even for users without extensive experience designing materials. In the following, we present a panel of results obtained with our system for artistic and scientific illustration. A detailed collection of BADFs, including reference artwork, is given in supplementary materials and videos.

INRIA

Artistic illustration

Figure 6 shows a series of results intended to reproduce artistic illustrations. The simplest one is shown at center. It is composed of only two layers: a black fill-in layer, and a white specular layer depicting convex shape features in a quite exaggerated way. Thanks to the warping technique, the resulting shading provides a compelling depiction of shape within a single layer. A pair of moderately more complex BADFs is shown at the top-left. In both cases, two opposing lights are used to depict all the details of the face. They only differ in their intensity and profile between the two figures. The comics-like example at bottom-left of Figure 6 shows a character with a single BADF, lit from three different directions. Shape features are used to align shading with anatomical details and darken colors in concavities. Simpler BADFs, more akin to manga styles, are shown at bottom-right. Here, screenaligned layers create effects such as contour lines and embossing. Glossy effects are easy to obtain as shown at top-right, where we illustrate two BADFs using a same number of layers, but with different colors and light directions. Here again, surface features help convey object shape regardless of chosen illumination. In the bottom-left of Figure 7, we focused on reproducing the brilliant sheen of a beetle, and added directional control to enhance concavities when the specimen is seen from behind. Finally, we have combined three different BADFs in the image at bottom-right, in order to create a scientific illustration in symbolic colors. Here we make use of a pastel BADF for the bones to focus on internal organs, providing a neutral context.

Scientific illustration

Discussion and future work

We designed our BADF model with a focus on interactive applications in mind. Its main idea is to encapsulate the behavior and effects of a few light sources to provide an integrated shading style in real-time, which is controlled at run-time INRIA by manipulating an illumination reference frame. This is adapted to interactive applications such as digital sculpture, scientific illustration, video games or even movie production (especially cartoons), in cases where end-users need not bother with the locations of individual light sources. The choice of directional lights follows this logic as it makes it easy to apply a same BADF to multiple 3D objects. Similarly, directional control permits richer variations without requiring any additional manipulation on the end-user side.

The extension to point light sources might not require many changes to our BADF model though. But the technique is also likely to lose some of its appeal since light positions would need to be adjusted for every new 3D object. In cases where a more advanced control on light positions is needed, a nodal BADF model might be more adapted. In such a scenario, each light source corresponds to a single node organized in a shading network, and is controlled independently. This is the conventional approach in movie production, and we believe our shading functions offer a contribution in this application domain as well. Artists might prefer to use key-framed instead of directional control in this case though.

There is room for improvements to our current BADF model. Shading functions are isotropic, and thus produce simple shading shapes. Although it is sufficient for creating a rich variety of shading styles, specific effects might require additional parameters. As shown in the supplemental materials, we experimented with simple transformations of the spherical parametrization, allowing us to reproduce more complex shading shapes such as in the work of Anjyo et al. [START_REF] Anjyo | Stylized Highlights for Cartoon Rendering and Animation[END_REF][START_REF] Anjyo | Tweakable light and shade for cartoon animation[END_REF] or Pacanowski et al. [START_REF] Pacanowski | Sketch and Paint-based Interface for Highlight Modeling[END_REF]. However, we believe that such parametrized models are limited because they become less intuitive with an increasing number of parameters. In the future, we would thus like to investigate example-based approaches, enabling painting-like interactions in the vein of the work of Okabe et al. [START_REF] Okabe | Illumination Brush: Interactive Design of All-Frequency Lighting[END_REF].

We presented two techniques to depict shape through shading, and combined them to produce shading styles that more resemble hand-made artistic and scientific illustrations. However, for extreme concave-to-convex transitions, both techniques tend to "break" the perception of material. This is explained by the fact that the gradient of shading is modified in a way that destroys basic material cues. A better understanding of material perception [START_REF] Adelson | On Seeing Stuff: The Perception of Materials by Humans and Machines[END_REF] might give rise to more intuitive BADF models in the future. Moreover, we only considered opaque materials and ignored visibility. The extension to transparent and translucent materials or cast shadows might require novel parametrizations which take a 3D object's volume into account. Again, understanding the perceptual cues triggered by such materials [START_REF] Fleming | Perceiving translucent materials[END_REF] is likely to provide for simple and intuitive manipulations.

A Profile function

In this appendix, we give the formula for constructing the profile function introduced in Section 3.2. Without loss of generality, we consider the profile function The cutoff parameter c corresponds to the angle at which the profile function reaches zero. The bias function β(c) gives the bias value that satisfies this condition: η(c) = 0. It is defined as:

β(c) = - cos(c) 1 -cos(c) (10)
The falloff parameter f marks the angle at which the profile function passes through an inflection point, which corresponds to the zero-crossing of its second derivative. The exponent function γ(f, c) gives the (strictly positive) exponent value that satisfies this condition: η ′′ (f ) = 0. The profile's second derivative is:

η ′′ (x) = (β -1)γ cos x(β + (1 -β) cos x) γ-1 + (1 -β) 2 (γ -1)γ(β + (1 -β) cos x) γ-2 sin 2 x
There is only strictly positive zero-crossing for η ′′ (f ):

γ(f, c) = β -β cos(f ) -1 (β -1) sin 2 (f )
Replacing β by the formula for the bias function β(c) gives:

γ(f, c) = 1 -cos(c) cos(f ) sin 2 (f ) (11)

B Biharmonic interpolation

We detail our least-squares technique used to interpolate weights in Section 3.4. Our approximation consists in solving a discretized version of the bi-Laplace equation on a regularly tessellated sphere, which yields a biharmonic interpolation. It is an adaptation of a technique used for 2D animation [START_REF] Baxter | N-way morphing for 2D animation[END_REF] from the planar to spherical domain. For each pair of adjacent triangles, the method minimizes the difference between the Jacobians of their vertex weights:

arg min w a,b∈E(T ) J a (w) -J b (w) 2 ,
where E(T ) is the set of triangle pairs that share an edge and J a (w) is the Jacobian of triangle a. The only difference with the planar approach is that we unfold adjacent triangles onto a common plane. This amounts to performing parallel transport of vectors along geodesics of the discrete sphere.

Figure 1 :

 1 Figure 1: Layered BADF model: each light source is assigned to a layer in a way similar to image editing software. Here we illustrate this organization on a simple BADF. From left to right: a fill layer, a diffuse layer, a specular layer, a back-lighting layer, and the composited result. More complex BADFs are shown in subsequent figures.

Figure 2 :

 2 Figure 2: Shading parametrization: Using a light direction oriented at top left with different parametrizations, a shading function achieves various behaviors on this blobby object. From left to right, we show the color code we used, then diffuse, specular and screen-aligned parametrizations.

Figure 3 :

 3 Figure 3: Shading profile: our generalized Phong shading function is able to create a variety of styles. Left: The profile is controlled by intuitive offset, fall-off and cut-off parameters. The red, green and brown curves correspond to Lambert, Half-Lambert and custom profiles. Right: profile parameters are directly controlled on top of the object surface via isocurve widgets (in blue).

Figure 4 :

 4 Figure 4: Shape depiction: we incorporate shape features into shading via two approaches. Left: we vary color in convexities and concavities as in Mean Curvature Shading. Middle: we modify shading extent in convex and concave regions. Right: applying both enhancement techniques.

Figure 5 :

 5 Figure 5: Directional control: Shading parameters are easily varied per layer to interpolate between a few key directions. Left: the original BADF illuminated from top left. Middle: when lit from behind, back-light and highlights are too bright. Right: we slightly change the diffuse profile, darken the back-light and highlight colors, and extend the back-light toward concavities for this key lighting direction.

Figure 6 :

 6 Figure 6: Artistic illustrations: Our BADF model provides for a vast range of artistic shading styles, such as minimal lighting (center), dual lighting (top-left), comics-like (bottom-left), manga-like (bottom-right), and glossy effects (top-right).

Figure 7

 7 Figure7demonstrates the subtle shading behaviors obtained with BADFs in scientific illustration scenarios. The pastel shading shown at top-left resembles sanguine drawings, and directional control permits to accurately position reddish and whitish shading functions around the hand silhouette when seen from behind. With a careful use of BADF layers, it also becomes possible to mimic existing materials as shown on the lemon and orange models at top-right.In the bottom-left of Figure7, we focused on reproducing the brilliant sheen of a beetle, and added directional control to enhance concavities when the specimen is seen from behind. Finally, we have combined three different BADFs in the image at bottom-right, in order to create a scientific illustration in symbolic colors. Here we make use of a pastel BADF for the bones to focus on internal organs, providing a neutral context.

Figure 7 :

 7 Figure 7: Scientific illustrations: Various kinds of scientific illustrations are simple to produce with our BADF model, including pastel sketches (top-left), mimicking real-world materials (top-right), creating a brilliant sheen in front while enhancing shapes from behind (bottom-left), and making use of symbolic and subdued colors (bottom-right).

  η(x) = |β + (1 -β) cos(x)| + γ that is related to α by α(u) = η • ζ(u -τ ), and where we have omitted the i subscript for clarity.

Table 1 :

 1 . List of layer parameters

	blending function b i	shading type	d i
	light direction	l i	shading offset	τ i
	shading fall-off	f i shading cut-off c i
	shape slope	s i shape offset	t i
	convex color	ĉi concave color	či

RR n°7195


Morgane Bergot
email: morgane.bergot@inria.fr

Gary Cohen
email: gary.cohen@inria.fr

Marc Duruflé
email: marc.durufle@math.u-bordeaux1.fr

HIGHER-ORDER FINITE ELEMENTS FOR HYBRID MESHES USING NEW NODAL PYRAMIDAL ELEMENTS *

Keywords: pyramidal element, higher-order finite element, hybrid mesh, conformal mesh, continuous finite element, discontinuous Galerkin method, error estimates, quadrature formula

We provide a comprehensive study of arbitrarily high-order finite elements defined on pyramids. We propose a new family of high-order nodal pyramidal finite element which can be used in hybrid meshes which include hexahedra, tetrahedra, wedges and pyramids. Finite elements matrices can be evaluated through approximate integration, and we show that the order of convergence of the method is conserved. Numerical results demonstrate the efficiency of hybrid meshes compared to pure tetrahedral meshes or hexahedral meshes obtained by splitting tetrahedra into hexahedra.

Introduction

Highly efficient finite element methods using hexahedral meshes have been developed by Cohen [START_REF] Cohen | Higher-Order Numerical Methods for Transient Wave Equations[END_REF] and his collaborators (Fauqueux [START_REF] Cohen | Mixed Finite Elements with Mass-Lumping for the Transient Wave Equation[END_REF], Pernet and Ferrières [START_REF] Cohen | A Spatial High-Order Hexahedral Discontinuous Galerkin Method to Solve Maxwell Equations in Time Domain[END_REF], [START_REF] Pernet | hp a-priori Error Estimates for a Non-Dissipative Spectral Discontinuous Galerkin Method to Solve the Maxwell Equations in the Time Domain[END_REF], Duruflé [START_REF] Duruflé | Intégration numérique et éléments finis d'ordre élevé appliqués aux équations de Maxwell en régime harmonique[END_REF], [START_REF] Duruflé | Influence of the Gauss and Gauss-Lobatto quadrature rules on the accuracy of a quadrilateral finite element method in the time domain[END_REF]) but currently the only systematic way to generate unstructured hexahedral meshes for a complex geometry is to generate a tetrahedral mesh, and split each tetrahedron into four hexahedra, which introduce needlessly substantial increase in the cost. However, some mesh generators are able to produce hexahedral-dominant meshes that include a minor number of tetrahedra, wedges and pyramids. The aim here is to study finite element methods on hybrid meshes in order to preserve the efficiency of the method developed for hexahedra.

Nodal finite elements are detailed in Hesthaven and Teng [START_REF] Hesthaven | Stable Spectral Methods on Tetrahedral Elements[END_REF] for tetrahedra, and Cohen [START_REF] Cohen | Higher-Order Numerical Methods for Transient Wave Equations[END_REF] for hexahedra. Wedge (or triangular prism) nodal finite elements are constructed as a tensor product between Legendre-Gauss-Lobatto (LGL) points on [0, 1] and electrostatic points on the triangle including LGL points on the edges [START_REF] Hesthaven | Stable Spectral Methods on Tetrahedral Elements[END_REF]. In this work, the main effort is devoted to the construction of pyramidal finite elements, preserving conformity with the other types of elements.

Since obtaining a proper base for nodal pyramidal elements is a tricky point, two approaches have been attempted. A first approach consists in using rational functions in order to obtain nodal shape functions.

First works about nodal pyramidal elements have been made by Bedrosian in [START_REF] Bedrosian | Shape Functions and Integration Formulas for Three-Dimensional Finite Element Analysis[END_REF] where he noticed the impossibility of choosing polynomial shape functions if we want to preserve the conformity with other elements. As a solution, he proposes what he calls "rabbitfunctions" for first-order and second-order approximations. But, the second-order approximation does not include a node at the center of the quadrilateral base, which prohibits the conformity with the second-order hexahedron.

Zgainski et al. [START_REF] Zgainski | A New Family of Finite Elements : The Pyramidal Elements[END_REF] perform numerical experiments with the basis functions given by Bedrosian, and propose a modified second-order set of shape functions by adding a node at the center of the quadrilateral base. However, the central basis function proposed does not satisfy the nodal condition ϕ i (M j ) = δ ij , and the modification does not improve the accuracy, since the finite element space generated by this set of basis functions does not contain P 2 . The same idea is taken back by Graglia et al. [START_REF] Graglia | Higher Order Interpolatory Vector Bases on Pyramidal Elements[END_REF] who achieve to improve the accuracy with their own second-order central basis function.

Chatzi and Preparata [START_REF] Chatzi | Using Pyramids in Mixed Meshes -Point Placement and Basis Functions[END_REF] introduce a generalization of Bedrosian basis functions at any order for nodes regularly distributed on the pyramid. Unfortunately, these basis functions are not consistent for order greater or equal to three since polynomials are not generated by these functions.

The second approach is to split the pyramid into tetrahedra to avoid the use of rational fractions, which have the debatable reputation to make the basis functions hard to manipulate, and instead use polynomial basis functions.

Wieners [START_REF] Wieners | Conforming discretizations on tetrahedrons, pyramids, prisms and hexahedrons[END_REF], Knabner and Summ [START_REF] Knabner | The Invertibility of the Isoparametric Mapping for Pyramidal and Prismatic Finite Elements[END_REF], and Bluck and Walker [START_REF] Bluck | Polynomial Basis Functions on Pyramidal Elements[END_REF] provide a consistent firstorder set of shape functions which ensures the conformity with tetrahedra and hexahedra, by splitting a pyramid into two tetrahedra. Second-order shape functions have been proposed by Wieners, and high order shape functions by Bluck and Walker. However, the finite element space of higher order does not contain the low order finite element space, which leads to a non-consistent method in the case of non-affine pyramids. Moreover, this method requires expensive quadrature on each tetrahedron.

Liu et al. [START_REF] Liu | On symmetric pyramidal Finite Elements[END_REF] propose to symmetrize shape functions of Wieners, but this modification barely improves the accuracy of the method.

An other popular alternative for finite element is the hp approach (Szabó and Babuška [START_REF] Szabó | Finite Element Analysis[END_REF]), e.g. with Šolín et al. [START_REF] Šolín | Higher-order finite elements methods[END_REF] for hexahedra, tetrahedra and wedges. Several papers extend the hp finite element to pyramidal elements.

Warburton [START_REF] Warburton | Spectral/hp Methods on Polymorphic Multi-Domains: Algorithms and Applications[END_REF], Sherwin [START_REF] Sherwin | Hierarchical hp Finite Element in Hybrid Domains[END_REF], Sherwin et al. [START_REF] Sherwin | Spectral/hp Methods For Elliptic Problems on Hybrid Grids[END_REF], and Karniadakis and Sherwin [START_REF] Karniadakis | Spectral/hp element methods for CF -Second Edition[END_REF] provide a tensorial set of basis functions for all types of elements based on the degeneration of a cube. For tetrahedra, hexahedra and wedges, the generated finite element spaces are standards. For pyramids, the proposed generated finite element space provides an optimal convergence for affine pyramids, but not for distorted pyramids for an order greater or equal to two. Moreover, the continuous transition between pyramids and tetrahedra is not achievable for general unstructured meshes. [START_REF] Nigam | Higher-Order Finite Elements on Pyramids[END_REF] propose an original finite element space by deriving pyramidal finite elements from a reference element that is the infinite pyramid. With the finite element space they obtain, the accuracy is preserved but the dimension of this space could be reduced. Demkowicz et al. [START_REF] Demkowicz | Computing with hp-adaptive finite elements[END_REF] and Zaglmayr [START_REF] Zaglmayr | High Order Finite Elements for Electromagnetic Field Computation[END_REF] give the construction of partial-orthogonal basis functions for tetrahedra, hexahedra and wedges, and exploit the use of a degenerated cube for pyramidal elements to get a finite element space that preserve the optimal accuracy, with a smaller dimension than Nigam and Phillips.

Nigam and Phillips

In this paper, the reference element is the symmetric unit pyramid (Fig. 1.1). Let P r be the polynomial space of degree r, we claim that if we choose the following finite element space

Pr = P r (x, ŷ, ẑ) ⊕ 0≤k≤r-1 xŷ 1 -ẑ r-k P k (x, ŷ),
we are able to produce optimal error estimates in

H 1 norm u -π r u 1,K ≤ Ch r u r+1,K
with the notations detailed in Section 4, for continuous finite elements.

In order to evaluate integrals, we propose in Section 3 to use the same technique as Bedrosian, detailed by Hammer, Marlowe and Stroud in [START_REF] Hammer | Numerical Integration Over Simplexes and Cones[END_REF], adapted to the pyramid and which does not deteriorate the accuracy, as it will be proved in Section 4. An extension of this work is proposed for discontinuous Galerkin formulation with the same finite element space Pr .

To validate this new pyramidal finite element, a dispersion analysis is carried out in the case of periodic meshes. We have observed an optimal dispersion error in O(h 2r ) as obtained for other element shapes. Furthermore these elements have been tested for the Helmholtz equation with the continuous Galerkin formulation, and for the unsteady wave equation with the discontinuous Galerkin method. The numerical experiments show that they are much more efficient than purely tetrahedral elements, or hexahedral meshes generated by splitting each tetrahedron into four hexahedra.

The outline of our paper is as follows:

In Section 1, following the classical notations of Ciarlet [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF], we define two pyramidal finite elements of order r, ( K, Pr , Σ) on the reference element, and (K, P r , Σ) for any pyramid in the mesh ;

A comparison to existing hp finite element spaces is given in Section 2, along with possible improvements of these spaces, in propositions 2.2 and 2.3 ;

The quadrature formula used to get exact integrals, whenever it is possible, for the basis functions constructed from the finite element space Pr are presented in Section 3, thanks to a change of variable from the unit cube ; Section 4 is devoted to the error analysis which is performed in a classical way ;

The case of a discontinuous Galerkin formulation is briefly treated in Section 5 ; Section 6 is devoted to numerical results: a dispersion analysis is performed on the wave equation in section 6.1, the stability condition (CFL) is computed on a periodic infinite mesh in section 6.2, and numerical experiments are performed in section 6.4 along with explanations about storage.

1. Arbitrary High-Order Pyramidal Element 1.1. Pyramidal Element Definition 1.1 A pyramid K(x, y, z) is the image of the reference pyramid K(x, ŷ, ẑ) taken as the unit symmetrical pyramid, centered at the origin by the transformation F given by Bedrosian [START_REF] Bedrosian | Shape Functions and Integration Formulas for Three-Dimensional Finite Element Analysis[END_REF] using rational fractions, as shown in Fig. 1.1

F = 1≤i≤5 S i φ1 i , (1)
where S i = (x i , y i , z i ) are the vertices of the pyramid K and φ1 i are the following mapping functions The case of a non-invertible transformation may occur when considering a degenerated element, e.g. when the five vertices are co-planar, but the characterization of pyramids for which F is invertible remains an open question, as for hexahedra (Duruflé et al. [START_REF] Duruflé | Influence of the Gauss and Gauss-Lobatto quadrature rules on the accuracy of a quadrilateral finite element method in the time domain[END_REF]). In the sequel, we assume that F is always invertible.

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ φ1 1 = 1 4 1 -x -ŷ -ẑ + xŷ 1 -ẑ φ1 2 = 1 4 1 + x -ŷ -ẑ - xŷ 1 -ẑ φ1 3 = 1 4 1 + x + ŷ -ẑ + xŷ 1 -ẑ φ1 4 = 1 4 1 -x + ŷ -ẑ - xŷ 1 -ẑ φ1 5 = ẑ when F is invertible. K Ŝ1 = (-1, -1, 0) Ŝ2 = (1, -1, 0) Ŝ5 = (0, 0, 1) Ŝ4 = (-1, 1, 0) Ŝ3 = (1, 1, 0) F ẑ x ŷ K S 5 S 1 S 2 S 3 S 4
The transformation F can be explicitly written as

4F = (S 1 + S 2 + S 3 + S 4 ) + x (-S 1 + S 2 + S 3 -S 4 ) + ŷ (-S 1 -S 2 + S 3 + S 4 ) + ẑ (4S 5 -S 1 -S 2 -S 3 -S 4 ) + xŷ 1 - ẑ (S 1 + S 3 -S 2 -S 4 ).
We notice that F is affine when

S 1 + S 3 = S 2 + S 4 ,
i.e. when the base of the pyramid is a parallelogram. Furthermore, F ensures the conformity with tetrahedra and hexahedra as the shape functions becomes a two-dimensional triangular or quadrilateral shape function, since adjacent tetrahedra, wedge and hexahedra have the same property. That would not be the case if F had been chosen to be polynomial.

Remark 1.3

The shape function of Bedrosian can be found by defining the transformation T from the unit cube Q to the reference pyramid K

T ∶ ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ x = (1 -z)(2x -1) ŷ = (1 -z)(2ỹ -1) ẑ = z.
(

) 2
For a basis function of the hexahedron ϕ(x, ỹ, z)

= (1 -x)(1 -ỹ)(1 -z), the transformation T gives indeed ϕ ○ T -1 (x, ŷ, ẑ) = 1 4 (1 -x -ẑ)(1 -ŷ -ẑ) 1 - ẑ = φ1 1 (x, ŷ, ẑ).
Similarly, we find the other functions of Bedrosian.

A Pyramidal Finite Element Space of Order r

We place ourselves in the most restrictive case, that is continuous finite elements. The finite element space V h on an open set Ω of R 3 is given by

V h = {u ∈ H 1 (Ω) u K ∈ P F r (K)},
where P F r is the real space of order r for an element K of the mesh defined by

P F r (K) = {u u ○ F ∈ Pr ( K)},
The finite element space Pr of order r on K is

P r (x, ŷ, ẑ) = xi ŷj ẑk , i + j + k ≤ r when K is a tetrahedron ; Q r (x, ŷ, ẑ) = xi ŷj ẑk , i, j, k ≤ r when K is a hexahedron ; P r (x, ŷ) ⊗ P r (ẑ) = xi ŷj ẑk , i + j ≤ r, k ≤ r when K is a wedge ;
and defined by identity (5) when K is a pyramid.

To use the Bramble-Hilbert's lemma and get optimal error estimates, the real space P F r for a pyramidal element K of the mesh must be such that P r (x, y, z) ⊂ P F r .

Theorem 1.4 When F is affine, the minimal space Pr of order r such that we have the inclusion (3) is

Pr = P r (x, ŷ, ẑ). ( 4
)
When F is not affine, the minimal space Pr of order r such that we have the inclusion (3) is

Pr = P r (x, ŷ, ẑ) ⊕ 0≤k≤r-1 xŷ 1 -ẑ r-k P k (x, ŷ). ( 5
)
Proof. When F ∈ P 1 , it is easy to see that

Pr ( K) = P r ( K) ⇐⇒ P F r (K) = P r (K),
which means that taking Pr = P r when the base of the pyramid is a parallelogram is necessary and sufficient to satisfy [START_REF] Bluck | Polynomial Basis Functions on Pyramidal Elements[END_REF].

For any base of the pyramid, we take f ∈ P r , i.e.

f = 0 ≤ i, j, k ≤ r, i + j + k ≤ r x i y j z k .
We study the case f = x n , n ≤ r. Using the transformation F , f can be written as

4 n f = 4 n x n = x 1 (1 -x -ŷ -ẑ) + x 2 (1 + x -ŷ -ẑ) + x 3 (1 + x + ŷ -ẑ) + x 4 (1 -x + ŷ -ẑ) + 4 ẑx 5 + xŷ 1 - ẑ (x 1 + x 3 -x 2 -x 4 ) n .
As the part

x 1 (1 -x -ŷ -ẑ) + x 2 (1 + x -ŷ -ẑ) + x 3 (1 + x + ŷ -ẑ) + x 4 (1 -x + ŷ -ẑ) + 4 ẑx 5 n
is in P n (x, ŷ, ẑ), it remains to handle the terms

(a + bx + cŷ + dẑ) k ( xŷ 1 - ẑ ) n-k k ≤ n -1.
Developing the first factor, we get terms of the form

ẑp (α + β x + γ ŷ) k-p ( xŷ 1 - ẑ ) n-k .
If p = 0, the factor belongs to xŷ 1ẑ n-k P k (x, ŷ). Otherwise, we decrease the power of ẑ, by

writing ẑ = 1 -ẑ + ẑ ẑp-1 (α + β x + γ ŷ) k-p ( xŷ 1 - ẑ ) n-k + ẑp-1 (α + β x + γ ŷ) k-p x ŷ ( xŷ 1 - ẑ ) n-k-1 .
Iterating this method, we erase all the powers of ẑ to obtain a term of higher degree

(α + β x + γ ŷ) k-p (x ŷ) p ( xŷ 1 - ẑ ) n-k-p .
However, when k + p ≥ n, the iterative procedure stops as we obtained the polynomial ẑp+k-n (α

+ β x + γ ŷ) k-p (xy) n-k , (6)
and the degree of this polynomial is equal to k ≤ r -1. Since k + p ≤ r -1,

(α + β x + γ ŷ) k-p (x ŷ) p ∈ P k+p (x, ŷ),
and the term is finally in

P m (x, ŷ) ( xŷ 1 - ẑ ) n-m with m = k + p, m ≤ n -1.
We let the reader convince himself that other cases can be treated similarly. At this point, we proved that it is sufficient to take Pr as specified by Theorem 1.4 to obtain the inclusion [START_REF] Bluck | Polynomial Basis Functions on Pyramidal Elements[END_REF]. ◻ Corollary 1.5

dim Pr = 1 6 (r + 1)(r + 2)(2r + 3).
Proof. We classically have

dim P r (x, ŷ, ẑ) = 1 6 (r + 1)(r + 2)(r + 3)
and, using the direct sums property,

dim 0≤k≤r-1 xŷ 1 -ẑ r-k P k (x, ŷ) = 0≤k≤r-1 dim P k (x, ŷ) = 0≤k≤r-1 (k + 1)(k + 2) 2 = r(r + 1)(r + 2) 6 , that is dim Pr (x, ŷ, ẑ) = 1 6 (r + 1) (r + 2) (r + 3) + 1 6 r (r + 1) (r + 2)
which provides the claimed result.

◻ Proposition 1.6 Pr x=1-ẑ or x=ẑ-1 = P r (ŷ, ẑ). Pr ŷ=1-ẑ or ŷ=ẑ-1 = P r (x, ẑ). Pr ẑ=0 = Q r (x, ŷ). (7)
Proof. Any function p ∈ Pr can be written as

p(x, ŷ, ẑ) = p r (x, ŷ, ẑ) + 0≤k≤r-1 p k (x, ŷ) xŷ 1 -ẑ r-k , with p r ∈ P r (x, ŷ, ẑ) and p k ∈ P k (x, ŷ).
On a triangular face, we replace x by ±(1ẑ) or ŷ by ±(1ẑ), according to the considered face. For example for the face x = (1ẑ), as p r ((1ẑ), ŷ, ẑ) obviously belongs to P r (ŷ, ẑ), rational parts become

p k (x, ŷ) xŷ 1 -ẑ r-k = p k ((1 -ẑ), ŷ) y r-k , 0 ≤ k ≤ r -1.
As p k ((1ẑ), ŷ) ∈ P k (ŷ, ẑ), we have p k ((1ẑ), ŷ)y r-k ∈ P r (ŷ, ẑ), and finally p ∈ P r (ŷ, ẑ). The same simplification can be done for the other faces.

On the quadrangular base, we replace ẑ by 0: p r (x, ŷ, 0) is obviously in Q r (x, ŷ), and the rational parts become

p k (x, ŷ) xŷ 1 -ẑ r-k = p k (x, ŷ) x r-k y r-k , 0 ≤ k ≤ r -1. As p k (x, ŷ) ∈ P k (x, ŷ), we have p k (x, ŷ)x r-k y r-k ∈ Q r (x, ŷ), and finally p ∈ Q r (x, ŷ).
The proposition is finally proved using a dimension argument. ◻ Proposition 1.7

P F r (K) ⊂ H 1 (K). Proof. For p ∈ P F r , p ∈ C ∞ ( K S 5
) as a the rational fraction when its pole is not in the domain. The continuity in S 5 is proved by considering four pseudo-faces

F i ε , 0 ≤ i ≤ 4, 0 ≤ ε ≤ 1 mapping a quarter Q i of the pyramid. We consider the face F 2 ε represented in Fig. 1.2 in red, Q 2 being represented in blue. F 2 ε is such that ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ x = (1 -z)(1 -ε) -(1 -z)(1 -ε) ≤ y ≤ (1 -z)(1 -ε) 0 ≤ z ≤ 1,
and we have

∀M = (x, y, z) ∈ Q 2 , ∃ε ∈ [0, 1], M ∈ F 2 ε . x ŷ ẑ ε Fig. 1.2. Pseudo face F 2 ε
All cases amount to study p = xy 1z which is the difficult case.

p F 2 ε = (1 -z)(1 -ε)y 1 -z = (1 -ε)y z→1 → 0.
The three other cases are similar by symmetry, and finally p ∈ C 0 ( K).

As K is bounded, we get p ∈ L 2 (K). As for p, we consider ∇p on a quarter of pyramid, for example Q 2 and we consider an ε such that M ∈ F 2 ε , and

-(1 -ε) 2 ≤ ∂p ∂z F 2 ε ≤ (1 -ε) 2 ,
that is ∂ z p is bounded. The same technique applied for ∂ x p and ∂ y p leads to conclude that ∇p in bounded inside K. As K is also bounded, we have ∇p ∈ L 2 (K). ◻ Proposition 1.6 ensures to a function u ∈ V h to be continuous across the interface between elements, whatever the type of the elements adjacent to the face, and therefore to belong to H 1 (Ω) due to proposition 1.7.

Proposition 1.8 The optimal finite element space of order r on the unit cube Q is

C r = Pr ○ T = 0≤k≤r Q k (x, ỹ)(1 -z) k .
Proof. Using the transformation T, the polynomial part of Pr becomes {x m ŷn ẑp , m

+ n + p ≤ r} ○ T = {(1 -2x) m (1 -2ỹ) n (1 -z) m+n zp , m + n + p ≤ r} = {x m ỹn zp (1 -z) m+n , m + n + p ≤ r} ⊂ C r ,
whereas the fractional part of Pr becomes

xi ŷj xŷ 1 -ẑ r-p , i + j ≤ p ≤ r -1 ○ T = (2x -1) r-p+i (2ỹ -1) r-p+j (1 -z) r-p+i+j , i + j ≤ p ≤ r -1 = xr-k+i ỹr-k+j (1 -z) r-k+i+j , 0 ≤ i + j ≤ k ≤ r -1 ⊂ C r , that is Pr ⊂ C r . We now notice that dim C r = 0≤k≤r (k + 1) 2 = 1 6 (r + 1)(r + 2)(2r + 3) = dim Pr = dim Pr ,
which proves the proposition. ◻

Location of the Degrees of Freedom

We wish to link continuously pyramidal elements with other elements of the mesh hexahedra with Legendre-Gauss-Lobatto (LGL) points ; tetrahedra with Hesthaven "electrostatic points" constructed with LGL points on edges (Hesthaven and Teng [START_REF] Hesthaven | Stable Spectral Methods on Tetrahedral Elements[END_REF]); wedges obtained by a tensorial product of a face of a tetrahedron of Hesthaven, that is a triangle of Hesthaven (Hesthaven [19]), with an edge with LGL points.

We place the degrees of freedom on LGL points on the quadrangular base of the pyramid, and on Hesthaven points on each triangular face. The number of degrees of freedom n f on the faces is then

n f = 3r 2 + 2.
We add n i degrees of freedom inside the pyramid

n i = 1 6 (r -1)(r -2)(2r -3) = 1≤k≤r-2 k 2 .
and we place them on (r -2) parallel planes of k 2 degrees of freedom, as shown in the Fig. 1.3. The total number of degrees of freedom is

n r = n i + n f = 1 6 (r + 1)(r + 2)(2r + 3).
which is precisely the dimension of Pr .

Degrees of freedom can then be placed systematically on the pyramid, at any order. Each category of point is represented by a color in the Fig. 1.4 for the pyramidal elements of order two to four.

Basis Functions

The basis functions on the reference pyramid K are obtained by inverting a Vandermonde system as follow.

Let ( Mi ) 1≤i≤nr the locations of the interpolation points on the pyramid, and ( ψi ) 1≤i≤nr a base of Pr , Definition 1.9 The Vandermonde matrix V DM ∈ M nr (R) is defined by

V DM i,j = ψi ( Mj ), 1 ≤ i, j ≤ n r , (8)
and the basis function φi linked to the interpolation point Mi is then defined as

φi = 1≤j≤nr (V DM -1 ) i,j ψj . ( 9
)
Proposition 1.10 The following set of basis functions is an orthogonal base of Pr

P 0,0 i x 1 -ẑ P 0,0 j ŷ 1 -ẑ (1 -ẑ) max(i,j) P 2max(i,j)+2,0 k (2ẑ -1), 0 ≤ i, j ≤ r, 0 ≤ k ≤ r -max(i, j) ,
where P i,j m (x) denotes the Jacobi polynomial of order m, orthogonal for the weight (1-x) i (1+x) j .

Proof. We note ψi,j,k (x, ŷ, ẑ) = P 0,0

i x 1 -ẑ P 0,0 j ŷ 1 -ẑ (1 -ẑ) max(i,j) P 2max(i,j)+2,0 k (2ẑ -1).
We first prove that the family is orthogonal by using the transformation (2) on Qr

K ψi,j,k (x, ŷ, ẑ) ψi ′ ,j ′ ,k ′ (x, ŷ, ẑ) dx dŷ dẑ = 1 0 P 0,0 i (2x -1)P 0,0 i ′ (2x -1) dx = C ii ′ δ i i ′ 1 0 P 0,0 j (2ỹ -1)P 0,0 j ′ (2ỹ -1) dỹ = C jj ′ δ j j ′ 4 1 0 (1 -z) max(i,j)+max(i ′ ,j ′ )+2 P 2max(i,j)+2,0 k (2z -1)P 2max(i ′ ,j ′ )+2,0 k ′ (2z -1) dz, with 0 ≤ i, j ≤ r, 0 ≤ k ≤ r -max(i, j
), and when i = i ′ and j = j

′ 1 0 (1 -z) 2max(i,j)+2 P 2max(i,j)+2,0 k (2z -1)P 2max(i,j)+2,0 k ′ (2z -1) dz = C kk ′ δ k k ′ .
We also have ψi,j,k

(x, ŷ, ẑ), 0 ≤ i, j ≤ r, k ≤ r -max(i, j) ○ T -1 = P 0,0 i (2x -1) P 0,0 j (2ỹ -1) (1 -z) max(i,j) P 2max(i,j)+2,0 k (2z -1), 0 ≤ i, j ≤ r, k ≤ r -max(i, j) ⊂ C r ,
that is, with an argument of dimension,

Span ψi,j,k (x, ŷ, ẑ), 0 ≤ i, j ≤ r, k ≤ r -max(i, j) ○ T -1 = C r ,
which proves the proposition. ◻

We compare the condition number of the Vandermonde matrix for monomial and orthogonal bases of Pr in the case of tetrahedral and pyramidal elements on Fig. 1.5. We notice that the condition number of the Vandermonde matrix is increasing faster for tetrahedra than for pyramids when using monomial base, whereas we observe the opposite for orthogonal base. Besides, the use of the orthogonal set of basis functions highly improves the condition number of the VDM matrix.

Remark 1.11 The characterization of the invertibility of the Vandermonde matrix is an open question, but we observed that the VDM matrix is invertible with our choice of position for the degrees of freedom, the element is therefore unisolvent.

Comparison with Existing Methods

The nodal basis functions we propose are the same as Bedrosian [START_REF] Bedrosian | Shape Functions and Integration Formulas for Three-Dimensional Finite Element Analysis[END_REF], Zgainski et al. [START_REF] Zgainski | A New Family of Finite Elements : The Pyramidal Elements[END_REF] and Chatzi and Preparata [START_REF] Chatzi | Using Pyramids in Mixed Meshes -Point Placement and Basis Functions[END_REF] for order one. They are identical to those of Graglia et al. [START_REF] Graglia | Higher Order Interpolatory Vector Bases on Pyramidal Elements[END_REF] for order two, and they are new for order greater or equal to three.

The finite element space C r of order r on the unit cube Q defined by in proposition 1.8 is the same as the one proposed by Zaglmayr, cited in [START_REF] Demkowicz | Computing with hp-adaptive finite elements[END_REF],

C 0 r (x, ỹ, z) = (1 -z) 2 x(1 -x) ỹ(1 -ỹ) z Cr-3 .
Proof. The basis functions obviously vanish on the boundary of Q and belongs to C r . The

dimension of the space is dim C r-3 = 1 6 (r -1)(r -2)(2r -3) = n i ,

which proves the proposition. ◻

We write transformation T from the infinite pyramid Q to the unit cube Q.

T ∶ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ x = x ỹ = ȳ z = z 1 + z . ( 10
)
Proposition 2.2 The finite element space U r proposed by Nigam and Phillips [START_REF] Nigam | Higher-Order Finite Elements on Pyramids[END_REF] on the infinite pyramid

Q satisfies U r ⊃ C r ○ T ,
and contains more degrees of freedom than C r since

dim U r = 1 + 3k + k 3 > dim C r .
The subspace U 0 r of U r whose trace is null on the boundary of the element is equal to

U 0 r (x, ȳ, z) = x(1 -x)ȳ(1 -ȳ)z (1 + z) r u(x, ȳ, z), u ∈ Q r-2 (x, ȳ, z) ,
and if we replace U r 0 by C 0 r ○ T , we get the optimal space Ūr = C r ○ T .

Proof. Using the transformation (10), we detail the following basis functions (the others can be treated similarly by symmetry)

For the vertex :

(

1 -x)(1 -ȳ) (1 + z) r ○ T -1 = (1 -x)(1 -ỹ)(1 -z) r ∈ C r .
For the apex : zr

(1 + z) k ○ T -1 = zr ∈ C r .
For the representative vertical edge :

(

1 -x)(1 -ȳ)z a (1 + z) r , 1 ≤ a ≤ r -2 ○ T -1 = {(1 -x)(1 -ỹ)(1 -z) r-a za , 1 ≤ a ≤ r -2} ⊂ C r .
For the representative base edge :

(1 -x)(1 -ȳ)x a (1 + z) r , 1 ≤ a ≤ r -2 ○ T -1 = {(1 -x)(1 -ỹ)x a (1 -z) r , 1 ≤ a ≤ r -2} ⊂ C r .
For the representative triangular face :

(1 -x)(1 -ȳ)x a zb (1 + z) r , a, b ≥ 0, a + b ≤ r -1 ○ T -1 = (1 -x)(1 -ỹ)x a (1 -z) r-b zb , 0 ≤ a + b ≤ r -1 ⊂ C r .
For the base face :

(

1 -x)(1 -ȳ)x a ȳb (1 + z) r , 1 ≤ a, b ≤ r -1 ○ T -1 = (1 -x)(1 -ỹ)x a ỹb (1 -z) r , 1 ≤ a, b ≤ r -1 ⊂ C r .
For the interior :

x(1 -x)ȳ(1 -ȳ)z (1 + z) r u(x, ȳ, z), u ∈ Q r-2 (x, ȳ, z) ○ T -1 = xi+1 (1 -x)ỹ j+1 (1 -ỹ)z k+1 (1 -z) r-k , 0 ≤ i, j, k ≤ r -2 ⊃ C 0
r . The subspace of U r whose trace is null on the boundary of the element is equal to 3 . Since there are n f = 3r 2 + 2 basis functions associated with the boundary, we have

U 0 r = x(1 -x)ȳ(1 -ȳ)z (1 + z) r u(x, ȳ, z), u ∈ Q r-2 (x, ȳ, z) whose dimension is dim U 0 r = dim Q r-2 = (r -1)
dim U r = 3r 2 + 2 + (r -1) 3 = 1 + 3r + r 3 > dim C r .
If we replace U 0 r by C 0 r ○ T , the new finite element space Ūr satisfies dim Ūr = dim C r , and Ūr ⊃ C r ○ T , that is we have the equality of these two spaces. ◻

We write transformation T from the cube

[-1, 1] 3 to the unit cube Q T ∶ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ x = 1 + a 2 ỹ = 1 + b 2 z = 1 + c 2 . ( 11
)
Proposition 2.3 The finite element space W r of order r introduced by Warburton [START_REF] Warburton | Spectral/hp Methods on Polymorphic Multi-Domains: Algorithms and Applications[END_REF] on the cube [-1, 1] 3 is not optimal. The subspace of W r whose trace is null on the boundary of the element is equal to

W 0 r ○ T -1 = x(1 -x)ỹ(1 -ỹ)z(1 -z) 2 u(x, ỹ, z), u ∈ P r-3 (x, ỹ, z) .
If we replace W r 0 by C 0 r ○ T , and the basis functions linked to the base face by the following set of functions

1 -a 2 1 + a 2 1 -b 2 1 + b 2 1 -c 2 max(i,j)+1 P 1,1 i-1 (a)P 1,1 j-1 (b), 1 ≤ i, j ≤ r -1 ,
we get the optimal space W r = C r ○ T .

Proof. Using the transformation [START_REF] Cohen | A Spatial High-Order Hexahedral Discontinuous Galerkin Method to Solve Maxwell Equations in Time Domain[END_REF], we detail the following basis functions (the others can be treated similarly by symmetry)

For the vertex : 1 -a 2 1 -b 2 1 -c 2 ○ T -1 = (1 -x)(1 -ỹ)(1 -z) ∈ C r .
For the apex :

1 + c 2 ○ T -1 = z ∈ C r .
For the vertical edge :

1 -a 2 1 -b 2 1 -c 2 1 + c 2 P 1,1 i-1 (c), 1 ≤ i ≤ r -1 ○ T -1 = (1 -x)(1 -ỹ)(1 -z)z P 1,1 i-1 (2z -1), 1 ≤ i ≤ r -1 ⊂ C r .
For the base edge :

1 -a 2 1 + a 2 1 -b 2 1 -c 2 i+1 P 1,1 i-1 (a), 1 ≤ i ≤ r -1 ○ T -1 = x(1 -x)(1 -ỹ)(1 -z) i+1 P 1,1 i-1 (2x -1), 1 ≤ i ≤ r -1 ⊂ C r .
For the base face :

1 -a 2 1 + a 2 1 -b 2 1 + b 2 1 -c 2 i+j+1 P 1,1 i-1 (a)P 1,1 j-1 (b), 1 ≤ i, j ≤ r -1 ○ T -1 = x(1 -x)ỹ(1 -ỹ)(1 -z) i+j+1 P 1,1 i-1 (2x -1)P 1,1 j-1 (2ỹ -1), 1 ≤ i, j ≤ r -1 ⊄ C r .
For the triangular face :

1 -a 2 1 + a 2 1 -b 2 1 -c 2 i+1 1 + c 2 P 1,1 i-1 (a)P 2i+1,1 j-1 (c), i + j ≤ r -1, i, j ≥ 1 ○ T -1 = (1 -x)x(1 -ỹ)(1 -z) i+1 z P 1,1 i-1 (2x -1)P 2i+1,1 j-1 (2z -1), i + j ≤ r -1, i, j ≥ 1 ⊂ C r .
For the interior :

1 -a 2 1 + a 2 1 -b 2 1 + b 2 1 -c 2 i+j+1 1 + c 2 P 1,1 i-1 (a)P 1,1 j-1 (b)P 2i+2j+1,1 k-1 (c), i + j + k ≤ r -1, i, j, k ≥ 1 ○ T -1 = x(1 -x)ỹ(1 -ỹ)z(1 -z) i+j+1 P 1,1 i-1 (2x -1)P 1,1 j-1 (2ỹ -1)P 2i+2j+1,1 k-1 (2z -1), i + j + k ≤ r -1, i, j, k ≥ 1 ⊂ C r .
The subspace W 0 r of W r whose trace is null on the boundary of the element is equal to

W 0 r ○ T -1 = x(1 -x)ỹ(1 -ỹ)z(1 -z) 2 u(x, ỹ, z), u ∈ P r-3 (x, ỹ, z) , whose dimension is dim W 0 r = dim P r-3 = (r -2)(r -1)r 6 .
Since there are 3r 2 + 2 basis functions associated with the boundary, we have

dim W r = (r -2)(r -1)r 6 + 3r 2 + 2 = (r + 1)(r + 2)(r + 3) 6 + r 2 < dim C r .
If we replace the proposed set of basis functions for the base face by the following one

1 -a 2 1 + a 2 1 -b 2 1 + b 2 1 -c 2 max(i,j)+1 P 1,1 i-1 (a)P 1,1 j-1 (b), 1 ≤ i, j ≤ r -1 ○ T -1 = x(1 -x)ỹ(1 -ỹ)(1 -z) max(i,j)+1 P 1,1 i-1 (2x -1)P 1,1 j-1 (2ỹ -1) 1 ≤ i, j ≤ r -1 ⊂ C r ,
and W 0 r by C 0 r ○ T , the new finite element space W r will satisfy

W r ⊂ C r ○ T and dim W r = dim C r ,
that is we have the equality of the two spaces. ◻ Remark 2.4 As W 1 = W 1 but W r ⊅ W 2 , using W r as a finite element space for pyramidal elements ensures not more than a first-order convergence in H 1 -norm.

Numerical study of the dispersion error has been conducted on periodic meshes containing non-affine pyramids in order to check these theoretical results (see Fig. 6.6 in section 6.3).

Quadrature Formula

To evaluate integrals, we use a quadrature rule defined over the reference pyramid K. A simple rule consists in taking Gauss points over the unit cube Q of coordinates (x, ỹ, z), and compute their image on the reference pyramid K of coordinates (x, ŷ, ẑ), via the change of variable T defined by equation [START_REF] Bedrosian | Shape Functions and Integration Formulas for Three-Dimensional Finite Element Analysis[END_REF], which is a diffeomorphism from the open Q to the open K.

For any function f , we denote f (x, ỹ, z) = f (x, ŷ, ẑ), and the change of variable provides

K f (x, ŷ, ẑ) dxdŷdẑ = Q 4 f (x, ỹ, z) (1 -z) 2 dxdỹdz. ( 12
)
Definition 3.1 Let M be the mass matrix for the pyramid K, defined by

M i,j = K ϕ i ϕ j dxdydz = K DF φi φj dxdŷdẑ (13)
and K the stiffness matrix such that

K i,j = K ∇ϕ i ⋅ ∇ϕ j dxdydz = K DF DF -1 DF * -1 ∇ φi ⋅ ∇ φj dxdŷdẑ. ( 14
)
Definition 3.2 We define the polynomial space

Q m,n,p = x i y j z k , 0 ≤ i ≤ m, 0 ≤ j ≤ n, 0 ≤ k ≤ p . Lemma 3.3 ∀i ∈ 1, n r , φi ∈ Q r (x, ỹ, z).
Proof. Using proposition 1.8, φi (x, ỹ, z) ∈ C r (x, ỹ, z), and we obviously have C r ⊂ Q r , which proves the lemma.

◻ Lemma 3.4 ∀i ∈ 1, n r , ∇ φi ∈ Q r-1,r,r-1 × Q r,r-1,r-1 × Q r,r,r-1 (x, ỹ, z).
Proof. We decompose φi (x, ŷ, ẑ) into the monomial base ψj (x, ŷ, ẑ) of Pr and we treat the different cases.

We first consider the derivative in x, the derivative in y being treated similarly by symmetry: either ψj (x, ŷ, ẑ) ∈ P r (x, ŷ, ẑ),

∂ ψj ∂ x (x, ŷ, ẑ) = xm-1 ŷn ẑp = (2x -1) m-1 (2ỹ -1) n (1 -z) m+n-1 zp , m + n + p ≤ r, or ∂ ψj ∂ x (x, ŷ, ẑ) = xr-p+i-1 ŷr-p+j (1 -ẑ) r-p = (2x -1) r-p+i-1 (2ỹ -1) r-p+j (1 -z) r-p+i+j-1 , i + j ≤ p ≤ r -1, that is ∂ ψj ∂ x (x, ỹ, z) ∈ Q r-1,
r,r-1 (x, ỹ, z) in both cases. Similarly, for the derivative in z, either

∂ ψj ∂ ẑ (x, ŷ, ẑ) = xm ŷn ẑp-1 = (2x -1) m (2ỹ -1) n zp-1 (1 -z) m+n , m + n + p ≤ r or ∂ ψj ∂ ẑ (x, ŷ, ẑ) = xr-p+i ŷr-p+j (1 -ẑ) r-p+1 = (2x -1) r-p+i (2ỹ -1) r-p+j (1 -z) r-p+i+j-1 , i + j ≤ p ≤ r -1, that is ∂ ψj ∂ ẑ (x, ỹ, z) ∈ Q r,r,r-1 (x, ỹ, z) in both cases. ◻ Lemma 3.5 DF = ∂ F ∂ x , ∂ F ∂ ŷ , ∂ F ∂ ẑ ∈ Q 3 0,1,0 × Q 3 1,0,0 × Q 3 1,1,0 (x, ỹ, z),
and

DF ∈ Q 1,1,0 (x, ỹ, z).
Proof. The derivatives of F can be written as

∂F ∂ x = 1 4 (-S 1 + S 2 + S 3 -S 4 ) + 1 4 (S 1 -S 2 + S 3 -S 4 ) ŷ 1 - ẑ = 1 4 (-S 1 + S 2 + S 3 -S 4 ) + 1 4 (S 1 -S 2 + S 3 -S 4 )(2ỹ -1) = A 1 + Cỹ, ∂F ∂ ŷ = 1 4 (-S 1 -S 2 + S 3 + S 4 ) + 1 4 (S 1 -S 2 + S 3 -S 4 ) x 1 - ẑ = 1 4 (-S 1 -S 2 + S 3 + S 4 ) + 1 4 (S 1 -S 2 + S 3 -S 4 )(2x -1) = A 2 + Cx, ∂F ∂ ẑ = 1 4 (4S 5 -S 1 -S 2 -S 3 -S 4 ) + 4 (S 1 -S 2 + S 3 -S 4 ) xŷ (1 -ẑ) 2 = 1 4 (4S 5 -S 1 -S 2 -S 3 -S 4 ) + 1 4 (S 1 -S 2 + S 3 -S 4 )(2x -1)(2ỹ -1) = A 3 + 2Cx ỹ.
The determinant is then

DF (x, ỹ, z) = det(A 1 + Cỹ, A 2 + Cx, A 3 + 2C x ỹ).
We can develop the expression in

DF (x, ỹ, z) = det(A 1 , A 2 , A 3 ) + x det(A 1 , C, A 3 ) + ỹ det(C, A 2 , A 3 ) + 2 x ỹ det(A 1 , A 2 , C),
which proves the lemma. ◻ Proposition 3.6 When F is affine, the quadrature formula must be exact for polynomials of (1 -z) 2 Q 2r for the mass matrix and (1 -z) 2 Q 2r,2r,2r-2 for the stiffness matrix, such that these matrices are exactly integrated. When F is not affine, for the mass matrix to be exactly integrated, the quadrature formula must be exact for polynomials of (1 -z) 2 Q 2r+1,2r+1,2r .

Proof. In the affine case, we can factorize the mass matrix by coefficient DF which is constant. For the mass matrix, the lemma 3.3 provides the term φi φj (1-z) 2 to be in (1-z) 2 Q 2r .

For the stiffness matrix, thanks to lemma 3.4, the term

∂ φi ∂x ∂ φj ∂x (1-z) 2 is in (1-z) 2 Q 2r-2,2r,2r-2 , the term ∂ φi ∂ỹ ∂ φj ∂ỹ (1 -z) 2 is in (1 -z) 2 Q 2r,2r-2,2r-2 and the term ∂ φi ∂z ∂ φj ∂z (1 -z) 2 is in (1 - z) 2 Q 2r,2r,2r-2 , that is the final term is in (1 -z) 2 Q 2r,2r,2r-2 .
In the non-affine case, using the result of the affine case and lemma 3.5, we can easily conclude that we need a quadrature rule exact for polynomials of (1-z) 2 Q 2r+1,2r+1,2r to exactly integrate the mass matrix. ◻ Remark 3.7 Because of the rational fraction due to DF -1 , the stiffness matrix can not be exactly integrated with a classical quadrature formula.

To integrate exactly the mass matrix, we choose the following quadrature formula :

(ξ G i , ξ G j , ξ HM k ), (ω G i ω G j ω HM k ),
where (ξ G i , ω G i ) is the classical Gauss quadrature rule, exact for polynomials of Q 2r+1 , and (ω HM k , ξ HM k ) a Gauss-Jacobi rule for the evaluation of (1 -x) 2 f (x), exact for polynomials of (1 -x) 2 Q 2r+1 . For this last rule, quadrature points and weights have been calculated in Hammer, Marlowe and Stroud [START_REF] Hammer | Numerical Integration Over Simplexes and Cones[END_REF]. Eventually, we need (r + 1) 3 integration points for an exact integration of the mass matrix.

Error Estimates

Functional Spaces and Basic Notations

Let Ω be an open set of R 3 .

Definition 4.1 We define

L 2 (Ω) = {u ∈ L 2 loc Ω u 2 < +∞} H m (Ω) = {u ∈ L 2 (Ω) ∂ α u ∂x α ∈ L 2 (Ω), α ≤ m} equipped with the usual norm in H m (Ω) u 2 m,Ω = α ≤m Ω ∂ α u ∂x α 2
and the usual semi-norm in H m (Ω)

u 2 m,Ω = α =m Ω ∂ α u ∂x α 2 .
We obviously have the inequality

∀u ∈ H m (Ω), u 2 m,Ω ≤ u 2 m,Ω . (15)
An approximate integral using a Gauss-type quadrature formula is denoted ∮ G , and π r will denote a projector on the polynomial space P r .

Presentation of the Problem

We consider a standard variational problem

Find u ∈ V such that ∀v ∈ V, a(u, v) = f (v), (16)
where the space V = H 1 (Ω), a(., .) a continuous bilinear and coercive form, and f (.) a continuous linear form. Then, given a finite-dimensional subspace V h of the space V , the discrete problem reads Find

u h ∈ V h such that ∀v h ∈ V h , a h (u h , v h ) = f h (v h ), (17)
where a h (., .) is a bilinear form defined over the space V h , uniformly V h -elliptic, and f h (.) is a linear form defined over the space V h . We shall consider the simple case where a(u, v)

= Ω uv + Ω ∇u ⋅ ∇v.

Abstract Error Estimate

We will consider the following version of Strang's lemma Lemma 4.2 (Strang's lemma). If u is the solution of ( 16) and u h the solution of ( 17), there exists a constant C > 0 which does not depend on the space step h such that

u -u h 1 ≤ C inf v h ∈V r h u -v h 1 interpolation error + sup w h ∈V r h a(v h , w h ) -a h (v h , w h ) w h 1 numerical integration error .
Proof. The proof of this version of Strang's lemma is similar to the proof proposed by Ciarlet [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] for his version, noticing that, in our case, V h ⊂ V .

◻

We now study separately the two terms of the right hand-side of the Strang's lemma, namely the interpolation error and the quadrature error. As P r ⊂ P F r , we first consider the case of v h = π r u h ∈ P r and get estimates in this case, and we will then take the infimum for v h ∈ P F r to get Strang's lemma. We also suppose that u is in H r+1 (Ω) for an Ω regular enough.

Interpolation error

Let Ω a set composed of n pyramids

K i Ω = ⋃ 1≤i≤n K i . We denote h Ki = diam (K i ) = sup x,y∈Ki x -y , h = sup 1≤i≤n h Ki .
Let us recall the Bramble-Hilbert's lemma (Ciarlet [7])

Lemma 4.3 (Bramble-Hilbert's lemma). Let r ≥ 1, m ≤ r + 1.
We denote by π r the projector to P r satisfying ∀p ∈ P r (K), π r p = p.

Then, there exists a constant C K > 0 which only depends on K and r such that

∀u ∈ H r+1 (K), u -π r u m,K ≤ C K h r+1-m K u r+1,K .
We consider a domain made of a single pyramid K, (Ω = K). For u in H r+1 (Ω), we apply Bramble-Hilbert lemma

u -π r u m,K ≤ C K h K r+1-m u r+1,K ,
where C K depends on the shape of the pyramid K.

We suppose that the constants C Ki for each element K i of Ω are bounded by a constant C,

∀K i ∈ Ω, C Ki ≤ C. ( 18
)
For m = 1, the sum for all the elements K i of Ω and the inequality (15) on norms give the error

estimates on Ω u -π r u 1,Ω ≤ Ch r u r+1,Ω . (19)
Remark 4.4 The condition C K ≤ C is satisfied in the case of periodic meshes of the following sections, since the number of pyramids of different shape is finite. In a more general case, we conjecture that the boundedness of C K is related to the existence of an upper bound for the inverse of jacobian matrix DF , as in the case of hexahedra (Girault and Raviart [START_REF] Girault | Finite Element Approximation of the Navier-Stokes Equations[END_REF]).

Quadrature Error

We now seek the minimal integration formula in order to obtain a quadrature error of order r.

Definition 4.5 For (v h , w h ) ∈ P F r and K ∈ Ω, we denote

E K (v h , w h ) = K v h w h dxdydz -∮ K v h w h dxdydz,
where the approximate integral

∮ K v h w h dxdydz is exact for (1 -z) 2 Q m,m,n .
We first get estimates for the mass term of the bilinear form a.

Lemma 4.6 ∀(v h , w h ) ∈ P F r , K ∈ Ω and 0 ≤ p, q ≤ m -r -1, we have E K (v h , w h ) = E K (v h -π p v h , w h -π q w h )
for a quadrature formula exact for polynomials of (1

-z) 2 Q m,m,m-1
Proof. Let us consider the integral

K π p v h w h dxdydz.
After change of variable, we get

Q 4 (1 -z) 2 (π p ṽh )( DF wh ) dxdỹdz.
For a quadrature rule exact for polynomials of (1-z) 2 Q m,m,m-1 , as π p ṽh DF wh ∈ Q p+r+1,p+r+1,p+r (x, ỹ, z), we have

E K (π p v h , w h ) = 0, (20)
for r + p + 1 ≤ m, that is p ≤ mr -1, using lemma 3.3. Similarly, we have

E K (v h , π q w h ) = 0
for q ≤ mr -1, and for p, q ≤ mr -1, we check that

E K (π p v h , π q w h ) = 0, which provides the claimed result. ◻ Proposition 4.7 For K ∈ Ω and v h ∈ P r (K), sup w h ∈P F r E(v h , w h ) w h 1,K ≤ C ′ K h r K v h r,K
for a quadrature formula exact for polynomials of (1

-z) 2 Q m,m,m-1 , with m ≥ 2r -1.
Proof. We apply lemma 4.6 for p = r -2 and q = 0

E(v h , w h ) = E(v h -π r-2 v h , w h -π 0 w h ).
Using the Cauchy-Schwarz inequality for both the exact and the approximate integral, which is possible since the norm defined by approximate integration is equivalent to the usual norm with a constant C N , we get

E(v h , w h ) ≤ C N v h -π r-2 v h 0,K w h -π 0 w h 0,K .
As v h ∈ P r (K) ⊂ H r (K) and w h ∈ H 1 (K), Bramble-Hilbert's lemma and the inequality on the norms [START_REF] Erlangga | Some numerical aspects for solving sparse large linear systems derived from the Helmholtz equation[END_REF] provides

v h -π r-2 v h 0,K ≤ C K h r-1 K v h r,K . w h -π 0 w h 0,K ≤ C K h K w h 1,K . Hence E(v h , w h ) ≤ C ′ K h r K v h r,K w h 1,K with C ′ K = C N C K independent of w h .
Therefore, we obtain the claimed result. ◻

We then find estimates for the stiffness term of a.

Lemma 4.8 (1 -z) 2 DF DF * -1 ∇ φi (x, ỹ, z) ∈ (1 -z) 2 Q r,r+1,r-1 × Q r+1,r,r-1 × Q r+1,r+1,r-1 (x, ỹ, z).
Proof. As DF DF * -1 is the comatrix of DF , using lemma 3.5, we get

DF DF * -1 = (A 2 ∧ A 3 , -A 1 ∧ A 3 , A 1 ∧ A 2 ) + (C ∧ A 3 , 0, A 1 ∧ C) x + (0, -C ∧ A 3 , C ∧ A 2 ) ỹ +2 (A 2 ∧ C, -A 1 ∧ C, 0) xỹ that is DF DF * -1 ∈ Q 1,1,0 (x, ỹ, z) 3 . As (∇ φi )(x, ỹ, z) ∈ Q r-1,r,r-1 × Q r,r-1,r-1 × Q r,r,r-1 (x, ỹ, z),
by summing degrees, we obtain the claimed result. ◻ Lemma 4.9 ∀(v h , w h ) ∈ P F r , K ∈ Ω and 0 ≤ p ≤ mr -1, we have

E K (∇v h , ∇w h ) = E K (∇v h -π p ∇v h , ∇w h )
for a quadrature formula exact for polynomials of (1

-z) 2 Q m,m,m-2
Proof. As for lemma 4.6, we prove that

E K (π p ∇v h , ∇w h ) = 0 by considering the integral K π p ∇v h ⋅ ∇w h dxdydz.
After the change of variable, we get

Q 4 (1 -z) 2 (π p ∇ṽ h ) ⋅ ( DF DF * -1
∇ wh ) dx dỹ dz.

For a quadrature rule exact for polynomials of (1 -z) 2 Q m,m,m-2 , using lemma 4.8, (π p ∇ṽ h )⋅ ( DF DF * -1

∇ wh ) ∈ Q p+r+1,p+r+1,p+r-1 (x, ỹ, z), that is E K (π p ∇v h , ∇w h ) is equal to zero if r + p + 1 ≤ m, i.e. p ≤ m -r -1. ◻ Proposition 4.10 For K ∈ Ω and v h ∈ P r (K), ∀w h ∈ P F r , E K (∇v h , ∇w h ) = 0
for a quadrature formula exact for polynomials of (1

-z) 2 Q m,m,m-2 , with m ≥ 2r.
Proof. For v h ∈ P r (K), that is ∇v h ∈ (P r-1 (K)) 3 , we apply lemma with p = r -1 to conclude. ◻

We can now get the final quadrature error estimate using propositions 4.7 and 4.10, by summing for all the elements K of Ω and using [START_REF] Hammer | Numerical Integration Over Simplexes and Cones[END_REF] sup

w h ∈V r h (a -a h )(v h , w h ) w h 1,Ω ≤ C ′ h r v h r+1,Ω , v h ∈ P r (Ω) (21)
with

C ′ = C N C.
Theorem 4.11 For a quadrature formula which is exact for polynomials of (1 -z) 2 Q m,m,m-2 , for m ≥ 2r, the final estimate from Strang's lemma is

inf v h ∈V r h ⎛ ⎝ u -v h 1,Ω + sup w h ∈V r h (a -a h )(v h , w h ) w h 1,Ω ⎞ ⎠ ≤ max(C, C ′ ) h r u r+1,Ω .
Proof. Summing equations ( 19) and ( 21) with v h = π r u, we get

u -π r u 1,Ω + sup w h ∈V r h (a -a h )(π r u, w h ) w h 1,Ω ≤ C ′ h r π r u r,Ω + Ch r u r+1,Ω , v h ∈ P r (Ω).
We then use the inequality π r u m,Ω ≤ u m,Ω to upper bound the right-hand side of the inequality with max(C, C ′ ), and we take the infimum for the left-hand side thanks to inclusion [START_REF] Bluck | Polynomial Basis Functions on Pyramidal Elements[END_REF]. ◻ Remark 4.12 We can take the points and weights of quadrature (ξ G k , ω G k ) to get a quadrature formula exact for Q 2r+1 , instead of (ξ HM k , ω HM k ) in the previous quadrature formula, without losing accuracy. Yet, this deteriorates the CFL as we shall see in section 6.2. are enforced, so that the plane wave is a solution of the continuous equations problem. The numerical eigenvalue closest to the wave number ⃗ k is denoted by ω h , so that we define q h as

q h = ω h ⃗ k .
Since q h should be close to 1, we can write

q h = 1 + C h p + o(h p )
where p denotes the dispersion order of the scheme.

A dispersion order of 2r is obtained for both continuous and discontinuous Galerkin method, for regular as well as for distorted meshes (see Babuska and Osborn [START_REF] Babuska | Eigenvalue Problems, Handbook of Numerical Analysis Vol II[END_REF] for the factor 2), which coincides with the theoretical error estimates results obtained. This order of dispersion is clearly shown on the log-log curves in Fig. 6.3 for the continuous hybrid elements on a distorted mesh. Dispersion curves are shown for regular elements of orders one to three in Fig. 6.4 for the continuous elements, and in Fig. 6.5 for the discontinuous elements. For the pyramids, exact and approximate integrations are giving very close results, and all elements exhibit similar dispersion properties. The less dispersive element is the pyramidal one in most of the cases. The same study has been performed for distorted meshes and leads to the same conclusion.

Stability

The stability condition of the standard leap frog scheme is also computed on a periodic infinite mesh.

The CFL number, for which we have the stability condition ∆t ≤ CFL h, is defined by

CF L = 1 max ⃗ k ≤π λ(M -1 ( ⃗ k) K( ⃗ k))
,

where M ( ⃗ k) and K( ⃗ k) are the mass and stiffness matrices associated with the periodic cell, and ⃗ k the wave vector.

For each type of element, the CFL number is given in Table 6.1 in the continuous case in Table 6.1, and in the discontinuous case in 6.2, up to order four. The stability criteria have been tested in the instationary case to check the correctness of the results. The same study has been performed for distorted meshes and leads to the same results, the CFL being smaller for the distorted meshes, and the CFL finally ranks as follows, in all the cases CF L Hexa > CF L W edge > CF L T etra > CF L P yr-ExactInt > CF L P yr-ApproxInt .

The CFL for a hybrid mesh is better than the tetrahedron's and pyramid's, which is a quite surprising result. Remark 6.1 A study of optimal location for the degrees of freedom inside the pyramid has been performed to get an optimal CFL, but this location appears to have almost no influence on the CFL.

Numerical Comparison with Other Existing Methods

We display the dispersion obtained for the existing methods of Sherwin et al. [START_REF] Sherwin | Spectral/hp Methods For Elliptic Problems on Hybrid Grids[END_REF], Nigam and Phillips [START_REF] Nigam | Higher-Order Finite Elements on Pyramids[END_REF], and Bluck and Walker [START_REF] Bluck | Polynomial Basis Functions on Pyramidal Elements[END_REF] on Fig. 6.6, for order two and three. For order one, all these methods provide the same accuracy. The dispersion obtained with the space proposed by Sherwin et al. is of order two, whatever the order of approximation, as the basis functions for the base face and inside the pyramid are not sufficient to contain the optimal space. However, in the affine case, we check that we have a dispersion of order 2r.

The dispersion obtained for the optimal space is equal to the space proposed by Nigam and Phillips, that is the degrees of freedom they add are not necessary as they do not increase the accuracy.

For order two, the dispersion obtained by Bluck and Walker is not consistent (dispersion of order 0) because the space proposed does not contain the space of order one. However, in a case of an affine pyramid, the dispersion obtained is in h 4 for order two. where

Ω = [-1, 1] 3 ω = 1.92 π,
and f is a Gaussian source centered at the origin. We study convergence on a hybrid mesh with a similar pattern as for the dispersion (see Fig. 6.2).

Displaying the L 2 error obtained versus the space step h in a log-log scale in Fig. 6.7, we observe an L 2 error in O(h r+1 ) as proved in the error estimates. For these experiments, we used exact integration with r + 1 HM points in the direction z and r + 1 Gauss points in directions x and ỹ.

Low Storage Matrix-Vector Product

Since the involved matrices can require a huge amount of memory, in particular for high order approximation like Q 5 , we use a low storage matrix vector product (i.e. a matrix-free implementation). This is quite classical in discontinuous Galerkin methods on hexahedra (Castel et al. [START_REF] Castel | Duruflé Discontinuous Galerkin method for hexahedral elements and aeroacoustic[END_REF]) or on tetrahedra (Hesthaven and Teng [START_REF] Hesthaven | Stable Spectral Methods on Tetrahedral Elements[END_REF]), and the extension to pyramid and prismatic elements is straightforward. For continuous formulation, we use similar techniques, i.e. that we exploit factorizations of elementary mass and stiffness matrix

K DF φi φj = Ĉ D h Ĉt K DF DF -1 DF * -1 ∇ φi ⋅ ∇ φj = R B h Rt ,
where matrices Ĉi,j = φi (ξ j ) Ri,j = ∇ φi (ξ j )

are independent on the geometry, so they are not stored for each element, and matrices D h and B h (D h ) i,j = ω j DF (ξ j )δ i,j (B h ) i,j = ω j ( DF DF -1 DF * -1 )(ξ j )δ i,j

are respectively diagonal and block-diagonal (each block being a symmetric 3x3 matrix). Such factorization is explained by Cohen and Fauqueux in [START_REF] Cohen | Mixed Finite Elements with Mass-Lumping for the Transient Wave Equation[END_REF], and shown to be very efficient for hexahedral elements, since R is very sparse. In the case of other elements, R is dense and does not induce any gain in computational time. Yet, we still use this factorization, since the storage induced is very low as we only store matrices D h and B h . To have a good approximation of the geometry, curved isoparametric elements are used. The implementation of such elements is explained in Šolín et al. [START_REF] Šolín | Higher-order finite elements methods[END_REF], except for pyramids for which the extension is straightforward. The reference solution is computed on a refined pure hexahedral mesh with Q 7 elements. On Fig. 6.9, three different meshes used for third order approximation are displayed.

The COCG solver (Clemens and Weiland [START_REF] Clemens | Iterative methods for the solution of very large complex symmetric linear systems of equations in electrodynamics[END_REF]) is implemented, and can be used with or without preconditioning. In Table 6.3, the required number of degrees of freedom necessary to reach an error between one and two percent (measured in L 2 norm) are displayed for each type of mesh and at orders two, three and four. The results obtained without preconditioning and for a p-multigrid preconditioning using a damped Helmholtz equation (see Erlangga [START_REF] Erlangga | Some numerical aspects for solving sparse large linear systems derived from the Helmholtz equation[END_REF] for finitedifference, and Duruflé [START_REF] Duruflé | Intégration numérique et éléments finis d'ordre élevé appliqués aux équations de Maxwell en régime harmonique[END_REF] for finite element) are also displayed, along with the computational time. The order of the coarsest mesh is set to 1 for orders two and three, and set to 2 for order five and P 4 .

The performance of hybrid mesh is quite similar to a purely hexahedral mesh while purely tetrahedral meshes result in much more expensive computations. The last row of this table concerns the use of a Gauss-Seidel smoother instead of the Jacobi smoother used in other rows. Jacobi smoother fails for hybrid meshes, but we don't have any explanation on this issue. In this case, a Q 5 approximation on hexahedral mesh is much faster than for tetrahedral elements. This is the reason why it is important to use hybrid meshes with a high percentage of hexahedra.

where Γ has the shape of the resonance cavity of a piano, and F is the surrounding parallelepiped box, as displayed in Fig. 6.10.

The source is chosen as f (x, t) = 1 r 2 0 e -13 r r 0 2 e -4(t-t0) 2 sin(2πf 0 t),

where r is the distance to the center of the source, r 0 the distribution radius of the Gaussian,

Fig. 1 . 1 .Remark 1 . 2

 1112 Fig. 1.1. Transformation of the reference pyramid K to the pyramid K via the transformation F

Fig. 1 . 3 .

 13 Fig. 1.3. Location of the degrees of freedom inside the pyramidal element of order 5

Fig. 1 . 4 .

 14 Fig. 1.4. Location of the degrees of freedom for the pyramidal elements of order 2, 3 and 4 (Color online)

  Fig. 1.5. Condition number of the VDM matrix versus the order for tetrahedral and pyramidal elements, for monomial and orthogonal basis functions Proposition 2.1 The subspace C 0 r of C r with zero trace on the boundary of Q is

Fig. 6 . 2 .

 62 Fig. 6.2. Periodic pattern for the hybrid case, with distorted pyramids (purple) and tetrahedra (gray) (Color online)

Fig. 6 .

 6 Fig. 6.3. log-log dispersion error for continuous finite elements of orders 1 to 4 for a hybrid distorted mesh

Fig. 6 . 5 .

 65 Fig. 6.5. Dispersion curves for discontinuous Galerkin method of orders 1 to 3 for regular meshes (K = 6kh 2πr )

Fig. 6 . 6 .

 66 Fig. 6.6. Dispersion errors for the different types of existing elements of orders 2 and 3 with a distorted pyramidal mesh

6. 4 .

 4 Numerical Experiments 6.4.1. Test Case on a Cube We first consider the Helmholtz equation on a cubic cavity with homogeneous Dirichlet boundary conditions -ω 2 u -∆u = f (x, y, z),

Fig. 6 . 7 .

 67 Fig. 6.7. L 2 error versus h for the cubic cavity at different order of approximations.

6. 4 . 3 .

 43 Test Case with Curved Isoparametric Elements on a Sphere Let us consider a sample test case of scattering by a sphere (see Fig 6.8)⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ -ω 2 u -∆u = 0 in Ω ∂u ∂n = -∂u incident ∂n on Γ ∂u ∂n -iωu = 0 on Σ,where Γ is a sphere of radius 3, ω = 2π, and Σ is the boundary of the cube [-5, 5] 3 .

Fig. 6 . 8 .

 68 Fig. 6.8. Real part of diffracted field for a sphere of radius 3 with Neumann condition.

Fig. 6 . 9 .

 69 Fig.6.9. Meshes used for third order approximation.

Table 6 .

 6 1: CFL for continuous finite elements with regular meshes.

		Order 1 Order 2 Order 3 Order 4
	Hexahedron	0.28868 0.11785 0.06697 0.04264
	Wedge	0.15794 0.07176 0.04344 0.02926
	Pyramid ExactInt	0.09682 0.04803 0.03083 0.02143
	Pyramid ApproxInt 0.07217 0.03335 0.01985 0.01316
	Tetrahedron	0.12247 0.06253 0.03919 0.02669
	Hybrid	0.15138 0.07245 0.04558 0.03191

Table 6 .

 6 2: CFL for discontinuous Galerkin method with regular meshes.

		Order 1 Order 2 Order 3 Order 4
	Hexahedron	0.14434 0.07144 0.04348 0.02934
	Wedge	0.11471 0.04348 0.03957 0.02717
	Pyramid ExactInt	0.07184 0.04058 0.02618	0.0184
	Pyramid ApproxInt 0.04811 0.02544 0.01566	0.0112
	Tetrahedron	0.07373 0.04467 0.03041 0.02173
	Hybrid	0.09283 0.05363 0.03527 0.02448
	For pyramidal elements, the CFL are computed with the quadrature formula (ξ HM k presented in paragraph 4.3.2 for exact integration, and with Gauss quadrature formula (ξ G , ω HM k k , ω G k ) )
	for approximate integration. The CFL of pyramidal elements are clearly lower for approximate
	integrals (60% lower on average) than for exact integrals. In this case, the CFL of pyramidal
	elements is close to the tetrahedral's.

Journal of Scientific Computing

Acknowledgments. We thank P. Ciarlet and S. Imperiale for useful suggestions and comments on the manuscript.

Discontinuous Galerkin Method

For any element K of boundary ∂K, we consider the Local Discontinuous Galerkin (LDG) formulation (Castillo et al. [START_REF] Castillo | An a-priori Error Analysis of the Local Discontinuous Galerkin Method for Elliptic Problems[END_REF])

, where u 1 is the value of u restricted to K, and u 2 the value of u restricted to an element adjacent to K.

We can take the space Pr identical to the one in the continuous case. The error estimates are quite different from the ones obtained for continuous elements as the norm used for DG formulations is the L 2 norm, instead of the H 1 norm.

Numerical Results

Dispersion

In order to study the pyramidal elements, a dispersion analysis is performed on the wave equation, relying on the computation of the phase error on infinite periodic meshes, as in Cohen [START_REF] Cohen | Higher-Order Numerical Methods for Transient Wave Equations[END_REF]. The periodic cell is a cube that can be composed of a single hexahedron; of two wedges; of two pyramids and two tetrahedra (hybrid); or of six pyramids; or of six tetrahedra as shown in Fig. 6.1. The analysis has also been carried out on periodic cells made up of distorted cubes in order to check the consistency of our method when the base of the pyramid is not a parallelogram as shown in Fig. 6.2 for the hybrid cell. , where ⃗ k is the wave vector, and ⃗ k the wave number. The computational domain is reduced to a periodic cell where quasi-periodic conditions f 0 the frequency, and t 0 a constant. Here we have taken

so that the parallelepiped box is as large as 32λ × 26λ × 10λ where λ = 1 f 0 is the wavelength.

We compute the solution from t = 0 until t = 6, and we obtain the result of Fig. 6.11.

A second-order leap frog scheme (Cohen and Fauqueux [START_REF] Cohen | Mixed Finite Elements with Mass-Lumping for the Transient Wave Equation[END_REF]) is used for the time discretization.

The reference solution is computed on a very fine mesh, and we compare two kind of meshes: a hybrid mesh and a hexahedral mesh obtained by splitting tetrahedra. Third order approximation (Q 3 ) is used. The results are given in Table 6.4. We have given the computational time we would have obtained on a single processor, this time is computed by summing the computational times for all the processors and we subtract the cost of communications. Since curved elements are not used, hexahedral meshes generated by splitting tetrahedra produce a bad approximation of the geometry, therefore they require a larger number of degrees of freedom.

Conclusion

Highly efficient pyramidal elements of any order have been constructed. Numerical experiments conducted with these elements (up to order six) exhibit a low phase error, a good CFL, and a very good behaviour in hybrid meshes.


Amandine Michelas
email: michelas@lpl-aix.fr

Mariapaola D'imperio
email: mariapaola.dimperio@lpl-aix.fr

Phonetic cues and phonological constraints in prosodic unit definition: evidence for the intermediate phrase in French

The studies referring to prosodic constituency in French generally agree on two levels of phrasing -the IP (Intonation Phrase) and the AP (Accentual Phrase) in [START_REF] Jun | The accentual phrase and the prosodic structure of French[END_REF][START_REF] Jun | A phonological model of French intonation[END_REF]-while the existence of an intermediate level of phrasing (ip) between these two units seems restricted to specific intonation structures observed in marked syntactic constructions [START_REF] Jun | The Realizations of the accentual phrase in French intonation[END_REF]. The two experiments reported here provide evidence for the existence of an ip which is not restricted to marked syntactic structures but which can appear within all focus utterances if both syntactic structure and prosodic constraints allow it. The analysis here is grounded in the framework of Generalized Alignment Theory [START_REF] Prince | Optimality Theory: Constraint Interaction in Generative Grammar[END_REF]. The analysis includes two important constraints, which reflect the prosodysyntax interface: 1) ALIGN XP,R,ip,R: align the right edge of a maximal projection with the right edge of an ip [START_REF] Prieto | Syntactic and eurhythmic constraints on phrasing decisions in Catalan[END_REF]; and 2) a prosodic constraint stating that non-final intermediates phrases consist of minimally two APs. We specifically predict that the interaction between these two constraints conspires to place an ip-right boundary made of minimally two APs in correspondence with a major syntactic break such as the boundary between a Subject NP and a VP.

In a first experiment, we focus on phonetic properties (duration and F0 values) of syllables occurring 1) within a prosodic word, 2) at an AP-boundary which is not aligned with the right edge of a syntactic XP (AP), 3) at an AP-boundary which is aligned with the right edge of a syntactic XP (AP/XP) and 4) at an IP-boundary. The corpus was read at two speech rates (normal/fast). The results of a mixed model show that the right edge of an AP/XP boundary appears to be marked by a significant syllable lengthening (Fig. 1) at normal speech rate (Fig. 2; t=2.83, p<0.05) whereas F0 values of AP/XP syllables are not significantly higher than AP final syllables (Fig. 3;t=0.85,p>0.05). While the alignment between syntactic and prosodic structure seems to reinforce preboundary lengthening, the ALIGN-XP,R constraint is not sufficient to account for the emergence of an the boundary itself, since the subject NP is always produced as a single Accentual Phrase.. In a second experiment, we investigated the duration and F0 values of all focus SVO utterances in which the structure of the subject NP was varied in order to obtain either a sequence of 2 or of 3 AP as exemplified in (Fig. 4). The corpus was read at two speech rates (normal/fast). We measured duration and F0 values of V2 (the AP 2 final vowel) in both 2APs and 3APs conditions. F0 values were included as ratios of the relevant H target divided by the F0 height of the first H in the utterance [START_REF] D'imperio | Interface entre structure syntaxique et structure prosodique: le syntagme intermédiaire en français. Interface, Discours & Prosodie[END_REF]. This was done in order to normalize for pitch range variation across repetitions. Note also that V2 is associated with an ip-boundary in the 2 APs condition while it is only associated with an AP-boundary in the 3APs condition. We also measured the height of the first LH* within the following ip (Hpb in Fig 4). Note that for both normal and fast speech rates, V2 is significantly longer in the 2APs condition than in the 3APs condition (normal speech rate t=-8.487, p<0.05; fast speech rate t=-3.250, p<0.05) and that H2 is significantly higher in the 2APs than in the 3APs condition (Fig. 5; normal speech rate: t=-3.67, p<0.05; fast speech rate: t=-9.43, p<0.05). We also found total pitch reset occurring before the ip-boundary, while we only found partial reset at the first LH* peak following the intermediate phrase boundary (Hpb; Fig. 6) as already observed for Germanic languages [START_REF] Ladd | Declination "Reset" and the Hierarchical Organization of Utterances[END_REF][START_REF] Van Den Berg | Downstep in Dutch: Implications for a model[END_REF][START_REF] Truckenbrodt | Embedded register levels[END_REF].

Hence, we propose that the intermediate phrase is signaled in French by (i) preboundary lengthening (ii) the presence of H-phrase accent which is responsible for blocking recursive downstep of subsequent AP final rises within the ip and (iii) partial reset across the ip boundary. The analyses show that both local phonetic cues and phonological constraints have great importance in the ip definition for French. An interesting interplay of duration and pitch height is responsible for signaling the boundary and will be investigated through perceptual experiments.

  Fig. 2 : Syllable lengthening in ms across boundary type for both speech rates for Experiment 1.

Fig. 1 :

 1 Fig. 1 : Target syllables in Experiment 1

Fig. 3 :

 3 Fig. 3 : F0 height across boundary type for both speech rates for Experiment 1.

Fig. 4a :

 4a Fig.4a: F0 curve for the sentence La mamie de Rémy demandait l'institutrice ("Remy's grandmother asked for the teacher") where the subject NP is made of 2 APs.