Théodore Cherriè

Luc Laurent

François Sami Hlioui

François Louf

Pierre Duysinx

Christophe Geuzain

Hamid Ben

Mohamed Gabsi

Eduardo Fernández

Théodore Cherrière
email: theodore.cherriere@ens-paris-saclay.fr

Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator

Keywords:

This work uses Multi-Material Topology Optimization (MMTO) to maximize the average torque of a 3-phase Permanent Magnet Synchronous Machine (PMSM). Eight materials are considered in the stator: air, soft magnetic steel, three electric phases, and their three returns. To address the challenge of designing a 3-phase PMSM stator, a generalized density-based framework is used. The proposed methodology places the prescribed material candidates on the vertices of a convex polytope, interpolates material properties using Wachspress shape functions, and defines Cartesian coordinates inside polytopes as design variables. A rational function is used as penalization to ensure convergence towards meaningful structures, without the use of a filtering process. The influences of different polytopes and penalization parameters are investigated. The results indicate that a hexagonal-based diamond polytope is a better choice than the classical orthogonal domains for this MMTO problem. In addition, the proposed methodology yields high-performance designs for 3-phase PMSM stators by implementing a continuation method on the electric load angle.

Since the seminal work of Bendsøe and Kikuchi [START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF], Topology Optimization (TO) has experienced considerable advances driven by the widespread interest in industry. The achievement owes to the high number of contributions that have improved computational efficiency [START_REF] Aage | Topology optimization using PETSc: An easy-to-use, fully parallel, open source topology optimization framework[END_REF] and expanded its applicability beyond solid mechanics to other physics [START_REF] Bendsøe | Topology optimization -Broadening the areas of application[END_REF]. Among different TO approaches [START_REF] Sigmund | Topology optimization approaches[END_REF], such as level-set the density method stands out. In the context of structural optimization, each mesh element is associated with a design variable (called density) that takes a value of 0 to represent a void phase or a value of 1 to represent a solid phase. After solving the optimization problem, solid elements define the optimized design. In order to obtain a well-posed problem, design variables are interpolated continuously between 0 and 1. The relaxation also enables sensitivities computation in order to use efficient optimization algorithms such as gradient descent, thanks to the adjoint variable method. To avoid intermediate values in the final design, these are penalized, for instance, by a power-law [START_REF] Bendsøe | Optimal shape design as a material distribution problem[END_REF] or a rational function [START_REF] Stolpe | An alternative interpolation scheme for minimum compliance topology optimization[END_REF]. A substantial extension of density-based methods called Multiple Materials in Topology Optimization (MMTO) appeared with the work of Thomsen [START_REF] Thomsen | Topology Optimization of Structures Composed of One or Two Materials[END_REF]. Notably, these methodologies have been widely used in compliance minimization problems that consider isotropic materials with different Young's moduli. Advances in MMTO have spanned different fields involving the design of compliant mechanisms [START_REF] Gaynor | Multiple-material topology optimization of compliant mechanisms created via PolyJet three-dimensional printing[END_REF], lattice structures [START_REF] Kazemi | Multi-material topology optimization of lattice structures using geometry projection[END_REF], thermomechanically constrained porous composites [START_REF] Takezawa | Design methodology for porous composites with tunable thermal expansion produced by multi-material topology optimization and additive manufacturing[END_REF], structural joints [START_REF] Jung | Multi-material topology optimization considering joint stiffness using a two-step filtering approach[END_REF], piezoelectric structures [START_REF] He | Multi-material topology optimization of piezoelectric composite structures for energy harvesting[END_REF], among many others. Other MMTO methodologies include level-set approaches [START_REF] Wang | Color" level sets: A multi-phase method for structural topology optimization with multiple materials[END_REF][START_REF] Allaire | Multi-phase structural optimization via a level set method[END_REF], phase field approaches, such as Cahn-Hilliard dynamics [START_REF] Zhou | Multimaterial structural topology optimization with a generalized Cahn-Hilliard model of multiphase transition[END_REF] or accelerated constrained Allen-Cahn Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator. Structural and Multidisciplinary Optimization volume, 65, 2022. doi: 10.1007/s00158-022-03460-1 dynamics [START_REF] Tavakoli | Multimaterial topology optimization by volume constrained Allen-Cahn system and regularized projected steepest descent method[END_REF][START_REF] Tavakoli | Alternating active-phase algorithm for multimaterial topology optimization problems: A 115-line MATLAB implementation[END_REF]. Recently, Chandrasekhar and Suresh [START_REF] Chandrasekhar | Multi-Material Topology Optimization Using Neural Networks[END_REF] used neural networks to represent the volume fraction of different materials.

In the magnetostatics field, MMTO emerged with the work of Dyck and Lowther [START_REF] Dyck | Automated design of magnetic devices by optimizing material distribution[END_REF]. The authors optimized a magnetic bearing, composed of steel, air and conductors, by interpolating each physical property (magnetic permeability and current density). They also eliminated intermediate materials with penalization. Concerning density methods, the most common applications are C-core electromagnets as in Wang, Park, and Kang [START_REF] Wang | Multi-domain topology optimization of electromagnetic systems[END_REF] and Wang et al. [START_REF] Wang | Topology optimization of electromagnetic systems considering magnetization direction[END_REF], which contain steel, air and a source -current density, Permanent Magnets (PM) or both. The application on which the highest number of materials are used is the design of Halbach arrays, where different PMs are distributed to maximize the flux density. This problem was tackled in 3D [START_REF] Choi | Optimization of magnetization directions in a 3-D magnetic structure[END_REF] and with arbitrary PM directions [START_REF] Lee | Topology optimization for three-dimensional design of segmented permanent magnet arrays[END_REF], which were included within the optimization variables. Steel was added in the optimization to design moving actuators. Lee, Dede, and Nomura [START_REF] Lee | Simultaneous Design Optimization of Permanent Magnet, Coils, and Ferromagnetic Material in Actuators[END_REF] placed several materials to maximize the force profile of a plunger. This work was extended in Lee et al. [START_REF] Lee | Multi-material topology optimization of magnetic actuator with segmented permanent magnets[END_REF] to include more PMs and in Jung, Lee, and Lee [START_REF] Jung | Design and Fabrication of Magnetic System Using Multi-Material Topology Optimization[END_REF] to produce manufacturable structures.

Concerning MMTO applied to electrical machines, many works use heuristic approaches as in Sato, Watanabe, and Igarashi [START_REF] Sato | Multimaterial topology optimization of electric machines based on normalized Gaussian network[END_REF] and Jung, Ro, and Jung [START_REF] Jung | A Hybrid Algorithm Using Shape and Topology Optimization for the Design of Electric Machines[END_REF], mainly to optimize only the rotor of electrical machines. Density-based methodologies are also used to optimize rotors of different kinds of actuators, such as wound field synchronous machines [START_REF] Guo | Multimaterial Magneto-Structural Topology Optimization of Wound Field Synchronous Machine Rotors[END_REF], permanent magnet synchronous reluctance motors [START_REF] Gauthey | Multi-Material Topology Optimization with Continuous Magnetization Direction for Permanent Magnet Synchronous Reluctance Motors[END_REF][START_REF] Lee | Multi -material topology optimization for the PMSMs under the consideration of the MTPA control[END_REF], or permanent magnet synchronous machines [START_REF] Park | Multi-material topology optimization of permanent magnet synchronous motors[END_REF], and can be hybridized with shape optimization [START_REF] Kuci | Combination of topology optimization and Lie derivative-based shape optimization for electro-mechanical design[END_REF]. By contrast, very few works tackle the stator optimization. Among them, Lee, Seo, and Kikuchi [START_REF] Lee | Topology optimization of switched reluctance motors for the desired torque profile[END_REF] and Choi et al. [START_REF] Choi | Topology optimization of the stator for minimizing cogging torque of IPM motors[END_REF] optimize the stator teeth without changing the conductor distribution, and Labbé and Dehez [START_REF] Labbé | Convexity-Oriented Mapping Method for the Topology Optimization of Electromagnetic Devices Composed of Iron and Coils[END_REF] includes only one phase of a reluctant machine within the optimization.

In order to consider many materials, the traditional density method can be extended. One solution is to interpolate M materials using M -1 auxiliary variables per mesh element [START_REF] Bendsøe | Topology optimization: theory, methods, and applications[END_REF]. Other methods require a single variable per element by using peak functions [START_REF] Yin | Topology optimization of compliant mechanisms with multiple materials using a peak function material interpolation scheme[END_REF], or by ordering several interpolation functions in ascending order [START_REF] Zuo | Multi-material topology optimization using ordered SIMP interpolation[END_REF]. Authors have proposed the use of M optimization variables (one per candidate layer) in a scheme denoted as Discrete Material Optimization [DMO, see 40]. More recently, Bruyneel [START_REF] Bruyneel | SFP-a new parameterization based on shape functions for optimal material selection: application to conventional composite plies[END_REF] proposed placing candidate materials at the vertices of a rectangular domain. The design variable is a point inside the domain that interpolates material properties. Thus, M = 4 candidate materials can be handled using two design variables: the Cartesian coordinates of the inner point. This idea, denoted as Shape Function with Penalization (SFP), was then extended to handle M = 8 candidate materials using hexahedral domains and three design variables per element [START_REF] Bruyneel | Extensions of the shape functions with penalization parameterization for composite-ply optimization[END_REF][START_REF] Gao | A bi-value coding parameterization scheme for the discrete optimal orientation design of the composite laminate[END_REF] and has been successfully applied, mainly in the design of fiber-reinforced composites [START_REF] Zhang | Topology optimization of hyperelastic structures with anisotropic fiber reinforcement under large deformations[END_REF].

