Hervé Guillard
email: herve.guillard@inria.fr

Jalal Lakhlili
email: jalal.lakhlili@ipp.mpg.de

Adrien Loseille
email: adrien.loseille@inria.fr

Alexis Loyer
email: alexis.loyer@inria.fr

Boniface Nkonga
email: ¶boniface.nkonga@inria.fr

Ahmed Ratnani ||

Ali Elarif

Tokamesh : A software for mesh generation in Tokamaks

Keywords: Tokamak fusion reactor, Automatic mesh generation, unstructured meshes, blockstructured meshes, spline tensor product Tokamak, Maillage automatique, maillage non-structurés, maillage bloc-structuré, produit tensoriel de fonctions splines

The generation flux aligned mesh for tokamaks is a difficult task that is not easily automated and very often ask for manual intervention and specific expertise. This report presents the work done in the framework workpackage of the WP 5.1 of the EoCoE action to design and build a software for this task.

Part I Introduction 0.1 Generalities

Small scale turbulence as well as large scale MHD instabilities in Tokamak plasmas are believed to be quasi bi-dimensional : the fluctuations scales in the parallel direction to the magnetic field are typically 3 to 4 orders of magnitude larger than those in the transverse plane. This property is used in many tokamak plasma simulation codes to have a different spatial resolution in the parallel and perpendicular direction. The necessity to consider a sufficiently fine spatial resolution in the transverse direction (that can be of the order of the Larmor radius) make particularly attractive the use of grids/meshes which are adapted to the magnetic flux surfaces in each 2D poloidal plane while the third periodic toroidal direction can be discretized with great accuracy with Fourier series. The shape of the poloidal cross-section of the averaged flux surfaces in the poloidal planes is governed by magnetic equilibrium. These magnetic equilibria are themselves computed by dedicated equilibrium codes. Having as input a prescribed magnetic field equilibrium, the purpose of the WP5.1 task of the EoCoE project is to develop a generic software to construct 2D triangular/quadrangular meshes of the poloidal plane of a tokamak, which will be aligned to magnetic flux-surfaces.

Design of the software and contents of the report

Tokamesh has been designed based on minimal assumptions about the user's input. To construct flux surface aligned meshes of the poloidal plane, the minimal inputs that the meshing software must require are the boundaries of the computational domain and the definition of a flux function inside this computational domain. Although, there is a growing interest for tokamak simulation to include the walls into the computational domain, the domains considered in tokamesh are restricted to the vacuum vessel of the Tokamaks or to a part of it. Specifically, edge plasma codes (SOLPS,FBGKI, TOKAM3X, SOLEDGE) as well as GYSELA exclude a central part of the core plasma around the magnetic axis (O-point). The boundary of this central part will be considered as a constant flux surface. We note also that in some machines, (for instance MAST) poloidal coils are located inside the vacuum vessel. In in this case, for the tokamesh software, the computational domain must exclude these parts of the machine. Thus instead of being defined by the walls of the vacuum vessel, the boundaries of the computational domain is rather defined by a (part of) constant flux surface. The definition of these constant flux surface boundaries cannot be automated as it depends on the user's interest. There are also a priori no reasons for the computational grid where the flux function is computed to contain the boundary of the vacuum vessel or the in-vessel flux surfaces that the user wishes to be a boundary of the computational domain. Therefore in tokamesh, the definition of the vacuum vessel is independent of the definition of the grid where the flux function is defined. Given these two inputs : a mesh where the magnetic flux function is defined and a definition of the vacuum vessel boundary, the tokamesh software will follow the following steps :

• Pre-processing of the data • Domain decomposition in logical squares • Meshing and gluing of the different subdomains.

Pre-processing of the data Pre-processing the data have been found necessary in the first stage of this work. This was not foreseen at the beginning of the project but the interactions with potential users has showed that the inputs that they were providing were not appropriate to generate flux aligned grids. Specifically, the equilibrium codes like EFIT [START_REF] Fitzgerald | Efit tokamak equilibria with toroidal flow and anisotropic pressure using the two-temperature guiding-centre plasma[END_REF] or cedres++ [START_REF] Heumann | Quasi-static free-boundary equilibrium of toroidal plasma with cedres++[END_REF] use low-order approximations and coarse grids to solve the Grad-Shafranov equation with the results that the inputs (mesh and flux solution) lack a sufficient quality and resolution to generate "good" flux aligned meshes. Therefore, we have developed in tokamesh, some tools to improve the quality of the data generated by these equilibrium solvers. This is described in section 2.1. In addition, the software contains an independent free equilibrium Grad-Shafranov solver working on unstructured meshes that allows to get a C 1 equilibrium solution provided that the user can give the domain boundary conditions and the functional form with respect to the magnetic flux of the pressure and toroidal current (see section 2.2).

Domain decomposition

Once, the data is judged to be of sufficient quality, a fundamental task in tokamesh consists to split the computational domain into physical subdomains that can be mapped onto computational square grids. This is a common task with existing tokamak mesh generator e.g [START_REF] Marchand | Carre: a quasi-orthogonal mesh generator for 2d edge plasma modelling[END_REF]. However, the approach used in tokamesh extends these previous approaches and does not rely on the physical interpretation of the subdomain decomposition. Instead we use the Morse theory for C 2 functions [START_REF] Morse | The Calculus of Variations in the Large[END_REF] [START_REF] John | Morse Theory[END_REF] as well as the companion concept of Reeb graph [START_REF] Reeb | Sur les points singuliers d'une forme de pfaff complètement intégrable ou d'une fonction numérique[END_REF] to show rigorously that the computational domain enclosed by a level set of the magnetic flux can be decomposed into subdomains that are either homeomorphic to a disk or a ring. The notion of Reeb graph allows to organize the mappings of these different subdomains and to construct a global mesh of the collection of these subdomains. These constructions are given in section 4.2.

Meshing and gluing of the different subdomains Up to now, the tasks in tokamesh are independent of the discretization techniques that are used in the simulation codes. From this points however, one must distinguish between the numerical techniques that use unstructured meshes and those that rely on block-structured meshes. We recall that an unstructured mesh is a tessellation of a domain of the plane by simple polygonal shapes in an irregular pattern. In practice, the polygons are very often triangles or quadrangles or a mixture of the two. On the opposite, block-structured meshes relies on the partition of the computational domain into a small number of sub-domains where each sub-domain is discretized by a structured grid. Unstructured meshes can easily handle very complex geometrical shapes but must maintain a list storing neighboring relations between the mesh elements. Parallelism are vectorization are less efficient with these type of meshes than with block-structured ones. On the opposite, the discretization of complex geometrical shape is extremely difficult with block-structured meshes. In tokamaks, since the grids must be flux aligned, if the vacuum vessel boundaries do not coincide with a flux surface, it is not possible to create block-structured meshes that are both aligned with the flux surfaces and that respect the vacuum chamber boundaries. Therefore, following a common practice in this domain, the computational domain is restricted to the interior of a flux surface close to the vacuum chamber boundary.

Unstructured meshes

Once the domain decomposition has been generated, for the simulation codes using unstructured triangulation, the strategy used in tokamesh consist to generate first a cloud of node lying on the level sets of the magnetic flux function in each subdomain. The generation of these sets of points is done using two different techniques, one of them gives Inria orthogonal grids but will cluster the grid points in regions of high curvature, in contrast the other technique will give a uniform distribution of the nodes along the isolines but the mesh will not be orthogonal. These two techniques are described in section 3.2.1. In practice, the actual distribution of the nodes can be a weighted average of the two methods. Once the cloud of node has been generated, a second step consists to eliminate all the nodes that are not inside the vacuum vessel and the final step creates the mesh using a constrained and anisotropic Delaunay technique as described in section 3.2.2 Block-structured meshes For the construction of block-structured meshes, the first steps of the strategy is almost identical to the previous case except that instead of generating in each subdomain a cloud of nodes, one generates a mapping between the physical domain and a logical square [0, 1] × [0, 1]. The construction of these mappings uses approximation of the level sets of the flux function by cubic b-splines and construction of the mapping by spline tensor product see section 4.4.1 and 4.4.2. The construction of these mapping uses heavily the fact that each subdomain is homeomorphic to a disk or annulus. Once done, these mapping have to be combined in order to have a global mesh of the computational domain. A global mapping involving only C 0 continuity only requires that the curves defining the boundaries of the subdomains are identical and this imposes some compatibility conditions between the discretization of the boundaries of the subdomain that are described in section 4.3. However, in addition to the previous requirement, the construction of a global mapping with C 1 continuity imposes also that the derivatives of the mapping are continuous when passing through a subdomain boundary. This is described in section 4.4.3.

The remaining of this report is as follows : In part I, we give the necessary theoretical concepts used in the algorithms developed in the software. This part is split in 4 sections developing successively the general tools used through all the algorithms, the pre-processing of the data provided by the users, the algorithms for the generation of unstructured grids and finally the ones developed for the creation of block-structured meshes. Part II then presents the installation of the software and how to use it. For the potential users of tokamesh, this last part have to be complemented by the html documentation present in the sphinx folder of the documentation and presenting the individual functions.

Part II

Theoretical part 1 General tools 1.1 Clough-Tocher Interpolation

The construction of a flux aligned mesh in tokamesh requires the knowledge of the flux function denoted ψ(x) in the sequel. For all modern tokamak configurations, this function cannot be specified analytically and is known as the numerical solution of the Grad-Shafranov (GS) equation [START_REF] Shafranov | Plasma equilibrium in a magnetic field[END_REF]. Solving this equation requires to define a mesh covering the domain of interest as well as a numerical approximation procedure. Although the methods to solve the GS equations are numerous, many of them, specially when coupled to real-time or near real-time identification methods use low -order numerical methods. In tokamesh, we have decided to use the Cedre++ software [START_REF] Heumann | Quasi-static free-boundary equilibrium of toroidal plasma with cedres++[END_REF]. This choice is guided by the fact that Cedre++ using a P1 triangular finite element method is very general and is not limited by the geometry in contrast to other GS solvers relying on more accurate methods. Tokamesh also offers an interface with the output of the EFIT software [START_REF] Fitzgerald | Efit tokamak equilibria with toroidal flow and anisotropic pressure using the two-temperature guiding-centre plasma[END_REF] that appears to be very often used in the fusion community. However, the experiments that we have done on several input provided by potential users of tokamesh show that the output produced by these software very often uses a too coarse spatial resolution to produce high quality and smooth outputs. As the construction of flux aligned meshes depends on the smoothness of the iso-contours computed by the GS equilibrium codes, we had to develop some generic tools to improve the quality of the fields provided by the GS solvers. For the sake of generality, since it is always possible to convert a mesh into an equivalent triangular one, these tools have to be effective on unstructured triangular meshes. In a first version of the software, we have used for this task Powell-Sabin finite element but the present version use instead the reduced Clough Tocher or (Heish Clough Tocher) interpolation [START_REF] Clough | Finite element stiffness matrices for analysis of plates in bending[END_REF] as this family of finite element is simpler to handle than Powell-Sabin and seems in practice to have the same accuracy. The following section is a brief description of Clough Tocher interpolation. More detailed and accurate description can be found in [START_REF] Boisserie | Sur l'implementation des éléments finis de hseieh-clough-tocher complet et reduit[END_REF], [START_REF] Alfeld | Smooth macro-elements based on clough-tocher triangle splits[END_REF] and for a general presentation of splines defined on triangulation in [START_REF] Lai | Spline Functions on Triangulations[END_REF].

Definition of Heish Clough Tocher (HCT) finite element

The HCT method is a finite element method that uses as basis function C 1 functions. In the sequel, we use the same notation as in [START_REF] Boisserie | Sur l'implementation des éléments finis de hseieh-clough-tocher complet et reduit[END_REF].

A reduced HCT finite element is the triple (K, P K , Σ K ) where the data K, P K and Σ K are defined as follows :

• The set K is a triangle with vertices a 1 , a 2 , a 3 and a the barycenter of K. Denoting by K i (cf. Figure 1 ) the triangles with vertices a, a i+1 , a i+2 , i ∈ {1, 2, 3[≡ 3]}.

• The space P K is given by :

P K = P ∈ C 1 ; P |Ki ∈ P 3 (K i ), i = 1, 2, 3 and ∂P ∂n ∈ P 1 (K i ) ∀ K ⊂ K
where ∂P ∂n is the outer normal derivative of the function P |K i and K i the side opposite to the vertex a i .

Σ K = {P (a i ), DP (a i )(a i+1 -a i ), DP (a i )(a i+2 -a i ), i = 1, 2, 3[≡ 3]} where DP (a i )(a i+1 -a i ) is the directional derivatives across the side [a i a i+1 ].
Note that in contrast to other finite element methods, the definition of the finite element space can be done on the triangles in the physical space without the use of a reference element and a mapping [START_REF] Boisserie | Sur l'implementation des éléments finis de hseieh-clough-tocher complet et reduit[END_REF]. This is one of the advantage of this method.

Passage from local to global degrees of freedom

Let v a fonction of class C 1 on the domain Ω. Denoting by DL(v) and DG(v) respectively the vector of local and global degrees of freedom, which are defined by :

