Charles Dapogny

Rafael Estevez

Alexis Faure

Georgios Michailidis

Shape and topology optimization

come   L'archive ouverte pluridisciplinaire

Introduction

Modern shape and topology optimization algorithms have by now widely penetrated the mechanical engineering community as valuable assets for industrial design [START_REF] Abballe | Rodin project, topology optimization 2.0?[END_REF][START_REF] Bendsoe | Topology optimization-broadening the areas of application[END_REF]. Not only do they help in improving the physical insight of engineers, but they are also well adapted to the design of structures satisfying the constraints imposed by the manufacturing process. To name a few, casting constraints [START_REF] Michailidis | Manufacturing Constraints and Multi-Phase Shape and Topology Optimization via a Level-Set Method[END_REF][START_REF] Xia | A level set based method for the optimization of cast part[END_REF][START_REF]Simultaneous optimization of cast part and parting direction using level set method[END_REF][START_REF] Zhou | Progress in topology optimization with manufacturing constraints[END_REF][START_REF] Zhou | Topology optimization with manufacturing constraints[END_REF], thermal constraints [START_REF] Ebrahimi | Sensitivity analysis and nonlinear programming applied to investment casting design[END_REF][START_REF] Michailidis | Manufacturing Constraints and Multi-Phase Shape and Topology Optimization via a Level-Set Method[END_REF][START_REF] Tavakoli | Optimal riser design in sand casting process by topology optimization with simp method i: Poisson approximation of nonlinear heat transfer equation[END_REF], and local size constraints (related e.g. to the thickness or the minimum distance between members) [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF][START_REF] Chen | Shape feature control in structural topology optimization[END_REF][START_REF] Guest | Achieving minimum length scale in topology optimization using nodal design variables and projection functions[END_REF][START_REF] Lazarov | Length scale and manufacturability in density-based topology optimization[END_REF] have been extensively discussed in the literature and are by now addressed by commercial softwares.

Parallel to this increasing enthusiasm for shape and topology optimization are the great headway in the development of additive manufacturing technologies, which herald revolutionary construction capabilities in the near future; see e.g. [START_REF] Barlier | Fabrication additive -Du Prototypage Rapide à l'impression 3D[END_REF][START_REF] Gibson | Additive manufacturing technologies: rapid prototyping to direct digital manufacturing[END_REF] and Chapter 1 in [START_REF] Clausen | Topology Optimization for Additive Manufacturing[END_REF]. The unprecedented complexity allowed by these techniques in terms of constructible designs when compared to traditional manufacturing methods makes it possible to process almost as is the intricated structures predicted by topology optimization. Beyond the exterior outline, additive technologies also offer the possibility to redesign (and optimize) the infill part of shapes, i.e. their interior, possibly porous regions, with the aim of reducing their total mass and improving their structural performance; see for instance [START_REF]Exploiting additive manufacturing infill in topology optimization for improved buckling load[END_REF].

Despite the opportunities they offer, additive manufacturing technologies suffer from limitations of a new nature, which have been the focus of multiple recent investigations. Not mentioning their potential lack of scalability and the significant amount of time needed by the construction process, the perhaps most critical issue is that of overhangs: additive manufacturing technologies generally experience difficulties in assembling large, nearly horizontal regions hanging over void without sufficient support from the lower structure. Several attempts have been made to understand these features and to prevent their emergence during the design optimization stage; see for instance [START_REF] Brackett | Topology optimization for additive manufacturing[END_REF][START_REF] Gaynor | Topology optimization considering overhang constraints: Eliminating sacrificial support material in additive manufacturing through design[END_REF][START_REF] Langelaar | An additive manufacturing filter for topology optimization of print-ready designs, Structural and Multidisciplinary Optimization[END_REF][START_REF]Topology optimization of 3d self-supporting structures for additive manufacturing[END_REF][START_REF] Mass | Topology optimization for additive manufacturing: Accounting for overhang limitations using a virtual skeleton[END_REF] in the framework of density-based topology optimization methods, and recently [START_REF] Allaire | Structural optimization under overhang constraints imposed by additive manufacturing technologies[END_REF][START_REF] Allaire | Shape optimization of a layer by layer mechanical constraint for additive manufacturing[END_REF][START_REF] Amir | Topology optimization for staged construction[END_REF], where a new formulation is proposed which appeals to the mechanical origin of overhangs.

Another major issue is to understand the influence of the peculiar, layer by layer additive construction process on the mechanical properties of the assembled structures, and notably the induced anisotropy of the constituent material [START_REF] Frazier | Metal additive manufacturing: a review[END_REF]. The phase transition phenomena governing the melting and cooling of the material during the assembly are very intricated; in particular, they depend on the selected additive technology and even possibly on the particular machine tool. When it comes to the design optimization stage of mechanical parts, the consideration of these anisotropic effects is key in the correct prediction of their structural response and in the accurate evaluation of their performance; see [START_REF] Guessasma | Challenges of additive manufacturing technologies from an optimisation perspective[END_REF][START_REF] Zhang | Role of anisotropic properties on topology optimization of additive manufactured load bearing structures[END_REF] about this point.

In this work, our contribution to shape and topology optimization in keeping with additive manufacturing processes is twofold. In a first time, we introduce a simplified, ad hoc model for the anisotropic features induced by several printing patterns for the two-dimensional layers featured by the additive construction process [START_REF] Liu | Concurrent deposition path planning and structural topology optimization for additive manufacturing[END_REF][START_REF] Liu | Deposition path planning-integrated structural topology optimization for 3d additive manufacturing subject to self-support constraint[END_REF]. In a nutshell, we assume different mechanical behaviors in the thin crust of material accounting for the structural boundary and in the interior infill region. The elastic material filling the crust is considered to be orthotropic, with principal axes depending on the local orientation of the boundary, while the properties of the constituent material of the infill region are driven by the pattern followed by the machine tool during the construction [START_REF] Ahn | Anisotropic material properties of fused deposition modeling abs[END_REF][START_REF] Lee | Measurement of anisotropic compressive strength of rapid prototyping parts[END_REF][START_REF] Ziemian | Anisotropic mechanical properties of ABS parts fabricated by fused deposition modelling[END_REF]. Moreover, inspired by the work presented in [START_REF] Liu | Concurrent deposition path planning and structural topology optimization for additive manufacturing[END_REF][START_REF] Liu | Deposition path planning-integrated structural topology optimization for 3d additive manufacturing subject to self-support constraint[END_REF], we extend this model to simulate the mechanical properties of parts printed along contour offsets.

In a second time, we consider the possibility of improving the structural performance of a shape by taking advantage of the non homogeneous material properties allowed by additive manufacturing technologies, namely by using a non-regular, optimized infill region: this region is seen as a whole new structure to be optimized by the same shape and topology optimization methods as those used for the exterior outline, under constraints on the available amount of material and on the maximum allowed structural thickness, along the lines of [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF][START_REF] Dapogny | Geometric constraints for shape and topology optimization in architectural design[END_REF][START_REF] Guest | Imposing maximum length scale in topology optimization[END_REF]. This echoes to the prior contributions [START_REF] Wu | Infill optimization for additive manufacturing-approaching bone-like porous structures[END_REF][START_REF] Wu | Minimum compliance topology optimization of shell-infill composites for additive manufacturing[END_REF], in which a local volume constraint is enforced so as to avoid large features and to obtain optimal non-regular infill structures. Let us mention another elegant and efficient way to get non-regular infill structures, introduced in [START_REF] Pantz | A post-treatment of the homogenization method for shape optimization[END_REF][START_REF]Construction of minimization sequences for shape optimization[END_REF], then used in [START_REF] Allaire | Topology optimization of modulated and oriented periodic microstructures by the homogenization method[END_REF][START_REF] Geoffroy-Donders | Optimization of oriented and parametric cellular structures by the homogenization method[END_REF][START_REF] Groen | Homogenization-based topology optimization for high-resolution manufacturable microstructures[END_REF], which is based on a post-treatment of optimized composite structures.

The remainder of this article is organized as follows. In Section 2, we provide a brief account of additive manufacturing technologies, and notably of the Fused Filament Fabrication methods, with a particular emphasis on why they induce anisotropy in the material properties of the constructed designs. In Section 3, we describe our theoretical and numerical framework for dealing with shape and topology optimization of elastic structures; next, in Section 4, we briefly recall elementary properties of the signed distance function to a domain, a key mathematical tool in our models. In Section 5, we present several anisotropic descriptions of the material properties of shapes produced by additive techniques; the question of how to optimize the (non-regular) infill region of structures is addressed in Section 6. In Section 7, the mathematical results needed in the device of our shape optimization algorithm are presented -namely the shape derivatives of the various objective functions at stake. Several numerical experiments are then conducted in Section 8. This article ends with three technical Appendices; Appendix A is a short reminder of the basic concepts of anisotropic linearized elasticity; Appendix B contains the change of basis formulae for elasticity tensors involved in the construction of our anisotropy models in Section 5. Finally, Appendix C outlines the main steps of the calculation of our shape derivatives -see Theorem 7.1 for the corresponding statement.

Let us already point out that the numerical experiments of this work are only two-dimensional, owing to unfortunate limitations imposed by our numerical environment; see Section 8.1 for more details. However, our models hold in greater generality; to emphasize this fact, the discussion takes place in the d-dimensional setting insofar as possible.

Additive manufacturing, Fused Filament Fabrication and anisotropy of the produced material

All additive manufacturing technologies assume the datum of the constructed shape under the form of a surface mesh (for instance in the popular STL format). This model is processed by a slicer, a software which converts it into a series of thin layers: a G-code is issued, containing the instructions for the machine tool to drive the successive construction of these layers, from the ground to the top of the structure (see Figure 1). Very different technologies are available to build these individual 2d layers, depending mainly on the nature of the processed material (metal, polymer, etc.); for example, some of them rely on a laser beam to melt a powder bed of metallic grains into the desired shape (Selective Laser Sintering, Electron Beam Melting), or to force a liquid photopolymer to harden (Stereolithography, Digital Light Processing). Another category of technologies, often gathered under the label of material extrusion methods, proceeds by direct deposition of a half-molten filament of material [START_REF] Turner | A review of melt extrusion additive manufacturing processes: I. process design and modeling[END_REF]. The so-called Fused Filament Fabrication methods (often abridged as FFF) fall into this class, which is of particular interest in the present article.

Fused Filament Fabrication methods were pioneered by the firm Stratasys under the name of Fused Deposition Modeling (FDM) as a means to process thermoplastic materials (such as ABS, PLA, PS,...).

FFF is a reliable and fairly cheap technology -an asset which makes it very popular for both domestic and industrial uses. While it was originally intended for rapid prototyping, it has recently aroused an increasing enthusiasm for the fabrication of end parts; for example, variations of the FFF strategy for assembling metallic parts have recently received much attention; see [START_REF] Mireles | Fused deposition modeling of metals[END_REF].