The literature proposes two possibilities to handle a different number of materials. The majority of works use redundant materials (mostly air) as Sigmund [START_REF] Sigmund | Design of multiphysics actuators using topology optimization -Part II: Two-material structures[END_REF] to fill the vertices of an orthogonal interpolation domain. However, no element indicates that orthogonal domains suit all MMTO problems. Jung and Min [START_REF] Jung | Material Interpolation in Multi-Material Topology Optimization for Magnetic Device Design[END_REF] have shown that the attribution of these redundant materials may affect the optimization results. A comparison of two different materials placements is given in Fig. 1a and Fig. 1b. Other authors, as Choi and Yoo [START_REF] Choi | Simultaneous structural topology optimization of electromagnetic sources and ferromagnetic materials[END_REF] use linear programming to handle a linear constraint per mesh element, as illustrated in Fig. 1c. As the obtained triangular domain is not equilateral as in Fig. 1d, a permutation of the materials' positions changes the optimization problem.

As an alternative, this paper presents a new method for MMTO based on SFP formalism, using polytopes as catalog domains that allow the inclusion of any number of candidate materials without repeating any of them. To interpolate material properties inside the polytope, we use generalized basis functions [START_REF] Wachspress | A Finite Element Rational Basis[END_REF]. Since the proposed method uses non-orthogonal domains, many linear constraints are required to keep the design variables within the polytope. However, since we consider low-dimensional polytopes, we observe that a straightforward Projecting Gradient Descent (PGD) algorithm is enough to tackle the 3-phase stator design problem.

The structure of this article is organized as follows. Section 2 presents the adopted topology optimization approach, the optimization problem to be solved, and the implemented algorithm. Section 3 details the novelties of this work related to the polytope domains for placing material candidates and the Wachspress functions for interpolating material properties. Section 4 presents design results for the stator of a Permanent Magnet Synchronous Machine (PMSM). Finally, conclusions and perspectives are (c) Adapted from [START_REF] Choi | Simultaneous structural topology optimization of electromagnetic sources and ferromagnetic materials[END_REF].

Source

Air Iron given in Section 5.

Topology optimization framework 2.1 Physical problem

Among the wide variety of electromechanical actuators, the 3-phase PMSM is one of the most popular in industrial applications, especially in electrical mobility. It comprises two parts separated by an airgap: a rotating element -the rotor -which contains one or several Permanent Magnets (PM), and a fixed part -the stator -which creates a rotating magnetic field generated by varying electric currents inside coils. For standard electric frequencies, the distribution of the magnetic flux density can be described by magnetostatics, which comes from Maxwell's equations:

∇ × (ν(|b|).∇ × a) = j + ∇ × m, (1)
where a = [a x a y a z ] T is the vector potential ensuring the divergence-free of flux density b = ∇ × a, ν is the magnetic reluctivity, j = [j x j y j z ] T is the current density, and m = [m x m y m z ] T is the PM magnetization. Note that ν depends on |b| in the ferromagnetic materials, as shown in Fig. 2. In a 2D problem within the plane (x, y), the x and y components of a and j vanish, as well as the z component of m. Thus, the magnetostatics equation in 2D becomes:

∇.[ν(|b|)∇a z ] = j z + ∇ R -π 2 m xy , (2)
with R -π

Equation ( 2) can be discretized using, for instance, the Finite Element Method (FEM), which gives the resulting non-linear system: where a, j and m, are the discretized a z , j z and m xy , respectively. K is the stiffness matrix, a the degrees of freedom, s j and s m compose the right-hand side, representing the discretized source terms of the system. The interaction between PM and stator magnetic fluxes creates a torque T . A formula to compute the torque of 2D rotating machines using spatial integration of Maxwell stress tensor was proposed by Arkkio [START_REF] Arkkio | Analysis of Induction Motors Based on the Numerical Solution of the Magnetic Field and Circuit Equations[END_REF]. It can be written as:

K(ν(a))a = s j (j) + s m (m), (4)
T = 2πL S airgap airgap x × (x) ds • e z , ( 5
)
where L is the axial length of the machine, S airgap the surface of the airgap between the stator and the rotor, x the position of the integration point, e z the axial unitary vector, and the magnetic part of the Maxwell stress tensor, which depends on the flux density b:

Σ ij = ν 0 b i b j - |b| 2 δ ij . (6)

Optimization problem

This work adopts the density approach. Namely, the design domain is discretized into N triangular finite elements, and the material properties of each element are defined by interpolating the properties of prescribed materials with n design variables. Since multiple materials are considered for the optimization problem and handled in an SFP approach, n depends on the total number of finite elements (N ) and the dimension of the interpolation domain. By contrast with the original SFP scheme, the interpolation domain is not limited to orthogonal ones and can be any convex polytope. From now on, this polytope is denoted by D for the sake of simplicity. To keep a general notation, we define x e by the vector containing the n design variables assigned to the mesh element e. That is, x e = [x e 1 x e 2 ... x e dim(D) ], where x e i denotes the i-Cartesian coordinate inside the polytope associated with the e-element, as shown in Fig. 3. In this work, only 2D and 3D polytopes are considered: dim(D) = 2 or 3.

The proposed methodology for MMTO is presented for designing a PMSM stator. The objective is to maximize the average torque given by Equation ( 5), considering multiple materials in the stator with different magnetic reluctivities and current densities, detailed in Section 4. The design problem is formulated as follows:

min f (x) = -T (a(x)) s.t.: x e ∈ D , e = 1, ..., N K(ν(a, x))a = s j (j(x)) + s m (m) (7)
Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator.

Structural and Multidisciplinary Optimization volume, 65, 2022. doi: 10.1007/s00158-022-03460-1

x 1

x 2 (a) Orthogonal domain with x the set of all optimization variables, i.e. x = [x e ] e∈ 1,N . Usually, a volume constraint is necessary to avoid the trivial solution, such as a fully solid beam in the compliance minimization problem. However, a fully solid stator cannot produce torque in the case of electrical machines, so a volume constraint is unnecessary and was not used. Moreover, global constraints are non-trivial to implement in the presented approach and thus will be tackled in future works.

Optimization algorithm

The optimization problem ( 7) is solved with a Projected Gradient Descent (PGD). The iterative optimization process consists of four steps illustrated in Fig. 4. The first step is to evaluate the material properties within each finite element from the optimization variables. This procedure will be detailed in Section 3 as it presents the main novelty of this work. The second step is to solve Equation (4) to evaluate the magnetic performance as defined previously (Section 2.1). The third step is to compute the sensitivities using the Adjoint Variable Method (AVM). To do so, the following adjoint problem must be solved [START_REF] Hinton | Analysis and optimization of Prismatic and Axisymetric Shell Structure[END_REF]:

∂(Ka -s j -s m ) ∂a T λ = ∂f ∂a , (8)
then the sensitivities to each component i of x e are computed for each element e with:

df dx e i = ∂f ∂x e i + λ T ∂s j ∂x e i - ∂K ∂x e i a . (9)
Our implementation of the AVM has been validated by comparison with finite differences. The fourth step of the optimization process is to update the design variables using the Projected Gradient Descent (PGD) algorithm. This algorithm consists of three phases, which are detailed below:

1. Compute the descent direction. As sensitivities magnitudes can be very different from one mesh element to another, which leads to slow convergence, each component i of the local descent direction vectors

d e = [d e 1 d e 2 ... d e dim(D)
] associated with a finite element e are normalized as follows: Interpolate material properties with (13) using Wachspress ( 16) & RAMP [START_REF] Tavakoli | Alternating active-phase algorithm for multimaterial topology optimization problems: A 115-line MATLAB implementation[END_REF] Solve FEM (4)