DL(v) = [v(a 1 ), v(a 2 ), v(a 3 ), Dv(a 1 )(a 3 -a 1 ), Dv(a 1 )(a 2 -a 1 ), Dv(a 2 )(a 1 -a 2 ), Dv(a 2 )(a 3 -a 2 ), Dv(a 3 )(a 2 -a 3 ), Dv(a 3 )(a 1 -a 3 )] DG(v) = [v(a 1 ), v(a 2 ), v(a 3 ), ∂ x v(a 1 ), ∂ y v(a 1 ), ∂ x v(a 2 ), ∂ y v(a 2 ), ∂ x v(a 3 ), ∂ y v(a 3 )]
The directional derivatives and the global derivatives are related by :

Dv(a i )(a i-1 -a i ) Dv(a i )(a i+1 -a i ) = x i-1 -x i y i-1 -y i x i+1 -x i y i+1 -y i ∂ x v(a i ) ∂ y v(a i )
So that, DL(v) and DG(v) are linked by :

DL(v) = DG(v)D RR n°9230
where D is a block diagonal matrix of size 9 × 9 given by :

D =         1 0 1 1 d 1 d 2 0 d 3        
where we denote d i the matrix defined by :

d i = x i-1 -x i y i-1 -y i x i+1 -x i y i+1 -y i , fori ∈ {1, 2, 3 mod 3}

Expression of basis functions

The basis functions value can be evaluated with the use of the barycentric coordinates of the points in K.

Denoting by [R i ], the vector of 9 basis functions associated with the subtriangle K i (i ∈ {1, 2, 3 mod 3}). From [START_REF] Boisserie | Sur l'implementation des éléments finis de hseieh-clough-tocher complet et reduit[END_REF], we have that :

[R i ] = [A i ][λ] (1)
where [λ] is the vector of Bernstein polynomials given by :

λ = λ 3 i λ 3 i+1 λ 3 i+2 λ 2 i λ i+2 λ 2 i λ i+1 λ 2 i+1 λ i λ 2 i+1 λ i+2 λ 2 i+2 λ i+1 λ 2 i+2 λ i λ i λ i+1 λ i+2
and [A i ] is a 9 × 10 matrix defined according to the excentricity parameters [START_REF] Boisserie | Sur l'implementation des éléments finis de hseieh-clough-tocher complet et reduit[END_REF].

The derivatives of these basis functions are therefore defined by :

[∂ j R i ] = [A i ][∂ j λ] j = 1, 2, 3 [∂ 2 jk R i ] = [A i ][∂ 2 jk λ] j, k = 1, 2 (2)
with

∂ j R = ∂R ∂λj and ∂ 2 jk R = ∂ 2 R ∂λj ∂λ k 1.1.

HCT interpolation of function and derivatives of interpolate

Let (K, P K , Σ K ) a given reduced HCT finite element. Given a function v defined over the set K, the HCT-interpolant Π K (v) of the function v is defined through the relations :

Π K (v) ∈ P K and R i (Π Ki (v)) = R i (v)
with R i the basis functions vector associated to K i . From this definition, the interpolation of a function v ∈ C 1 (Ω) on the subtriangle

K i of K ⊂ Ω is equivalent to : Π Ki (v) = [DL(v)][R i ]
with [R i ] and DL(v) defined in [START_REF] Alfeld | Smooth macro-elements based on clough-tocher triangle splits[END_REF].

Inria

Then, from (2) and this equivalence, the derivatives of the interpolating polynomials Π Ki (v) are given by the follows expressions :

               ∂ x Π Ki v = DL(v)[A i ][∂ x λ] ∂ y Π Ki v = DL(v)[A i ][∂ y λ] ∂ x xΠ Ki v = DL(v)[A i ][∂ x xλ] ∂ x yΠ Ki v = DL(v)[A i ][∂ x yλ] ∂ y yΠ Ki v = DL(v)[A i ][∂ y yλ] (3)
with DL, [A i ] and λ defined previously.

Regular family of reduced HCT Let K a triangle of reference and B a compact in K. Let (K) a family of reduced HCT. We denote by : h = diamK and ρ = diam(registred cercle in K)

The family (K) is regular if :

1. There is a constant α > 0 such that α ≤ ρ h for all element of family (K). 2. Let F denote the affine transformation of the plane that verifies

F (ā i ) = a i , 1 ≤ i ≤ 3. Then ā = F -1 (a) is in the compact B.
with āi and a i respectively the vertices of K and K and ā, a some points of K and K respectively.

Definition 1.1. Interpolation error

Given a regular family of reduced HCT triangles, there exists a constant C such that

|v -Π K (v)| m,K ≤ Ch 3-m K |v| 3,K ∀v ∈ H 3 (K) and m = 0, 1, 2 (4)
where,

|v| m,K =   K |α|=m |∂ α v| 2   1 2
and h K is the diameter of K Proof : The proof of this result is given in [START_REF] Ciarlet | Interpolation error estimates for the reduced hsieh-clough-tocher triangle[END_REF] The figure 2 illustrates this result. On this figure, the function v(x, y) = x 4 + y 4 is interpolated on meshes successively refined and the L 2 discrete norm of the interpolation error |v -Π K (v)| versus the mesh parameter is plotted. On can check on this figure that in agreement with (4) the error is of order 3 on the function and 2 on its derivatives.

Stream-line integration

Constructing a structured grid on a domain Ω is equivalent to identify a mapping

(s, t) ∈ [0, 1] × [0, 1] → x(s, t) ∈ Ω (5)
Since we want to construct flux aligned grids, the isolines of the flux function correspond to one of the coordinate lines, say t = for instance. Therefore the streamlines of the ode where v(s) is any scalar velocity will cross the isolines. In tokamesh, this result is used many times to identify seed points to generate isolines. More precisely, in a connected subdomain corresponding to an edge of the Reeb graph (see section 4.2.2 for these notions), the isolines are closed and unique. One way, to generate a particular isoline ψ(x) = ψ 0 is thus to identify a seed point x 0 where ψ(x 0 ) = ψ 0 and then to follow the isoline ψ 0 by computing the intersection of this isoline with the elements. If the mesh where the function ψ(x) is defined, is a P 1 triangular mesh, this procedure is not sensible to round-off error and the generated isoline is closed and will come back exactly on the seed point x 0 . To identify the seed point, we will use streamline integration of the ode [START_REF] Czarny | Bézier surfaces and finite elements for mhd simulations[END_REF]. More precisely, we used the following result that allow to control accurately the length of the integration step :

   dx ds = v(s)∇ψ x (0) = x 0 (6) RR n°9230
Lemma : Let x 0 ∈ Ω, we set ψ 0 = ψ(x 0 ) then on the curve S (s) defined by the solution of the ordinary differential system :

   dx ds = ψ M -ψ 0 ∇ψ 2 ∇ψ x (0) = x 0 (7)
the function Ψ(s) = ψ(x(s)) is linear and

Ψ(s) = ψ M s + (1 -s)ψ 0 . Proof : We have dψ(x(s))/ds = (ψ M -ψ 0 )∇ψ.dx(s)/ds = ψ M -ψ 0 = const, thus ψ(s) = ψ 0 + (ψ M -ψ 0 )s.

Inria

This result therefore allows to find easily a point where a specific isoline passes through : Choose a random point in a connected component of Ω and then integrate the ode (7) until s = 1. The solution x M = x(s = 1) will verify ψ(x M ) = ψ M . This procedure is illustrated Figure 3: Streamline integration from an isoline on figure 1.2 that shows the integration from several points on the same isoline (green curve) toward another isoline (green dashed curve). The blue curves represent an spline interpolation of the curves generated through the integration of (6). All the curves end on the same isoline. Although any ode solver can used to integrate (6), we have found that this task requires accurate solver. In practice, the evaluation of the gradient ∇ψ is done with the Clough-Tocher C1 interpolation described in section 1.1 and we use dopri5, [START_REF] Hairer | Solving Ordinary Differential Equations I. Nonstiff Problems[END_REF], an explicit Runge-Kutta method of order (4)5 due to Dormand and Prince that use step size control to monitor the accuracy and is available in the scipy library [START_REF] Hairer | Solving Ordinary Differential Equations I. Nonstiff Problems[END_REF]. This procedure is quite costly and can be simplified for many cases. However, it has been found necessary for some difficult cases and it is the default procedure used in tokamesh 1 .

Note also that the result [START_REF] Deboor | A practical guide to splines[END_REF] allows for a simple construction of a one-to-one mapping from [0, 1] × [0, 1] on Ω such that the curves t = constant are isolines of ψ : Let t ∈ [0, 1] and let C 0 (t) be a parametric representation of the isoline ψ = ψ 0 , then x(s, t) is defined as the solution at s ∈ [0, 1] of the ode 7 with initial condition x (0) = C 0 (t). This mapping is obviously one-toone since being streamlines of the gradient of ψ, for s = s , the curves S (s) and S (s ) cannot cross. Moreover for any s 0 ∈ [0, 1], the curve x(s 0 , t) is the isoline ψ 1 s 0 + (1 -s 0 )ψ 0 . Such a mapping is orthogonal since isolines and gradient streamlines are orthogonal by definition. While orthogonality is usually a desirable feature for a grid, it can also create extreme clustering of the grid nodes in the regions where the isolines are convex and possess a high curvature. Thus this procedure cannot be used alone but must be mixed with the methods described in section (4.4.1) if the mesh extends in the plasma core. However, the use of this mapping can be an option for edge simulation codes.

Spline approximation 1.3.1 Generalities on splines

This section recalls some basic properties of B-splines curves and surfaces. We also recall some fundamental algorithms (knot insertion and degree elevation). For a basic introduction to the subject, we refer to the book [START_REF] Tiller | The NURBS Book[END_REF]. A B-Splines family, (N i ) 1 i n of order k, can be generated using a non-decreasing sequence of knots T = (t i ) 1 i n+k .

Definition 1.1 (B-Splines 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) = χ [tj ,tj+1[ (x)
for k ≥ 1 and 1 ≤ j ≤ n.

We note some important properties of a B-splines basis:

• B-splines are piecewise polynomial of degree p = k -1,

• 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,

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

The vectorial space spanned by these B-splines, which we denote S k (T, I), where I denotes the interval [t 1 , t n+1 ], is called the Schoenberg space.

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), i.e. open knot vector,

• 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 the i th control point P i affects C(t), only in the interval [t i , t i+k ],

• The control polygon approaches the behavior of the curve.

Remark 1.3. In order to model a singular curve, we can use multiple knots or control points, i.e. P i = P i+1 .

Generation of the B-spline functions N k j

In order to generate the B-spline functions N k j , we need to give the appropriate knot vector:

T = {t 1 = t 2 = • • • = t k = 0 < • • • < t n+1 = t n+2 = • • • = t n+k = 1}
Classical choices are:

• Equally spaced knots:

t i = i-k n+1-k , ∀j ∈ {k + 1, . . . , n},
• Average knots:

t i = i-1 i-k tj p
, ∀j ∈ {k + 1, . . . , n}.

Choice of the parameterization tj

Usually, three methods are used to find the values of tj :

• Equally spaced parametrization: tj = j-1 n-1 , ∀j ∈ {1, . . . , n}.

• Chordal parametrization: tj = tj-1 +

Pj -Pj-1 d
, ∀j ∈ {3, . . . , n}.

Where we defined:

d = n j=2
P j -P j-1 , and t0 = 0

• Centripetal parametrization: tj = tj-1 + √ Pj -Pj-1 d
, ∀j ∈ {3, . . . , n}.

Where we defined:

d = n j=2
P j -P j-1 , and t0 = 0

In the case of tokamaks, for the core region where isolines are closed, an alternate parametrization that have been found useful is to parametrize the data point P j by the angular value of the segment OP j where O denotes the magnetic axis :

t j = angle(OP 0 OP j ) 2π
where P 0 is an arbitrary point on the curve.

Fundamental geometric operations

By inserting new knots into the knot vector, we add new control points without changing the shape of the B-Spline curve. This can be done using the DeBoor algorithm [START_REF] Deboor | A practical guide to splines[END_REF]. We can also increase the degree of the B-Spline family and keep unchanged the curve [START_REF] Huang | Fast degree elevation and knot insertion for b-spline curves[END_REF]. In (Fig. 5), we apply these algorithms to a quadratic B-Spline curve and we show the position of the new control points.

Deriving a B-spline curve

The derivative of a B-spline curve is obtained as:

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 (8) where Q i = p Pi+1-Pi
ti+1+p-ti+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. More generally, by introducing the B-splines family {N k-j i *

, 1 ≤ i ≤ n -j} generated by the knot vector T j * obtained from T by removing the first and the last knot j times, we have the following result: Proposition 1.4. The j th derivative of the curve C is given by

C (j) (t) = n-j i=1 N k-j i * (t)P (j) i
where, for j > 0,

P (j) i = p -j + 1 t i+p+1 -t i+j P (j-1)
i+1 -P (j-1) i and P (0) i

= P i .
By denoting C and C the first and second derivative of the B-spline curve C, it is easy to show that: Proposition 1.5. We have,

• C (0) = p tp+2 (P 2 -P 1 ),

• C (1) = p 1-tn (P n -P n-1 ),

• C (0) = p(p-1) tp+2

1 tp+2 P 1 -{ 1 tp+2 + 1 tp+3 }P 2 + 1 tp+3 P 3 , • C (1) = p(p-1) 1-tn 1 1-tn P n -{ 1 1-tn + 1 1-tn-1 }P n-1 + 1 1-tn-1 P n-2 .
Example: Let us consider the quadratic B-spline curve associated to the knot vector T = {000 2 5 3 5 111} and the control points {P i , 1 ≤ i ≤ 5}:

C(t) = 5 i=1 N 3 i (t)P i
we have,

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 5 3 5 11}.

