Asma Gdhami

Devant Le

Amel Ben Abda

Aux Dérivées

Partielles Hyperboliques

Professeur Mr Christophe Chalons

Mrs Amel

Professeur Ben Abda

Mrs Claire Scheid

Mr Maher

Petrov-Galerkin Rk

Runge-Kutta Petrov-Galerkin

Méthodes isogéométriques pour les équations aux dérivées partielles hyperboliques

Keywords: problèmes hyperboliques, méthode des éléments finis, méthode de Galerkin discontinue, analyse isogéométrique, fonctions B-splines, extraction hyperbolic problems, Finite Element method, discontinuous Galerkin method, Isogeometric analysis, B-spline functions, Bézier extraction, curve fitting, least squares method Ordinary differential equation. PDE Partial differential equation

HYPERBOLIQUES

R Esumé: L'nalyse isogéométrique (AIG) est une méthode innovante de résolution numérique des équations différentielles, proposée à l'origine par Thomas Hughes, Austin Cottrell et Yuri Bazilevs en 2005. Cette technique de discrétisation est une généralisation de l'analyse par éléments finis classiques (AEF), conçue pour intégrer la conception assistée par ordinateur (CAO), afin de combler l'écart entre la description géométrique et l'analyse des problèmes d'ingénierie. Ceci est réalisé en utilisant des B-splines ou des B-splines rationnelles non uniformes (NURBS), pour la description des géométries ainsi que pour la représentation de champs de solutions inconnus.

L'objet de cette thèse est d'étudier la méthode isogéométrique dans le contexte des problèmes hyperboliques en utilisant les fonctions B-splines comme fonctions de base. Nous proposons également une méthode combinant l'AIG avec la méthode de Galerkin discontinue (GD) pour résoudre les problèmes hyperboliques. Plus précisément, la méthodologie de GD est adoptée à travers les interfaces de patches, tandis que l'AIG traditionnelle est utilisée dans chaque patch. Notre méthode tire parti de la méthode de l'AIG et la méthode de GD.

Les résultats numériques sont présentés jusqu'à l'ordre polynomial p = 4 à la fois pour une méthode de Galerkin continue et discontinue. Ces résultats numériques sont comparés pour un ensemble de problèmes de complexité croissante en 1D et 2D.

ISOGEOMETRIC METHODS FOR HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS

A Bstract: Isogeometric Analysis (IGA) is a modern strategy for numerical solution of partial differential equations, originally proposed by Thomas Hughes, Austin Cottrell and Yuri Bazilevs in 2005. This discretization technique is a generalization of classical finite element analysis (FEA), designed to integrate Computer Aided Design (CAD) and FEA, to close the gap between the geometrical description and the analysis of engineering problems. This is achieved by using B-splines or non-uniform rational B-splines (NURBS), for the description of geometries as well as for the representation of unknown solution fields.

The purpose of this thesis is to study isogeometric methods in the context of hyperbolic problems using B-splines as basis functions. We also propose a method that combines IGA with the discontinuous Galerkin (DG) method for solving hyperbolic problems. More precisely, DG methodology is adopted across the patch interfaces, while the traditional IGA is employed within each patch. The proposed method takes advantage of both IGA and the DG method.

ACKNOWLEDGMENTS

T Hese PhD years have been a rich experience from the professional as well as the personal points of view.

First of all I would like to express my gratitude to my supervisors, Regis Duvigneau & Maher Moakher for inspiring me, introducing me to new challenges and for helping me carry out this work. I would like to thank them for showing much interest in my work, and for their help and support at the frequently advising sessions. At your sides, I have learnt a lot. You shared with me your technical knowledges and above all your scientific rigorous. In particular, thanks to you, I got the opportunity to attend international congresses.

It represents a lot to me. I would like to express my sincere thanks for your patience, your help, your regular availability and your encouragement and advice. I dedicate this work to you reflecting my deep respect.

I express my heartfelt thanks to Mrs. Amel BEN ABDA & Mrs. Claire SCHEID for being my thesis examiner and providing valuable suggestions and corrections. Also, I present my highly express of thanks to Sir Christophe CHALONS & Sir Maatoug HASSINE for agreeing to be the reviewers of my work and for your constructive suggestions.

Then, I address a special appreciation for Sir. Mekki Ayadi that helped me to accomplish this work by extensive discussions. I would like to express my sincere thanks for your patience, your help, your encouragement and advice.

My next big thank you is to my families Gdhami, Essid and Azaouzi. To my parents Mohamed and Kaouther, my brother and my sister for their unconditional love and support all though my life. To them I owe all that I am and all that I have ever accomplished and it is to them I dedicate this thesis. Moreover, I will never forget my cousins Sana, Ichraf, Boutheina, Soumaya & Olfa who have always encouraged me to accomplish my studies.

vii On a more personal note, I would like to express my deep gratitude to my husband Aymen for his support. This work is dedicated to him and to my son with my deepest love. I address a special message for my second family: my friends from the Modeling Laboratory in Engineering Sciences (LAMSIN) National Engineering School of Tunis (ENIT), Boutheina, Imen, Rabeb, Maroua, hamouda, Anis, ... and my friends of the National Institute for Research in Computer Science and Automatic of Nice (INRIA), I would like to express my gratitude for their continuous encouragement.

In the life path, we meet people that can change our life in a different measure -some change it forever. I am sure I would have forgotten some people in this acknowledgement section, that is why for the sake of completion I am thanking all the people who shared with me parts of these years. Thanks for being part of my path. In a certain way, that you maybe do not measure, you contribute to allow me to arrive here today.

N el

Finite number of cells of DG method.

N el

Finite number of patchs of IGDG method.

N Dimension of the polynomial approximation.

B k p k-th Bernstein polynomial of degree p in parameter space.

B k p k-th Bernstein polynomial of degree p in physical space.

N i ,p i -th B-spline function of degree p in parameter space.

N i ,p i -th B-spline function of degree p in physical space.

u k + h
The field on the exterior of the element boundary Γ k .

u k - h
The field on the interior of the element boundary Γ k .

C C F L

Courant-Friedrichs-Lewy (CFL) number.

c Advection velocity.

-→ c = (c x , c y ) Velocity field.

-→ n

The outward unit normal.

-→ n e

The outer unit normal to Γ e of the element Ω e . T Final time.

M

Mass matrix for the P 1 FEM.

M -1

Inverse of mass matrix.

M s L SUPG lumped mass matrix.

R

Stiffness matrix for the P 1 FEM.

M s

Mass matrix for the P 1 SUPG-FEM.

M s

Mass matrix for the P 2 SUPG-FEM.

R s

Stiffness matrix for the P 1 SUPG-FEM.

R s

Stiffness matrix for the P 2 SUPG-FEM.

M B

Mass matrix for the B-spline IGFEM.

xxiv (M B ) T Transpose of the mass matrix for the B-spline IGFEM.

R B

Stiffness matrix for the B-spline IGFEM.

Q

The matrix defined by the least square method.

M k

Local mass matrix.

R k

Local stiffness matrix.

M k

Local mass matrix for IGDGM.

R k

Local stiffness matrix for IGDGM.

J

The Jacobian matrix.

| J | Determinant of the Jacobian.

J e

The elemental Jacobian matrix in the physical domain Ω e .

∂ x u = ∂u ∂x Partial derivative of u with respect to space x.

∂ t u = ∂u ∂t
Partial derivative of u with respect to time t .

⊗

Tensor product.

C k

The set of functions with kt h order continuous derivatives.

L Differential operator.

T Transformation of the parametric domain to the physical domain.

T -1

Inverse of transformation T.

F e

Numerical flux in the patch Ω e .

f C E N Central flux.

f G Godunov flux.

INTRODUCTION

H Yperbolic systems of partial differential equations (PDEs) are mathematical mod- els expressing the conservation of a physical quantity, as for instance mass, energy, etc. They arise naturally from the conservation laws in physics. In particular, they describe a wide variety of phenomena that involve wave motion (such as acoustic, elastic, electromagnetic) or the advective transport of substances.

The wide range of applications of hyperbolic PDEs led to a very intense research activity in this field. It allowed to develop very early a set of numerical methods for accurate and computationally efficient approximations to the solutions of such problems. There are three major families of methods which are widely used:

the finite difference method, the finite volume method and the finite element method. These methods have proved to be extremely useful in modeling a broad set of phenomena. To keep this thesis self-contained, we briefly introduce each of these three methods in the context of hyperbolic PDEs.

Historically, the finite difference method (FDM) was the first method used to produce approximations of the solutions of hyperbolic PDEs. They were introduced by Euler in the 18t h century and represent the easiest method to solve problems on simple geometries [START_REF] Morton | Numerical solution of partial differential equations: an introduction[END_REF]. The main idea of this method, is to replace the functional derivatives of the unknown by their FD approximations. The FDM is notable for the large variety of schemes that can be used to approximate a given PDE; e.g. explicit schemes (Forward Euler, Upwind, Lax-Friedrichs, Lax-Wendroff, Leapfrog, ...) and implicit schemes (Backward Euler, Crank-Nicolson,...) [START_REF] Trefethen | Finite difference and spectral methods for ordinary and partial differential equations[END_REF].

Although this method can be easily formulated and implemented, its application to problems with realistic geometries is rather cumbersome, thus making the method not very attractive for industrial problems.

This fact urged the need for other methods with more flexibility, such as finite volume and finite element methods.

The finite volume method (FVM) can handle complex geometries which makes it more attractive for complex problems than the FDM. It is based on the conservative form instead of the differential form to estimate the values of unknown fields. We split the domain into grid cells and approximate the total integral over each grid cell, or actually the cell average, which is this integral divided by the volume of the cell. The links between the cell quantities in the FVM rely on the flux between neighbouring control volumes. This means that the FVM represents the flux of information within the structure of the mesh in a conservative way [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. It allows the use of unstructured grids to handle complex geometries.

The finite element method (FEM), had its origins in the early 1960s and is nowadays the predominating method in analysis of elliptic or parabolic problems due to its flexibility to represent complex geometric domains and its strong theoretical basis. The idea consists in decomposing the domain into many small, "finite" elements which are defined by a set of nodal points and interpolating basis functions.

Although the FEM has been used widely in simulating many physical phenomena due to its flexibility to represent complex geometric domains, it is well known in the FE literature that numerical difficulties arise when solving hyperbolic PDEs. Indeed, when using the standard Galerkin FE method applied to hyperbolic PDEs, unwanted spurious (non-physical) oscillations (Gibbs phenomenon) are frequently detected in the numerical solutions. A cure to this drawback, widespread in the literature, is to add some "artificial" viscosity to a standard (unstable) numerical scheme. On the one hand, this artificial viscosity should damp the oscillations but, on the other hand, it should not smear the numerical solution. In the late 1970s and early 1980s, a large number of so-called stabilized methods have been developed with different ideas [START_REF] Brooks | Streamline Upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF] [32] [START_REF] John | On spurious oscillations at layers diminishing (SOLD) methods for convection-diffusion equations: Part I I -Analysis for P 1 and Q 1 finite elements[END_REF]. This is achieved through the use of a Petrov-Galerkin formulation [START_REF] Griffiths | An analysis of the Petrov-Galerkin finite element method[END_REF], where the test functions are modified such that they weight the upstream node more than the downstream node [17] [38]. Among them, the most popular, so called Streamline-Upwind Petrov-Galerkin (SUPG) method, was introduced by Brooks and Hughes. It was first proposed in the context of advection-diffusion equations and incompressible Navier-Stokes equations [START_REF] Brooks | Streamline Upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF], and then extended to various other problems, e.g., coupled multidimensional advective-diffusive systems [START_REF] Hughes | A new finite element formulation for computational fluid dynamics: v. Circumventing the Babuška-Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations[END_REF], first-order linear hyperbolic systems [START_REF] Johnson | Finite element methods for linear hyperbolic problems[END_REF] or first-order hyperbolic systems of conservation laws [START_REF] Hughes | Finite element methods for first-order hyperbolic systems with particular emphasis on the compressible Euler equations[END_REF].

Later, Galerkin/least-squares (GLS) has emerged as a generalization of the SUPG method, developed by Hughes et al. for convective transport problems [START_REF] Billaud | Eléments finis stabilisés pour des écoulements diphasiques compressible-incompressible[END_REF] [30] [START_REF] Shakib | A new finite element formulation for computational fluid dynamics: x. the compressible Euler and Navier-Stokes equations[END_REF], in which residuals of the equations in leastsquares form are added to the standard Galerkin formulation. GLS has been successfully employed in a wide variety of applications where enhanced stability and accuracy properties are needed, including problems governed by Navier-Stokes and the compressible Euler equations in fluid mechanics [START_REF] Shakib | A new finite element formulation for computational fluid dynamics: x. the compressible Euler and Navier-Stokes equations[END_REF].

We can also mention the pressure-stabilizing/Petrov-Galerkin (PSPG) [START_REF] Wervaecke | Simulation d'écoulements turbulents compressibles par une méthode d'éléments finis stabilisée[END_REF] formulation which has been introduced for the stabilization of the Stokes equations [START_REF] Hughes | A new finite element formulation for computational fluid dynamics: v. Circumventing the Babuška-Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations[END_REF] and incompressible Navier-Stokes equations [START_REF] Tezduyar | Incompressible flow computations with stabilized bilinear and linear equal-order-interpolation velocity-pressure elements[END_REF].

The main idea of all these methods is to transform the original Galerkin method into a Petrov-Galerkin formulation adapted to the physics considered. In fact, these formulations stabilize the method without introducing excessive numerical dissipation. Because its symptoms are not necessarily qualitative, excessive numerical dissipation is not always easy to detect. This concern makes it desirable to seek and employ stabilized formulations developed with objectives that include keeping numerical dissipation to a minimum. In these stabilized formulations, judicious selection of the stabilization parameter, which is almost known as τ, plays an important role in determining the accuracy of the formulation. This yielded a significant amount of attention and research [14] [30] [82] [START_REF] Tezduyar | Finite element formulations for convection dominated flows with particular emphasis on the compressible Euler equations[END_REF]. Typically this stabilization parameter involves a measure of the local length scale (also known as "element length") and other parameters such as the local Reynolds and Courant numbers. However, this stabilization parameter requires special attention, as it strongly depends on the problem under consideration and the chosen numerical method.

More recently, an alternative approach has emerged, the discontinuous Galerkin method (DGM), which shares some features with both the FVM and the FEM. Indeed, discontinuous polynomial functions are used and a numerical flux is defined at the interface between cells to reconstruct the solution. It has been proved

very useful in solving a large range of problems. It was first introduced in 1973 by Reed and Hill for solving a time-independent linear hyperbolic equations [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF] and, later on, it has been extended for solving nonlinear time-dependent equations. Subsequently, during the nineties of the last century the DGM experienced a series of developments by Cockburn and Shu, where numerical schemes for hyperbolic problems were proposed by combining DG approximation in space with Runge-Kutta time stepping strategies [START_REF] Cockburn | The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. I V . The multidimensional case[END_REF] [19] [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws v: multidimensional systems[END_REF].

In fact, the DG method combines the advantages of stability of FVM and the accuracy of continuous FEM.

Within each element, the solution is approximated by a polynomial of degree p ≥ 0 (as in FEM), while the continuity conditions applied to the solution are relaxed at the boundaries of elements (as in FVM), however, motivated by the FVM, interface terms of the problem are approximated by a consistent, monotone and Lipschitz continuous numerical flux. This ensures that the scheme obtained is conservative, which does not hold in case of the classical FEM. In particular, the increasing interest in these kind of formulations are due to the following interesting features: they have good stability properties, they offer flexibility in the mesh construction (irregular meshes are admissible) and in the handling of boundary conditions (Dirichlet boundary conditions are weakly imposed), the accuracy is obtained by means of high-order polynomials within elements, without any regularity constraint at element interfaces. Furthermore, they are locally conservative.

The fundamental difference between the DGM and the classical FEM relies on the continuity of basis functions. In comparison with the classical FEM, in DGM the basis functions are completely discontinuous across each element interface and they consist of local piecewise polynomials. Due to the fact that of basis function have compact support, integration can be achieved locally in each element. This simplifies the implementation of the method, since the mass matrix becomes block diagonal and the solution of a large system is avoided. In addition, the discontinuity across each element allows the use of different degrees of freedom in each element independently, which is not allowed in classical finite element method. Consequently, we can easily apply adaptivity strategies by increasing the degrees of freedom near phenomena of interest to obtain better approximations to the solution.

As explained above, FEM decomposes the computational domain into many small, "finite" elements with simple shapes. A drawback is that with such elements there is usually no continuity higher than C 0 between elements. Even with higher-order polynomials it is difficult to guarantee C 1 continuity for arbitrarily shaped elements. With the advancement in design technology a more accurate and flexible handling of the geometry becomes necessary.

The design of free-from shapes by mathematical methods is a discipline, named computer-aided-design (CAD). It had its origins slightly later than the computer-aided-engineering (CAE). In fact, CAD is the use of computer technology for design: it allows the creation, modification, analysis and optimization of drawings and geometric modeling. The Bézier curve was the first method used to construct free-form curves and surfaces, and is named according to its inventor, Dr. Pierre Bézier. ) which are basic elements in geometric modeling, but also allows very flexible modeling of free-form surfaces [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF].

Today, there is a strong need for reducing the gap between CAD and FEM in terms of geometric representations to gain in accuracy, flexibility and ease of interaction. A new form of analysis, named isogeometric analysis (IGA), tries to close this gap between CAD and FEA in such a way that both disciplines work on the same geometric models.

IGA is an extension of the FEM for solving PDEs. It was first introduced in 2005 by Hughes, Cottrell,

and Bazilevs [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF], and expanded in 2006 [START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF] in an effort to bridge the gap between FEM and CAD. The key idea is to use for analysis the same geometry used for geometric modeling. In fact, we use the same basis functions, which are used for the representation of the geometry in computer aided design (CAD) models, also for the approximation of the solution of the PDE or the system of PDEs describing the physical phenomenon. The idea of using Bézier, B-splines or NURBS as basis functions is driven by the desire to integrate CAD within FEM, and to have a strategy to replace a huge number of little cells (the FEs) by a reduced set of larger patches covering the entire domain. Moreover, there are several advantages of this approach over the FEM: easily control of the continuity, as C p-1 -continuity is obtained using p-th order NURBS, exact representation of the underlying NURBS geometry on the coarsest level of discretisation, as well as exact representation of the geometry as the mesh is refined [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF].

IGA has been applied to a wide variety of different physical phenomena, including computational solid dynamics problems, computational fluid dynamics [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF], coupled solid-fluid interaction problems [START_REF] Bazilevs | Computational fluid-structure interaction: methods and applications[END_REF] and the diffusion equation [START_REF] Hall | The application of isogeometric analysis to the neutron diffusion equation for a pincell problem with an analytic benchmark[END_REF]. In last years, there has been an increasing interest in DG-IGA for the numerical solution of elliptic PDEs [START_REF] Hofer | Dual-primal isogeometric tearing and interconnecting solvers for large-scale systems of multipatch continuous Galerkin IgA equations[END_REF] [54] [69] [START_REF] Zhang | Discontinuous Galerkin methods for isogeometric analysis for elliptic equations on surfaces[END_REF]. The advantages of the local approximation spaces without continuity requirements that DG methods offer [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF] [28] is thus employed to manage multi-patch computations.

The main purpose of this thesis is to study the use of IGA to solve some hyperbolic problems. In particular, we describe the continuous and discontinuous Galerkin method using a B-spline basis. Special emphasis is on the discontinuous Galerkin method, since it is considered as one of the most powerful and fastest growing methods with applications in various problems, not necessarily hyperbolic. The discontinuity of basis functions, which provides more flexibility in analysis, makes the method tedious however for handling realistic geometries from CAD.

The thesis is structured in four main parts: the first gives particular focus on the Bernstein and B-splines basis functions used in CAD. It is devoted to giving their definitions and basic properties. We present in the second and third parts the extension from classical analysis to IGA for the FE and DG methods, for the onedimensional advection problem. In the last part, we deal with two dimensional hyperbolic problems by combining the IGA with the DG method. It should be mentioned, in all this work, that the discretization of equations in time is done by means of high-order explicit Runge-Kutta methods [START_REF] Cockburn | The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. I V . The multidimensional case[END_REF] [19] [34] [80].

More precisely, chapter 2 and chapter 3 provide a comprehensive introduction to the main ideas and properties of the Bernstein, B-splines and NURBS, which form the basis for the IGA. In the same context, an analysis of fitting B-spline curve and surface in the least squares sense is presented in chapter 4. The analysis is illustrated by examples of univariate and bivariate problems. Since IGA is an extension of FEM, we start by revisiting the original analysis framework in chapter 5, i.e. FEM. The need for stabilization is outlined and stabilization ideas based on the Petrov-Galerkin concept are discussed. We focus on the SUPG stabilization method and a special attention is given to the study of the stabilization parameter τ.

In chapter 6, the various computational procedures for IGA are reviewed in the context of FEM, by revisiting the one-dimensional advection problem that is given in the previous chapter. While in this chapter we use B-splines (due to the simplicity of the domain) as a basis function, it is not hard to generalize it to other splines such as NURBS. Detailed comparisons between both IGA and classical FEM are discussed.

In this context of IGA, we consider then the application of DG methods. Indeed, the major argument for using DG methods lies with their ability to provide stable numerical methods for hyperbolic PDE problems, for which classical FEM is well known to perform poorly. Therefore, in chapter 7, we deal with one-dimensional advection problem by combining IGA method with the DG method. We note that the DG methodology is adopted at patch level, i.e., we employ the classical IGA within each patch, and employ the DG method across the patch interfaces. Moreover, a transformation of the B-spline basis is necessary to introduce discontinuities at the interfaces, without modifying the geometry of the domain. The advantageous features of both IGA and DG method enable us to design a promising formulation.

With some adjustments, chapter 8 and chapter 9 are devoted to the study of two numerical examples in 2D. The advection problem is first presented, followed by the acoustic wave equations, where both systems are solved over several domains (Cartesian, linear and curvilinear).

Finally, in chapter 10 we end with some concluding remarks and outlooks. The results of the various studies performed are summarized and discussed, and ideas for future research are proposed.

INTRODUCTION

L Es systèmes d'équations aux dérivées partielles (EDPs) hyperboliques sont des mod- [START_REF] Brooks | Streamline Upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF], puis elle a été étendue à divers autres problèmes, par exemples les systèmes advectifs-diffusifs multidimensionnels couplés [START_REF] Hughes | A new finite element formulation for computational fluid dynamics: v. Circumventing the Babuška-Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations[END_REF], les systèmes hyperboliques linéaires du premier ordre [START_REF] Johnson | Finite element methods for linear hyperbolic problems[END_REF], les systèmes hyperboliques de lois de conservation [START_REF] Hughes | Finite element methods for first-order hyperbolic systems with particular emphasis on the compressible Euler equations[END_REF].

Plus tard, la méthode de Galerkin/moindres carrés (GLS) est apparue comme une généralisation de la méth- représentation exacte de la géométrie en utilisant une NURBS au niveau de discrétisation le plus grossier, ainsi que la représentation exacte de la géométrie lorsque le maillage est affiné [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF].

L'AIG a été appliquée à une grande variété de phénomènes physiques, y compris la dynamique des fluides computationnelle [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF], les problèmes d'interaction couplé fluide-structure [START_REF] Bazilevs | Computational fluid-structure interaction: methods and applications[END_REF] et l'équation de diffusion [START_REF] Hall | The application of isogeometric analysis to the neutron diffusion equation for a pincell problem with an analytic benchmark[END_REF].

Part I CAD REPRESENTATIONS

C H A P T E R 2 BÉZIER CURVES
B Ézier curves are parametric curves commonly used in computer graphics and re- lated fields. They are named after their inventor, Dr. Pierre Bézier, an engineer from the Renault car company who developed in the early 1960 s a curve formulation for use in shape design. The main interest of Bernstein-Bézier patches is that they lend to an easy geometric understanding of the underlying mathematical concepts. Some basic properties and a brief discussion of Bernstein polynomials and Bézier curves [START_REF] Bézier | Essai de définition numérique des courbes et des surfaces expérimentales: contribution à l'étude des propriétés des courbes et des surfaces paramétriques polynomiales à coefficients vectoriels[END_REF] are presented in the present chapter.

Bernstein basis

Bézier curves are expressed in terms of Bernstein polynomials which where introduced by Sergei Bernstein in order to formulate a constructive proof of the Weierstrass approximation theorem.

Definition 2.1.1. (Univariate Bernstein)

The Bernstein polynomials of degree p over the interval [0, 1] are defined explicitly by:

B k p (ζ) = C k p ζ k (1 -ζ) p-k ∀ k = 0, ..., p,
with the binomial coefficients C k p given by:

C k p =          p! k!(p-k)! if 0 ≤ k ≤ p, 0 otherwise.

BERNSTEIN BASIS

An example of constant, linear, quadratic and cubic Bernstein polynomials are presented in Fig. 2.1:

B k p (ζ) = B k 1 p 1 (ζ 1 ) ⊗ B k 2 p 2 (ζ 2 ) ⊗ ... ⊗ B k d p d (ζ d ).
The bivariate Bernstein polynomials are illustrated in Fig. 2.2 for the linear, quadratic and cubic cases.

Properties of the Bernstein polynomials

The Bernstein polynomials B k p of degree p, have several important properties [START_REF] Kelisky | Iterates of Bérnstein polynomials[END_REF]:

1.

Recursion

The Bernstein polynomials satisfy the following recursion relation: for p > 0, we have:

         B 0 p (ζ) = 1 -ζ p if k = 0, B k p (ζ) = 1 -ζ B k p-1 (ζ) + ζB k-1 p-1 (ζ) if k = 1, ..., p -1, B p p (ζ) = ζ p if k = p.
Proof.

∀1 ≤ k ≤ p -1, we have 1 -ζ B k p-1 (ζ) = 1 -ζ C k p-1 ζ k 1 -ζ (p-1)-k = C k p-1 ζ k 1 -ζ p-k = p -1 ! k! p -1 -k ! ζ k 1 -ζ p-k = p -k p B k p (ζ).
Similarly,

ξB k-1 p-1 (ζ) = ζC k-1 p-1 ζ k-1 1 -ζ (p-1)-(k-1) = C k-1 p-1 ζ k 1 -ζ p-k = p -1 ! k -1 ! p -1 -k -1 ! ζ k 1 -ζ p-k = k p B k p (ζ).
Therefore,

1 -ζ B k p-1 (ζ) + ζB k-1 p-1 (ζ) = p -k p B k p (ζ) + k p B k p (ζ) = B k p (ζ). ■ 2.

Non-negativity

The Bernstein polynomials are positive everywhere in [0, 1].

B k p (ζ) ≥ 0 ∀ζ ∈ [0, 1] 0 ≤ k ≤ p.
Note that the Lagrange interpolating polynomials, commonly used as basis for the numerical solution of partial differential equation (PDE), do not satisfy this property.

Partition of unity

The Bernstein polynomials of degree p forms a partition of unity, that is:

p k=0 B k p (ζ) = 1 ∀ζ ∈ [0, 1]. Proof. p k=0 B k p (ζ) = p k=0 C k p ζ k 1 -ζ p-k = ξ + (1 -ξ) p = 1. ■ 4.

Unique maximum

The Bernstein polynomial B k p has a unique maximum at ζ = k p on [0, 1].

Symmetry

The Bernstein polynomials are symmetric in the sens:

B k p (1 -ζ) = B p-k p (ζ) ∀0 ≤ k ≤ p.

Basis for polynomials of degree less or equal p

The Bernstein polynomials of degree p form a basis for the space of polynomials of degree less than or equal to p.

Degree elevation

Any Bernstein polynomial of degree less than p can be expressed as a linear combination of Bernstein polynomials of degree p. In particular, any Bernstein polynomial of degree p -1 can be written as a linear combination of Bernstein polynomials of degree p.

B k p-1 (ζ) = p -k p B k p (ζ) + k + 1 p B k+1 p (ζ). 2.3. DERIVATIVES Proof. p -k p B k p (ζ) + k + 1 p B k+1 p (ζ) = p -k p C k p ζ k 1 -ζ p-k + k + 1 p C k+1 p ζ k+1 1 -ζ p-(k+1) = C k p ζ k 1 -ζ p-k p -k p + k + 1 p p -k k + 1 ζ 1 -ζ = C k p ζ k 1 -ζ p-k p -k p 1 1 -ζ = (p -1)! k!(p -k -1)! ζ k 1 -ζ (p-1)-k = B k p-1 (ζ).

■

This property will play an important role, in the context of numerical solution of PDE, by allowing prefinement process.

Derivatives

The derivative of the k-th Bernstein polynomial of degree p is given by:

d d ζ B k p (ζ) = p B k-1 p-1 (ζ) -B k p-1 (ζ) .
Proof.

d d ζ B k p (ζ) = d d ζ C k p ζ k 1 -ζ p-k = C k p kζ k-1 1 -ζ p-k -p -k ζ k 1 -ζ p-k-1 = k ζ C k p ζ k 1 -ζ p-k - p -k 1 -ζ C k p ζ k 1 -ζ p-k = k ζ p k C k-1 p-1 ζ k 1 -ζ p-k - p -k 1 -ζ C k p ζ k 1 -ζ p-k = pC k-1 p-1 ζ k-1 1 -ζ p-k -p -k C k p ζ k 1 -ζ p-k-1 = pB k-1 p-1 (ζ) -p -k p! k!(p -k)! ζ k 1 -ζ p-k-1 = pB k-1 p-1 (ζ) -p (p -1)! k!(p -k -1)! ζ k 1 -ζ p-k-1 = pB k-1 p-1 (ζ) -pB k p-1 (ζ).

■

This formula is used to evaluate the derivatives of basis functions in the variational formulation, in a recursive way.

Bézier curves

Definition 2.4.1. Given (p + 1) distinct points P 0 , P 1 , ..., P p in space, the Bézier curve (of degree p) defined from these points is the parametric curve C p defined by:

C p (ζ) = p k=0 B k p (ζ)P k ∀ζ ∈ [0, 1].
The points P i i =0,...,p are called the control points and the line segments P 0 P 1 , P 1 P 2 , ..., P p-1 P p form in this order, the control polygon.

Properties

Let us recall some important properties of a Bézier curve [START_REF] Jaxon | Isogeometric analysis on triangulations[END_REF] [60], that will be useful for the numerical solving PDE systems.

Interpolation at the extremities

A Bézier curve C p of degree p, always starts at the first control point P 0 and ends at last control point This property is very important because, even if it is not necessary to fully control the curve in the middle, it is essential to know where it starts and where it ends. If we want to connect several Bézier curves, it is mandatory to know where the ends are. Moreover, this is important for applying Dirichlet boundary conditions when solving PDE systems.

P p .

Invariance under affine transformation