d e i = - df dx e i df dx e . ( 10

Compute sensitivities

with ( 8) and ( 9)

Update design variables • descent direction [START_REF] Takezawa | Design methodology for porous composites with tunable thermal expansion produced by multi-material topology optimization and additive manufacturing[END_REF] • compute step size • projection [START_REF] He | Multi-material topology optimization of piezoelectric composite structures for energy harvesting[END_REF] Report solution counter k is heuristically adapted according to the quality of the linearized model. A quality indicator q is defined as:

no x k κ (x k ) f (x k ) df (x k ) dx x k x k+1
q = f (x k ) -f (x k-1 ) df dx T (x k -x k-1 ) . ( 11
)
If the linearized model is considered as bad (when q is too low), the iteration is rejected and the step size is reduced. If the model is considered good enough, the step size is increased to accelerate the convergence. Algorithm 1 details the whole procedure. The initial step size α 0 is set to 0.1.

Projection.

After applying the descent, some points may lie outside the polytope domain D. In that case, they should be projected back onto it using the min-norm operator denoted as P. The PGD update reads:

x k+1 = P(x k + α k d k ). ( 12
)
Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator.

Structural and Multidisciplinary Optimization volume, 65, 2022. doi: 10.1007/s00158-022-03460-1

Algorithm 1 Control of the step size α

1: Compute q 2: if q < 5 × 10 -3 then
The linear model is bad.

3:

Reject the iteration

4: α ⇐ α/2
The step size is decreased.

5: else if q ∈ [5 × 10 -3 , 5 × 10 -2 ] then 6:
Accept the iteration 7: else if q > 5 × 10 -2 then

The linear model is good.

8:

Accept the iteration

9: α ⇐ 1.2α
The step size is increased. 10: end if 11: Continue the optimization process For orthogonal domains as in Fig. 3a, the application of P is a simple truncation of x-components. For other types of interpolation domains, it implies considering many linear constraints and finding the active one to project the point onto, which can be a facet, an edge, or a vertex of the polytope. This task can demand an extensive computation time on high-dimensional polytopes without efficient projection algorithms (Haddock 2018). Fortunately, in the present case, P can be decomposed in N projections of x e onto a low dimensional convex polytope D (2D or 3D). If x e is located outside D, it belongs to the normal cone associated with the face i, onto which x e should be projected, as shown in Fig. 5.

x 1 P (x 1 ) P (x 2 )
x 2 A vectorized Matlab 2020b implementation of the projection algorithm is available on Github and archived in Zenodo [START_REF] Cherrière | [END_REF].

Interpolation scheme

The expressions of the interpolated material properties are essential in density-based methods to obtain realistic design solutions. This work proposes to extend the SFP formalism introduced by Bruyneel [START_REF] Bruyneel | SFP-a new parameterization based on shape functions for optimal material selection: application to conventional composite plies[END_REF] to a broader class of interpolation domains. In order to build such interpolations, four choices have to be made:

1. The interpolation domain D with M vertices.

2. The materials placement on its vertices.

Note that it is possible to associate a different penalization function for each vertex. However, we consider the same penalization function P for every vertex for simplicity. The derivative of (13) with respect to the i-component of x e is necessary to compute the sensitivities with [START_REF] Kazemi | Multi-material topology optimization of lattice structures using geometry projection[END_REF], and reads:

d κe dx e i (x e , |b e |) = M v=1 κ v (|b e |)P (ω D v (x e )) dω D v (x e ) dx e i . (14)

Interpolation domain

The first step consists in choosing a suitable domain D for the material interpolation, with each prescribed material associated with one of its vertices. In the literature, only hypercubes of [0, 1] dim are chosen.

Hypercubes are orthogonal domains that can contain 2 dim materials. Their main advantage is their boundary constraints, which can be easily handled by a simple truncation algorithm. However, compared to other polytopes, hypercubes may not be suitable for some MMTO problems because, as shown in our numerical examples (Section 4), they promote local optima even when the number of prescribed materials is a power of two. This work considers 2D and 3D convex polytopes as admissible domains. Some examples are presented in Table 1.

Materials placement

The next step consists of associating each material considered in the optimization to a vertex of D.

As shown by [START_REF] Jung | Material Interpolation in Multi-Material Topology Optimization for Magnetic Device Design[END_REF], the results of the optimization depend on the materials placement. An intuitive choice for the 3-phase PMSM stator optimization problem is to place the electric sources on the same plane, according to their phases, in orthogonality with steel. This can be made with a hexagonal-based diamond but not with a classical orthogonal domain such as a cube. The methodology presented in this paper can efficiently handle any convex polytope, which is essential for complex MMTO problems such as a 3-phase PMSM stator optimization as shown by our numerical experiments (Section 4.2).

Basis functions

Once the domain D and the materials placement are chosen, an interpolation can be defined thanks to a set of differentiable basis functions ω D v : D → [0, 1] associated with each vertex v, such as ω D v = 1 on the vertex v and 0 on the others. An interesting property is the partition of unity:

Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator. Structural and Multidisciplinary Optimization volume, 65, 2022. doi: 10.1007/s00158-022-03460-1

∀x e ∈ D, M v=1 ω D v (x e ) = 1, (15)
which allows the interpretation the basis functions as material volume fractions. However, partition of unity is required only at the end of the optimization process, so the penalization of the intermediate values of the basis functions is possible [START_REF] Bruyneel | SFP-a new parameterization based on shape functions for optimal material selection: application to conventional composite plies[END_REF] and will be detailed in Section 3.4. Several different basis functions meet the partition of unity, such as Wachspress, mean-value, or Malsch interpolant, see Sukumar and Malsch [START_REF] Sukumar | Recent advances in the construction of polygonal finite element interpolants[END_REF] and Kraus, Rajagopal, and Steinmann [START_REF] Kraus | Investigations on the polygonal finite element method: Constrained adaptive Delaunay tessellation and conformal interpolants[END_REF]. In this paper, we use [START_REF] Wachspress | A Finite Element Rational Basis[END_REF] basis functions, which read in 2D:

ω v (x e ) = ϕ v (x e ) M j=1 ϕ j (x e ) , with ϕ v (x e ) = A v A + v (x e ).A - v (x e ) , ( 16
)
where A - v (x e ), A + v (x e ), A v represent the areas of the triangles defined, respectively, by vertices

[x e v-1 , x e v ,
x e ], [x e v , x e v+1 , x e ] and [x e v-1 , x e v , x e v+1 ] as shown in Fig. 7a. An example of basis function is plotted in Fig. 7b. In higher dimensions, A v becomes a simplex hypervolume related to the adjacent facets. A code to compute the basis functions values as well as their gradients for any convex 2D and 3D polytopes is given by Floater, Gillette, and Sukumar [START_REF] Floater | Gradient bounds for Wachspress coordinates on polytopes[END_REF]. A vectorized version of this code was used to handle many points [START_REF] Cherrière | [END_REF].

x v-1 x v x v+1 x e A v A ¡ (x e ) v A + (x e ) v (

Penalization

An interpolation built with basis functions matching the partition of unity may not lead to meaningful designs. Some intermediate materials often remain, and a penalization P is necessary to avoid them. A wide variety of such functions exists. The most popular is the power-law from the SIMP method [START_REF] Bendsøe | Optimal shape design as a material distribution problem[END_REF]. In this paper, we chose to use the RAMP [START_REF] Stolpe | An alternative interpolation scheme for minimum compliance topology optimization[END_REF] applied to each basis function:

P : [0, 1] → [0, 1] ω → ω 1 + p(1 -ω) , ( 17
)
where p is the penalization factor. The RAMP was preferred to the classical SIMP scheme because it preserves the non-zero gradient of the basis functions on all the vertices of the polytope, as shown in Fig. 8.

Numerical examples and discussion

The presented methodology can handle numerous materials. As an example, we applied it to the MMTO of a complete three-phase PMSM stator. Although three-phase machines are widely used in industrial applications, only their rotor is generally optimized as in Sato and Igarashi [START_REF] Sato | Deep learning-based surrogate model for fast multi-material topology optimization of IPM motor[END_REF], and sometimes with the stator teeth [START_REF] Lee | Topology optimization of switched reluctance motors for the desired torque profile[END_REF][START_REF] Choi | Topology optimization of the stator for minimizing cogging torque of IPM motors[END_REF]. When the magnetic source distribution has been addressed, only one electric phase had been considered, as in Labbé and Dehez [START_REF] Labbé | Convexity-Oriented Mapping Method for the Topology Optimization of Electromagnetic Devices Composed of Iron and Coils[END_REF]. This numerical example is the MMTO of a 3-phase PMSM, containing all electrical phases, steel, and air. To the authors' knowledge, this problem has not been addressed in the literature, perhaps because of the high number of materials and the lack of appropriate interpolation domains.

Presentation of the problem

We consider Equation ( 7) as the optimization problem. It contains M = 8 materials to distribute within the stator, summarized in Table 2:

• Air, which is a linear material with a magnetic reluctivity of ν 0 = 1 4π × 10 -7 m/H.

• Steel, which has a non-linear behavior. Its reluctivity is given in Fig. 2b and follows the anhysteretic behavior of standard FeSi.

• Six different types of copper conductors, which have the same magnetic behavior as air. They carry current densities j α with the same amplitude but different electrical phases:

j α (θ e ) = J cos(p p θ -απ/3 + ψ optim ), α ∈ 0, 5 (18)
where J is the amplitude of the current density set to 10 A/mm 2 , θ is the mechanical angle varying from 0 • to 30 • , p p = 6 is the number of poles pairs, and ψ optim is the fixed load angle set to 288 • , which corresponds to the load angle where the torque is maximum on the reference design in Fig. 9. The load angle ψ can be interpreted as either a magnetic rotation of the magnetic fields created from the stator without mechanical rotation or a mechanical rotation of the stator without magnetic rotation, as shown in Fig. 10.

The torque is computed with (5) on 60 angular positions of the rotor and averaged. The rotor is inspired by the one of the BMW i3 [START_REF] Staton | Open Source Electric Motor Models for Commercial EV & Hybrid Traction Motors[END_REF] shown in Fig. 9. The problem is discretized on a mesh generated with Gmsh [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF] shown in Fig. 11. As the rotor structure is periodic, only one magnetic pole is considered in the FEM analysis. Therefore, only this pole of the obtained designs will be shown to save space in the rest of the article. The numerical examples are used to assess three main aspects of the proposed methodology: (i) the polytope domain D, (ii) the penalization parameter p, and (iii) the load angle ψ. In the following section, the first aspect is discussed.

0 Steel (FeSi) ν(|b|) -see Fig. 2b 0 Phase A + Phase A + Phase A + Phase A + Phase A + Phase A + Phase A + Phase A + Phase A + Phase A + Phase A + Phase A + Phase A + Phase A + Phase A + Phase A + Phase A + ν 0 J cos(θ e + ψ) Phase A - Phase A - Phase A - Phase A - Phase A - Phase A - Phase A - Phase A - Phase A - Phase A - Phase A - Phase A - Phase A - Phase A - Phase A - Phase A - Phase A - ν 0 J cos(θ e -π + ψ) Phase B + Phase B + Phase B + Phase B + Phase B + Phase B + Phase B + Phase B + Phase B + Phase B + Phase B + Phase B + Phase B + Phase B + Phase B + Phase B + Phase B + ν 0 J cos(θ e -2π/3 + ψ) Phase B - Phase B - Phase B - Phase B - Phase B - Phase B - Phase B - Phase B - Phase B - Phase B - Phase B - Phase B - Phase B - Phase B - Phase B - Phase B - Phase B - ν 0 J cos(θ e -5π/3 + ψ) Phase C + Phase C + Phase C + Phase C + Phase C + Phase C + Phase C + Phase C + Phase C + Phase C + Phase C + Phase C + Phase C + Phase C + Phase C + Phase C + Phase C + ν 0 J cos(θ e -4π/3 + ψ) Phase C - Phase C - Phase C - Phase C - Phase C - Phase C - Phase C - Phase C - Phase C - Phase C - Phase C - Phase C - Phase C - Phase C - Phase C - Phase C - Phase C -

Interpolation domain

The choice of the interpolation domain is essential and beyond the number of materials M . Indeed, it exists several possible polytopes which have M vertices. Since this work considers eight candidate materials, it is interesting to compare the performance of the proposed MMTO approach using different 8-vertex polytopes. The polytopes under analysis are a 2D octagon (Fig. 12a), a 3D cube (Fig. 12b) and a 3D diamond (Fig. 12c). The 3D cube is representative of the literature [START_REF] Bruyneel | Extensions of the shape functions with penalization parameterization for composite-ply optimization[END_REF], but the octagon and the diamond are new domains in MMTO, to the best of our knowledge.

We place the conductors consistently with their electric phases, which is easy for the diamond domain, but not for the others. At the beginning of the optimization process, the design variables are placed at the centroid of the polytopes. All the other optimization parameters are strictly the same. The material properties ν and j are interpolated by using Equation ( 17), with p ν = 4 and p j = 3. These penalization values are justified in Section 4.3. The PGD was applied during 500 iterations, and the final structures are plotted in Fig. 12.

We note that no design contains air, which is logical because air has the same magnetic behavior as conductors but carries no current. Therefore, they cannot increase the torque more than conductors. Visually, the best-defined design is given by the diamond domain in Fig. 12c. For that design, the Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator. optimized design variables are displayed in Fig. 13 as points inside the diamond polytope. There, it can be observed that most design variables have converged to the vertices of the polytope. Notably, a large number of variables that point to iron lie in the vicinity of the vertex. Indeed, at high flux densities, intermediate materials may carry more flux than iron because of penalization, which is analyzed in detail in Section 4.3.

Dirichlet Sliding band Antiperiodic

Regarding the optimized design, one can identify small design features where portions of materials appear isolated from the main bulk. Fig. 14 decomposes the topology by material to facilitate the visualization of these zones. According to our numerical experiments, these small design features are caused by a poor convergence of the optimization problem to a local optimum and not by numerical instabilities such as a checkerboard. This can be seen in Figure 15, where two results obtained with the same set of parameters but with different discretization (coarse and refined) are presented. Both topologies are the same, i.e., the global shape and arrangement of the materials are the same. Similar to Choi and Yoo [START_REF] Choi | Simultaneous structural topology optimization of electromagnetic sources and ferromagnetic materials[END_REF], neither checkerboard nor mesh-dependence was found, so no filtering process was used for the rest of the numerical experiments. This could be the subject of future work, such as when incorporating manufacturing constraints.

The optimized design in Fig. 12c contains three conductor slots, while the reference design contains six slots. The topology change between the reference and the optimized structures demonstrates the usefulness of a MMTO, as the optimized design is unreachable with parametric or shape optimization approaches. However, the numerical results given in Table 3 indicates that the cube (Fig. 12b) returns a structure associated with higher torque. These obtained structures are then fed with currents with Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator. different load angles ψ, as for a real machine where the optimal feeding is determined from its design. The resulting torques plotted in Fig. 17 show clearly that the diamond gives a structure that can reach higher torque than the two other designs, and even higher than the reference machine. Therefore, the structure obtained with a diamond-based catalog is better than the two others, as ψ represents only a rotation between the same stator and the rotor. Moreover, it seems closer to convergence than the other domains as the relative variation of the objective function is 1000 times lower, as shown in the 4 th row of Table 3. Note that when the step size becomes too big, some iterations may be rejected by Algorithm 1, so Table 3 also gives the number of accepted iterations. For the rest of the numerical experiments, we thus choose to use the diamond domain. Note that a diamond-type domain can be used with this approach to address stator configurations that have a number of materials other than 8, as shown in Fig. 16.

-f end-1 | f end 2.0 × 10 -5 2.3 × 10 -5 3.4 × 10 -8
Nevertheless, Fig. 17 shows that the load angle where the torque is maximum ψ max is shifted from the angle ψ optim imposed during the optimization process. It means that a simple rotation applied to paper, we use the RAMP interpolation scheme proposed on Equation ( 17) directly on ν. While the obtained design does not seem sensitive to p j , the ν-interpolation strongly influences the result, which should be investigated. Its value is tested between 0 (which corresponds to a linear interpolation) and 10. We found that for low p ν values, the final designs tend to have very little steel, as shown in Fig 18b.

The steel quantity in the stator teeth and the average torque increase with p ν , until the steel is replaced by gray material for a high p ν value, as shown in Fig. 18. The presence of intermediate materials in converged design can be explained because they carry a higher flux density |b| than the steel for high magnetic field |h|. Indeed, the physical limit of a ferromagnetic material's reluctivity is:

lim |b|→∞ ν(|b|) = ν 0 , (19)
while the limit reluctivity ν of the interpolated materials with a strictly convex penalization scheme such as [START_REF] Tavakoli | Alternating active-phase algorithm for multimaterial topology optimization problems: A 115-line MATLAB implementation[END_REF] with p ν > 0 is:

lim |b|→∞ ν(x e , |b|) = ν 0 M v=1 P (ω D v (x e )) < ν 0 . (20)
Therefore, intermediate materials can carry more flux than steel and are better magnetic materials than steel for high flux densities |b|, which justifies their existence in the converged designs. In practice, the larger p, the lower the |b| limit under which the steel is better than intermediate materials. To quantify the presence of intermediate materials, we use the following discreteness indicator in the optimization zone: with ω v the Wachspress basis function verifying the partition of unity, V e the surface of the e-mesh element, and P a convex penalization function. Equation ( 21) is 0 when all the x e are in the center x 0 of D, and 1 if they are all on its vertices (prescribed materials). We arbitrary chose Equation ( 17) with p = 2 for P in [START_REF] Wang | Topology optimization of electromagnetic systems considering magnetization direction[END_REF]. A compromise must be made to get a realistic design with a high torque, as shown in Fig. 18a. The relative variations of the objective function becomes lower than 10 -6 for all designs, which indicates gray materials exists in converged results for high p ν . We chose p ν = 4 and p j = 3 for the remaining numerical tests, which gives good performances with a discreteness around 0.99.

D(x) = N e=1 V e M v=1 P (ω v (x e )) -P (ω v (x 0 )) N e=1 V e 1 - M v=1 P (ω v (x 0 )) , (21)

Discussion on the local optima issue

As shown in Fig. 17, the average torque is not maximum at ψ = ψ optim = 288 • . The situation is illustrated in Fig. 19. The optimum load angle shift means that we can rotate the optimized structure, shown in Fig. 19a to obtain a better design, shown in Fig. 19b under the same electric conditions. Therefore the optimized design is necessarily a local optimum. At least two causes may lead to local optima: the PM induced a field in the initial situation, which

Influence of the remanence

Only two different components can produce a non-zero average torque. The hybrid component is due to the interaction between the winding flux and the magnet flux, while the reluctant component is due to the non-uniformity of the magnetic flux paths. The initial situation is a uniform material with no current. The only magnetic flux comes from the permanent magnets, which may promote too much the hybrid torque against the reluctant one at the beginning of the optimization process. This may lead to the shift between ψ optim and ψ max shown in Fig. 17.

This influence can be verified numerically by removing the permanent magnet remanence B r from the optimization, so the concurrence between the two torque components vanishes. In practice, it is impossible to start the optimization with no source because the physical field a is null, and so are the sensitivities -see Equation [START_REF] Kazemi | Multi-material topology optimization of lattice structures using geometry projection[END_REF]. Therefore, we started with B r = 1 µT. The influence of B r on the final design is plotted in Fig. 20. We note that the bigger the remanence, the bigger the shift between ψ max and ψ optim , which highlights the role of B r in getting local optima.

Influence of the penalization

Another difficulty for avoiding local optima is the penalization. It is known that for simple problems, such as compliance minimization, a penalized material interpolation function breaks the convexity of the problem [START_REF] Abdelhamid | Revisiting non-convexity in topology optimization of compliance minimization problems[END_REF].

To illustrate the non-convexity induced by the penalization, we use the following procedure. The penalization parameters of the magnetic reluctivity (p ν ) and the current density (p j ) are initialized as p ν = 4 and p j = 3. After 20 iterations, the design is almost converged and similar to Fig. 12c, which is a local optimum, so the penalization is removed by setting p ν = 0 and p j = 0. Next, 10 more iterations are carried out, then the penalization of the current density is reactivated as p j = 3 to remove intermediate materials, and the optimization problem is run until convergence. The average torque and some intermediate designs during optimization are plotted in Fig. 22 and 23, respectively. Fig. 23 shows that during the optimization process, materials that were committed for a specific penalization can be decommitted and associated to another material when the penalization is changed. This method was tested for several n ψ . The results are given in Fig. 25 and the structures are shown in Fig. 26. The initial ψ value is ψ 0 = 288 • . As expected, the more often ψ 0 is adjusted, the higher the torque. Moreover, the maximum torque with ψ adjustment is significantly higher than the maximum torque without, even with a high n ψ value.

Algorithm 2 Control of ψ

Require: k, x k , n ψ k is the number of the current iteration of the optimizer.

1: if modulo(n ψ , k) = 0 then 2: for each ψ i ∈ 0 • , 359 • do 3: Compute T i (x k ) with ψ = ψ i 4:
Store T i (x k ) and ψ i 5:

end for 6:

Find the max. torque and associated ψ max 7:

ψ ⇐ ψ max 8: end if 9: Continue the optimization Concerning the load angle, ψ max drifts more from the initial ψ when n ψ is small. Concerning the obtained structures, the final designs become more symmetric with frequent ψ adjustments, which indicates that the asymmetry is also a local optimum. The existence of this local optimum may be due to the PM flux as Fig. 21 seems to show a correlation between symmetry and remanence.

Moreover, low n ψ gives "cleaner" designs, as all elements are the same conductor in a same slot for n ψ = 1 and n ψ = 5: the number of isolated islands of materials decreases with the frequency of the load angle adjustment. Therefore, we recommend adjusting ψ during the optimization process to avoid local optima. The computing time in Algorithm 2 can be reduced by considering a smaller angle interval Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator.

Conclusion

This paper presented a general density-based topology optimization methodology designed to support many different materials. It extends the SFP scheme to non-orthogonal interpolation domains using Wachspress basis functions with a projected gradient algorithm. We applied this methodology to the challenging 8-materials optimization of a 3-phase PMSM stator. Performing and meaningful designs were obtained, demonstrating this methodology's potential. In addition, we highlighted the importance and difficulties of choosing both a suitable interpolation domain and an appropriate material penalization to obtain satisfactory designs through numerical examples. We show that classic orthogonal domains are not suitable for all MMTO problems. In fact, a diamond domain appears to be better suited for this application and may also apply to a different number of electric phases. However, with a load-angle analysis, we found that this methodology could return local optima, as in all gradient-based optimization methods. In this application, they are caused at least by the presence of an external field source -identified as the permanent magnets -and the penalization of the material properties. Adjusting the load angle during optimization solves this issue and leads to better designs, transitioning from asymmetric to almost symmetric structures. This methodology can be improved by adding global constraints to address industrial problems. For instance, imposing the conductors' volume associated with a given current density is directly related to controlling the total electric current. A solution is to use popular optimizers such as MMA or GCMMA [START_REF] Svanberg | A Class of Globally Convergent Optimization Methods Based on Conservative Convex Separable Approximations[END_REF]. It requires including the projection algorithm within the optimizer, which will be done in future work. This methodology could then be applied to other physics and objective functions. As this work extends the SFP scheme designed for composite materials in mechanical engineering, it could be used to optimize different fiber orientations, similar to the different electric phases within a PMSM.

Figure 1 :

 1 Figure 1: Different possible materials placements for a C-core electromagnet optimization.

  and m xy = [m x m y ] T . The 2D flux density reads: b = R -π 2 ∇a z .

Figure 2 :

 2 Figure 2: FeSi anhysteretic behavior: h is the magnetic field, b the flux density and ν = |h| |b| the reluctivity.

Figure 3 :

 3 Figure 3: Examples of interpolation domains. Black dots represent prescribed materials, gray dots represent interpolated materials.

4 Figure 4 :

 44 Figure 4: Flowchart of the optimization algorithm.

Figure 5 :

 5 Figure 5: Min-norm projection onto a pentagon: x 1 is associated to an edge and x 2 to a vertex. It follows a simple procedure to compute P(x e ): a) Compute the normal fan of the polytope -defined as the collection of all the normal cones from each of its faces, as shown in Fig. 6. This needs to be done once, at the beginning of the optimization process. b) Find the normal cone where x e is located. The closest face is associated with this cone. It requires testing each cone. c) Project the point x e onto the closest face with the standard orthogonal projector.

  Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator. Structural and Multidisciplinary Optimization volume, 65, 2022. doi: 10.1007/s00158-022-03460-1 (a) Pentagon (b) Tetrahedron (c) Diamond