Intersection of two spline curves

The problem of finding the intersections between two parametric curves is fundamental in Geometric Modeling. In this work, we need to compute the intersection of the wall with a given magnetic surface. Among the available algorithms, we present in the sequel, the one developped by Mørken et al. [START_REF] Mørken | Computing intersections of planar spline curves using knot insertion[END_REF], which is based on knot insertion. The underlying idea is based on the fact that we can approach the B-splines coefficients, using the value at some knots. Let us consider the Greville abscissa, which are successive averages of knots, as

t j = t j+1 + • • • + t j+k-1 k -1 , ∀j ∈ {1, • • • , n + 1}. (9)
We have the following result Proposition 1.6.

For each S = n j=1 s j N k j ∈ S k (T, I) ∩ C 1 we have, |S(t j ) -s j | C(k)h 2 D 2 S [tj+1,t j+k-1 ] (10)
where

h := max 1≤i≤n (t i+1 -t i ). Moreoer, if we define V S = n i=1 S(t i )N i , we have S -V S h 2 , ∀S ∈ S k (T, I) ∩ C m , m 2 (11)

Inria

The simplest version of the Intersection Algorithm may be stated as the following

• Refine the two B-splines curves, for a given mesh resolution,

• Compute the intersection of their control polygons,

• Use the (Eq. 10) to compute knots approximation for all intersections of the two curves,

In (Fig. 6), we show the intersection of two quadratic B-splines curves. By inserting more knots, the intersection of the two control polygons converge to the desired points. Once the intersections points, and their corresponding knots are computed, we can split the two curves with respect to these knots as shown in (Fig. 6(d)).

(a) (b)

Figure 6: (a) Two quadratic B-spline curves and their control polygons. The knot vector for the blue curve is {000, 1 2 , 3 4 , 3 4 , 111} while for the red curve it is {000, 1 2 , 111}. (b) Curves and their control polygons after inserting the knots {0.1, 0.2, 0.3, 0.4, 0.8} in the blue curve and {0.4, 0.7, 0.9} in the red curve. Intersections of the two control polygons approximate the the intersection of the two curves.

Multivariate tensor product splines

Let us consider d knot vectors

T = {T 1 , T 2 , • • • , T d }.
For simplicity, we consider that these 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, S ki (T i , I). The tensor product of all these spaces is also a Schoenberg space, namely S k (T ),

where k = {k 1 , • • • , k d }.
The cube 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 k1 i1 ⊗ N k2 i2 ⊗ • • • ⊗ N k d i d where, i = {i 1 , • • • , i d }.
A typical cell from P is a cube of the form :

Q i = [ξ i1 , ξ i1+1 ] ⊗ • • • ⊗ [ξ i d , ξ i d +1 ].
We are looking for a B-spline curve of order k (with a degree p = k -1), such that it interpolates the points P i :

n j=1 N k j ( ti )c j = P i , ∀i ∈ {1, . . . , n} (12)
As noticed, in order to solve this problem, we need to know how to generate the B-spline functions N k j and the parameters ti , then we will have to solve the following linear system:

M C = P (13)
M i,j = N k j ( ti ), C i = c i , and

P i = P i (14)
Remark 1.7. Here, we have considered the case of interpolation without end conditions. Additional conditions can be added on the two extremeties of the curve.

Approximation of iso-lines by BSplines

Generally in CAD/CAM systems, one has to reconstruct smooth surfaces from discrete data. Let us consider a set of n points V i , i ∈ {1, • • • , n}. In order to have a well-posed problem, it is necessary to have an additional information about the topology of the points set. This can be done by defining the set of neighbors for each entry point. The approximation problem writes:

Find a smooth parameterized surface F such that minimize i=1,n dist(F(u i ) -V i ) ( 15
)
where dist is a distance function to be defined, as well as the n parameters u i . Many distance functions can be used. In this work, we are interested in least square distance. Moreover, we consider the parametric distance, which means that the parameter u i has been defined previously. Another approach, is to consider the geometry distance which is defined as inf u F(u) -V i .

Approximation by a B-spline curve

The approximation problem can be formulated as the following: Given n points (V k ) 1≤k≤n in the plane and a parametrization

0 = u 1 ≤ u 2 ≤ • • • ≤ u k ≤ u n = 1, k ∈ {1, . . . , n}, find a B-spline curve C(t) = n j=1 N k j (t)P j of order k (with a degree p = k-1), such that
• it approximates the points V k , in the least squares sense:

F(c) = n-1 k=2 C(u k ) -V k 2 → min (16)

Inria

• it interpolates the start and end points (i.e V 1 = C(0) and V m = C(1)).

We can easily show that the minimum is achieved if and only if

D T DP δ = D T V δ , δ ∈ {x, y, z} (17)
D i,j = N k i ( tj ) (18)
Where V δ contains the δ ∈ {x, y, z} coordinate of the points V k , ( tj , 1 ≤ j ≤ n) are the Greville knots, while the matrix D will be defined in 2D in the next section. More details can be found in [START_REF] Tiller | The NURBS Book[END_REF]. It is interesting to see the impact of the parametrization on the constructed curve (c.f. Fig 7).

As we can see, this may lead to high variations of the first derivative. In [START_REF] Pottmann | A concept for parametric surface fitting which avoids the parametrization problem[END_REF], the authors show how to avoid the parametrization problem but the construction is much more complicated.

Approximation by a B-spline surface

Let (V k ) 1≤k≤n be a sequence of points in the plane. In this section, we will construct a mapping F, as a B-spline surface, that maps the unit square (i.e. patch) onto an approximation of the given set of points.

The pure approximation problem can be sketched as follows:

• For each point (V k ) 1≤k≤n we associate a parameter (site) u k . This means also that we have endowed the patch with a specific topology. Hence the mapping F will be such that

F(u k ) V k .
• Given 2 knot vectors, which will define the B-splines family that we will use later.

Construct the mapping F by minimizing the following functional:

J 0 [F] = n k=1 F(u k ) -V k 2 (19)
RR n°9230

Here again, the unknowns are the B-splines coefficients. In general, the constructed mapping given by (Eq. 19), may be not regular. One needs to add a constraint on the curvature. In this case, the general approximation problem writes:

minimize J 0 [F] + µ J 1 [F] (20)
and is called approximation with fairing. The choice of the fairness functional J 1 must be a compromise between low numerical complexity and high quality of the resulting surfaces [START_REF] Greiner | Interpolating and approximating scattered 3D-data with hierarchical tensor product B-splines[END_REF]. This can be adjusted using the constant µ which measures the fairing weight with respect to fitting. High quality can be achieved by using better approximation of the curvature functionals which can be given by the exact thin plate energy. For low numerical complexity, the simple thin plate energy can be used

J 1 [F] = P F ξξ 2 + 2 F ξη 2 + F ηη 2 dP (21)
A good compromise between the simple and exact thin energies is the data dependent thin plate energy [START_REF] Greiner | Data dependent thin plate energy and its use in interactive surface modeling[END_REF]. More details can be found in [START_REF] Greiner | Data dependent thin plate energy and its use in interactive surface modeling[END_REF][START_REF] Greiner | Interpolating and approximating scattered 3D-data with hierarchical tensor product B-splines[END_REF].

Example

Let us consider for the moment the approximation problem using one patch. In order to simplify the text, we will only consider one coordinate (x). We will also use c ij instead of P x ij . The patch coordinates will be denoted by ξ = (ξ, η). A parameter u will the two components

(u 1,k , u 2,k ). Let us define, F[F] = n k=1 F(u k ) -V k 2 + µJ 1 [F] (22)
with,

F(ξ) = n1 i=1 n2 j=1 N (1)
i (ξ)N

(2)

j (η)c ij (23)
and J 1 the simple thin plate energy.

Remark 1.8. In the general case, we shall write c ij instead of c ij . We will consider the vector representation of the matrix c. The corresponding 1D index for (i b , j b ) will be denoted by b.

Simple computations lead to the following expression of J 1 [F]:

J 1 [c] = c T Ec (24)
where the matrix E can be written as

E = A + 2B + C, (25)
with the matrices A, B and C are defined as

A b,b = 1 0 N i b N i b 1 0 N j b N j b (26) B b,b = 1 0 N i b N i b 1 0 N j b N j b (27) C b,b = 1 0 N i b N i b 1 0 N j b N j b (28)

Inria

The minimum is achieved if,

(D T D + µE)c = D T x (29)
where the matrix D is defined as

D b,k = N (1) i b (u 1,k )N (2) j b (u 2,k ) (30)
the matrix D T D can be easily computed

M b,b = D T D b,b = n k=1 N (1) i b (u 1,k )N (2) j b (u 2,k )N (1) i b (u 1,k )N (2) j b (u 2,k )
which can be seen as a mass matrix with weights equal to one (but not evaluated on quadrature points).

Global regularity

For our applications, we need to use multi-patch description in order to take into account the topology of the magnetic field. Hence, we will have to stick these patches together ensuring a C 0 or C 1 regularity. This can be done easily with B-splines surfaces (c.f. subsection 1.3.3).

As an example, let us consider the case of sticking two patches, namely P m and P s , on the faces f m and f s respectively, ensuring the C 1 continuity condition. This can be written as follows:

c :,fm = c :,fs (31)
and,

α m (c :,fm -c :,δfm ) = β s (c :,fs -c :,δfs ) (32)
where the coefficients α m and β s are computed using the knot vectors. The first condition enforces the C 0 continuity be requiring that the curves separating the two patches are identical while [START_REF] Si | Tetgen, a delaunay-based quality tetrahedral mesh generator[END_REF] ensures the C 1 continuity. These conditions can be enforced explicitly or can be added when solving the problem (Eq. 29).

Spline approximation and singular points

In section 4.4, we will construct C 1 meshes by gluing together spline tensor-product patches. A minimal requirement is therefore that the spline curves are well-behaved and do not contain loops or cusps. Moreover, in the process of gluing together two patches, we will impose that the end derivatives of the curves are colinear. We therefore have to take care that this process do not introduce singular points on the spline curve. Thus, in this section, we consider a cubic Bezier curve defined by the control points x k , k = 0, • • • , 3 and the parametric expression :

C k (s) = 3 k=0 x k N k (s) (33)
the basis functions N k (s) are displayed in figure 1.3.7. The expression of these functions and of their derivatives are given below :

N 3 1 (s) = (1 -s) 3 (N 3 1 ) (s) = -3(1 -s) 2 (N 3 1 ) (s) = 6(1 -s) N 3 2 (s) = 3s(1 -s) 2 (N 3 2 ) (s) = 3(1 -s)(1 -3s) (N 3 2 ) (s) = -6(2 -3s) N 3 3 (s) = 3s 2 (1 -s) (N 3 3 ) (s) = 3s(2 -3s) (N 3 3 ) (s) = 6(1 -3s) N 3 4 (s) = s 3 (N 3 3 ) (s) = 3s 2 (N 3 3 ) (s) = 6s (34)
Instead of the 4 control points, x k , the curve [START_REF] Sovinec | Nonlinear magnetohydrodynamics with high-order finite elements[END_REF] is usually defined by their end points x 0 and x 3 and its derivatives d 0 and d 1 . Since we are interested by G -1 continuity, the norm of the end derivatives is not relevant and we can assume that d 0 and d 1 are unit vectors. The control points x 1 and x 2 are then given by 2 :

x 1 = x 0 + αd 0 /3 x 2 = x 3 -αd 1 /3 ( 35
)
where α is a positive parameter whose expression has to be found. The curve then has the following expression :

C k (s) = x 0 (1 -s) 3 + 3 (α d 0 + x 0 ) s (1 -s) 2 + 3 (-α d 1 + x 3 ) s 2 (1 -s) + x 3 s 3 (36)
Introducing the notation u = 1 -s it will be useful to write it as :

C k (s) = x 0 + (s 3 + 3us 2 )(x 3 -x 0 ) + αu 2 sd 0 -αus 2 d 1 ( 37
)
Its derivative is given by :

C k (s) = α[(u 2 -2us)d 0 + (s 2 -2us)d 1 ] + 6su(x 3 -x 0 ) (38)
while its second derivative is

C k (s) = α[(2s -4u)d 0 + (4s -2u)d 1 ] + 6(u -s)(x 3 -x 0 ) (39)
Inria We will examine under what conditions the curve (33) has cusps or loops. Figure 1.3.7 shows a Bezier curve defined by 4 control point with a loop. The red and blue points are the end points of the curve while the green and yellow control the end derivatives. We want to give sufficient conditions on the position of these points in order to avoid loops or cusps on the curve. The curve C k (s) is given by 4 vectors among which only two of them are independent. Our discussion on the existence of singular points for the curve will thus depends on which of them are independent.

The case d 0 = d 1 = d We first examine the case where the two end tangents are identical. In this case we have :