An affine transformation of a Bézier curve is obtained by applying the transformation to the control points.

Proof. Let ψ be an affine transformation in R p :

ψ(X ) = AX + b with A ∈ M p (R) and b ∈ R p .
The affine transformation of a Bézier curve C p of degree p is:

ψ(C p (ζ)) = ψ p k=0 B k p (ζ)P k = A p k=0 B k p (ζ)P k + b = p k=0 AB k p (ζ)P k + p k=0 bB k p (ζ) = p k=0 AP k + b B k p (ζ) = p k=0 ψ(P k )B k p (ζ). So ψ(C p (ζ)
) is a Bézier curve. ■

Convex hull

An important property of Bézier curves is that they always lies within the convex hull of their control points. To explain this property we need to define the convex hull of a set of points. The convex hull of a set P = {P 0 , P 1 , ..., P p } of control points is the smallest convex polygon that contains all the control points of P . In Fig. 2.5, the convex hull of the 4 control points is shown. Again, this property can be exploited when solving PDE, in particular in case of discontinuity capturing.

Local control

A Bézier curve is not interpolating the control points and moreover a Bézier curve changes globally when a control point is modified.

An example for a quadric Bézier curve is presented in Fig.

Degree elevation

The degree p of a Bézier curve C p may be increased without changing the curve geometrically or parametrically. The new p + 2 control points P 0 , P 1 , ..., P p , P p+1 are formed from the original p + 1 control points P 0 , P 1 , ..., P p by:

P k = k p + 1 P k-1 + p + 1 -k p + 1 P k ∀ 0 ≤ k ≤ p + 1.
Proof.

C p (ζ) = ζ + 1 -ζ C p (ζ) = ζC p (ζ) + 1 -ζ C p (ζ) = ζ p k=0 B k p P k + 1 -ζ p k=0 B k p P k = p k=0 k + 1 p + 1 B k+1 p+1 (ζ)P k + p k=0 p + 1 -k p + 1 B k p+1 (ζ)P k (see page 13) = p+1 k=1 k p + 1 B k p+1 (ζ)P k-1 + p k=0 p + 1 -k p + 1 B k p+1 (ζ)P k = p+1 k=0 k p + 1 B k p+1 (ζ)P k-1 + p+1 k=0 p + 1 -k p + 1 B k p+1 (ζ)P k = p+1 k=0 k p + 1 P k-1 + p + 1 -k p + 1 P k B k p+1 (ζ). ■
This process can be used to increase the order of the solution, in the context of PDE solving, without modifying the solution itself.

An example of degree elevation is depicted in Fig.

Derivatives of a Bézier Curve

The derivative of a Bézier curve of degree p, is another Bézier curve of degree p -1 given by:

d d ζ C p (ζ) = p p-1 k=0 B k p-1 (ζ) P k+1 -P k . (2.1)
More generally, derivatives of higher order are given by:

d r d ζ r C p (ζ) = p! (p -r )! p-r k=0 B k p-r (ζ) P k+1 -P k r ∀r ≤ p. (2.2)
This property has important consequences. In particular the curve at extremities is tangent to the first and last control points lines. This could be used to apply Neumann boundary conditions for instance.

Proof.

For the first derivative, we have:

d d ζ C p (ζ) = d d ζ p k=0 B k p (ζ)P k = p k=0 d B k p (ζ) d ζ P k = p k=0 p(B k-1 p-1 (ζ) -B k p-1 (ζ)) P k = p p k=0 B k-1 p-1 (ζ)P k -p p k=0 B k p-1 (ζ)P k note that B -1 p-1 (ζ) = B p p-1 (ζ) = 0 = p p-1 k=0 B k p-1 (ζ)P k+1 -p p-1 k=0 B k p-1 (ζ)P k = p p-1 k=0 B k p-1 (ζ) P k+1 -P k .
This can be written as:

d d ζ C p (ξ) = p p-1 k=0 B k p-1 (ζ)P k = p! (p -1)! p-1 k=0 B k p-1 (ζ)P k .
Now assume that (2.2) is true up to the order r . Let us prove it for r + 1:

d r +1 d ζ r +1 C p (ζ) = d d ζ d r d ζ r C p = d d ζ p! (p -r )! p-r k=0 B k p-r (ζ)(P k+1 -P k ) r = p! (p -r )! p-r k=0 d d ζ B k p-r (ζ)(P k+1 -P k ) r = p! (p -r )! p-r k=0 (p -r )(B k-1 p-r -1 (ζ) -B k p-r -1 (ζ)) P k+1 -P k r +1 = p! (p -r -1)! p-r k=0 B k-1 p-r -1 (ζ) -B k p-r -1 (ζ) P k+1 -P k r +1 = p! (p -r -1)! p-r k=1 (P k+1 -P k ) r +1 B k-1 p-r -1 (ζ) - p-r -1 k=0 (P k+1 -P k ) r +1 B k p-r -1 (ζ) = p! (p -r -1)! p-r -1 k=0 (P k+1 -P k ) r +1 B k p-r -1 (ζ) - p-r -1 k=0 (P k+1 -P k ) r +1 B k p-r -1 (ζ) = p! (p -(r + 1))! p-(r +1) k=0 (P k+1 -P k ) r +1 B k p-(r +1) (ζ) .
Therefore, the result is true for (r + 1). ■

Subdivision of Bézier curves

The De Casteljau algorithm is probably the most fundamental algorithm in the field of curve and surface design. This algorithm was devised in 1959 by Paul De Casteljau, a French mathematician from the Citroen automobile company. A main interest of this algorithm is to subdivide a Bézier curve C p of degree p, de-

fined in [ζ 1 , ζ l ], into two Bézier curves of degree p, C 1 p defined in [ζ 1 , ζ l 1 ] and C 2 p defined in [ζ l 1 , ζ l ]
whose union is equivalent to the original curve.

Consider now a Bézier curve of degree p defined by the control points P 0 , P 1 , ..., P p . The De Casteljau algorithm is used to obtain a point on the curve at a parameter value ζ ∈ [0, 1] from the control polygon composed of the points P i . We construct P j i : the new control point during the subdivision step j . Formally, the algorithm of De Casteljau can be written as [START_REF] Nava-Yazdani | De Casteljaus algorithm on manifolds[END_REF]: Thus, it is possible to subdivide a Bézier curve of degree p, into several curves which describe the same curve. Because the resulting Bézier curves must have their own new control points, the original set of control points is discarded. Moreover, since the original Bézier curve C p is cut into several pieces, each of which is a subset of the original degree p Bézier curve, the resulting Bézier curves must be of degree p.

   P 0 i = P i i = 0, ..., p -r, P j i = (1 -ζ)P j -1 i + ζP j -1 i +1 j = 1, ..., p i = 0, ..., p -j.

RATIONAL BÉZIER CURVES

Note that this algorithm is used in practice to evaluate the curve at a given parameter ζ (for visualization for example).

Rational Bézier curves

Given (p +1) control points P 0 , P 1 , ..., P p and associated nonnegative weights w 0 , w 1 , ..., w p we can define, for ζ ∈ [0, 1], the rational Bézier curve of degree p as:

R p (ζ) = p i =0 w i B i p (ζ)P i p i =0 w i B i p (ζ)
.

Note that if all weights are 1 (or if all weights are simply equal), a rational Bézier curve reduces to a polynomial Bézier curve. However, almost all properties of Bézier curves still hold for rational Bézier curves (details can be found in [START_REF] Floater | Derivatives of rational Bézier curves[END_REF]). In particular, rational Bézier curves have end-point interpolation, prescribed tangent lines at the endpoints. They also satisfy the convex hull property. Moreover, degree elevation and subdivision can be extended to rational curves. Further, rational Bézier curves have several advantages over polynomial Bézier curves. They provide more control over the shape than does a polynomial Bézier curve.

Elsewhere, it is possible to reparametrize the curve by simply changing the weights in a specific manner.

Moreover, rational Bézier curves are needed to exactly represent conic sections [START_REF] Lee | Approximate conversion of rational Bézier curves[END_REF]. Therefore, their use is fundamental in Computer-Aided-Design.

Bézier surface

One can easily extend the previous concepts to define Bézier surfaces, by using a second-order tensor product of Bernstein polynomials. More precisely, a Bézier surface of degree p 1 × p 2 can be formulated as:

S(ζ 1 , ζ 2 ) = p 1 i =0 p 2 j =0 B i p 1 (ζ 1 )B j p 2 (ζ 2 )P i j ∀ζ 1 , ζ 2 ∈ [0, 1],
where the P i j define a control net. The degrees of Bernstein polynomials p 1 and p 2 do not necessarily have to be the same in the two parameter directions. The corresponding properties of the Bézier curve apply to the Bézier surface:

• Invariance under affine transformations.

• Subdivision produces smooth continuous surfaces.

• The Bézier surface is contained within the convex hull of the control nets.

• The Bézier surface can be exactly represented as a Bézier surface of higher degree.

• The Bézier surface does not in general pass through the control nets except for the corners of the control net grid.

• Closed surfaces can be formed by setting the last control net equal to the first.

CONCLUSION

L Es courbes de Bézier ont été inventées dans les années 1960 par l'ingénieur Pierre Bézier. Il a étudié le problème de conception de surfaces 3D (carrosseries d'automobiles, fuselages d'avion, etc.) pour les premiers programmes de CAO (Conception Assistée par Ordinateur). Le but était de trouver un moyen pour définir des courbes paramétriques de manière précise et simple.

C H A P T E R

B-SPLINES CURVES

B -Spline and Non-Uniform-Rational B-spline (NURBS) (a generalization of B-spline) are the de facto standard for geometric representation in computer aided design (CAD) systems [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF]. A brief discussion of B-spline and NURBS is presented in this chapter. We begin by introducing necessary background concepts, defining commonly used notations and introducing B-spline/NURBS curves and surfaces [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF] [61] [START_REF] Régis | An introduction to isogeometric analysis with application to thermal conduction[END_REF].

B-spline functions

Univariate B-spline functions are piecewise polynomial functions with compact support. They are defined in parametric space using a so-called knot vector denoted Ξ, in one-dimensional space (1D), Ξ is a set of m non-decreasing coordinates:

Ξ = ξ 1 , ξ 2 , ..., ξ m . Definition 3.1.1.
The univariate B-spline function N i ,p of degree p is defined according to the Cox-de Boor recursion formula [START_REF] Boor | On calculating with B-splines[END_REF]:

For p = 0 :

N i ,0 (ξ) =    1 if ξ i ≤ ξ < ξ i +1 ∀i = 1, ..., m -1, 0 otherwise. (3.1)
For p ≥ 1: We use the convention that fraction in front of the basis functions is set equal to zero in the case of the denominator being zero. That is:

N i ,p (ξ) = ξ -ξ i ξ i +p -ξ i N i ,p-1 (ξ) + ξ i +p+1 -ξ ξ i +p+1 -ξ i +1 N i +1,p-1 (ξ). ( 3
ξ -ξ i ξ i +p -ξ i ≡ 0 if ξ i +p -ξ i = 0, ξ i +p+1 -ξ ξ i +p+1 -ξ i +1 ≡ 0 if ξ i +p+1 -ξ i +1 = 0.
We see that there are exactly n basis functions. This is a direct consequence of the knot vector consisting of n +p +1 knots. If we increase the number of elements in the knot vector, we will also increase the number of basis functions. Note that the degree and the number of basis functions are independent, contrary to Bézier functions.

With Eq. (3.1) and Eq. (3.2), it can be noted that for p = 0 and p = 1, the basis functions of isogeometric analysis are identical to those of the standard piecewise constant and linear finite elements, respectively. In fact, the greatest wealth of B-splines comes from the case p = 2.

In order to define multivariate B-splines functions in higher dimensions, we make use of the tensor product.

Definition 3.1.2. Let p = (p 1 , p 2 , ..., p d ) be a vector in N d and for all j = 1, ..., d , let Ξ j be a 1D knot vector defined by:

Ξ j = ξ j 1 , ξ j 2 , ..., ξ j n 1 +p 1 +1 .
Furthermore, if we denote the i j univariate B-spline of degree p j defined on the knot vector Ξ j by N i j ,p j (ξ j ), then, with the multi-indices i = (i 1 , i 2 , ..., i d ), p = (p 1 , p 2 , ..., p d ) and n = (n 1 , n 2 , ..., n d ), the d-dimensional tensor product B-spline is defined by:

N i ,p (ξ) = N i 1 ,p 1 (ξ 1 ) ⊗ N i 2 ,p 2 (ξ 2 ) ⊗ ... ⊗ N i d ,p d (ξ d ).

Knot Vectors

A knot vector Ξ in one-dimensional space is a set of finite, real-valued, monotonically increasing sequence of real numbers written,

Ξ = ξ 1 , ξ 2 , ..., ξ n+p+1 ,
where

ξ i ∈ R is the i -t h knot, i ∈ {1, 2, ..., n + p + 1}
is the knot index, p is the polynomial degree and n is the number of basis functions which define the B-spline. The interval [ξ 1 , ...,

ξ i ] is called a patch ∀1 < i ≤ n + p + 1, whereas the interval between two knots [ξ i , ξ i +1 ] is called knot span. A knot span is called empty if ξ i = ξ i +1 and is called interior if ξ i < ξ i +1 ∀1 ≤ i ≤ n + p.
There are different types of knot vectors:

• Periodic knot vector: a knot vector Ξ is periodic if there exists an integer I and a real T such that for every i ∈ Z,

ξ i +I = ξ i + T.
• Uniform knot vector: if knots are equally-spaced in the parametric space, they are said to be uniform and non-uniform otherwise.

• Open knot vetor: a knot vector for a B-spline basis of degree p is said to be open if the first and last knots are repeated p + 1 times.

Ξ = ξ 1 , ..., ξ 1 (p+1)times , ξ p+2 , ..., ξ n+p+1 , ..., ξ n+p+1 (p+1)times .
Open knot vectors are standard in the CAD literature. In one dimension, basis functions formed from open knot vectors are interpolatory at the ends of the parametric space interval, [ξ 1 , ξ n+p+1 ], and at the corners of patches in multiple dimension but they are not, in general, interpolatory at interior knots. That is a distinguishing feature between "knots" and "nodes" in FEA. Note that the basis functions are interpolatory at the ends of the interval and also at ξ = 1, the location of a repeated knot, where only C 0 -continuity is attained (see Fig. 3.4). This property will be used later, to generate a computational domain suitable to Discontinuous Galerkin methods.

Properties of the B-spline functions

The B-spline functions N i ,p defined using the procedure described above are polynomials of degree p.

There are several important features of the basis functions that are pointed out by Hughes et al. [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF], in the perspective of PDE analysis.

1. B-spline functions form a partition of unity, i.e.

n i =1 N i ,p (ξ) = 1 ∀ξ ∈ Ξ = ξ 1 , ξ 2 , ..., ξ n+p+1 .
2. The basis functions are interpolatory at the end points of the knot vector, such that:

N i ,p (ξ 1 ) = δ i ,1 and N i ,p (ξ n+p+1 ) = δ i ,n
where, δ is the Kronecker symbol.

DERIVATIVES OF B-SPLINE FUNCTIONS

3. The support of each N i ,p is compact and contained in [ξ i , ξ i +p+1 ]. As seen in Fig. 3.4, the supports of the functions are growing with increasing polynomial degree. In fact, the support will cover exactly p + 2 knots. Note however that some of these knots may be equal and thus have knot multiplicity greater than one. In these cases, the support will not go over p + 1 knot spans, which is the maximum support it can have.

4. Each B-spline function is nonnegative over the entire domain, that is:

N i ,p (ξ) ≥ 0 ∀ξ ∈ Ξ = ξ 1 , ξ 2 , ..., ξ n+p+1 ,
which means that all the coefficients of a mass matrix computed from a B-spline basis are also nonnegative, which can be useful for mass lumping schemes [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF].

5. The B-spline functions indeed form a basis for the space of polynomials of degree less than or equal to p, P p . That is, they are all linearly independent i.e.

n i =1 α i N i ,p (ξ) = 0 ⇐⇒ α i = 0 ∀i = 1, 2, ..., n.
6. The B-spline function N i ,p (ξ i ) is of regularity C p-r at each knot of multiplicity r . When the multiplicity of a knot is exactly p, the basis function is interpolatory.

Derivatives of B-spline functions

The derivatives of B-spline functions are represented in terms of lower order B-spline basis, as it comes directly from the recursive definition given in equations (3.1) and (3.2). Thus, the first derivative of the i -th B-spline basis function of degree p is given by:

d d ξ N i ,p (ξ) = p 1 ξ i +p -ξ i N i ,p-1 (ξ) - 1 ξ i +p+1 -ξ i +1 N i +1,p-1 (ξ) . (3.3) Proof. d d ξ N i ,0 (ξ) = 0 N i ,0 (ξ) and N i +1,0 (ξ) are constants . d d ξ N i ,1 (ξ) = d d ξ ξ -ξ i ξ i +1 -ξ i N i ,0 (ξ) + ξ i +2 -ξ ξ i +2 -ξ i +1 N i +1,0 = 1 ξ i +1 -ξ i N i ,0 (ξ) - 1 ξ i +2 -ξ i +1 N i +1,0 .
We assume that for all 0 ≤ q ≤ p, we have:

d N i ,q (ξ) d ξ = q ξ i +q -ξ i N i ,q-1 (ξ) - q ξ i +q+1 -ξ i +1 N i +1,q-1 (ξ).
Let us prove that (3.3) is true for p + 1.

d d ξ N i ,p+1 (ξ) = d d ξ ξ -ξ i ξ i +p+1 -ξ i N i ,p (ξ) + ξ i +p+2 -ξ ξ i +p+2 -ξ i +1 N i +1,p (ξ) = 1 ξ i +p+1 -ξ i N i ,p (ξ) + ξ -ξ i ξ i +p+1 -ξ i d N i ,p (ξ) d ξ - 1 ξ i +p+2 -ξ i +1 N i +1,p (ξ)
+ ξ i +p+2 -ξ ξ i +p+2 -ξ i +1 d N i +1,p (ξ) d ξ = 1 ξ i +p+1 -ξ i N i ,p (ξ) + ξ -ξ i ξ i +p+1 -ξ i p ξ i +p -ξ i N i ,p-1 (ξ) - p ξ i +p+1 -ξ i +1 N i +1,p-1 (ξ) - 1 ξ i +p+2 -ξ i +1 N i +1,p (ξ) + ξ i +p+2 -ξ ξ i +p+2 -ξ i +1 p ξ i +p+1 -ξ i +1 N i +1,p-1 (ξ) - p ξ i +p+2 -ξ i +2 N i +2,p-1 (ξ) = 1 ξ i +p+1 -ξ i N i ,p (ξ) + p ξ i +p -ξ i ξ -ξ i ξ i +p+1 -ξ i N i ,p-1 (ξ) - 1 ξ i +p+2 -ξ i +1 N i +1,p (ξ)
+ -p ξ i +p+1 -ξ i +1 ξ -ξ i ξ i +p+1 -ξ i + p ξ i +p+1 -ξ i +1 ξ i +p+2 -ξ ξ i +p+2 -ξ i +1 N i +1,p-1 (ξ) + -p ξ i +p+2 -ξ i +2 ξ i +p+2 -ξ ξ i +p+2 -ξ i +1 N i +2,p-1 (ξ). p ξ i +p+1 -ξ i +1 ξ i +p+2 -ξ ξ i +p+2 -ξ i +1 - ξ -ξ i ξ i +p+1 -ξ i N i +1,p-1 = p ξ i +p+1 -ξ i +1 ξ i +p+2 -ξ ξ i +p+2 -ξ i +1 - ξ i +p+2 -ξ i +1 ξ i +p+2 -ξ i +1 N i +1,p-1 + p ξ i +p+1 -ξ i +1 ξ i +p+1 -ξ i ξ i +p+1 -ξ i - ξ -ξ i ξ i +p+1 -ξ i N i +1,p-1 = p ξ i +p+1 -ξ i +1 ξ i +1 -ξ ξ i +p+2 -ξ i +1 + ξ i +p+1 -ξ ξ i +p+1 -ξ i N i +1,p-1 .

B-SPLINE CURVES

Then, we get:

d d ξ N i ,p+1 (ξ) = 1 ξ i +p+1 -ξ i N i ,p (ξ) - 1 ξ i +p+2 -ξ i +1 N i +1,p (ξ) + p ξ i +p+1 -ξ i ξ -ξ i ξ i +p -ξ i N i ,p-1 (ξ) + ξ i +p+1 -ξ ξ i +p+1 -ξ i +1 N i +1,p-1 (ξ) + p ξ i +p+1 -ξ i +1 ξ i +1 -ξ ξ i +p+2 -ξ i +1 N i +1,p-1 (ξ) + -p ξ i +p+2 -ξ i +2 ξ i +p+2 -ξ ξ i +p+2 -ξ i +1 N i +2,p-1 (ξ) = p + 1 ξ i +p+1 -ξ i N i ,p (ξ) - 1 ξ i +p+2 -ξ i +1 N i +1,p (ξ) + -p ξ i +p+2 -ξ i +1 ξ -ξ i +1 ξ i +p+1 -ξ i +1 N i +1,p-1 (ξ) + ξ i +p+2 -ξ ξ i +p+2 -ξ i +2 N i +2,p-1 (ξ) = p + 1 ξ i +p+1 -ξ i N i ,p (ξ) - p + 1 ξ i +p+2 -ξ i +1 N i +1,p (ξ).
■

B-spline curves

B-spline curves are defined as a linear combination of control points and B-spline basis functions.

Given n basis functions N i ,p , i = 1, ..., n and corresponding control points P i ∈ R, i = 1, ..., n, a piecewisepolynomial B-spline curve is obtained as:

C p (ξ) = n i =1 N i ,p (ξ)P i .
The linear interpolation of the control points is called the control polygon. It can be seen that a Bézier curve of order n + 1 (degree n) is a special case of B-spline curve with no internal knots and the ends knots are repeated n + 1 times. The knot vector is thus:

Ξ = 0, 0, ..., 0 (n+1)t i mes , 1, 1, ..., 1 (n+1)t i mes .
An example is shown in Fig. 3.5 for the quadratic basis functions considered previously. Note that the curve is interpolatory at the first and last control points P 1 and P 5 , due to the fact that the knot vector is uniform and open, Ξ = 0, 0, 0, 1, 2, 3, 3, 3 . The curve C 2 is tangent to the control polygon at the first and last control points. The curve is C 1 continuous everywhere. B-spline curves are generalization of Bézier curves and share many important properties with them. Moreover, B-spline curves have additional useful properties [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF]. We list below some of the most important properties of B-spline curves:

• The curve C p is C p-1 continuous everywhere except at the knot or control point of multiplicity r , where it is C p-r .

• Invariance with respect to affine transformations. In fact, an affine transformation of a B-spline curve is obtained by applying the transformation directly to the control points.

• If the knot vector is open, the curve starts at the point P 1 and ends at P n . Moreover, it is tangent to (P 0 P 1 )

and (P n-1 P n ) at the extremities.

• Changing the position of control point P i affects C p only in the interval [ξ i , ξ i +p ], which is a significant difference with respect to Bézier curves.

HIERARCHICAL REPRESENTATION

• B-spline curve is contained in the convex hull of its control point. More specifically, if ξ ∈ [ξ i , ξ i +1 ), then C p (ξ) is in the convex hull of control points P i -p , P i -p+1 , ..., P i .

• The derivative of a B-spline curve is given by:

d d ξ C p (ξ) = p n-1 i =0 P i +1 -P i ξ i +p+1 -ξ i +1 N i +1,p-1 (ξ).

Hierarchical representation

The B-spline curves can be enriched by three types of refinements, which are termed h-, p-and krefinements, without changing the shape of the geometry. For further details, we refer the reader to [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF].

Knot insertion

The first mechanism by which one can enrich the basis is knot insertion. Knots may be inserted without changing a curve geometrically or parametrically. More precisely, given a B-spline curve with m control points P i , the same curve can be obtained with m + 1 control points P i by inserting ξ ∈ [ξ k , ξ k+1 [ given by:

P i = α i P i + (1 -α i )P i -1 ,
where,

α i =          1 if i ≤ k -p, ξ-ξ i ξ i +p -ξ i if k -p + 1 ≤ i ≤ k, 0 if i ≥ k + 1.
An example of knot insertion is presented in Fig. 3.7. The original curve, shown on the left, consists of a cubic B-spline curve with 9 control points and has 9 + 3 + 1 = 13 knot values. The new curve, shown on the right, is geometrically and parametrically identical to the first one, but the basis functions and control points are changed, there is one more of each. This process may be repeated to enrich the solution space by adding more basis functions of the same order while keeping the curve unchanged. Insertion of new knot values clearly has similarities with the classical h-refinement strategy in FEA. As explained in section 3.1.1, inserting a knot on an existing knot results in a decrease of the curve regularity (multiple knots).

Order elevation

The second mechanism by which one can enrich the basis is order elevation (sometimes also called "degree elevation") can be thought of as p-refinement in FEA. During this process, the multiplicity of each knot is increased by one but no new knots are added. Let Ξ = {ξ 1 , ..., ξ m } be a knot vector, the degree p of a B-spline curve C p defined for the knot vector Ξ may be increased without changing the geometry or parametrization. It is possible to define another B-spline curve of degree p +1 that is identical to the original one. The number of new control points depends on the multiplicities of existing knots. The process for order elevation is the following:

1. Begins by subdividing the curve into many Bézier curves of degree p by knot insertion (see section 3.7).

2. The next step is to elevate the order of the polynomial on each of these individual segments (see section 2.5.1).

3. Last, excess knots are removed to combine the segments into a B-spline curve of degree p + 1.

The basic idea of degree raising and knot insertion is to achieve the flexibility without changing the shape of the curve or surface. We refer to [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF] for the mathematical details.

k-refinement

k-refinement refers to the process in which order elevation is followed by knot insertion. It has no analogous in FEA. It is important to point out that the order elevation and knot insertion do not commute.

This process results in a higher order and higher continuity basis than the process of knot insertion followed by order elevation. We refer to [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF] for a thorough treatment and application examples.

B-spline surfaces and volumes

Given a control net (P i , j ) i =1,...,n j =1,...,m and two knot vectors, Ξ 1 = {ξ 1 , ..., ξ n+p+1 } and Ξ 2 = {η 1 , ..., η m+q+1 }, a tensor-product B-spline surface is defined as:

S (ξ, η) = n i =1 m j =1
N i ,p (ξ)N j ,q (η)P i , j .

V (ξ, η, ζ) = n i =1 m j =1 l k=1 N i ,p (ξ)N j ,q (η)N k,r (ζ)P i , j ,k .
It is important to note that many properties of a B-spline surface and volume are the results of the tensor product nature. Multivariate B-splines basis functions are nonnegative, have local support, are invariant to affine transformations and form a partition of unity.

Non-Uniform Rational B-spline (NURBS)

B-splines are convenient for free-form modelling, but they lack the ability to exactly represent some simple shapes such as circles and ellipsoids. NURBS are non-rational functions of B-splines, they allow for the exact parametrization of common curves and surfaces such as circles, cylinders and spheres. They are also extremely flexible and intuitive when dealing with more complex shape creation and deformation. A NURBS entity in R d is obtained by the projective transformation of a B-spline entity in R d +1 [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF].

NURBS basis functions

Let (w i ) 1≤i ≤n be a sequence of non-negative reals (weights for control point). The i -th NURBS function of degree p, associated to the knot vector Ξ and the weights w, is given by:

R p i (ξ) = N i ,p (ξ)w i n j =1 N j ,p (ξ)w j ,
where N i ,p (ξ) denotes the i -th B-spline basis function of degree p.

The first derivative of a NURBS basis function is given by:

d d ξ R p i (ξ) = w i d d ξ N i ,p (ξ)W (ξ) -N i ,p (ξ) d d ξ W (ξ) W 2 (ξ) ,
where,

W (ξ) = n j =1 N j ,p (ξ)w j , d d ξ W (ξ) = n j =1 d d ξ N j ,p (ξ)w j .

NURBS curves and surfaces

The NURBS curve of degree p associated to the knot vector Ξ, the control points (B i ) 1≤i ≤n and the weights (w i ) 1≤i ≤n is defined as:

C p (ξ) = n i =1 R p i (ξ)B i .
Similarly, the NURBS surfaces of degree p 1 × p 2 associated to the knot vectors Ξ 1 and Ξ 2 , the control nets (B i , j ) 1≤i ≤n 1 ,1≤ j ≤n 2 and the weights (w i j ) 1≤i ≤n 1 ,1≤ j ≤n 2 , are defined by:

S(ξ,η) = n 1 i =1 n 2 j =1 R i j (ξ, η)B i , j .
With,

R i j (ξ, η) = w i j N i ,p 1 (ξ)N j ,p 2 (η) n 1 i =1 n 2 j =1 w i j N i ,p 1 (ξ)N j ,p 2 (η)

EXTRACTING BÉZIER CURVES FROM B-SPLINES

We recall some important properties of NURBS (we refer to [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF] for more details):

• NURBS basis functions forms a partition of unity.

• The continuity and support of NURBS basis functions are the same as for B-splines.

• Affine transformations in physical space are obtained by applying the transformation to the control points, that is, NURBS possess the property of affine invariance.

• If the weights are equal, NURBS become B-splines.

• NURBS surfaces and solids are the projective transformations of tensor product, piecewise polynomial entities.

Extracting Bézier curves from B-splines

Let us consider the knot vector Ξ = ξ 1 , ..., ξ n+p+1 where, ξ 1 = ... = ξ p+1 and ξ n+1 = ξ n+p+1 .

To decompose a B-spline (or NURBS) curve to its Bézier elements, called Bézier extraction, a straightforward approach consists in using the knot-insertion procedure several times, for each of the existing interior knots (ξ p+2 , ..., ξ n ) until interior knots have a multiplicity p + 1. Then, the original B-spline curve is separated in independant Bézier elements. The curve is geometrically unmodified, but its representation is split in a set of discontinuous elements. This is the key procedure to build a computational domain suitable to Discontinuous Galerkin method.

It is important to point out that the Bézier patch is a particular case of B-spline patch, for which the number n of functions (and control points) is equal to p +1. An example illustrating Bézier extraction from quadratic B-spline is shown in Fig. 3.9.

CONCLUSION

L Es courbes B-splines ont été définies dans les années C Urve or surface fitting, also known as regression analysis, is used to approximate a curve or surface for a series of data points. Most of the time, the curve/surface fitting is used to find points and derivatives anywhere along the curve/surface. This procedure is attractive because of its ability to produce an interpolating curve/surface that retains extremely high accuracy with a minimal number of data to represent the curve/surface, so that the error is minimized in the least squares sense. The purpose of this chapter is to present the problem of fitting a given ordered set of data with a B-spline curve or surface in the least squares sense.

In the context of isogeometric analysis, this technique will be used in the following two cases:

• to define the B-spline computational domain if the geometry of the problem is not defined originally in terms of B-spline,

• to define the B-spline initial conditions.

Curve fitting 4.1.1 Basic concepts

Curve or surface fitting is a fundamental problem in many fields, such as computer graphics, image processing, shape modeling and data mining. Depending on applications, different types of curves such as parametric curves, implicit curves and subdivision curves are used for fitting [53] [96]. In this chapter, we discuss the problem of B-spline curve and surface fitting [10] [57]. Due to the fact that the curve and surface fitting are closely related topics, we define and discuss the curve fitting problem first and show how to generalize it to surface reconstructions later on.

CURVE FITTING

The general objective of curve fitting is to theoretically describe experimental data with a model (function or equation) and to find the parameters associated with this model.

We recall that a B-spline curve of degree p is defined for a collection of n control points (P i ) i =1,...,n by:

C p (ξ) = n i =1 N i ,p (ξ)P i .
The functions N i ,p (ξ) are the B-spline basis functions of degree p, which are defined recursively and require the selection of a sequence of scalars The basic aim is to fit a B-spline curve that will approximate n ev al measured data points in a least squares sense (with n ev al n). This leads to find an optimal set of control point {P i } 1≤i ≤n producing an optimal approximating B-spline C p with minimal distances to the points Q k . We seek therefore to minimize the least-squares error, defined as the sum of squares of the points distance expressed as:

ξ i ∈ Ξ = ξ 1 ,
1 2 n ev al k=1 n i =1 N i ,p (ζ k )P i -Q k 2 ,
where

n i =1 N i ,p (ζ k )P i is the B-spline curve point at ζ k and Q k is the corresponding measured data point.

Description of the least squares method

The least-squares estimation procedure is a mathematical tool that was developed independently by Carl Friedrich Gauss in 1795 and Adrien-Marie Legendre who published it first in 1805 [START_REF] Krarup | A contribution to the mathematical foundation of physical geodesy[END_REF]. This theory describes a frequently used approach evolved from statistical methods to estimate values of parameters of a mathematical model from measured data, which are subject to errors.

For a specified set of control points, the least-squares error function between the B-spline curve and sample points is the scalar-valued function:

E (P ) = 1 2 n ev al k=1 n j =1 N j ,p (ζ k )P j -Q k 2 . ( 4.1)
As it has been pointed previously, the goal is to find values of the control points that minimize the error.

Therefore, we find the values of {P i } 1≤i ≤n in such a way that:

∂E ∂P i = 0.
Thus, carrying out the differentiation leads to:

∂E ∂P i = ∂ ∂P i 1 2 n ev al k=1 n j =1 N j ,p (ζ k )P j -Q k 2 = n ev al k=1 n j =1 N j ,p (ζ k )P j -Q k N i ,p (ζ k ) = n ev al k=1 n j =1 N i ,p (ζ k )N j ,p (ζ k )P j - n ev al k=1 N i ,p (ζ k )Q k = n ev al k=1 n j =1
a ki a k j P j -

n ev al k=1 a ki Q k ,
where

a ki = N i ,p (ζ k ).
Setting the partial derivatives equal to the zero vector (0 R n ) leads to the system of equations:

0 R n = n ev al k=1 n j =1
a ki a k j P j -

n ev al k=1 a ki Q k = A t AP -A t Q,
where a matrix A ∈ R n ev al ×n given by: A = a ki 1≤k≤n ev al ,1≤i ≤n .

This system of equations is a least-squares formulation for an over determined problem.

Hence, the minimization leads to the linear system:

A t AP = A t Q. (4.2)
As consequence, the equation 4.2 is equivalent to solving:

P = A t A -1 A t Q. (4.3)

B-spline curve fitting -example

In this section, an example of Least-squares B-spline curve fitting is given to illustrate the performance of the method.

Consider the sinusoidal function:

f (x) = sin(2πx) ∀x ∈ Ω = [-1, 1].
In our example we consider fitting a data set of n ev al uniformly distributed parameters (ζ k ) 1≤k≤n ev al and

points Q k = sin(2πζ k ), ∀1 ≤ k ≤ n ev al .
In this case, we choose n ev al = 80 n = 8. As seen in Fig. 4.1 an accurate approximation is obtained using only 8 control points.

LEAST-SQUARES B-SPLINE SURFACE FITTING

Least-squares B-spline surface fitting

A B-spline tensor product surface is defined for a bi-dimensional array of n × m control points P i j with

1 ≤ i ≤ n and 1 ≤ j ≤ m, S (ξ, η) = n i =1 m j =1 N i ,p (ξ)N j ,q (η)P i , j . (4.4)
The numbers p and q are the degrees for the surface, N i ,p (ξ) and N j ,q (η) are the B-spline basis functions.

The sample data points are (

ζ k 1 , θ k 2 ,Q k 1 k 2 ) with 1 ≤ k 1 ≤ n ev al and 1 ≤ k 2 ≤ m ev al .
It is assumed that:

ζ 1 < ζ 2 < • • • < ζ n ev al and θ 1 < θ 2 < • • • < θ m ev al .
The control points may be arranged formally as an n × m matrix.

P =      P 11 • • • P 1m . . . . . . . . . P n1 • • • P nm      .
Similarly, the samples Q k 1 k 2 may be arranged formally as an n ev al × m ev al matrix:

Q =      Q 11 • • • Q 1m ev al . . . . . . . . . P n ev al 1 • • • P n ev al m ev al      .
For a specified set of control points, the least-squares error function between the B-spline surface and sample points is the scalar-valued function:

E (P ) = 1 2 n ev al k 1 =1 m ev al k 2 =1 n j 1 =1 m j 2 =1 N j 1 ,p (ζ k 1 )N j 2 ,q (θ k 2 )P j 1 j 2 -Q k 1 k 2 2 . (4.5)
We determine as previously the control points by minimizing the error function.

The first-order partial derivatives are written in terms of the control points P j 1 j 2 rather than in terms of the components of the control points:

∂E ∂P i 1 i 2 = ∂ ∂P i 1 i 2 1 2 n ev al k 1 =1 m ev al k 2 =1 n j 1 =1 m j 2 =1 N j 1 ,p (ζ k 1 )N j 2 ,q (θ k 2 )P j 1 j 2 -Q k 1 k 2 2 = n ev al k 1 =1 m ev al k 2 =1 n j 1 =1 m j 2 =1 N j 1 ,p (ζ k 1 )N j 2 ,q (θ k 2 )P j 1 j 2 -Q k 1 k 2 N i 1 ,p (ζ k 1 )N i 2 ,q (θ k 2 ) = n ev al k 1 =1 m ev al k 2 =1 n j 1 =1 m j 2 =1 N j 1 ,p (ζ k 1 )N j 2 ,q (θ k 2 )N i 1 ,p (ζ k 1 )N i 2 ,q (θ k 2 )P j 1 j 2 - n ev al k 1 =1 m ev al k 2 =1 N i 1 ,p (ζ k 1 )N i 2 ,q (θ k 2 )Q k 1 k 2 = n ev al k 1 =1 m ev al k 2 =1 n j 1 =1 m j 2 =1 a k 1 j 1 b k 2 j 2 a k 1 i 1 b k 2 i 2 P j 1 j 2 - n ev al k 1 =1 m ev al k 2 =1 a k 1 i 1 b k 2 i 2 Q k 1 k 2 ,
where a k 1 j 1 = N j 1 ,p (ζ k 1 ) and b k 2 j 2 = N j 2 ,q (θ k 2 ). These equations may be written in matrix notation as:

∂E ∂P = A t AP B t B -A t QB, ( 4.6)
where, A = a k j 1≤k≤n ev al ,1≤ j ≤n is a n ev al × n matrix, and B = b k j 1≤k≤m ev al ,1≤ j ≤m is a m ev al × m matrix.

Setting the partial derivatives equal to the zero matrix 0 R n×m leads to the matrix system of equations,

A t AP B t B -A t QB = 0 R n×m . (4.7)
Then, we get: cases will be used in forthcoming chapters for isogeometric analysis applications. Note that the leastsquares fitting method presented previously is applied patch by patch when the computational domain is composed by a set of patches.

P = A t A -1 A t Q B B t B -1 = A t A -1 A t Q B t B -1 B t t .

B-spline surface fitting -examples

The first example is the two-dimensional Gaussian function.

f (x, y) = a exp - (x -x 0 ) 2 2σ 2 x + (y -y 0 ) 2 2σ 2 y ∀(x, y) ∈ Ω = [-1, 1] × [-1, 1].
Here the coefficient a is the amplitude, x 0 , y 0 is the center and σ 2 x , σ 2 y are the x and y spreads of the blob. Figure 4.2 created using a = 1, (x 0 , y 0 ) = (0, 0), and

σ x = σ y = 1/ 2.
The physical domain Ω is a set of quadratic Bézier patch with 5 × 5 patches. The mesh type considered is a simple Cartesian grid. In each patch, we define local data points ( Bessel equation (named according to the astronomer Friedrich Wilhelm Bessel) is a second-order differential equation with two linearly independent solutions: a Bessel function of the first kind of order ν, J ν and a Bessel function of the second kind of order ν, Y ν , where ν is a non-negative real number. We refer the reader to [START_REF] Dattoli | Theory and applications of generalized Bessel functions[END_REF].

ζ k 1 , θ k 2 ,Q k 1 k 2 ) ∀1 ≤ k 1 ≤ n ev al , 1 ≤ k 2 ≤ m ev
The Bessel function of the first kind J ν (x) can be written as an infinite polynomial:

J ν (x) = ∞ i =0 (-1) i i !Γ(ν + i + 1) x 2 ν+2i ,
where Γ is the gamma function, satisfy: Γ(i ) = (i -1)! for i a positive integer.

Some Bessel functions are plotted in Fig. 4.3.

CONCLUSION

L 'ajustement des courbes et des surfaces est un problème fondamental dans de nom- breux domaines, tels que l'infographie, le traitement d'images, la modélisation de formes et l'exploration de données. Selon les applications, différents types de courbes, telles que les courbes paramétriques, les courbes implicites et les courbes de subdivision, sont utilisées pour ajuster. Dans ce chapitre, nous avons discuté le problème de l'ajustement de courbes et surfaces B-splines en se basant sur l'approximation par moindres carrés.

La méthode d'approximation par moindres carrés est un outil mathématique, développé indépendamment par Carl Friedrich Gauss en 1795 et Adrien-Marie Legendre qui l'a publié en 1805. Cette théorie décrit une approche fréquemment utilisée à partir de méthodes statistiques pour estimer les valeurs des paramètres d'un modèle mathématique à partir de données mesurées, qui sont sujettes à des erreurs.

Dans le cadre de l'analyse isogéométrique, cette technique est utilisée dans deux cas suivants:

• pour définir le domaine de calcul B-spline si la géométrie du problème n'est pas définie à l'origine en termes de B-spline,

• pour définir les conditions initiales dans l'espace de représentation B-spline.

Part II ISOGEOMETRIC ANALYSIS -FINITE ELEMENT FRAMEWORK (ILLUSTRATION FOR A 1D PROBLEM)

C H A P T E R

5

SUPG -FINITE ELEMENT METHOD

T He aim of IGA is to generalize and improve upon classical FEA. In this chapter, we start by giving an introduction to IGA by revisiting the original analysis, i.e. FEA, in the context of hyperbolic PDEs.

The stabilized method Streamline Upwind/Petrov Galerkin (SUPG) [START_REF] Brooks | Streamline Upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF] [48] [START_REF] Mizukami | A Petrov-Galerkin finite element method for convection-dominated flows: an accurate upwinding technique for satisfying the maximum principle[END_REF] had its origin in the late 1970s and early 1980s. In the present chapter, the capability of this stabilized Finite Element Method (FEM) is illustrated by means of the linear advection problem. A special attention is given to the identification of the stabilization parameter τ of this method which weights the stabilization terms. The need for stabilization is outlined and the basic idea of the Petrov-Galerkin (PG) concept are discussed.

Preliminaries

For the arguments to follow we would like to introduce some notations and definitions.

The first space we need is the space of square-integrable functions L 2 (Ω) defined by:

L 2 (Ω) = u | Ω | u | 2 d Ω < ∞
where Ω an open internal of R.

There is a norm related to this space, denoted by ∥ . ∥ L 2 (Ω) which is defined by:

∥ u ∥ L 2 (Ω) = Ω u 2 d Ω 1 2 .
We also will need one of the wide family of Sobolev spaces H 1 (Ω) which is defined by:

H 1 (Ω) = u ∈ L 2 (Ω) such that d u d x ∈ L 2 (Ω) . H 1 (Ω) is a Hilbert space with inner product (v, w) 1,Ω = Ω v(x)w(x)d x + Ω v (x)w (x)d x,
and the associated norm:

∥ v ∥ 1,Ω =∥ v ∥ 2 L 2 (Ω) + ∥ v ∥ 2 L 2 (Ω) .

Standard Galerkin FEM

Given a one-dimensional domain Ω =]a, b[, we consider the one-dimensional hyperbolic model problem of a scalar conservation law, with boundary and initial conditions, which can be written as:

         ∂u(x,t ) ∂t + ∂ ∂x f (u(x, t )) = 0 ∀(x, t ) ∈]a, b[×[0, T ], u(x, 0) = u 0 (x) ∀x ∈]a, b[, u(a, t ) = u a (t ) ∀t ∈ [0, T ], (5.1)
where u denotes an unknown scalar variable, while f (u) is called flux function and ∂u ∂t indicates the partial derivative of u with respect to time t , and t 0 = 0 indicates initial time. In addition, we represent by ∂u ∂x the partial derivative of u with respect to space.

The classical variational approach associated with (5.1) is obtained by multiplying this equation by a test function v supposed to be sufficiently regular and by integrating over the domain Ω. We denote by V 1 such space of functions, verifying the boundary conditions:

V 1 = w ∈ H 1 (Ω) such that w(a) = u a .
We introduce the space H 1 ([0, T ],V 1 ) of functions v: from [0, T ] to V 1 , such that ∂v ∂t ∈ L 2 ([0, T ]). The weak form of the state equation (5.1) is given by: Find u ∈ H 1 ([0, T ],V 1 ) such that:

Ω ∂u(x, t ) ∂t v(x)d x + Ω ∂ ∂x f (u(x, t ))v(x)d x = 0 ∀t ∈ [0, T ] ∀v ∈ H 1 (Ω). (5.2)
By integrating by parts and using Dirichlet boundary conditions, the variational formulation reads:

Find u ∈ H 1 ([0, T ],V 1 ) such that: Ω ∂u(x, t ) ∂t v(x)d x = Ω f (u(x, t ))v (x)d x + f (u a )v(a) -f (u b )v(b) ∀v ∈ H 1 (Ω). (5.3)
As consequence, the weak formulation of the problem (5.1) may be rewritten as:

   Find u ∈ H 1 ([0, T ],V 1
) such that :

a(u, v) = L(v) ∀v ∈ H 1 (Ω), (5.4)
where,

a(u, v) = Ω ∂u(x, t ) ∂t v(x)d x - Ω f (u(x, t ))v (x)d x,
and

L(v) = f (u a )v(a) -f (u b )v(b).

Lagrange P 1 elements

We first consider a mesh of the 1D computational domain Ω =]a, b[, where we want to compute the solution. A mesh is simply a set of points (x i ) 1≤i ≤N or intervals The cell size or space step is defined by

Ω i = [x i , x i +1 ], ∀1 ≤ i ≤ N -1 such that: a = x 1 < x 2 < ... < x N = b.
h 1 = b-a N -1 ,
where N -1 is the number of cells in the mesh. The coordinates of the grid points are then defined by

x i +1 = x 1 + i h 1 , ∀1 ≤ i ≤ N -1.
The FEM for Lagrange P 1 elements involves the space of globally continuous affine functions on each interval:

V 1 h = w h ∈ C 0 (Ω) w h | Ω i ∈ P 1 (Ω i ), ∀i = 1, 2, ..., N -1 such that w h (a) = u a , V 1 h is a subspace of V 1 of dimension N . Moreover, every function u h ∈ V 1
h is uniquely determined by:

u h (x, t ) = N j =1 u h (x j , t )ϕ j (x) ∀x ∈ Ω,
where ∀2 ≤ j ≤ N -1, ϕ j is the basis of the shape functions defined as:

ϕ j (x) =          x-x j -1 h 1 if x ∈ [x j -1 , x j ], x j +1 -x h 1 if x ∈ [x j , x j +1 ], 0 otherwise.
These functions are shown in Fig. 5.2. Their construction involves satisfies: By choosing the test functions v h equal to the basis functions of V 1 h , one obtains:

ϕ j (x i ) = δ i j ∀i , j = 1, ..., N .
Find u h ∈ V 1 h such that for all i = 1, ..., N , N j =1 ∂u j (t )
∂t Let us consider the problem:

Ω ϕ i (x)ϕ j (x)d x = Ω f (u h (x, t ))ϕ i (x)d x + f (u a )ϕ i (a) -f (u b )ϕ i (b) ∀t ∈ [0, T ]. ( 5
         ∂u(x,t ) ∂t + c ∂u(x,t ) ∂x = 0 ∀(x, t ) ∈] -1, 1[×[0, T ], u(x, 0) = u 0 (x) = sin(2πx) ∀x ∈] -1, 1[, u(-1, t ) = u a (t ) ∀t ∈ [0, T ], (5.6)
where, u a (t ) = sin(2π(-1ct )).

(5.7)

Note that no boundary condition is prescribed at x = 1 due to the hyperbolicity of the problem.

The exact solution to this problem is:

u ex (x, t ) = u 0 (x -c t ) = sin(2π(x -c t )).
We present in Fig. 5.3 the exact solution of the problem (5.6) for a transport speed c = 1. Thus, we get the discrete weak formulation given by: Find

u 1 (t ) = u h (x 1 , t ), ..., u N (t ) = u h (x N , t ) such that for all i = 1, ..., N N j =1 ∂u j (t )
∂t

Ω ϕ i (x)ϕ j (x)d x = c N j =1 u j (t ) Ω ϕ j (x)ϕ i (x)d x + cu a ϕ i (a) -cu(b)ϕ i (b) ∀t ∈ [0, T ], (5.8)
which can be rewritten in the form of linear system:

M ∂ t U = RU + cU a -cU b , (5.9)
where U (t ) is the vector whose components are u j (t ), the unknown values at the grid point at time t, M and R are the mass and stiffness matrices whose coefficients are:

M i j = Ω ϕ i (x)ϕ j (x)d x ∀i , j = 1, ..., N , R i j = c Ω ϕ i (x)ϕ j (x)d x ∀i , j = 1, ..., N , and
U =      u 1 . . . u N      ∈ R N , U a =          sin(2π(-1 -c t )) 0 . . . 0          ∈ R N , U b =          0 0 . . . u N          ∈ R N .
In the FEM, the matrices M and R are computed from the corresponding elementary matrices which are obtained by applying a local change of variables to map each cell Ω onto the reference element.

After integration, we get:

M = h 1 6                2 1 0 • • • • • • 0 1 4 1 0 • • • 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 • • • 0 1 4 1 0 0 • • • 0 1 2                and R = c 2                1 -1 0 • • • • • • 0 1 0 -1 0 • • • 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 • • • 0 1 0 -1 0 • • • 0 0 1 1                .
Therefore, the semi-discrete formulation in space is given by: ∀2 ≤ i ≤ N -1,

h 1 6 ∂ t u i -1 + 4∂ t u i + ∂ t u i +1 + c u i +1 -u i -1 2 ,
where, u i -1 , u i and u i +1 are the values of u at nodes i -1, i and i + 1 respectively.

We can conclude that classical Galerkin FEM has similarities with the use of central differenciation in the Finite Difference Method (FDM) for the advection problem for which the advection term writes c U i +1 -U i -1

2

(identical to -RU ).

The coefficients u i are time-dependent while the basis and test functions depend just on spatial coordinates. Further, the time derivative is not discretized in the time domain. One approach would be to use FEM for the time domain as well, but this can be rather computationally expensive. Alternatively, an independent discretization of the time domain is often applied using the method of lines. For example, it is possible to use the explicit Euler scheme:

Given the initial value problem:

   ∂U (x,t ) ∂t = L t ,U (x, t ) , U (x, t 0 ) = U 0 (x), (5.10)
we approximate the partial differential equation (5.10) by the finite-difference formulation:

U (x, t + ∆t ) -U (x, t ) ∆t = L t ,U (x, t ) .
This is exactly one step of the explicit Euler method. Introducing the notation:

   t ň+1 = t ň + ∆t , U ň = U (t ň ),
we have,

U ň+1 = U ň + ∆t L t ň ,U ň .
Euler explicit method is employed to integrate the equations in time with the timestep ∆t chosen small enough to ensure that time discretization errors can be neglected. Assuming that M is invertible, we have:

U ň+1 = U ň + ∆t M -1 RU ň + c∆t M -1 U ň a -U ň b .
(5.11)

This example demonstrates the limitations of the standard Galerkin method for solving advection problem, we observe in Fig. 5.4 that the approximate solution of the problem oscillates. Therefore, standard Galerkin FEM applied to these problems is far from "optimal" and gives unphysical oscillatory solutions, well known as Gibbs phenomena [START_REF] Brooks | Streamline Upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF]. The attenuation of these oscillations has been the subject of extensive research for several decades during which a huge number of so-called stabilized methods have been developed. The stabilizing effect can be often interpreted as the addition of some artificial diffusion to a standard (unstable) numerical scheme. On the one hand, this artificial diffusion should damp down the oscillations. On the other hand, it should not smear the numerical solution so that the design of a proper stabilization is a very difficult task.

In the late 1970s and early 1980s, a number of methods have been proposed to remove or, at least, to diminish these oscillations (a good summary of the very early literature is contained in [START_REF] Brooks | Streamline Upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF] [32] [START_REF] John | On spurious oscillations at layers diminishing (SOLD) methods for convection-diffusion equations: Part I I -Analysis for P 1 and Q 1 finite elements[END_REF]). Thus, the attention of FE researchers was turned to the development of Petrov-Galerkin methods (PGM), which are based on searching the test functions that provide exact nodal values for a selected class of solutions.

Using these test functions in the general case induces a stabilizing effect which removes the wiggles obtained with the classical Galerkin method [START_REF] Mizukami | A Petrov-Galerkin finite element method for convection-dominated flows: an accurate upwinding technique for satisfying the maximum principle[END_REF]. In the context of FEM, there are several approaches, among the most popular techniques, we can name the so called Streamline-Upwind Petrov-Galerkin (SUPG) [START_REF] Brooks | A Petrov-Galerkin finite element formulation for convection dominated flows[END_REF] [84], Galerkin Least-Squares (GLS) [START_REF] Billaud | Eléments finis stabilisés pour des écoulements diphasiques compressible-incompressible[END_REF] [30] [START_REF] Shakib | A new finite element formulation for computational fluid dynamics: x. the compressible Euler and Navier-Stokes equations[END_REF] and Pressure-Stabilizing Petrov-Galerkin (PSPG) [START_REF] Wervaecke | Simulation d'écoulements turbulents compressibles par une méthode d'éléments finis stabilisée[END_REF]. All these methods are based on a PG FEM. The main idea of all these methods is to add products of suitable perturbation terms and the residuals, thereby maintaining consistency. In these stabilized formulations, a judicious selection of the stabilization parameter, which is often denoted as τ, plays an important role in determining the accuracy of the formulation [1] [15]. This stabilization parameter requires special attention, as it strongly depends on the problem under consideration and the chosen numerical method.

We here give a brief outline of the theoretical basis for PG FEM. Given a differential operator L and a function G, we consider the problem:

L u(x) = G(x) ∀x ∈ Ω, (5.12)
where Ω ⊂ R (without loss of generality on R d ). The weak form of (5.12) is given by:

Ω (L u -G)v * d Ω = 0.
The problem is discretised based on the following formula:

u = N j =1 u j ϕ j , (5.13)
where {ϕ j } 1≤ j ≤N are mesh-based shape functions.

Choosing v * = ϕ i leads to the classical Galerkin method, whereas for the PGM we can have v * = ϕ i .

In a stabilized FEM, a disturbance is added to the test function of the Galerkin FEM which is given by:

v * = v + τL * (v),
where τ is a nonnegative stabilization parameter and L * is a differential operator which may or may not coincide with L . Let us examine some formulations:

Streamline-Upwind/Petrov-Galerkin (SUPG) method: In the context of FEM, a very popular stabilization technique is the SUPG method. This method, developed by Brooks and Hughes [START_REF] Brooks | Streamline Upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF], can be considered as the first successful stabilization technique to prevent oscillations in convection-dominated problems in the FEM. The basic idea of this method is to add diffusion (or viscosity) which acts only in the flow direction.

Extended to a Petrov-Galerkin formulation, the standard Galerkin test functions are modified by adding a streamline upwind perturbation, which again acts only in the flow direction. The modified test function is applied to all terms in the equation, resulting in a consistent weighted residual formulation.

Galerkin/Least-Squares (GLS) method: An alternative stabilization technique, known as Galerkin/Least Squares (GLS) formulation, was introduced in 1988 by Hughes, Franca and Hulbert [START_REF] Wervaecke | Simulation d'écoulements turbulents compressibles par une méthode d'éléments finis stabilisée[END_REF]. It can be interpreted as a generalization of the SUPG method. It is similar to the SUPG in some aspects, and for the hyperbolic equations and/or piecewise linear interpolation functions in the general case, the two methods become identical. The way GLS works is as follows: Least-squares forms of residuals are added to the Galerkin method. These terms enhance the stability of the Galerkin method without degrading accuracy.

The result is that practically convenient interpolations, which are unstable within the Galerkin framework, become convergent.

Pressure-Stabilizing

SUPG FEM for one-dimensional linear advection problem

Selection of the SUPG stabilization parameter

Like many other stabilized methods, the SUPG method contains a stabilization parameter, τ, for which a general "optimal" choice is not known. Since the SUPG method attracted a considerable attention, many research works had also been devoted to the choice of the parameter τ [1] [START_REF] Carbonel | Numerical study of Petrov-Galerkin formulations for the shallow water wave equations[END_REF]. Theoretical investigations of model problems only provide asymptotic behavior of this parameter (with respect to the local mesh size)

and certain bounds for which the SUPG method is stable and leads to (quasi-)optimal convergence of the discrete solution u h . The choice of τ may dramatically influence the accuracy of the discrete solution and therefore it has been a subject of an extensive research. However, the stabilization parameter τ depends on the problem under consideration and unfortunately, a general optimal definition of τ is still not known.

Note that for 1D advection-diffusion problems, an optimal value of the stabilization parameter can be defined.

In our computations, for the one-dimensional linear advection problem, we define τ, on any element Ω i , by the formula:

τ = α h i c , (5.14)
where, h i is the length of Ω i , in the present case is simply h i = h 1 and α is a parameter to be determined in [0, 1]. We then try to look for the optimal choice of the coefficient α ∈ [0, 1].

SUPG finite element approximation

The standard Galerkin FEM produces non-physical oscillations for the advection problem (5.6), that pollute the whole computational domain. Because of this undesirable feature of the Galerkin method, several approaches have been proposed to cure this problem within the framework of FEMs. In this thesis, we investigate the most favorite one: the SUPG stabilisation method, which adds an additional term to the Galerkin FEM to control the derivatives in the streamline direction. The SUPG weak form of (5.6) can then be written as follows:

Ω ∂u(x, t ) ∂t + c ∂u(x, t ) ∂x v(x)d x + N -1 k=1 x k+1 x k τc ∂v(x) ∂x ∂u(x, t ) ∂t + c ∂u(x, t ) ∂x d x = 0 ∀v ∈ V 1 . (5.15)
For the sake of generality, we adopt here a classical formulation, where stabilization terms are integrated only inside the elements to avoid possible regularity problems at the interfaces.

After integration by parts, one obtains:

Ω ∂u(x, t ) ∂t v(x)d x + N -1 k=1 τc x k+1 x k ∂u(x, t ) ∂t ∂v(x) ∂x d x -c Ω u(x) ∂v(x) ∂x d x + N -1 k=1 τc 2 x k+1 x k ∂u(x, t ) ∂x ∂v(x) ∂x d x -cu a v(a) + cu b v(b) = 0 ∀v ∈ V 1 .
Therefore, the weak formulation reads:

   Find u ∈ H 1 ([0, T ],V 1
) such that :

a SU PG (u, v) = L(v) ∀v ∈ V 1 (Ω), (5.16)
where

a SU PG (u, v) = Ω ∂u(x, t ) ∂t v(x)d x-c Ω u(x, t ) ∂v(x) ∂x d x+τc N -1 k=1 x k+1 x k ∂u(x, t ) ∂t ∂v(x) ∂x d x+τc 2 N -1 k=1 x k+1 x k ∂u(x, t ) ∂x ∂v(x) ∂x d x,
and

L(v) = cu a v(a) -cu b v(b).
Let now the space H 1 (Ω), in which the solution of (5.1) is sought, be approximated by a conforming FE subspace

V 1 h .
Thus, the discritized SUPG method reads as follows:

             Find u h ∈ H 1 ([0, T ],V 1 ) such that : Ω ∂u h (x, t ) ∂t v h (x)d x -c Ω u h (x, t ) ∂v h (x) ∂x d x + τc N -1 k=1 x k+1 x k ∂u h (x, t ) ∂t ∂v h (x) ∂x d x + τc 2 N -1 k=1 x k+1 x k ∂u h (x, t ) ∂x ∂v h (x) ∂x d x -cu a (t )v(a) + cu b (t )v(b) = 0.
It is important to point out that the SUPG solution with τ = 0 corresponds to the standard semi-discrete Galerkin approximation.

By choosing the test functions v h as the basis functions of V 1 h , the approximate problem can be written as follows:

                 Find u 1 , u 2 , ..., u N such that : N j =1 Ω ϕ j (x)ϕ i (x) ∂u j (t ) ∂t d x + τc N -1 k=1 x k+1 x k ϕ j (x)ϕ i (x)u j (t )d x -c Ω ϕ j (x)ϕ i (x)u j (t ) +τc 2 N -1 k=1 x k+1 x k ϕ i (x)ϕ j (x)u j (t )d x -cu a + cu b = 0.
(5.17)

A consistent SUPG-FE spatial discretization (5.17) leads to the following set of ordinary differential equations:

M s 1 ∂ t U = R s 1 U + cU a -cU b . (5.18)
The elements of mass matrix M s 1 and the stiffness matrix R s 1 are given by:

M s 1 i j = Ω ϕ j (x)ϕ i (x)d x + τc N -1 k=1 x k+1 x k ϕ i (x)ϕ j (x)d x ∀i , j = 1, ..., N , (5.19)
and