Figure 6 :

 6 Figure 6: Normal fans of different polytopes. Red cones are associated with vertices, green cones with edges and blue cones with facets.

3 .

 3 The basis functions set {ω D v } v∈ 1,M . Each basis function is associated with a vertex v of D. 4. The penalization function P . After that, the material property κe of a mesh element e filled with an intermediate material can be constructed from the properties of the prescribed materials κ v , which possibly depends on the local flux density b e : κe (x e , b e ) = M v=1 κ v (|b e |) • P (ω D v (x e )).

Théodore Cherrière , 1 dim

 Cherrière1 Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator. Structural and Multidisciplinary Optimization volume, 65, 2022. doi: 10.1007/s00158-022-03460-

Figure 7 :

 7 Figure 7: Illustration of a Wachspress basis function associated with a vertex of a pentagon.

Théodore Cherrière ,

 Cherrière Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator. Structural and Multidisciplinary Optimization volume, 65, 2022. doi: 10.1007/s00158-022-03460-1 (a) SIMP -p = 3. The gradient is zero in the circled zone. (b) RAMP (17) -p = 4.The gradient does not vanish anywhere.

Figure 8 :

 8 Figure 8: Examples of penalized Wachspress functions.

ν 0 JFigure 9 :

 09 Figure 9: Reference design adapted from the BMW i3 [57].