Roughly speaking, Fused Filament Fabrication techniques act by heating a filament of material up to its melting point, and by extruding it through a nozzle; the construction of each two-dimensional layer is then guided by the printing path followed by the machine tool. Generally, the boundary of the 2d layer is carefully printed in a first time, in a rather slow and accurate way; it forms a thin crust of material. Then, the bulk, or infill region is constructed; several strategies are available to this end, some of which being exemplified on Figure 2: • One possibility consists in printing the infill region in a fast way, according to a predefined pattern: often, rasters of material are deposited along a given orientation. We hereafter refer to this model as the crust-pattern model ; see Figure 3 (left). • Another possibility is to construct the whole 2d layer by following the contour offsets of the structural boundary. This is usually the way the upper and lower surfaces of three-dimensional structures are printed, so as to fully hide and protect the infill part. We shall refer to this model as the offset model ; see Figure 3 (right). • Last but not least, the infill design of the 2d layer (or even that of the whole three-dimensional structure under construction) may be subject to optimization, with the goal to enhance the structural performance: the infill design becomes an optimization variable on its own. We shall refer to this possibility as the crust-personalized infill model.

Constructed shape Printing contour

Anisotropic crust, isotropic bulk anisotropic bulk

Anisotropic crust, Offset model Isotropic crust, Infill bulk personalized infill (bottom, from the left to the right) assembly of the bulk of the layer using an isotropic infill, an anisotropic (horizontally oriented) infill, along contour offsets, and using a personalized infill.

The material properties of structures assembled by Fused Filament Fabrication methods are typically anisotropic; indeed, the bonds between regions that melt and solidify together are stronger than those between regions which solidify aside from one another; see for instance [START_REF] Ahn | Anisotropic material properties of fused deposition modeling abs[END_REF][START_REF] Lee | Measurement of anisotropic compressive strength of rapid prototyping parts[END_REF][START_REF] Ziemian | Anisotropic mechanical properties of ABS parts fabricated by fused deposition modelling[END_REF]. Therefore, the material properties of three-dimensional structures constructed by FFF are strongly influenced by the build direction, and by the particular model retained for the construction of each two-dimensional layer. For instance, in the aforementioned 'crust-pattern' model, the boundary of the structure has distinct properties from those of the bulk phase: it arises as a thin crust with a laminated structure (several strata are melted and cooled one atop the other), whose elastic behavior differs between the tangential and normal directions (see Figure 4).

The anisotropy in the properties of structures assembled by additive techniques has a strong impact on their mechanical performance; it is therefore of utmost importance to understand this behavior, and, if possible, to take advantage of it. In the following, we first propose a model for the anisotropic behavior predicted by the crust-pattern and offset models (Section 5), then we turn to the question of optimizing the infill design (Section 6).

Shape and topology optimization framework

Shape and topology optimization strategies mainly differ from one another by the considered category of shapes (e.g. 'black and white' designs, density functions,...), the evaluation of the sensitivity of the optimized criterion, and the numerical representation of designs and their evolution. In the present section, we briefly outline the main ingredients of our setting.

n ⌧ n E n E ⌧ Figure 4.
Orthotropic model for the elastic properties of the crust: the Young's modulus E is larger in the tangential direction τ to the structural boundary than in the normal direction n.

Shape optimization of linear elastic structures.

In our context, shapes are bounded Lipschitz domains Ω ⊂ R d . They are clamped on a subset Γ D of their boundary, and they are submitted to surface loads g ∈ L 2 (Γ N ) d , applied on a disjoint subset Γ N ⊂ ∂Ω. Body forces are omitted for simplicity. The boundary of any shape Ω is required to contain both regions Γ D and Γ N , so that the complementary free boundary Γ := ∂Ω \ (Γ D ∪ Γ D ) is the only region of ∂Ω which is subject to optimization.

In such circumstances, the elastic displacement u Ω of the structure Ω belongs to H 1 Γ D (Ω) d , where

H 1 Γ D (Ω) := u ∈ H 1 (Ω)
, u = 0 on Γ D , and it is the unique solution in this space to the linearized elasticity equations:

(3.1)        -div(σ Ω (u Ω )) = 0 in Ω, u Ω = 0 on Γ D , σ Ω (u Ω )n = g on Γ N , σ Ω (u Ω )n = 0 on Γ.
In the above system, σ Ω (u) is the stress tensor associated to a displacement u. The relation between σ Ω (u) and the linearized strain tensor e(u) := 1 2 (∇u + ∇u T ) is described by the compliance tensor S Ω ∈ L(S d (R), S d (R)), where L(S d (R), S d (R)) contains the linear mappings from the set S d (R) of real, d × d symmetric matrices into itself. An equivalent connection involves the Hooke's tensor A Ω , which is the inverse of S Ω : e(u) = S Ω σ Ω (u) and σ Ω (u) = A Ω e(u); see Appendix A for additional details about the physical information encoded in S Ω and A Ω . The Ω subscript in the above notations reflects the dependence of the material properties with respect to the domain Ω itself, which it is our aim to model in the next sections.

Our purpose is to optimize the external design of shapes Ω with respect to some performance criterion; to be quite precise, we solve the problem:

(3.2) min Ω∈U ad J(Ω),
where J(Ω) is a given function of the domain. In the present work, we typically consider J(Ω) = C(Ω), the structural compliance of a shape Ω:

(3.3) C(Ω) = Ω A Ω e(u Ω ) : e(u Ω ) dx = Γ N g • u Ω ds,
which measures the energy stored in Ω in its deformed configuration, or equivalently the work of external loads acting on Ω. Other objective functions could be considered -such as least-square functionals over the elastic displacement u Ω , stress-based functionals, etc. Admittedly, compliance is a particularly simple performance criterion to be optimized in shape and topology optimization, and dealing with more intricated objective functions would cause concerns related to e.g. the existence of lots of local minima; nevertheless, the analyses carried out in the present article would be essentially similar.

The optimization problem (3.2) takes place over a set U ad of admissible shapes Ω, which are smooth and enclose the fixed regions Γ D and Γ N as parts of their boundaries.

In practice, (3.2) is often supplemented with constraints which are modelled by shape functionals. Two simple examples of such are the volume Vol(Ω) = Ω dx and the perimeter Per(Ω) = ∂Ω ds. In our context where shapes are filled with material in an inhomogeneous way (in particular, the material density is not uniform within the structure), a more relevant constraint is related to the mass M (Ω) of the structure. The mathematical definition of M (Ω) depends on the particular model retained for the distribution of the material inside Ω, and this issue is discussed in Section 5.4.

Shape sensitivity using Hadamard's method.

Most optimization algorithms for problems of the form (3.2) rely on the knowledge of the derivatives of the objective and constraint functionals. Since the optimization variable is the shape itself in our applications, this requires a notion of differentiation with respect to the domain. We rely on the classical Hadamard's boundary variation method [START_REF]Conception optimale de structures[END_REF][START_REF] Henrot | Variation et optimisation de formes: une analyse géométrique[END_REF][START_REF] Simon | Sur le contrôle par un domaine géométrique[END_REF][START_REF] Soko | Introduction to shape optimization[END_REF], whereby variations of a shape Ω are considered under the form:

Ω θ := (Id + θ)(Ω), θ ∈ W 1,∞ (R d , R d ), ||θ|| W 1,∞ (R d ,R d ) < 1.
Definition 3.1. A function J(Ω) of the domain is shape differentiable at Ω if the underlying mapping θ → J(Ω θ ), from W 1,∞ (R d , R d ) into R, is differentiable at θ = 0 in the sense of Fréchet; denoting by J (Ω)(θ) the corresponding derivative, the following expansion holds in the neighborhood of θ = 0:

J(Ω θ ) = J(Ω) + J (Ω)(θ) + o(θ) where lim θ→0 |o(θ)| θ W 1,∞ (R d ,R d ) = 0.
According to the Structure theorem (see e.g. [START_REF] Henrot | Variation et optimisation de formes: une analyse géométrique[END_REF], §5.9), the shape derivative J (Ω)(θ) depends only on the normal trace θ • n of the deformation θ on the boundary ∂Ω. Actually, for a wide class of shape functionals, J (Ω)(θ) has the more precise structure:

(3.4) J (Ω)(θ) = ∂Ω w Ω θ • nds,
where the integrand w Ω depends on the function J(Ω), on the state u Ω , and possibly on an adjoint state.

It is easy to infer a descent direction for J(Ω) from (3.4); indeed, letting

(3.5) θ = -w Ω n,
one has, for a small enough descent step t > 0:

J(Ω tθ ) = J (Ω) -t ∂Ω w 2 Ω ds + O(t 2 ) < J(Ω),
so that Ω tθ is 'better' than Ω in terms of J(Ω). More elaborated strategies, relying on a change of inner products to identify a shape gradient from (3.4) are described in e.g. [START_REF] Burger | A framework for the construction of level set methods for shape optimization and reconstruction[END_REF][START_REF] Gournay | Velocity extension for the level-set method and multiple eigenvalues in shape optimization[END_REF]. See also Remark C.1 below about this point.

Remark 3.1. So that variations Ω θ of admissible shapes Ω ∈ U ad stay admissible, the considered deformations θ are in practice restrained to a set Θ ad of admissible deformations; in the setting of Section 3.1, these are smooth and they vanish on the non optimizable region Γ D ∪ Γ N .

3.3. Evolution of shapes using the Level Set method.