R s 1 i j = c Ω ϕ i (x)ϕ j (x)d x -τc 2 N -1 k=1 x k+1 x k ϕ j (x)ϕ i (x)d x ∀i , j = 1, ..., N . (5.20)
By evaluating the integrals, we get:

M s 1 = 1 6                2h 1 + 3τc h 1 -3τc 0 • • • • • • 0 h 1 + 3τc 4h 1 h 1 -3τc 0 • • • 0 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 • • • 0 h 1 + 3τc 4h 1 h 1 -3τc 0 • • • • • • 0 h 1 -3τc 2h 1 + 3τc                .
The interpretation of the stabilization term as an additional diffusion is clear here.

R s 1 =                -τc 2 h 1 + c 2 τc 2 h 1 -c 2 0 • • • • • • 0 τc 2 h 1 + c 2 -2τc 2 h 1 τc 2 h 1 -c 2 0 • • • 0 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 • • • 0 τc 2 h 1 + c 2 -2τc 2 h 1 τc 2 h 1 -c 2 0 • • • • • • 0 τc 2 h 1 + c 2 -τc 2 h 1 + c 2                .
It is noteworthy that M s 1 is not diagonal. It is common in the literature to approximate (5.18) in time by means of explicit time stepping. To avoid having to solve linear systems involving the mass matrix at each time step, it also common to simplify (5.18) by lumping the mass matrix.

Mass lumping

Mass lumping is a numerical technique employed in FEM that has been widely used in different applications (like heat equation, wave equation and time-dependent transport equation). This technique consists

of replacing the consistent mass matrix by a diagonal matrix whose entry in row i is the sum of all the entries of the consistent mass matrix in row i , ∀1 ≤ i ≤ N , usually referred to as the lumped mass matrix (for more details see [START_REF] Guermond | A correction technique for the dispersive effects of mass lumping for transport problems[END_REF] [90] [START_REF] Wervaecke | Simulation d'écoulements turbulents compressibles par une méthode d'éléments finis stabilisée[END_REF]).

Mass lumping can be shown in one space dimension to be equivalent to approximate the consistent mass matrix by using the following trapezoidal quadrature rule:

x j +1 x j f (x)d x ≈ (x j +1 -x j ) f (x i ) + f (x i +1 ) 2 .
(5.21)

This quadrature is exact for linear polynomials. Using this quadrature, the mass matrix coefficients can be approximated as follows: ∀2 ≤ i , j ≤ N -1,

x j +1 x j -1 ϕ i (x)ϕ j (x)d x = x j x j -1 ϕ i (x)ϕ j (x)d x + x j +1 x j ϕ i (x)ϕ j (x)d x ≈ h 1 2 ϕ i (x i )ϕ j (x i )δ i , j + ϕ i (x i -1 )ϕ j (x i -1 )δ i -1, j + h 1 2 ϕ i (x i )ϕ j (x i )δ i , j + ϕ i (x i +1 )ϕ j (x i +1 )δ i +1, j = h 1 2 x i -x i -1 h 1 =1 x i -x j -1 h 1 δ i , j + x i -1 -x i -1 h 1 =0 x i -1 -x j -1 h 1 δ i -1, j + h 1 2 x i +1 -x i h 1 =1 x j +1 -x i h 1 δ i , j + x i +1 -x i +1 h 1 =0 x j +1 -x i +1 h 1 δ i +1, j = h 1 δ i j .
Where δ i j is the Kronecker symbol.

τc x j +1 x j -1 ϕ i (x)ϕ j (x)d x = τc x j x j -1 ϕ i (x)ϕ j (x)d x + τc x j +1 x j ϕ i (x)ϕ j (x)d x = τc h 1 2 1 h 1 x j -1 -x j -1 h 1 + 1 h 1 x j -x j -1 h 1 + τc h 1 2 -1 h 1 x j +1 -x j h 1 - 1 h 1 x j +1 -x j +1 h 1 = τc h 1 2 ( 1 h 1 - 1 h 1 ) = 0.
Hence, the coefficients of the lumped mass matrix are:

M s L i j = h 1 δ i j ∀2 ≤ i , j ≤ N -1
This technique of mass lumping presents a computational advantage as we need mass matrix inversion.

Upon replacing the consistent mass matrix M s 1 by the lumped mass matrix M s L , we obtain a new matrix form of advection equation which is fully explicit and writes as follows:

M s L ∂ t U = R s 1 U + cU a -cU b .
(5.22)

∂ t U = (M s L ) -1 R s 1 U + c(M s L ) -1 U a -U b ,
where, the so-called lumped mass matrix M s L thus computed is diagonal, is given by:

M s L = h 1 2             1 0 0 • • • 0 0 2 0 • • • 0 . . . . . . . . . . . . 0 • • • 0 2 0 0 • • • 0 0 1             .

Runge-Kutta time discretization

The space semidiscrete problem (5.22) represents a system of ordinary differential equations (ODEs).

This system can be solved by any of the available ODE solvers. In this thesis we focus on the Runge Kutta schemes. These methods have been developed [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF] [80] for solving:

∂ t U = L t ,U (x, t ) ,
where L t ,U (x, t ) is a spatial discretization operator.

We divide the time interval (0, T ) into ň time steps [t k , t k+1 ] ∀k = 0, ..., ň -1, where t 0 = 0 and t ň = T and denote the step length of interval [t k , t k+1 ] by ∆t = t k+1t k .

First-order Runge-Kutta formula (RK 1):

U ň+1 = U ň + ∆t L t ,U ň .

Second-order Runge-Kutta formula (RK 2):

We define an intermediate estimate U ň+ 1 2 :

U ň+ 1 2 = U ň + ∆t 2 L t ,U ň , U ň+1 = U ň + ∆t L t + ∆t 2 ,U ň+ 1 2 .
Third-order Runge-Kutta formula (RK 3):

U ň+1 = U ň + ∆t 6 K 1 + 4K 2 + K 3 ,
where,

K 1 = L t ň ,U ň , K 2 = L t ň + ∆t 2 ,U ň + ∆t 2 K 1 , K 3 = L t ň + ∆t ,U ň -∆t K 1 + 2∆t K 2 .
Fourth-order Runge-Kutta formula (RK 4):

U ň+1 = U ň + ∆t 6 K 1 + 2K 2 + 2K 3 + K 4 ,
where,

K 1 = L t ň ,U ň , K 2 = L t ň + ∆t 2 ,U ň + ∆t 2 K 1 , K 3 = L t ň + ∆t 2 ,U ň + ∆t 2 K 2 , K 4 = L t ň + ∆t ,U ň + ∆t K 3 .
In all these formula, U ň+1 is an approximation to U (t ň+1 ) and K 1 , K 2 , K 3 and K 4 are the intermediate evaluations.

Courant-Friedrichs-Lewy (CFL) condition

During the time integration cycles, the length of the time step needs to be chosen, according to a stability criterion ruled by the C C F L number. The concept of the C C F L number was originally published in 1928, the aim was to prove the existence of solutions of some PDEs. Consequently, while proving the existence, Courant, Friedrichs and Lewy found the necessary condition to stabilize the numerical methods. For more detailed discussion of the C C F L number see [START_REF] Hughes | Finite element methods for first-order hyperbolic systems with particular emphasis on the compressible Euler equations[END_REF]. The Courant Friedrichs-Lewy (CFL) condition is given by:

| c | ∆t h 1 ≤ C C F L ,
where ∆t is the time step, c is the speed and h 1 is the length for the computational element.

Note that this condition may not be sufficient in the presence of the SUPG stabilization term, for which a diffusion stability criterion has to be taken into account.

Numerical results

In this section, the performance of the SUPG method will be illustrated for the linear one-dimensional advection problem. We use equidistant grids with mesh size h 1 = 2 N -1 . It is known that standart Galerkin discretization results in a strongly non-stable scheme which then leads to a numerical solution exhibiting non-physical phenomena such as spurious numerical oscillations. We get such results on the left in Fig. 5.5.

On the right, the SUPG numerical solution is stable thanks to the additional diffusion term. The time stepping is done with the standard RK 2 method; this ensures that the error induced by the time approximation is small compared to the spatial error. The solution is computed at T = 0.4. Special attention is given to the role of the stabilization parameter τ (more precisely of the coefficient α) given by the formula:

τ = α h 1 c .

Influence of the SUPG parameter

The SUPG FEM contains a stabilization parameter τ for which a general "optimal" choice is not known.

It is pointed out that the optimal choice of this parameter is still an open question. The aim of this section is to describe how the SUPG stabilization parameter impacts the solution for the one-dimensional advection problem. The influence of this parameter on the numerical results is illustrated in Fig. 5.6. We can see clearly in Fig. 5.6 important oscillations for 0 ≤ α ≤ 0.01 and α ≥ 0.6. For α too small, the stabilization effect is not strong enough to counterbalance the instability of Galerkin scheme. When α is too large, oscillations are due to the violation of the stability criteria in RK time integration for diffusive terms. For 0.01 < α ≤ 0.5, we do not observe oscillations. However, a too large additional diffusion (α = 0.5) impacts the accuracy of the solution, even if the scheme is stable. In what follows, we propose to determine the optimal choice by studying the L 2 -error.

Error Estimates for the SUPG FE method

We want to evaluate the numerical behavior of the SUPG linear FEM for the problem (5.1), where an analytic solution is explicitly known. The interest here is in the evaluation of the errors committed for different mesh steps h 1 as well as the order of convergence of the method. A mesh convergence study allows to test the accuracy of the numerical method by numerically evaluating an order of accuracy. The error between the approximate solution u h and the analytic solution u ex for a final time T is determined in the standard

L 2 -norm: ∥ e(T ) ∥ L 2 (Ω) = ∥ u ex (x, T ) -u h (x, T ) ∥ L 2 (Ω) = Ω u ex (x, T ) -u h (x, T ) 2 d x 1 2 = N -1 k=1 x k+1 x k sin(2π(x -cT )) -u h (x, T ) 2 d x 1 2
.

Before we proceed further, let us denote by r the convergence rate of the SUPG FEM. We assume that a norm ∥ e(T ) ∥ L 2 (Ω) of the computational error behaves according to the formula:

∥ e(T ) ∥ L 2 (Ω) = C r h r 1 ,
where C r > 0 is a constant, h 1 the mesh size and r ∈ R.

Let ∥ e 1 (T ) ∥ L 2 (Ω) and ∥ e 2 (T ) ∥ L 2 (Ω) be computational errors of the numerical solutions obtained on two different meshes h 1 and h 2 , respectively. Then we obtain:

r = ln ∥ e 2 (T ) ∥ L 2 (Ω) -ln ∥ e 1 (T ) ∥ L 2 (Ω)
ln h 2ln h 1 . We observe that the convergence rates obtained are lower than the optimal rate expected (of value 2). As

shown in [START_REF] Liou | Finite element solution techniques for large-scale problems in computational fluid dynamics[END_REF] [91], this reduction of the convergence rate is due to the use of the mass lumping procedure, in addition to the use of the stabilization term (for some τ values).

SUPG FE method for high-order elements

In some applications, the affine approximation on each element of the mesh Ω i = [x i , x i +1 ] can be considered as not good enough in the sense that it provides an approximate function that is too far from the exact function u. To overcome this problem, we can approximate u on each mesh by polynomials of higher degree. We now extend the above considerations to higher-order finite-elements. We particularly focus our attention in this section on the FE approximation P 2 , which consists in approaching the solution u by a continuous function on Ω and a polynomial of degree 2 on each element Ω i . Now, let us revisit the one-dimensional advection problem that is given in (5.6) and we examine the quadratic Lagrange bases functions. These are constructed by adding an extra node x i + 1 2 at the midpoints of each The FEM for Lagrange P 2 elements involves the discrete space:

Ω i = [x i , x i +1 ], x i + 1 2 = x i + h 1 2 , ∀1 ≤ i ≤ N -1.
V 2 h = w h ∈ C 0 ([-1, 1]) w h | Ω i ∈ P 2 ∀1 ≤ i ≤ N -1 such that w h (a) = u 2a .
These spaces are composed of piecewise polynomials functions (polynomials of degree less than or equal to 2). The P 2 FEM consists in applying the internal variational approximation approach to these spaces.

V 2 h is a subspace of H 1 (Ω) of dimension 2N -1.
As with the piecewice-linear basis, one basis function is associated with each node. Those associated with vertices are:

φ j (x) =                        2 (x-x j -1 )(x-x j -1 2 ) h 2 1 if x j -1 ≤ x ≤ x j , 2 (x j +1 -x)(x j + 1 2 -x) h 2 1 if x j ≤ x ≤ x j +1 , 0 otherwise, (5.23)
and those associated with element midpoint are:

φ j + 1 2 (x) =          4 (x j +1 -x)(x-x j ) h 2 1 if x j ≤ x ≤ x j +1 , 0 otherwise.
(5.24)

These functions are shown in Fig. 5.9. Their construction (to be described) involves satisfying:

   φ j (x i ) = δ i j φ j (x i + 1 2 ) = 0 and    φ j + 1 2 (x i ) = 0 φ j + 1 2 (x i +1 ) = δ i j Figure 5
.9: Global shape functions for the space V 2 h .

The SUPG variational formulation of the problem (5.6) consists in finding

u h ∈ H 1 ([0, T ],V 2 h ). Every function u h ∈ V 2
h defined by its values at the mesh vertices (x i ) 1≤i ≤N and at the midpoints (x i + 1 2 ) 1≤i ≤N -1 :

u h (x, t ) = N j =1 u j (t )φ j (x) + N -1 j =1 u j + 1 2 (t )φ j + 1 2 (x) ∀x ∈ Ω.
Furthermore, since it is enough to use basis functions of V 2 h as test functions v, the discrete weak formulation can be rewritten in form of linear differential system, wich can be described as follows:

M s 2 ∂ t U = R s 2 U + c U 2a -U 2b .
(5.25)

Since the shape functions φ i have a compact support, the matrices M s 2 and R s 2 are mostly composed of zeros. In contrast to the Lagrange P 1 FEM, the matrices M s 2 and R s 2 are no longer tridiagonal matrices. Their elements are given by:

M s 2 i j = Ω φ j (x)φ i (x)d x + Ω φ j + 1 2 (x)φ i + 1 2 (x)d x + τc N -1 k=1 x k+1 x k φ i (x)φ j (x)d x + N -1 k=1 x k+1 x k φ i + 1 2 (x)φ j + 1 2 (x)d x , and R s 2 i j = c Ω φ i (x)φ j (x)d x + Ω φ i + 1 2 (x)φ j + 1 2 (x)d x -τc 2 N -1 k=1 x k+1 x k φ j (x)φ i (x)d x + N -1 k=1 x k+1 x k φ j + 1 2 (x)φ i + 1 2 (x)d x ,
where,

U 2a =          sin(2π(-1 -ct )) 0 . . . 0          ∈ R 2N -1 and U 2b =          0 0 . . . u N          ∈ R 2N -1
In this case, the elementary contributions of the element Ω i = [x i , x i +1 ] to the stiffness matrix and to the mass matrix are given by the 3 × 3 elementary matrices R s 2e | Ω i and M s 2e | Ω i :

R s 2e | Ω i =            -7τc 2 3h 1 -2c 3 + 8τc 2 3h 1 c 6 -τc 2 3h 1 2c 3 + 8τc 2 3h 1 -16τc 2 3h 1 -2c 3 + 8τc 2 3h 1 -c 6 -τc 2 3h 1 2c 3 + 8τc 2 3h 1 -7τc 2 3h 1            , and
M s 2e | Ω i =            2h 1 15
h 1 15 -2τc 3 -h 1 30 + τc 6 h 1 15 + 2τc 3 8h 1 15 h 1 15 -2τc 3 -h 1 30 -τc 6 h 1 15 + 2τc 3 2h 1 15            .
In the present case of P 2 elements, Simpson's rule can be derived. In particular, let the basis function φ j be tabulated at points x j , x j +1 ∀1 ≤ j ≤ N -1 (equally spaced by distance h 1 ), and midpoint

x j + 1 2 = x j +x j +1 2 .
Then Simpson's rule states that:

x j +1 x j f (x)d x ≈ x j +1 -x j 6 f (x j ) + 4 f (x j + 1 2 ) + f (x j +1 ) . (5.26)
Since it uses quadratic polynomials to approximate functions, Simpson's rule actually gives exact results when approximating integrals of polynomials up to cubic degree.

Hence, the so-called lumped mass matrix M s 2 thus computed is diagonal and is given by:

M s 2L | Ω i = 1 6            h 1 -3τc 0 0 0 4h 1 0 0 0 h 1 + 3τc            .
Upon replacing the consistent mass matrix M s 2 by the lumped mass matrix M s 2L , we rewrite the matrix form 5.25 as follows:

∂ t U = M s 2L -1 R s 2 U + c M s 2L -1 U 2a -U 2b .
(5.27)

Matrix assembly

The construction of the global mass and stiffness matrices, M s 2L and R s 2 is often referred to as "assembling" due to its method of construction. The assembly is obtained algorithmically using a loop over all mesh elements Ω i , ∀1 ≤ i ≤ N -1 (which is composed of three nodes x j , x j + 1 2 and x j +1 ) and adding their contributions to the corresponding coefficients of the global system. It is more convenient to loop over each element instead of looping over all basis functions at the outer loop.

The algorithm of constructing the global mass and stiffness matrices M s 2L and R s 2 can be described as follows:

Pseudo-code of matrix assembly:

1. For k = 1 : N -1 2.
For i = 1 : 3

3. i g = 2(k -1) + i 4.
For j = 1 : 3 5.

j g = 2(k -1) + j 6. M s 2L (i g , j g ) = M s 2L (i g , j g ) + m s 2e (i , j ) 7. R s 2 (i g , j g ) = R s 2 (i g , j g ) + r s 2e (i , j ) 8.
End (loop j )

9.

End (loop i )

End (loop k)

The space semidiscrete problem (5.25) represents a system of ODEs, which has to be solved with the RK schemes. The choice of the timestep that ensures a stable scheme, for the quadratic Lagrange P 2 SUPG FE method is based on the condition:

∆t ≤ h 1 C C F L 2 | c | .

Numerical results

To give an illustration of the SUPG quadratic Lagrange FEM in conjunction with RK 2 method for the temporal discretization, some numerical results are presented in this section. The exact and numerical In Fig. 5.11 we show the L 2 -error norm depending on the numbers of d.o.f for the optimal choice of α.

Conclusion

In the present chapter we have focuced our attention to the SUPG FEM. As mentioned before, no standard scheme has been formulated to select the stabilisation parameter τ, which depends on the problem and the grid size. In our work, we have assumed that the stabilisation parameter τ takes the form:

τ = α h c ,
and thus we have focused our attention on the selection of α ∈ [0, 1].

In order to have a precise idea for the linear and quadratic Lagrange SUPG FEM in conjunction with RK 2

for time discretization, we represented the L 2 -error norm for an optimal choice for each method for a successive refinement of the mesh. Even if the convergence rate does not seem to be strongly impacted by the choice of α, the influence of the stabilization on the error level was reported. Une attention particulière a été accordée à l'identification du paramètre de stabilisation τ de cette méthode qui pondère les termes de stabilisation, et son influence sur l'erreur de la solution.

C H A P T E R

ISOGEOMETRIC ANALYSIS: B-SPLINE AS A FEM BASIS

I Sogeometric Analysis (IGA) is a generalization of classical FEA. It was first introduced by T.J.R. Hughes, J.A. Cottrell and Y. Bazilevs [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF] with the main aim of closing the gap between the geometrical description and the analysis of engineering problems. The isogeometric paradigm consists of using basis functions commonly found in CAD geometries such as B-spline, to represent both the geometry and the physical fields in the solution of problems governed by partial differential equations (PDEs) [43] [67]. However, Hughes and co-authors considered only elliptic or parabolic problems so far. The objective of this chapter is to consider IGA in the hyperbolic context. This chapter reviews the various computational procedures for IGA, by revisiting the one-dimensional advection problem that is given in (5.6). While in this chapter we use B-splines (due to the simplicity of the domain), it is not hard to generalize it to other splines such as NURBS. Detailed comparisons between IGA and FEA approaches are carried out.

IGA: a B-spline based approach

IGA employs the same mathematical foundations as FEA to obtain the numerical solutions of differential equations, but the idea behind IGA method is to use the same functions that define the physical domain Ω ⊂ R 2 , to approach the solution. Although the main difference between FEA and IGA lies compactly in the set of basis functions used, this change influences all steps of traditional FEA: preprocessing, solving and postprocessing.

The basis functions of IGA B-splines (more generally NURBS) are defined on a parameter space (or parameter domain) that we denote Ω. It is defined by control points, and not nodes like discretized domain in FEA. The input of node coordinates is replaced by coordinates of control points. Before we proceed further, let us notice that the B-spline parameter space Ω corresponds to "patches" instead of elements. A patch can be seen as a "macro-element", whereas knot intervals can be seen as traditional elements. Thus, the IGA allows to control the geometry and variables at the control points, in contrast to FEA which writes these quantities at the nodes. The transformation of the parametric domain Ω to the physical domain Ω is introduced:

T : Ω -→ Ω (ξ, η) -→ x(ξ, η), y(ξ, η) .
Note that this transformation is not linear and is simply defined by the parametric representation of the computational domain. Unless otherwise specified, we assume that the mapping T to be invertible and its inverse:

T -1 : Ω -→ Ω,
transforms points in the physical domain to their corresponding parameter values. Any point of coordinate (x, y) in the physical domain Ω is mapped to a point (ξ, η) in the parametric domain Ω. The mapping from the parametric domain to the physical domain is defined by associating a control net P i j to each basis function in such a way that:

x(ξ, η), y(ξ, η) = n i =1 m j =1
N i ,p (ξ)N j ,q (η)P i j .

Thus, the basis functions in the physical domain Ω are defined by:

N i j ,pq (x, y) = N i j ,pq T(ξ, η) = N i j ,pq (ξ, η) = N i ,p (ξ)N j ,q (η).
It is important to point out that in some applications (which uses the simple geometries such as square, circle/ellipse), the computational domain can be modeled straightforwardly with only a single patch. But for more complex geometries, it is not possible to describe the physical computational domain with just one geometrical mapping T, we need to use multi-patch geometries (we refer the reader to [26] [40]). In this thesis we focus on single-patch geometries for the sake of simplicity but the extension is quite simple.

In an isoparametric formulation, the variable field is approximated by the same shape functions:

u(ξ, η) = n i =1 m j =1 N i ,p (ξ)N j ,q (η)u i j .
Hence, the variable u is obtained as functions of ξ and η. Here u i j denotes the value of the variable field u corresponding to the control net P i j . It is therefore referred to as a control variable or more generally a degree of freedom (d.o.f.). Note that B-spline representations are not interpolatory, so u i j is not the solution value at the location P i j .

Computational procedures for IGA

The fact that the functions are defined in the parametric space while the space derivatives are with respect to the coordinates of the physical space makes it necessary to compute the Jacobian matrix of the geometric mapping defined as: We also highlight that for the one-dimensional case, the Jacobian matrix reduces to:

J ξ = ∂x ∂ξ = n i =1 ∂N i ,p (ξ) ∂ξ P i ,
and we denote, J -1 ξ = ∂ξ ∂x . We will try to give an idea of the calculation of the integrals. Integrals over the entire geometry (physical domain) are split into elementary integrals over a domain denoted by Ω e . Let f be a function of the two variables x and y. Then,

Ω f(x, y)d Ω = N el e=1 Ω e f(x, y)d Ω e = N el e=1 Ω e f(T(ξ, η))d Ω e = N el e=1 Ω e f(ξ, η) | J | d Ω e = N el e=1 Ω f ξ, η | J || J | d Ω (via the mapping T(F i g . 6.1)).
We also highlight that ξ and η on the parent domain Ω are given as functions of ξ and η on parametric

domain Ω i j = [ξ i , ξ i +1 ] × [η j , η j +1 ]: ξ = 1 2 (ξ i +1 -ξ i ) ξ + (ξ i +1 + ξ i ) , η = 1 2 (η j +1 -η j ) η + (η j +1 + η j ) .
Therefore, the Jacobian of this transformation reads:

| J |= 1 4 ξ i +1 -ξ i η j +1 -η j .
For the one-dimensional case | J | reduces to:

| J |= 1 2 ξ i +1 -ξ i .

Isogeometric FE formulation

After a general description of the physical domain, geometrical mappings, and an introduction of notations for the IG method, we will discusses in this section specific details of the implementation of IGFE method.

The IGM based on B-splines is not much different than the classical FEM. Subtle differences are introduced due to the non-interpolatory character of B-splines and the definition of an element.

We focus now on the one-dimensional advection problem that is given in the previous chapter by:

         ∂u(x,t ) ∂t + c ∂u(x,t ) ∂x = 0 ∀(x, t ) ∈ [a, b] × [0, T ], u(x, 0) = u 0 (x) ∀x ∈ [a, b], u(a, t ) = u a (t ) ∀t ∈ [0, T ], (6.1)
with,

   u 0 (x) = sin(2πx) ∀x ∈ [a, b], u a (t ) = sin 2π(-1 -ct ) ∀t ∈ [0, T ].
In order to solve Eq. ( 6.1), we consider a spatial discretization of its domain by means of B-spline functions IGA in the framework of the Galerkin method and we need to write the weak form. Before we proceed further, let us define:

V N := span N i , for i = 1, .., n , V N := span N i , for i = 1, .., n .
According to the IGA concept, these spaces will be used to build the test function spaces for the approximation of (6.1). However, we define a subspace V B which only span a finite number of B-spline basis functions.

That is:

V B = V ∩ V N ,
where,

V := w B ∈ H 1 (Ω) such that w B (a) = u a .
For IGA in the framework of the Galerkin method, we start by multiplying the strong form (6.1) by a B-spline test function v B and integrate over the domain Ω, then using Green's first identity and applying the boundary conditions we finally get:

Ω ∂ t u(x, t )v B (x)d x -c Ω u(x, t )∂ x v B (x)d x = cu a v B (a) -cu b v B (b).
As in the previous chapter, we introduce a SUPG stabilization term yielding:

               Find u ∈ H 1 ([0, T ],V B ) such that : Ω ∂u(x, t ) ∂t v B (x)d x -c Ω u(x, t ) ∂v B (x) ∂x d x + τc n-p k=1 x k+1 x k ∂u(x, t ) ∂t ∂v B (x) ∂x d x +τc 2 n-p k=1 x k+1 x k ∂u(x, t ) ∂x ∂v B (x) ∂x d x = cu a v B (a) -cu b v B (b) ∀v B ∈ H 1 ([a, b]).
The solution u is obtained by solving the finite-dimensional problem:

   Find u ∈ H 1 ([0, T ],V B )
such that :

a sp (u, v B ) = L sp (v B ) ∀v B ∈ H 1 ([a, b]). (6.2)
Where, a sp :

V B × H 1 ([a, b]) -→ R is given by: a sp (u, v B ) = Ω ∂u(x, t ) ∂t v B (x)d x -c Ω u(x, t ) ∂v B (x) ∂x d x + τc n-p k=1 x k+1 x k ∂u(x, t ) ∂t ∂v B (x) ∂x d x + τc 2 n-p k=1 x k+1 x k ∂u(x, t ) ∂x ∂v B (x) ∂x d x,
and

L sp : H 1 ([a, b]) -→ R
, that contains the right hand side term of (6.2) and

L sp (v B ) = cu a v B (a) -cu b v B (b).
Galerkin's method consists in constructing finite-dimensional approximations of V B , denoted V B h . Strictly speaking, these will be subsets such that: V B h ⊂ V B . Therefore, the SUPG formulation of this problem can be written:

             Find u h ∈ H 1 ([0, T ],V B h ) such that : Ω ∂u h (x, t ) ∂t v B h (x)d x -c Ω u h (x, t ) ∂v B h (x) ∂x d x + τc n-p k=1 x k+1 x k ∂u h (x, t ) ∂t ∂v B h (x) ∂x d x +τc 2 n-p k=1 x k+1 x k ∂u h (x,t ) ∂x ∂v B h (x) ∂x d x -cu a (t )v B h (a) + cu b (t )v B h (b) = 0.
According to the isogeometric paradigm, the solution u h from the IGA space H 1 ([0, T ],V B h ) can be represented using B-spline basis functions for the one-dimensional case, in the form:

u h (x, t ) = n j =1 N j ,p (x)u j (t ),
where n is the number of basis functions which is equal to the number of control points and the degrees of freedom (d.o.f.) of u h associated with the control points, respectively. This representation is similar to that used in classical FEM. By choosing the test function v B h equal to the B-spline basis function N i ,p , the Galerkin IGA scheme reads as follows:

                 Find u 1 , u 2 , ..., u n such that : n j =1 Ω N j ,p (x)N i ,p (x) ∂u j (t ) ∂t d x + τc n-p k=1 x k+1 x k N j ,p (x)N i ,p (x)u j (t )d x -c Ω N j ,p (x)N i ,p (x)u j (t ) +τc 2 n-p k=1 x k+1 x k N i ,p (x)N j ,p (x)u j (t )d x = cu a N i ,p (a) -cu b N i ,p (b) ∀i = 1, ..., n.
Finally, we obtain the linear system which is:

M B ∂ t U = R B U + cU B a -cU B b , (6.3)
similar to that resulting from the FEM in conjunction with SUPG.

The elementary mass and stiffness matrices are defined by:

M B i j = Ω N j ,p (x)N i ,p (x)d x + τc n-p k=1 x k+1 x k N j ,p (x)N i ,p (x)d x ∀i , j = 1, ..., n,
and

R B i j = c Ω N j ,p (x)N i ,p (x)d x -τc 2 n-p k=1 x k+1 x k N j ,p (x)N i ,p (x)d x ∀i , j = 1, ..., n,
and the right-hand side given by:

U B a =          sin(2π(-1 -ct )) 0 . . . 0          ∈ R n , and
U B b =          0 0 . . . u n          ∈ R n .
We note that on the element Ω j , only p +1 B-splines functions N j ,p , N j +1,p , ..., N j +p,p non zero and only p +1 derivatives N j ,p , N j +1,p , ..., N j +p,p are non zero. The correspond matrices are therefore sparse.

By performing the integration in the parametric space and using the Jacobian of the B-spline mapping, we get:

M B i j = Ω N j ,p (ξ)N i ,p (ξ)J ξ d ξ + τc n-p k=1 Ω k N j ,p (ξ)N i ,p (ξ)d ξ,
and

R B i j = c Ω N j ,p (ξ)N i ,p (ξ)d ξ -τc 2 n-p k=1 Ω k N j ,p (ξ)N i ,p (ξ)J -1 ξ d ξ.
By interpreting the integration in the parent element, we get:

M B i j = Ω N j ,p ( ξ)N i ,p ( ξ) Ĵ ξd ξ + τc n-p k=1 Ω k N j ,p ( ξ)N i ,p ( ξ)d ξ,
and