1 ψ=0(

 1 Structural and Multidisciplinary Optimization volume, 65, 2022. doi: 10.1007/s00158-022-03460ψ seen as a mechanical angle.

Figure 10 :

 10 Figure 10: Equivalence between load angle ψ interpretations. Arrows represent magnetic field reference.The rotor is represented in orange and the stator in blue.

Figure 11 :

 11 Figure 11: Mesh with boundary conditions obtained (11331 triangles, 5890 nodes).

  Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator. Structural and Multidisciplinary Optimization volume, 65, 2022. doi: 10.1007/s00158-022-03460-1

Figure 12 :

 12 Figure 12: Different interpolation domains associated with their final design.

Figure 13 :

 13 Figure 13: Repartition of the optimization variables for the design plotted in Fig. 12c. Each point represents the material of a mesh element from the stator.

Figure 15 :

 15 Figure 15: Optimized stator designs with a diamond domain and different meshes.

Théodore Cherrière ,Figure 16 :igure 17 :

 Cherrière1617 Figure 16: Examples of diamond-type domains as in Fig. 12c for other multiphased stator. The numbers indicate the electric phases.

  Final torque and discreteness (21) with pν . (b) pν = 0 (c) pν = 4 (d) pν = 8

Figure 18 :

 18 Figure 18: Effect of p ν on the converged designs.

  Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator. Structural and Multidisciplinary Optimization volume, 65, 2022. doi: 10.1007/s00158-022-03460-1 (a) Obtained design with a load angle shift (local optimum). (b) Potential global optimum (no load angle shift).

Figure 19 :

 19 Figure 19: Illustration of the effect of the load angle shift on the design.

Figure 21 :

 21 Figure 21: Final designs obtained for different remanences.

  Fig.23b

Fig.23e

  Fig.23e

Fig.23f

  Fig.23f

Fig.23a

  Fig.23a

Fig.23c

  Fig.23c

Fig.23d

  Fig.23d

Figure 22 :

 22 Figure 22: Average torque during topology optimization with penalization changes. The prescribed penalization is shown in red.

  Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator. Structural and Multidisciplinary Optimization volume, 65, 2022. doi: 10.1007/s00158-022-03460-1 (a) Initial situation (b) it. 8 (c) it. 21 (d) it. 25 (e) it. 31 (f) Final design

Figure 23 :igure 24 :

 2324 Figure 23: Design snapshots during the optimization process with penalization changes. The average torque for each figure is illustrated in Fig. 22.

  Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator. Structural and Multidisciplinary Optimization volume, 65, 2022. doi: 10.1007/s00158-022-03460-1

	Air	(0,1)	(1,1)	Source	Iron	(0,1)	(1,1)	Air
	Air	(0,0)	(1,0)	Iron		Air	(0,0)	(1,0)	Source
	(a) RM placement adapted	(b) DMO placement adapted
		from [46].				from [46].
	Source	(0,1)	(1,1)
		Air	(0,0)	(1,0)	Iron

  The step size at iteration k is denoted as α k and is adjusted by a simplified trust-region algorithm inspired from Nocedal and Wright[START_REF] Nocedal | Numerical optimization[END_REF]. The step Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator. Structural and Multidisciplinary Optimization volume, 65, 2022. doi: 10.1007/s00158-022-03460-1

	Discretization
	Initialize
	design variables
	Convergence?

) 2. Estimate the appropriate step size.

Table 1 :

 1 Examples of admissible interpolation domains.

  Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator. Structural and Multidisciplinary Optimization volume, 65, 2022. doi: 10.1007/s00158-022-03460-1

Table 2 :

 2 Materials used in the optimization

	Material	Reluctivity	Current density
	Air	ν 0

Table 3 :

 3 Optimization results with different domains.

	Domains	Octagon	Cube	Diamond
	Torque (Nm/m)	1703	1836	1720
	Total iter.	500	500	500
	Accepted iter.	401	400	395
	|f end

  Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator.

	Structural and Multidisciplinary Optimization volume, 65, 2022. doi: 10.1007/s00158-022-03460-1
	0	2	4	6	8	10

Penalization parameter

Once the domain and materials placement are chosen, which are global and apply to all physical properties, the choice of the penalization is tricky. There are two material coefficients to interpolate in the partial differential equation (2): the magnetic reluctivity ν and the current density j. The associated penalization coefficients are p ν and p j , respectively. Their role is to avoid intermediate materials to have interesting properties simultaneously (i.e. high current density and good magnetic property). In this Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator. It can also be observed that such a scheme presents in its final design a smaller number of isolated materials (material zones that are represented by only one element) compared to the designs obtained with a fixed penalization.

Regarding the performance, Fig. 23 shows that the final design artificially changes the magnetic phase by mixing some electrical phases within a single slot instead of mechanically rotating the whole structure. By doing so, the load angle shift is highly reduced, as shown in Fig. 24, and the torque is increased. This result confirms that the penalization leads to a local optimum shown in Fig. 23c, and the modification of the penalization parameters allows to get better designs such as the structure shown in Fig. 23f.

Parametric ψ adjustement

The shift between the optimization load angle ψ optim and the maximum torque load angle ψ max is an indicator of local optimum. A simple way to avoid this problem is to control ψ optim to ensure that it is always equal to ψ max with a parametric study. To save computation time, we can update ψ all the n ψ iterations with Algorithm 2. n ψ = 1 means ψ is ajusted every iterations, n ψ = ∞ means it is never adjusted.

Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi & Eduardo Fernández. Multi-material topology optimization using Wachspress interpolations for designing a 3-phase electrical machine stator.

Conflict of interest

On behalf of all authors, the corresponding author states that there is no conflict of interest.

Replication of results

Our work relies on the Matlab R2020b programming language. An implementation of the projection algorithm presented in Section 2 and the Wachspress interpolants computation presented in Section 3 are freely available on GitHub [START_REF] Cherrière | [END_REF]. All other information required to replicate the results of this work is clearly detailed in the manuscript. The model, as well as the optimization algorithm, are given in Section 2. Interested readers may send a request to the corresponding author to obtain the .mat files of results.