C k (s) = x 0 + (s 3 + 3us 2 )(x 3 -x 0 ) + α(u 2 s -us 2 )d (40) C k (s) = -6 s 2 + 6 s (x 3 -x 0 ) + α 6 s 2 -6 s + 1 d (41) and C k (s) = 6(1 -2s)((x 3 -x 0 ) -αd) (42)
Let us first study when the curve can have a cusp where C k (s) = 0 or a loop. We have two cases: (x 3 -x 0 ) and d are not colinear A loop will be possible if there exits two distinct reals a, b ∈]0, 1[ such that C k (a) = C k (b). Since (x 3 -x 0 ) and d are independent this implies that we must have :

a 2 (3 -2a) = b 2 (3 -2b) a(1 -a)(1 -2a) = b(1 -b)(1 -2b) ( 43
)
and this is not possible since in [0, 1], the function x 2 (3 -2x) is monotone. For the same reason the existence of a cusp where C k (s) = 0 is not possible. (x 3 -x 0 ) and d are colinear

In this case, we note λ = |x 3 -x 0 | and a necessary and sufficient condition for

C k (s) = 0 is that the equation f (s) = 6 s 2 -6 s + 1 α -6 λ s 2 + 6 λ s = 0 (44)
has a solution for s ∈]0, 1[. But equation (42) that writes now :

C k (s) = 6(1 -2s)(λ -α)d (45)
shows that s = 1/2 is always an extremum for the quadratic curve f (s). On s = 1/2, the value of f (s) on the extremum is

f (1/2) = 3λ -α 2 Since in any case f (0) = f (1) = α, we can conclude that f (s) is always in [min(α, f (1/2)), max(α, f (1/2)]. Thus, if α ≤ λ, f (s) > α and is strictly positive while for λ ≤ α < 3λ, f (s) > f (1/2) = 3λ-α 2
and is again always positive. Therefore for α < 3λ, the derivative is always strictly positive and the curve cannot have a cusp. In the same way, it cannot have a loop since Rolle's theorem will imply in this case the existence of a solution in ]0, 1[ of the equation C k (s) = 0. For α > 3λ, the curve admits a cusp.

We summarize these results : if the two end derivatives are colinear, d 0 = d 1 for α ≤ 3|x 3 -x 0 |, the curve C k (s) has no singular point. d 0 and d 1 are independent Thus we can develop (x 3 -x 0 ) on the basis formed by these two vectors :

(x 3 -x 0 ) = λd 0 + µd 1 (46)
where

λ = (x 3 -x 0 ) × d 1 d 0 × d 1 µ = (x 3 -x 0 ) × d 0 d 1 × d 0 (47)
and the expression of the derivative of the curve becomes :

C k (s) = [6suλ + α(u 2 -2us)]d 0 + [6suµ + α(s 2 -2us)]d 1 (48)
The existence of a cusp or a loop 3 therefore means that it exists

(s, u) in ]0, 1[×]0, 1[ such that 6suλ + α(u 2 -2us) = 0 6suµ + α(s 2 -2us) = 0 (49)
3 for a cusp, we must have u = 1 -s while a loop can correspond to different abscissa Inria i.e : 6sλ + α(u -2s) = 0 6uµ + α(s -2u) = 0 (50) whose solutions are

s = α/λ 3(α/λ -2) u = α/µ 3(α/µ -2) (51) Figure 1.3.7 shows the function f (x) = x 3(x -2)
. It can be seen that the region of the plane If µ > 0 and α < 3µ then the curve C k (s) has no loop or singular points. To summarize the results of this section : Let λ > 0 or µ > 0 be the coordinates of (x 3 -x 0 ) in the basis d 0 = d 1 if these vectors are independent, then if α < 3 min(µ, λ, |x 3 -x 0 |, the curve (36) has no cusps or loops.

Pre-processing of the data

The construction of flux aligned meshes depends on the smoothness of the iso-contours computed by the equilibrium codes. In particular, the inputs (mesh and flux solution) must have a sufficient quality and resolution to generate "good" flux aligned meshes. The problem is particularly important near the X point since in many cases, the iso-contours of the flux function in this region are not very regular. Figure 11 shows the level set passing through the x-point as computed by Cedre++. It can be seen that the curve oscillates in this region. This characteristic Figure 11: Level set passing through the x-point for GS equilibrium computed by Cedre++ of the solution is not specific to the Cedre++ software. Figure 12 shows the same curve for an equilibrium computed in the jet tokamak with the EFIT software where the same problem can be noticed.

Data regularization

This absence of smoothness of the equilibrium computed by the equilibrium solvers leads us to develop some tools to regularize the data produced by these codes. The first one is a local refinement tool that selectively allows a mesh refinement in the vicinity of the x-points. Note that a global refinement tool is also available in tokamesh. However, it has been found that an aggressive refinement dividing the spatial step size by factor of order 2 4 or 2 5 is necessary to yield smooth contours near the x-point and a global refinement by these factors is too costly in practice and produce too large meshes (of the order of several M points). Figure 13 shows an example of this refinement process for a double disconnect null equilibrium in the West tokamak.

However, this refinement process does not change by itself the data, we have therefore complemented it by a smoothing using cubic interpolation on this refined mesh. More precisely, given the initial mesh where the equilibrium is computed and the flux function defined on this Inria Figure 13: Local refinement near x-points mesh, we replace this function by its Clough Tocher interpolant (see section 1.1) defined by the value of the function on the node and derivatives at the node defined as an average value of the gradient on the triangles surrounding this node. This CT-interpolant is then evaluated on the refined triangulation and use to produce smooth iso-contours. Although there is no theory on this technique, it has been found quite effective to regularize the data in such a way that it can be used for the computation of smooth iso-contours. Figure 14 shows the result of this technique applied to the data of figure 12. Note that in general this technique produces a slight displacement of the x-points. This is usually, barely noticeable except in the case of a connected double null (CDN) where in agreement to the result given [START_REF] Ma | Structural classification and stability of divergence-free vector field[END_REF] (see section 4.2), these slight displacement of the x-points produce two disconnected x-points lying on two distinct but extremely close isolines. Thus for this specific case, we had to implement a function to connect again the x-points.

Re-computing a Grad-Shafranov equilibrium

One of the drawback of the previous technique is that it gives a regularization of a GS equilibrium computed with a low-order method. Therefore even if the results are smooth, their accuracy is not greater than the accuracy of the original data. For this reason, it can be interesting to provide in tokamesh a GS solver allowing the use of high-order C 1 method to compute the solution of the equilibrium problem expecting that the use of high order C 1 method will provide solutions that will be both smoother and of greater accuracy than the original input. Note however that it is not guaranteed that the use of high-order elements alone can yield an improved regularity as remarked in [START_REF] Heumann | Quasi-static free-boundary equilibrium of toroidal plasma with cedres++[END_REF] : "We are solving here a nonlinear elliptic problem with discontinuous coefficients(in the case of iron-transformer tokamaks) and discontinuous righthand side. The standard convergence theory for finite elements and elliptic regularity theory does not yield improved approximation results for polynomials of degree higher than 1".

Nevertheless we have developed a free boundary Grad-Shafranov solver using Clough-Tocher reduced cubic finite element (see section 1.1). The interest of this finite element family is that it can use the same mesh that the one where the initial data has been computed. Moreover, since the Grad-Shafranov equation is a non-linear one that must be solved with an iterative procedure, the initial data provides a convenient initial guess for this procedure. Nonhomogeneous Dirichlet boundary conditions coming from the equilibrium data computed on the initial mesh are used and the non-linear system is solved by Picard iterations. Note that to take into account non-homogeneous Dirichlet boundary conditions, it has been necessary to develop a penalization method of Nischte type. Figure 15 compare the solution obtained for the computation an equilibrium in the JET tokamak with three different numerical methods : The original equilibrium computed with EFIT, a solution computed with a P1 finite element solver and the one computed with the C 1 Clough-Tocher finite element method. This Grad-Shafranov Figure 15: Comparison of magnetic equilibrium computed by EFIT, a P1 method and a C 1 method using Clough-Tocher finite elements solver can only be used if the user possesses information about the non-linear functions defining the right-hand side of the GS equation, namely the functional form of the pressure and toroidal current with respect to the magnetic flux. In case where these information are missing, the only possibility will be to use the data smoothing procedure described in the previous section 2.1.

Unstructured meshes

Introduction

The use of unstructured meshes is the method of choice when one has to deal with complex geometries. In the tokamak fusion community, the number of codes using this technology is at present limited. A preeminent example is EIRENE [START_REF] Reiter | The eirene and b2-eirene codes[END_REF] devoted to the interactions of the plasma and the wall that require an accurate description of the wall geometry, description that can only be addressed easily with non-structured meshes. Another examples for MHD studies are the Inria the Nimrod [START_REF] Sovinec | Nonlinear magnetohydrodynamics with high-order finite elements[END_REF] and M3DC1 [START_REF] Jardin | A triangular finite element with first-derivative continuity applied to fusion mhd applications[END_REF] codes that uses triangular high-order finite elements. However, the use of triangular unstructured meshes is actually growing since a realistic description of the interactions between the plasma and the vacuum chamber walls becomes mandatory. Some new edge simulation codes are using this approach [START_REF] Giorgiani | A new high-order fluid solver for tokamak edge plasma transport simulations based on a magnetic-field independent discretization[END_REF]. In consequence mesh generation software dedicated to tokamak simulation are currently being developed [START_REF] Zhang | Mesh generation for confined fusion plasma simulation[END_REF]. The previous examples, however does not seem to fully exploit the use of adaptive meshes that can easily be generated with the modern simplicial mesh generation techniques. We describe in the following the techniques and algorithms that have been developed in Tokamesh for this purpose.

Triangulation

A triangular mesh consists of a set of nodes and a set of simplicial elements (triangle in 2D, tetrahedron in 3D) whose edges connect the nodes. To define a triangular mesh, it is then necessary to define these two sets.

Choice of the nodes

There are at least two ways to define the set of nodes required to built a triangular mesh : One can start from a small set of elements, and introduce the nodes in an incremental way, one by one maintaining at each step of the procedure, a triangulation of the set of nodes. In this way, the decision to add a node and where to add it, is governed by the overall quality of the triangulation at the previous step. On the other hand, one can also generate once for all, a cloud of unconnected nodes and then built the triangulation of this set of nodes. In tokamesh, we have chosen the second way and we first generate a cloud of nodes. The reason for this choice is that we want to have a flux aligned mesh where the elements have at least one edge on an isoline of the magnetic flux function. It is therefore natural to create a set of nodes that consists of parallel collections of nodes sampling the iso-lines. In addition since orthogonality is a desirable property, it can be interesting to sample the iso-lines to obtain an "orthogonal" mesh where the triangles are close to right-angled triangles. However, in case of a strong curvature of the isolines, orthogonal meshes produces a non uniform distribution of the nodes along the iso-lines with very small steps near the magnetic axis. The most satisfactory strategy than we have found therefore consists to use the domain segmentation that will be described in section 4.2. Then in the regions corresponding to leaves of the Reeb graph (see section 4.2 or in more physical terms in the plasma core and private flux domain, a set of node based on spline approximation of level sets is generated. Figure 3.2.1 shows this distribution of nodes for an equilibrium in jet where 31 isovalues have been generated in this region. The red lines indicates the P1 isovalues computed on the background mesh while the green dots show the distribution of nodes on the spline approximation of these isovalues. On the opposite, in the region corresponding to edges of the Reeb graph (the SOL region), we generate the set of nodes based on streamline integration of the ode

   dx ds = v(s) ∇ψ ∇ψ x (0) = x i 0 (52)
where v(s) is a scalar velocity and where x i 0 are a set of nodes lying on the separatrix. In this way, we generate a mesh that will be locally orthogonal near the separatrix while being characterized by a uniform distribution near the magnetic axis. In the SOL region, this generation of nodes, is done independently of the existence of the vacuum chamber walls, the nodes exterior to the vacuum chamber are then eliminated in a second Figure 16: Distribution of nodes in the core region for a jet equilibrium step. From a practical point, it is much simpler to proceed in this way than to stop the integration when the vacuum chamber walls are crossed by a streamline as the algorithm to decide if a 2D point is inside or outside a closed curve relies only on simple operations. The triangulation algorithm will then use the nodes defining the vacuum chamber and the interior nodes generated in the SOL to construct a triangulation of the vacuum chamber. Figure 3.2.1 illustrates the procedure for an equilibrium in JET where the integration of equation 52 has been stopped at seven different values of the flux function.

Triangulation of a set of nodes

An overview of the modern techniques for non-structured mesh generation can be found in [START_REF] Loseille | Chapter 10 -unstructured mesh generation and adaptation[END_REF]. We give below the principal concepts used in this field and applied in the algorithms used in tokamesh. The construction of an adapted triangular mesh relies on the notion of the "quality" of an element. In 2d, the "best" possible element is the equilateral triangle where the three edges of the triangle are of equal length. An optimal triangulation is therefore a triangulation where all triangles are equilateral. As detailed for instance in [START_REF] Frey | Mesh Generation, 2nd Edition[END_REF], one way to adapt this notion of optimal triangulation to adaptive meshing is to define a metric M (x) represented as a symmetric positive definite matrix allowing to define the distance between two points a and b as :

d(a, b) = (b -a)M (x)(b -a)
The notion of equilateral triangles then refers to lengths evaluated with this metric. In two dimensions, a positive definite symmetric matrix can be written as the product

M = RΛR t
where R is an orthogonal matrix that contains the eigenvectors n 1 and n 2 while Λ is the diagonal matrix of the positive eigenvalues. Therefore a useful geometrical interpretation of the Inria Riemanian structure defined on the domain, is to consider that on each point of the domain, the unit ball is transformed into an ellipsis whose axis are given by the two directions n 1 and n 2 and whose length h 1 and h 2 on the axis n 1 and n 2 are given by the inverse of the square root of the eigenvalues of Λ. The metric M (x) depends on the point x of the domain. Since in general, it is not possible to define analytically the relation

x → M (x)
It is then necessary to find a way to bring this information to the mesh generator. Tokamesh therefore uses the notion of background mesh to carry this information. In tokamesh, the necessary information to construct an adapted mesh is contained in three files : a background mesh in gmf format, that extends on a spatial domain containing the vacuum chamber, a solution file that contains on the nodes of the background mesh the solution of a Grad-Shafranov equilibrium and a third file containing a polygonal description of the vacuum chamber geometry as a closed list of 2D points (see section 4.2 and 5.1). The Pyamg software will use the information contained in these files to construct the desired metric and triangulation.

Delaunay and constrained Delaunay triangulation

Once the generation of iso-lines mesh is done (iso-line extractions and distribution of points along the iso-lines), a surface mesh is generated to fill the domain within the edge of the Reeb graph. The iso-line mesh then becomes an input but also a constraint as all the input edges have to match an edge of the triangular mesh. Two different approches have emerged and have proved to be robust to the complexity of the geometry: the frontal and the Delaunay methods. We focus in this report on the Delaunay approach.

The constrained Delaunay approach starts from an initial simple mesh of a box surrounding the surface mesh (composed of six tetrahedra). We then have the following steps:

1. Insert the points of composing the iso-line mesh in the current mesh; 2. Recover the boundary corresponding to the initial edge mesh (list of edges );

H k H k -C p H k+1 = H k -C p + B P
3. Fill the interior of the domain by inserting internal points;

4. Optimize the mesh with the smoothing of points and the swap of edges and faces.

Contrary to the frontal approach, a valid 2D mesh is always kept through the entire process. This is due to the insertion procedure based on a iterative process, see Figure 18. Once Step (i) is completed, some faces or edges of the initial mesh may not be present in the current mesh, a boundary recovery is used. It is generally used in enforcing these entities by applying successively or randomly standard optimization operators as the swap of edges [START_REF] George | Back to edge flips in 3 dimensions[END_REF]. In addition, some theoretical and constructive proofs exist to show that this procedure can succeed to generate a mesh, see [START_REF] George | ultimate' robustness in meshing an arbitrary polyhedron[END_REF][START_REF] Si | Tetgen, a delaunay-based quality tetrahedral mesh generator[END_REF]. The most critical step is the second one. However, if we accept to modify the initial surface mesh, this procedure can always succeed to output a surface mesh with a (slightly) modified surface mesh. Consequently, this approach is more robust than the frontal approach. Note that Step 3. and 4. are not performed in this specific setting as only the node composing the iso-lines are kept in the final mesh. Note that a lot of hybrid approaches are a combinaison of both. The frontal creation of points can be used with Delaunay insertion, or the closure of the front can used a complete constrained Delaunay approach.

Pyamg software Pyamg relies on the Feflo.a and AMG-Lib INRIA peaces of software. Feflo.a is an robust anisotropic local remeshing software. Surface and volume mesh adaptation are handled in a coupled-way. It also includes boundary layers mesh generation for RANS simulations. It includes:

• CAD re-projection and discrete surface remeshing • Hybrid mesh generation for boundary-layers • High-quality quasi-structured grids for complex geometries and complex corners: multinormals, normals deactivation, ...

• Highly anisotropic mesh adaptation, ratios: O(1:1000)

• Anisotropic/Boundary-layer coupling for supersonic shock/boundary layer interaction Inria AMG-Lib is the Adaptive Mesh Generation Library based on the Feflo.a technology. Its intent is to ease coupling with external flow solvers. In this project, Pyamg is used. It is a Python interface to the FEFLO.A/AMG-Library. However, only the 2D constrained mesh generation functionality is used. Pyamg uses Python (2.7.x or 3.x) dictionary to store mesh, solution, metric and remeshing options. For instance, the following script defines a simple 2D mesh and performs a uniform refinement:

mesh = {} mesh['xy'] = [ [0,0], [1,0], [1,1], [0,1] ] mesh['Triangles'] = [ [1,2,3,0],[3,4,1,0]]
Then, to refine this mesh with a constant size of 0.1, we have to declare a dictionary of options:

remesh_options = {} remesh_options['hmax'] = 0.1 remesh_options['logfile'] = "remesh-unif.log"
To perform the remeshing, just do: mesh_adap = pyamg.adapt_mesh(mesh, remesh_options)

For the remeshing options, a generic key exists options. It is then possible to pass a list of options directly. For instance, the previous options can be declared as: remesh_opt = {} remesh_opt["options"] = " -hmax 0.1 -logfile remesh-unif.log " Example We illustrate in this paragraph the different steps involved in the generation of the mesh of a poloidal section of a tokamak, namely the JET tokamak. The starting point is the segmentation of the domain that create a collection of subdomains. Figure 19 illustrates the construction of these subdomains (see section 7.3 in part II of this report). Inside each of the The second step of the algorithm then consists to choose inside each subdomain a collection of isolines and to discretize each isolines with a constant number of mesh points (per subdomain) as explained in section 3.2.1. For each subdomain, the segments connecting the points belonging to the same isoline are constrained to belong to the triangulation. Once all the nodes and all the constrained edges have been generated, the following Python script is used to generate the mesh for this example where cloud_mesh denote the set of generated nodes and edges.

import pyamg remesh_options = {} remesh_options['logfile'] = "gen2d.log" remesh_options["gen2d"] = [] remesh_options["MeshAllSubDomains"] = [] pymsh = pyamg.adapt_mesh(cloud_mesh,remesh_options)

Figure 20 shows the final mesh generated by this procedure containing a total of 801 nodes. In practice, the meshes used for simulations will be much more refined. This example presents intentionally a coarse mesh in order to be able to distinguish its different features.

Block-structured meshes 4.1 Domain segmentation and Morse function

The generation of block structured grids can be divided into two different tasks :

• First, one have to identify a partition of the computational domain into sub-domains such that the number of singular connection between the sub-domains is minimal and such that each sub-domain can be mapped to the square [0,1]x[0,1]. We will call this part, the domain segmentation problem.

• Second, one has to mesh in a structured way, each sub-domain ensuring some compatibility between the sub-domain (continuity of the finite element)

In the general case, the first part is by far, the most difficult part of the problem. It usually requires some manual input from the users relying on the knowledge that they have from the physical problem at hand. In the tokamak plasma community this is also the case for the simulation codes and meshing tools in use. The meshing tool CARRE [START_REF] Marchand | Carre: a quasi-orthogonal mesh generator for 2d edge plasma modelling[END_REF] relies on a a priori segmentation of the domain into several sub-domain corresponding to the known expected magnetic configurations (single null geometry, double null geometry, disconnected double null). This is also the case in the Jorek code [START_REF] Czarny | Bézier surfaces and finite elements for mhd simulations[END_REF], where the meshing tool contains different subroutines corresponding to known configurations : grid_xpoint.f90, grid_double_xpoint.f90, etc or for the SOLEDGE code [START_REF] Bufferand | Development of a fluid code for tokamak edge plasma simulation[END_REF] where a graphical user interface has been designed to help the user to define the sub-domains. While in the general case, the segmentation problem appears as a very difficult one, the construction of flux surface aligned block-structured grids can benefit from some peculiarities of the problem. Actually, it appears that all the partition strategies of the computational domain used for the generation of 2D meshes in tokamaks plasma modeling rely on some assumptions on the flux function and that it is possible to unify these approaches by a preliminary analysis of the equilibrium fields without having to use some pre-defined configurations. This analysis uses some notion of topology that we recall below.

Morse functions

Morse theory relates the study of the critical points of a function f defined on a smooth manifold Ω to the topology of this space. Introduced by M.Morse in the twenties, this theory etc In this work, we will restrict ourselves to 2D smooth manifolds and consider that Ω is a 2-D compact manifold diffeomorphic to a submanifold of the unit sphere S 2 . In order to avoid technical difficulties, we will assume that the boundary of Ω is an iso-contour of f and considering the embedding of S 2 into IR 3 we introduce a virtual point located

We begin to state some definitions and classical results.

Some definitions and results

Critical points

Let C r be the space of r differentiable scalar field defined on Ω. For r ≥ 2, a point p ∈ Ω is a critical or singular point of f if ∇f = 0. Considering a (possibly local) coordinate system around p, we therefore have :

∂f ∂x i (p) = 0 for i = 1, 2.
A critical point is regular (or non-degenerate) if the Hessian

H f = ∂ 2 f ∂x i ∂x j (p) of f at p is invertible.
A function f is a Morse function if all its critical points are regular.

The index λ(p) of a regular critical point p is the number of negative eigenvalues of the Hessian H f . A Morse function can only have isolated critical points, this is a consequence of the Morse lemma ( [START_REF] John | Morse Theory[END_REF] p 6 ) that states that it exists a local coordinate system y 1 , •, y n with y i (p) = 0 such that in the neighborhood of p,

f = f (p) -(y 1 ) 2 -• • • -(y λ ) 2 + (y λ+1 ) 2 + • • • + (y n ) 2 .
For n = 2, if we note µ 1 , µ 2 the two eigenvalues of H f the only possible critical points of a Morse function are therefore :

Maxima index = 2 µ 1 < 0 µ 2 < 0 Saddles index = 1 µ 1 < 0 µ 2 > 0 Minima index = 0 µ 1 > 0 µ 2 > 0
The following result allows to check that the identification of critical point is correct : Theorem 4.1. Let f be a Morse function defined on Ω, a region defined by a closed iso-contour, then the number of critical points (counting the virtual minimum) verify the relation :

C M -C S + C m = 2
Now let v be a vector field defined on Ω, and x 0 ∈ Ω to this vector field can be associated the ordinary differential equation

dx dt = v(x) x(0) = x 0
an orbit, streamline or trajectory (these notions are identical for autonomous odes) is a curve

x(t) = x 0 + t 0 v(u)du
These curves are everywhere tangent to the vector v.

If we associate to the scalar field f the vector field v = ∇f or the rotated vector field v ⊥ = k×∇f where k is a unit vector field of IR 3 orthogonal to Ω, we see that the iso-contours γ(s) of f defined by f (γ(s)) = C where C is constant are the orbits of the vector field v ⊥ . Conversely, the orbits of ∇f are orthogonal to the iso-contour of f .

The critical points of f are also the critical points of the vector fields v = ∇f and v ⊥ = k × ∇f . The vector field v is curl free while v ⊥ is divergence free and we have the following correspondence between the critical points Inria f maximum saddle minimum v source saddle sink v ⊥ center saddle center

Since the orbits of the divergence free vector field v ⊥ are the iso-contours of f , we use the results given in [START_REF] Ma | Structural classification and stability of divergence-free vector field[END_REF] to give a structural classification of the iso-contours of f . Proof. see [START_REF] Ma | Structural classification and stability of divergence-free vector field[END_REF] The following result also proven in [START_REF] Ma | Structural classification and stability of divergence-free vector field[END_REF] gives a refined result on the possible saddle connections.

Theorem 4.2. Let f be a Morse function defined on Ω. This field is structurally stable if and only if all saddle points are self-connected.

This last result establishes that the situation depicted in figure 4.2.1 is not structurally stable. This is intuitively easy to understand : An arbitrary small perturbation near any of the two connected saddles will change the value of f at this point and break the saddle connection. In the plasma physics application where the function f is the magnetic flux, the situation depicted in figure 4.2.1 is known as a connected double null (CDN) while the one of figure 4.2.1 is called a disconnected double null (DDN).

As remarked for instance in [START_REF] Marchand | Carre: a quasi-orthogonal mesh generator for 2d edge plasma modelling[END_REF], the CDN is only an idealization and in real experiments, the two saddles are never on exactly the same iso-contour. For the meshing algorithm, we are considering in this work, there is no necessity to consider the CDN pattern and it will be sufficient to consider the figure 4.2.1 pattern as the generic one : if the values of the magnetic flux in the two saddles points are too close, we will simply fuse the two iso-values and there will be no need to mesh the inter-region between the saddle iso-contours since its surface will be zero.

To sum up the results of this section in a concrete way, we have established that Ω consists of connected regions that contains only closed orbits. These regions are either circle cells containing an extremum critical point or circle bands separated from the other connected regions by self connected saddles. In the next section, we will give a concrete way to construct these regions and to organize the neighboring relationship between them.

Reeb graph

In the previous section, we have seen that it is possible to split the domain Ω into connected component that are either homeomorphic to a disk or to an annulus. The Reeb graph [START_REF] Reeb | Sur les points singuliers d'une forme de pfaff complètement intégrable ou d'une fonction numérique[END_REF] gives a concrete way to construct and store this decomposition of the domain. To define formally the Reeb graph, we first define an equivalence relation between points of Ω. Definition 4.1. Given a topological space Ω and a continuous function f : Ω → IR, two points p and q are equivalent p ∼ q if they belong to the same connected component of a single isocontour f -1 (c) for some c ∈ IR.

The formal definition of the Reeb graph is then : Definition 4.2. Reeb Graph : The Reeb graph is the quotient space Ω/∼ endowed with the quotient topology.

Loosely speaking, the Reeb graph concatenates all the points belonging to the same connected component of a level set into a single representative. In the case of a Morse function that have only isolated critical points, the construction of the Reeb graph can be visualized in the following way : Let us consider the graph F of the function f defined by x = (x 1 , x 2 , x 3 ) = (p, f (p)) ∈ IR 3 for p ∈ Ω. We "slice" F by the plane x 3 = C and evolve C in the positive direction starting from from -∞. Each intersection of F with the plane x 3 = C corresponds to a set of connected iso-contour f -1 = C. The change in topology of this set occurs at critical points of f . As C goes from -∞ (the virtual pit) to the global maximum of f , a new connected level set will appear when C passes trough a minimum. Conversely a connected iso-contour will disappear when C goes beyond a maximum. When C passes through a saddle, two components of the iso-contour f -1 = C will merge giving birth to two new connected components as depicted in Figure 4.2.2. The nodes of the Reeb graph then corresponds to critical points of f and each edge corresponds to a change in the number of connected components that occur when the slicing plane passes through this node. Consequently an edge originates (resp. terminates) at the nodes corresponding to a minimum (resp. maximum) and these nodes have degree 1. If the node corresponds to a saddle, two components of the level set merge and two new closed component appear. The corresponding node of the Reeb graph is therefore of degree 3 and looks like the letter "Y". For Morse functions defined on a orientable surface of genus 0, the Reeb Graph contains no loop and on a flat Ω, it is a tree. Reeb Graph are therefore also called sometimes contour tree. Due to their numerous applications in computational geometry and computer graphics, the efficient construction of Reeb graphs has been well studied. In this work, we have used the algorithm described in [START_REF] Shigeo | Algorithms for extracting correct critical points and constructing topological graphs from discrete geographical elevation data[END_REF]. Several efficient alternative exists but this one has been found sufficient for our purpose.

Inria

Meshing algorithm for C 0 meshes

We are now in position to describe our meshing algorithm. Given a domain Ω and a Morse function f : Ω → IR we first identify the critical points of f and construct its Reeb Graph. Each edge of the Reeb Graph corresponds to a connected subdomain of Ω that contains only closed level sets of f . Thus let e be an edge of the Reeb graph and e 1 , e 2 be its two extremities.

Let Ω e be the subdomain defined by e. In Ω e , the value of f evolves monotically from f m = min(f (e 1 ), f (e 2 )) to f M = max(f (e 1 ), f (e 2 )). Moreover, the domains Ω e are of only two different types : either the vertices e 1 and e 2 are two saddle points or one and only one of these two vertices, say e 1 is an extremum of f . We therefore have the following result : But Ω e can only be of two different types :

1. e is an extremal edge of the graph and one of its extemities, say e 1 is an extremum of f then Ω e is homeomorphic to a disk and we can associate to any s ∈ [0, 1] the unique p = f -1 (c) ∩ D where D is the straight line passing through e 1 and making an angle 2π × s with the x-axis.

2. e is internal edge of the graph and its two extremities are saddle points of f then Ω e is Inria a circle band and we can identify a point of f -1 (c) by its scaled curvilinear coordinate s 0 |dl|/L where dl is a differential displacement vector along the curve f -1 (c) and L its total length.

Remark :

In the first case, the mapping G is singular since on the extremum the contour line is reduced to a single point. On Ω e , the generated mesh is a polar mesh that contains a singular point at its center. In practical applications (see the example section) we either assume that the numerical method is able to handle the degenerate case of the center of a polar grid or we will modify locally the mesh in the vicinity of the center to avoid this situation.

The previous result shows that it is possible to construct a block structured grid Ω h e , aligned with the iso-contours of f into each of the subdomain Ω e . We now show that it is possible to glue these individual grids to obtain a single mesh covering all the domain. In the following, we will use a lagrange type discretization for concretness, but other type of discretization can be used. Since inside each Ω e the contour lines are closed and their values vary monotically in

[f e m , f e M ], let us choose {f e 1 , f e 2 , • • • , f e m } an increasing sequence of values ∈ [f e m , f e M ]
. This set can be choosen in a totally arbitrary manner with no influence on the discretization of the other subdomains Ω e with e = e and therefore it is sufficient to examine the discretization in the s coordinate. Let {s e j : i = j, • • • , n e } be a sequence of nodes choosen ∈ [0, 1]. Using the mapping G e , it is clear that the set Ω h e = G e (f e i , s e j ) defines a structured grid of Ω e .

Then let us consider the Reeb graph, and assume that this graph contains K extremal edge or alternatively K vertices of degree 1 (the leaves). We now construct the K paths P k in the Reeb graph connecting the K extremal vertices to the virtual pit. These K paths are composed of a unique extremal edge and a sequence of edges whose vertices are saddle points. We will call

{[m k , S k 1 ], [S k 1 , S k 2 ] • • • , [S k m , p]
} the sequence of edges composing the k th path with m k the extremal vertex and p the virtual pit. Note that this implies that the edges are ordered and that in this way we have an order relation between edges of the same path. We adopt the notation Ω [v1,v2] to designate the subdomain corresponding to the edge [v 1 , v 2 ] in the Reeb graph. To initialize the grid generation algorithm, we begin by choosing K sequences {s k

i : i = 1, • • • , n k } of nodes ∈ [0, 1] and construct the K grids Ω h [m k ,S k 1 ] = G [m k ,S k 1 ] (f k i , s k j )
corresponding to the leaves of the Reeb graph. In this construction, the only requirement imposed to these K grids will be that the saddle S k 1 belongs to the grid Ω h

[m k ,S k 1 ]
. Then the construction of the grids Ω h [Sj ,Sj+1] will be done according to the following rules :

1. if two extremal paths P k and P k meet at the saddle S j , let [S k j-1 , S j ] and [S k j-1 , S j ] be the two edges of the Reeb graph distinct of [S j , S j+1 ] that contain S j , then construct first the grids Ω h

[S k j-1 ,S-j] and Ω h [S k j-1 ,S-j]
2. The saddles S j and S j+1 belong to the grid Ω h [Sj ,Sj+1] .

the number of nodes n

[Sj ,Sj+1] used to discretize the s-coordinate in Ω h [Sj ,Sj+1] verifies : n [Sj ,Sj+1] = n [S k j-1 ,Sj ] + n [S k j-1 ,Sj ] RR n°9230

C 1 meshes

Consider a spatial domain composed of a collection of patches or subdomains corresponding to edges of a Reeb graph. The boundary of each patch is therefore a closed isoline (or a part of it) and the generic case is that a closed isoline is the boundary of two different patches. One can construct a grid or mapping of each subdomain in an independent way and "glue" the different patches together. If S i (s, t) is the mapping defining the patch i, a C 0 gluing only requires that the two curves defining the boundaries of the patches are identical, for instance :

S j (0, t) = S i (0, t) = C (t) (53)
if C (t) is the curve defining the boundaries of the patches i and j. On the contrary a C 1 gluing will require not only that the relation ( 53) is verified but also that the s-derivatives to this curve are proportional.

∂S j ∂s (0, t) = α(t) ∂S i ∂s (0, t) (54)
where α is a scalar factor. We emphasize that (54) is sufficient for C 1 continuity between the patches when only two patches are involved. If more than two patches meet at one point, additional constraints are required. However, this is the general situation for our problem, except on the x-points where more than two patches are present. Thus in the next section, we will describe how the mapping of the patches and the gluing of two patches is done when the curves separating the patches are defined by splines functions.

Grid construction of a patch enclosed by isolines

We consider a sub-domain bounded by two level sets ψ(x) = ψ m and ψ(x) = ψ M . This section explains the construction of the mapping of this region using a spline tensor product. For simplicity, we consider only two elements in the radial (orthogonal to the level sets) direction but the procedure can be easily generalized to an arbitrary number of elements. Thus, let S = [s 1 , s 2 , s 3 , s 4 , s 5 , s 6 , s 7 , s 8 , s 9 ] = [0, 0, 0, 0, 0.

Inria

Define n as the number of "elements" and let

T = [ -3 n , -2 n -1 n , 0, 1 n , • • • , n-1 n , 1, 1 + 1 n , 1 + 2 n , 1 + 3
n ] be a periodic knot vector that defines a periodic spline vector space S 3 per of dimension n with generating basis function N j (t), j = 1, • • • , n. Taking the tensor product of the two spaces S 3 ⊗ S 3 per , we can construct bivariate functions by :

S (s, t) = 5 i=1 n j=1 C i,j N i (s)N j (t) (55)
Now, let C k (t), k = 1, • • • , 5 be 5 curves defined by

C k (t) = n j=1 x k,j N j (t) k = 1, • • • , 5 (56)
We ask under what conditions the 5 curves S (s k , t) for s k ∈ {0, 0.25, 0.5, 0.75, 1} can be identical to the curves

C k (t), k = 1, • • • , 5.
That is we want to solve for the 5 × n unknowns C i,j the system:

5 i=1 n j=1 C i,j N i (s k )N j (t) = n j=1 x k,j N j (t) (57)
since the N j (t) form a basis, this gives for each j = 1, • • • , n, the 5 × 5 system :

5 i=1 C i,j N i (s k ) = x k,j for s k ∈ {0, 0.25, 0.5, 0.75, 1} (58)
Actually, (see figure 4.4.1) since N i (0) = 0 except for i = 1 (resp. N i (1) = 0 except for i = 5), we get C 1,j = x 1,j and C 5,j = x 5,j and we just have to solve the 3 equations :

5 i=1 C i,j N i (s k ) = x k,j for s k ∈ {0.25, 0.5, 0.75} (59)
Since the basis functions are of compact support, these equations gives a 3 × 3 system that writes explicitly (see figure 4

    C 2,j C 3,j C 4,j   =   x 2,j -0.125x 1,j x 3,j x 4,j -0.125x 5,j   (60)
However, this procedure does not give any means to control the curves S (s, t * ) with t * fixed. Obviously, one can reverse the role of s and t and instead of the

5 curves C i (t), i = 1, • • • , 5 consider n curves T j (s), j = 1, • • •
, n in the s direction. However this time, the procedure will not control the curves in the t direction.

Thus instead, we propose a procedure where one control some of the s = constant curves and some of the t = constant curves. More precisely, we will impose that the three curves S (s k , t) for s k ∈ {0, 0.5, 1} are identical to the three curves C k (t) for k = 1, 3, 5. As seen previously, this requirement sets the control points C 1,j and C 5,j and we are left with the system (60) for s k = 0.5 that writes :

C 2,j 2 + C 3,j + C 4,j 2 = 2 x 3,j ∀j = 1, • • • , n (61)
This system gives explicitly C 3,j if C 2,j and C 4,j are known. To compute these last control points, we will require that for each t the derivatives at s = 0 and s = 1 of the curves S (s, t) have known values, say x (0, t) and x (1, t). The rationale of such a requirement is to insure the C 1 continuity of the curves t = constant when passing from one patch to a neighboring one.

The system that we have now to solve is thus

∂ s (S (s, t))(0, t) = x (0, t) ∀t ∂ s (S (s, t))(1, t) = x (1, t) ∀t (62)
Computing the derivatives of the basis functions N i (s) in 0 and 1 gives explicitly :

n j=1 6(C 2,j -C 1,j )N j (t) = x (0, t) ∀t n j=1 6(C 5,j -C 4,j )N j (t) = x (1, t) ∀t (63)
To solve (63)we will develop x (0, t) and x (1, t) on the N j (t) basis. Thus, let x 0,j and x 1,j be the control points defining the functions x (0, t) and x (1, t) :

x (0, t) = n j=1 x 0,j N j (t) ∀t x (1, t) = n j=1 x 1,j N j (t) ∀t (64)
The equation (62) gives :

6(C 2,j -C 1,j ) = x 0,j ∀j 6(C 5,j -C 4,j ) = x 1,j ∀j (65)
and can be solved for C 2,j , C 4,j . Now solving for C 3,j thanks to (61) defines completely the spline surface S (s, t). On this surface, the curves S (s k , t)), s k ∈ {0, 0.5, 1} exactly coincide with the three curves C k (t), k ∈ {1, 3, 5} and the definition of the set of control points C 2,j , C 4,j ensures that the derivatives of the curves S (s, t)) in s = 0 and s = 1 are controlled. This procedure is illustrated in figure 4.4.1. The black lines represents the curves S (s * , t) for s * fixed while the blue ones are the curves S (s, t * ) for t * fixed. Observe that the three controlling curves C k (t) (in red) are exactly superposed with the curves S (s, t)) (in black) for s = 0, 0.5, 1 defined on the surface. Here the derivatives at s = 0 and s = 1 have been chosen to be proportional to the gradient of the function ψ. This makes the grid locally orthogonal in the vicinity of the two level sets S (0, t)) S (1, t)) as can be seen on figure 4.4.1. As already mentioned, this procedure can be easily generalized to an arbitrary number of elements in the s direction. The only change is that instead of the equation (61), we will have to solve a tridiagonal linear system. In this way, one can obtain a grid whose S (s * , t)) curves are exactly aligned on an arbitrary number of level sets.

