Ahmed Ratnani
email: ahmed.ratnani@cea.fr

Πgasus : Python for IsoGeometric AnalysiS and Unified Simulations

B-splines and NURBS (Non Uniform Rational B-splines) are widely used in CAD (Computer Aider Design) models. IGA (IsoGeometric Analysis) consists of using these functions to both define the geometry and represent the unknowns that are solution of a Partial Differential Equation, using the Finite Element principle. In this paper we present a new library, namely Πgasus , that was developped in order to bring a common framework between the users (especially physicists) and mathematicians. We want to provide a stable and robust framework, that handles complex geometries and models as it is the case in Plasma Physics. Physicists will be able to use the recent works and results obtained by mathematicians. Πgasus is a 1D, 2D and 3D Fortran code, interfaced with Python. It provides a Geometry module, a FEM (Finite Element Method) computational engine and a Visualization module.

Introduction

The IsoGeometric Analysis has been introduced by Hughes et al [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF], in order to bring a common framework between CAD (Computer Aided Design) and Numerical Simulations. During the whole process, including refinement and analysis, the geometry is maintained exact. In addition to the classical refinement strategies, hp-refinement, IGA introduces a new one, namelly the k-refinement. The latter offers the possiblity to control the regularity of the basis functions, by increasing or decreasing the multiplicities of knots. From this point of view, IGA can be seen as a generalization of standard FEMs where more regular functions are used.

In addition to handle complex and CAD geometries, it has been shown [START_REF] Beirao Daveiga | Some estimates for h-p-k refinement in isogeometric analysis[END_REF] that, the use of regular elements reduces the dimension of the Finite Element spaces, while keeping the desired precision. It allows to deal with higher order differential operators [START_REF] Gomez | Isogeometric analysis of the isothermal navierøstokesøkorteweg equations[END_REF]. Moreover, it has been noticed [START_REF] Ratnani | An arbitrary high-order spline finite element solver for the time domain maxwell equations[END_REF], that the use of regular elements tends to give better CFL numbers.

For a comprehensive introduction on the subject, we may refeer the reader to the recent book by Cottrel et al [START_REF] Cottrell | Isogeometric Analysis, toward Integration of CAD and FEA[END_REF]. For an algorithmic overview on B-splines we refeer to the book by DeBoor [START_REF] Deboor | A practical guide to splines[END_REF] or Piegl and Tiller [START_REF] Piegl | The NURBS book[END_REF]. For an extensive overview on the approximation theory, we may suggest the books [START_REF] Devore | Constructive Approximation[END_REF][START_REF] Schumaker | Spline Functions: Basic Theory[END_REF].

IGA is getting a large success inside the Engineering and Mathematical communities. The increasing number of publications shows the interest of these communities to this new approach. It has been used in electromagnetism [START_REF] Buffa | Isogeometric analysis in electromagnetics: B-splines approximation[END_REF][START_REF] Buffa | Isogeometric discrete differential forms in three dimensions[END_REF][START_REF] Ratnani | An arbitrary high-order spline finite element solver for the time domain maxwell equations[END_REF], in incompressible fluid dynamics [START_REF] Buffa | Isogeometric analysis: Stable elements for the 2d stokes equation[END_REF], in fluid-structure interaction [START_REF] Hsu | Blood vessel tissue prestress modeling for vascular fluidøstructure interaction simulation[END_REF][START_REF] Bazilevs | Isogeometric fluidøstructure interaction analysis with emphasis on non-matching discretizations, and with application to wind turbines[END_REF] or using a mapped Finite Volume method [START_REF] Ch | A finite volume method on nurbs geometries and its application in isogeometric fluidøstructure interaction[END_REF], in structural and contact mechanics [START_REF] Shojaee | Free vibration analysis of thin plates by using a nurbs-based isogeometric approach[END_REF][START_REF] Konyukhov | Geometrically exact theory for contact interactions of 1d manifolds. algorithmic implementation with various finite element models[END_REF]. IGA has been used in Plasmas Physics, in MHD (Magneto-HydroDynamics) problems [START_REF] Ratnani | Isogeometric analysis in reduced magnetohydrodynamics[END_REF], and the Kinetic approach [START_REF] Crouseilles | An isogeometric analysis approach for the study of the gyrokinetic quasi-neutrality equation[END_REF][START_REF] Abiteboul | Solving the vlasov equation in complex geometries[END_REF][START_REF] Back | An axisymmetric pic code based on isogeometric analysis[END_REF]].

An active area of research is the study of local refinement. It seems that T-splines are getting more popular. Many studies are under investigation to understand the behavior of the discrete spaces [START_REF] Berdinsky | On the problem of instability in the dimension of a spline space over a t-mesh[END_REF][START_REF] Beirào Da Veiga | Analysis-suitable t-splines are dualcompatible[END_REF] or to glue 2 patchs [START_REF] Beirào Da Veiga | Isogeometric analysis using t-splines on two-patch geometries[END_REF], which is of a big importance for parallelization.

Another active area is how to construct a volume description of a computational domain, given by its boundary, so we can directly use all existing CAD-models. For those interested in this subject, we may refeer to the works [START_REF] Martin | Volumetric parameterization and trivariate b-spline fitting using harmonic functions[END_REF][START_REF] Martin | Volumetric parameterization of complex objects by respecting multiple materials[END_REF][START_REF] Martin | Mixed-element volume completion from nurbs surfaces[END_REF][START_REF] Cohen | Analysis-aware modeling: Understanding quality considerations in modeling for isogeometric analysis[END_REF][START_REF] Xu | Optimal analysis-aware parameterization of computational domain in 3d isogeometric analysis[END_REF][START_REF] Xu | Parameterization of computational domain in isogeometric analysis: Methods and comparison[END_REF].

Some useful softwares can be found ( [START_REF] Vuong | Isogat: A 2d tutorial matlab code for isogeometric analysis[END_REF], [START_REF] De Falco | Geopdes: A research tool for isogeometric analysis of pdes[END_REF]) to accompany researchers and introduce them to this new subject. However, for real and complex applications, they present some disadvantages. They are inflexible for Numerical Simulations combined with Modeling. In the latter, the user may want to see the effect of adding a new term, or remove it, test a new time scheme, etc .... In order to solve partial differential equations (PDEs), one needs to copy the old written code, and then change it. Hence, we will end up with as many as codes, as the problems we have treated, even if the principal changes occure while assembling the element matrices.

All these points, motivated me to developp a new library, Πgasus , written in Fortran and interfaced with Python. The idea was to simplify as much as possible this interface, so we can recover the classical language of the Finite Element Method, while maintaining the advantages of a scriptedlanguage. By sharing the library between different projects, we will keep the same computational engine, and so reduce redundancy. In addition, the architecture of Πgasus offers the user the ability to manipulate the geometry and refine it (using the hpk-refinement strategies), use high order and regular elements, and to manipulate the differential operators even for non-linear problems. The design of Πgasus gives the user some freedom to define the spaces and manipulate the grids (thus relax the isoparametric approach, by considering edge, surface, . . . elements).

Πgasus is not intend (for the moment) to be a concurrent of automated FEM softwares like Freefem++ or Feelpp [START_REF] Pironneau | Freefem++[END_REF][START_REF] Prudh Homme | Feel++: A computational framework for galerkin methods and advanced numerical methods[END_REF].

This work was motivated by the increasing interactions between physicists and mathematicians in (computational) plasma physics, trying to answer the question : how can a mathematician study and improve physicists works and challenges, and make his developpements promptly available? It was initiated during my Phd-thesis [START_REF] Ratnani | Isogeometric analysis in plasma physics and electromagnetism[END_REF] at the INRIA, and pursued at the CEA. We wanted to provide both physicists and mathematicians a common framework where they can developp their researchs, take advantage from the IGA approach, but also in the future for production cases.

As the reader will notice it, the user does not need to have a deep knowledge on CAD or Bsplines. The important point is to know how to derive a variational (weak) formulation, and then ask Πgasus to assemble the different parts of the formulation. Πgasus is delivered with a (2D) simple and easy GUI interface to create and manipulate the computational domains. The user can move and manipulate the control points in order to generate the geometry. The basic distribution of Πgasus does not need any third package or library.