Jiajia Dong
email: jdong2@mnhn.fr

Natallia Vicente

Ioana C Chintauan-Marquier

Cahyo Ramadi

Agnès Dettai

Tony Robillard

Complete mitochondrial genome and taxonomic revision of Cardiodactylus muiri Otte, 2007 (Gryllidae: Eneopterinae: Lebinthini)

Keywords: mitogenome, long-PCR approach, next generation sequencing technology, New Guinea, Indonesia

In the present study, we report the high-coverage complete mitochondrial genome (mitogenome) of the cricket Cardio-dactylus muiri Otte, 2007. The mitogenome was sequenced using a long-PCR approach on an Ion Torrent Personal Ge-nome Machine (PGM) for next generation sequencing technology. The total length of the amplified mitogenome is 16,328 bp, representing 13 protein-coding genes, 22 transfer RNA genes, two ribosomal RNA genes and one noncoding region (D-loop region). The new sets of long-PCR primers reported here are invaluable resources for future comparative evolu-tionary genomic studies in Orthopteran insects. The new mitogenome sequence is compared with published cricket mito-genomes. In the taxonomic part, we present new records for the species and describe life-history traits, habitat and male calling song of the species; based on observation of new material, the species Cardiodactylus buru Gorochov & Robillard, 2014 is synonymized under C. muiri.

Introduction