C k (t) = n j=1 x k,j N j (t) k = 1, 2 (68)
we require the 2 curves S (0, t) and S (1, t) to be identical to the curves C k (t), k = 1, 2. This immediately gives :

C 1,j = x 1,j C 4,j = x 2,j ∀j
To defines the surface, again as in the previous section, we will require that the derivatives in s = 0 and s = 1 of the curves S (s, t) match two prescribed directions d 1 (t) and d 2 (t) defined by two sets of control points d 1,j and d 2,j defined on the basis N j (t). Given the expression of the derivatives of the basis functions given in 66, we obtain :

C 2,j = C 1,j + d 1,j /3 C 3,j = C 4,j -d 2,j /3 ∀j (69)
However, in section 4.4.1, the direction fields d 1,j , d 2,j can be defined by the gradient of the magnetic flux field. This is no more possible on the C 1 (t) curve since the gradient is not defined on the x-point and is a highly discontinuous quantity. We have therefore found convenient to define the d 1,j control points as an average between the normal vectors of the two curves

C k (t), k = 1, 2 : d 1,j ← αd 1,j + (1 -α)d 2,j (70)
where d 1,j and d 2,j are the control points of the gradient field along the curves C k (t), k = 1, 2 and α is a weighting factor null on the x-point and close to 1 far away from this point. In this way, we smooth out the discontinuity on the x-point and construct a surface where only the elements directly touching the x-point are not C 1 . This can be verified on figure 4.4.2 that display the grid between the isoline ψ = -1 and ψ = -0.954 . The element is quite narrow in order to confine the width of the discontinuity to a very small extent in the s direction. As can be seen on the enlarged view of figure 4.4.2, in the t direction, the derivative discontinuity is confined to the first element near the x-point.