Throughout this article, we will present 2D studies even if the code treats also the 1D and 3D cases. The current paper is structured as follows. In section 2, we give a brief overview on B-splines, NURBS and IsoGeometric Analysis approach. In section 3, we present Πgasus on a simple example (Poisson's equation). We have made an effort to link each of Πgasus 's functions to the FEM language. In section 4, we give much more details on Πgasus : its architecture, the important concepts, utilities and features. Finally, in section 5, we a non-linear 2D pde as a simple application. This article is not a tutorial, but we have tried to present as mush as possible the different and interesting parts.

An overview on IsoGeometric Analysis

More details on the subject can be found in [START_REF] Devore | Constructive Approximation[END_REF][START_REF] Schumaker | Spline Functions: Basic Theory[END_REF][START_REF] Piegl | The NURBS book[END_REF][START_REF] Deboor | A practical guide to splines[END_REF][START_REF] Cottrell | Isogeometric Analysis, toward Integration of CAD and FEA[END_REF].

Basic properties of B-splines

Definition 2.1 (B-Spline) Let X = {x 0 , • • • , x p } a non-decreasing sequence of p + 1 points such that x 0 = x p . The B-Spline is defined by the following reccurence formula:

N (x; x 0 , • • • , x p ) = x -x 0 x p-1 -x 0 N (x; x 0 , • • • , x p-1 ) + x p -x x p -x 1 N (x; x 1 , • • • , x p ) (2.1)
the initialization is given by : N (x; x 0 , x 1 ) = 1, if x 0 = x 1 , and 0 otherwise.

In order to construct a B-spline of degree p, we need p+1 points which are called knots in the spline terminology. Then, to create a family of B-splines, we will need to have a non-decreasing sequence of knots, also called knot vector.

Let T = (t i ) 1 i N +k be a non-decreasing sequence of knots, with k = p + 1. Each set of p + 1 consecutive knots T j = {t j , • • • , t j+p } will generate a B-spline N j . This leads to the following definition: Definition 2.2 (B-Spline series) The j-th B-Spline of order k is defined by the recurrence relation:

N k j = w k j N k-1 j + (1 -w k j+1 )N k-1 j+1
where,

w k j (x) = x -t j t j+k-1 -t j N 1 j (x) = χ [t j ,t j+1 [ (x)
for k ≥ 1 and 1 ≤ j ≤ N .

In figures 1 and 2, we give examples of the generated B-splines families depending on the knot vectors (figure 1) and the order (figure 2). We note some important properties of a B-splines basis:

• B-splines are piecewise polynomial of degree p = k -1, 0 0.5 • Compact support; the support of

N k j is contained in [t j , t j+k ], • if x ∈ ]t j , t j+1 [, then only the B-splines {N k j-k+1 , • • • , N k j } are non vanishing at x, • Positivity: ∀j ∈ {1, • • • , N } N j (x) > 0, ∀x ∈]t j , t j+k [, • Partition of unity : N i=1 N k i (x) = 1, ∀x ∈ R,
• Local linear independence, and finally

• If a knot t i has a multiplicity m i then the B-spline is C (p-m i ) at t i

Multivariate tensor product splines

Let us consider d knot vectors

T = {T 1 , T 2 , • • • , T d }.
For simplicity, we consider that those knot vectors are open, which means that k knots on each side are duplicated so that the spline is interpolating on the boundary, and of bounds 0 and 1. In the sequel we will use the notation I = [0, 1]. Each knot vector T i , will generate a basis for a Schoenberg space [START_REF] Devore | Constructive Approximation[END_REF], S k i (T i , I). The tensor product of all those spaces is also a Schoenberg space, namely S k (T ),

where k = {k 1 , • • • , k d }.
The hypercube P = I d = [0, 1] d , will be referred to as a patch.

The basis for S k (T ) is defined by a tensor product :

N k i := N k 1 i 1 ⊗ N k 2 i 2 ⊗ • • • ⊗ N k d i d where, i = {i 1 , • • • , i d }.
A typical cell from P is a cube of the form :

Q i = [ξ i 1 , ξ i 1 +1 ] ⊗ • • • ⊗ [ξ i d , ξ i d +1 ].
To any cell Q, we will associate its extension Q, which is the union of the supports of basis functions, that intersects Q.

Splines in CAD

In order to have a control on the regularity of the curve, we need to use a piecewise-polynomial form. This is why the use of B-splines has known a large success. Moreover, the control points (and polygone) have many geometric interpretations. Let (P i ) 1 i N ∈ R d be a sequence of control points, forming a control polygon.

Definition 2.3 (B-Spline curve)

The B-spline curve in R d associated to T = (t i ) 1 i N +k and (P i ) 1 i N is defined by :

C(t) = N i=1 N k i (t)P i
We have the following properties for a B-spline curve:

• If N = k, then C is just a Bézier-curve,
• C is a piecewise polynomial curve,

• The curve interpolates its extremas if the associated multiplicity of the first and the last knot are maximum (i.e. equal to k),

• Invariance with respect to affine transformations,

• Strong convex-hull property:

if t i ≤ t ≤ t i+1 , then C(t) is inside the convex-hull associated to the control points P i-p , • • • , P i ,
• Local modification : moving P i affects C(t), only in the interval [t i , t i+k ],

• The control polygon approaches the behavior of the curve.

Remark 2.4

Remark that there is a kind of duality between knots and control points. We can use multiple control points : P i = P i+1 , instead of multiple knots.

Deriving a B-spline curve:

We have:

C (t) = n i=1 N k i (t)P i = n i=1 p t i+p -t i N k-1 i (t)P i - p t i+1+p -t i+1 N k-1 i+1 (t)P i = n-1 i=1 N k-1 i * (t)Q i (2.
2) where Q i = p P i+1 -P i t i+1+p -t i+1 , and {N k-1 i * , 1 ≤ i ≤ n -1} are generated using the knot vector T * which is obtained from T by reducing by one the multiplicity of the first and the last knot (in the case of open knot vector), i.e. by removing the first and the last knot.

Example: T = {000 2 5 3 5 111}, p = 2, n = 5. We have C(t) = 5 i=1 N 3 i (t)P i , then

C (t) = 4 i=1 N 2 i * (t)Q i where Q 1 = 5{P 2 -P 1 }, Q 2 = 10 3 {P 3 -P 2 }, Q 3 = 10 3 {P 4 -P 3 }, Q 4 = 5{P 5 -P 4 }.
The B-splines {N 2 i * , 1 ≤ i ≤ 4} are associated to the knot vector T * = {00 2

(P i,j ) 1 i N 1 ,1 j N 2 , is defined by M(ξ, η) = N 1 i=1 N 2 j=1 N i,j (ξ, η)P i,j with N i,j (ξ, η) = N (1) i (ξ)N (2) j (η)

Fundamental geometric operations

We can change the curve's parameters without changing the curve. This refers as to geometric operations. For example, we can elevate the spline's degree or insert new knots. Many algorithms are available from the CAD community, and have proved their efficiency [START_REF] Deboor | A practical guide to splines[END_REF][START_REF] Goldman | Knot Insertion and Deletion Algorithms for B-Spline Curves and Surfaces[END_REF][START_REF] Huang | Fast degree elevation and knot insertion for b-spline curves[END_REF][START_REF] Prautzsch | A fast algorithm to raise the degree of spline curves[END_REF][START_REF] Tiller | The NURBS Book[END_REF].

In figure 4, we show the use of some of these geometric operations (elevation degree and knot insertion) on a quadratic B-Spline curve. We show also how one can split a domain by raising the multiplicity of knots to the spline's degree. This is an important strategy that may help us for domain decomposition and parallelization.

In figure 5, we show the use of the elevation degree and insertion knots algorithms on a 2D domain.

As we can see, the geometry is kept unchanged.

NURBS

Let ω = (ω i ) 1 i N be a sequence of non-negative reals. The NURBS (Non-Uniform Rational Bsplines) functions are defined by a projective transformation:

Definition 2.6 (NURBS)

The i-th NURBS of order k, associated to the knot vector T and the weights ω, is defined by

R k i = ω i N k i N j=1 ω j N k j . (2.3)
Notice that when the weights are equal to 1 the NURBS are B-splines.

Definition 2.7 (NURBS curve)

The NURBS curve of order k associated to the knot vector T , the control points (P i ) 1 i N and the weights ω, is defined by

M(t) = N i=1 R k i (t)P i .

Modeling conics using NURBS

In this section, we will show how to construct an arc of conic, using rational B-splines. Let us consider the following knot vector : T = {000 111}, the generated B-splines are Bernstein polynomials. According to 2.3, the general form of a rational Bézier curve of degree 2 is:

C(t) = ω 1 N 2 1 (t)P 1 + ω 2 N 2 2 (t)P 2 + ω 3 N 2 3 (t)P 3 ω 1 N 2 1 (t) + ω 2 N 2 2 (t) + ω 3 N 2 3 (t) (2.4)
Let us consider the case ω 1 = ω 3 = 1. Because of the multiplicity of the knots 0 and 1, the curve C is linking the control point P 1 to P 3 . Depending on the value of ω 2 , we get different type of curves (Table 1). 4 , 1 2 , 3 4 }, this leads to a Bézier description. We can then, split the curve into 4 pieces (sub-domains), each one will corresponds to a quadratic Bézier curve.

Multivariate tensor product NURBS

As for splines, one can define multivariate tensor product NURBS. For surfaces, we have the following definition.

Definition 2.8 (NURBS surface)

The NURBS surface of order k associated to the knot vectors {T (1) , T (2) }, the control points (P i,j ) 1 i N 1 ,1 j N 2 and the weights (ω i,j ) 1 i N 1 ,1 j N 2 , is defined by M(t (1) , t (2) (1) , t (2) )P i,j with R i,j (t (1) , t (2)

) = N 1 i=1 N 2 j=1 R i,j (t
) = ω i,j N 1 i N 2 j r,s ωr,sN 1 r N 2 s
Remark 2.9 NURBS functions inherit most of B-splines properties. Remark that in the interior of a knot span, all derivatives exist, and are rational functions with non vanishing denominator. We present here the definition of the perspective mapping. We construct the weighted control points P ω i = (ω i x i , ω i y i , ω i z i , ω i ), then we define the B-spline curve in four-dimensional space as

nature of the curve

ω 2 = 0 line 0 < ω 2 < 1 ellipse arc ω 2 = 1 parabolic arc ω 2 > 1 hyperbolic arc
M ω (t) = N i=1 N k i (t)P ω i . (2.5)
For fundamental geometric operations on NURBS curves, we use the latest transformation and algorithms on B-spline curves.

Remark 2.10 NURBS functions allow us to model, exactly, much more domains than B-splines. In fact, all conics can be exactly represented with NURBS. For more details, see [START_REF] Tiller | The NURBS Book[END_REF].

IsoGeometric Analysis

The idea behind the IGA method is to use the same functions that define the physical domain, to approach the solution of a partial differential equation.

In the sequel, we consider 2 knot vectors

T ξ = {ξ 1 , • • • , ξ N 1 +p 1 +1 } and T η = {η 1 , • • • , η N 2 +p 2 +1 }. Let W ξ = {ω ξ 1 , • • • , ω ξ N 1 } and W η = {ω η 1 , • • • , ω η N 2
} be two weight sequences, and (P ij ) 1 i N 1 ,1 j N 2 a sequence of control points. This defines a mapping

F(ξ, η) = 1 i N 1 ,1 j N 2 R ξ i (ξ)R η j (η)P ij (2.6) that maps the rectangular patch [ξ 1 , ξ N 1 ] × [η 1 , η N 2 ] onto the physical domain Ω.
Where R ξ and R η are NURBS functions defined by knot vectors T ξ and T η , and weights W ξ and W η .

As said before, we consider only open knot vectors. Without loss of generality, we shall consider knot vectors of the form:

ξ 1 = • • • = ξ p 1 +1 = η 1 = • • • = η p 2 +1 = 0, and ξ N 1 +1 = • • • = ξ N 1 +p 1 +1 = η N 2 +1 = • • • = η N 2 +p 2 +1 = 1.
Let K be a cell in the physical domain. Q is the parametric associated cell and such that K = F(Q). Let J F be the Jacobian of the transformation F, that maps any parametric domain point (ξ, η) into the physical domain point (x, y) (figure 6).

Refinement strategies

Refining the grid can be done in 3 different ways. This is the most interesting aspects of B-splines basis, • using the patch parameter h, by inserting new knots. This is the h-refinement, it is the equivalent to mesh refinement of the classical Finite Element Method.

• using the degree p, by elevating the B-spline degree. This is the p-refinement, it is equivalent to use higher finite element order in the classical FEM.

• using the regularity of B-splines, by increasing / decreasing the multiplicity of inserted knots. This is the k-refinement. This new strategy does not have an equivalent in the classical FEM.

An active area of research is the study of local refinement. It is important to notice that the use of tensor products leads to the existence of a lot of superfluous control points, that might exist because of the presence of a cartesian grid in the parametric domain. Sederberg et al. [START_REF] Sederberg | T-spline simplification and local refinement[END_REF] defined the notion of T-splines that allows us to reduce the number of those control points. In [START_REF] Michael | Adaptive isogeometric analysis by local hrefinement with t-splines[END_REF] D örfel et al. use T-splines for local h-refinement in isogeometric analysis. It seems that T-splines are getting more popular. Many studies are under investigation to understand the behavior of the discrete spaces [START_REF] Berdinsky | On the problem of instability in the dimension of a spline space over a t-mesh[END_REF][START_REF] Beirào Da Veiga | Analysis-suitable t-splines are dualcompatible[END_REF] or to glue 2 patchs [START_REF] Beirào Da Veiga | Isogeometric analysis using t-splines on two-patch geometries[END_REF], which is of a big importance for parallelization. For more information on this subject, we refer to [START_REF] Schillinger | An isogeometric design-through-analysis methodology based on adaptive hierarchical refinement of nurbs, immersed boundary methods, and t-spline cad surfaces[END_REF][START_REF] Schillinger | An unfitted hp-adaptive finite element method based on hierarchical b-splines for interface problems of complex geometry[END_REF][START_REF] Nguyen-Thanh | Isogeometric analysis using polynomial splines over hierarchical t-meshes for two-dimensional elastic solids[END_REF][START_REF] Scott | Local refinement of analysis-suitable tsplines[END_REF][START_REF] Buffa | Characterization of t-splines with reduced continuity order on t-meshes[END_REF] or [START_REF] Wang | Adaptive isogeometric analysis using rational pht-splines[END_REF][START_REF] Tian | Adaptive finite element methods for elliptic equations over hierarchical t-meshes[END_REF] for PHT-splines. For details on implementation see [START_REF] Rypl | Object oriented implementation of the t-spline based isogeometric analysis[END_REF]. In [START_REF] Vuong | A hierarchical approach to adaptive local refinement in isogeometric analysis[END_REF], authors propose the use of a hierarchical local refinement method. THB-splines are used in [START_REF] Giannelli | Thb-splines: The truncated basis for hierarchical splines[END_REF]. Dokken, Kvamsdal and their team from SINTEF are currently developing another approach for local refinement, based on LR-splines [START_REF] Dokken | Non-Standard Numerical Methods for PDE's[END_REF].

Patch

Starting from the geometric description of the domain, which gives a coarse mesh, we can then, use h/p/k refinements to create the computational grid. We can also use multiple patches to describe more complex domains [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF][START_REF] Cottrell | Isogeometric Analysis, toward Integration of CAD and FEA[END_REF]. There are many ways to stick those patches together.

Grid generation

For this purpose, we use alternatively h and p-refinement. The minimal degree of basis functions is imposed by the domain design. When inserting knots, we can use uniformly-spaced knots or non uniformly-spaced ones. In figure 7, we give an example of such a refinement.

Introducing Πgasus

In this section, we will consider the resolution of an elliptic partial differential equation on a square domain Ω. We will see how we can recover the classical FEM language, for the IGA approach. We also will separate the different differential operators in the equation and how we can construct them using Πgasus . Let us consider the following problem:

For given functions A, f, b, find u such that:

-∇ • (A∇u) + bu = f , Ω u = 0 , ∂Ω (3.7)
Introducing the matrix :

Σ = (Σ i,j ) 1 i,j n
and the vectors :

L = (L i ) T 1 i n-n D [u h ] = ([u h ] i ) T 1 i n
where, we denote for i, j ∈ {1, .., n}

Σ i,j = Ω ∇ϕ i A ∇ϕ j dΩ + Ω b ϕ i ϕ j dΩ L i = Ω f ϕ i dΩ
here, n denotes the dimension of the discrete space and n D the number of Dirichlet elements.

Using a variational formulation and the Green's formulae, we know that after discretization, we have to solve the linear system:

Σ [u h ] = L .
The matrix Σ is nothing else but the sum of the Stiffness matrix S and Mass matrix M .

S = Ω ∇ϕ i A ∇ϕ j dΩ M = Ω b ϕ i ϕ j dΩ
In practice, the user will need to define, discretize and assemble the two differential operators M and S. To do so, we must start by defining the discrete FE space (see subsection 3.1), with the appropriate boundary conditions and set the computational grid. Then, he needs to define two fields, one for the right hand side term, and the other for the solution of the pde (see subsection 3.2). Finally, he defines the discretized differential operators M and S (see subsection 3.3). In the following subsections, we will detail and explain each step of a typical Πgasus script. For instance, we will solve the Poisson's equation on the unit square domain with the source term :

f (x, y) = 8π 2 sin(2πx) sin(2πy)
and the analytical solution u(x, y) = sin(2πx) sin(2πy)

Spaces Definition

The Galerkin Finite Element Method relies on the fact that we approach the inifinite dimensional spaces, where the unknowns may live, by a sequence of finite dimensional subspaces. With Πgasus we recover this approach and start by creating these finite dimensional subspaces.

# d e f i n e t h e d i s c r e t e space V = space ( as_file=ls_domain ) # s e t boundary c o n d i t i o n s # here we use Homogeneous D i r i c h l e t on whole t h e boundary

V . dirichlet ( faces= [ [ 1 , 2 , 3 , 4 ] ] ) V . set_boundary_conditions ( )
The discrete space needs a geometry that may be given by in XML or HDF5 file. Once the space was defined, we have to associate a grid (quadrature points).

# c r e a t e a g r i d ( using Gauss-Legendre quadrature p o i n t s ) V . create_grids ( type=" l e g e n d r e " , k=lpi_glorder )

Fields Definition

The Field object can be used to define either the right hand side term or the unknown.

Matrices Definition

Here, we show how we can construct and assembly the differential operators involved in our PDE. The Mass matrix operator is of the form Mass

S_V = matrix ( spaces= [ V , V ] , ai_type=STIFFNESS , func=func_stiff )

Initialization

Once we have defined all the operators and fields that we need to model our PDE, we have to initialize Πgasus .

fe . initialize ( )

Assembling Process

Using the function assembly, we can assemble the operators that we need. We can imagine, that for a given problem, one of the operators will remain unchanged while the others may change at each iteration. In this case, the user can specify which matrices or fields to assemble, by given it in the matrices or fields lists.

fe . assembly ( matrices= [ M_V , S_V ] , fields= [ F_V ] )

Solving the linear system

You can either use the integrated linear solver for Πgasus , or export matrices to python and use your favorite solver. Next, we show how the later one can be done using spsolve and SuperLU.

Computing the error norm

In Πgasus , the norms are modeled as operators on fields. The user can compute the classical L 2 , H 1 , etc, . . . , norms for a given field. An example of declaration is:

N_U = norm ( field=U_V , type=NORM_L2 )
Here we show how to compute the L 2 norm.

# s e t t h e NURBS/B-s p l i n e s c o e f f i c i e n t s , # i n t h e c a s e o f using e x t e r n a l s o l v e r U_V . set ( lpr_u ) # assemble t h e norm o p e r a t o r fe . assembly ( norms= [ N_U ] ) lr_norm = N_U . get ( ) As seen, through the previous example, Πgasus is designed following the classical FEM language. In the next section, we shall give much more detail about the architecure of Πgasus and the offered notions and objects.

Dive into Πgasus

In this section, we will introduce some advanced utilities and give much more precision on the design of the library. We start by introducing the Geometry module, the architecture and then present the fundamentals concepts/objects.

Geometry Module

An important difference between CAD and the numerical simulation worlds, is the nature of a given domain. In CAD, a circle is a curve, but in order to solve the Poisson's equation (for example) on a circular domain (for example) we need a 2D description. Hence, even if there are many powerful tools for CAD, they present at least this disadvantage. It is also complicated for mathematicians and physicists to use such tools. Remember that the user wants a volume description of its domain. For more details, we refer to [START_REF] Martin | Volumetric parameterization and trivariate b-spline fitting using harmonic functions[END_REF][START_REF] Martin | Volumetric parameterization of complex objects by respecting multiple materials[END_REF][START_REF] Martin | Mixed-element volume completion from nurbs surfaces[END_REF][START_REF] Cohen | Analysis-aware modeling: Understanding quality considerations in modeling for isogeometric analysis[END_REF]. For all these reasons, we prefered to developp our own module for geometry, with a simple interface for the user and a classical XML format. We have also added HDF5 format to handle heavy data, as the user may want to store his domain after the refinement process.

Formats

There are two ways to store the geometry. For heavy data (especially in 3D), you have to use the HDF5 format. Weights are stored here as the d + 1 component of the control point (where d is the space dimension). In (Appendix, source code 1) we show the XML-file of the linear description of the unit square.

Geometry utilities

In (Appendix, source code 2), we show how to use a p-refinement. Now let's take the linear description of the unit square, and perform 1. a p-refinement: we will elevate the spline degree : +1 in the ξ-direction and +1 in the ηdirection,

2. an h-refinement: we will insert {0.25, 0.5, 0.75} in the ξ-direction and {0.3, 0.7} in the ηdirection,.

The detail of the script is given in (Appendix, source code 3). In order to do a k-refinement, you can either do a h-refinement then a p-refinement, or do directly a hp-refinement with duplicated knots (Appendix, source code 4). In (Appendix, source code 5) we show a simple example of how we can duplicate a given geometry and modify it to have a multi-patchs description.

Remark 4.1

We have developed a 2D-CAD designer tool, in order to manipulate and generate directly, the volume description of the computational domain [START_REF] Ratnani | Caid : A cad tool for isogeometric computational domains[END_REF].

Πgasus architecture

In the following subsection, we present the architecture behind Πgasus . Although there is no Singleton notion in Python, we can define an alternative using the following classical idea. Then, we can define the common obj class, which will contain all spaces, fields, matrices, . . . , declarations. Each time the user will declare a new space (for example), it will be automatically added into the common obj.spaces list.

Remark 4.2

All objects (Fields, Matrices and Norms) are stored with a sub-domain (patch) id. As seen before, the user will use the fundamental geometric operations, to split the computational domain into many subdomains while keeping the geometry exact. Then he will provide this new description to construct the Finite Element space.

Spaces

The notion of spaces has been introduced in the previous section. Let us just give some additional remarks. The user can define a space using an exterior mapping (not defined by the IGA approach). For example, he can give an analytic mapping (or a new metric) or even a CAD-description using splines or nurbs.

Remark 4.3

We have choosed to use a tensorial approach for Πgasus . This means that inside each element, the grid is viewed as n-arrays rather than an array of R n . A patch is a collection of elements. The user can stick patchs together by using the duplicate function. [START_REF] Back | An axisymmetric pic code based on isogeometric analysis[END_REF] The user can provide how to calculate the basis (and their derivatives) on each element. For example, he can give the values of the basis (and their derivatives) over the 1D-quadrature points. Remark 4.5 For (non-homogeneous) boundary conditions, the user can create an adequate space using a grid with a boundary profile.

The dirichlet function is used to set Dirichlet boundary condition on the appropriate domain faces. This is given as a list over patchs (see figure 10). which will set the boundaries defined by {ξ = 0.0, ∀η} and {ξ = 1.0, ∀η}, in the 2D case, to Dirichlet condition.

In the case where we use a multi-patch description, the user can specify the Dirichlet boundary condition as shown in figure 11. The duplicate function can be used to set periodic boundary conditions or to stick patchs together with a C 0 condition. It tells Πgasus that all basis defined by "face-base" have the same ID (see [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF] for details on FEM connectivities structures) as those by "face". Both "face-base" and "face" must be lists of couples [patch id, num face]. To set a periodic boundary condition:

V . dirichlet ( faces= [ [ 1 , 2 , 3 ] , [ 1 , 3 , 4 ] ] ) 1 1 2 3 3 4 4 2 Ω 1 Ω 2
V . duplicate ( faces_base= [ [ 0 , 1 ] ] , faces= [ [ 0 , 3 ] ] )
To stick two patchs together:

V . duplicate ( faces_base= [ [ 0 , 4 ] ] , faces= [ [ 1 , 2 ] ] )

Vectorial spaces

The basic definition of spaces, leads to scalar functions. In order to define a vectorial space, the user must proceed as on paper, he starts by defining spaces for each coordinate, on the patch and without a mapping, then he defines the whole space as a vect space. In the following example, we show a typical construction, for the Maxwell's time domain problem [START_REF] Ratnani | An arbitrary high-order spline finite element solver for the time domain maxwell equations[END_REF]. In 2D domains, Maxwell's equations can be decoupled into two systems. The first involving the (E x , E y , H z ) components is called the Transverse Electric (TE) mode, and the second, involving the (H x , H y , E z ) components is called the Transverse Magnetic (TM) mode. As both modes can be discretized in the same manner, The TE mode reads

∂E ∂t -rot H = -J, (4.8
)

∂H ∂t + rot E = 0, (4.9
)

div E = ρ, (4.10)
where the components are defined by E = E x E y , H = H z . Let us define the following scalar spaces

: V = S p,p α,α , W 1 = S p,p-1 α,α-1 , W 2 = S p-1,p α-1,α
and the vectorial space

W = W 1 × W 2
The discrete spaces V and W are involved in the DeRham sequence [START_REF] Ratnani | An arbitrary high-order spline finite element solver for the time domain maxwell equations[END_REF].

p = 3 n = 63 ls_etiq = " p "+str ( p ) +"x"+str ( p ) +" n "+str ( n ) +"x"+str ( n ) ls_domain_V = "domain p"+str ( p ) +"x"+str ( p ) +" n "+str ( n ) +"x"+str ( n ) +" . xml" ls_domain_W1 = "domain p"+str ( p ) +"x"+str ( p -1 ) +" n "+str ( n ) +"x"+str ( n ) +" . xml" ls_domain_W2 = "domain p"+str ( p -1 ) +"x"+str ( p ) +" n "+str ( n ) +"x"+str ( n ) +" .

Fields

As said previously, the notion of Field is intended to model a right hand side or an unknown. Next we list some of its attributs:

• space : This is the discrete space on which the field is defined.

• func : This is a given function for the field. It can be used either for the source term or the exact solution.

• type : This is the projection/interpolation type.

• operator : This is used to apply an operator to a field, defined previously.

• field : This is the operande for the operator (must be given if the type is FIELD OPERATOR).

• func arguments : This is a list of fields used for a non-linear field (i.e. f ([g], x)).

Remark 4.6

The user can define a field as a function of a list of fields. In the following example, we show how to model a term of the form :

f (u, x, y) = u 2 + x 2 ∂ x u + (x + y)∂ y u
In fact, we will need to define a new field u2 = v • ∇u which is a FIELD OPERATOR of type GRAD applied on the fields u with an argument function v = [x

[ list_F [ 0 ] ˆ2 + list_F [ 1 ] ] F = field ( space=V , func = func_F , func_arguments= [ U , U2 ] ) # . . .

Matrices

A matrix ise stored naturally in a compact format thanks to its profile. The latter is easily computed thanks to the spaces connectivities. In order to construct a matrix, you need to specify the following parameters:

• space : This is a couple of the discrete spaces on which the Differential Operator is defined.

• ai type : This is the type of the Differential Operator (MASS, STIFFNESS, ADVECTION, ...)

• func : This is a given parameter function for the matrix.

• matrices : This is used to construct a composed (block) matrix.

• transpose : If we want to transpose the matrix during the construction.

• func arguments : This is a list of fields used for a non-linear field (i.e. f ([u], x)).

• addto : if we want to add the local element matrix directly to a global matrix previously defined. This will reduce the memory cost.

In table 2, we give the list of the implemented Differential Operators. The user can also create Block-matrices or import a matrix. In the following example we show how to construct a Block-Matrix. The user can choose which part to update during the assembling process.

Operator Contribution

Mass Ω f ([u], x) ϕ b ϕ b dΩ Stiffness Ω (A([u], x) ∇ϕ b ) • ∇ϕ b dΩ Advection Ω v([u], x) • ∇ϕ b ϕ b dΩ Second derivatives Ω v([u], x) • D 2 ϕ b dΩ
M_V = matrix ( matrices= [ [ A_1 , B_1 ] , [ C_1 , D_1 ] ] )
and then, he can change the blocks by calling

M_V . set ( matrices= [ [ A_2 , B_2 ] , [ C_2 , D_2 ] ] )
Where each two matrices X 1 and X 2 (X ∈ {A, B, C, D}) are defined on the same spaces, and have the same profile. The user can choose which matrices to update, by calling

Grids

Πgasus offers different types of grids. The user can associate to the discrete space, a volume, a surface, an edge or a boundary grid. The volume, surface, edge profiles are related to the elements, while the boundary profile is for the whole patch (sub-domain). Here is an example of how to create a grid for a given discrete space. The user can also create a grid without associating it to any space. This may be interesting if he wants to take advantage of the B-splines and NURBS tranformations, in order to developp, for example, a Finite Difference code. In addition, Πgasus offers the possibility to create, manipulate a metric object and associate it to a grid. Remember that when using a CAD-description to manipulate the metric, we never need to store the whole grid. On each element, we only store the 1D array, the grid will be constructed, on the fly, at each assembling step.

F = lambda r , t : [ r * np . cos ( 2 . * np . pi * t ) , r * np . sin ( 2 . * np . pi * t ) ] DF = lambda r , t : [ np . cos ( 2 . * np . pi * t ) \ , -2 . * np . pi * r * np . sin ( 2 . * np . pi * t ) \ , np . sin ( 2 . * np . pi * t ) \ , 2 . * np . pi * r * np . cos ( 2 . * np . pi * t ) ] M2 = metric ( analytic= [ F , DF ] ) # . . . # . . . # using a l i s t o f p o i n t

Diagnostics

The user can apply different operators to a given field. Once for example, you have the Finite Element description, you can ask Πgasus to evaluate the field or its derivatives. Now let's take the example of Anisotropic Diffusion, which writes :

∂ t u -∇ • (K∇u) = f, Ω (4.11)
u = 0, ∂Ω (4.12)
In general, u denotes the temperature ( inside the plasma ), and K the conductivity. In this example, we consider the evolution of a Gaussian pulse, where:

K = Dµ sin 2 (θ) + µ ⊥ cos 2 (θ) D(µ -µ ⊥ ) sin(θ) cos(θ) D(µ -µ ⊥ ) sin(θ) cos(θ) Dµ cos 2 (θ) + µ ⊥ sin 2 (θ) (4.13)
We would like to visualize the perpendicular diffusion effect.

Visualization using Pylab

Let us start by recalling a very important result for splines curves.

The following result shows how we can approach the values at some knots using the B-splines coefficients. Let us define the average knots as, t j =

t j+1 +•••+t j+k-1 k-1 , for each j ∈ {1, • • • , N + 1}.
|S(t j ) -[S] j | C(k)h 2 D 2 S [tj+1,tj+k-1] (4.14)
where h := max i∈Λ (t i+1 -t i ).

If we define V S = i∈Λ S(t i )N i , then we have

S -V S h 2 , ∀S ∈ S k (T, I) ∩ C m , m 2 (4.15)
As we can notice, the results of this form are important; we do not need to evalute the spline surface for visualization. We will only need to associate at each average knot, the correspondant coefficient.

In the next example, we show how to do a fast-plot of a given field:

1 # D e f a u l t c a l l 2 F_V . fast_plot ( )

The function fast plot has the following arguments:

ai patch id the current patch, useControlPoints is True if the B-slines coefficients will be associated to the control points as a first approximation, savedPoints used if useControlPoints is False. If it is True, the user must provide list P.

list P a list of the associated points to B-slines coefficients.

The user still has the possibility to create his own grid and then evaluate the field.

Parallelization

Πgasus was designed in order to be parallelized, even if the current version is sequential. In the future it will be possible to use generic frameworks as Murge [START_REF] Hénon | PaStiX: A High-Performance Parallel Direct Solver for Sparse Symmetric Definite Systems[END_REF], PetSc [START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF], Hips [START_REF] Gaidamour | HIPS : a parallel hybrid direct/iterative solver based on a Schur complement approach[END_REF] and have a direct link with some linear solvers (Pastix [START_REF] Hénon | PaStiX: A High-Performance Parallel Direct Solver for Sparse Symmetric Definite Systems[END_REF], Mumps [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF][START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF]).

Application to a non-linear equation

In this section, we show as an application of IGA, the resolution of a non-linear partial differential equation. We will present two methods: Picard and Newton's algorithms.

In the sequel, we shall consider the following problem: Find u such that:

-∇ • (A∇u) + Bu = F (x, u) , Ω u = 0 , ∂Ω (5.16)
Let V h be the discrete space, such that V h = span{ϕ b , b ∈ {1, • • • , n}}, then the variational formulation of (5.16) is :

Ω (A∇u) • ∇ϕ b + Ω Buϕ b = Ω F (x, u)ϕ b , ∀b ∈ {1, • • • , n} (5.17)
thus, by expanding u h over V h , using 5.18) this leads to solve the problem :

u h = b∈{1,••• ,n} [u] b ϕ b , we get : b∈{1,••• ,n} [u] b { Ω (A∇ϕ b ) • ∇ϕ b + Ω Bϕ b ϕ b } = Ω F (x, u h )ϕ b , ∀b ∈ {1, • • • , n} (
S[u] = F([u]) (5.19)
where,

S b,b = Ω (A∇ϕ b ) • ∇ϕ b + Ω Bϕ b ϕ b , ∀b, b ∈ {1, • • • , n} (5.20)
F([u]) b = Ω F (x, u h )ϕ b , ∀b ∈ {1, • • • , n} (5.21)

Numerical results : Example from combustion theory

In this section, we shall solve the equation :

-∆u = -ae βu (5.29) This example occurs in combustion theory, but also models the electrostatic potential in a charged body.

The general form of solutions is :

u(x, y) = 1 β ln 8C aβ - 2 β ln |(x + A) 2 + (y + B) 2 -C| (5.30)
for more solutions, we refer to [START_REF] Polyanin | Handbook of Nonlinear Partial Differential Equations[END_REF]. In order to have the function u, vanishing at the boundary, we shall take the following values of parameters:

C = - 1 2 , A = B = 0, aβ = -4 which gives, u(x, y) = - 2 β ln |x 2 + y 2 + 1 2 | (5.31)
One can easily check that u verifies:

-∆u = 4 β e βu (5.32)
In the following test, we took β = -1.

To have homogeneous Dirichlet boundary condition, the domain will be a circle of radius √ 2 2 , centered at 0.

In figure 13, we plot the numerical solution. In figure 14, we plot the evolution of the error with respect to the number of iterations for Picard and Newton algorithms. As expected, we recover the good behavior of the Newton's algorithm.

In the Appendix, we give the detailed codes for Picard's (source code 6) and Newton's methods (source code 7).

Conclusions

In this article, we present a new library Πgasus developed to use the new paradigm of the IsoGeometric Analysis. Πgasus was designed in order to give a comfortable framework for both mathemacians and users. The parallelization of the library is under developpement. The library will also offer local refinements and multiple entry-points for linear solvers. Another active aera is to provide fast solvers and preconditionners, for a specific class of domains and problems, especially in plasmas physics problems and electromagnetism.

For the moment, the library is mainly used at the C.E.A. Cadarache, and shortly at the IPP Garching and FOM Institute for Plasma Physics. People interested by the library, can contact the author. In the future, it will be available for downloads, subjected to a a specific licence. np . savetxt ( " runs/minerror . t x t " , np . array ( list_minerror ) ) p r i n t ( " S o l v i n g done using P i c a r d ' s method" ) p r i n t ( " a f t e r "+str ( i ) +" i t e r a t i o n s , with a t o l e r a n c e o f "+str ( tol ) ) # p r i n t " l i s t m i n e r r o r = " , l i s t m i n e r r o r # p r i n t " l i s t n o r m = " , l i s t n o r m Source Code 6: Non-elliptic equation using Picard's algorithm ) and ( i < niter ) ) : p r i n t ( " S o l v i n g done using Newton ' s method" ) p r i n t ( " a f t e r "+str ( i ) +" i t e r a t i o n s , with a t o l e r a n c e o f "+str ( tol ) ) p r i n t " f i n a l e r r o r i s " , str ( minerror ) Source Code 7: Non-elliptic equation using Newton's algorithm

Figure 1 :

 1 Figure 1: B-splines functions associated to the knot vector T = {000 111}, of order k = 3. These are Bernstein polynomials

Figure 2 :

 2 Figure 2: B-splines functions associated to the knot vector T = {000 1 2 3 44 555}, of order k = 1, 2, 3 (left to right).

5 3 5

 55 11}.In figure 3, we show an example of a B-spline curve. In this case, the B-splines generated are simply Berstein polynomials, and the obtained curve is a Bézier curve.

Figure 3 :

 3 Figure 3: (left) A B-spline curve and its control points, (right) B-splines functions used to draw the curve. N = 9, p = 2 , T = {000, 1 4

Figure 4 :

 4 Figure 4: First line : (left) A quadratic B-spline curve and its control points. The knot vector is T = {000, 1 4 , 1 2 , 3 4 , 111}. (right) The curve after a p-refinement, the degree was raised by 2. Second line : (left) The curve after a hrefinement by inserting the knots {0.1, 0.2, • • • , 0.9} while the degree is kept equal to 2. (right) The curve after duplicating the multiplicity of the internal knots { 1 4 , 1 2 , 3 4 }, this leads to a Bézier description. We can then, split the curve into 4 pieces (sub-domains), each one will corresponds to a quadratic Bézier curve.

Figure 5 : 3 ,

 53 Figure 5: First line : (left) A quadratic NURBS description of a coil and its control points. The knot vector is T = {000, 1 3

Figure 6 :

 6 Figure 6: Mapping from the patch to the physical domain: (left) initial patch, (right) patch after h-refinement in the η direction. Here, we have K = F(Q)

Figure 7 :

 7 Figure 7: Grid generation: (left) The coarsest mesh, (right) Domain after h-refinement. The minimal degree of basis functions is 2 in this example.

# d e f

  i n e t h e r i g h t hand s i d e func_f = lambda x , y : [ 2 . 0 * ( ( 2 * pi ) * * 2 ) * sin ( 2 * pi * x ) * sin ( 2 * pi * y ) ] F_V = field ( space=V , func = func_f ) # t h i s i s t h e unknown func_u = lambda x , y : [ sin ( 2 * pi * x ) * sin ( 2 * pi * y ) ] U_V = field ( space=V , func = func_u )

Ωf

  ([u], x) ϕ b ϕ b dΩ while the Stiffness matrix is Ω (A([u], x) ∇ϕ b ) • ∇ϕ b dΩ In this example, we have f ([u], x) = 1.0 and A([u], x) = Id R 3 is the identity matrix in R 3 . Using Πgasus , this can be defined as # d e f i n e t h e mass matrix func_mass = lambda x , y : [ 1 . 0 ] M_V = matrix ( spaces= [ V , V ] , ai_type=MASS , func=func_mass ) # d e f i n e t h e s t i f f n e s s matrix func_stiff = lambda x , y : [ 1 . 0 , 0 . 0 , 0 . 0 , 1 . 0 ]

# e x pRemark 3 . 1

 31 o r t m a t r i c e s t o scipy-c s r format Mass_V = M_V . to_csr ( ) Stiffness_V = S_V . to_csr ( ) Sigma = Mass_V + Stiffness_V # e x p o r t t h e r i g h t hand s i d e lpr_rhs = F_V . get ( ) # s o l v i n g using s p s o l v e from scipy . sparse . linalg import spsolve lpr_u = spsolve ( Sigma , lpr_rhs ) # s o l v i n g using splu from scipy . sparse . linalg import splu op_Sigma = splu ( Sigma . tocsc ( ) ) # F a c t o r i z a t i o n lpr_u = op_Sigma . solve ( lpr_rhs ) # import t h e r i g h t hand s i d e i n t o P i ga s u s U_V . set ( lpr_u ) The user can also manipulate the objects Fields and matrices using the operators *=, += or .dot(. . . )In (figure8), we plot the numerical solution of Poisson's equation on the unit square.

Figure 8 :

 8 Figure 8: Poisson's equation on the unit square : plot of the numerical solution

  def singleton ( cls ) : instances = {} def getinstance ( ) : i f cls not in instances : instances [ cls ] = cls ( ) r e t u r n instances [ cls ] r e t u r n getinstance

  @singleton c l a s s common_obj ( object ) : def __init__ ( s e l f ) : # . . . # importing t h e F o r t r a n module # . . . import pyfem as py s e l f . pyfem = py . pyfem # . . . # . . . # d e f i n i n g o b j e c t s l i s t s # . . . s e l f . fields =

Figure 9 :

 9 Figure 9: An UML class diagram for Πgasus

2 Figure 10 :

 210 Figure 10: Faces numerotation for boundary conditions. (left) 1D case, (right) 2D case.

Figure 11 :

 11 Figure 11: Faces numerotation for boundary conditions for a multi-patchs domain in the 2D case.

  xml" ls_domain_W = "domain W . xml" # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # D e f i n i t i o n o f t h e space V # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * V = space ( as_file=ls_domain_V ) V . set_boundary_conditions ( ) V . create_grids ( type=" l e g e n d r e " , k=lpi_ordregl ) # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # D e f i n i t i o n o f t h e space W1 # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * W1 = space ( as_file=ls_domain_W1 ) W1 . set_boundary_conditions ( ) W1 . create_grids ( space=V ) # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # D e f i n i t i o n o f t h e space W2 # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * W2 = space ( as_file=ls_domain_W2 ) W2 . set_boundary_conditions ( ) W2 . create_grids ( space=V ) # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # D e f i n i t i o n o f t h e space W as W1 x W2 # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * W = space_vect ( spaces= [ W1 , W2 ] ) W . set_boundary_conditions ( ) W . create_grids ( space=V ) # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

#

  M i s a Block-matrix = [ [ A, B ] , [ C ,D] ] d e f i n e d p r e v i o u s l y by # M = matrix ( m a t r i c e s = [ [A, B ] , [ C ,D ] ] ) # i f f o r example , only A and D must be re-e v a l u a t e d then t h e u s e r must c a l l M . assembly ( [ [ True , F a l s e ] , [ F a l s e , True ] ] )

#

  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # D e f i n i t i o n o f t h e space V # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * V = space ( as_file=ls_domain_V ) V . set_boundary_conditions ( ) V . create_grids ( type=" l e g e n d r e " , k=lpi_ordregl , profile=" boundary " ) # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * By grids, we mean, by default, quadrature ones. The user can use the predefined ones, or gives his own grids, as follows: # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # D e f i n i t i o n o f t h e space V # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * V = space ( as_file=ls_domain_V ) V . set_boundary_conditions ( ) V . create_grids ( list_nodes=list_x , list_weights=list_w ) # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

#

  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # D e f i n i t i o n o f t h e g r i d G # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * geo = cg . cad_geometry ( ls_domain ) G = grids ( profile="volume" , api_k=lpi_ordregl , ao_geometry=geo . list [ : ] , as_type=" l e g e n d r e " ) # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

#

  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # D e f i n i t i o n o f M e t r i c s # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # . . . # using a CAD-d e s c r i p t i o n # . . . M1 = metric ( geometry=geo . list [ : ] ) # . . . # . . . # using an a n a l y t i c a l f u n c t i o n # . . .

  s and d e r i v a t i v e s # . . . M3 = metric ( points=lpr_points ) # . . . # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Then we can associate the metric to a grid, as the following. # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # D e f i n i t i o n o f t h e g r i d G # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * geo = cg . cad_geometry ( ls_domain ) G = grids ( type=" l e g e n d r e " , k=lpi_ordregl , metric=Met ) # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

  def func_b ( x , y ) : lr_tetha = get_tetha ( x , y ) r e t u r n [ cos ( lr_tetha ) , sin ( lr_tetha ) ] U_V . set_func ( func_b ) lpr_values = U_V . eval ( 0 , operator=GRAD_S )In figure12, we plot the numerical solution and the transverse part.

Figure 12 :

 12 Figure 12: Anisotropic Diffusion on the [-1, 1] × [-1, 1] square : plot of the solution (right) and the normal (transverse) part (left)

Theorem 4 . 7

 47 With the above notations, we have, for each S ∈ S k (T, I) ∩ C 1 ,

Figure 13 :

 13 Figure 13: Non-elliptic equation using Picard's algorithm : plot of the numerical solution

Figure 14 :

 14 Figure 14: Non-elliptic equation: evolution of the error during the (left) Picard's and (right) Newton's algorithms.

  # ! /usr/bin/python from test_params import * from scipy . io import mmread , mmwrite from pigasus . fem . constants import * from pigasus . fem . field import * from pigasus . fem . norm import * from pigasus . fem . grids import * from pigasus . fem . matrix import * from pigasus . fem . space import * import pigasus . fem . fem as fem fe = fem . fem ( stdoutput=True , ai_detail=0) # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # D e f i n i t i o n o f t h e space V # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * V = space ( as_file=ls_domain ) V . dirichlet ( faces = [ [ 1 , 2 , 3 , 4 ] ] ) V . set_boundary_conditions ( ) V . create_grids ( type=" l e g e n d r e " , k=lpi_ordregl ) # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # D e f i n i t i o n o f f i e l d s # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * from numpy import log , exp func_u0 = lambda x , y : [-2 . 0 * log ( x * * 2 + y * * 2 + 0 . 5 ) ] U0_V = field ( space=V , func = func_u0 ) func_u = lambda x , y : [-2 . 0 * log ( x * * 2 + y * * 2 + 0 . 5 ) ] U_V = field ( space=V , func = func_u ) # . . . t h i s i s t h e non-l i n e a r p a r t def func_F ( list_F , x , y ) : r e t u r n [ 4 . 0 * exp ( list_F [ 0 ] ) ] F_V = field ( space=V , func = func_F , func_arguments=[U_V ] ) def func_dF ( list_F , x , y ) : r e t u r n [ 4 . 0 * exp ( list_F [ 0 ] ) ] dF_V = field ( space=V , func = func_dF , func_arguments=[U_V ] ) # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # D e f i n i t i o n o f norms # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * N_U = norm ( field=U_V , type=NORM_L2 ) # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # D e f i n i t i o n o f m a t r i c e s # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * func_mass = lambda x , y : [ 1 . 0 ] Ma_V = matrix ( spaces=[V , V ] , ai_type=MASS , func=func_mass ) Mn_V = matrix ( spaces=[V , V ] , ai_type=MASS , func=func_dF , func_arguments=[U_V ] ) func_stiff = lambda x , y : [ 1 . 0 , 0 . 0 , 0 . 0 , 1 . 0 ] S_V = matrix ( spaces=[V , V ] , ai_type=STIFFNESS , func=func_stiff ) # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * fe . initialize ( ) fe . assembly ( matrices=[Ma_V , S_V ] , fields=[U0_V , U_V ] ) Mass_V = Ma_V . to_csr ( ) Stiffness_V = S_V . to_csr ( ) from scipy . sparse . linalg import spsolve # i n i t i a l i z a t i o n # l p r u = s p s o l v e ( Mass V , U0 V . g e t ( ) ) lpr_u = np . zeros ( U_V . size ) # p r i n t " s i z e : " , U V . s i z e U_V . set ( lpr_u ) ; np . savetxt ( " runs/u 0 . t x t " , lpr_u ) list_norm = [ ] fe . assembly ( norms=[N_U ] ) list_norm . append ( N_U . get ( ) ) # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * i = 0 miniter = i minerror = error list_minerror = [ ] # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * while ( ( error > tol

Table 1 :

 1 Modeling conics using NURBS.

Table 2 :

 2 List of the implemented Differential operators

Acknowledgement : We warmly thank G. Latu, V. Grandgirard, G. Dif-Pradalier (CEA-Cadarache), B. Nkonga (University of Nice) and E. Sonnendr ücker (University of Strasbourg and IPP-Garching), for useful discussions and advices. This project was founded by the ANYMOS-ANR.

Picard's algorithm

To solve iteratively 5.19, let us start with the Picard algorithm, which is the simplest one but also the less accurate.

1. X 0 is given, 2. knowing X n , we solve :

Newton's algorithm

Let us define the function :

thus [u] is a zero of the function g. To solve 5.19, we use Newton's method. As J g(X) = S -J F (X) , the Newton's method is:

• knowing X n , we solve :

The algorithm is the following:

1. we compute the mass matrix associated to the function : ∂ u F, i.e : Source Code 2: p-refinement example: Starting from the latter description, we use the following script to elevate the spline degree : +2 in the ξ-direction and +3 in the η-direction.


Nicola Gigli
email: gigli@unice.fr

Michel Ledoux
email: ledoux@math.univ-toulouse.fr

FROM LOG SOBOLEV TO TALAGRAND: A QUICK PROOF

Keywords: 2000 Mathematics Subject Classification. Primary: 51F99; Secondary: 39B62 Metric measure spaces, log-Sobolev inequality, Talagrand inequality

We provide yet another proof of the Otto-Villani theorem from the log Sobolev inequality to the Talagrand transportation cost inequality valid in arbitrary metric measure spaces. The argument relies on the recent development [2] identifying gradient flows in Hilbert space and in Wassertein space, emphasizing one key step as precisely the root of the Otto-Villani theorem. The approach does not require the doubling property or the validity of the local Poincaré inequality.

1. Introduction. The Otto-Villani theorem [START_REF] Otto | Generalization of an inequality by Talagrand and links with the logarithmic Sobolev inequality[END_REF] indicates that, on a Riemannian manifold, the validity of the log Sobolev inequality implies the Talagrand transportation cost inequality (preserving moreover optimal constant). The original proof is inspired by the Riemannian formalism on the space of probability measures established by Otto in [START_REF] Otto | The geometry of dissipative evolution equations: the porous medium equation[END_REF]. An alternate more direct proof relying on infimum convolutions and Hamilton-Jacobi equations was presented later in [START_REF] Bobkov | Hypercontractivity of Hamilton-Jacobi equations[END_REF]. Since then, the relation between the two families of inequalities, and its extension to abstract metric measure spaces, has been investigated by a number of authors, e.g. [START_REF] Lott | Hamilton-Jacobi semigroup on length spaces and applications[END_REF], [START_REF] Gozlan | A characterization of dimension free concentration in terms of transportation inequalities[END_REF] and the recent [START_REF] Gozlan | Characterization of Talagrand's transport-entropy inequalities in metric spaces[END_REF] (see also [START_REF] Villani | Optimal transport, Grundlehren der Mathematischen Wissenschaften[END_REF], [START_REF] Gozlan | Transport inequalities. A survey, Markov Process[END_REF]). In particular, the latter authors address the implication in a rather general metric space context by means of dimension free measure concentration tools [START_REF] Gozlan | A characterization of dimension free concentration in terms of transportation inequalities[END_REF], further investigated in the recent [START_REF] Gozlan | Characterization of Talagrand's transport-entropy inequalities in metric spaces[END_REF] which develops the use of Hamilton-Jacobi equations and the Herbst argument (and covers extended versions for modified log Sobolev inequalities and general transportation costs). In [START_REF] Lott | Hamilton-Jacobi semigroup on length spaces and applications[END_REF], the Otto-Villani theorem is established in metric spaces satisfying doubling property and local Poincaré inequality assumptions. These conditions actually allow the authors to rely on the existence of a weak differentiable structure on the metric measure space (built by Cheeger in [START_REF] Cheeger | Differentiability of Lipschitz functions on metric measure spaces[END_REF]) which in turn gives the possibility to mimic the arguments valid in a Riemannian context.

By the dimension free nature of both the log Sobolev and the Talagrand inequalities, it is quite unnatural to impose either the doubling condition on the measure or the validity of a local Poincaré inequality, as these assumptions are finite dimensional in nature. Recent investigations [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF] (see also [START_REF]Density of lipschitz functions and equivalence of weak gradients in metric measure spaces[END_REF], and [START_REF] Ambrosio | Heat Flow and calculus over spaces with Ricci curvature bounded from below -the compact case[END_REF] for a simplified exposition in the compact case) achieved to develop a differential calculus in metric measure spaces, without the doubling property and the local Poincaré inequality conditions, towards the equivalence of the heat flow generated by a suitable Dirichlet energy and the Wasserstein gradient flow of the relative entropy functional in the space of probability measures. As one key step in this identification process, the heat flow is shown to have a metric speed in the Wasserstein space controlled by the Fisher information, which is precisely the root of the Otto-Villani theorem as emphasized in the formalism of [START_REF] Otto | Generalization of an inequality by Talagrand and links with the logarithmic Sobolev inequality[END_REF]. This key argument has been put forward first in [START_REF] Gigli | Heat flow on Alexandrov spaces[END_REF] in the context of Alexandrov spaces, and its proof is presented below as 6 and 7 in the smooth case. It has been then extended in [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF] to general metric measure spaces.

It is the purpose of this short note to use this observation to provide a simple and direct proof of the Otto-Villani theorem in a general metric measure theoretic context. While the result itself does not necessarily go strictly beyond the framework of [START_REF] Gozlan | A characterization of dimension free concentration in terms of transportation inequalities[END_REF] and [START_REF] Gozlan | Characterization of Talagrand's transport-entropy inequalities in metric spaces[END_REF], the emphasis is perhaps more on methodology, the differentiable structure on metric measure space of [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF] providing the suitable and coherent circle of ideas leading to the Otto-Villani theorem as already expected by the formalism of [START_REF] Otto | Generalization of an inequality by Talagrand and links with the logarithmic Sobolev inequality[END_REF].

We deal below with standard Polish metric measure spaces. Actually, the most natural setting would be that of extended Polish spaces as developed in [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF], because this framework covers the Wiener space, which is the standard infinite dimensional space where log Sobolev inequalities hold. The approach we propose, based on the calculus tools developed in [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF], covers also this general framework: proofs are verbatim the same. We however preferred to state everything in the more familiar context of complete and separable metric spaces just to avoid dealing with the uncommon technology presented in [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF].

The new argument towards the Otto-Villani theorem put forward in this note is already of interest in the smooth case so that we briefly present it below in Section 3. The rest of the paper will be devoted to recall and organize the main conclusions of the recent investigation [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF] linking the gradient flows in Hilbert space and in the Wasserstein space in metric measure spaces relevant to the Otto-Villani theorem. In particular, it is of importance to suitably describe the family of log Sobolev inequalities under consideration. The main result is contained in the last section. We start by collecting standard notions in metric (measure) spaces.

2. Metric notions. This section recalls the necessary metric notions towards the subsequent investigation. We refer to [START_REF] Villani | Optimal transport, Grundlehren der Mathematischen Wissenschaften[END_REF], [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF] for complete details. The metric spaces (X, d) will always be complete and separable. Dealing with metric measure space (X, d, m), the measure m will always be a Borel probability measure. A curve γ : [0, 1] → X is said to be absolutely continuous if there exists g ∈ L 1 (0, 1) such that

d(γ t , γ s ) ≤ s t g(r) dr, ∀t < s ∈ [0, 1]. (1)
In this case, the metric speed of γ is well-defined for a.e. t by

| γt | := lim h→0 d(γ t+h , γ t ) |h| .
It turns out that | γ| ∈ L 1 (0, 1) and this is the minimal L 1 function for which 1 holds.

Given f : X → R, the local Lipschitz constant |∇f | : X → [0, ∞] is defined by |∇f |(x) := lim y→x |f (x) -f (y)| d(x, y) ,
and the ascending/descending slopes

|∇ + f |, |∇ -f | : X → [0, ∞] by |∇ + f |(x) := lim y→x (f (y) -f (x)) + d(x, y) , |∇ -f |(x) := lim y→x (f (y) -f (x)) - d(x, y) ,
where (•) + , (•) -are the positive and negative part, respectively. One says that

G : X → [0, ∞] is an upper gradient of f : X → R provided for any absolutely continuous curve γ : [0, 1] → X it holds f (γ 1 ) -f (γ 0 ) ≤ 1 0 G(γ t )| γt | dt. ( 2
)
It is not difficult to check (see for instance Proposition 2.5 of [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF] for the simple proof) that if f is Lipschitz, then |∇f |, |∇ + f | and |∇ -f | are all upper gradients. Denote by P(X) the space of Borel probability measures on X and by P 2 (X) ⊂ P(X) the space of probability measures with finite second moment. The Wasserstein distance W 2 (µ, ν) between two measures µ, ν ∈ P(X) is defined by

W 2 2 (µ, ν) := inf γ X×X d 2 (x, y) dγ(x, y)
where the infimum is taken among all plans γ ∈ P(X × X) such that π 1 γ = µ and π 2 γ = ν. Notice that we are defining W 2 between arbitrary measures in P(X), so that it is possible that the value +∞ is attained. The Wasserstein distance W 2 is a metric on P 2 (X), and (P 2 (X), W 2 ) is again complete and separable.

The relative entropy functional Ent m : P(X) → [0, ∞] with respect to a given m ∈ P(X) is defined by

Ent m (µ) :=    X ρ log ρ dm, if µ = ρ m, +∞, otherwise.
Recall that sublevels of the entropy are tight. Indeed, using first the bound z log(z) ≥ -1 e and then Jensen's inequality we get that for any µ = ρ m such that Ent m (µ) ≤ C and for any Borel set E, it holds

1 e + C ≥ m(X \ E) e + Ent m (µ) ≥ E ρ log ρ dm ≥ µ(E) log µ(E) m(E)
so that the claim follows from the tightness of m. We will denote by D(Ent m ) ⊂ P(X) the domain of the entropy, i.e. the set of measures µ ∈ P(X) such that Ent m (µ) < ∞. Under these notations, the Talagrand transportation cost inequality (proved in [START_REF] Talagrand | Transportation cost for Gaussian and other product measures[END_REF] in Euclidean space for m the Gaussian measure) states that for some constant C > 0 and every µ ∈ P(X),

W 2 2 (µ, m) ≤ C Ent m (µ). (3)
The Otto-Villani theorem expresses that this inequality holds under a log Sobolev inequality (to be recalled below), preserving the constants.

3. The smooth case. In this short section, we present the argument of [START_REF] Gigli | Heat flow on Alexandrov spaces[END_REF] and the principle of proof of the Otto-Villani theorem which we would like to emphasize in the case of a smooth Riemannian manifold (X, g) with the standard differential calculus (as presented e.g. in [START_REF] Villani | Optimal transport, Grundlehren der Mathematischen Wissenschaften[END_REF]). The log Sobolev inequality for a probability measure m on the Borel sets of X indicates that there is a constant α > 0 such that for every smooth function

f : X → [0, ∞) with X f dm = 1, 2α X f log f dm ≤ {f >0} |∇f | 2 f dm. ( 4
)
The task is thus to establish, under the log Sobolev inequality 4, the transportation cost inequality 3 (with constant C = 2 α ; for the standard Gaussian measure in Euclidean space, α = 1 is optimal for both the log Sobolev and the Talagrand inequalities).

Given µ ∈ P(X) with (smooth) density f with respect to m, denote by (f t ) the heat flow with respect to the Dirichlet form X ∇u, ∇v dm starting from f . Set ϕ(t) = Ent m (µ t ), t ≥ 0, where

µ t = f t m. The standard calculus (Boltzmann formula) yields d dt Ent m (µ t ) = - {ft>0} |∇f t | 2 f t dm, (5)
so that the log Sobolev inequality 4 expresses equivalently that ϕ ≤ -2αϕ.

Following the strategy of [START_REF] Otto | Generalization of an inequality by Talagrand and links with the logarithmic Sobolev inequality[END_REF] through the formal Riemannian structure on (P 2 (X), W 2 ), the aim is to show that the curve t → µ t = f t m in (P 2 (X), W 2 ) has a metric speed satisfying

| μt | 2 ≤ {ft>0} |∇f t | 2 f t dm (a.e. t). (6)
Once this holds, the conclusion easily follows. Namely, by 5, the latter amounts to | μt | ≤ -ϕ (t). Since by the log Sobolev inequality ϕ ≤ -2αϕ and thus

-ϕ ≤ - 2 α ( √ ϕ) ,
it follows by integration (according to 1 in (P 2 (X), W 2 )) that, for all T > 0,

W 2 (µ, µ T ) ≤ T 0 | μt | dt ≤ - 2 α T 0 ϕ(t) dt = 2 α ϕ(0) -ϕ(T ) . Since Ent m (µ T ) → 0 as T → +∞, W 2 (µ, m) ≤ lim T →+∞ W 2 (µ, µ T ) ≤ lim T →∞ 2 α ϕ(0)-ϕ(T ) = 2 α ϕ(0) = 2 α Ent m (µ),
which amounts to the transportation cost inequality 3 with constant C = 2 α . The central point of the approach is thus to establish 6. The following simple and direct proof has been thus put forward in [START_REF] Gigli | Heat flow on Alexandrov spaces[END_REF] and shown in [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF] to be flexible enough to extend to the case of general metric measure spaces. The argument relies on a suitable interlacing of the heat flow and the Hamilton-Jacobi equation.

For a continuous Lipschitz and bounded function g on X, and t, s > 0, write

X (Q 1 g)f t+s dm - X gf t dm = 1 0 d dr X (Q r g)f t+rs dm dr where Q r g(x) = inf y∈X g(y) - d(x, y) 2 2r
is the infimum convolution of g with the quadratic cost. The semigroup (Q r ) r>0 represents the Hopf-Lax solution of the Hamilton-Jacobi equation

d dr Q r f (x) + |∇Q r f | 2 (x) 2 = 0
(for almost every t, x). Hence,

X (Q 1 g)f t+s dm- X gf t dm = 1 0 X - |∇Q r g| 2 2 f t+rs -s ∇Q r g, ∇f t+rs dm dr.
By the Cauchy-Schwarz inequality

-∇Q r g, ∇f t+rs ≤ 1 2s |∇Q r g| 2 f t+rs + s 2 |∇f t+rs | 2 f t+rs so that X (Q 1 g)f t+s dm - X gf t dm ≤ s 2 2 1 0 X |∇f t+rs | 2 f t+rs dm dr.
By the dual Kantorovich representation of the Wasserstein distance W 2 (cf. e.g. [START_REF] Villani | Optimal transport, Grundlehren der Mathematischen Wissenschaften[END_REF]), taking the supremum on g on the left-hand side of the preceding inequality yields, for all t, s > 0,

W 2 2 (µ t+s , µ t ) ≤ s 2 1 0 X |∇f t+rs | 2 f t+rs dm dr, (7)
from which 6 follows. The Otto-Villani theorem is thus established in this way.

4. The general metric measure space case. To address the preceding proof of the Otto-Villani theorem in the context of a general metric space, we first recall the necessary differential calculus in non-smooth spaces developed in the recent [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF] (and to which we refer for further details) which provides the key step 6. We repeat that metric measure spaces (X, d, m) studied in this paper consist of a complete and separable metric space (X, d) equipped with a Borel probability measure m.

4.1.

The Sobolev space W 1,2 (X, d, m). There are several equivalent definitions of the Sobolev space W 1,2 (X, d, m), here we follow the approach -and use the results -of [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF]. Define the functional Ch :

L 2 (X, m) → [0, ∞] by Ch(f ) := inf lim n→∞ 1 2 X |∇f n | 2 dm,
where the infimum is taken among all sequences (f n ) of Lipschitz functions converging to f in L 2 (X, m). Then the space W 1,2 (X, d, m) is defined as the set of f 's such that Ch(f ) < ∞ endowed with the norm d,m) there exists a map |∇f | * ∈ L 2 (X, m) -which plays the role of the modulus of distributional gradient for Sobolev functions on R d -which is characterized by the following two properties:

f 2 W 1,2 := f 2 L 2 + 2Ch(f ). It turns out that for f ∈ W 1,2 (X,
Ch(f ) = 1 2 X |∇f | 2 * dm, |∇f | * ≤ G, m -a.e.
∀G such that there exists a sequence (f n ) of Lipschitz functions converging to f in L 2 such that

|∇f n | G in L 2 (X, m).
Two theorems are worth to mention.

Theorem 4.1 (Relation with upper gradients). Let f ∈ L 2 (X, m) and G an upper gradient for it as defined in 2. Assume that G ∈ L 2 (X, m).

Then f ∈ W 1,2 (X, d, m) and |∇f | * ≤ G, m -a.e. ( 8
)
Theorem 4.2. For any f ∈ W 1,2 (X, d, m), there exists a sequence of Lipschitz functions

(f n ) converging to f in L 2 (X, m) such that |∇f n | converges to |∇f | * in L 2 (X, m).
For the proof of the first theorem, see Theorem 6.2 of [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF] (or [START_REF] Ambrosio | Heat Flow and calculus over spaces with Ricci curvature bounded from below -the compact case[END_REF] for the compact case). The second is obvious.

Let us briefly explain the importance of these results. The typical definition of W 1,2 (X, d, m), see for instance [START_REF] Cheeger | Differentiability of Lipschitz functions on metric measure spaces[END_REF], is to proceed by relaxing upper gradients of L 2 functions, rather than the local Lipschitz constant of Lipschitz functions. With this approach, Theorem 4.2 is totally non trivial, while Theorem 4.1 is straightforward (we refer to the survey [START_REF] Heinonen | Nonsmooth calculus[END_REF] for an overview of the standard approach to Sobolev spaces over metric measure spaces). A step forward made in [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF] has been to show that to relax the local Lipschitz constant produces the same space given by the relaxation of upper gradients, which therefore yields to both Theorem 4.1 and Theorem 4.2, 4.2. Properties of the gradient flow of Ch. Notice that the functional Ch : L 2 (X, m) → [0, ∞] is convex, lower semicontinuous and with dense domain, so that the standard theory of gradient flows on Hilbert space applies. This means that for any f ∈ L 2 (X, m) there exists a unique locally absolutely continuous curve

(f t ) ⊂ L 2 (X, m) on [0, ∞) such that f 0 = f and it holds d dt f t ∈ -∂Ch(f t ), a.e. t.
In [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF] (with an approach inspired by [START_REF] Gigli | Heat flow on Alexandrov spaces[END_REF]) properties of this gradient flow have been studied extensively, especially in connection with concepts coming from optimal transport theory like Wasserstein distance and relative entropy. In particular, the following main results of [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF] are critical in the identification of the gradient flows for the Cheeger energy in L 2 (X, m) and for the entropy functional in (P 2 (X), W 2 ). As announced and illustrated in Section 3, they actually contain the key step towards the Otto-Villani theorem (described here in item (ii) corresponding to 6). Theorem 4.3. Let (X, d, m) be a metric measure space, µ ∈ D(Ent m ) and let f be its density with respect to m. Let (f t ) be the gradient flow of Ch starting from f . Then the following are true:

(i) f t m ∈ P(X) for any t ≥ 0 ;

(ii) the curve t → µ t := f t m is locally absolutely continuous on (0, ∞) with respect to the Wasserstein distance W 2 and for its metric speed | μt | it holds

| μt | 2 ≤ {ft>0} |∇f t | 2 * f t dm, a.e. t;
(iii) the map t → Ent m (µ t ) = X f t log f t dm is locally absolutely continuous and it holds

d dt Ent m (µ t ) = - {ft>0} |∇f t | 2 * f t dm, a.e. t.
Items (i) and (iii) are proved in (c), (d) of Theorem 4.16 in [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF] while (ii) is the content of Lemma 6.1 therein. For the simplified approach in the case of compact spaces, we refer to [START_REF] Ambrosio | Heat Flow and calculus over spaces with Ricci curvature bounded from below -the compact case[END_REF].

5. The Otto-Villani theorem. On the basis of the previous differential calculus and the main Theorem 4.3, we address in this section the Otto-Villani theorem in the general context of a metric measure space with the scheme of proof put forward in the smooth case (Section 3).

We start by noticing that various typical formulations of the log Sobolev inequality are all equivalent. The result of the following proposition is already well-known if m is doubling and (X, d, m) supports a local Poincaré inequality (see [START_REF] Lott | Hamilton-Jacobi semigroup on length spaces and applications[END_REF]). However, as remarked in the introduction of [START_REF]Density of lipschitz functions and equivalence of weak gradients in metric measure spaces[END_REF], there is no need for these hypotheses, and the statement below is just an exemplification of what is expressed in [START_REF]Density of lipschitz functions and equivalence of weak gradients in metric measure spaces[END_REF] for the particular case of the log Sobolev inequality. Technically speaking, the proof is based on the density of Lipschitz functions in W 1,2 (X, d, m) proved in [START_REF]Calculus and heat flows in metric measure spaces with Ricci curvature bounded from below[END_REF].

Proposition 1 (Equivalent formulations of the log Sobolev inequality). Let (X, d, m) be a metric measure space and let α > 0. Then the following are equivalent:

(i) For any Lipschitz function f :

X → [0, ∞) with X f dm = 1, it holds 2α X f log f dm ≤ {f >0} |∇f | 2 f dm.
(ii) For any Lipschitz function f :

X → [0, ∞) with X f dm = 1, it holds 2α X f log f dm ≤ {f >0} |∇ -f | 2 f dm.
(iii) For any continuous f : X → [0, ∞) with X f dm = 1 and any upper gradient

G of f it holds 2α X f log f dm ≤ {f >0} G 2 f dm. (iv) For any non-negative f ∈ W 1,2 (X, d, m) with X f dm = 1, it holds 2α X f log f dm ≤ {f >0} |∇f | 2 * f dm.
Proof. (ii) ⇒ (i) is trivial, and given that for Lipschitz functions slopes are upper gradients (iii) ⇒ (ii) is also trivial. The property of weak gradients expressed in inequality 8 shows that also (iv) ⇒ (iii) is true. Thus everything boils down in proving that (i) implies (iv). For this we use the density Theorem 4.2. Pick f ∈ W 1,2 (X, d, m) and assume without loss of

generality that {f >0} |∇f | 2 * f dm < ∞. For c > 0, let f c := a c min{c -1 , max{f, c}}, where a c is such that X f c dm = 1. Notice that a c → 1 as c ↓ 0. Also, observe that X |∇fc| 2 * fc = 4Ch( √ f c ) < ∞ and use Theorem 4.2 to find a sequence (g n c ) of Lipschitz functions such that g n c → √ f c and |∇g n c | → |∇ √ f c | * in L 2 (X, m).
Up to passing to a subsequence we can also assume that the convergence of (g n c ) to √ f c is also m-a.e.. Then the Lipschitz functions f n c := (g n c ) 2 converge to f c pointwise a.e. and thus by dominated convergence it holds

X f c log f c dm = lim n→∞ X f n c log f n c dm.
Since f n c is Lipschitz for any n, c by hypothesis we know that 2α

X f n c log f n c dm ≤ {f n c >0} |∇f n c | 2 f n c dm = 4 X |∇g n c | 2 dm.
The right-hand side of this equality converges, by construction, to 4Ch(

√ f c ) = X |∇fc| 2 * fc dm. Hence we proved that 2α X f c log f c dm ≤ X |∇f c | 2 * f c dm,
for any c > 0. Letting c ↓ 0 we complete the proof.

As the proof shows it would be equivalent to state (ii) with the ascending slope

|∇ + f | in place of |∇ -f |.
The formulation here is the way the log Sobolev inequality is usually written (see e.g. Theorem 3.22 of [START_REF] Villani | Optimal transport, Grundlehren der Mathematischen Wissenschaften[END_REF]). Also, we point out that the equivalence between (i) and (iv) relies only on Theorem 4.2 (and not on the delicate Theorem 4.1), which in turn is a direct consequence of the way |∇f | * is defined.

We then say that (X, d, m) supports the log Sobolev inequality with constant α provided any of the inequalities in the Proposition 1 is true. Lemma 5.1. Assume that (X, d, m) supports the log Sobolev inequality with constant α > 0 and let µ ∈ D(Ent m ). Let f be the density of µ with respect to m and (f t ) be the gradient flow of Ch starting from f . Set µ t = f t m, t ≥ 0. Then as t → +∞, we have Ent m (µ t ) → 0 and µ t m in duality with C b (X).

Proof. Let ϕ(t) := Ent m (µ t ), t ≥ 0, and recall that (iii) of Theorem 4.3 yields that ϕ is locally absolutely continuous with ϕ (t) = -X |∇ft| 2 * ft dm. Then the fact that (X, d, m) supports the log Sobolev inequality with constant α implies (see formulation (iv) of Proposition 1) that ϕ ≤ -2αϕ, which by the Gronwall lemma further gives ϕ(t) ≤ e -2αt ϕ(0) = e -2αt Ent m (µ), [START_REF] Gozlan | Characterization of Talagrand's transport-entropy inequalities in metric spaces[END_REF] so that the first claim is proved. To prove the second one, recall that the sublevels of Ent m are tight, and so for any sequence t n → ∞ there exists a subsequence, not relabeled, such that (µ tn ) weakly converges to some ν as n → ∞ in duality with C b (X). By 9 we then have Ent m (ν) ≤ lim n→∞ Ent m (µ tn ) = 0, hence ν = m. Since this is independent on the chosen sequence (t n ), the claim is proved.

Provided with these results, the Otto-Villani theorem is now established exactly as in the smooth case (Section 3). Theorem 5.2 (Otto-Villani theorem). Let (X, d, m) be a metric measure space supporting the log Sobolev inequality with constant α > 0 (Proposition 1). Then it also supports the Talagrand inequality with constant 2 α , i.e. it holds

W 2 2 (µ, m) ≤ 2 α Ent m (µ)
for all µ ∈ P(X).

The theorem immediately leads to the following consequence.

Corollary 1. Let (X, d, m) be a metric measure space supporting the log Sobolev inequality for some α > 0. Then m has finite second moment.

Proof. Thanks to the Talagrand inequality, it is sufficient to prove that D(Ent m ) ∩ P 2 (X) is non empty. But this is obvious, as any measure of the form m(A) -1 m | A , with A ⊂ X Borel, bounded and with positive m-measure is in both D(Ent m ) and P 2 (X).