The Level Set method was pioneered by S. Osher and J. Sethian in [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[END_REF], then introduced in the shape optimization context in [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF][START_REF] Osher | Level set methods for optimization problems involving geometry and constraints: I. frequencies of a two-density inhomogeneous drum[END_REF][START_REF] Sethian | Structural boundary design via level set and immersed interface methods[END_REF][START_REF] Wang | A level set method for structural topology optimization[END_REF]. It relies on an implicit representation of an arbitrary shape Ω ⊂ R d , as the negative subdomain of an auxiliary 'level set' function φ : D → R, defined on a larger, fixed computational domain D ⊂ R d ; more precisely,

(3.6)    φ(x) < 0 if x ∈ Ω, φ(x) = 0 if x ∈ ∂Ω, φ(x) > 0 if x ∈ D \ Ω.
If Ω ≡ Ω(t) is evolving in time according to a velocity field with (scalar) normal component V (t, x), the motion translates in terms of an associated level set function φ(t, •) -i.e. (3.6) holds at each time t -as the following Hamilton-Jacobi equation:

(3.7) ∂φ ∂t (t, x) + V (t, x)|∇φ(t, x)| = 0, t > 0, x ∈ D.
In our applications, the velocity V (t, x) is the shape gradient of the optimized functional J(Ω) (the integrand w Ω in (3.4), (3.5)) and the (pseudo-)time t > 0 is interpreted as a descent step for the optimization algorithm. As far as the numerical approximation of (3.7) is concerned, one may use an explicit second-order upwind scheme when D is equipped with a Cartesian grid [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF][START_REF] Sethian | Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science[END_REF], or a semi-Lagrangian algorithm based on the method of characteristics when the computational support is a simplicial (triangular in 2d) mesh [START_REF] Bui | An accurate anisotropic adaptation method for solving the level set advection equation[END_REF][START_REF] Strain | Semi-lagrangian methods for level set equations[END_REF].

Eventually, let us underline a classical difficulty arising from the use of the Level Set method in the shape optimization context. The velocity field V (t, x) guiding the motion of the shape is calculated from the solution u Ω to (3.1), a partial differential equation posed on Ω. However, no mesh of Ω is available to calculate u Ω with e.g. the Finite Element method; to circumvent this drawback, we rely on the so-called 'ersatz' material approach [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF], whereby the void D \ Ω is filled with a very soft elastic material (i.e. with a small, yet positive Young's modulus); this allows to approximate (3.1) by a system posed on D; see also Section 5.1.3.

'Classical' facts around the signed distance function

As we have mentioned in Section 2, the local material properties of structures assembled via additive processes strongly depend on the distance to the structural boundary and on the orientation of the latter. For this reason, the signed distance function to a domain plays a key role in our models and we recall the basic notions, referring to [START_REF] Allaire | Mutli-phase structural optimization via a level-set method[END_REF][START_REF] Delfour | Shapes and geometries: metrics, analysis, differential calculus, and optimization[END_REF] for details.

Let us start with a definition: When it comes to the study of the signed distance function, the following notions are particularly relevant; see Figure 5 (left) for an illustration. Definition 4.2. Let Ω ⊂ R d be a bounded Lipschitz domain.

(4.1) ∀x ∈ R d , d Ω (x) =    -d(x, ∂Ω) if x ∈ Ω, 0 if x ∈ ∂Ω, d(x, ∂Ω) if x ∈ c Ω,
• Let x ∈ R d . A point p ∈ ∂Ω is a projection of x onto ∂Ω if d(x, ∂Ω) = |x -p|, i.e. p is a minimizer in the definition (4.2). • The set Π ∂Ω (x) of all such points is the set of projections of x onto ∂Ω. When Π ∂Ω (x) is a singleton, its unique element is denoted p ∂Ω (x) and is called the projection of x onto ∂Ω.

• Let p ∈ ∂Ω; the ray emerging from p is the set ray ∂Ω (p) defined by:

ray ∂Ω (p) = x ∈ R d , Π ∂Ω (x)
is a singleton and p ∂Ω (x) = p .

• The skeleton Σ of Ω is the set of points in R d where the squared distance function x → d(x, ∂Ω) 2 is not differentiable.

⌦ ⌃ • • x p @⌦ (x) x 0 • • • y 1 y 2 z • ray @⌦ (z) ⌦ ⌃ x • • p @⌦ (x) rd ⌦ (x)
n(p @⌦ (x)) We now recall a classical result about the differentiability of the signed distance function, in a deliberately loose manner as far as smoothness assumptions are concerned; see [START_REF] Delfour | Shapes and geometries: metrics, analysis, differential calculus, and optimization[END_REF][START_REF] Henrot | Variation et optimisation de formes: une analyse géométrique[END_REF] and Figure 5 (right).

Theorem 4.1. Let Ω ⊂ R d be a smooth bounded domain; there exists a tubular neighborhood U of ∂Ω such that d Ω is smooth on U. Moreover, for any point x ∈ R d ,

• Either x ∈ ∂Ω; then d Ω is differentiable at x, with ∇d Ω (x) = n(x), the unit normal vector to ∂Ω at x, pointing outward Ω; • Or x ∈ R d \ ∂Ω; then d Ω is differentiable at x if and only if x belongs to the complement of the skeleton Σ of Ω. In such a case, there exists a unique projection p ∂Ω (x) of x onto ∂Ω and the gradient of d Ω at x reads

(4.3) ∇d Ω (x) = x -p ∂Ω (x) d Ω (x) .
In particular ∇d Ω takes constant values on each set ray ∂Ω (z), z ∈ ∂Ω.

As a consequence of Theorem 4.1, the normal vector field n : ∂Ω → R d and any tangential vector field τ : ∂Ω → R d may be extended to U, and actually to R d \ Σ via the formulae:

(4.4) n(x) := n(p ∂Ω (x)), and τ (x) := τ (p ∂Ω (x)), x ∈ R d \ Σ.
This convention is used throughout the article.

The signed distance function is key in the description of many geometric features of shapes, such as the notion of maximum thickness, as introduced in [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF].

(4.5) ∀x ∈ Ω, d Ω (x) ≥ -δ/2.
Before ending this section, let us mention that it is possible to calculate shape derivatives of functionals involving the signed distance function d Ω -for instance the derivative of the mapping Ω → d Ω (x), for a given point x ∈ R d . This comes in handy in the next sections, since the elasticity tensors associated to the models for the material properties of shapes discussed in Section 2 bring d Ω into play. We omit the technical details for brevity and refer to [START_REF] Allaire | Mutli-phase structural optimization via a level-set method[END_REF][START_REF] Delfour | Shape identification via metrics constructed from the oriented distance function[END_REF] about this issue.

Anisotropic modeling of the mechanical properties of shapes constructed by additive processes

We now discuss the anisotropic material properties of a shape Ω ⊂ D induced by the 'crust-pattern' and the 'offset' models introduced in Section 2. Sections 5.1, 5.2, 5.3 and 5.4 below deal with two-dimensional shapes. From the physical point of view, this amounts to focusing on the assembly of one individual 2d layer of material. The extension of these concepts to the three-dimensional case is considered in Section 5.5.

In the remainder of this article, we denote by

A ref (resp. E ref , ν ref , B ref
) the Hooke's tensor (resp. the Young's modulus, the Poisson's ratio and the bulk modulus) of the isotropic constituent material of shapes when they are constructed by the molding process, the reference situation in our study. Consistently, we use a crust (resp. bulk ) subscript for all the elastic properties within the crust (resp. bulk) region of shapes, and we indicate with a cp (resp. off ) superscript that these are attached to the 'crust-pattern' (resp. 'offset') model.

Material properties of shapes using the 'crust-pattern' model.

In the 'crust-pattern' model, the computational domain D is decomposed into three regions, as exemplified on Figure 6, (left):

• The crust of Ω is the band {x ∈ D, -ε < d Ω (x) < ε} with width 2ε around ∂Ω -note that depending on the particular application, ε may not be 'small' (although it is in ours);

• The interior of Ω, or 'infill' region is the set {x ∈ D, d Ω (x) < -ε}; • The void region is {x ∈ D, d Ω (x) > ε}.
In keeping with this decomposition, the entries of the elasticity tensor A cp Ω (x) at a particular point x ∈ D depend on the distance from x to the boundary ∂Ω, and on the local orientation of ∂Ω, say:

(5.1)

A cp Ω (x) = A cp (d Ω (x), n(x)
), where we recall the convention (4.4). The tensor A cp (d, n) is constructed by gluing together three Hooke's tensors A cp bulk , A crust (n) and A void accounting for the local material properties in the bulk, crust and void regions respectively, and whose expressions are discussed in the next sections:

(5.2)

A cp (d, n) = h m (d) A cp bulk + (1 -h m (d) -h p (d))A crust (n) + h p (d)
A void , where the smooth interpolation functions h m , h p : R → R are defined by:

(5.3) h m (s) =    1 if s < -ε, -s ε + 1 2π sin( 2πs ε ) 1 β if -ε ≤ s ≤ 0, 0 if s > 0, and
(5.4) h p (s) =    0 if s < 0, s ε -1 2π sin( 2πs ε ) n ⌦ " " " 0 n " P bulk P crust P void

Material properties

A crust (n) in the crust.
As we have mentionned in Section 2, it is somewhat natural to consider the elastic properties of the crust to be orthotropic, with stronger stiffness in the tangential direction to the structural boundary ∂Ω. More precisely, let τ stand for the 90 degree clockwise rotate of n; the material properties at a point x in the crust are given by the following compliance tensor, expressed using Voigt's notation, in the local orthonormal frame (τ (x), n(x)) of the plane:

(5.5)

S crust =    1 E c τ - ν c nτ E c n 0 - ν c τ n E c τ 1 E c n 0 0 0 1 µ c    ,
where E c τ , E c n are the Young's moduli in the tangential, and normal directions respectively, ν c τ n and ν c nτ are the Poisson's ratios, and µ c is the shear modulus; see Appendix A for details. Details about the choice of these parameters are presented in Section 5.3; for now, let us simply mention that typically, E c τ is expected to take larger values than E c n . The expression S crust (n) of the compliance tensor S crust in the canonical basis (e 1 , e 2 ) of the plane then reads:

S crust (n) = R e (n) -1 S crust R σ (n),
where the rotation matrices R e (n) and R σ (n) are defined by (B.3) and (B.4) respectively. Eventually, the Hooke's tensor A crust (n) in the crust region is defined by:

(5.6)

A crust (n) = (S crust (n)) -1 .

5.1.2.

Material properties A cp bulk in the bulk of the structure.

The anisotropy in the elastic properties of the infill of shapes is dictated by the pattern adopted by the machine tool to assemble this region. In our case, the material is deposited as rasters of material with orientation τ 0 = (τ 0 1 , τ 0 2 ). Letting n 0 := (-τ 0 2 , τ 0 1 ), we write, as in Section 5.1.1:

S cp bulk = R e (n 0 ) -1 S cp bulk R σ (n 0 ), where (5.7) S cp bulk =     1 E b τ - ν b nτ E b n 0 - ν b τ n E b τ 1 E b n 0 0 0 1 µ b     .
Here, E b τ and E b n are the Young's moduli in the respective directions τ 0 and n 0 , the Poisson ratios ν b nτ and ν b τ n are relative to these directions, and µ b is the shear modulus. Again, see Section 5.3 about the practical choice of these parameters. We finally define:

(5.8)

A cp bulk = (S cp bulk ) -1 . 5.1.3. Material properties A void in the void.

In the void region D \ Ω, according to the 'ersatz' material approximation (see Section 3.3), the elasticity tensor reads:

(5.9)

A void = ηA ref ,
where η 1 is a given parameter.

Remark 5.1.

• The whole Hooke's tensor is interpolated in the definition (5.1), (5.2) of A cp Ω , and not the different material parameters independently (the Young's moduli, Poisson's ratios). The latter option indeed proves difficult when the combined materials have different principal directions.

• Although we have hitherto focused on anisotropic elastic properties in the crust and bulk of shapes, the above models obviously make it possible to describe isotropic properties in these regions, a situation that we also consider in the numerical examples of Section 8.

Material properties of structures using the 'offset' model.

A similar mathematical representation to that of Section 5.1.1 can be applied to construct the Hooke's tensor A off Ω accounting for the material properties of a 2d layer printed along contour offsets. The procedure is actually easier in this case since we only need to glue together the tensors A off bulk (n) and A void for the bulk and void parts respectively; more precisely:

(5.10)

A off Ω (x) = A off (d Ω (x), n(x)), where A off (d, n) = h o (d) A off bulk (n) + (1 -h o (d))
A void , the interpolation profile h o : R → R is defined by: (5.11)

h o (s) =    1 if s < -ε, 1 2 -s 2ε -1 2π sin πs ε if -ε ≤ s ≤ ε, 0 if s > ε,
where in this case, the parameter ε > 0 is 'small', and represents the half-thickness of the transition region between the shape Ω and the void region. The tensor A off bulk (n) reads exactly as that A crust (n) given in Section 5.1.1 (see (5.5) Remark 5.2. Our formulation for the anisotropic mechanical properties of shapes constructed according to the 'offset' model is similar to that proposed in [START_REF] Liu | Concurrent deposition path planning and structural topology optimization for additive manufacturing[END_REF][START_REF] Liu | Deposition path planning-integrated structural topology optimization for 3d additive manufacturing subject to self-support constraint[END_REF]. In those works however, the Level Set method is used for the numerical representation of shapes, but also in the definition of their elastic properties when they are printed along contour offsets: using the notations of the present article, the definition of A off bulk (n) featured in [START_REF] Liu | Concurrent deposition path planning and structural topology optimization for additive manufacturing[END_REF][START_REF] Liu | Deposition path planning-integrated structural topology optimization for 3d additive manufacturing subject to self-support constraint[END_REF] is based on the level set function φ and its gradient ∇φ instead of d Ω and n Ω respectively. The sensitivity of this tensor with respect to the domain is evaluated by taking derivatives with respect to φ. On the contrary, our model relies solely on intrinseque geometric quantities of shapes, (the normal vector, the signed distance function), which are independent of the particular numerical representation of shapes. From the theoretical point of view, this allows for a rigorous mathematical analysis of the shape differentiability of the considered objective functionals (see Section 7), while from the practical viewpoint, it prevents the undesired trend that the numerical algorithm may take advantage of the optimization of the values of φ far from its 0 level set, which are physically irrelevant; see [START_REF] Allaire | Mutli-phase structural optimization via a level-set method[END_REF] for a detailed discussion of these aspects.

Material parameters of the effective elasticity tensors in the crust and bulk regions.

In this section, we discuss the values of the effective Young's moduli, Poisson's ratios, and shear moduli involved in the expressions of A crust (n), A cp bulk and A off bulk (n) in the basis of their principal axes -see (5.5) and (5.7). These values depend on the density ρ crust , ρ cp bulk or ρ off bulk of the crust and bulk regions of shapes, which are in turn related to the size of the air gap between rasters of deposited material.

An intuitive idea, inspired from the geometric arrangement of the deposited material, is to describe the crust or the bulk region as a rank 1 laminated structure, i.e. a mixture of the reference material A ref and void, arranged as bands at the microscopic scale (see e.g. [START_REF] Allaire | Shape optimization by the homogenization method[END_REF] or [START_REF] Milton | The theory of composites, The Theory of Composites[END_REF], Chap. 22). Unfortunately, applying the homogenization theory to this crude model results in effective elasticity tensors A crust (n), A cp bulk and A off bulk (n) with non physical behaviors; their shear components are all equal to 0 (see e.g. [START_REF] Allaire | Shape optimization by the homogenization method[END_REF][START_REF]Conception optimale de structures[END_REF] about this point).

In reality, the material filling each 2d slice of a 3d structure retains a certain shear resistance from the upper and lower layers -an effect which is difficult to predict with a purely two-dimensional model. For this reason, we use in practice an ad hoc description for the elastic properties in the 'crust-pattern' and 'offset' models, based on the following considerations:

• The crust region in the 'crust-pattern' model is assumed to be fully dense, i.e. ρ crust = 1, which reflects the fact that it is printed with particular care; see Section 2. The entries of the tensor A crust (n) in the basis of its principal axes (i.e. the entries of S crust as defined in (5.5)) are based on experiments; see Section 8 for numerical values. • The material properties of the tensors A cp bulk and A off bulk (n), which may not be fully dense, are inferred from their fully-dense counterparts, relying on the homogenization theory (again, see [START_REF] Allaire | Shape optimization by the homogenization method[END_REF] or [START_REF] Milton | The theory of composites, The Theory of Composites[END_REF]):

-When the material filling the bulk region is assumed to be isotropic with density ρ = ρ cp bulk or ρ off bulk , we regard it as a composite structure formed by mixing the reference material A ref and void in proportions ρ and (1ρ) respectively. More precisely, we describe this material as one attaining the upper Hashin-Shtrikman bounds [START_REF] Hashin | A variational approach to the theory of the elastic behaviour of multiphase materials[END_REF] for its shear modulus µ bulk = µ cp bulk or µ off bulk and bulk modulus B bulk = B cp bulk or B off bulk . In other terms, these two quantities are the largest possible among all isotropic materials obtained by mixing A ref and void in proportions ρ and (1ρ):

µ bulk = µ ref - 2µ ref (µ ref + B ref )(1 -ρ) 2µ ref (1 -ρ) + (2 -ρ)B ref and B bulk = B ref - B ref (B ref + µ ref )(1 -ρ) (1 -ρ)B ref + µ ref ; see [3], Theorem 2.3.2.
3 for a proof of existence of materials attaining these bounds. -When the material filling the bulk region is supposed to be orthotropic, we regard it as a composite structure formed by mixing the fully dense material described by A crust (n) and void in proportions ρ and (1ρ). More precisely, we assume that this material attains the upper Voigt-Reuss bound, i.e. that it is the stiffest possible among all composite mixtures of A crust (n) and void in proportions ρ and (1ρ):

A cp bulk = ρA crust (n 0 ) and A off bulk (n) = ρA crust (n),
where n 0 is defined in Section 5.1.2. As far as attainability of these bounds is concerned, the situation is far more complicated than in the isotropic case, and we refer to [START_REF] Milton | On the possible effective elasticity tensors of 2-dimensional and 3-dimensional printed materials[END_REF], Section 2, for a recent overview.

Remark 5.3. In our applications, the density ρ = ρ crust , ρ cp bulk or ρ off bulk of the crust or bulk phases, and the orientation n 0 of the rasters in the 'crust-pattern' model are user-defined parameters, accounting respectively for the air gap between rasters of material, and the orientation of the machine tool; they are fixed throughout the optimization process of the shape Ω. However, they could be optimized as well, relying on similar techniques to those presented in this article (in a simpler framework of parametric optimization, however).

Mass of the structure.

In the 'crust-pattern' model, the mass of the structure Ω depends on the densities ρ crust and ρ cp bulk of the crust and bulk regions respectively. Assuming that the thickness 2ε of the crust is 'small', the volume of this region is approximately 2εPer(Ω), and that of the bulk is (Vol(Ω) -εPer(Ω)), so that the mass of the structure Ω reads in this case:

(5.12)

M cp (Ω) = 2ερ crust Per(Ω) + ρ cp bulk (Vol(Ω) -εPer(Ω)), where we recall from Section 5.3 that ρ crust = 1.

In the 'offset' model, the only parameter at play is the constant density of material ρ off bulk , so that the mass M off (Ω) of the structure is simply:

M off (Ω) = ρ off bulk Vol(Ω).
Remark 5.4. When the 'crust-pattern' model is used and the mechanical properties of the crust are better than those of the bulk with respect to the optimized structural criterion, the numerical shape optimization problem becomes ill-posed. Severe oscillations appear on the boundary of the shape as the algorithm creates artificial interface regions to take advantage of their better properties. Constraining the volume of the crust, and thereby the perimeter of the shape, is a classical remedy to this artifact. Another solution consists in including a minimum length scale to the optimization problem, either by imposing feature-size constraints [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF] or by using projection methods [START_REF] Faure | Design of thermoelastic multi-material structures with graded interfaces using topology optimization[END_REF].

Extension to three space dimensions.

The models presented in Section 5.1 and 5.2 account for the material properties of 2d structures, and they may in particular be applied to each horizontal layer involved in the fabrication of a 3d structure [START_REF] Liu | Deposition path planning-integrated structural topology optimization for 3d additive manufacturing subject to self-support constraint[END_REF]. They are also the building blocks for a similar model in the general case of whole 3d shapes, as we now outline. The orthotropy of the material properties of a three-dimensional structure Ω built by an additive technology is dictated by the build direction, and by the path followed by the machine tool for the assembly of the horizontal 2d layers. Denoting by (e 1 , e 2 , e 3 ) the canonical basis of R 3 (hoping that this does not induce confusion with our notation for strain tensors), we assume this build direction to be the third coordinate vector e 3 . Let h min (resp. h max ) be the minimum (resp. maximum) height of a point in Ω, and for h ∈ (h min , h max ), let Ω h := {x = (x 1 , x 2 , x 3 ) ∈ Ω, x 3 = h} be the horizontal layer at height h appearing in the construction of Ω; we also denote D h := {x ∈ D, x 3 = h}. Elaborating on Sections 5.1 and 5.2, the material properties inside a horizontal slice Ω h are described in terms of the signed distance function d Ω h : D h → R to Ω h in this layer and of the two horizontal vectors τ H and n H defined by:

n H = n -(n • e 3 )e 3 |n -(n • e 3 )e 3 |
, and τ H = n H × e 3 ;

see Figure 8 for an illustration. To be precise:

• If each 2d layer Ω h is constructed according to the 'crust-pattern' model, the bulk, crust and void regions of the 3d structure Ω are respectively defined by:

x ∈ D, d Ωx 3 (x) < -ε , x ∈ D, -ε < d Ωx 3 (x) < ε , and x ∈ D, d Ωx 3 (x) > ε .
The elastic tensor A cp Ω is then constructed by gluing together contributions from these regions (i.e. via a formula of the form (5.1), (5.2) where d Ω and n are replaced by d Ω h and n H respectively) in such a way that:

-In the crust, A cp Ω is orthotropic with principal axes τ H , n H and e 3 . The Young's modulus is weaker along the build and normal directions e 3 and n H than along τ H . -In the bulk, A cp Ω is orthotropic with principal axes τ 0 , n 0 and e 3 , where τ 0 is the (horizontal) deposition direction of the rasters of materials and n 0 = e 3 × τ 0 . Again, the Young's modulus is weaker along e 3 and n 0 than along τ 0 . -In the void, the properties of A cp Ω are still given by (5.9). • Likewise, if the 'offset' model of Section 5.2 is retained for the construction of each 2d layer, the orthotropic Hooke's tensor A off Ω (x) at a given point x ∈ D has for principal axes the build direction e 3 and the horizontal tangent and normal vectors τ H (x) and n H (x). These 3d expressions of the Hooke's tensors A cp Ω and A off Ω are admittedly more complicated than the two-dimensional counterparts under scrutiny in this article (notably because of the structure of the rotation matrices transforming the canonical frame of R 3 into (τ H , n H , e 3 )), but let us emphasize that there is no conceptual obstruction to the treatment of this three-dimensional extension.

Modeling and optimization of the infill design

We have hitherto focused on the description of the elastic properties of shapes when the infill design follows a (uniform) predefined pattern. In practice, patterns are mainly used for the bulk region when the structural efficiency of the assembled part is not the main target, but one is rather interested in the aesthetics of the shape or in building a prototype.

On the contrary, when the structural efficiency is the priority, significant improvements may be achieved on the total mass or on the mechanical performance of the shape by optimizing not only its exterior outlook, but also the infill design.

In this direction, one possibility, investigated in [START_REF] Geoffroy-Donders | Optimization of oriented and parametric cellular structures by the homogenization method[END_REF][START_REF] Groen | Homogenization-based topology optimization for high-resolution manufacturable microstructures[END_REF], is to consider a regular, lattice-shaped infill pattern, parametrized by its local thickness and by the local orientation of bars. These parameters are then the optimization variables.

Alternatively, one could regard the boundary of the constructed shape Ω as a new working domain, inside which the non regular infill design is sought as a general subdomain ω of Ω (i.e. which is not necessarily parametrized as a set of bars) to be optimized with respect to a criterion related e.g. to the mass of the structure. In fact, using such a non regular bulk region offers several advantages: beyond alleviating the burden of printing an infill structure with variable density and orientation, it enhances the robustness with respect to manufacturing errors that are likely to occur in additive manufacturing, especially at regions with low density [START_REF] Wu | Infill optimization for additive manufacturing-approaching bone-like porous structures[END_REF].

Such a joint optimization of the exterior boundary ∂Ω of a shape Ω and of its non regular infill ω ⊂ Ω may follow at least three different trails:

(1) One idea is to search for an optimal infill design ω inside the fixed computational domain D, under a constraint on the minimum length-scale allowed by the machine-tool, or under a maximum lengthscale constraint so as to drive ω towards a lightweight, porous structure. In a post-processing stage, the structural shape Ω is revealed as the exterior envelope of the region where this infill has high density. Our experience suggests that doing so results in local minima of poor structural performance: in a nutshell, the optimization of ω generally leads to a design filling almost all the computational domain D (with a non uniform density though), so that the exterior structure Ω resulting from the post-processing operation differs little from D as a whole; in other terms, it is difficult to interprete such an infill design ω ⊂ D in terms of a structure composed of a crust and a personalized infill; see the numerical illustrations in Section 8.3 about this point. (2) A second idea consists in optimizing simultaneously the exterior outline ∂Ω of the shape Ω and its infill design ω ⊂ Ω. In this spirit, in the work [START_REF] Wu | Minimum compliance topology optimization of shell-infill composites for additive manufacturing[END_REF] -arising in the context of density-based method -, the authors consider two distinct density functions, one for the exterior shape Ω, and the other representing the infill design ω. The material properties inside the computational domain D are modelled by relying on adapted interpolation schemes depending on Ω and ω, and the considered performance criterion is then minimized with respect to both designs at the same time. However relevant this approach, it seems a little difficult to implement in our context. (3) One last possibility -that retained in the present article -is to carry out a two-step optimization process: the infill design ω is only optimized after the exterior design of the shape Ω is set.

Let us describe this last alternative with a little more details: the original shape optimization problem (3.2) is solved for the exterior structure, assuming for instance that the constituent material is homogeneous (isotropic or anisotropic). Then, the resulting shape Ω becomes a new design space, and the infill region is a structure ω ⊂ Ω to be optimized in itself, using the same shape and topology optimization techniques as for Ω. In other words, we solve the new problem: min ω⊂Ω J(ω), imposing an additional constraint on ω so as to drive it toward an elaborate structure mimicking a network of bars, rather than 'bulky', thick regions. Namely, we impose a constraint on the maximum thickness in the spirit of [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF][START_REF] Dapogny | Geometric constraints for shape and topology optimization in architectural design[END_REF][START_REF] Michailidis | Manufacturing Constraints and Multi-Phase Shape and Topology Optimization via a Level-Set Method[END_REF]:

∀x ∈ ω, d ω (x) ≥ -δ/2,
where δ > 0 is the desired, user-defined maximum thickness; see Definition 4.3. This pointwise constraint is aggregated into a global one as: (6.1)

P MT (ω) ≤ δ 2 , where P MT (ω) :=     1 Ω h(d ω )dx ω h(d ω ) (d ω ) 2 dx     1 2
, and h : R → R is a regularized version of the characteristic function of the set {s ∈ R, |s|> δ/2}:

h(s) = 1 2 1 + tanh |s| -δ/2 η δ/2 ,
η > 0 being a parameter that controls the regularization of the constraint. In this work, we have chosen η = 2 5 ∆x, where ∆x is the maximum size of an element in the computational mesh.

Remark 6.1. A completely different and quite interesting method for creating non regular infills is the so-called deshomogenization method, introduced by O. Pantz and K. Trabelsi in [START_REF] Pantz | A post-treatment of the homogenization method for shape optimization[END_REF][START_REF]Construction of minimization sequences for shape optimization[END_REF], then used in [START_REF] Allaire | Topology optimization of modulated and oriented periodic microstructures by the homogenization method[END_REF][START_REF] Geoffroy-Donders | Optimization of oriented and parametric cellular structures by the homogenization method[END_REF][START_REF] Groen | Homogenization-based topology optimization for high-resolution manufacturable microstructures[END_REF]. A composite homogenized structure is projected onto a classical shape that resembles the original microstructure at a macroscopic (yet thin) scale. One great benefit of this approach is its computational efficiency, since fine enough details can be produced even when the shape optimization stage takes place on a relatively coarse mesh.

Shape derivative of the optimized criterion

In this section, we provide the shape derivative of the compliance C(Ω) defined by (3.3), in the general d-dimensional setting. For the sake of completeness, a formal proof, relying on the method of Céa [START_REF] Céa | Conception optimale ou identification de formes: calcul rapide de la dérivée directionnelle de la fonction coût[END_REF] is presented in Appendix C. So as to deal with the various possibilities of Section 5 about the nature of A Ω in a unified fashion, we consider a generic Hooke's tensor of the form:

A Ω (x) = A(d Ω (x), n(x)),
C (Ω)(θ) = Γ ray ∂Ω (x ) q Ω (x, d Ω (z)) ∂A ∂s (d Ω , n)e(u Ω ) : e(u Ω ) (z) d (z) (θ • n)(x) ds(x) - Γ div ∂Ω ray ∂Ω (x ) Q Ω (x, d Ω (z)) ∂A ∂n (d Ω , n)e(u) : e(u) (z) d (z) (θ • n)(x)ds(x) + Γ ray ∂Ω (x ) Q Ω (x, d Ω (z)) ∂A ∂n (d Ω , n)e(u) : e(u) (z) d (z) • n κ(x)(θ • n)(x) ds(x),
where d denotes the line integral over the sets ray ∂Ω (x) for x ∈ ∂Ω, div ∂Ω v := divv -(∇vn) • n is the tangential divergence of a smooth enough vector field v : ∂Ω → R d , and for x ∈ ∂Ω and s ∈ R, the quantities q Ω (x, s) and Q Ω (x, s) are defined by:

(7.2) q Ω (x, s) = d-1 i=1
(1 + sκ i (x)), and Q Ω (x, s) = -I + sII Ω (x)(I + sII Ω (x)) -1 .

In (7.2), II Ω (x) is the d × d matrix of the second fundamental form of ∂Ω at x ∈ ∂Ω; its eigenvalues are the principal curvatures κ i (x), i = 1, ..., d -1 of ∂Ω at x (the associated eigenvectors lie in the tangent plane to ∂Ω at x) together with the value 0, whose associated eigenvector is n(x); eventually, κ := κ 1 + ... + κ d-1 is the mean curvature of ∂Ω.

Remark 7.1.

• When the 'crust-pattern' model is used, an analysis in the spirit of [START_REF] Dapogny | Optimisation de formes, méthode des lignes de niveaux sur maillages non structurés et évolution de maillages[END_REF], Chap. 4, reveals that q Ω (x, s) ≈ 1, Q Ω (x, s) ≈ -I and the first term in the right-hand side of (7.1) is of the order O(1) as the half-thickness ε of the crust tends to 0, whereas the last two terms there are of the order O(ε). These may therefore be neglected in (7.1), thus alleviating the burden of calculating numerically the tangential divergence and the mean curvature κ. These simplifications are consistently used in the numerical results of Section 8. • We refer to [START_REF] Dapogny | Geometric constraints for shape and topology optimization in architectural design[END_REF] for the shape derivative of the maximum thickness constraint functional (6.1).

Numerical results

We eventually report on several numerical experiments about the impact of the models of Section 5 for the construction of shapes, and about the two-step numerical strategy of Section 6 for the non regular infill design.

Details about the numerical implementation.

In the examples of Section 8.2 concerning the influence of the 'crust-pattern' and 'offset' models on the elastic properties of structures, the computational domain D is equipped with a Cartesian grid; the level set function φ accounting for shapes (see Section 3.3) is discretized as a Q 1 Finite Element function (in particular, its values are stored at the vertices of the grid); the Finite Element analyses involved in the resolution of the linearized elasticity equations (3.1) and the level set calculations are also performed on the grid, using a Scilab implementation [START_REF] Campbell | Modeling and Simulation in SCILAB[END_REF].

As far as the examples of Section 8.3 are concerned, the first stage of the proposed methodology consists in the resolution of a shape optimization problem similar to those tackled in Section 8.2; the numerical framework used to this end is identical. The second stage of the process takes place in a triangular mesh of the resulting design: the shape optimization problem involved there is solved within the FreeFem++ environment [START_REF] Pironneau | [END_REF], using P 1 Finite Elements when it comes to mechanical analyses; the signed distance function is calculated by using an algorithm from our previous work [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF] and the level set Hamilton-Jacobi equation (3.7) featured in the Level Set method is solved using an algorithm based on the method of characteristics [START_REF] Bui | An accurate anisotropic adaptation method for solving the level set advection equation[END_REF].

In all cases, the constrained optimization problems under scrutiny are solved with an SLP-type algorithm, similar to the one presented in [START_REF] Dunning | Introducing the sequential linear programming level-set method for topology optimization[END_REF].

One feature of our implementation is not totally classical and deserves a little attention, namely the way we calculate integrals of the form

I := ray ∂Ω (x) w(z)d (z),
for a given function w : D → R (in practice, w is generally defined as a constant per element of the grid of D, see (7.1)), and a given point x ∈ ∂Ω, identified as a point x ∈ D belonging to an element cut by the 0 isoline of the level set function φ for Ω or by that of its signed distance function d Ω .

• The set ray ∂Ω (x) is the segment with direction ∇φ(x) = ∇d Ω (x) and whose endpoints lie on the skeleton Σ of Ω (which is practically detected by a change of monotonicity of d Ω along the line with direction ∇φ(x)). • A quadrature rule, using points {p i } i=1,...,N , is chosen for integrating along this segment;

• The integrand w of I is interpolated at the nodes of the computational grid of D, then at the points p i , and I is calculated according to the selected quadrature rule.

Let us finally discuss the numerical values used in our examples. In all cases, the isotropic reference material A ref is characterized by E ref = 1.93 and ν ref = 0.2. The parameter η for the Ersatz material in (5.9) equals 10 -3 . When the mechanical properties of shapes are described according to the 'crust-pattern' model, the orthotropic behavior of the crust is based on experiments performed in the SIMaP laboratory, and we use the following Hooke's tensor: In particular, the crust is almost twice stiffer in the tangential direction, which is in good agreement with experimental measurements reported in [START_REF] Ahn | Anisotropic material properties of fused deposition modeling abs[END_REF].

Unfortunately, this numerical framework does not allow us to use very fine numerical grids, out of intrinseque limitations of the Scilab environment in terms of memory management (similar limitations related to Matlab were reported in [START_REF] Andreassen | Efficient topology optimization in matlab using 88 lines of code[END_REF]). Such fine grids would be necessary for an accurate modelling of a thin crust in the examples of Section 8.2, or for the creation of a non-regular infill design with fine details in the discussions of Section 8.3. As we have already mentionned, this is the main reason why the presented examples in this work are only two-dimensional.

Comparison of optimized shapes using different material models.

In this section, we conduct experiments about the influence on optimized shapes of the 'crust-pattern' and 'offset' models proposed in Section 5 for describing their material properties.

To ease the visualisation of the results, the different regions of optimized designs are displayed with textures depending on their elastic properties, as depicted on Figure 9.

Two-dimensional cantilever.

Our first example is the classical two-dimensional cantilever, enclosed in a rectangular computational domain D with size 1 × 0.5, meshed with 160 × 80 Q 1 elements. The shape is clamped on its left-hand side and a vertical load g = (0, -1) is applied at the middle of its right-hand side (see Figure 10). The mass of the structure is minimized under a constraint on its compliance, i.e. we solve: (f) The 'crust-pattern' model is used with an anisotropic crust and an isotropic bulk.

To better assess the influence of each model on material properties, the values of the elastic coefficients A 1111 and A 2222 of the material filling the shape of Figure 10 (right) are displayed on Figures 11 and12 respectively, in the above six cases (notice that the clamping region Γ D is not considered to be part of the crust). In particular, the various regions of the computational domain (bulk, crust and void) and the principal directions of the elasticity tensors characterizing them are easily identified. cases, the compliance constraint is active at convergence: the compliance of all the optimized designs equals α c . The CPU time of one such optimization process is approximately 40 min on a standard laptop.

The obtained results suggest that when the crust region of the shape is endowed with anisotropic properties, the algorithm tries to take advantage of the stiffness of the crust in the tangential direction: near the region where loads are applied, the structural boundary tends to align with the direction of the load; the portion of the boundary which is aligned with the load is wider as the infill region is less dense (and so less efficient from the mechanical point of view). On a different note, as the density of the infill region decreases, shapes 'inflate' so as to increase their moment of inertia. This especially happens when an isotropic infill is used, since for a given value of the density ρ cp bulk , the Young's modulus of the isotropic infill defined in Section 5.3 is weaker than the maximum Young's modulus of its anisotropic counterpart; in the present case, the latter is associated to the horizontal direction (since τ 0 = e 1 ), which is close to the principal stress direction near the clamping region Γ D (recall that the upper Voigt-Reuss bound is stiffer than the upper Hashin-Shtrikman bound). Let us now take a closer look at the mechanical efficiency of the obtained designs. Our comparisons are based on the ratio

(8.3) r(Ω) := M (Ω) M (Ω * )
between the mass of the various optimized shapes Ω and that of the optimized shape Ω * in the 'molded' case. This ratio is presented on Figure 17 in the cases (c-f); in the case (b) of the 'offset' model, its value is r(Ω) = 3.46. As a general trend, for a given, common infill design (be it isotropic or anisotropic), shapes having an isotropic crust perform better than the ones with an anisotropic crust. This is to be expected, since the latter are almost as stiff in the tangential direction to the structural boundary, and twice weaker in the normal direction. Moreover, comparing cases (c) and (d) or (e) and (f) suggests that for a given mechanical behavior of the crust (isotropic or anisotropic), those with an anisotropic infill perform better than those with an isotropic hatch. Again, this is justified by the facts that the maximum Young's modulus predicted by the Voigt-Reuss bounds used for the anisotropic modelling of the bulk is stiffer than the Young's modulus given by the Hashin-Shtrikman bounds used for the isotropic modelling, and that the maximum Young's modulus in the presented case featuring an anisotropic bulk is attained in the horizontal direction, which is aligned with the principal stress direction of the structure close the the clamping region. Last but not least, let us remark that the optimized design of Figure 13 (right) obtained using the 'offset' model (b) has a dramatically poor performance in this example, even when compared with structures obeying the 'crust-pattern' model with a much lower density for the infill region. This is why we have only considered the case ρ off bulk = 1 in this example, whereas lower values could be used in principle. 8.2.2. Two-dimensional L-shaped beam.

Our second example is a two-dimensional L-shaped beam, clamped at its upper side and subjected to a vertical load g = (0, -1) applied at the middle of its right-hand side; see Figure 18 (left).

Our purpose is to compare the following six models for the construction of shapes: (a) The 'molded' case, where the shape is filled with the isotropic reference material A ref ;

(b) The 'offset' model; (c) The 'crust-pattern' model with an isotropic crust and an isotropic bulk;

(d) The 'crust-pattern' model with an anisotropic crust and an isotropic bulk;

(e) The 'crust-pattern' model with an anisotropic crust and an anisotropic bulk with horizontal rasters (τ 0 = e 1 with the notation of Section 5.1.2); (f) The 'crust-pattern' model with an anisotropic crust and an anisotropic bulk with vertical rasters (τ 0 = e 2 ). Starting from the initial design of Figure 18 (right), we first solve (8.2) for α c = 60 in the 'molded' case, i.e. by assuming isotropic mechanical properties A ref inside shapes. This yields the optimized shape Ω * of Figure 19 (left) whose mass M (Ω * ) equals 0.21.

Starting from Ω * , the optimized shape using the 'offset' model (b) with ρ off bulk = 1 is displayed on Figure 19 (right), and those associated to the 'crust-pattern' models (c-f) with densities ρ off bulk = 0.90, 0.75 and 0.55 In Figure 23, the mass ratio r(Ω) (see (8.3)) of the optimized shapes using the 'crust-pattern' models (c-f) are reported; in the case of the 'offset' model (b), which is not represented on that graph, we get r(Ω) = 2.07.

The results obtained in this section partially confirm the trends observed in the previous example. One observes that the vertical bars close to the clamping region are thicker when the anisotropic infill is composed of vertical rasters, whereas when it shows horizontal rasters, more material is put at the lower part. Grossly speaking, in order to fulfill the compliance constraint featured in (8.2), the shape concentrates more material in regions where rasters are aligned with the principal stress directions. When it comes to mechanical efficiency, the shapes obtained using models with anisotropic infills (cases (e) and (f)) lead to designs with worse performance, in sharp contrast with the results of Section 8.2.1. The reason is that, in the present example, the direction of principal stresses change significantly from one part of the structure to the other. Therefore, both horizontal and vertical orientations for the rasters composing the infill region prove dramatically inefficient in this case (their inefficiency coming from different regions of the structure). Finally, one notices that some designs for model (c) seem to perform slightly better than the reference 'molded' case (i.e. r(Ω) < 1.0) when ρ bulk > 0.7. This fact may either be attributed to the mechanical fact that shape and topology optimization problems generally favor 'grayscale', intermediate regions between material and void over clear 'black-and-white' designs (see again [START_REF] Allaire | Shape optimization by the homogenization method[END_REF] for instance), or to the numerical approximations used for the mass computation; see (5.12). 8.2.3. Optimization of a two-dimensional MBB Beam.

Our third example is a two-dimensional MBB-beam with size 6 × 1, as depicted in Figure 24: the structure is clamped on its lower left corner, and its displacement in the vertical direction is prevented at the lower Again, we remark that when an anisotropic model is used for the crust region, the structural boundary tends to get aligned with the direction of principal stresses (which in this case does not coincide with the direction of the loads where they are applied). Besides, the structural boundary is also aligned with the support reaction force (the two clamping corners).

Regarding now the performance of optimized designs, one observes that the most efficient design (i.e. that with the lowest mass) does not correspond to the 'molded' case (a), but to the case where an anisotropic infill is used with horizontal rasters (case (d)). This is not so surprising, since in this example, shapes are submitted to pure flexion in the interval between the two applied loads, and so an anisotropic infill with stronger stiffness in the horizontal direction is expected to perform well. Also, in case (d), the relation between density and performance is non-monotonic. This confirms that the performance of an optimized printed part is subtly linked to both the modelling and the density of the bulk region, and that the choice of a raster pattern for printing the infill region of structures cannot be decoupled from the particular application.

8.3. Optimization of the non-regular infill design.

In this section, we evaluate the approach presented in Section 6 for the optimization of the non regular infill of structures.

8.3.1. Non-regular infill for the two-dimensional cantilever.

Our first example takes place in the setting of the two-dimensional cantilever of Section 8.2.1; see again Figure 10 (left) for the details of the test case.

As we have hinted at in Section 6, there are mainly three possibilities for the optimal design of a shape and its non regular infill: (1) optimizing only the infill region and revealing the exterior of the structure as the region with vanishing infill density, (2) optimizing jointly the exterior shape and the infill design, which are described as two nested shapes, and (3) optimizing in a first step the exterior design (assuming for instance a homogeneous, isotropic infill), then optimizing the infill design in a second step. Our preliminary aim in this section is to illustrate why the first of these processes is not retained. Leaving the second one aside owing to its implementation complexity, we shall then examine the third one with further scrutiny. Starting from the design of Figure 27 (a), we solve (8.2) on a Cartesian grid with size 400 × 200, assuming shapes to be filled with a fully dense isotropic material and using the value α c = 40. This yields the optimized shape of Figure 27 (b).

Let us consider the design of a structure with a non-regular infill; as discussed in Section 6, to achieve this in our setting, we add a constraint on the maximum thickness of shapes to (8.2); starting again from the design of Figure 27 (a), we solve:

(8.4) min Ω M (Ω), s.t. C(Ω) ≤ α c , P MT (Ω) ≤ δ 2 , for δ = 10∆x
, where ∆x stands for the uniform mesh size. This 'one-step' procedure results in the optimized shape Ω os of Figure 27 (c). Admittedly, a design with an organic-like infill region has been achieved, but the exterior shape barely differs from the boundary of the computational domain D.

Let us now evaluate the proposed two-step approach proposed in Section 6 for the non regular infill design; at first, we solve (8.2) assuming shapes with an isotropic crust and an isotropic uniform infill, with density ρ cp bulk = 0.6. This yields the structure of Figure 16 (left). Then, as a second stage, we search to trade the uniform infill of this shape for a non-regular infill with better performance.

To this end, at the end of the first step, the Level Set function is defined on a uniform regular grid. Using these values, one constructs a mesh fitted to the actual shape boundary. Different approaches are available to deal with this issue; see for instance [START_REF] Allaire | Shape optimization with a level set based mesh evolution method[END_REF][START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF][START_REF] Xia | A level set solution to the stress-based structural shape and topology optimization[END_REF]. Applying this methodology in our example results in the new working domain of Figure 28 (right).

Then, to obtain non-regular infill structures, we start from the quite arbitrary initialization ω 0 of Figure 29 (left) and solve the optimization problem: 1. We observe that the resulting shapes from this process are significantly lighter than the initial infill design ω 0 for the same mechanical performance. Also, the shapes resulting from this two-step optimization process generally perform better than that Ω os obtained by directly solving (8.4).

Case

Reference design ω 0 The same two-step approach is applied in the case of the L-shaped beam of Section 8.2.2. In a first time, assuming shapes with an isotropic crust and an isotropic infill with density ρ cp bulk = 0.55, the mass of shapes is minimized under a constraint on their compliance; i.e. Problem (8.2) is solved using α c = 0.6. This results in the optimized shape of Figure 30 (left). Thence, a triangular mesh of this structure is obtained, which is represented on Figure 30 (right). Starting from the arbitrary initialization ω 0 of Figure 31 (left), we solve the new optimization problem (8.4) for the infill structure ω ⊂ Ω * , for several values of the parameter δ for the maximum allowed thickness of the infill design. This results in the shapes displayed on Figure 31. The data associated to the computation are reported in Table 2, and they lead to the same conclusions as in the cantilever example of Section 8.3.1.

Ω os δ = 3∆x δ = 4∆x δ = 5δx M (Ω) 0.
Remark 8.1. In both Sections 8.3.1 and 8.3.2, the performance of the resulting infill designs ω does not decrease monotonically as the thickness constraint is strenghtened (see Tables 1 and2). Let us emphasize that maximum thickness constraints may be modelled and implemented in a wide variety of manners, and that they generally induce the existence of many local minima in the associated shape optimization problem (see e.g. [START_REF] Michailidis | Manufacturing Constraints and Multi-Phase Shape and Topology Optimization via a Level-Set Method[END_REF]). Our conclusions may therefore be biased by our implementation choices. We now recall the expression of the compliance and Hooke's tensors in the two-dimensional plane stress situation. Being symmetric, σ(u) and e(u) have in total 3 independent components; in practice, once an orthonormal frame of R 2 is chosen, they are represented by three-dimensional vectors, using Voigt's notation: In this article, we consider a particular class of elastic materials, called orthotropic (or sometimes orthorombic): their mechanical properties are symmetric with respect to 2 orthogonal planes, which entails simplifications in the expressions (A.2) of the general, fully anisotropic tensor A and S: denoting by by e 1 , e 2 the corresponding unit normal vectors -or principal directions -, then expressing e(u) and σ(u) in the direct orthonormal frame ( e 1 , e 2 ) of R 2 , the compliance tensor S is of the form:

(A.3) S =   1 E1 -ν21 E2 0 -ν12 E1 1 E2 0 0 0 1 µ   .
In this formula,

• E i (i = 1, 2) is the Young's modulus in direction e i ; it is a measure of the resistance of the material to a compression in this direction, • the ν ij (i, j = 1, 2) are the Poisson's ratios; they appraise the decline in the direction e j of a piece of material which is stretched in the direction e i , • µ is the shear modulus; it describes the response of a portion of material when a deformation in the e 1 -direction is applied to the face oriented by e 2 . Recall in addition that in the particular case of isotropic materials, all the Young's moduli (resp. the Poisson's ratios) share the same value E (resp. ν), and the shear modulus µ equals µ = E 2(1+ν) , so that the Hooke's law A rewrites under the simple form (again, in 2d): Let (e 1 , e 2 ) be the canonical basis of R 2 , and let ( e 1 , e 2 ) be a reference basis (i.e. a basis in which the material properties take a simple form such as (A.3)). We assume that ( e 1 , e 2 ) is obtained from (e 1 , e 2 ) by a rotation of angle α, whose matrix in the canonical basis of R 2 is: R(α) = cos αsin α sin α cos α .

If the displacement u : R 2 → R 2 of a structure is expressed in the canonical basis, the corresponding displacement u : R 2 → R 2 in the rotated basis reads: Notice that an easy calculation reveals that R e (α) -1 = R e (-α), and likewise for R σ (α). Eventually, the relation between the compliance tensors S and S (resp. the Hooke's tensors A and A) in the rotated and canonical bases reads: S = R e (-α) S R σ (α) and A = R σ (-α) A R e (α).

B.2. Expression of the elasticity tensor in the local basis (τ, n) of the plane.

In this section, we now assume that the reference basis ( e 1 , e 2 ) for the properties of the elastic material is the local basis (τ, n) of the plane defined by the tangent and normal vectors to the boundary ∂Ω of a structure Ω at a given point x ∈ ∂Ω.

Letting n = (n 1 , n 2 ), the rotation matrix of the mapping from the basis (e 1 , e 2 ) into (τ, n) reads:

R(n) = n 2 n 1 -n 1 n 2 .
Like in the previous section, we then obtain, using Voigt's notation: We provide in this appendix a formal calculation of the shape derivative featured in Theorem 7.1, based on Céa's method [START_REF] Céa | Conception optimale ou identification de formes: calcul rapide de la dérivée directionnelle de la fonction coût[END_REF] (see also [START_REF]Conception optimale de structures[END_REF]). Let us define the Lagrangian L : So far, the shape derivative (C.4) does not satisfy the convenient structure (3.4). To achieve this feature, some additional technical work is needed, which we only outline, referring to [START_REF] Allaire | Mutli-phase structural optimization via a level-set method[END_REF][START_REF] Dapogny | Optimisation de formes, méthode des lignes de niveaux sur maillages non structurés et évolution de maillages[END_REF][START_REF] Michailidis | Manufacturing Constraints and Multi-Phase Shape and Topology Optimization via a Level-Set Method[END_REF] for full details:

U ad × H 1 Γ D (R d ) d × H 1 Γ D (R d
• For a given point x ∈ D \ Σ, the derivative of the mapping θ → d Ω θ (x) reads:

d Ω (θ)(x) = -θ(p) • n(p), where p = p ∂Ω (x).

• In the same spirit, one has; for x ∈ D \ Σ:

∇d Ω (θ)(x) = -I + d Ω (x)II ∂Ω (p)(I + d Ω (x)II ∂Ω (p)) -1 ∇ ∂Ω (θ • n)(p), p = p ∂Ω (x).

• We then use the coarea formula, as in [START_REF] Allaire | Mutli-phase structural optimization via a level-set method[END_REF], as a 'curved' version of the Fubini theorem, in order to transform integration on D in (C.4) into nested integrations over the surface ∂Ω and on the rays ray ∂Ω (x), x ∈ ∂Ω.

Using these facts, and under the previous assumptions, the shape derivative (C.4) becomes:

(C.5)

C (Ω)(θ) = ∂Ω ray ∂Ω (x )
q Ω (x, d Ω (z)) ∂A ∂s e(u) : e(u) (z) d (z) θ(x) • n(x) ds(x)

+ ∂Ω ray ∂Ω (x )
Q Ω (x, d Ω (z)) ∂A ∂n e(u) : e(u) (z) d (z) • ∇ ∂Ω (θ • n) (x) ds(x), in which d stands for the one-dimensional, line integration, and we have defined q Ω (x, d) and Q Ω (x, d) by (7.2). Eventually, integration by parts on ∂Ω yields (see e. Remark C.1. Although the expression (C.6) has the convenient structure (3.4), allowing for the easy identification of a shape gradient as outlined in Section 3.2, the tangential divergence and mean curvature featured in there are difficult to calculate in practice, leading to large numerical errors. Another possibility to obtain a descent direction for C(Ω) which is not so prone to numerical errors relies on the alternative expression (C.5) for C (Ω)(θ), and on a change of inner products for the identification of a shape gradient; see e.g. [START_REF] Burger | A framework for the construction of level set methods for shape optimization and reconstruction[END_REF][START_REF] Gournay | Velocity extension for the level-set method and multiple eigenvalues in shape optimization[END_REF]. For example, considering vector fields of the form θ = wn for some scalar field w : ∂Ω → R, the shape derivative (C.5) rewrites: Then, searching for V ∈ H 1 (D) such that the following variational problem is satisfied:

∀φ ∈ H 1 (D), D α 2 ∇V • ∇φ + V φ dx = C (Ω)(φ)
for a given scalar parameter α (usually chosen of the order of the mesh size) and letting w = -V in (C.7) guarantees a descent direction for C(Ω) since

C (Ω)(-V ) = - D α 2 |∇V | 2 +V 2 dx < 0.

Figure 1 .

 1 Figure 1. Rough sketch of the slicing procedure, initiating the additive manufacturing processes.

Figure 2 .

 2 Figure 2. Construction of one 2d layer of a 3d structure with the FFF technology. (Topleft) outline of the layer; (top-right) construction of the contour;(bottom, from the left to the right) assembly of the bulk of the layer using an isotropic infill, an anisotropic (horizontally oriented) infill, along contour offsets, and using a personalized infill.

Figure 3 .

 3 Figure 3. Examples of two-dimensional layers constructed according to (left) the 'crustpattern' model, and (right) the 'offset' model (Photo courtesy of Richard Horne † ).

Definition 4 . 1 .

 41 Let Ω ⊂ R d be a bounded Lipschitz domain. The signed distance function d Ω : R d → R to Ω is defined by:

  where

  is the usual, Euclidean distance from a point x ∈ R d to the boundary ∂Ω.

Figure 5 .

 5 Figure 5. (Left) Projection point p ∂Ω (x) of x ∈ Ω, projection set Π ∂Ω (x ) = {y 1 , y 2 } of another point x ∈ Ω and ray emerging from z ∈ ∂Ω; (right) the gradient of d Ω (x), if it is defined, equals the normal vector to ∂Ω at p ∂Ω (x).

Definition 4 . 3 .

 43 The domain Ω ⊂ R d has maximum thickness lower than δ > 0 if:

Figure 6 .

 6 Figure 6. Interpolation of any (scalar) entry P of the tensor A cp Ω using the 'crust-pattern' model (5.2).

Figure 7 .

 7 Figure 7. Interpolation of any of the (scalar) entries P of the Hooke's tensor A off Ω using the 'offset' model (5.10).

3 D x 3 Figure 8 .

 338 Figure 8. In three space dimensions, the natural frame for the orthotropy of the material in the crust region is (τ H , n H , e 3 ).

Figure 9 .

 9 Figure 9. Schematic representation used in the presentation of optimized shapes in Section 8.2.

Figure 10 . 1 .

 101 Figure 10. (Left) Setting and (right) initial design in the two-dimensional cantilever test case of Section 8.2.1.Based on the discussions of Section 5, we compare six situations for the elastic properties of shapes:(a) Shapes are filled with the (isotropic) reference material A ref ; since this situation typically arises when structures are constructed by molding, we refer to it as the 'molded' case; (b) The 'offset' model of Section 5.2 is retained; (c) The 'crust-pattern' model of Section 5.1 is used with isotropic (yet different) properties in the crust and bulk regions; (d) The 'crust-pattern' model is used, with still an isotropic crust, and an anisotropic bulk filled with horizontal rasters of material (τ 0 = e 1 , with the notations of Section 5.1.2); (e) The 'crust-pattern' model is used with both anisotropic crust and bulk (τ 0 = e 1 for the latter);

Figure 11 .Figure 12 .

 1112 Figure 11. Distribution of the elastic coefficient A 1111 using the six considered models for the elastic properties of shapes in the cantilever test case of Section 8.2.1.

Figure 13 .

 13 Figure 13. Optimized shapes Ω * and Ω for the two-dimensional cantilever test case of Section 8.2.1 in (left) the 'molded' case (a), and (right) the 'offset' model (b).

Figure 14 .

 14 Figure 14. Optimized shapes in the cantilever example of Section 8.2.1 obtained with the four 'crust-pattern' models considered for the assembly of shapes, using an infill density ρ cp bulk = 0.90.

Figure 15 .

 15 Figure 15. Optimized shapes in the cantilever example of Section 8.2.1 obtained with the four 'crust-pattern' models considered for the assembly of shapes, using an infill density ρ cp bulk = 0.75.

Figure 16 .

 16 Figure 16. Optimized shapes in the cantilever example of Section 8.2.1 obtained with the four 'crust-pattern' models considered for the assembly of shapes, using an infill density ρ cp bulk = 0.60.

Figure 17 .

 17 Figure 17. Performance of the optimized cantilevers of Section 8.2.1 for different densities and mechanical modelings.

D • g 1 2Figure 18 .

 118 Figure 18. (Left) Setting and (right) initial shape in the L-shaped beam test case of Section 8.2.2.

Figure 19 .

 19 Figure 19. Optimized shapes Ω * and Ω obtained in the L-shaped beam test case of Section 8.2.2 (left) in the 'molded' situation, (right) when the material properties of shapes are described using the 'offset' model.

Figure 20 .

 20 Figure 20. Optimized shapes in the L-beam example of Section 8.2.2 obtained with the four considered 'crust-pattern' models for the assembly of shapes, using the infill density ρ cp bulk = 0.90.

Figure 21 .

 21 Figure 21. Optimized shapes in the L-beam example of Section 8.2.2 obtained with the four considered 'crust-pattern' models for the assembly of shapes, using the infill density ρ cp bulk = 0.75.

2 ,

 2 . C(ω) ≤ C(ω 0 ), P MT (ω) ≤ δ This problem is solved for various values δ of the allowed maximum thickness; see Figure29. The data of these calculations are reported in Table

Figure 22 .

 22 Figure 22. Optimized shapes in the L-beam example of Section 8.2.2 obtained with the four considered 'crust-pattern' models for the assembly of shapes, using the infill density ρ cp bulk = 0.55.

Figure 23 . 1 6Figure 24 .

 23124 Figure 23. Performance of the optimized L-shaped beams of Section 8.2.2 for different densities and mechanical modelings.

Figure 25 .

 25 Figure 25. Optimized shapes in the MBB beam example of Section 8.2.3 obtained with the four considered models for the assembly of shapes, using the infill density ρ cp bulk = 0.6.

Figure 26 .

 26 Fig. 25

Figure 27 .Figure 28 .

 2728 Figure 27. (a) Initial design, and (b) optimized two-dimensional cantilever obtained in Section 8.3.1 without imposing any thickness constraint (i.e. (8.2) is solved); (c) optimized two-dimensional cantilever for (8.4) with δ = 10∆x. where A : S d (R) → S d (R) is the Hooke's tensor. From the physical point of view, this behavior is often reflected by the inverse relation: e(u) = Sσ(u),

Figure 29 .

 29 Figure 29. (a) initial infill design ω 0 for the second stage of the two-step process for non regular infill design in the cantilever example of Section 8.3.1; (b-d) optimized shapes ω having the same compliance as ω 0 setting (b) δ = 3∆x, (c) δ = 4∆x and (d) δ = 5∆x.

Figure 30 .A

 30 Figure 30. First step of the two-step process for non regular infill design in the L-beam example of Section 8.3.2; (left) optimized shape considering shapes with an isotropic crust and an isotropic infill with density 0.6; (Right) exported unstructured mesh, used as a new design domain for the second step.

Ae

  = 2µe + (Bµ)tr(e)I, where B := E 2(1 + ν)(1 -2ν) is the bulk modulus of the material. Appendix B. Change of orthonormal bases and elasticity tensors in two space dimensions B.1. Rotation of elasticity tensors.

Figure 31 .

 31 Figure 31. (a) initial infill design ω 0 for the second stage of the two-step process for non regular infill design in the L-beam example of Section 8.3.2; (b-d) optimized shapes ω having the same compliance as ω 0 setting (b) δ = 4∆x, (c) δ = 5∆x and (d) δ = 8∆x.

u cos 2 α sin 2 α -sin α cos α sin 2 α cos 2 α cos 2 α sin 2 α -2 sin α cos α sin 2 α cos 2 α 2

 22 ( x) = R(α)u(x), and ∇ u( x) = R(α)∇u(x)R(α)T , where x = R(α)x.Hence, the strain tensor e in the basis ( e 1 , e 2 ) is connected to its counterpart e in the canonical basis (e 1 , e 2 ) via the following formula: e = R(α)eR(α) T , and a similar transformation rule applies for the stress tensors σ and σ. Expressing this identity using Voigt's notation (A.1) yields: (B.1) e = R e (α)e, where R e (α) :=  sin α cos α 2 sin α cos α -2 sin α cos α cos 2 αsin 2 α σ = R σ (α)σ, where R σ (α) :=  sin α cos α sin α cos αsin α cos α cos 2 αsin 2 α   .

(B. 3 ) 2 -2n 1 n 2 -n 1 n 2 n 1 n 2 n 2

 3222 e = R e (n)e, where R e (n) := σ = R σ (n)σ, where R σ (n) := Sketch of the proof of Theorem 7.1

  ) d by: (C.1) L(Ω, u, p) = D A(d Ω , ∇d Ω )e(u) : e(u) dx + D A(d Ω , ∇d Ω )e(u) : e(p) dx -D f • p dx -Γ N g • p ds.For a fixed shape Ω ∈ U ad , we now search for the critical points (u, p)∈ H 1 Γ D (R d ) d × H 1 Γ D (R d ) d of the functional L(Ω, •, •), and to this end, we equate its partial derivatives to 0:• At first, the relation∀φ ∈ H 1 Γ D (R d ) d , ∂L ∂p (Ω, u, p)(φ) = 0, characterizes u as the unique solution u Ω of (3.1), in which A Ω is replaced by A(d Ω , ∇d Ω ).• Likewise, cancelling the partial derivative of L(Ω, •, •) with respect to u at (u, p) yields:(C.2) ∀φ ∈ H 1 Γ D (R d ) d , 2 D A(d Ω , ∇d Ω )e(u) : e(φ) dx + D A(d Ω , ∇d Ω )e(φ) : e(p) dx = 0, whence we readily identify p = -2u Ω .We now rely on the fact that, for arbitrary φ ∈ H 1 Γ D (R d ) d , we have by definition:C(Ω) = L(Ω, u Ω , φ);differentiating with respect to Ω in the above relation and using chain rule, we therefore obtain:(C.3) C (Ω)(θ) = ∂L ∂Ω (Ω, u Ω , φ)(θ) + ∂L ∂u (Ω, u Ω , φ)(u Ω (θ)),where u Ω (θ) is the shape (or Eulerian) derivative of Ω → u Ω . Now, since (C.3) holds for any φ ∈ H 1 Γ D (R d ) d , we may take φ = -2u Ω so that the last term in the right-hand side vanishes (see (C.2)), leaving:C (Ω)(θ) = ∂L ∂Ω (Ω, u Ω , p Ω )(θ),and we are left with the calculation of the latter (partial) shape derivative.To this end, we remark that the only explicit dependence with respect to the shape Ω in the expression (C.1) of L is in the elasticity tensor A (d Ω , ∇d Ω ); it follows: (C.4) C (Ω)(θ) = -D ∂A ∂Ω (θ)e(u) : e(u) dx = -D ∂A ∂s e(u) : e(u) d Ω (θ) dx -D ∂A ∂n e(u) : e(u) • ∇(d Ω (θ)) dx.

  g.[START_REF] Henrot | Variation et optimisation de formes: une analyse géométrique[END_REF], Prop. 5.4.9):(C.6) C (Ω)(θ) = ∂Ω ray ∂Ω (x ) q Ω (x, d Ω (z)) ∂A ∂s e(u) : e(u) (z) d (z) (θ • n)(x)ds(x) -∂Ω div ∂Ω ray ∂Ω (x ) Q Ω (x, d Ω (z)) ∂A ∂n e(u) : e(u) (z)d (z) (θ • n)(x) ds(x) + ∂Ω ray ∂Ω (x ) Q Ω (x, d Ω (z)) ∂A ∂n e(u) : e(u) (z) d (z) • n(θ • n)(x) ds,which is the desired expression.

(C. 7 )

 7 C (Ω)(wn) = ∂Ω ray ∂Ω (x )q Ω (x, d Ω (z))∂A ∂s e(u) : e(u) (z) d (z) w(x) ds(x)+ ∂Ω ray ∂Ω (x ) Q Ω (x, d Ω (z)) ∂A ∂n e(u) : e(u) (z) d (z) • ∇ ∂Ω w(x) ds(x),

  for some mapping A :R s × R d n → L(S d (R), S d (R)).One word about notations is in order: we denote by ∂A ∂s (s,n) ∈ L(S d (R), S d (R)) and ∂A ∂n (s, n) ∈ L(S d (R), S d (R))d the partial derivatives of A with respect to the variables s and n respectively, at a particular point (s, n) ∈ R× R d . Thus, for arbitrary e ∈ S d (R), ∂A Theorem 7.1. The compliance C(Ω) is shape differentiable at any admissible shape Ω ∈ U ad ; its shape derivative in a direction θ ∈ Θ ad reads:(7.1)

∂s (s, n)e : e ∈ R is the derivative of the mapping R s → A(s, n)e : e ∈ R and ∂A ∂n (s, n)e : e ∈ R d is the gradient of R d n → A(s, n)e : e ∈ R.

Table 1 .

 1 Values of the mass and of the maximum thickness for the optimized cantilevers with non regular infill obtained in Section 8.3.1.

		180	0.172	0.165	0.172	0.165
	P MT (Ω)	-	0.0538	0.0183	0.0248	0.0301

8.3.2. Non-regular infill for the two-dimensional L-shaped beam.

Table 2 .

 2 Values of the mass and of the maximum thickness for the optimized L-shaped beams with non regular infill obtained in Section 8.3.2.

	M (Ω)	0.234	0.178	0.175	0.181	0.175
	P MT (Ω)	-	0.0237	0.0331	0.0402	0.0515

† http://richrap.blogspot.fr

β if 0 ≤ s ≤ ε, 1 if s > ε,where β > 0 is a user-defined parameter. Grossly speaking, h m and h p are smoothed versions of the characteristic functions of the sets (-∞, -ε) and (ε, +∞) respectively, which converge to their exact, discontinuous counterparts as β tends to 0. In our applications, we consistently used the value β = 1. The interpolation procedure (5.2) has already been presented in[START_REF] Faure | Design of thermoelastic multi-material structures with graded interfaces using topology optimization[END_REF][START_REF] Vermaak | Material interface effects on the topology optimizationof multi-phase structures using a level set method[END_REF] (see also[START_REF] Clausen | Topology optimization of coated structures and material interface problems[END_REF] for non monotonic projection operations in the context of density-based topology optimization), and it is illustrated on Figure6, (right).

Acknowledgements. The authors are very grateful to Prof. G. Allaire and Dr. N. Vermaak for stimulating and insightful discussions. They also warmly thank R. Horne for providing some of the pictures of Section 2. C.D. is partially supported by the IRS-CAOS grant from Université Grenoble-Alpes.

In this short section, we briefly recall the main concepts and notations about linearly elastic structures which are used in this article, emphasizing on the case of two space dimensions (a similar representation being valid in the 3d situation); we refer to [START_REF] Gould | Introduction to linear elasticity[END_REF] for further details.

The properties of a linearly elastic material are characterized by the (linear) relation between the stress tensor σ(u) and the (linearized) strain tensor e(u) = 1 2 (∇u + ∇u T ) associated to a displacement field u : R d → R d : σ(u) = Ae(u),