C 1 gluing of two patches

We consider in this section the problem of constructing a global C 1 mapping of a domain described by a collection of patches. This is in general a difficult problem. However for domains described by a Reeb graph, except on saddle points, we need only to stick two patches since a patch has only one direct neighbor. This sticking of two patches can be realized using parametric continuity (C 1 ) or a geometric continuity (G 1 ). For the C 1 condition, the derivative vectors Inria are in opposite directions with the same magnitude, while for the G 1 the derivative vectors are in opposite directions but can be of different magnitude. As explained in section 1.3.6, these conditions translate on geometric conditions on the control points requiring the alignment of three control points. When only two patches have to be glued, this condition is therefore easy to satisfy. As a first example, consider the construction of a regular grid inside a region limited by the separatrix. As explained in the previous sections (4.4.1 and 4.4.2) since we want the derivative discontinuity generated by the x-point to extend in the smallest region possible, we will generate two patches : S c (s, t) an inner surface for the core region with C 1 regularity and S e (s, t) a one-element surface containing the x-point where the derivative discontinuity is limited to the 2 elements containing the x-points. The common boundary of the two surfaces is a spline curve C ( t) defined by the control points P 0,j (purple curve on figure 4.4.3). Given for instance, the control points P c 1,j that define the s-derivative on the boundary of S c (s, t), it is easy to define by the equation [START_REF] Si | Tetgen, a delaunay-based quality tetrahedral mesh generator[END_REF] the P e 1,j control points in the patch S e (s, t) using for instance a G 1 condition. patches is the curve approximating the separatix (cf. the 1 st step). Note that here, to define the normal derivative to the separatrix, we have chosen to use the gradient of the flux function. This create a local orthogonal grid in this region.