In the era of next generation sequencing technologies, complete mitochondrial genomes (mitogenomes) are becoming a common and valuable resource for phylogenetic studies and molecular systematics [START_REF] Crampton-Platt | Soup to Tree: The Phylogeny of Beetles Inferred by Mitochondrial Metagenomics of a Bornean Rainforest Sample[END_REF], 2016;[START_REF] Gomez-Rodriguez | Validating the power of mitochondrial metagenomics for community ecology and phylogenetics of complex assemblages[END_REF]. In a recent phylogeny of the diverse and ancient insect order Orthoptera, [START_REF] Song | 300 million years of diversification: elucidating the patterns of orthopteran evolution based on comprehensive taxon and gene sampling[END_REF] used both mitogenomes and traditional molecular markers as sources of phylogenetic information. They showed that despite the long evolutionary history of these insects, their mitogenome has remained relatively stable, making it a powerful tool to resolve phylogenetic relationships. However, mitogenome sequences have been little investigated in Orthopteran insects compared with other animal groups [START_REF] Zhou | The phylogeny of Orthoptera inferred from mtDNA and description of Elimaea cheni (Tettigoniidae: Phaneropterinae) mitogenome[END_REF], and large areas of their phylogenetic tree are barely or not documented when considering these data. This is particularly true in the superfamily Grylloidea. Despite the diversity of crickets and their allies, only a handful of mitogenomes are available today in databases (four in [START_REF] Song | 300 million years of diversification: elucidating the patterns of orthopteran evolution based on comprehensive taxon and gene sampling[END_REF]. The families Gryllotalpidae, Myrmecophilidae and Schizodactylidae are each documented by one species, while the "true crickets", comprising three large families, namely Gryllidae, Phalangopsidae and Trigonidiidae [START_REF] Chintauan-Marquier | Laying the foundations of evolutionary and systematic studies in crickets (Insecta, Orthoptera): a multilocus phylogenetic analysis[END_REF], are represented only by five species belonging to a single subfamily (Gryllinae): Teleogryllus commodus [START_REF] Wolff | Selective enrichment and sequencing of whole mitochondrial genomes in the presence of nuclear encoded mitochondrial pseudogenes (numts)[END_REF], Teleogryllus oceanicus [START_REF] Zhou | The complete mitochondrial genome of the black field cricket, Teleogryllus oceanicus[END_REF], Loxoblemmus equestris, Teleogryllus emma and Velarifictorus hemelytrus (Yang et al. 2016a). The direct consequence of this is that the phylogenetic information of mitogenomes between and within cricket families has never been evaluated.

In this study, we present a new set of primers to obtain mitogenomes from long-range PCRs in next generation sequencing (NGS) (e.g., [START_REF] Jia | Long-range PCR in next-generation sequencing: comparison of six enzymes and evaluation on the MiSeq sequencer[END_REF][START_REF] Hinsinger | Fishing for barcodes in the Torrent: from COI to complete mitogenomes on NGS platforms[END_REF] and report the complete mitogenome of Cardiodactylus muiri Otte, 2007. This species belongs to the subfamily Eneopterinae, which is characterized by diverse shapes and original features related to acoustic communication (e.g., Robillard & Desutter-Grandcolas 2004a, b;[START_REF] Robillard | Mechanisms of high-frequency song generation in brachypterous crickets and the role of ghost frequencies[END_REF][START_REF] Ter Hofstede | Evolution of a Communication System by Sensory Exploitation of Startle Behavior[END_REF]. In the taxonomic part of the study, we revise C. muiri based on recent collections in New Guinea (Papua province, Indonesia) and C. buru Gorochov & Robillard, 2014 is synonymized under C. muiri. We also present new geographic records for the species and describe life-history traits, habitats and male calling song.

Materials and methods

Materials. The new material was collected in the "Lengguru 2014" expedition in Papua, Indonesia and will be deposited in the collections of the Muséum national d'Histoire naturelle Paris, France (MNHN) and in the Museum Zoologicum Bogoriense, Bogor, Indonesia (MZB). The crickets were collected mostly by sight, night and day, in order to observe the habitat and activity of the species. Square brackets in the list of specimens are used for additional information not mentioned on the specimen labels or translated from French labels.

Molecular techniques. Total genomic DNA was extracted from the median femoral muscles using the NucleoSpin ® 96 Tissue Kit (Macherey-Nagel, Germany) and the automatic pipetting robot Eppendorf epMotion 5075 TMX following the manufacturer's manual. The complete mitochondrial genome sequence of C. muiri was obtained by overlapping three DNA fragments (> 5.5kb) amplified by long-range PCRs and the HotStart LongAmp ® Taq DNA polymerase (New England BioLabs ® Inc., Ipswich MA). The primers used for these PCRs are presented in Table 1. The amplicons were sequenced using the Ion Torrent Personal Genome Machine (PGM) sequencing platform with two level multiplexing [START_REF] Hinsinger | Fishing for barcodes in the Torrent: from COI to complete mitogenomes on NGS platforms[END_REF] at the "Service de Systématique Moléculaire" of the MNHN. The reads obtained from each long-range PCR were separately trimmed by error probability algorithm and assembled using the mapping to reference algorithm in GENEIOUS (8.1.7 version) (Table 1) [START_REF] Kearse | Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data[END_REF]. The references used for mapping each long-range PCR were the partial mitochondrial gene markers of C. muiri obtained from a phylogenetic study of Cardiodactylus in preparation molecular data library (Dong et al. in prep., see details in Table 1). After this first assembly, the three contigs obtained were re-assembled into a circular mitochondrial genome with high coverage by the de novo assembly algorithm in GENEIOUS. Sequence analysis. Annotation of the mitogenome was performed using the MITO webserver with invertebrate genetic code [START_REF] Bernt | MITOS: Improved de novo metazoan mitochondrial genome annotation[END_REF] and modified after comparisons with other mitogenomes from Grylloidea species. The validation of tRNA sequences was performed in tRNAscan-SE (http://lowelab.ucsc.edu/ tRNAscanSE/) using the invertebrate mitogenome genetic codon [START_REF] Schattner | The tRNAscan-SE, snoscan and snoGPS web servers for the detection of tRNAs and snoRNAs[END_REF]. The nucleotide base composition of the complete mitogenome was calculated in GENEIOUS.

Phylogenetic analysis. To validate the reliability of the mitogenome sequence and its identification, the cytb gene and rrnS sequences of this mitogenome were extracted and compared with other species of Cardiodactylus. In addition, a phylogenetic analysis of the family Gryllidae based on the sequences of cytb and rrnS was conducted using a Maximum Likelihood approach. The phylogenetic trees were reconstructed on the IQ-TREE web server [START_REF] Nguyen | IQ-TREE: A Fast and Effective Stochastic Algorithm for Estimating Maximum-Likelihood Phylogenies[END_REF].

The whole mitogenomes nucleotide base composition of C. muiri was calculated in GENEIOUS. The Loxoblemmus equestris (KU562919), Teleogryllus emma (KU562917) and Velarifictorus hemelytrus (KU562918)). Some of these sequences however include nuclear mitochondrial DNA sequences (numts), as reported from [START_REF] Wolff | Selective enrichment and sequencing of whole mitochondrial genomes in the presence of nuclear encoded mitochondrial pseudogenes (numts)[END_REF] and Yang (2016) and were thus discarded here. We only considered the sequences of L. equestris (KU562919), T. emma (KU562917) V. hemelytrus (KU562918) and C. muiri which was obtained in this study. The nucleotide base composition comparison analysis for the whole mitogenomes, two ribosomal genes (the large and small ribosomal subunit, rrnL and rrnS) and the A+Trich regions in grylline crickets were also calculated.

The percentage of pairwise identical bases (% identity) for the whole mitogenomes, each protein coding genes (PCGs) and two ribosomal genes (rrnL and rrnS) of grylline crickets were calculated in GENEIOUS. The average percentage of identical bases were then calculated between two subfamilies (Eneopterinae and Gryllinae) and within Gryllinae.

Besides, the comparison analysis of synonymous and non-synonymous substitutions (Ka/Ks) for all PCGs in grylline crickets were calculated using MEGA 7.0 [START_REF] Kumar | MEGA7: Molecular Evolutionary Genetics Analysis Version 7.0 for Bigger Datasets[END_REF].

Acoustic analysis. The basic cricket song terminology follows [START_REF] Ragge | The songs of the grasshoppers and crickets of Western Europe[END_REF]. One song unit is called a syllable and corresponds to one opening-closure cycle of the male forewings. Three males of C. muiri have been recorded in the field, in their natural habitat or in semi-captivity in a net cage hanged to a low tree. The recordings were made with a modified Condenser Microphone Capsule CM16 (Avisoft Bioacoustics, Berlin), with a flat frequency response from 3 to 150 kHz (R. Specht pers. comm.), connected to a Fostex FR-2LE digital field recorder (96 kilo-samples per second sampling frequency, 16 bit). Acoustic analyses were performed using the computer software Avisoft-SASLab Pro version 4.40 [START_REF] Specht | Avisoft-SASLab Pro: Sound Analysis and Synthesis Laboratory[END_REF]. All recording files are deposited in the Sound Library of the MNHN.

Taxonomic part

Subfamily Eneopterinae Saussure, 1874 Tribe Lebinthini Robillard, 2004 Genus Cardiodactylus Saussure, 1878 Life history traits and habitats. C. muiri lives in dense forested habitats on tree trunks of various sizes or on dead trees still standing (Fig. 1). Groups of singing males (three or more) distant by ca. 50 cm to one meter, are commonly found during afternoon on large tree trunks above three meters high (Fig. 2). They sit on the bark surface and sometimes in crevices of the bark or of epiphytes. Songs were not heard during night, even when habitats with singing males heard during the day were visited after sunset. It could indicate that the species is diurnal or crepuscular and that the formed couples hide at night for mating. However, the species seems to have a low-light preference, given that light intensity under forest cover is limited during afternoon and in the confined habitats (Fig. 1). Females are found ovipositing during early night in the leaf litter or in roots near the ground (Fig. 3).

Calling song. The calling song of C. muiri (Fig. 4) consists of only one long syllable. The song bouts are rather irregular, which disqualifies them as echemes. At 26 °C (MNHN-EO-ENSIF4405, measurement of 28 syllables) the call duration is 38 ± 11.4 ms (mean ± SD), with a period of 4.78 s ± 3.83 s, giving a syllable duty cycle of 8%. The dominant frequency is 11.98 ± 0.24 kHz, which is rather low for the genus and corresponds to the third peak of the frequency spectrum, the two first peaks being little marked.

Results of sequencing and phylogenetic analyses

The mitogenome is 16,328 bp long (Fig. 5), organized in the typical set of 37 genes (13 protein-coding genes, 22 transfer RNA genes and two ribosomal RNA genes) and one A+T-rich control region including two 279 bp repeated regions (Table 2). The ML gene trees, based on cytb and rrnS markers, are congruent and showed that the sequences issued from the present study clustered with published sequences of Cardiodactylus species and other representatives of the subfamily Eneopterinae (Fig. 6).

The nucleotide composition of Cardiodactylus muiri is typically insect A+T biased, with 76.5%, and is slightly A skewed (AT-skew=0.08) and strongly C skewed (GC-skew= -0.29). The nucleotide compositional skew of the whole mitogenomes, rrnL genes, rrnS genes and the A+T-rich regions of Gryllinae crickets are also represented (Table 3, Fig. 7).

The conventional start/stop codons and incomplete stop codons are represented in the PCGs (see details in Table 2). All the protein-coding genes (PCGs) in the mitogenome start with a typical ATN codon but end with different stop codons. The complete stop codon TAA are determined in genes cox2, atp8, atp6, nd3, nd4l, nd6 and cytb, while the incomplete stop codons (T) are found in genes nd2, cox1, cox3, nd5, nd4, nd1. The lengths of 22 tRNA genes range from 62 to 71 bp. Two rRNA genes (rrnL and rrnS), separated by tRNA-Val, are located between tRNA-Leu and A+T-rich region. The lengths of rrnL and rrnS are 1,362 and 771 bp, respectively. Their A+T contents are 78.5% and 74.9%, respectively. The A+T-rich region of 1,472 bp, with an A+T content of 80.1%, is located between rrnS and trnI. The percentage of pairwise identical bases (% identity) for the whole mitogenomes, each PCG and two ribosomal genes, and the percentage of identity between the two subfamilies (Eneopterinae and Gryllinae) and within Gryllinae are presented in Table S1 and Fig. 8. Except for gene nd3, the percentage identities and the pairwise genetic distance within Gryllinae were higher than between the two subfamilies. Pairwise comparison of Ka/Ks for each PCG are also presented in Table S2 and Fig. 9. The overall Ka/Ks values are less than 1.0 between the two subfamilies, and higher than 1.0 in genes nd1, cox1, cox2, nd3 and cytb within the subfamily Gryllinae. The values for Cytochrome Oxidase complex (cox1, cox2 and cox3) and cytb are less than 1.0, but higher than for other genes between the two subfamilies. The smallest value between the two subfamilies is that for the atp8. The variability of Ka/Ks ratio within Gryllinae is similar with that between the two subfamilies.

Discussion

We generated the first mitogenome for the cricket subfamily Eneopterinae, which allows comparisons with existing mitogenomes from the subfamily Gryllinae. The gene order and orientation within the mitogenome of C.muiri are identical with those described for grylline crickets (Loxoblemmus equestris, Teleogryllus emma and Velarifictorus hemelytrus). The length of the whole mitogenome and the length range for tRNAs in C. muiri are also very close to Gryllinae (Yang et al. 2016b). The whole mitogenomes of these four crickets are weakly AT-skewed (0.08 to 0.10) and strongly GC-skewed (-0.35 to -0.27). The average AT-and GC-skew in the whole mitogenome of the Gryllinae are close to those of C. muiri. However, nucleotide compositional skew in rrnL, rrnS and A+T-rich region are stronger in Gryllinae than in C. muiri, with the exception of the GCskew in A+T-rich region. It is logical that the nucleotide compositional skew in different genes of different subfamilies are slightly divergent, but more mitogenomes of other Eneopterinae crickets are necessary to test whether this tendency is consistent between the subfamilies.

There is no significant difference of the percentage identity, between the two subfamilies and within Gryllinae, for the mitogenome and the genes considered separately, except for the atp8 and nd6. The percentage of pairwise identical bases for the mitogenome and for the individual genes are higher than 70%, both between the two subfamilies and within Gryllinae. The level of identity between the two subfamilies is higher or almost equal to that within Gryllinae in the genes nd3 and nd6, respectively. According to the Ka/Ks ratio, the evolutionary rate of Cytochrome Oxidase complex (cox1, cox2 and cox3) and cytb genes is high between the subfamilies and within Gryllinae. There is an obvious evolutionary rate discrepancy in NADH dehydrogenase complex genes and ATP synthase, since the genes nd2 and nd6 are less divergent than the other NADH dehydrogenase complex genes. In addition, the rate of evolution of nd6 is almost the same, both between the two subfamilies and within Gryllinae. The gene atp8 is found to have the slowest evolution rate between the two subfamilies, while the gene cox2 has the highest rate of evolution within Gryllinae. Thus, the present study revealed these two genes as potential markers for future phylogenetic study at family and subfamily level in crickets. Our results lead to the conclusion that the mitogenomes in crickets have a relatively stable structural and functional evolution, as suggested by previous studies made in Orthoptera [START_REF] Song | 300 million years of diversification: elucidating the patterns of orthopteran evolution based on comprehensive taxon and gene sampling[END_REF] and in insects in general (e.g. [START_REF] Simon | Incorporating Molecular Evolution into Phylogenetic Analysis, and a New Compilation of Conserved Polymerase Chain Reaction Primers for Animal Mitochondrial DNA[END_REF]; [START_REF] Cameron | Insect Mitochondrial Genomics: Implications for Evolution and Phylogeny[END_REF]; [START_REF] Konstantinov | DNA import into mitochondria[END_REF]).

With the quick development of next generation sequencing (NGS) methods, the new sets of long-PCR primers defined here will allow to easily and cheaply acquire large amount of new data for large taxonomic samples. The new mitogenome described in the present study, combined with the existing datasets, arises as an invaluable resource for future comparative evolutionary genomic studies in crickets and insects, in general. TABLE S1. Genetic distance for complete mitogenome, all PCGs and two ribosomal genes showed by the percentage of pairwise identical bases (% identity). The average genetic distances and their standard deviation between two subfamilies (Eneopterinae and Gryllinae) and within the subfamily Gryllinae are also shown in the table. The abbreviations of species as following: C. muiri (C), L. equestris (L), T. emma (T) and V. hemelytrus (V).

Percentage of pairwise identical bases between subfamilies (%)

Percentage of pairwise identical bases within Gryllinae (%)

C / L C /T C / V Mean ± SD L / T L / V T / V Mean ±

FIGURE 2 .

 2 FIGURE 2. Tree trunk with males of Cardiodactylus muiri Otte, 2007 engaged in singing activity during late afternoon in Kumawa forest (Papua, Indonesia). Red circles mark males' positions; height ca. 4 m.

FIGURE 3 .

 3 FIGURE 3. Females of Cardiodactylus muiri Otte, 2007. A, walking on leaf litter during early night; B, ovipositing in a mossy tree root at night; in Kumawa (Papua, Indonesia).

FIGURE 4 .

 4 FIGURE 4. Calling song of Cardiodactylus muiri Otte, 2007. A, oscillogram of 19 syllables (= echemes); B-C, detailled oscillogram (B) and sonogram (C) of five syllables; D, detailed oscillogram of one syllable; E, frequency spectrum of one syllable.

FIGURE 5 .

 5 FIGURE 5. Map of the mitochondrial genome of Cardiodactylus muiri Otte, 2007.

FIGURE 7 .

 7 FIGURE 7. Comparison of AT-skews and GC-skews of Eneopterinae and Gryllinae. TABLE3. General nucleotide composition of complete mitogenome and nucleotide compositional skew of the Gryllidae mitogenomes, rrnL genes, rrnS genes and the A+T-rich regions in grylline crickets. A (%), T (%), G (%) and C (%) mean the percentage of adenine, thymine, guanine and cytosine. The number with asterisk (*) means the average nucleotide compositional skew in the subfamily Gryllinae.

FIGURE 8 .

 8 FIGURE 8. Percentage of pairwise identical bases (% identity) of mitogenome and each gene between the two subfamilies and within Gryllinae.

TABLE 1 .

 1 List of primers used to amplify the mitogenome and mapping reference in assembly procedure of Cardiodactylus muiriOtte, 2007.

	R1281	Reverse	AGTAATGATNTAGACTGCAATTCTAAAG	partial 12S

  Cardiodactylus muiriOtte, 2007aOtte, : 349-2007b: 30 : 30 (confirmation of depository)-Robillard et al. 2014: 68 (redescription).Synonym name. Cardiodactylus buru[START_REF] Robillard | Revision of the cricket genus Cardiodactylus (Orthoptera, Eneopterinae, Lebinthini): the species from both sides of the Wallace line, with description of 25 new species[END_REF], in Robillard et al. 2014: 25. New synonymy. The differences between C. buru and C. muiri in male genitalia described in[START_REF] Robillard | Revision of the cricket genus Cardiodactylus (Orthoptera, Eneopterinae, Lebinthini): the species from both sides of the Wallace line, with description of 25 new species[END_REF] are clearly present in specimens of the same populations in the newly examined material from New Guinea and Kai Islands. These differences must then be considered as intraspecific variation of C. muiri; C. buru is consequently a junior synonym of C. muiri.

	Cardiodactylus muiri Otte, 2007a
	(Figs 1-4)
	Type material. Male holotype: Indonesia: Malaka Province, Ceram, Piroe [Seram I., Piru], II.1909, F. Muir (BPBM) [not
	examined].
	Type locality. Indonesia, Seram Island, Piru.
	Distribution. Indonesia, Maluku islands (Ambon, Seram, Buru), Kai Islands, South-west part of New Guinea (surroundings of
	Kumawa and Kaimana).

Type species: Cardiodactylus novaeguineae

(Haan, 1842)

TABLE 2 .

 2 Annotation of the mitochondrial genome of Cardiodactylus muiri Otte, 2007.

	Gene	Strand	Anticodon	Start codon	Stop codon	Position
	tRNA-Ile	+	AUC			1-63 (0)
	tRNA-Gln	-	UUG			61-129 (0)
	tRNA-Met	+	AUG			143-211 (0)
	nd2	+		ATT	T	212-1229 (0)
	tRNA-Trp	+	UCA			1230-1294 (-8)
	tRNA-Cys	-	GCA			1287-1348 (0)
	tRNA-Tyr	-	GUA			1361-1425 (-8)
	cox1	+		ATC	TAA	1418-2965 (0)
	tRNA-Leu 2	+	UUA			2972-3036 (0)
	cox2	+		ATG	T	3042-3720 (+3)
	tRNA-Lys	+	AAG			3718-3788 (0)
	tRNA-Asp	+	GAC			3789-3854 (0)
	atp8	+		ATT	TAA	3855-4016 (+7)
	atp6	+		ATG	TAA	4010-4693 (+1)
	cox3	+		ATG	T	4693-5479 (0)
	tRNA-Gly	+	GGA			5480-5542 (0)
	nd3	+		ATA	TAA	5543-5896 (0)
	tRNA-Ala	+	GCA			5904-5968 (0)
	tRNA-Arg	+	CGA			5979-6041 (-1)
	tRNA-Glu	-	UUC			6041-6105 (0)
	tRNA-Ser 1	-	GCU			6110-6176 (0)
	tRNA-Asn	-	GUU			6177-6244 (0)
	tRNA-Phe	-	GAA			6246-6307 (0)
	nd5	-		ATT	T	6308-8021 (0)
	tRNA-His	-	GUG			8040-8103 (-1)
	nd4	-		ATG	T	8103-9445 (+7)
	nd4l	-		ATG	TAA	9439-9735 (0)
	tRNA-Thr	+	ACA			9747-9810 (0)
	tRNA-Pro	-	UGG			9811-9877 (0)
	nd6	+		ATC	TAA	9880-10404 (+1)
	cytb	+		ATG	TAA	10404-11537 (0)
	tRNA-Ser 2	+	UCA			11562-11624 (0)
	nd1	-		ATA	T	11625-12594 (+6)
	tRNA-Leu 1	-	UAG			12589-12656 (+41)
	rrnL	-				12616-13977 (+18)
	tRNA-Val	-	UAC			13960-14027 (0)
	rrnS	-				14028-14798 (0)
	A+T-rich region	+	N/A			14857-16328 (0)

TABLE S2 .

 S2 Comparison of Ka/Ks ratios for all PCGs. The average Ka/Ks ratio and their standard deviation between subfamilies (Eneopterinae and Gryllinae) within the subfamily Gryllinae are also shown in the table. The abbreviations of species as following: C.

									SD
	Mitogenome	71.80	72.01	71.80	71.87±0.12	74.82	73.76	78.62	75.73±2.56
	nd2	66.76	66.08	67.36	66.73±0.52	68.44	69.35	74.40	70.73±3.21
	cox1	83.1	83.63	83.27	83.34±0.22	84.41	82.65	87.15	84.74±2.27
	cox2	78.96	82.01	79.94	80.30±1.27	81.04	79.70	81.42	80.72±0.90
	atp8	59.88	64.15	57.86	60.63±2.62	64.15	64.15	69.93	66.08±3.34
	atp6	76.51	77.09	76.36	76.65±0.31	80.47	75.77	0.816	79.29±3.11
	cox3	79.77	80.15	80.03	79.98±0.16	80.03	81.17	83.21	81.47±1.61
	nd3	75.21	74.93	76.35	75.50±0.61	73.50	72.65	76.64	74.26±2.10
	nd5	73.73	75.60	76.42	75.25±1.13	77.51	77.21	81.09	78.60±2.16
	nd4	75.02	76.88	75.39	75.76±0.80	78.82	76.96	84.04	79.94±3.67
	nd4l	72.11	72.45	75.17	73.24±1.37	79.25	79.93	82.99	80.72±1.99
	nd6	68.01	66.29	67.82	67.37±0.77	69.19	64.72	69.19	67.70±2.58
	cytb	77.90	77.28	76.83	77.34±0.44	79.54	79.89	80.16	79.86±0.31
	nd1	76.19	73.55	78.70	76.15±2.10	74.22	78.60	0.798	77.55±2.95
	rrnL	71.73	73.13	69.90	71.59±1.32	75.78	72.06	78.71	75.52±3.33
	rrnS	69.50	71.66	69.67	70.28±0.98	76.28	74.04	81.04	77.12±3.57
		Ka/Ks ratio pairwise comparison			Ka/Ks ratio pairwise comparison within
		between subfamilies			Gryllinae
	Gene	C / L	C / T	C / V	Mean ± SD	L / T	L / V	T / V	Mean ± SD
	nd2	0.47	0.79	0.57	0.61±0.17	1.01	0.97	1.05	1.01±0.04
	cox1	1.00	0.91	0.95	0.95±0.05	1.01	1.25	0.73	1.00±0.26
	cox2	1.02	0.70	1.24	0.99±0.27	1.10	1.41	1.32	1.28±0.16
	atp8	0.07	0.15	0.34	0.19±0.13	0.77	0.32	0.75	0.61±0.25
	atp6	0.67	0.81	0.76	0.75±0.07	0.71	1.18	0.90	0.93±0.23
	cox3	0.72	0.79	0.74	0.75±0.04	1.11	0.88	0.91	0.97±0.13
	nd3	0.57	0.70	0.73	0.67±0.09	1.27	1.28	1.28	1.27±0.00
	nd5	0.56	0.41	0.33	0.43±0.12	0.66	0.65	0.61	0.64±0.03
	nd4	0.50	0.42	0.49	0.47±0.04	0.64	0.84	0.42	0.63±0.21
	nd4l	0.67	0.61	0.28	0.52±0.21	0.96	0.71	0.63	0.77±0.17
	nd6	0.68	0.72	0.61	0.67±0.05	0.76	1.35	0.81	0.97±0.32
	cytb	0.75	0.99	1.13	0.96±0.19	0.93	0.83	1.26	1.01±0.22
	nd1	0.73	0.81	0.56	0.70±0.13	0.66	0.66	0.58	0.63±0.05

muiri (C), L. equestris (L), T. emma (T) and V. hemelytrus (V).

Acknowledgements

Fieldwork was supported by the Project Lengguru 2014 (www.lengguru.org), conducted by the French Institut de Recherche pour le Développement (IRD), the Indonesian Institute of Sciences (LIPI), the University of Papua (UNIPA), the University of Cendrawasih (UNCEN), the University of Musamus (UNMUS) and the Sorong Fisheries Academy (APSOR) with corporate sponsorship from COLAS Group, Veolia Water and the Total Foundation. Laboratory access and assistance was provided by the "Service de Systematique Moleculaire" of the Muséum national d'Histoire naturelle, Paris (CNRS UMS 2700). This study was conducted in the context of the PhD thesis of JD, who was funded by China Scholarship Council (CSC), the Innovation Funds of Graduate Programs of Shaanxi Normal University (2012CXB019) and a supporting grant from "La Société des Amis du Muséum", and of NV, who was supported by the National Council for Scientific Development and Technological (CNPq-Brazil, PROTAX [158726/2013-8]), the Coordination for the Improvement of High Education Personnel (CAPES; PDSE-[BEX 2609/14-2]), and benefitted from financial support provided by CNPq and FAPEMIG, Process N. 563360/2010-0, Edital N. 47/2010 of the Ministério de Ciência e Tecnologia/CNPq/Ministério de Meio Ambiente/ Ministério da Educação/CAPES/Fundo Nacional de Desenvolvimento Científico e Tecnológico, Ação Transversal/Fundações de Apoio à Pesquisa, Sistema Nacional de Pesquisa em Biodiversidade, SISBIOTA Brazil. This work was also partially benefitted from the project of ICCM, who was supported by a grant from Agence Nationale de la Recherche under the LabEx ANR-10-LABX-0003-BCDiv, in the program "Investissements d'avenir" n ANR-11-IDEX-0004-02.