R B i j = c Ω N j ,p ( ξ)N i ,p ( ξ)d ξ -τc 2 n-p k=1 Ω k N j ,p ( ξ)N i ,p ( ξ) Ĵ -1 ξ d ξ.
Contrary to Lagrange elements, the element integrals appearing are evaluated using Gauss points via Gaussian quadrature which can be described as follows:

Ω N j ,p ( ξ)N i ,p ( ξ)d ξ n G k=1 N j ,p (X G k )N i ,p (X G k )ω G k ,
where n G are the number of integration points X G k , ω G k is the weight corresponding to the k-th integration point (see Appendix A). This approach is more flexible than the use of analytical integration.

The global mass and stiffness matrices are obtained using the assembly technique.

As mentionned before, we focus in this thesis on the RK method for the time integration. The ODE (6.3) is integrated in time by means of a second and fourth-order RK schemes. In one spatial dimension, a space discretization using polynomials basis functions of degree p, associated to a (p +1)-stage RK time integrator of order p + 1, the stability limit for the CFL number is defined by:

∆t | c | h ≤ C C F L .
where h is the size of the image of the knot interval considered.

Numerical results

The standard Galerkin formulation for the advection problem produces unstable discretizations. This is well known in FEA for the Lagrange P 1 and P 2 elements, and unfortunately, is also observed for the linear and quadratic B-splines-based approach as we can see in this section. We investigate the ability of the isogeometric approach based on B-splines of arbitrary degree, in conjunction with SUPG, to solve advection problem. In all calculations the mesh is uniform with element size length h = 2 n-p . The solutions are calculated from p = 1 to p = 4. Notice that in the case of linear and quadratic B-spline, the standard SUPG formulation is used with the choice of the stabilization parameter τ = αh c while for the cubic and quartic B-spline it is important to point out that the method has been found stable, and therefore τ = 0. Due to the equivalence of hat shape functions and linear B-spline functions, the linear FEM and IGA give the same results (as shown in Fig. 7.3). However, this is not the case for representations of higher degree.

If one compares quadratic B-spline basis functions with basis functions based on Lagrange interpolating polynomials, one can notice some critical differences.

CONCLUSION

L 'analyse isogéométrique (AIG) trouve ses origines depuis 2005 dans les travaux de Hughes, Cottrell, et Bazilevs [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF] avant d'être détaillée et formalisée dans le livre de Cottrell et al. [START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF] dans un contexte de développement de nouveaux outils pour le calcul numérique en ingénierie basée sur la simulation, celle-ci introduit un nouveau paradigme pour tenter de combler le fossé entre la MEF et la CAO.

L'idée principale de l'AIG est de modéliser exactement la géométrie avec des fonctions qui vont servir à approximer la solution, ces fonctions permettant une description paramétrique d'un domaine. Elles sont définies dans un espace appelé espace paramétrique; la géométrie est définie dans un espace nommé espace physique. La paramétrisation de la géométrie est obtenue par un morphisme non linéaire entre l'espace paramétrique et l'espace physique défini à l'aide des fonctions B-splines. S'il est nécessaire d'effectuer des intégrations lors du processus de résolution, elles sont effectuées élément par élément sur un élément de référence obtenu par une transformation linéaire à partir de l'espace paramétrique.

On a introduit brièvement dans ce chapitre le contexte de l'AIG. Pour cela, on a commencé par rappeler les fondements et les premiers objectifs de la méthode. On se restreint aux B-splines mais on ne doute pas de l'applicabilité des méthodes développées ici à d'autres représentation de l'AIG (comme NURBS), ces fonctions possèdent en fait une continuité de classe supérieure, ce qui permet d'obtenir plus de précision pour un même nombre de degrés de libertés comparé aux EF classiques.

Les travaux présentés dans ce chapitre ont pour objectif de mettre au point une méthode d'EF stabilisée, la méthode SUPG pour un problème d'advection dans le cadre de l'AIG. L'idée principale est donc d'utiliser les fonctions de base B-spline représentant la géométrie pour générer l'espace de recherche de la solution au problème hyperbolique souhaité.

.

Part III

ISOGEOMETRIC DISCONTINUOUS GALERKIN METHOD (IGDGM)

In fact, the FV method, which is well suited to hyperbolic conservation laws, can only use low-degree polynomials to locally represent the solution. In contrast, FEM are unstable for hyperbolic problems and, as seen in previous chapters, stabilization relies on tedious choice of a parameter. Therefore, the idea of this method is to decompose the original problem into a set of subproblems, solved by using a FEM approach, that are connected using an appropriate transmission condition (known as the numerical flux).

Although DG methods have gained increasing attention in large-scale modeling applications, a shortcoming of the conventional DG methodology is the inability to fully recover complex underlying geometries in the meshing domain. To overcome this problem, we combine IGA method with the DG method. As mentioned before, IGA is a computational technique that improves and generalizes the classical FE method.

The main benefit of this method is the exact representation of the geometry in the language of computer aided design (CAD) tools. This simplifies the meshing as the computational mesh is directly created by the engineer using the CAD tools.

The proposed IGDG method is the DG method formulated on elements that exactly preserve the geometry generated by CAD tools while the PDE solution exhibits discontinuities at element interfaces. An important property of B-spline in the context of DG is the ability to perform Bézier extraction. Bézier extraction provides the capability of recovering a local Bernstein-Bézier representation of the geometry from the global B-spline CAD. In this chapter, we will discuss specific details of the implementation of IGDG method for the one-dimensional linear advection problem in contrast to the DGFE method.

DGFE framework for one-dimensional scalar conservation law

The DGFE method was first designed as an effective numerical method for solving hyperbolic conservation laws [START_REF] Proft | Discontinuous Gelerkin methods for convection-diffusion equations for varying and vanishing diffusivity[END_REF] [START_REF] Xu | Discontinuous Galerkin method for fractional convection-diffusion equations[END_REF]. In this section, we will present the details of the DGFE method, the stability analysis, and the error estimates for the one-dimensional scalar conservation law given by:

         ∂u(x,t ) ∂t + ∂ f (u(x,t )) ∂x = 0 ∀(x, t ) ∈ Ω × [0, T ], u(x, 0) = u 0 (x) ∀x ∈ Ω, u(a, t ) = u a (t ) ∀t ∈ [0, T ]. (7.1)
(with Ω = [a, b] and f (u) > 0 ∀x ∈ Ω).

Discontinuous Galerkin-space discretization

In the DG method, the domain Ω = [a, b] is subdivided into a union of finite number N el of cells {D k } N el k=1 , each element being delimited by two nodes of coordinates x l k and x r k , such that:

Ω = [a, b] = N el k=1 D k with D k D l = ∀1 ≤ k = l ≤ N el
We note that: x l 1 = a and x r N el = b.

We define uniform N el cells D k of length h by:

D k = [x l k , x r k ] and h = x r k -x l k , ∀1 ≤ k ≤ N el .
We denote by T the subdivision of Ω into N el elements D k ,

T = D k , 1 ≤ k ≤ N el .
The starting point for the derivation of a DG scheme is obtained by multiplying Eq. ( 7.1) by a polynomial test function v on each cell and then integrating over an arbitrary subset D k of T :

D k ∂u(x, t ) ∂t v(x)d x + D k ∂ f (u(x, t )) ∂x v(x)d x = 0, ∀t ∈ [0, T ].
We underline that this integration is carried-out over a single element D k and not over the whole computational domain Ω, as done in a classical FE approach. After integration by parts, one obtains:

D k ∂u(x, t ) ∂t v(x)d x = D k f (u(x, t )) ∂v(x) ∂x d x + f u(x l k , t ) -f u(x r k , t ) , ∀t ∈ [0, T ].
The DG method represents the unknowns like the FEM by piecewise polynomial functions, but unlike FEM the polynomials are discontinuous at the cell interfaces. A numerical flux is defined at the cell interface in the same way as for Finite Volume (FV) methods. So, on each cell D k , the discrete unknown u k h is represented as a linear combination of well chosen basis functions of the space of polynomials of degree p. Then, the finite-dimensional subspace V p h is defined as:

V p h = v ∈ L 2 (Ω) | v | D k ∈ P p (D k ) ∀1 ≤ k ≤ N el , D k ∈ T ,
where P p (D k ) represents the space of polynomials of degree up to p defined on the element D k . Notice that functions in V p h are discontinuous across cell interfaces.

As a consequence, the flux f evaluated at the interfaces of the element D k is replaced by a numerical flux function f , due to the fact that u h is a priori discontinuous at the interfaces. Therefore, the local approximate solution u k h is then determined as the unique solution of the following weak formulation:

For each element D k ∈ T :

D k ∂u k h (x, t ) ∂t v h (x)d x = D k f u k h (x, t ) ∂v h (x) ∂x d x + f l k -f r k , ∀t ∈ [0, T ], ∀v h ∈ V p h . (7.2)
Note that, for the boundary element D 1 , the numerical flux for the left edge is defined using the given boundary condition:

f l 1 = f (u a (t )).
For the element D N el , the numerical flux for the right edge is evaluated using interior solution:

f r N el = f (u N el h (t )).
For interior interfaces, a basic idea is to use a numerical flux f that is defined according to the left and right values of the solution at the interface, as in a classical FV approach:

f l k = f u k-1 h (x l - k , t ), u k h (x l + k , t ) , f r k = f u k h (x r - k , t ), u k+1 h (x r + k , t ) , with, u k h = u h | D k .
Naturally, the choice of the flux is important.

Numerical flux

To complete the definition of the approximate solution u h , it only remains to choose the numerical flux f . An approximation to the true flux [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF] at an element's interface is defined by considering the following conditions:

• A numerical flux is defined using interface solution values regardless of the polynomial space chosen for the solution,

f = f u k + h , u k - h ,
where u k + h and u k - h are the fields on the two sides of the element interface Γ k .

• Consistency: For continuous solutions, the numerical flux must be equivalent to the normal flux at the interface:

f (u k h , u k h ) = f (u k h ).
• Conservation: If a piecewise constant approximation is used, the discretization results in a monotone FV scheme. This is ensured if we have a conservative flux,

f u k + h , u k - h -f u k - h , u k + h = 0.
There are many possible choices of the numerical flux [START_REF] Proft | Analytical and numerical study of diffusive fluxes for transport equations with near-degenerate coefficients[END_REF] satisfying the above properties:

Central flux:

f C E N = 1 2 f (u k- h ) + f (u k+ h ) .
Upwind flux:

f U p =          f (u k- h ) if f (u k h ) ≥ 0, f (u k+ h ) else.
Lax-Friedrichs flux:

f LF = 1 2 f (u k- h ) + f (u k+ h ) + C LF 2 u k- h -u k+ h = f C E N + C LF 2 u k- h -u k+ h ,
where C LF is the maximum wave speed at the interface:

C LF = max | f (u h ) | for min(u k- h , u k+ h ) ≤ u h ≤ max(u k- h , u k+ h ).
The Lax-Friedrichs flux adds an extra diffusive term to the central flux in an attempt to smear out instabilities. We will focus on the Lax-Friedrichs flux, in the following work, which is a classical choice in DG methods.

Elementary linear system

Returning now to Eq. (7.2) and substituting the local solution approximation in element D k , there are (p + 1) equations to be solved for each component of the field corresponding to the (p + 1) degrees of freedom. Indeed, if a local basis of P p (D k ) is chosen and denoted as ϕ k j (x) j =1,...,p+1 , for x ∈ D k , we can express the local numerical solution u k h as:

u k h (x, t ) = u h | D k (x, t ) = p+1 j =1 u k j (t )ϕ k j (x), ∀x ∈ D k , ∀t ∈ [0, T ].
If Lagrange polynomials are used as local basis, u k j corresponds to the solution value in element k at interpolation point j .

Since the support of the basis functions is restricted to the element D k , the left-hand side of equation 7.2 is simply the product of a local mass matrix M k of size (p + 1) × (p + 1) with the vector of the time-derivative of the degrees of freedom ∂ t u k for the element D k . The right-hand side represents the residual for the element D k , composed of a volume integral R k and fluxes at the interfaces F l and F r , which ensures the coupling with solution in neighboring elements.

Therefore, we should solve the following local system for each element k:

M k ∂ t u k = R k (u k ) + f l k -f r k ∀t ∈ [0, T ] k = 1, ..., N el , (7.3)
for the coefficients:

u k =          u k 1 u k 2 . . . u k p+1          .
The global solution u(x, t ) is then assumed to be approximated by the piecewise p-th order polynomial approximation u h (x, t ):

u(x, t ) u h (x, t ) = N el k=1 u k h (x, t ),
defined as the direct sum of the N el local polynomial solutions u k h (x, t ).

Computation of residual and mass matrix

Contrary to the FEM, only local matrices on each element (in practice only the elementary matrices on the reference interval [-1, 1]) need to be assembled. The coefficients of the local mass matrix for the element D k are:

(M k ) i j = D k ϕ k i (x)ϕ k j (x)d x i , j = 1, ..., p + 1.
The first contribution to the residual is the integral over each element D k of the flux multiplied by the test function gradient:

R k (u k ) i j = D k f u k h (x, t ) j ∂ x ϕ k i (x)d x i , j = 1, ..., p + 1.
The numerical integration is performed using the Gauss-Legendre quadrature rules described above.

CFL condition for DG Method

An important feature of the above-presented method is the use of a local mass matrix, which results from the local support of basis functions. As a consequence, the local mass matrix can be inverted easily before time integration and the method is well suited to high-order explicit time integration, like RK methods, and is highly parallelizable. Therefore, in the present work, RK 2 and RK 4 are used for time integration.

Because we are focusing on DG schemes, we discuss the limitation on the C C F L number when the DG method is used in conjunction with the RK time integration approach. A stability condition on the size of the timestep must indeed be satisfied. It corresponds to the Courant Friedrichs-Lewy (CFL) condition:

λ ∆t h ≤ C C F L , where λ = max u | f (u) |, h
is the smallest element width, and ∆t is the length of the time step.

Physically this condition bounds the size of the timestep to ensure the physical features of the solution are resolved over the mesh. It is also noteworthy to realise that numerical stability is ensured by bounding the CFL number by (2p + 1) -1 , i.e.:

λ ∆t h ≤ 1 2p + 1 ,
where p is the degree of the approximating polynomial (when polynomial of degree p is used, a RK of order (p + 1) must be used to recover optimal convergence rate). This condition has been proven for the polynomial order p = 0 and p = 1, there is no analytical proof for higher order polynomials as for as we know [START_REF] Cockburn | The Runge-Kutta local projection P 1 -discontinuous-Galerkin finite element method for scalar conservation laws[END_REF].

Numerical results

We study the DGFE method in conjunction with RK for the time disretization for the one-dimensional advection problem:

   ∂u(x,t ) ∂t + c ∂u(x,t ) ∂x = 0 ∀(x, t ) ∈ [-1, 1] × [0, T ] c > 0, u(x, 0) = u 0 (x) = sin(2πx) ∀x ∈ [-1, 1], (7.4)
and left boundary condition,

u(-1, t ) = sin(2π(-1 -c t )).
The effect of 1D advection is thus to move an initial distribution with speed c > 0. To complete the numerical scheme, the Lax-Friedrichs flux was chosen, where C LF = c. We get:

f LF (u k- h , u k+ h ) = 1 2 f (u k- h ) + f (u k h ) + C LF 2 u k- h -u k+ h = cu k- h ∀2 ≤ k ≤ N el -1.
which is equivalent to the upwind flux in this case.

Setting c = 1, an exact solution can be found,

u ex (x, t ) = sin(2π(x -t )).
To test the validity of the classical DG method in conjunction with the RK method for the time discretisation, the L 2 -error of the difference between the numerical and exact solution across the domain was computed.

Hesthaven et al. [START_REF] Hesthaven | Nodal discontinuous Galerkin methods: algorithms, analysis, and applications[END_REF] show that at time T the error measured in the L 2 -norm should be of the form:

∥ u ex (T ) -u h (T ) ∥ L 2 (Ω) ≤ C (N p )(h) N p (1 +C 1 (N p )T ),
where h is the smallest elements size, N p (= p+1) is the number of nodes per element and C , C 1 are constants independent of h.

Table 7.1 displays the L 2 -error as a function of the number of elements N el , and polynomial degree, p = N p -1.

Isogeometric -discontinuous Galerkin framework (IGDG)

In this section, we present a method that combines isogeometric analysis (IGA) with the discontinuous Galerkin (DG) method for solving hyperbolic equations. The basis functions are continuous within each patch, and discontinuous only on patch boundaries. We have already mentioned that the IGA space is defined using patches rather than elements as in FEA spaces. Therefore, the DG application in IGA is a patch-to-patch relation instead of an element-to-element one. This fact is important to remember, since every time we refer to partitions in the domain, we are in fact referring to patches and not to elements.

As in the previous part of the present chapter we have already seen the technical part of DGFEM, it is now sufficient to focus on the necessary adjustments for IGDG.

Construction of the DG basis

In order to apply the IGDG methodology for the problem (5.1), we assume that the physical domain Ω is represented by a set of non overlapping B-splines patches D l :

Ω = l D l • D l • D l = ∀ l = l .
As it is common in the IGA analysis, we assume a parametric domain D l composed of the knots (ξ 1 , ..., ξ n+p+1 ).

We propose to consider as DG patches the intervals delimited in the parametric domain D l by all the knots (ξ 1 , ..., ξ n+p+1 ).

Therefore, in the framework of the DG method, we can define a set of "elements" delimited by the interfaces

x l k and x r k (assuming p + 1 equal knots at domain extremities) for each patch D l :

x l k = x(ξ k+p ) and x r k = x(ξ k+p+1 ), k = 1, ..., n -p = N el .
Thus, this baseline discretization fits the definition of the computational domain. If a finer discretization is required, additional knots can be inserted locally in the representation, without any modification of the geometry, using the knot-insertion procedure described in chapter 3.

To complete the construction of the DG framework, it only remains to define for each interval

D k = [x l k , x r k ] the basis functions (Φ k j ) j =1,.
..,p+1 . The B-spline basis functions cannot be used directly, because they do not exhibit discontinuities at the interfaces. A modification of the basis should therefore be achieved before.

A straightforward approach consists in using again the knot-insertion procedure p times, for each of the existing interior knots (ξ p+2 , ..., ξ n ) as illustrated in section 3.7.

By doing so, the computational domain is split into a set of Bézier patches, without modification of the geometry. A Bézier patch is a particular case of B-spline patch, for which the number n of functions (and control points) is equal to p +1. As a consequence, p +1 basis functions (Φ k j (x)) j =1,...,p+1 are defined in each interval D k , which can be identified with Bernstein polynomials of degree p in the parametric domain (with a change of parameter from [0, 1] to [ξ k+p , ξ k+p+1 ]).

The resulting representation has several suitable properties:

• (Φ k j ) j =1,...,p+1 are non-zero for x ∈]x(ξ k+p ), x(ξ k+p+1 )[. • (Φ k j ) j =1,...,p+1 are C ∞ for x ∈]x(ξ k+p ), x(ξ k+p+1 )[. • Φ k
1 and Φ k p+1 are equal to one at x + (ξ k+p ) and x -(ξ k+p+1 ).

• p+1 j =1 Φ k j (x) = 1 ∀x ∈]x(ξ k+p ), x(ξ k+p )[ (partition of unity). • u h | D k (x + (ξ k+p )) = u k 1 (local solution interpolates first degree of freedom). • u h | D k (x -(ξ k+p+1 )) = u k p+1 (local solution interpolates last degree of freedom). • min j u k j ≤ u h | D k ≤ max j u k j (convex hull property)
. These properties will be exploited to define an efficient numerical procedure, as described below.

Isogeometric discontinuous Galerkin approximation spaces

Let us now consider the physical element D k , that can be exactly parametrized with a mapping T:

T : D k -→ D k , ξ -→ x(ξ).
Thus, any point of coordinate x in the physical domain D k is mapped to a point of parameter ξ in the parametric domain D k . The transformation is defined by associating a control point to each basis function:

x(ξ) = p+1 j =1
B k j ,p (ξ)P j .

P j , for j = 1, ..., p + 1 are defined by the knot insertion procedure as described above.

We will revisit the one-dimensional advection problem that is given in (5.1). Using the Bernstein basis functions in the element D k , u h | D k is described as:

u h | D k (x, t ) = u k h (x, t ) = p+1 j =1 B k j ,p (x)u k j (t ).
For the derivation of a IGDG scheme, equation (5.1) is multiplied by a polynomial test function on each local patch D k , ∀1 ≤ k ≤ N el which can be choosen equal to the Bernstein basis function and after integration by parts and introducing the numerical flux f , we get:

p+1 j =1 ∂ t u k j (t ) D k B k i ,p (x)B k j ,p (x)d x = c p+1 j =1 u k j (t ) D k B k j ,p (x)∂ x B k i ,p (x)d x+ f (x l k , t )-f (x r k , t )∀t ∈ [0, T ] k = 1, ..., N el .
Since the support of the basis functions is restricted to the element D k , the left-hand side of the previous equation is simply the product of a local mass matrix M k of size (p + 1) × (p + 1) with the vector of the timederivative of the degrees of freedom ∂ t u k for the element D k . The right-hand side represents the residual for the element D k , composed of a volume integral R k and fluxes at the interfaces f l k and f r k , which ensures the coupling with solution in neighboring elements. Therefore, the problem can be expressed in local matrix form:

M k ∂ t u k = R k u k + f l k -f r k ∀t ∈ [0, T ] k = 1, ..., N el . (7.7)
As for the DGFEM, for the first and last elements, the boundary fluxes f (u a ) and f (u b ) replace f l 1 and f r N el respectively. Note that, if one uses a piecewise constant representation p = 0, the mass matrix is just a scalar equal to the element volume, the volumic residual is zero and one recovers a classical first-order FV method.

Computation of residual and mass matrix

The coefficients of the local mass matrix and stiffness matrix for the element D k are:

(M k ) i j = D k B k i ,p (x)B k j ,p (x)d x i , j = 1, ..., p + 1, (R k ) i j = c D k B k j ,p (x)∂ x B k i ,p (x)d x i , j = 1, ..., p + 1.
By performing the integration in the parametric cell D k , we obtain:

(M k ) i j = D k B k i ,p (ξ)B k j ,p (ξ) J k ξ d ξ i , j = 1, ..., p + 1, and
(R k ) i j = c D k B k j ,p (ξ)∂ ξ B k i ,p (ξ)d ξ i , j = 1, ..., p + 1,
where we denote

J k ξ = ∂x ∂ξ | D k .
The computation is achieved using Gauss-Legendre quadrature rules:

(M k ) i j = n G l =1 B k i ,p X G (l ) B k j ,p X G (l ) ω G l J k X G (l ) i , j = 1, ..., p + 1,
where (X G (l )) l =1,...,n G and (ω G l ) l =1,...,n G are quadrature abscissae and weights. Its inverse is computed numerically, in a pre-processing phase.

(R k ) i j = c n G l =1 B k j ,p X G (l ) ∂ ξB k i ,p X G (l ) ω G l i , j = 1, ..., p + 1.
The solution in each element D k interpolates the local first and last degrees of freedom u k 1 and u k p+1 . Therefore, the flux computation only depends on the two degrees of freedom located at each interface:

f (x l k , t ) = f u k-1 p+1 , u k 1 , f (x r k , t ) = f u k p+1 , u k+1 1 .
The solutions obtained numerically are compared with the exact solutions in terms of L 2 norm, for basic functions of degrees 0, 1, 2, 3 and 4. The plots of the numerical solutions obtained using only 4 elements and a quadratic, cubic and quartic basis are depicted in Figure 7.3.

We underline the accuracy of the solution obtained with a so small number of elements. Note also that the discontinuities at the interfaces are decreasing as p increases.

One method to test the convergence of a spatial discretisation is by computing the error between the numerical solution and the analytical one. The numerical approximation error in L 2 -norm:

∥ e ∥ L 2 (Ω) =∥ u ex (T ) -u h (T ) ∥ L 2 (Ω) ,
for various element sizes h and degrees p of the Bernstein basis and it is depicted in Fig. 7.3. The convergence rates are shown in Table 7.2 and Table 7.

3. An optimal convergence rate is observed, the method being of order p + 1 with respect to L 2norm.

We notice that, as expected, the convergence rate is limited when using RK 2 time integrator.

Conclusion and comparison

In this chapter a new family of discontinuous Galerkin methods which combines the IGA with the DG method, called IGDG method, has been developed for the one-dimensional advection problem. Our method takes advantage of both IGA and DG methods. In fact, DG methodology is adopted at Bézier patch level, i.e., we employ the traditional IGA within each Bézier patch, and employ the DG method across the patch interfaces to glue the multiple patches. Bézier patches, considered as elements, are constructed by transformation of the initial B-spline domain. Due to IGA, all conic sections can be represented exactly, thus eliminating the geometrical errors by the construction. Obviously, this property will be more evident for 2D problems.

As mentioned before, the major reason for using DG methods lies with their ability to provide stable numerical methods for first order PDE problems, for which classical FEM is well known to perform poorly. Due to the piecewise discontinuity of basis functions, the DG method can be applied locally in each element. This simplifies the implementation of the method, since the mass matrix becomes block diagonal and the solution of a large system is avoided. The solution for the whole computational domain is achieved by summing over all the elements of the mesh.

Compared to SUPG-FEM for solving hyperbolic problems, an attractive feature of the DG method is that it still uses the basic Galerkin method for the volume integrals, although this is extended by element boundary integrals to achieve the 'upwinding'. Because of the simple concept of the method it is basically the same for any space dimension. The treatment of boundary conditions is also relatively easy.

Now we focus our attention on the comparison of the numerical methods combining IGFEM, DGFEM and IGDG space discretization and explicit RK time integration, for the one-dimensional advection problem. I N the previous chapter, we have described the various computational procedures for IGDG analysis especially for the one-dimensional case. This chapter is devoted to the extension to the two-dimensional case. The geometrical representation will obviously play a more critical role in 2D cases and we will especially underline the tretment of the geometry. Numerical experiments validate the presented methodology. Furthermore, the boundary of the domain Ω is denoted by ∂Ω.

Computational procedures in two dimensions

Isogeometric analysis (IGA): physical domain and geometrical mappings

The main idea behind the isogeometric approach is to discretize the unknowns of the problem with the same set of basis functions that CAD employs for the construction of the geometries.

COMPUTATIONAL PROCEDURES IN TWO DIMENSIONS

We use the definitions from IGA that we introduced in chapter 6. In our specific case the physical domain Note that this transformation is non linear and is simply defined by the parametric representation of the computational domain. Unless otherwise specified, we assume this mapping to be invertible. Its inverse:

Ω ⊂ R 2 is
T -1 : D e -→ D e ,
takes points in the physical domain to their corresponding parameter values. Figure 8.1 gives a schematic overview of our proposed approach.

Basic function space for the parametric domain and physical domain

In this section we focus on the construction of the basis functions. A C p-1 continuous basis Φ(ξ, η) is constructed over the parameteric domain D e . The basis is used to construct a geometric map T so that it maps a point (ξ, η) ∈ R 2 in the parametric domain D e to a point in D e .

We have already stated that we consider elements to be the images of knot spans under the Bézier mapping.

We will denote these knot spans in the parameter space by D e , and their image in the physical space as D e , where e = 1, ..., N el , with N el being the total number of elements.

To implement the DG method in the isogeometric framework, i.e. based on a computational domain defined from a B-spline representation, we must first define a set of elements, which are the supports of a polynomial representation with discontinuities at each interface between elements. Given a B-spline surface defining the computational domain, the insertions of knots p times are used for each of the existing interior knot for each parametric direction sequentially. In doing so, the computational domain is divided into a set of Bézier patches, without modification of the geometry. We remind that a Bézier patch is a special case of B-spline patch, for which the number n of functions (and control points) is equal to p + 1 (with p is the degree of basis function).

Finally, all the Bézier patches created by the insertion process are considered as elements. Each element D e is therefore defined by (p + 1) × (q + 1) basis functions, B p,q (x, y) e , which can be identified with Bernstein polynomials of degrees p and q:

B p,q (X ) q j e j =1,...,q+1 are respectively the i -th and the j -th Bernstein polynomials of degree p and q, defined on the interval [0, 1] by:

B p i (ξ) e = C i p ξ i 1 -ξ p-i ∀t ∈ [0, 1] i = 0, ..., p.
A given control point P e will have local indices associated to the Bézier patch D e . Thus, the geometrical mapping local to element e can be defined as:

X e (ξ, η) = x e (ξ, η), y e (ξ, η) = p+1 i =1 q+1 j =1 B p i (ξ) e B q j (η) e P e i j .
Before we proceed further, we denote by J e the elementary Jacobian matrix of this transformation defined in D e by: Thus, we can also calculate the inverse of the elementary Jacobian matrix (J e ) -1 given by: Then, we define the test functions in the physical domain D e by using the same representation as the geometry:

Φ p,q (x, y) | D e = Φ p,q (x, y) e = B p,q (x, y) e = B p,q T(ξ, η) e = B p,q (ξ, η)

e = B p (ξ) e ⊗ B q (η) e = Φ p (ξ) e ⊗ Φ q (η) e .
We also highlight that, ∀1 ≤ i ≤ p + 1 and ∀1 ≤ j ≤ q + 1,

Φ p,q (x, y) e i , j = B p (ξ) e i B q (η) e j = Φ p (ξ) e i Φ q (η) e j .
In addition, we need also to compute ∇Φ e , which is the gradient in the Cartesian form:

∇Φ e =      ∂Φ p,q (x,

Numerical integration

In this thesis we use Gaussian quadrature as a choice of numerical integration. To do this, we will need to map our functions over to a reference square ( ξ, η)

∈ D = [-1, 1] × [-1, 1].
While this is common for most FEM as well, we have the additional mapping from the parametric space to the physical space. The first one is an affine mapping (see Fig. 8.1):

T : D -→ D e .
The numerical integration has to take place in the reference square, while all basis functions are defined over the parametric domain and the differential equation is formulated in the physical space. Therefore, a typical integration writes:

2D advection problem: IGDG space semi-discretization

In the following, we describe the discretization of the transient, bi-dimensional, linear conservation law or transport equation over a domain Ω ⊂ R 2 with periodic boundary conditions by the IGDG method:

   ∂ t u + ∇ • - → c u(x, y, t ) = 0 ∀(x, y, t ) ∈ Ω × [0, T ], u(x, y, 0) = u 0 (x, y) ∀(x, y) ∈ Ω, (8.1)
where u(x, y, t ) is a scalar quantity transported by a continuous velocity field -→ c = (c x , c y ) t .

Similarly as in the previous chapter, we derive the IGDG space semidiscretization leading to a system of ordinary differential equations. As a matter of fact, the IGDG method for multi-dimensional conservation law has the same structure it has for one-dimensional scalar conservation laws, we only need to describe the DG-space discretization.

Applying a IGDG method, the solution u is approximated by u h ∈ V p , which we assume to have the following form:

u e h (x, y, t ) = p+1 i =1 q+1 j =1 B i ,p (ξ) e B j ,q (η) e u e i j (t ),
where u e i j : [0, T ] -→ R, ∀ 1 ≤ i ≤ p + 1, 1 ≤ j ≤ q + 1 are local unknown coefficients. Before we proceed further, let us put v e h ∈ V p h , the approximate solution that will be sought for each t ∈ [0, T ] in the finite-dimensional space:

V p h = v := v(x(ξ, η), y(ξ, η)) ∈ L 2 (Ω), v | D e ∈ span B e ,
where B e is the set of (p + 1) × (q + 1) Bernstein polynomials defined over D e .

2D ADVECTION PROBLEM: IGDG SPACE SEMI-DISCRETIZATION

Isogeometric discontinuous Galerkin space semi-discretization

The starting point for a DG discretization is the weak formulation, which is obtained by multiplying Eq. (8.1) by a local arbitrary test function v e (x, y) ∈ V p , and then integrating on each patch D e separately. Note that, in this framework, no continuity on the state vector u e and the test function v e is enforced along the interfaces between patches. The variational formulation is given for e = 1, ..., N el by: Therefore, the local problem takes the form of a linear system of size (p +1) 2 ×(q +1) 2 , which can be written in the following matrix form:

M e ∂ t u e = R e u e -F e
∀t ∈ [0, T ] ∀1 ≤ e ≤ N el .

(8.5)

Elementary linear system

The coefficients of the local mass and stiffness matrix for the element D e are written as: The computation is achieved using Gauss-Legendre quadrature rules. To complete the scheme (8.5), we choose the Lax-Friedrichs flux. Obviously, these vectors are normalized for the computations.

The RK time discretization

The space semidiscrete problem (8.5) represents a system of ordinary differential equations (ODEs), which has to be solved with the Runge Kutta schemes of order 2 and 4. As mentioned before, for the 1D case the time step ∆t is strongly restricted by the Courant-Friedrichs-Lewy (CFL) stability condition. In the context of a general multi-dimensional hyperbolic system, the combination of RK of order p + 1 with DG at order p, the CFL condition can be written as: In what follows, we consider a partition 0 = t 0 < t 1 < ... < t ň = T of the time interval [0, T ] and set ∆t = t k+1t k for k = 1, ..., ň -1.

∆t ≤ C C F L min | d e | ∥ - → c (2p + 1) ∥ , ( 8

Numerical results

In order to demonstrate the performance of the present method, we consider an example of 2D advection problem given above by Eq. (8.1) whose initial solution is given by: u 0 (x, y) = exp -5(x 2 + y 2 ) .

(8.9)

The analytical solution to this problem is:

u analytic (x, y, t ) = exp -5(x -c x t ) 2 -5(y -c y t ) 2 ) .
We present in Fig. taken to be of bi-degree (p, q). We focus for the case p = q which will be specified in each example. Before we proceed further, let us recall that we denote by ∥ u ∥ L 2 (Ω) , the L 2 (Ω)-norm of a function u:

∥ u ∥ L 2 (Ω) = Ω u 2 d Ω 1 2 .
We measure the convergence of the numerical methods in the L 2 -norm for the Cartesian, linear and curvilinear grids. Our aim is to identify the numerical order of convergence. We define the error of the numerical solution as:

e h = u h -u ex , (8.10)
and its L 2 -norm is:

∥ e h ∥ L 2 (Ω) = Ω u h (x, y) -u ex (x, y) 2 d Ω 1/2 . ( 8

Linear grids

Linear patches are still employed in this case. The initial configurations of N el = N 1 el × N 2 el patches and degrees (p, q) of bivariate Bernstein function are the same with those described in the previous case. Convergence results for the 2D advection problem with different choices of N el are given in Tab. (8.3) for the linear and quadratic Bernstein and in Tab. (8.4) Table 8.5 and table 8.6 show the quality of the numerical approximations in L 2 norm for various size of element h and order of polynomial (p and q). The error of the numerical approximations is depicted in Figure 8.12. The convergence rates are 2, 3, 4 and 5 using bases of degrees p = q = 1, p = q = 2, p = q = 3 and p = q = 4, respectively.

Conclusion

As it has been pointed out previously, we can conclude that IGDG is a powerful tool. In this chapter, we have introduced a new analysis framework, called isogeometric discontinuous Galerkin (IGDG) method for bi-dimensional hyperbolic problem, which is based on Bézier extraction. We have confined our attention to the development of the RK-IGDG methods for bi-dimensional advection problem. The resulting IGDG method in conjunction with RK method is stable, high-order accurate, and highly parallelizable. The scheme can easily handle complicated geometries and boundary conditions. The flexibility of the method to handle different geometries and to work with different elements has been shown.

As consequence, this method can be easily formulated and implemented. Thus, the numerical behavior of the method is evaluated and it has shown an optimal convergence rate.

For (x, y, t ) ∈ Ω × 0, T , the linear acoustic wave equation can be formulated as a first-order hyperbolic system in terms of the pressure field p and velocity field -→ U = u, v :

         ∂ -→ U ∂t + ∇p = 0 (x, y, t ) ∈ Ω × [0, T ], ∂p ∂t + ∇ • -→ U = 0 (x, y, t ) ∈ Ω × [0, T ], -→ U • - → n = 0 on ∂Ω × [0, T ]. (9.1)
Here, p represent the pressure perturbation and -→ U the velocity perturbation with respect to a reference state at rest. Moreover, the equations have been adimensionalized by assuming that the sound speed is equal to unity.

At boundary, we assume perfect wall conditions, all waves being ideally reflected.

We supplement the system (9.1) with the initial condition:

         u(x, y, 0) = u 0 (x, y) (x, y) ∈ Ω, v(x, y, 0) = v 0 (x, y) (x, y) ∈ Ω, p(x, y, 0) = p 0 (x, y) (x, y) ∈ Ω. (9.2)
An equivalent formulation of the system (9.1) is the following:

               ∂u ∂t + ∂p ∂x = 0 (x, y, t ) ∈ Ω × [0, T ], (1)
∂v ∂t + ∂p ∂y = 0 (x, y, t ) ∈ Ω × [0, T ], (2)
∂p ∂t + ∂u ∂x + ∂v ∂y = 0 (x, y, t ) ∈ Ω × [0, T ], (3)
-→ U • - → n = 0 on ∂Ω × [0, T ]. (9.3)

IGDG approximation of the acoustic wave equations 9.2.1 Spatial discretization

This section introduces the weak formulation of the linear acoustic wave equations given above by equation (9.3) using IGDG method. As highlighted in the previous chapter, in DG method the basis functions are allowed to be discontinuous at the element boundaries, the integrals are performed element-wise, and the compling of the wave field across the elements is weakly imposed using fluxes. The first step is to derive the weak formulations of equation (9.3). This weak formulation is then discretized by introducing an approximation for the pressure and velocity in a finite-dimensional subspace to obtain a linear system of ordinary differential equations.

In order to apply the IGA methodology, we highlight that the physical domain Ω is subdivided into Bspline patches D e , S (Ω) := D e N pa e=1 .

We suppose that the computational domain Ω can be exactly represented by the union of non-overlapping patches D e :

Ω = D e with • D e • D l = ∀ e = l .
Each patch D e defined by:

D e = X e = (x e , y e ) ∈ R 2 | X e = T(ξ, η) such that (ξ, η) ∈ D e ,
where the transformation T is defined for all (ξ, η) ∈ D e by:

T : D e -→ D e (ξ, η) -→ x e (ξ, η), y e (ξ, η) .

As in the previous chapter, a set of Bézier elements is obtained from the B-spline patches by multiple knot insertion. Within each element D e we assume that the local solution is well approximated by twodimensional Bernstein polynomials of degrees p and q:

     u h (x, y, t ) v h (x, y, t ) p h (x, y, t )      e = p+1 i =1 q+1 j =1 Φ p,q i , j (x, y) e      u i j (t ) v i j (t ) p i j (t )      e = p+1 i =1 q+1 j =1 B p,q i , j (x, y) e      u i j (t ) v i j (t ) p i j (t )      e = p+1 i =1 q+1 j =1 B p,q i , j (T(ξ, η)) e      u i j (t ) v i j (t ) p i j (t )      e = p+1 i =1 q+1 j =1 B i ,p (ξ) e B j ,q (η) e      u i j (t ) v i j (t ) p i j (t )      e
where u e i j : [0, T ] -→ R, v e i j : [0, T ] -→ R and p e i j : [0, T ] -→ R, ∀1 ≤ i ≤ p + 1, ∀1 ≤ j ≤ q + 1 are local unknown coefficients.

The number of degree of freedom inside the element D e is (p + 1) × (q + 1). Φ p,q i , j e are the bivariate Bernstein polynomials of degree p × q. The global solution can then be approximated by the direct sum of the local solutions:

     u(x, y, t ) v(x, y, t ) p(x, y, t )           u h (x, y, t ) v h (x, y, t ) p h (x, y, t )      = N el e=1      u h (x, y, t ) v h (x, y, t ) p h (x, y, t )      e .
Before we proceed further, let us put w e h ∈ V p h , the approximate solution will be sought for each t ∈ [0, T ] in the finite-dimensional space:

V p h = w := w x(ξ, η), y(ξ, η) ∈ L 2 (Ω), w | D e ∈ span B e ,
where B e is the set of (p + 1) × (q + 1) Bernstein polynomials defined over D e .

The starting point for a DG discretization is the weak formulation, which is obtained by multiplying each equations of the system (9.3) by a local arbitrary test function w e (x, y) ∈ V p , and then integrating on each Bézier patch D e separately. Note that, in this framework, no continuity on the state vectors u e , v e , p e and the test function w e is enforced along the interfaces between Bézier patches. The weak formulation of the acoustic wave equation is given for e = 1, ..., N el for each equations of the system (9.3) by the following statement:

First variational equation

Applying a IGDG method to equation (1) of the system (9.4), the weak form of the problem can be written for each element D e as: By discretizing the problem on the Bézier basis associated with the element and using B e k,l as test function, the problem is written:

p+1 i =1 q+1 j =1 ∂ t u e i , j (t ) D e B e i , j (x, y)B e k,l (x, y)d D e = p+1 i =1 q+1 j =1 p e i , j (t ) D e B e i , j (x, y) - → e x • ∇B e k,l (x, y)d D e - Γ e B e k,l (x, y) - → e x p e • - → n e d Γ e ∀1 ≤ k ≤ p + 1, ∀1 ≤ l ≤ q + 1.
Therefore, the local problem takes the form of a linear system of size (p +1) 2 ×(q +1) 2 , which can be written in the following matrix form:

M e ∂ t u e = R e x p e -F e x ∀t ∈ [0, T ], ∀1 ≤ e ≤ N el .
(9.5)

Second variational equation

In the same way, the weak form of the equation (2) of the system (9.4) can be described for each element D e as follows: By discretizing the problem on the Bézier basis associated with the element and using B e k,l as test function, the problem is written:

p+1 i =1 q+1 j =1 ∂ t v e i , j (t ) D e B e i , j (x, y)B e k,l (x, y)d D e = p+1 i =1 q+1 j =1 p e i , j (t ) D e B e i , j (x, y) - → e y • ∇B e k,l (x, y)d D e - Γ e B e k,l (x, y) - → e y p e • - → n e d Γ e , ∀1 ≤ k ≤ p + 1, ∀1 ≤ l ≤ q + 1.
Therefore, the local problem takes the form of a linear system of size (p + 1) 2 × (q + 1) 2 , which can be written in the following matrix form:

M e ∂ t v e =
≤ p + 1, ∀1 ≤ l ≤ q + 1.
Therefore, this local problem leads to a linear system of size (p + 1) 2 × (q + 1) 2 , which can be written in the following matrix form: (9.9)

M e ∂ t p e =

Elementary linear system

The coefficients of the local mass and stiffness matrix for the element D e are written as: Achieving the integration in the unit square D, we get:

M e kl ,i j = D B i ,p ( ξ) e B j ,q (η) e B k,p ( ξ) e B l ,q (η) e | Ĵ e (ξ, η) || J e (ξ, η) | d D, ∀i , k = 1, ..., p + 1 j , l = 1, ..., q + 1.
The integrals are evaluated numerically via Gaussian quadrature. = R e x kl ,i j + R e y kl ,i j .

By expressing the integration in the local parametric space D e , we get: By performing the integration in the reference square D, we get:

R e x kl ,i j = D B i ,p ( ξ) e B j ,q (η) e ∂B e k,p ( ξ) ∂ ξ B e l ,q (η) ∂ ξ ∂x + B e k,p ( ξ) ∂B e l ,q (η) ∂η ∂η ∂x | Ĵ e ( ξ, η)| | J e (ξ, η)| d D, R e y kl ,i j = D B i ,p ( ξ) e B j ,q (η) e ∂B e k,p ( ξ) ∂ ξ B e l ,q (η) ∂ ξ ∂y + B e k,p ( ξ) ∂B e l ,q (η) ∂η ∂η ∂y | Ĵ e ( ξ, η)| | J e (ξ, η)| d D.
At the end, the computation is achieved using Gauss-Legendre quadrature rules.

Numerical Lax-Friedrichs fluxes

The selection of the numerical flux is of utmost importance when formulating a DG method. The numerical flux must accurately couple the neighboring elements, while yielding a stable scheme. There are several flux methods. In our formulation, we only discuss the Lax-Friedrich flux [START_REF] Qiu | A numerical study for the performance of the Runge-Kutta discontinuous Galerkin method based on different numerical fluxes[END_REF], given by:

f e n (u l , u r ) = 1 2 f n (u e l ) + f n (u e r ) - ρ 2 u e r -u e l , (9.10)
where u l and u r are the left and right limits of the discontinuous solution u h . For the Lax-Friedrichs flux, ρ is taken as an upper bound for | f (u) | in the scalar case, or for the absolute value of eigenvalues of the Jacobian for the system case. For the acoustic wave equation, ρ is the sound speed assumed ρ = 1.

Since the weak form of the DG method is written elementwise, the numerical flux between adjacent elements must be defined. For this purpose, it is possible to write: Obviously, these vectors are normalized for the computations. The initial conditions are:

F e = 4 k=1 Γ e k f e n B e k,l (x, y)d Γ e k = [0,1] f e n | Γ e 1 B k,p (ξ) e B l ,q (0) =1 e | J e (ξ, 0) || J e ( ξ, 0) | d ξ + [0,1] f e n | Γ e 2 B k,p (1) =1 e B l ,q (η) e | J e (1, η) || J e (1, η) | d η + [0,1] f e n | Γ e 3 B k,p (ξ) e B l ,q (1)
   p | t =0 = J 0 α 1 + r -r 1 r 2 -r 1 (α 2 -α 1 ) , -→ U | t =0 = u, v | t =0 = (0, 0).
This solution is provided where α 1 = 3.8317, α 2 = 10.1735 (α 1 and α 2 : roots of J 0 = J 1 ) (see appendix C ) , r 1 = 0.25, r 2 = 1 and J 0 the Bessel function of the first kind.

The exact solutions is given in cylindrical coordinates (r, θ) by:

         p ex = J 0 α 1 + r -r 1
r 2 -r 1 (α 2 -α 1 ) cos α 2 -α 1 r 2 -r 1 t , u ex = J 1 α 1 + r -r 1 r 2 -r 1 (α 2 -α 1 ) sin α 2 -α 1 r 2 -r 1 t cos θ , v ex = J 1 α 1 + r -r 1 r 2 -r 1 (α 2 -α 1 ) sin α 2 -α 1 r 2 -r 1 t sin θ . The bivariate Bernstein functions are taken to be of bi-degree (p, q). We focus for the case p = q which will be specified in each example.

Rectilinear grids

The first mesh type considered is a simple rectilinear grid. As can be seen on Figures 9.7 and 9.8, the use of a rectilinear grid strongly impacts the accuracy of the solution near the boundaries of the domain. In particular, one can notice the presence of fictious variations of the solution at each boundary vertex, which are due to the fact that the boundary conditions are not prescribed at the true location and normals are not well approximated. When the rectilinear grid is refined, these effects are reduced but are still dommageable. In Figures 9.9 to 9.11, one can observe that increasing.

The approximation order of the solution only does not provide any remedy to this issue. 9.12.

Figure 9.12: L 2 -error for the 2D acoustic problem using the IGDG method in conjunction with RK 4.

Curvilinear grids

We consider now the case of curvilinear grids. Figure 9.14 depicts the IGDG numerical solutions for the bivariate quadratic Bernstein case for the initial mesh with N el = 4 × 4 patches as shown on the left in Fig. 9.13. We subsequently add a reffinement of ( Contrary to the rectilinear case, the use of curvilinear grids yields a far better accuracy of the solution near the boundaries, as shown in figures 9.14 and 9.15. In particular, one can notice that the boundary isovalues are correctly captured, even on the coarsest grid. Figures 9.16 to 9.18 show that the increase of the approximation order of the solution yields a more accurate solution inside the domain while preserving the boundary accuracy. Note that, in this case, the circle is not exactly represented because B-spline representations are employed, instead of NURBS. Nevertheless, this does not seem to be critical in this example. The impact of the boundary curvature on the solution accuracy, already observed in the solution plots, is confirmed when the convergence rates are computed. For the isogeometric approach using curvilinear grids, quasi-optimal convergence rates close to order p +1 are obtained as the B-spline grids are refined.

However, when using rectilinear grids, a far lower convergence rate is observed, close to second-order (globally), whatever the degree of the solution employed (see Figure 9.12). One can even notice that the accuracy is degraded as the degree of the solution increases, for a fixed number of degrees of freedom. The solver tends to a solution corresponding to the uncorrect boundary geometry (location and normals). Therefore, the use of an accurate geometry description, exact if possible, is mandatory for the efficiency of high-order discretization schemes

Conclusion

In this chapter we gave a description of the ongoing development of a IGDG method to solve the firstorder acoustic wave equation in 2D, expressed as a linear hyperbolic PDE system. As it has been pointed out previously, the computational domain is divided into non-overlapping sub-domains, composed of B-spline patches. The DG approach was applied on element level, each element being a Bézier patch constructed from initial B-spline patches. The solution of the problem is approximated in every element without any continuity requirements for the discrete solution on the interfaces. As a result, IGDG method can be easily formulated and implemented. Finally, the numerical behavior of the method is evaluated and it has shown an optimal convergence rate.

CONCLUSION

D Ans le chapitre 7, nous avons décrit les différentes procédures de calcul pour la méthode de GD dans le cadre de l'AIG, spécialement pour le cas du problème hy- I Sogeometric Analysis can be considered as a Finite Element method that general- izes the set of basis functions from polynomials to B-splines or more generally to non-uniform-rational B-splines. As shown in the literature, this choice guarantees several advantages, from the exact parametrization of geometries defined via CAD to a higher accuracy per-degrees-of-freedom. It also allows for solution fields with higher smoothness. These reasons have made IGA a successful topic in recent years.

In this thesis we applied the IGA method to some hyperbolic problems. We considered standard Galerkin methods as well as stabilized methods, with a special emphasis on B-splines. IGA possesses a set of attractive features from the view point of accuracy and implementational convenience not present in standard FE discretizations. In fact, IGA allows exact representation of a wide class of geometries even on very coarse meshes. In particular, geometry domains having conic sections like circles, cylinders, spheres, ellipsoids, etc... can be represented exactly using NURBS. Refinements can be performed by subdivision of the grid (by inserting knots) or by elevation of the polynomial order of the basis functions in the same way as using the traditional finite element method.

Due to the fact that the aim of IGA is to generalize and improve upon classical FEA, we started this thesis by giving an introduction to IGA by revisiting the original analysis, i.e. FEA, in the context of hyperbolic PDE with SUPG stabilization.

CONCLUSION GÉNÉRALE & PERSPECTIVES

L 'analyse isogéométrique peut être considérée comme une méthode par éléments fi- When Ω = [-1, 1], the roots of the Legendre polynomials and their corresponding weights have been extensively tabulated, so we can simply use these tables without redoing the calculations.

n G Weights ω G Nodes X G 1 2 0 2 -1 3 , 1 3
1, 1 3 -3 5 , 0, The Runge-Kutta method is used to solve a system of ODEs given by: ∂ t u(x, t ) = L t , u(x, t ) .

(B.1)

The s-stage Runge-Kutta (RK) method for (B.1) is written in the form:

∂ t u (1) = L t ň , u ň ∂ t u (i ) = L t ň + c i ∆t , u ň + ∆t i -1 j =1
a i j ∂ t h ( j ) i = 2, ..., s

u ň+1 = u ň + ∆t s i =1 b i ∂ t u (i ) ,
where ∆t = t ň+1t ň is the time step, and u ň and u ň+1 represent the values of u at time t ň and t ň+1 respectively. The coefficients a i j , b i and c i can be summarized in matrix/vector form by the Butcher tableau.

We have a i j = 0 ∀ j ≥ i . Moreover, the coefficients c i and a i j are connected by the condition:

c i = s j =1
a i j i = 1, ..., s.

7. 2

 2 L 2 -error for the IGDG method in conjunction with RK 2 time discretisation for various element sizes and degree of Bézier basis p = 0, 1, 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 L 2 -error for the IGDG method in conjunction with RK 4 time discretisation for various element sizes and degree of Bézier basis p = 2, 3, 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.1 L 2 -error for the 2D advection problem and convergence order for the IGDG method for the linear (left) and quadratic (right) Bernstein bases in conjunction with RK 4 time discretisation. . 8.2 L 2 -error for the 2D advection problem and convergence order for the IGDG method for the cubic (left) and quartic (right) Bernstein bases in conjunction with RK 4 time discretisation. . . . 8.3 L 2 -error for the 2D advection problem and convergence order for the IGDG method for the linear (left) and quadratic (right) Bernstein bases in conjunction with RK 4 time discretisation. . 8.4 L 2 -error for the 2D advection problem and convergence order for the IGDG method for the cubic (left) and quartic (right) Bernstein bases in conjunction with RK 4 time discretisation. . . . 8.5 L 2 -error for the 2D advection problem and convergence order for the IGDG method for the linear (left) and quadratic (right) Bernstein bases in conjunction with RK 4 time discretisation. . 8.6 L 2 -error for the 2D advection problem vs. mesh parameter and convergence order for the IGDG method for the cubic (left) and quartic (right) Bernstein bases in conjunction with RK 4 time discretisation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv List of Tables 9.1 L 2 -error for the 2D acoustic problem and convergence order for the IGDG method for the quadratic Bernstein bases in conjunction with RK 4 time discretisation. . . . . . . . . . . . . . . . . . . . . . 9.2 L 2 -error for the 2D acoustic problem and convergence order for the IGDG method for the cubic (left) and quartic (right) Bernstein bases in conjunction with RK 4 time discretisation. . . . . . . . 9.3 L 2 -error for the 2D acoustic problem and convergence order for the IGDG method for the quadratic Bernstein bases in conjunction with RK 4 time discretisation. . . . . . . . . . . . . . . . . . . . . . 9.4 L 2 -error for the 2D acoustic problem and convergence order for the IGDG method for the cubic (left) and quartic (right) Bernstein bases in conjunction with RK 4 time discretisation. . . . . . . . A.1 Gauss-Legendre nodes and coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv LIST OF FIGURES 2.1 Constant, linear, quadratic and cubic Bernstein polynomials. . . . . . . . . . . . . . . . . . . . . . 2.2 Linear, quadratic and cubic bivariate Bernstein polynomials. . . . . . . . . . . . . . . . . . . . . . 2.3 Bézier curves of various degrees and their control polygons. . . . . . . . . . . . . . . . . . . . . . . 2.4 Bézier curves with endpoint interpolation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Cubic Bézier curve, its control polygon and the convex hull. . . . . . . . . . . . . . . . . . . . . . . 2.6 Quadric Bézier curve and repositioning of the control point P 2 . . . . . . . . . . . . . . . . . . . . . 2.7 Degree elevation of a quadratic Bézier Curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Geometric construction according to De Casteljau's algorithm for p = 3 and ζ = 2/3. . . . . . . . . 2.9 Quadratic rational Bézier curves. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Tensor product Bézier patch of degree 3 × 3 and its control net. . . . . . . . . . . . . . . . . . . . . 3.1 Basis functions of degrees 1, 2 and 3 for uniform knot vector Ξ = {0, 1, 2, 3, ...}. . . . . . . . . . . . . 3.2 Bivariates quadratic and cubic B-spline basis functions [16]. . . . . . . . . . . . . . . . . . . . . . . 3.3 Quadratic basis functions for the open-uniform knot vector Ξ = 0, 0, 0, 1, 2, 3, 3, 3 . . . . . . . . . . 3.4 Quadratic basis functions with reduced continuity at ξ = 1, Ξ = 0, 0, 0, 1, 1, 2, 3, 3, 3 . . . . . . . . . 3.5 Quadratic B-spline curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 A quadratic B-spline curve and its control points. In the right, the curve after moving the control point P 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Before and after knot insertion (cubic B-spline curve). . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 B-spline surface example. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Bézier decomposition (bottom) from a quadratic B-spline basis (top) by knot insertion. . . . . . . 4.1 The least-squares quadratic B-spline curve fitting. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 The least-squares quadratic B-spline Gaussian surface fitting. . . . . . . . . . . . . . . . . . . . . . 4.3 Bessel functions of the first kind-1D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Representation of the physical domain Ω. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 The least-squares quadratic B-spline surface fitting. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Uniform P 1 mesh of [a, b]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Global shape functions for the space V 1 h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 The exact sine wave solution for the one-dimensional advection problem. . . . . . . . . . . . . . . 5.4 Exact and standard Galerkin FEM P 1 solution for the one-dimensional linear advection problem at T = 0.4s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.5 Exact, Galerkin and SUPG solutions at T = 0.4s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.6 SUPG FE P 1 solution for different values of α ∈ [0, 1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.7 Convergence in number of d.o.f. for different choices of α. . . . . . . . . . . . . . . . . . . . . . . . 82 5.8 Uniform P 2 mesh of [a, b]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.9 Global shape functions for the space V 2 h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.10 SUPG quadratic Lagrange P 2 FEM in conjunction with RK 2 for the 1D advection problem. . . . . 88 5.11 L 2 -error for the advection problem with the linear and quadratic Lagrange FEM. . . . . . . . . . . 90 6.1 An example of a B-spline patch in physical space Ω, parametric space Ω, and the reference element Ω used to perform numerical integration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 6.2 (a) SUPG B-spline linear solution for the advection problem (α = 0.1). (b) SUPG FEM P 1 for the advection problem (α = 0.1) at T = 0.4s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 SUPG quadratic B-spline solutions in conjunction with RK 2 for the advection problem. . . . . . 6.4 Convergence rates in the L 2 -norm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Error in the L 2 -norm of IGFEM and classical FEM vs. number of d.o.f. . . . . . . . . . . . . . . . . 7.1 L 2 -errors for the 1D advection problem using the DGFE method in conjunction with the RK method for a sinusoidal initial condition and Lax-Friedrichs flux. . . . . . . . . . . . . . . . . . . . 7.2 Bézier decomposition (bottom) from a quadratic B-spline basis (top) by knot insertion. . . . . . . 7.3 IGDG solution of the 1D advection problem for quadratic (a), cubic (b) and quartic (c) basis and exact solution with 4 Bézier elements. RK 4 time discretization and Lax-Friedrichs flux were used. 7.4 L 2 -errors for the 1D advection problem with a sinusoidal initial condition, RK 2 and RK 4. . . . . 7.5 Convergence rates in IGDG method as a function of the Bernstein function is degree p for the finest grid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.6 Error in the L 2 -norm combining IGFEM, DGFEM, IGDG space discretization and explicit RK time integration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.1 An example of a B-spline patch in physical space, parametric space, and the parent element used to perform numerical integration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.2 Element D e , its faces Γ e k k=1,...,4 and the corresponding normals -→ n e | Γ e k . . . . . . . . . . . . . . . . . 8.3 Analytical solution for the bi-dimensional advection problem, for N el = 4 × 4 at T = 0.5s. . . . . . 8.4 Plots and contour plots of numerical results for bivariate quadratic Bernstein basis with (a) N el = 4 × 4 patches and (b) N el = 8 × 8 patches at T = 0.05s. . . . . . . . . . . . . . . . . . . . . . . . . . . 8.5 IGDG solution for N el = 8 × 8 patches for different degrees (p, q). . . . . . . . . . . . . . . . . . . .8.6 L 2 -error for the 2D advection problem using the IGDG method in conjunction with RK 4. . . . . 8.7 Contour plots of numerical results for bivariate quadratic Bernstein basis at T = 0.05s. . . . . . . 8.8 IGDG solutions for different degrees for N el = 8 × 8 uniform elements. . . . . . . . . . . . . . . . . 8.9 L 2 -error for the 2D advection problem using the IGDG method in conjunction with RK 4. . . . . xviii 8.10 Contour plots of numerical results for bivariate quadratic Bernstein basis at T = 0.05s. . . . . . . 8.11 IGDG solutions for different bivariate degrees (p, q) for N el = 8 × 8. . . . . . . . . . . . . . . . . . . 8.12 L 2 -errors for the 2D advection problem using the IGDG method in conjunction with RK 4. . . . 9.1 Control point lattice for quadratic rectilinear grid (on the right) and curvilinear grid (on the left). 9.2 Rectilinear grid on the right and curvilinear grid on the left (4 × 4 elements). . . . . . . . . . . . . 9.3 Plots and contour plots of the exact pressure p ex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.4 Plots and contour plots of u ex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.5 Plots and contour plots of v analytic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.6 Rectilinear patches. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.7 Plots and contour plots of numerical results for bivariate quadratic Bernstein basis with N el = 4 × 4 patches at T = 0.1s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.8 Plots and contour plots of numerical results for bivariate quadratic Bernstein basis N el = 8 × 8 patches at T = 0.1s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.9 IGDG solution u for different degrees p for N el = 4 × 4 elements at T = 0.1s. . . . . . . . . . . . . . 9.10 IGDG solution v for different degrees p for N el = 4 × 4 elements at T = 0.1s. . . . . . . . . . . . . . 9.11 IGDG solution p for different degrees p for N el = 4 × 4 elements at T = 0.1s. . . . . . . . . . . . . . 9.12 L 2 -error for the 2D acoustic problem using the IGDG method in conjunction with RK 4. . . . . . 9.13 Curvilinear patches. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .9.14 Plots and contour plots of numerical results for bivariate quadratic Bernstein basis with N el = 4 × 4 patches at T = 0.1s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.15 Plots and contour plots of numerical results for bivariate quadratic Bernstein basis N el = 8 × 8 patches at T = 0.1s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.16 IGDG solution u for different degrees p for N el = 4 × 4 elements at T = 0.1s. . . . . . . . . . . . . . 9.17 IGDG solution v for different degrees p for N el = 4 × 4 elements at T = 0.1s. . . . . . . . . . . . . . 9.18 IGDG solution p for different degrees p for N el = 4 × 4 elements at T = 0.1s. . . . . . . . . . . . . . 9.19 L 2 -error for the 2D acoustic problem using the IGDG method in conjunction with RK 4. . . . . . xix LIST OF SYMBOLS p Degree of the polynomial approximation. n Number of B-spline functions of degree p. n G Number of nodes for Gauss quadrature methods. ň Dimension of the time steps.

  The weights of nodes for Gauss quadrature methods.vTest function.

Figure 2 . 1 :

 21 Figure 2.1: Constant, linear, quadratic and cubic Bernstein polynomials.

Figure 2 . 2 :

 22 Figure 2.2: Linear, quadratic and cubic bivariate Bernstein polynomials.

Figure 2 . 3 :

 23 Figure 2.3: Bézier curves of various degrees and their control polygons.

Figure 2 .

 2 Figure 2.3 depicts three different Bézier curves, with their corresponding control polygons. Each control polygon is composed of its control points that are connected with line segments. Note that these control polygons are not necessary closed.

Figure 2 . 4 :

 24 Figure 2.4: Bézier curves with endpoint interpolation.

Definition 2 . 5 . 1 .

 251 (Convex hull)

Figure 2 . 5 :

 25 Figure 2.5: Cubic Bézier curve, its control polygon and the convex hull.

2 . 6 .

 26 The original curve is shown on the left. The new Bézier curve after repositioning of the control point P 2 is shown on the right. As we can see, the curve is globally modified. This property can be a drawback and will justify the construction of B-spline curves (see chapter 3).

Figure 2 . 6 :

 26 Figure 2.6: Quadric Bézier curve and repositioning of the control point P 2 .

2 . 7 .

 27 The original quadratic Bézier curve C 2 is shown in the beginning. The numbers of control points and basis functions increase simultaneously. The locations of the control points change, whereas the elevated curve is geometrically and parametrically identical to the first one.

Figure 2 . 7 :

 27 Figure 2.7: Degree elevation of a quadratic Bézier Curve.

Fig. 2 .

 2 Fig. 2.8 illustrates the De Casteljau evaluation of a point on a cubic Bézier curve at ζ = 2/3. This algorithmcan also be represented as a triangular scheme, starting with four control points P 0 , P 1 , P 2 and P 3 , and eventually reducing them to a single point P 3 0 .

Figure 2 . 8 :

 28 Figure 2.8: Geometric construction according to De Casteljau's algorithm for p = 3 and ζ = 2/3.

Fig. 2 .

 2 Fig.2.9 shows several conic curves, which all share the same Bézier points, only the inner weight w changes.

Figure 2 . 9 :

 29 Figure 2.9: Quadratic rational Bézier curves.

Fig. 2 .

 2 Fig. 2.10 gives an example of a tensor product Bézier surface of degree 3 × 3 and the corresponding control net.

Figure 2 . 10 :

 210 Figure 2.10: Tensor product Bézier patch of degree 3 × 3 and its control net.

. 2 )Figure 3 . 1 :

 231 Figure 3.1: Basis functions of degrees 1, 2 and 3 for uniform knot vector Ξ = {0, 1, 2, 3, ...}.

Figure 3 . 2 :

 32 Figure 3.2: Bivariates quadratic and cubic B-spline basis functions [16].

The equations ( 3 . 1 )

 31 and (3.2) show clearly that the choice of the knot vector has a significant influence on the B-spline functions. An example of quadratic basis functions for an open, uniform knot vector is presented in Fig. 3.3.

Figure 3 . 3 :

 33 Figure 3.3: Quadratic basis functions for the open-uniform knot vector Ξ = 0, 0, 0, 1, 2, 3, 3, 3 .

Figure 3 . 4 :

 34 Figure 3.4: Quadratic basis functions with reduced continuity at ξ = 1, Ξ = 0, 0, 0, 1, 1, 2, 3, 3, 3 .

Figure 3 . 5 :

 35 Figure 3.5: Quadratic B-spline curve.

Figure 3 . 6 :

 36 Figure 3.6: A quadratic B-spline curve and its control points. In the right, the curve after moving the control point P 4 .

Figure 3 . 7 :

 37 Figure 3.7: Before and after knot insertion (cubic B-spline curve).

Figure 3 .

 3 Figure 3.8: B-spline surface example.

Figure 3 . 9 :

 39 Figure 3.9: Bézier decomposition (bottom) from a quadratic B-spline basis (top) by knot insertion.

Figure 4 . 1 :

 41 Figure 4.1: The least-squares quadratic B-spline curve fitting.

Figures 4 .

 4 Figures 4.2 and 4.5 show two examples of least squares quadratic B-spline fitting surfaces. These two

  al , with n ev al = m ev al = 30.

Figure 4 . 2 :

 42 Figure 4.2: The least-squares quadratic B-spline Gaussian surface fitting.

Figure 4 . 3 :

 43 Figure 4.3: Bessel functions of the first kind-1D.

Figure 4 . 4 :

 44 Figure 4.4: Representation of the physical domain Ω.

Figure 4 . 5 :

 45 Figure 4.5: The least-squares quadratic B-spline surface fitting.

Figure 5 . 1 :

 51 Figure 5.1: Uniform P 1 mesh of [a, b].

Figure 5 . 2 :

 52 Figure 5.2: Global shape functions for the space V 1 h .

. 5 )

 5 Now we rewrite equation (5.1) with Ω =] -1, 1[ and f (u(x, t )) = cu(x, t ) such that c > 0 is the advection velocity. For simplicity in what follows, we assume c to be constant.

Figure 5 . 3 :

 53 Figure 5.3: The exact sine wave solution for the one-dimensional advection problem.

Figure 5 . 4 :

 54 Figure 5.4: Exact and standard Galerkin FEM P 1 solution for the one-dimensional linear advection problem at T = 0.4s.

  Petrov-Galerkin (PSPG) method: Motivated by mathematical analysis, another type of stabilization scheme has been established: the pressure-stabilizing Petrov-Galerkin. It was introduced by Hughes and his collaborators in 1986 for the Stokes problem and incompressible Navier-Stokes equations.

Figure 5 . 5 :

 55 Figure 5.5: Exact, Galerkin and SUPG solutions at T = 0.4s.

Figure 5 . 6 :

 56 Figure 5.6: SUPG FE P 1 solution for different values of α ∈ [0, 1].

Figure 5 . 7 :

 57 Figure 5.7: Convergence in number of d.o.f. for different choices of α.

  α

Figure 5 . 8 :

 58 Figure 5.8: Uniform P 2 mesh of [a, b].

  solutions of the problem (5.6) investigated by these methods, are depicted in Fig. 5.10 for different choices of the coefficient α. The numerical behavior of the method is evaluated on the spacial domain [-1, 1], which is subdivided into N -1 uniformly distributed subintervals with a number of d .o. f . = 2N -1 , and time domain [0, T ] where T > 0 is considered as the final time.

Figure 5 . 10 :

 510 Figure 5.10: SUPG quadratic Lagrange P 2 FEM in conjunction with RK 2 for the 1D advection problem.

Figure 5 .CONCLUSIOND

 5 Figure 5.11: L 2 -error for the advection problem with the linear and quadratic Lagrange FEM.

6. 2 .

 2

Figure 6 . 2 :

 62 Figure 6.2: (a) SUPG B-spline linear solution for the advection problem (α = 0.1). (b) SUPG FEM P 1 for the advection problem (α = 0.1) at T = 0.4s.

Figure 7 . 4 :

 74 Figure 7.4: L 2 -errors for the 1D advection problem with a sinusoidal initial condition, RK 2 and RK 4.

Figure 7 . 5 :

 75 Figure 7.5: Convergence rates in IGDG method as a function of the Bernstein function is degree p for the finest grid.

Figure 7 . 6 :

 76 Figure 7.6: Error in the L 2 -norm combining IGFEM, DGFEM, IGDG space discretization and explicit RK time integration.

Figure 7 .

 7 Figure 7.6 shows that DG methods are far more accurate than FEM. The comparison between the classical DG and IGDG methods is less clear, because both method gives almost the same accuracy with optimal convergence results with respect to the L 2 -norm.

8. 1 . 1

 11 Preliminaries -IGDG notationIn order to apply the IGA methodology, the physical domain Ω is subdivided into B-spline patches D e , S (Ω) := D e N pa e=1 , denote the interpatch boundary between the two patches D e and D l by Γ e,l and the collection of all interfaces by Γ, i.e., Γ e,l = D e D l and Γ e := l >e Γ e,l .

  represented by a set of B-spline patches. Each B-spline patch denoted by D e , is an image under a B-spline mapping of a parametric domain D e . D e = X e = (x e , y e ) ∈ R 2 | X e = T(ξ, η) such that (ξ, η) ∈ D e , where the transformation T is defined for all (ξ, η) ∈ D e by: T : D e -→ D e (ξ, η) -→x e (ξ, η), y e (ξ, η) .

Figure 8 . 1 :

 81 Figure 8.1: An example of a B-spline patch in physical space, parametric space, and the parent element used to perform numerical integration.

e=

  B p,q (x, y) e = B p,q T(ξ, η) e = B p,q (ξ, η) ...,p+1 and B

(

  e T ∇Φ e . Finally, one obtains:∇ Φ e = J e T ∇Φ e , ∇Φ e = (J e ) T -1 ∇ Φ e .

  p (ξ, η) e B j ,q (ξ, η) e | J e (ξ, η)| d D e = D B i ,p ( ξ) e B j ,q ( η) e | J e ( ξ, η)| | J e (ξ, η)| d D,

∂( 8 . 4 )ii

 84 t u e (x, y, t )v e (x, y)d D e + D e ∇ • -→ c u e (x, y, t ) v e (x, y)d D e = 0. (8.2)By applying Green's formula and introducing the numerical flux f , the weak formulation can be written as:D e ∂ t u e (x, y, t )v e (x, y)d D e = D e u e (x, y, t ) -→ c • ∇v e (x, y)d D e -Γ e v e (x, y) -→ c u e • -→ n e d Γ e . (8.3)We denote by -→ n e the outer unit normal to Γ e for the element D e . We define the restriction of the approximate solution u e h ≈ u e to D e via:u e h (x, y, t ) | D e = u e h (x, y, t ) = B e (x, y)u e (t ).By discretizing the problem on the Bézier basis associated with the element, the problem can be written: , j (x, y)B e k,l (x, y)d D e = , j (x, y) -→ c • ∇B e k,l (x, y)d D e -(x, y)( -→ c u e ) • -→ n e d Γ e ∀1 ≤ k ≤ p + 1 1 ≤ l ≤ q + 1.

  j (x, y)B e k,l (x, y)d D e = D e B e i , j (T(ξ, η))B e k,l (T(ξ, η))d D e , e i , j (x, y) .∇B e k,l (x, y)d D e (T(ξ, η))d D e .By performing the integration in the local parametric space D e , we get:M e kl ,i j = D e B i ,p (ξ) e B j ,q (η) e B k,p (ξ) e B l ,q (η) e | J e (ξ, η) | d D e ∀i , k = 1, ..., p + 1, j , l = 1, ..., q + 1. (ξ, η)| d D e ,where | J e | is the local Jacobien determinant which is described in section 8.1.3.By performing the integration in the reference square D, we get:M e kl ,i j = D B i ,p ( ξ) e B j ,q (η) e B k,p ( ξ) e B l ,q (η) e | Ĵ e (ξ,η) || J e (ξ, η) | d D ∀i , k = 1, ..., p + 1 j , l = 1, ..., q + 1, e ( ξ, η)| | J e (ξ, η)| d D.

  It is also necessary to define for the patch D e , the normal vectors on the four interfaces Γ e 1 ,

. 8 )

 8 where | d e | is a typical length scale for the Bézier element D e . Obviously, 0 < C C F L ≤ 1 is the Courant-Friedrichs-Lewy (CFL) coefficient and ∥ -→ c ∥= c 2 x + c 2 y .

Figure 8 . 3 : 2 N 1 el and h y = 2 N 2 el.

 832122 Figure 8.3: Analytical solution for the bi-dimensional advection problem, for N el = 4 × 4 at T = 0.5s.

Figure 8 .

 8 Figure 8.5 shows the numerical solutions from IGDG space discretization and explicit RK 4 time integration for the linear, quadratic, cubic and quartic Bernstein cases. We can see the effect of the degree elevation of the bivariate Bernstein basis function on the accuracy.

Figure 8 . 5 :

 85 Figure 8.5: IGDG solution for N el = 8 × 8 patches for different degrees (p, q).

Figure 8 . 7 :

 87 Figure 8.7: Contour plots of numerical results for bivariate quadratic Bernstein basis at T = 0.05s.

4 :

 4 L 2 -error for the 2D advection problem and convergence order for the IGDG method for the cubic (left) and quartic (right) Bernstein bases in conjunction with RK 4 time discretisation.

Figure 8 . 9 :

 89 Figure 8.9: L 2 -error for the 2D advection problem using the IGDG method in conjunction with RK 4.

Figure 8 .

 8 Figure 8.11 shows the numerical solutions from IGDG space discretization and explicit RK 4 time integration for the linear, quadratic, cubic and quartic Bernstein cases. We can see the effect of the degree elevation of the bivariate Bernstein basis function on the accuracy.

Figure 8 . 11 :

 811 Figure 8.11: IGDG solutions for different bivariate degrees (p, q) for N el = 8 × 8.

Figure 8 . 12 :

 812 Figure 8.12: L 2 -errors for the 2D advection problem using the IGDG method in conjunction with RK 4.

  D e∂ t u e (x, y, t ) + ∂ x p e (x, y, t ) w e (x, y)d D e = 0.(9.4)By applying Green's formula, the weak formulation can be written as:∂ ∂t D e u e (x, y, t )w e (x, y)d D e + w e (x, y)d Γ ee (x, y)d D e = 0. the outer unit normal to Γ e of the element D e , so we get: ∂ ∂t D e u e (x, y, t )w e (x, y)d D e = D e p e (x, y, t ) -→ e x • ∇w e (x, y)d D e -Γ e p e (x, y, t ) -→ e x • -→ n e w e (x, y)d Γ e .

∂ 1 

 1 t v e (x, y, t ) + ∂ y p e (x, y, t ) w e (x, y)d D e = 0. (9.6) By applying Green's formula, the weak formulation can be written as:∂ ∂t D e v e (x, y, t )w e (x, y)d D e + w e (x, y)d Γ ee (x, y)d D e = 0. , so we get: ∂ ∂t D e v e (x, y, t )w e (x, y)d D e = D e p e (x, y, t ) -→ e y • ∇w e (x, y)d D e -Γ e p e (x, y, t ) -→ e y • -→ n e w e (x, y)d Γ e .

  R e x u e + R e y p e -F e x + F e y ∀t ∈ [0, T ] ∀1 ≤ e ≤ N el .

i

  , j (x, y)B e k,l (x, y)d D e = D e B e i , j (T(ξ, η))B e k,l (T(ξ, η))d D e .By performing the integration in the local parametric space D e , we get:M e kl ,i j = D e B i ,p (ξ) e B j ,q (η) e B k,p (ξ) e B l ,q (η) e | J e (ξ, η) | d D e , ∀i , k = 1, ..., p + 1, j , l = 1, ..., q + 1,where | J e | is the local Jacobien determinant defined in the previous chapter by:

  B e i j (x, y) • ∇B e kl (x, y)d D e T(ξ, η) d D e , x, y) • ∇B e kl (x, y)d D e T(ξ, η) d D e , R e kl ,i j

=1e|+

  J e (ξ, 1) || J e ( ξ, 1)d ξ ,q (η) e | J e (0, η) || J e (0, η)d η.It is also necessary to define for the patch D e , the normal vectors on the four interfaces Γ e 1 , Γ e 2 ,

Figure 9 . 2 :

 92 Figure 9.2: Rectilinear grid on the right and curvilinear grid on the left (4 × 4 elements).

Figure 9 . 3 :

 93 Figure 9.3: Plots and contour plots of the exact pressure p ex .

Figure 9 . 4 :

 94 Figure 9.4: Plots and contour plots of u ex .

Figure 9 . 5 :

 95 Figure 9.5: Plots and contour plots of v analytic .

  Figure 9.7 depicts the IGDG numerical solutions for the bivariate quadratic Bernstein case for the initial mesh with N el = 4×4 patches (as shown on the left in Fig. 9.6). We subsequently add a refinement (on the right of Fig. 9.6) of ( the IGDG numerical solutions are shown in Fig. 9.8.

Figure 9 . 6 :

 96 Figure 9.6: Rectilinear patches.

Figure 9 . 7 :

 97 Figure 9.7: Plots and contour plots of numerical results for bivariate quadratic Bernstein basis with N el = 4 × 4 patches at T = 0.1s.

Figure 9 . 8 :

 98 Figure 9.8: Plots and contour plots of numerical results for bivariate quadratic Bernstein basis N el = 8 × 8 patches at T = 0.1s.

Figure 9 . 9 ,

 99 Figure 9.9, Fig. 9.10 and Fig. 9.11 shows the numerical solutions from IGDG space discretization and explicit RK 4 time integration for the quadratic, cubic and quatric Bernstein cases respectively. We can see the effect of the degree elevation of the bivariate Bernstein basis function on the accuracy.

Figure 9 . 9 :

 99 Figure 9.9: IGDG solution u for different degrees p for N el = 4 × 4 elements at T = 0.1s.

2 )

 2 as shown on the right in Fig.9.13.

Figure 9 . 13 :

 913 Figure 9.13: Curvilinear patches.

Figure 9 . 14 :

 914 Figure 9.14: Plots and contour plots of numerical results for bivariate quadratic Bernstein basis with N el = 4 × 4 patches at T = 0.1s.

Figure 9 . 15 :

 915 Figure 9.15: Plots and contour plots of numerical results for bivariate quadratic Bernstein basis N el = 8 × 8 patches at T = 0.1s.

Figure 9 .

 9 Figure 9.16, Fig. 9.17 and Fig. 9.18 shows the numerical solutions from IGDG space discretization and explicit RK 4 time integration for the quadratic, cubic and quatric Bernstein cases respectively. We can see the effect of the degree elevation of the bivariate Bernstein basis function on the accuracy.

Figure 9 .

 9 Figure 9.16: IGDG solution u for different degrees p for N el = 4 × 4 elements at T = 0.1s.

Figure 9 . 19 :

 919 Figure 9.19: L 2 -error for the 2D acoustic problem using the IGDG method in conjunction with RK 4.

  perbolique unidimensionnel. Le présent chapitre et le chapitre précédent sont consacrés à l'extension au cas des problèmes hyperboliques bidimensionnels. Nous avons décrit le développement d'une méthode de GD dans le cadre AIG pour résoudre l'équation d'advection en 2D dans le chapitre précédent et le problème acoustiques de premier ordre en 2D dans le présent chapitre, qui est exprimé sous la forme d'un système EDP linéaire hyperbolique. Comme cela a été souligné précédemment, le domaine de calcul est divisé en sous-domaines, composés de patches B-splines. L'approche GD a été appliquée au niveau des éléments, chaque élément étant un patch de Bézier construit à partir de patches B-splines initiaux. La solution du problème est approximée dans chaque élément sans aucune exigence de continuité pour la solution discrète sur les interfaces. Par conséquent, la méthode de GD dans le cadre AIG peut être facilement formulée et mise en oeuvre, elle est capable de gérer facilement des géométries complexes. La flexibilité de la méthode pour gérer différentes géométries et pour travailler avec différents éléments a été montrée. Ainsi, le comportement numérique de la méthode est évalué et il a été montré un taux de convergence optimal selon la norme L 2 .

A. 1

 1 nis qui généralise l'ensemble des fonctions de base pôlynomiales aux B-splines ou plus généralement aux NURBS. Comme le montre la littérature, ce choix garantit plusieurs avantages de la paramétrisation exacte des géométries définies par CAO à une plus grande précision par degré de liberté. Il permet également des champs de solution réguliers. Ces raisons ont fait de l'AIG un sujet émergeant ces dernières années.Dans cette thèse, nous avons appliqué la méthode de l'AIG à certains problèmes hyperboliques. Nous avons considéré les méthodes standards de Galerkin ainsi que les méthodes stabilisées, en mettant particulièrement l'accent sur les B-splines. L'AIG possède un ensemble de caractéristiques attractives du point de vue de la précision et de la commodité de mise en oeuvre qui ne sont pas présentes dans les discrétisations d'EF standards. En fait, l'AIG permet une représentation exacte d'une large classe de géométries, même sur des maillages très grossiers. En particulier, les domaines géométriques ayant des sections coniques comme les cercles, les cylindres, les sphères, les ellipsoides, etc. peuvent être représentés exactement à l'aide de NURBS. Les raffinements peuvent être effectués par subdivision de la grille (en insérant des noeuds) ou par élévation de l'ordre polynomial des fonctions de base de la même manière que la méthode des EF traditionnels.Vu que l'objectif principal de l'AIG est de généraliser et d'améliorer l'AEF classique, nous avons commencé cette thèse en revisitant l'analyse originale, à savoir l'AEF, dans le contexte d'EDP hyperboliques avec la stabilisation SUPG.Des exemples numériques de problèmes d'advection ont été donnés pour les bases classiques deLagrange et B-splines. Il a été montré que la méthode d'AIG peut être bien adaptée aux problèmes hyperboliques si la stabilisation est correctement calibrée. Les taux de convergence ont été mesurés dans la norme standard L 2 . La méthode d'AIG a été testée sur différents ordres polynomiaux. On a constaté que, pour un ordre polynomial donné, une régularité élevée donne une erreur plus faible par degré de liberté, par rapport à une faible régularité, sans exception. Par ailleurs, un inconvénient possible de la méthode SUPG est la sensibilité de la solution au paramètre de stabilisation, dont la valeur optimale n'est pas déterminée précisément par la théorie disponible et qui est fastidieuse à sélectionner dans la pratique. Par conséquent, nous avons proposé dans cette thèse une nouvelle méthode qui combine l'AIG avec la méthode de GD, appelée méthode IGDG, pour résoudre des problèmes hyperboliques. La principale raison de l'utilisation des méthodes de GD dans cette thèse réside dans leur capacité à fournir des méthodes numériques stables pour les problèmes hyperboliques, pour lesquels la MEF classique est bien connue pour ses performances médiocres. Notre méthode tire parti des méthodes de l'AIG et de GD. En fait, la méthodologie de GD est adoptée au niveau des patches de Bézier, c'est-à-dire que nous utilisons l'AIG traditionnelle dans chaque patch de Bézier et la méthode de GD sur les interfaces de patch pour assembler les multiples patches. Les patches de Bézier, considérés comme des éléments, sont construits par transformation du domaine B-spline initial (l'extraction de Bézier est une technique classique de CAO). Nous considérons alors des cas scalaires et des cas systèmes de lois de conservation, comme cas tests. Ces exemples de tests montrent que la méthode IGDG obtenue en conjonction avec la méthode de RK pour la discrétisation temporelle est stable, de haute précision, et peut facilement gérer les géométries courbes et les conditions aux limites. Par conséquent, cette méthode peut être facilement formulée et mise en oeuvre. Ainsi, le comportement numérique de la méthode a été évalué et on a montré un taux de convergence optimal pour l'advection 2D et les problèmes acoustiques. La plus grande précision obtenue en utilisant des éléments curvilignes a été démontrée, en particulier lorsque des grilles grossières sont utilisées. Une précision de l'ordre de cinq a été atteinte. La méthode actuelle s'est avérée très efficace pour les problèmes hyperboliques linéaires et les systèmes à une et deux dimensions. Maintenant, la méthode IGDG proposée devrait être étendue à des lois de conservation non linéaires plus complexes, comme les équations d'Euler compressible, par exemple. Des méthodes de GD efficaces pour ces systèmes non linéaires ont été construites récemment pour les grilles rectilignes classiques. L'extension à des bases de Bézier ne devrait pas présenter de difficultés spécifiques, sauf pour la capture des discontinuités de solution pouvant apparaître pour les lois de conservation non linéaires. Dans cette perspective, on pourrait étendre les limiteurs généralisés proposés par Cockburn et ses coauteurs (voir annexe) aux représentations de Bézier locales. Une autre extension nécessaire concernerait la construction du domaine de calcul initial B-spline. Les géométries considérées dans ce travail sont plutôt simples et représentées par un seul patch B-spline, avant l'extraction de Bézier. Pour prendre en compte des géométries plus complexes, il convient de proposer un cadre général permettant de gérer facilement un ensemble de patches B-spline. Dans cette perspective, nous soulignons que l'approche proposée, basée sur les concepts GD, est plus flexible que la méthode d'AIG d'origine car aucune contrainte de régularité n'est nécessaire à l'interface entre les patches. Enfin, il convient d'exploiter la grande flexibilité offerte par les représentations B-splines en termes de raffinement p-et h-afin de proposer des stratégies de raffinement locales automatisées. A P P E N D I X A Gaussian quadrature Let Ω = [a, b] ⊂ R. A quadrature rule is defined by a set of nodes (X G (k)) 1≤k≤n G ∈ Ω and a set of weights (ω G (k)) 1≤k≤n G ∈ R. The fundamental result of Gaussian quadrature states that the optimal abscissas of the n G -point Gaussian quadrature formulas are precisely the roots of the orthogonal polynomial for the same interval and weighting function. Gaussian quadrature is optimal because it integrates all polynomials up to degree 2n G -1 exactly.

1 :ω 2 A

 12 Gauss-Legendre nodes and coefficientsWe can transform any given integral on the interval [a, b] into an integral on the interval [-1, 1], simply use:G (k) f (ba)X G (k) + (b + a)

  Bézier was an engineer in the Renault car company and developed this method in 1966. Actually, another French engineer, Paul de Casteljau at

Citroën developed the same technology some years earlier. A further development to Bézier's method were B-splines which provide more flexibility in the modeling of free-form curves and surfaces. Since 1975 nonuniform rational B-splines (NURBS) have been used in CAD programs, as a generalization of B-splines. The development of NURBS provided a technology that can exactly describe circular shapes (cylinders, spheres, etc.

  La méthode des éléments finis a ses origines au début des années 1960, elle est aujourd'hui la méthode prédominante dans l'analyse des problèmes elliptiques ou paraboliques en raison de sa flexibilité à représenter des domaines géométriques complexes et sa base théorique solide. L'idée consiste à décomposer le domaine en plusieurs éléments «finis», définis par un ensemble de points nodaux et des fonctions de base interpolantes.Bien que la MEF a été largement utilisée pour simuler de nombreux phénomènes physiques en raison de sa flexibilité à représenter des domaines géométriques complexes, il est bien connu dans la littérature de la MEF que des difficultés numériques se posent lors de la résolution des EDPs hyperboliques.

	En effet, lors de l'utilisation de la méthode d'EF standard appliquée à de tels EDPs, des oscillations para-
	sites (non physiques) indésirables (phénomène de Gibbs) sont fréquemment détectées dans les solutions
	numériques. Un remède à cet inconvénient, répandu dans la littérature, consiste à ajouter une viscosité "ar-
	tificielle" à un schéma numérique standard (instable). D'une part, cette viscosité artificielle doit amortir les
	oscillations mais, d'autre part, elle ne doit pas entacher la solution numérique. À la fin des années 1970 et au
	début des années 1980, un grand nombre de méthodes dites stabilisées ont été développées avec des idées
	différentes [14] [32] [48]. Ceci est réalisé grâce à l'utilisation d'une formulation Petrov-Galerkin [35], où les
	fonctions test sont modifiées de telle sorte qu'elles pondèrent le noeud en amont plus que le noeud en aval
	[17] [38]. Parmi eux, la plus populaire, dite Streamline-Upwind Petrov-Galerkin (SUPG), a été introduite
	par Brooks et Hughes. Elle a été proposée d'abord dans le contexte des équations d'advection-diffusion
	et les équations de Navier-Stokes incompressible

èles mathématiques permettant d'exprimer la conservation d'une quantité physique, comme par exemple la masse, l'énergie, etc. Ils découlent naturellement des lois de conservation en physique. En particulier, ils décrivent une grande variété de phénomènes impliquant le mouvement des ondes (acoustiques, élastiques, électromagnétiques) ou le transport advectif de substances.

Le large éventail d'applications des EDPs hyperboliques a conduit à une activité de recherche très intense dans ce domaine. Cela a permis de développer un ensemble de méthodes numériques pour des approximations précises et efficaces du point de vue du calcul des solutions à ces problèmes. Il existe trois grandes familles de méthodes qui sont largement utilisées: la méthode des différences finies (MDF), la méthode des volumes finis (MVF) et la méthode des éléments finis (MEF). Ces méthodes se sont révélées extrêmement utiles pour modéliser un large éventail de phénomènes. Pour garder cette thèse autonome, nous présentons brièvement chacune de ces trois méthodes dans le contexte des EDPs hyperboliques.

Historiquement, la méthode des différences finies (MDF) était la première méthode utilisée pour produire des approximations des solutions des EDPs hyperboliques. Elle a été introduite par Euler au 18ème siècle et représente la méthode la plus simple pour résoudre des problèmes sur des géométries simples

[START_REF] Morton | Numerical solution of partial differential equations: an introduction[END_REF]

.

L'idée principale de cette méthode est de remplacer les dérivées partielles de l'inconnu par leurs approximations par des différences finis. La méthode DF est remarquable par la grande variété de schémas qui peuvent être utilisés pour approcher une EDP donnée, par exemples, schémas explicites (Forward Euler, Upwind, Lax-Friedrichs, Lax-Wendroff, Leapfrog, ...) et schémas implicites (Backward Euler, Crank-Nicolson, ...)

[START_REF] Trefethen | Finite difference and spectral methods for ordinary and partial differential equations[END_REF]

. Bien que cette méthode puisse être facilement formulée et mise en oeuvre, son application à des problèmes avec des géométries réalistes est assez lourde, ce qui rend la méthode peu attrayante pour les problèmes industriels. Cela a nécessité d'autres méthodes plus flexibles, telles que les méthodes de FV et d'EF.

La méthode des VF peut gérer des géométries complexes ce qui la rend plus attrayante que la méthode de DF pour les problèmes complexes. Elle est basée sur la forme conservative au lieu de la forme différentielle pour estimer les valeurs des champs inconnus. On divise le domaine en cellules et on approxime l'intégrale totale sur chaque cellule de la grille. Le liens entre les quantités de cellules dans la MVF dépendent du flux entre les volumes de contrôle voisins. Cela signifie que la MVF représente le flux d'informations dans la structure du maillage de manière conservative

[START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]

. Elle permet l'utilisation de grilles non struc-CHAPTER 1. INTRODUCTION turées pour gérer des géométries complexes.

  Au cours des dernières années, la méthode de GD dans le cadre IG a manifesté un intérêt croissant pour

	L'AIG étant une extension de la MEF, nous commençons par revoir le cadre de l'analyse originale au cin-
	quième chapitre, à savoir: la MEF. Le besoin de stabilisation est souligné et les idées de stabilisation basées
	sur le concept Petrov-Galerkin sont discutés. Nous nous concentrons sur la méthode de stabilisation SUPG
	et une attention particulière est accordée à l'étude du paramètre de stabilisation τ.
	Au niveau du sixième chapitre, les différentes procédures de calcul pour l'AIG sont passées en revue
	dans le contexte de la MEF, en revisitant le problème d'advection unidimensionnel qui est présenté dans le
	chapitre précédent. Dans ce chapitre nous utilisons les B-splines (en raison de la simplicité du domaine)
	comme fonctions de base, il n'est pas difficile de le généraliser pour d'autres splines telles que les NURBS.
	Des comparaisons détaillées entre l'AIG et la MEF classiques sont discutées. Dans ce contexte de l'AIG, nous
	considérons alors l'application des méthodes GD. En effet, l'argument majeur pour utiliser les méthodes
	de GD réside dans leur capacité à fournir des méthodes numériques stables pour les EDPs hyperboliques,
	pour lesquelles la MEF classique est bien connue pour ses performances médiocres. Par conséquent, au
	septième chapitre, nous traitons le problème d'advection unidimensionnelle en combinant la méthode de
	la solution numérique des EDPs elliptiques [40] [54] [69][95]. Les avantages des espaces d'approximation l'AIG avec la méthode de GD. Nous notons que la méthodologie de GD est adoptée au niveau du patch,
	locaux sans exigences de continuité offerts par les méthodes de GD [5] [28] sont alors utilisés pour gérer les c'est-à-dire que nous employons l'AIG classique dans chaque patch et utilisons la méthode de GD à travers
	calculs multi-patch. les interfaces de patch. De plus, une transformation de la base B-spline est nécessaire pour introduire des
	L'objectif principal de cette thèse est d'étudier l'utilisation des AIG pour résoudre certains problèmes discontinuités aux interfaces, sans modifier la géométrie du domaine. Les caractéristiques avantageuses
	hyperboliques. En particulier, nous décrivons la méthode de Galerkin continue et discontinue en utilisant des deux méthodes AIG et GD nous permettent de concevoir une formulation prometteuse.
	la base des B-splines. Un accent particulier est mis sur la méthode de GD, car elle est considérée comme Avec quelques ajustements, les huitième et neuvième chapitres sont consacrés à l'étude de deux ex-
	l'une des méthodes les plus efficaces et à la croissance la plus rapide, avec des applications dans divers emples numériques en 2D, le problème d'advection est d'abord présenté, suivi par les équations d'ondes
	problèmes, pas nécessairement hyperboliques. Les discontinuités des fonctions de base, qui offrent une acoustiques, où les deux systèmes sont résolus sur plusieurs domaines (cartésien, linéaire et curviligne).
	plus grande souplesse d'analyse, rendent cependant la méthode délicate pour gérer des géométries réal-Enfin, au niveau du dernier chapitre, nous terminons avec quelques remarques et perspectives finales.
	istes à partir de la CAO. Les résultats des différentes études réalisées sont résumées et discutées et des idées de recherches futures
	La thèse est divisée en quatre parties principales: la première met l'accent sur les fonctions de base de sont proposées.
	Bernstein et B-splines utilisées en CAO. Elle est consacré à donner leurs définitions et propriétés de base.
	Nous présenterons dans la deuxième et troisième parties l'extension de l'analyse classique à l'AIG pour
	les méthodes d'EF et GD, pour le problème d'advection unidimensionnel. Dans la dernière partie, nous
	traitons des problèmes hyperboliques en deux dimensions en combinant l'AIG avec la méthode de GD. Il
	convient de mentionner que dans tout ce travail, la discrétisation des équations dans le temps se fait au
	moyen de méthodes de Runge-Kutta explicites d'ordre élevé [18] [19] [34] [80].
	Plus précisément, les deuxième et troisième chapitres fournissent une introduction complète aux prin-
	cipales idées et propriétés des fonctions Bernstein, B-splines et NURBS, qui constituent la base de l'AIG.
	Dans le même contexte, la construction des courbes B-spline et des surfaces au sens des moindres carrés
	est présentée au quatrième chapitre. L'analyse est accompagnée d'exemples de problèmes univariés et bi-
	variés.

Definition 2.1.2. (Multivariate Bernstein) In order to define Bernstein basis in higher dimensions, we make use of the tensor product construction. Let

	p = (p 1 , p 2 , ..., p d ) be a vector in N d . The d -dimensional Bernstein polynomials are defined by a tensor product
	of d univariate Bernstein polynomials with possibly different degrees p = (p 1 , p 2 , ..., p d ) and multi-indices
	k = (k 1 , k 2 , ..., k d ).
	Therefore, ∀ζ = (ζ 1 , ζ 2 , ..., ζ d ) ∈ [0, 1] d we get:

  1970 par Cox et DeBoor, pour remédier à l'inconvénient de la globalité des courbes de Bézier, à travers un algorithme efficace car hiérarchique (analogue à celui de Casteljau), stable numérique-

	C H A P T E R	4
	CURVE AND SURFACE FITTING

ment (coefficients multiplicatifs toujours positifs) et interprétable géométriquement: le déplacement d'un point de contrôle de la courbe n'affecte ainsi plus qu'une partie limitée de la courbe, ce qui amène un plus grand confort dans la Conception Assistée par Ordinateur (CAO). Par conséquent, les courbes B-splines permettent d'approcher des points de manière lisse, comme les coubes de Bézier. Leur avantage est qu'elles sont plus lisses et plus faciles à contrôler. Dans ce chapitre, nous avons présenté des définitions formelles des fonctions et des courbes B-splines, ainsi qu'une étude détaillée des propriétés les plus fondamentales des B-splines, illustrées par des exemples et des figures.

  ξ 2 , ..., ξ n+p+1 . Let (ζ k ) 1≤k≤n ev al such that ζ 1 < ζ 2 < ... < ζ n ev al and let Q be a set of ordered and distinct points Q k , ∀k = 1, .., n ev al to be fitted.

Table 5

 5

	.1 shows the error of the SUPG-P 1 FEM numerical approximations in L 2 -norm for various number of
	degrees of freedom (d.o.f.) (notice that the number of d.o.f. is equal to N ) and coefficient of the stabilization

parameter α ∈ [0, 1]. The L 2 -errors of the numerical approximations are depicted in Fig.

5

.7 for different values of α in log scale. ∥ e(T ) ∥ L 2 α d .o. f . = 16 d .o. f . = 32 d .o. f . = 64 d .o. f . = 128 d .o. f . = 256 d .o. f . = 512 d .o. f . =

Table 5 .

 5

1: L 2 -errors of the SUPG FE P 1 method for the one-dimensional advection problem.

Table 5 . 2

 52

	0.6	-	0.36	0.59	0.76	0.87	0.93	0.97
	0.5	-	0.3	0.65	0.8	0.91	0.93	0.97
	0.4	-	0.51	0.7	0.84	0.91	0.95	0.98
	0.3	-	0.61	0.77	0.87	0.93	0.97	0.98
	0.2	-	0.74	0.84	0.91	0.94	0.97	1
	0.1	-	1.03	0.92	0.94	0.97	1	1
	0.09	-	1.12	0.91	0.95	0.77	1.2	1
	0.08	-	1.2	0.92	0.95	0.98	1	1.02
	0.07	-	1.3	0.94	0.96	0.98	1	1
	0.06	-	1.5	0.96	0.96	0.98	1	1
	0.05	-	1.6	1	0.98	0.98	1	1
	0.04	-	1.7	1.07	0.98	0.99	0.98	1
	0.03	-	2	1.22	1	1	0.99	1

d .o. f . = 16 d .o. f . = 32 d .o. f . = 64 d .o. f . = 128 d .o. f . = 256 d .o. f . = 512 d .o. f . = 1024 : Convergence rates.

Table 5 .

 5 These results have been obtain without stabilization. e(T ) d .o. f . = 16 d .o. f . = 32 d .o. f . = 64 d .o. f . = 128 d .o. f . = 256 d .o. f . = 512 d .o. f . = 1024 1.420E -02 3.406E -03 8.476E -04 2.1154E -04 5.285E -05 1.321E -05 3.302E -06 4: The L 2 -error norm for RK 4 time discretization.

Table 7 .

 7 Therefore, we employ RK 4 method to recover optimal convergence rate for p ≥ 2. 2: L 2 -error for the IGDG method in conjunction with RK 2 time discretisation for various element sizes and degree of Bézier basis p = 0, 1, 2.

	p	h	h 2	h 4	h 8	h 16	rate
	0	4.649E -01 3.054E -01 1.699E -01 8.841E -02 4.477E -02	1
	1	5.137E -02 1.322E -02 3.255E -03 8.176E -04 2.024E -04	2
	2	4.255E -03 7.179E -04 1.594E -04 3.581E -05 8.652E -06 2.5
	p	h	h 2	h 4	h 8	h 16	rate
	2	3.692E -03 4.671E -04 5.758E -05 7.083E -06 8.709E -07	3
	3	9.299E -05 5.724E -06 3.322E -07 1.849E -08 1.062E -09	4
	4	8.069E -06 2.497E -07 7.767E -09 2.359E -10 7.241E -12	5

Table 7 .

 7 3: L 2 -error for the IGDG method in conjunction with RK 4 time discretisation for various element sizes and degree of Bézier basis p = 2, 3, 4.

Table 8 .

 8 .11) Tab. 8.1 and Tab. 8.2 summarize the convergence results in the L 2 -norm for the bivariate linear,quadratic, cubic and quartic Bernstein basis functions, from which we can see that the L 2 -convergence rate is approx-1: L 2 -error for the 2D advection problem and convergence order for the IGDG method for the linear (left) and quadratic (right) Bernstein bases in conjunction with RK 4 time discretisation.

	imately p + 1.
	Mesh	L 2 -error	rate	Mesh	L 2 -error	rate
	h	2.235E -01	-	h	1.136E -01	-
	h 2 h 4 h 8	1.012E -01 1.14 2.586E -02 1.96 5.807E -03 2.15	h 2 h 4 h 8	1.476E -02 2.94 1.713E -03 3.10 2.047E -04 3.06
	Mesh	L 2 -error	rate	Mesh	L 2 -error	rate
	h	3.170E -02	-	h	1.706E -02	-
	h 2 h 4 h 8	5.090E -03 2.63 2.6488E -04 4.26 1.650E -05 4.00	h 2 h 4 h 8	3.639E -04 5.55 9.763E -06 5.22 2.851E -07 5.09

Table 8 .

 8 2: L 2 -error for the 2D advection problem and convergence order for the IGDG method for the cubic (left) and quartic (right) Bernstein bases in conjunction with RK 4 time discretisation.

Table 8 .

 8 for the cubic and quartic Bernstein. 3: L 2 -error for the 2D advection problem and convergence order for the IGDG method for the linear (left) and quadratic (right) Bernstein bases in conjunction with RK 4 time discretisation.

	Mesh	L 2 -error	rate	Mesh	L 2 -error	rate
	h	2.248E -01	-	h	1.222E -01	-
	h 2 h 4 h 8	1.043E -01 1.10 2.931E -02 1.83 6.818E -03 2.10	h 2 h 4 h 8	2.220E -02 2.46 3.410E -03 2.70 3.874E -04 3.13
	Mesh	L 2 -error	rate	Mesh	L 2 -error	rate
	h	4.397E -02	-	h	1.845E -02	-
	h 2 h 4 8 h	5.362E -03 3.03 3.473E -04 3.94 1.915E -05 4.18	h 2 h 4	8.373E -04 4.46 2.871E -05 4.86

h 8 7.553E -07 5.19

Table 8 .

 8

  ∂ t p e (x, y, t ) + ∂ x u e (x, y, t ) + ∂ y v e (x, y, t ) w e (x, y)d D e = 0.

						R e y p e -F e y	∀t ∈ [0, T ], ∀1 ≤ e ≤ N el .	(9.7)
	9.2.4 Third variational equation
		In order to write the last variational equation of the system (9.4), we multiply the equation (3) with the
	test function w e and integrate over D e :
												(9.8)
			D e
	By applying Green's formula, the weak formulation can be written as:
	∂ ∂t D e	p Γ e	 	u e (x, y, t ) v e (x, y, t )	  • -→ n e w e (x, y)d D e -	D e	 	u e (x, y, t ) v e (x, y, t )	  • ∇w e (x, y)d D e = 0.
	so we get:
		∂ ∂t D e	p k,l as test function,
	the problem is written:
	p+1 i =1	q+1 j =1	∂ t p e i , j (t )	D e	B e i , j (x, y)B e k,l (x, y)d D e	=	p+1 i =1	q+1 j =1	u e i , j (t ) -→ e x + v e i , j (t ) -→ e y	D e	B e i , j (x, y) • ∇B e k,l (x, y)d D e
											-	Γ e	B e k,l (x, y) -→ e x u e + -→ e y v e • -→ n e d Γ e ,
										∀1 ≤ k

e (x, y, t )w e (x, y)d D e + e (x, y, t )w e (x, y)d D e = D e u e (x, y, t ) -→ e x + v e (x, y, t ) -→ e y • ∇w e (x, y)d D e -Γ e u e (x, y, t ) -→ e x + v e (x, y, t ) -→ e y • -→ n e w e (x, y)d Γ e .

By discretizing the problem on the Bézier basis associated with the element and using B e

Table 9 .

 9 The L 2 -error is now presented for rectilinear patches. Tables 9.1 and 9.2 summarize the convergence results in the L 2 -norm for the bivariate quadratic, cubic and quartic Bernstein basis functions. 1: L 2 -error for the 2D acoustic problem and convergence order for the IGDG method for the quadratic Bernstein bases in conjunction with RK 4 time discretisation.

			Mesh	L 2 -error	rate
			h	1.951E -01	-
			h 2 h 4 h 8 h 16	5.241E -02 1.89 7.224E -03 2.85 1.928E -03 1.90 2.524E -04 2.93
	Mesh	L 2 -error	rate		Mesh	L 2 -error	rate
	h	1.324E -01	-		h	9.624E -02	-
	h 2 h 4 h 8 h 16	2.622E -02 2.33 6.244E -03 2.07 1.889E -03 1.72 2.480E -04 2.92		h 2 h 4 h 8 h 16	2.579E -02 1.78 6.220E -03 2.05 1.890E -03 1.71 2.474E -04 2.93

Table 9 .

 9 2: L 2 -error for the 2D acoustic problem and convergence order for the IGDG method for the cubic (left) and quartic (right) Bernstein bases in conjunction with RK 4 time discretisation.The corresponding convergence data for uniform refinement h, h

2 , h 4 , h 8 and h 16 are shown in Figure

Table 9 .

 9 Tables 9.3 and 9.4 summarize the convergence results in the L 2 -norm for the bivariate quadratic, cubic and quartic Bernstein basis functions. 3: L 2 -error for the 2D acoustic problem and convergence order for the IGDG method for the quadratic Bernstein bases in conjunction with RK 4 time discretisation.

			Mesh	L 2 -error	rate
			h	1.983E -01	-
			h 2 h 4 h 8 h 16	4.541E -02 2.13 3.436E -03 3.72 4.066E -04 3.07 5.084E -05 2.99
	Mesh	L 2 -error	rate		Mesh	L 2 -error	rate
	h	9.907E -02	-		h	1.527E -02	-
	h 2 h 4 h 8 h 16	2.740E -03 5.17 3.173E -04 3.11 2.019E -05 3.97 1.310E -06 3.94		h 2 h 4 h 8 h 16	1.046E -03 3.86 2.875E -05 5.18 9.704E -07 4.88 3.056E -08 4.98

Table 9 .

 9 4: L 2 -error for the 2D acoustic problem and convergence order for the IGDG method for the cubic (left) and quartic (right) Bernstein bases in conjunction with RK 4 time discretisation.The corresponding convergence data for uniform and refinement h, h

	2 ,	h 4 ,

Dans ce chapitre, nous avons examiné les polynômes de Bernstein et les courbes de Bézier détaillées dans les deux cas, rationnel et non rationnel, et les surfaces par produits tensoriels ont été introduites. Diverses exemples de courbes de Bézier ont été générées.Nous avons discuté des propriétés des courbes de Bézier, ces propriétés sont directement impliquées dans l'algorithme de Casteljau ou dans les propriétés des polynômes de Bernstein. Nous avons mentionné les propriétés les plus fondamentales et nous avons fourni les preuves de ces propriétés brièvement.En conclusion, nous avons vu que les courbes de Bézier ont un certain nombre de caractéristiques qui les rendent intéressantes pour l'exploration du paramétrage géométrique. Ces courbes sont faciles à décrire paramétriquement. Elles sont faciles à représenter graphiquement elles ont aussi l'avantage d'être incroyablement pratiques, comme en témoigne leur utilisation dans la conception géométrique assistée par ordinateur et infographie depuis les année 1960.

 Figure 5.10 shows that the SUPG solution obtained for τ defined by α ≥ 0.2 contains large spurious oscillations whereas they are less important for α ≤ 0.01. These oscillations disappear if α ∈ [0.05, 0.1]. The results shown in Table 5.3 confirm what has been observed in Fig. 5.10. These results are similar as those obtained with SUPG P 1 method.

Accuracy study

To assess the quality of the numerical approximations, we compare the approximate solutions with the exact solution in terms of L 2 -norm. The errors of the numerical approximations for various number of d.o.f. and choice of the coefficient α are shown in Table 5.3 for the RK 2 time discretization and in Table 5 The L 2 -error norm in function of the choice of the stabilization parameter α.

IGA: A B-SPLINE BASED APPROACH

B-spline mapping transforms a patch of multiple elements in the parameter space into the physical space.

Each element in the physical space is the image of a corresponding element (knot interval) in the parameter space, but the mapping itself is global to the whole patch, rather than to the elements themselves. We refer to Fig. 6.1 for a schematic overview of this approach. Figure 6.1: An example of a B-spline patch in physical space Ω, parametric space Ω, and the reference element Ω used to perform numerical integration.

Isogeometric discretisation

The parameter domain Ω is defined by the knot vectors Ξ 1 and Ξ 2 :

and

where, p and q are prescribed degrees, ξ i and η j are the i -th and j -th knots, i , j are the knots indices, i = 1, 2, ..., n + p + 1, j = 1, 2, ..., m + q + 1 and n, m equals the number of basis functions. These results are very similar to those obtained with Lagrange quadratic elements in Fig. 5.10.

Error estimates for the quadratic B-spline

We define the error of the numerical solution as e(T ) = u ex (x, T )-u h (x, T ). The L 2 (Ω)-norm of the error is depicted in the Table 6.1 which is defined by: Convergence results in the L 2 -norm are shown in Tab. 6.1 and Tab. 6.2 for the quadratic B-spline. These results are also visualized in Fig. 6.4. As a result, the continuity of the basis is C 1 everywhere, for the quadratic B-spline case. As can be seen, the L 2 -convergence for quadratic B-spline in conjunction with the RK 4 for the temporal discretization is approximately 3. This is the best one could reasonably hope for.

Higher order B-spline

Table 6.3 shows the L 2 -error of the numerical approximations for various number of control point n and order of B-spline. Those results are depicted in Figure 6.4. These results have been obtain without stabilization.

Table 6.4 shows the convergence rates observed for different degrees. As can be seen, a sub-optimal rate of the value 2 is obtained for quadratic basis in conjunction with RK 2, certainly due to the stabilization. For quadratic basis in conjunction with RK 4, the optimal rate of the value 3 is obtained. For the degree p = 4, the rate is limited by the use of RK 4 time integrator. For p = 3 a sub-optimal rate of value 3.5 is also observed.

IGFEA and classical FEA: a comparisons

In this section we present a comparison of the performance of IGFEM and the classical FEM for the numerical simulation of the one-dimensional advection problem.

These two methods employ the same mathematical foundations. Therefore, they have many similarities.

However, some important differences lie in the choice of basis functions. Obviously, in classical FEM the basis which is chosen to approximate the unknown field is interpolatory. This often takes the form of polynomial functions and the geometry is in most cases only approximated.

Whereas the IGA approach has an advantage that the basis is chosen to exactly capture the geometry and this is also used to approximate the field of unknown quantities.

Another benefit is that the approximation is smooth. In fact, when the multiplicity at a knot is m, for 1 ≤ m ≤ p in IGA, the basis functions are then C p-m at the interfaces of the involved elements, whereas basis functions based on Lagrange interpolating polynomials are only C 0 .

Also, it is important to point out that, due to the support of a B-spline, function of order p is always p + 1 knot spans. Therefore, a higher-order B-spline function has support over much larger portion of the domain when compared to classical FEM. Therefore, the computational efficiency is reduced.

We proceed now further to numerically compare the IGFEM and FEM. A good comparison we might perform is to compare the numerical approximations error with respect to the number of d.o.f. However, a special attention should be paid to the accuracy of the SUPG stabilization method.

As mentioned before, the linear FEM and IGFEM give the same results (due to the equivalence of hat shape functions and linear B-spline functions).

Note that the error for the one-dimensional advection problem is plotted in Fig. 6.5 as a function of the global d.o.f. in conjunction with SUPG stabilization (if necessary), to be able to compare IGFEM and FEM. As one can see, results of this simple numerical test allow to make some conclusions. At first sight, the graphs of Fig. 6.5 indicate that the slope of the error lines are not exactly 1, 2, 3 and 4 for 1st , 2nd , 3r d and 4t h order elements, respectively. Moreover, we highlight that the values of the errors from the IGFEM are far lower than those of the FEM.

We also note that, increasing the order of Lagrangian polynomials may increase the amplitude of oscillations in the FEM. This problem is eliminated in IGA as a result of non-negativity and non-interpolatory nature of the B-splines shape functions. In conclusion, we believe that the isogeometric approach has considerable potential in practical problem solving and is a promising alternative to current analysis procedures.

However, a possible drawback of the SUPG method is the sensitivity of the solution to the stabilization parameter, whose optimal value is not determined precisely by the available theory.

In the context of IGA, another approach we could go for is the application of DG method. Indeed, the major argument for using DG methods lies in their ability to provide stable numerical methods for first order PDE problems, for which classical FEM is well known to perform poorly. Therefore, we will introduce in the next chapter a new DG method in the IG context, called IGDGM.

DISCONTINUOUS GALERKIN METHOD (DGM): FROM CLASSICAL TO ISOGEOMETRIC

W E propose a method that combines isogeometric analysis (IGA) with the discontin- uous Galerkin (DG) method for solving hyperbolic problems, namely the isogeometric discontinuous Galerkin (IGDG) method that merges exact geometry with high-order solution accuracy [54] [61]. In this chapter we formulate and analyze this method for the onedimensional advection problem. The solution of the problem is approximated in every sub-domain without any continuity requirement for the discrete solution at the interfaces. Finally, we numerically compare the performance of the IGDG method with the DGFE method.

Introduction and background

The discontinuous Galerkin finite element (DGFE) method was originally introduced in 1973 by Reed

and Hill [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF], for the numerical solution of the nuclear transport PDE problem. Subsequently, the method has found broad applications in large-scale data intensive science and engineering problems. DG is a class of FEM that uses completely discontinuous basis functions. Thanks to their flexibility in local approximation, they offer good stability properties when approximating convection dominated problems [79] [93].

In contrast to the stabilized continuous Galerkin FEM, DG method produces stable discretizations for hyperbolic problems without the need for stabilization parameters, stabilization resulting from the use of upwind fluxes. Therefore, this method combines the best properties of the finite volume (FV) method and continuous Galerkin FEM. As we can see, Table 7.1 demonstrates that greater accuracy can be achieved by either increasing the polynomial order or increasing the number of elements in the domain. The convergence rates were estimated by fixing the polynomial degree, and measuring the L 2 -error as a function of the number of elements.

Assuming a relationship of the form: In practice, basis functions used to represent the solution in the physical domain are defined in the parametric domain by:

Numerical studies

Now, we will revisit the one-dimentional advection problem (5.1) by using the developed IGDG method.

We will give the numerical results demonstrating the performance of the IGDG method in conjunction with RK time discretization.

We start with an initial B-spline patch which is formed on Ω = [-1, 1]. This uniform patch is used to solve the test problem. The solution is computed up to T = 0.4. The initial patch is composed of equal knot intervals of size h, and is split into a set of Bézier elements as explained in section 7.6.1 to apply the DG formulation.

CONCLUSION

L A méthode de Galerkin discontinue a initialement été introduite par Reed et Hill [74] en 1973, pour la discrétisation des équations caractérisant le transport de neutrons.

Numerical Lax-Friedrichs fluxes

In the IGDG method, continuity is not enforced between elements. The flux f n (u e ) = f e n = -→ c u e • -→ n e along the boundaries Γ e must be approximated by a numerical flux f e n . In the present case for the Bézier elements, we use the local Lax-Friedrichs flux which can be defined by:

Due to the fact that the weak form of the DG method is written elementwise, the numerical flux between adjacent elements must be defined. For this purpose, it is possible to write for each element D e :

Cartesian grids

The first mesh type considered are simple cartesian grids. We consider p = q ∈ {1, 2, 3, 4} and asses the quality of numerical approximations through the L 2 -norm.

Curvilinear grids

In this case, we consider N el curvilinear patches. Contrary to the previous cases, the geometry curvature is taken into account here. We underline that the boundary geometry is identically maintained whatever the basis degree (except for p = 1) and d.o.f number. In the following we investigate the numerical order of convergence in the the L 2 -norm of the IGDG discretizations, the problem is solved for several different values of bivariate Bernstein polynomial degrees (p, p) and numbers of patches N el . The L 2 -norm of the error and convergence rate are shown in Table 8.5 for the linear and quadratic case and in Table 8 Again, as can be seen, we obtain that the L 2 convergence rate is approximately p + 1.

C H A P T E R

2D ACOUSTIC WAVE EQUATIONS

T He purpose of this chapter is to develope and analyze the new IGDG method, that uses the IGA discretization concept combined with the DG technique, for solving the first-order acoustic wave equation in 2D, modelling sound propagation phenomena. The computational domain is divided into non-overlapping sub-domains, composed of B-spline patches. The DG approach was applied on element level, each element being a Bézier patch constructed from initial B-spline patches. The solution of the problem is approximated in every element without imposing any continuity requirements for the discrete solution on the interfaces. Basic tests of accuracy and stability are demonstrated, including optimal convergence rates with respect to L 2 -norm .

Introduction and basic theory

Acoustic equations model acoustic wave propagation in a medium. Several application fields are covered by such a model, like elastic wave propagation in the ground, or sound propagation in the air [START_REF] Rienstra | An introduction to acoustics[END_REF]. The form itself is usually considered as a linear problem governed by the compressible linearized Euler equations, in order to describe mean flow effects on sound propagation, such as refraction.

The pressure-velocity formulation of the acoustic wave equations is expressed as a linear PDE system, in which the acoustic pressure and velocity interact with one another to propagate waves through materials.

Let Ω ⊂ R 2 be a two dimensional domain with boundary ∂Ω. On each point on ∂Ω we denote by n the outward normal vector. Let T > 0 be a fixed time.

Numerical results

We consider the example of an ideal acoustic resonator between two cylinders solved using the acoustic wave equation (9.3) over the physical domain Ω represented in Fig. 9.1, for which an analytical solution exists.

Two types of computational domains are considered: a curvilinear one constructed by least squares approximation and knot insertion (described in chapter 4) and a rectilinear one, that corresponds to classical grids with straight element interfaces. Numerical examples for advection problem were given for both classical Lagrange and B-spline bases.

It was shown that the IG method may be well suited to the hyperbolic problems if stabilization is correctly calibrated. The convergence rates were measured in the L 2 -norm. The IG method was tested on different polynomial orders. It was found that, for a given polynomial order, high regularity gives lower error versus degrees of freedom compared to low regularity, without exceptions. Elsewhere, a possible drawback of the SUPG method is the sensitivity of the solution to the stabilization parameter, whose optimal value is not determined precisely by the available theory and is tedious to select in practice.

Therefore, we proposed in this thesis a new method which combines the IGA with the DG method, called IGDG method, for solving hyperbolic problems. The major reason for using DG methods in this thesis lies in their ability to provide stable numerical methods for hyperbolic problems, for which classical FEM is well known to perform poorly.

Our method takes advantage of both IGA and DG methods. In fact, DG methodology is adopted at Bézier patch level, i.e., we employ the traditional IGA within each Bézier patch, and employ the DG method across the patch interfaces to glue the multiple patches. Bézier patches, considered as elements, are constructed by transformation of the initial B-spline domain (Bézier extraction is a classical CAD technique).

We consider then scalar and system of conservation laws, as test cases. These test examples show that the resulting IGDG method in conjunction with RK method is stable, high-order accurate, and can easily handle curved geometries and boundary conditions. As consequence, this method can be easily formulated and implemented. Then, the numerical behavior of the method has been evaluated and it was shown an optimal convergence rate for 2D advection and acoustic problems. The higher accuracy obtained by using curvilinear elements has been demonstrated, in particular when coarse grids are employed. An accuracy up to order five has been reached.

The present method has shown to be very effective for linear hyperbolic problems and systems in one and two dimensions. Now, the proposed IGDG method should be extended to more complex non-linear

B.2 1D slope limiting

If the solution of the problem exhibits discontinuities, the proposed scheme generates oscillations (Gibbs phenomenon). To overcome this difficulty, a particular treatment is required, such as filtering or limiting. In the present work, we envisage to apply the generalized limiting approach. The idea is to modify locally the solution, if a discontinuity is detected, in order to satisfy the following conditions:

i ) Maintain the mass conservation principle in each element.

i i ) Satisfy the Total Variation Diminishing in the Means (TVDM) property.

i i i ) Do not degrade accuracy of the method.

B.2.1 TVDM limiter

The proposed approach, relies on the minmod function defined as:

Let u j denote the mean value of the solution over the element Ω j . Thanks to the partition of unity property of the Bézier representation, it can be simply evaluated by:

The minmod function is employed to extrapolate linearly the mean value u j at the element extremities, according to a limited slope estimation, using:

B.2.2 TVBM limiter

The TVDM generalized slope limiter described above yields a loss of accuracy in the vicinity of local extrema of the solution. To prevent this effect, one must construct a Total Variation Bounded in the Means (TVBM) limiter, instead of the TVDM one.

This can easily be achieved by introducing a modified minmod function (less restrictive):