Conversion to Hermite-Bézier mesh

The previous construction uses spline tensor patches defining a global surface. Once this is done, one can convert these spline surfaces into an Hermite-Bézier mesh using an arbitrary number of elements defined by the user. This transformation does not change the surface and therefore the G 1 character of the global ssurface is maintained. In particular, this construction implies that the control points corresponding to mixed derivatives are aligned across the elements corners as seen in figure 4.4.4

Figure 35: Hermite-Bézier mesh (for CORE region) with control points of mixed derivarives.

Inria pip install Scipy

• Matplotlib : a python 2D plotting library pip install Matplotlib

• IgaKit : a toolkit for IsoGeometric Analysis that implements many of the NURBS routines from Piegl's book using Fortran and Python.

pip install https://bitbucket.org/dalcinl/igakit/get/default.tar.gz

Installing the library

Getting the source code

The source code is located on the Inria GitLab platform. Thus first create an Inria Gitlab account as an external user. Then search for the EoCoE/tokamesh project and as External user you will be allowed to download the project, have access to the documentation and the error tracking information. You can clone the tokamesh project using: git clone git@gitlab.inria.fr:EoCoE_Mesh/tokamesh.git cd tokamesh Or download it form the project page on the git platform : https://gitlab.inria.fr/ users/sign_in.

Customization

Next, define a custom Path Variable enabling your system to locate the the source code : Add the following line to your * .bashrc * , * .bash_profile * or * .profile * file: export TOKAMESH_DIR=MY_SOURCE_CODE_PATH where MY_SOURCE_CODE_PATH is the path where you have clone the project.

Pyamg

After downloading the tokamesh library, you will find in the tokamesh folder an 'external' folder, which contains the archive of pyamg an external package for mesh generation. To install pyamg : Select the compressed folder according to your Operating System, unzip it and install pyamg :

• For Linux:

• Development mode :

python -m pip install --user -e .

Uninstall:

To uninstall the library, in a terminal, go to the main folder of the library, and type in : python -m pip uninstall tokamesh 6 Organization of the software tokamesh relies mainly on two different tools. The first one is Pyamg, a generic software developed by Adrien Loseille at Inria and built upon Feflo.a. Feflo.a is an robust anisotropic local remeshing software written in C that can handle highly anisotropic mesh adaptation (anisotropy ratio of more that O(1 : 1000). Pyamg provides a python interface with Feflo.a. The second is IgaKit written by Lisandro Dalcin at KAUST that includes many algorithms to create and handle spline and NURBS curves and surfaces. Beside these basic tools, tokamesh is composed of a collection of python modules representing approximately 12k lines of codes. In order to ease the construction of meshes and grids, tokamesh have been organized in different subsets contained in distinct directories.

The following sections give a presentation of the organization of the library. Its main parts are:

• doc • data • examples • externals • tokamesh
Figure 36: Tokamesh library main folders

Doc

This folders contains the documentation about the library and includes the present report. The auto-documentation of the code uses sphinx (http://www.sphinx-doc.org/en/master/). The subfolder sphinx present in the doc folder contains the html code documentation generated by sphinx using the doc_string of individual python modules and functions.

Data

In tokamesh, the inputs to the python scripts are gathered in a folder data/inputs that contains subfolders. Each subfolder corresponds to the definition of a different case. As can be seen in figure 37, data/inputs already contains several examples corresponding to different tokamaks and cases. To treat a new problem, the user then must create a new folder in data/inputs.

Then loading a problem can be done with the single instruction line :

filename="data/inputs/ * /.../name_of_case" data=load_data(filename)

provided the last folder contains three files respectively called name_of_case.mesh, name_of_case.sol and name_of_case.vac whose definition is as follows :

• a .mesh file containing the list of points ( by coordinates), edges (doublet of points) and elements of the mesh • a .sol file containing the list of,value of a function Z(x,y) on each vertices of the mesh Inria and the 'End' keyword at the bottom end of the file. Note that to use the tokamesh library, you do not need to download the LibMesh library, just to adhere for the definition of the .mesh and .sol files to the convention used in the LibMesh library. The . vac is specific to this library. Here is an example of .vac file: The first line indicate the number of vertices defining the boundary of the vacuum chamber. The following lines gives the coordinates (two reals) of the vertices defining the boundary plus an additional integer flag that the user can use to distinguish different parts of the boundary if needed. The list of coordinates present in a .vac have to be organized in such a way the vertices are consecutive : the segment a i and a i+1 is assumed to be a part of the vacuum chamber boundary as well as segment a n and a 1 (closed curve). It is assumed that the coordinates in a .vac file are unique and that there is no repeated nodes. In order to ensure that this is the case, the vacuum_utils module contains a python script : check_vac_bdry_double allowing to verify this point.

Note also that if the data corresponding to a case is in .eqdsk format, the EQDSKtoMesh Fortran tool allows you to generate the .mesh, .sol and .vac files from it. This tool can be found in the topology folder.

External

In this folder are stored the binaries of the external tools used by the tokamesh library : Pyamg, EqdskToMesh and transmesh. For more information about these tools, the first two have been previously described, the last one is described in section 7.4.2

tokamesh

The tokamesh folder contains the core of the library, the different subfolders contained in tokamesh are :

• preprocess : tools for initial data analysis (meshes and solutions),

• bsplines : functions for b-splines generation • meshgen : functions for mesh generation • hbezier : functions for Hermite-Bézier splines generation • utilities : various functions used by other tools

Examples

A small set of examples showing how to use the different part of the library. For instance, the analysis and segmentation of a problem named "pb" is demonstrated by running the python script Inria python nameoffile.py in the corresponding folder "pb". It is advisable to do not erase the folders present in Examples as they are used to provide reference results.

Pre-processing tools

Various preprocessing tools are included in the library, which you can find in the tokamesh/preprocess folder of it. They can be divided in two parts: data analysis and preparation, tools related to the definition of the vacuum chamber as well as function for domain decomposition. We present here the main functions, which use several smaller routine , for more information you can check the documentation generated with sphinx in the doc folder .

Data analysis

load_data

The first par of an usual Tokamesh routine is to check if the given input data are conform. The load_data(filename) routine contained in the load_data.py file, read the input data and return them in a data format. If they are not conform to the format it return an error. Example:

filename="/data/inputs/WEST/west_gmf" data=load_data(filename) mesh=data[0] z=data [START_REF] Alfeld | Smooth macro-elements based on clough-tocher triangle splits[END_REF] vac_nodes= data [START_REF] Boisserie | Sur l'implementation des éléments finis de hseieh-clough-tocher complet et reduit[END_REF] Here we load the data contained in west_gmf.mesh, west_gmf.sol and west_gmf.vac and stack them in data variable.

data_preparation

The data can be refined using the data_preparation function ( located in the preparation.py module. This function works as follow: First the critical points of the magnetic flux function are identified, then the x-points that are present inside the vacuum chamber are detected. This step take into account that some data sets use meshes that cover not only the vacuum vessel but also a large part of the tokamak. At present, since tokamesh have been designed for x-point configurations, if there are no saddle point inside the vacuum chamber, the function returns an error and stops. Then the mesh around the saddle points which are inside the vacuum chamber is refined, and finally the solution is interpolated on the new refined mesh using Clough-Tocher C 1 interpolation. This phase of preparation of the data has been found necessary in many cases to obtain better and cleaner data around the saddle points.

Here is an example of call: mesh_r,z_r=data_preparation(mesh,z,vac_nodes) """ writing the refined . mesh and .sol """ filename_r=filename+'prepared' write_mesh_sol_from_tri_mesh(filename_r,mesh_r,z_r) WriteVacuum(filename_r,vac_nodes) It return a new mesh, and new solution. Here we write down the new mesh, solution and vacuum border in new files, by using two routines: write_mesh_sol_from_tri_mesh and WriteVacuum

Vacuum chamber tools

Being able to describe the vacuum chamber and to know whether we are inside or outside of it is an important part of the Tokamesh library functions. You will find described here some of the useful routines for it.

is_in_vaccum_pt

For a random point of the 2d space, says if it is inside the vacuum chamber Arguments :

• pnt : Pt2D a 2D set of coordinates (x,y)

• bdry_list : an ordered list of 2D point defining the vacuum vessel.

Return : true if i is inside the vacuum, false otherwise.

interpolate_on_mesh

Read a function defined on a triangulation and return the value of the function on a list of points inside the triangulation . It can be used to define the value of the function on the points defining the vacuum chamber border if they aren't triangulation nodes, or to calculate the value at newly inserted points. Arguments:

• pt : list of mesh nodes • el : list of mesh triangles (triplets of int)

• z : values of a function on the point pt • pt_int : list of points where we want to interpolate the function • kind : interpolation method possible values : linear or cubic Return: the list of value of z on pt_int points given by interpolation

curve_refinement

In some data sets, the vacuum chamber is defined by a very coarse discretization containing very large as well as very small edges. This situation might be problematic for the generation of a new mesh using the vacuum chamber as external boundary. To solve this problem, the curve_refinement function compute a new distribution of nodes on the vacuum vessel boundary. Given an ordered list of points coordinates (in the Pt2D format), this function returns a new list of points in the same format, containing the old list plus new points inserted so that the edges of the boundary are of a regular size. Figure 38 show an example on the JET vacuum chamber .

Segmentation

The function segmentation present in the module preprocess/segmentation.py create a new mesh containing only the vacuum chamber boundary and the isolines passing through the saddle points present in the vacuum vessel. In addition, these isolines are sampled uniformly and the function also returns a list of critical points present in the vacuum vessel. This is necessary to handle the case where the data extends far beyond the vacuum vessel and contains also critical points related to the presence of external coils. This function has to be called to generate the different subdomains that will be meshed by the different meshing functions. Here is an example of call: sampling_coeff = 0.1 min_pt = 8 """ compute the domain segmentation """ seg_mesh, seg_sol,x_pt_i,o_pt = Segmentation(mesh,z,vac_nodes,\\ sampling_coeff , min_pt )

Where sampling_coeff and min_pt are parameters used for an uniform sampling of the isolines. isoline, isolineedge, l_newbdry = Gen_iso_from_ver(node,mesh, z,\\ iso_val, new_bdry_pt) new_bdry_pt = l_newbdry list_iso_saddle.append(isoline) list_iso_sad_edge.append(isolineedge) print("the isoline passing through x point of index",node) print("has", len(isolineedge), "separated branches") print(" -----------------------------------------------------------" Using those isolines, the border and the center point ( corresponding to the absolute min/max in the vacuum chamber), we call the routine get_void_mesh to generate a triangular mesh, with several ref corresponding to the different branch of the Reeb graph of our problem, as in figure 40 """creating a mesh structured around the vaccum border + the saddle isolines""" seg_mesh, seg_sol = get_void_mesh(l_nodes,l_edges,center_pt) flag_t=[t [START_REF] Bufferand | Development of a fluid code for tokamak edge plasma simulation[END_REF] pyamg.adapt_mesh : This function generate a mesh in the form of a python dictionary which contains its vertices, edges, triangles and corners. That mesh is in a compatible format with the libmesh data format. Arguments :

• mesh : a python dictionary consisting of :

-"xy" , a list of doublets corresponding to the 2D coordinates of the mesh vertices -"edges", a list of doublets corresponding to the edges between vertices we want the generated mesh to have : "edges" = [..., [START_REF] Bufferand | Development of a fluid code for tokamak edge plasma simulation[END_REF][START_REF] Clough | Finite element stiffness matrices for analysis of plates in bending[END_REF],... ] force the pyamg.adapt_mesh routine to generate a mesh where vertices 3 and 5 are linked by an edge.

-"corner", a list of vertices index which indicate the mandatory vertices of the mesh. Which mean that any mesh optimization or adaptation routines will be forced to keep those points.

• remesh_options

Inria

Output: pymsh : a mesh in pyamg mesh data storage format pyamg.write_mesh: This function write down a mesh in a .mesh(b) file format . Arguments :

• mshtab : a mesh in pyamg mesh data storage format.

• filename : the name of the file where it is going to be written, including the filetype (.mesh or .meshb)

Transmesh

This tool developed by Loic Maréchal allows you to transform libmesh binary format to ASCII format and the other way round. it works as follow :

./transmesh source_name destination_name

EqdskToMesh

This Fortran program transforms .eqdsk files into mesh files. From one Eqdsk file it produces three files:

• a .mesh file containing the list of points ( by coordinates), edges and elements of the mesh • a .sol file containing the list of,value of a function Z(x,y) on each vertices of the mesh • a .vac file containing the list of coordinates of the vertices describing the vacuum chamber boundary.

To compile it use the following line :

gfortran EqdskToMesh.f90 -o EqdskToMesh

To use it on a case, named mycase. In the same repository where the EqdskToMesh programm is, create a folder named "mycase" : mkdir mycase run the program on your case :

./EqdskToMesh mycase 7.5 Mesh generation

patch_connectivity

The function patch_connectivity present in the module meshgen/patch_definition takes as argument the mesh resulting from running the segmentation function. It returns the connectivity between the different subdomains (or patches) and the definition of the curves separating the different patches either as a list of points or as cubic spline curves.

Figure 1 :

 1 Figure 1: Triangle HCT réduit

Figure 2 :

 2 Figure 2: Interpolation for the function v(x, y) = x 4 + y 4

Inria

Figure 4 :4

 4 Figure 4: (left) A quadratic B-Spline curve and its control points (the knot vector is {000 1 2 3 4 3 4 111}), (right) the corresponding B-Splines.

Inria

Figure 5 : 4 ,

 54 Figure 5: (a) A quadratic B-spline curve and its control points. The knot vector is T = {000, 1 4 , 1 2 , 3 4 , 111}. (b) The curve after a h-refinement by inserting the knots {0.15, 0.35} while the degree is kept equal to 2.

Figure 7 :

 7 Figure 7: Impact of the parametrization: the discrete data, the approximated curve depending on the parameters (u k ) -(a): using uniform distribution, (b): using Gauss-Legendre points.

Figure 8 :

 8 Figure 8: b-spline basis function associated to the knot vector S = [0, 0, 0, 0, 1, 1, 1, 1]

Figure 9 :

 9 Figure 9: Bezier curve with a loop.Left : curve, middle : derivative, right : the curves 50 for large α showing that 50 have roots in [0, 1].

Figure 10 :

 10 Figure 10: f (x) = x 3(x -2)

Figure 12 :

 12 Figure 12: Level set passing through the x-point for GS equilibrium computed by EFIT

Figure 14 :

 14 Figure 14: Level set passing through the x-point after local refinement and CT interpolation

Figure 17 :

 17 Figure 17: left : Distribution of nodes in the SOL region for a jet equilibrium, right : zoom on the lower region.

Figure 18 :

 18 Figure 18: Illustration of the incremental Delaunay insertion of a point in a mesh.

Figure 19 :

 19 Figure 19: Segmentation of the vacuum chamber of the JET tokamak subdomains indicated on figure 19 by a different color, the level sets of the flux function are

Figure 20 :

 20 Figure 20: Illustration of the constrained Delaunay mesh generation from an initial set of isoline mesh.

Theorem 4 . 1 .

 41 Let f be a Morse function defined on Ω. Then the topological set of the iso-contours of f consists of finite connected components that are either • Circle cells which are homeomorphic to open disks • Circle bands which are homeomorphic to open annulus • Saddle connections

Figure 21 :

 21 Figure 21: A connected double null in the MAST tokamak

Figure 22 :

 22 Figure 22: A disconnected double null in the DIIID tokamak

Figure 23 :Figure 24 :Theorem 4 . 1 .

 232441 Figure 23: Behavior of the level sets when the plane x 3 = C passes through a saddle : For C < saddle value, the level set is the single blue line, at C = saddle value, the two components merge (green line) giving rise to two connected components for C > saddle value (red line)

5 , 1 , 1 , 1 , 1 ]

 51111 be an open-knot vector that defines a spline vector space S 3 of cubic splines with two elements. The b-spline basis functions N i (s), i = 1, • • • , 5 generating this space are represented in figure 4.4.1.

Figure 25 :

 25 Figure 25: b-spline basis function associated to the knot vector S = [0, 0, 0, 0, 0.5, 1, 1, 1, 1]

Figure 26 :

 26 Figure 26: Tensor product surface of a patch enclosed by two level sets.

Figure 28 :Figure 29 :

 2829 Figure 28: One-element surface containing an-x-point

Figure 30 :

 30 Figure 30: G 1 regular grid inside the Separatix.

Figure 31 :

 31 Figure 31: C 1 continuity of splines curves across the X-point.

Figure 32 :

 32 Figure 32: A regular grid of the SOL region.

Figure 33 :

 33 Figure 33: C 1 stick of CORE/Edge and SOL grids.

  (a) Partially inside vaccum vessel. (b) Fully inside vaccum vessel.

Figure 34 :

 34 Figure 34: Full mesh of CORE/Edge/SOL regions.

  cp pyamg-inria.linux64bits.tgz $MY_INSTALL_PATH cd $MY_INSTALL_PATH tar zxvf pyamg-inria.linux64bits.tgz • For macOS: Inria cp pyamg-inria.macosx.tgz $MY_INSTALL_PATH cd $MY_INSTALL_PATH tar zxvf pyamg-inria.macosx.tgz where MY_INSTALL_PATH is the path where you want the PyAMG package to be installed. You can use for example: " $TOKAMESH_DIR/usr " (after creating usr folder). Finally, add the following line to your *.bashrc*, *.bash_profile* or *.profile* file: export PYTHONPATH=$PYTHONPATH:MY_INSTALL_PATH export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:MY_INSTALL_PATH/pyamg/pyamg/dep 5.3.4 Installation: The last step installs the library : open a terminal, go to the main folder of the library , and type in one of the following command line: • Standard mode : python -m pip install .

Figure 37

 37 Figure 37: data folder

Figure 38 :

 38 Figure 38: left : initial JET border, right : refined JET border.

Figure 39 :

 39 Figure 39: Isolines branches as generated by Segmentation

  ) z_min=np.amin(z) z_max=np.amax(z) delta_z=z_max-z_min x_value=[c.valeur for c in x_pt_i] if len(x_value) >1 : if abs(x_value[0]-x_value[1]) < 1.e-5 * delta_z : list_iso_saddle,list_iso_sad_edge=glue_isolines(list_iso_saddle,\\ list_iso_sad_edge)

  Figure 40: Basic mesh generated by Segmentation

RR n°9230

Actually, the use of accurate ode solver is required when one need to reach an isoline close to the magnetic axis. Since in this region, ∇ψ becomes small, the integration of (6) requires a great accuracy. In practice, for all the other cases the integration of (6) can be done with cheaper methods RR n°9230

For the sake of simplicity, we have considered here the case where the multiplicative coefficients of the derivatives are equal but in practice they can be different if needed RR n°9230

The flux function is normalized between -1 on the x-point and 0 on the magnetic axis RR n°9230

Inria

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgments

This work has been partly realized in the framework of the H2020, Energy Oriented Center of Excellence for computing applications, EoCoE project -http://www.eocoe.eu/ under grant agreement EINFRA-676629. In particular Alexis Loyer and Jalal Lakhlili have been directly supported by EoCoE. Ali Elarif has been supported by a PhD grant from Inria. The authors also acknowledge the support of Inria through the Inria IPL program FRATRES https://team.inria.fr/ipl-fratres/

Grid construction of a patch enclosed by non C 1 isolines

The above procedure is appropriate for the case where the two isolines ψ(x) = ψ m and ψ(x) = ψ M have C 1 regularity. However when one of the two isolines enclosing the patch passes through an x-point, the previous strategy have to be adapted since the curve ψ(x) = ψ m is no more C 1 and cannot be expanded on a basis of periodic splines. Moreover, the normal to ψ(x) = ψ m is not defined on the x-point and has a discontinuity on this point making the spline approximation of the gradient defined in 64 problematic. We therefore have to modify the previous algorithm to cope with this case.

A second problem, is that we would like the loss of C 1 regularity to be as localized as possible. Ideally, we would like this loss of regularity to be confined only in the elements surrounding the x-point. Thus, we will in this section consider the construction of a surface composed by a unique element in the s direction enclosing a region defined by 2 isolines ψ(x) = ψ m and ψ(x) = ψ M where one of them ,say for instance ψ(x) = ψ m ) is not C 1 . The vector space of cubic splines defined by a unique element and the open-knot vector[0, 0, 0, 0, 0, 1, 1, 1, 1] reference, the expression of these functions and of their derivatives are given below :

In the t direction, we cannot use periodic splines, thus we use clamped splines defined on the open-knot vector T = [0, 0, 0, 0,

where n is the number of elements. As in the previous section, we consider the surface

Inria

Part III

User's manual 5 Installation

The tokamesh library uses a few external packages written in C and Fortran as well as a large number of modules written in Python. To use the tokamesh Library, you first need to install Python on your machine, and these few external packages. The following describes steps by steps this installation.

Installing Python

Python is installed by default on all major distributions. If this is not the case on your machine proceed as follows, otherwise go to step 5.2 Then you can use homebrew to install python and the package easily brew install python@2 export PATH="/usr/local/opt/python@2/libexec/bin:$PATH"

Package installation

Installing python also install automatically pip, the python package manager. Except the last one (IgaKit) the following packages are usually available by default. If they are not present on your machine use pip to install them.

• Numpy : the fundamental package for scientific computing with Python.

pip install Numpy

• Scipy : a Python-based ecosystem of open-source software for mathematics, science, and engineering.

First it detect the critical points inside the vacuum chamber, and store those who are inside the vacuum chamber in two list : one for the min/max and one for the saddle points :

""" calculate critical points in mesh""" cp= CriticalPoints(pt,el,z) """ get the x-points and extremal points inside the vaccum chamber """ M_pt_i.append(c)

The calculation of the critical points is done using the CriticalPoints routine, which can be found in the preprocess/morse_lib.py file . To know if one of critical points is inside the vacuum chamber, we use the is_in_vaccum_pt routine which given a contour described by a series of ordered points can return if a random point is inside it or not . This routine can be found in the preprocess/vaccum_utils file.

Then for each saddle point, it generate the isoline going through it, restricted to the vacuum chamber border. For this we used the routine Gen_iso_from_ver from the file preprocess/ret_iso. From a node of the mesh, this routine generate the isoline corresponding to the value of the solution field at this point. An isoline is defined by two list:

• isoline containing the coordinates of the points describing it • isolineedge containing doublet of index (called "edges") , which describe how those points are connected to each other. If the isoline is split in several part by the border of the vacuum chamber, the edges are stored in sublist inside isolineedge corresponding to each branch.

Example: in the figure 39, each color correspond to a different "branch" of an isoline .

If two saddle point have a value different from less than 0.00001 of the amplitude of the flux function, we decide that they form a double null configuration, and fuse their isoline, using the glue_isolines routine.