Charles Dapogny

Alexis Faure

Georgios Michailidis

Grégoire Allaire

Agnes Couvelas

Rafael Estevez

Geometric constraints for shape and topology optimization in architectural design

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Designing shapes in an 'optimal' way is a task of utmost interest in a wide variety of applications. Architects, designers, artists and engineers from various fields and with diverse sensitivities endeavour to create products and structures that suit their purposes in the best possible way, while meeting some specifications at the same time. Consciously or not, they try to solve an 'optimization problem', relying on their knowledge, experience and intuition.

Shape and topology optimization [START_REF] Allaire | Conception optimale de structures[END_REF][START_REF] Bendsoe | Topology optimization: theory, methods and applications[END_REF][START_REF] Christensen | An introduction to structural optimization[END_REF] is a discipline which combines mathematical and computational techniques in order to find the shape of a body or a structure as the solution to an optimization problem. It has aroused a considerable enthusiasm in engineering, where the design criteria usually lend themselves to mathematical formulations, and where engineers are already quite familiar with computational methods. The great headway made in computer science during the last decades and the subsequent soar in computational capacity made it possible to apply such techniques to the typical design problems in that field. Applications ranging from nanoscience to airplanes [START_REF] Sigmund | Topology optimization: from airplanes to nano-optics, Bridging from technology to society[END_REF] have been reported in the literature, where complex problems are tackled, in which intuition and experience are limited. Nowadays, several dedicated commercial softwares exist (OptiStruct, TOSCA, GENESIS, etc.), which are extensively used in industry to accelerate the design process.

On the contrary, shape and topology optimization techniques have hardly ever been applied in architecture. Excluding the fact that scientific computing has only recently started being acknowledged by architects as a valuable help in design, there are mainly three sources of explanation for this fact:

• Architects generally follow construction rules, based on the long-standing intuition and mechanical understanding of the discipline, to apprehend the physical situation of a problem, and the multiple manufacturing constraints at stake. Except in 'exceptional cases', where, for instance, the aesthetic quality of the designed structure is the main priority, or where this structure is 'unique', by nature or by size (e.g. skyscrapers), challenging these construction rules proves costly, lengthy, and eventually inefficient. • Although shape and topology optimization results are often directly exploitable in the field of mechanical engineering, this is untrue when it comes to architecture. Indeed, most of the available methods deal with continuous media, while architects and civil engineers typically simulate their structures using beams, plates, shells, etc. • Perhaps the main point lies in that architectural design is mostly driven by aesthetics and is influenced by the designer's inspiration, which is not always easily amenable to a mathematical formulation. As a result, the architect's intervention is limited to the mechanical inputs of the optimization problem, which has a major drawback: the results can be reproduced by anyone making the same choices, i.e. the notion of personal creation disappears! The main purpose of this article is to propose a shape and topology optimization framework which is oriented towards applications in architecture. Depending on the relevance in the particular situation, the material of this work may serve to improve the mechanical performance of the designed structure, or to influence the optimization process with an information about the personal taste of the architect.

From the mathematical viewpoint, structures are described as linearly elastic shapes, and their mechanical performance is evaluated in terms of their compliance, which is to be minimized under a volume constraint. This minimization uses a numerical representation of shapes via the level set method [START_REF] Allaire | A level-set method for shape optimization[END_REF][START_REF]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]; the clear description of shapes allowed by this technique is pivotal in accounting for their geometric attributes (normal vector, principal curvatures, etc.), as opposed to density-based methods [START_REF] Bendsoe | Topology optimization: theory, methods and applications[END_REF] where the same can be done only at the price of an additional approximate reconstruction of the shape from the material density (see [START_REF]Topology optimization of coated structures and material interface problems[END_REF] and references therein).

In this situation, we introduce three geometric constraints; two of them are aimed at alleviating the prejudicial lack of information related to any user-defined notion of aesthetics in shape and topology optimization problems, while the third one allows among other things to control the visibility of structures, and to penalize the length of their bars. The common point between these shape functionals is that they bring into play the notions of signed distance function and maximum thickness, either in the usual Euclidean setting (as in, e.g. [START_REF] Delfour | Shapes and geometries[END_REF] of [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF]), or in the anisotropic one, which is new to the best of our knowledge. The impact of these constraints is discussed in the perspective of computational mechanics, and especially architectural design, even though they could (and for some of them, already have) find an interest in totally different contexts.

The remainder of this article is organized as follows: in Section 2 we explain the sense in which we intend to apply shape and topology optimization techniques to conceptual architectural design. After discussing the limitations of the existing approaches, we present in Section 3 the shape and topology optimization framework adopted in this work to deal with civil structures: in particular, the mechanical model for structures, the mathematical method for evaluating the shape sensitivity of functions and the numerical representation of shapes are detailed. In Section 4, we recall some fairly classical material about the signed distance function and the unit normal vector to a domain; these notions are then extended to the case where the ambient space is equipped with a (constant), anisotropic metric tensor. In Section 5, we formulate and analyze three shape functionals of a geometric nature, which allow to add information about the personal taste of the user regarding aesthetics into the optimization problem, and to model important constraints on civil structures, e.g. constraints over the visibility of shapes, or over the length of their bars in some particular directions. In Section 6, numerical results are presented in two and three space dimensions to appraise the main features of the proposed shape and topology optimization framework. This article ends with a general conclusion in Section 7, and with a technical appendix where all the mathematical details underlying the main issues of this article are concentrated.

Shape and topology optimization for conceptual architectural design

Motivations of this work.

Admittedly, one should not expect such a universal applicability of shape and topology optimization in the field of architectural design as in that of mechanical parts in high-technology industry: as we have mentioned, the nature of the mechanical simulations involved when dealing with civil structures, and the generally overwhelming amount of manufacturing constraints limit the range of application of shape and topology optimization to exceptional structures, where in particular, the total construction budget is not the main concern. However, there are certainly two purposes for which shape and topology optimization techniques may be a valuable help during the stage of conceptual design: exploring new aesthetic concepts which are mechanically viable, and enforcing geometric constraints on the devised structures (thus accelerating the interactions between architects and engineers).

Let us comment on the first point. Shape optimization of structures may help in several directions to get interesting designs from an aesthetical point of view. Several architects indeed argue that an optimized form, e.g. a body bearing a load with the minimum possible material or maximizing some performance criterion for a certain material volume, enhances its proper aesthetic qualities, as it reveals the trajectory of loads in the structure, eliminating any underperforming part [START_REF] Vazquez | Design, analysis and fabrication of expressive, efficient shell structures: a prototype exploring synergy between architecture, engineering and manufacture[END_REF]. The visual impact of such a shape acts in an educational manner, helping the observer to understand how the structure works mechanically.

Shapes obtained as the result of a shape and topology optimization procedure are frequently characterized by a strong organic nature [START_REF] Dombernowsky | Three-dimensional topology optimisation in architectural and structural design of concrete structures[END_REF]. This is not a coincidence, since such shapes exhibit an optimal mechanical response to the applied loads, in a similar way as those found in nature. One famous example of shape optimization in nature concerns the evolution of trees [START_REF] Baumgartner | Sko (soft kill option): the biological way to find an optimum structure topology[END_REF][START_REF] Mattheck | Design and growth rules for biological structures and their application to engineering[END_REF], which is guided by the necessity to withstand wind and to maximize their capacity to collect nutrients. Another well-known example is bone remodelling [START_REF] Wolff | The law of bone remodelling[END_REF], i.e. the process under which bones change their shape and constitution according to the loads they have to bear. Although forces in nature act in a complex way, the principle remains the same [START_REF] Sasaki | Morphogenesis of flux structure[END_REF] and boils down to the minimization of some potential energy. Therefore, shape and topology optimization methods turn out to be a natural way to mimic the rules of biological growth and produce non-intuitive shapes. Let us mention that famous architects, such as Antoni Gaudi and Frei Otto, have extensively used such models in their works, either by directly taking inspiration from natural shapes or by performing simplified experiments under gravitational loads [START_REF] Larena | Shape design methods based on the optimisation of the structure. historical background and application to contemporary architecture[END_REF].

In addition to this 'natural' aesthetic behavior of mechanically optimal shapes, shape and topology optimization allows to devise optimization criteria to help in improving their geometric features (for instance the curvature of their boundaries), see Section 5.2.

Beyond aesthetic motivations, shape and topology optimization allows to model important geometric constraints, for instance visibility constraints, which are pivotal in the study of structures exhibiting a non-standard mechanical behavior, such as skyscrapers [START_REF] Besserud | Structural emergence: architectural and structural design collaboration at som[END_REF].

All in all, our belief is that both engineers and architects may benefit from such an approach, converging faster in a proposal of high aesthetic quality and firm structural coherence: as hinted at in [START_REF] Dombernowsky | Three-dimensional topology optimisation in architectural and structural design of concrete structures[END_REF], an optimized shape can be regarded as a convergence point between architecture and engineering, helping designers to shift from a typological to a performance-driven design [START_REF] Bradner | Parameters tell the design story: ideation and abstraction in design optimization[END_REF][START_REF] Feringa | An integral approach to structural optimization and fabrication[END_REF].

2.2. Shape and topology optimization studies in architecture: pioneer and modern contributions.

Although modern techniques of optimal design were not available at that time, perhaps the first famous example of shape and topology optimization in architecture dates back to the French architect Gustave Eiffel [START_REF] Sundaram | Gustave eiffel and his optimal structures[END_REF]: indeed, the outline of the Eiffel tower in Paris is similar to the optimal shape of a simple cantilever when it is subjected to gravity and wind loads. Moreover, its topology, ordered over three different levels, is strongly inspired by patterns commonly found in nature.

The first applications of modern shape and topology optimization techniques in architecture have been presented by Sasaki [START_REF] Sasaki | Morphogenesis of flux structure[END_REF]. Sasaki used the so-called BESO (Bi-directional Evolutionary Structural Optimization) method in an international competition for the design of the Florence station in 2002. Grossly speaking, he searched for the design which would achieve maximum mechanical efficiency and minimum use of material, imposing two non-optimizable regions for the slab and the entrance. The resulting optimal design avoids bending forces and achieves an iso-stress distribution in the structure, where members work either in tension or in compression [START_REF] Larena | Shape design methods based on the optimisation of the structure. historical background and application to contemporary architecture[END_REF]. The same method was used in the design of the entrance of the Qatar National Convention Centre, in Doha (see Figure 1). The final design is 250m long by 30m wide by 20m high; it shows a dominant tree-like structure constructed via double tubes supported by steel bars, which is strongly reminiscent of typical outlines obtained in shape and topology optimization. Sasaki further explored geometric optimization techniques in his work, i.e. techniques for optimizing the shape of a structure, starting from an initial design, without affecting its topology. Two famous examples presented in [START_REF] Sasaki | Morphogenesis of flux structure[END_REF] are the Island City Central Park GrinGrin and the Crematorium in Kakamigahara, which achieve optimal mechanical efficiency under stress and displacement constraints. Since then, several studies have relied on shape and topology optimization techniques for the design of structures, most of them employing density-based methods. In [START_REF] Dombernowsky | Three-dimensional topology optimisation in architectural and structural design of concrete structures[END_REF], the authors use a commercial software to optimize pre-fabricated, post-tensioned fiber-reinforced concrete elements. In that work, the optimization problem only considers the mechanical behavior of the structure -the deformation energy is minimized under manufacturing constraints linked to the CNC-milling process -and aesthetic aspects are considered implicitly as a post-processing issue, by imposing a subsequent constraint on the already optimized design, so as to drive it towards an organic shape. Similar techniques are used in [START_REF] Ponginan | Danish team uses hyperworks to prove the value of topology optimization for concrete architectural structures[END_REF] to optimize the topology of concrete shelters, in [START_REF] Vazquez | Design, analysis and fabrication of expressive, efficient shell structures: a prototype exploring synergy between architecture, engineering and manufacture[END_REF] for concrete shell structures, as well as in [START_REF] Feringa | An integral approach to structural optimization and fabrication[END_REF] where the so-called Robotic Hotwire Cutting (HWC) technique is selected for processing the structure. Finally, in [START_REF] Frattari | THE STRUCTURAL FORM: Topology Optimization in Architecture and Industrial Design[END_REF] three examples of topology optimization using a commercial software are presented. The optimization employs a density-based method and the various possible interpretations of the optimized shapes are of particular interest.

Eventually, let us mention the works [START_REF] Besserud | Structural emergence: architectural and structural design collaboration at som[END_REF] and [START_REF] Stromberg | Topology optimization for braced frames: combining continuum and beam/column elements[END_REF], where examples using topology optimization in the initial conceptual phase are presented and the importance of these techniques in obtaining knowledge about the structural behaviour of optimized systems is highlighted.

Limitations.

Let us now discuss the potential limitations of the use of shape and topology optimization in architecture, some of which have already arised from the above bibliography. These can be roughly classified into three categories with limitations related to:

(1) The practical realization of the obtained results.

(2) The simplified mechanical framework considered.

(3) The shape and topology optimization method.

The first class of limitations mainly comes up when the result of the optimization process is intended to be produced and used as such, or at least with minimal modifications. Then, the utility of shape and topology optimization methods is mostly restricted to art works, since it requires the use of costly advanced fabrication techniques. Although significant progress has been made in this direction, via techniques such as CNC-milling, HWC, or even additive manufacturing methods, they are still far from being practical for common applications. On the contrary, when shape and topology optimization is applied for conceptual design, i.e. when its results serve as a guide for novel shapes or efficient structural systems, this limitation is of minor concern.

Concerning the second class of limitations, the use of a simplified mechanical description of the real mechanical problem is legitimated by several arguments. Firstly, when shape and topology optimization is intended for conceptual design, the obtained results are likely to undergo a serious post-treatment. Thus, considering more sophisticated mechanical descriptions, as in [START_REF] Amir | Topology optimization for conceptual design of reinforced concrete structures[END_REF][START_REF] Amir | Reinforcement layout design for concrete structures based on continuum damage and truss topology optimization[END_REF] (for different applications however), would only complicate the optimization problem unreasonably. Furthermore, when aesthetics is the principal interest, it is quite reasonable to slip in the intuitive framework where simple gravity loads occur and to consider the compliance as the objective function. Admittedly when the performance of the structural system is at the heart of the problem, a more complete model, including e.g. a possible material anisotropy, i.e. different mechanical behaviors under tension and compression, or even multiple phases could be beneficial [START_REF] Aage | Advanced topology optimization methods for conceptual architectural design[END_REF][START_REF] Dombernowsky | Three-dimensional topology optimisation in architectural and structural design of concrete structures[END_REF].

Last but not least, let us comment on the shape and topology optimization methodologies used so far in the field of architecture. They either feature a density-based description of structures, like the SIMP method [START_REF] Bendsoe | Topology optimization: theory, methods and applications[END_REF], or rely on an evolutionary approach, like the BESO method [START_REF] Querin | Evolutionary structural optimisation (eso) using a bidirectional algorithm[END_REF], in which elements of the Finite Element mesh are removed or added according to a (usually heuristic) mechanical criterion.

In our opinion, a significant flaw of density-based methods as regards architectural applications is the loss of almost any description of the shape, and handling geometric criteria is then far from obvious. Recently, some works have been presented in this direction, such as the consideration of flexible void areas in [START_REF] Clausen | Topology optimization with flexible void area[END_REF] and the application of projection filters in [START_REF]Topology optimization of coated structures and material interface problems[END_REF] to simulate an interface. However, such methods are expected to be less specific and thus less efficient than traditional geometric methods for shape optimization, in which an explicit description of structures is brought into play. Another significant drawback of density-based methods is the well-known difficulty of interpreting the final result: a classical shape is revealed only in case the resulting density field from the optimization process takes values close to 0 or 1. Although there is usually an inherent mechanism (as in the SIMP formulation) that favors this tendency, numerical results very often contain a lot of intermediate densities, as shown in [START_REF] Feringa | An integral approach to structural optimization and fabrication[END_REF][START_REF] Vazquez | Design, analysis and fabrication of expressive, efficient shell structures: a prototype exploring synergy between architecture, engineering and manufacture[END_REF]. Hence, an interpretation step is necessary, and beyond the difficulty of such a step, it can cause in a significant loss of optimality for the final shape.

Finally, in BESO methods, the heuristic criterion used for the shape change can be proved efficient only for very simple criteria. In case of complicated functionals, where the gradient does not have a constant sign (as in the case of the compliance or the volume), such methods may not provide satisfying results.

In a nutshell, we believe that there is room for improving the existing structural optimization methods in the perspective of applications in architecture, notably by bringing into play a precise geometric description of shapes. The framework and geometric constraints introduced in this work are proposals in this direction.

An architecture-oriented approach

In this section, we detail the basic theoretical and numerical ingredients of the proposed shape and topology optimization method for architecture: in Section 3.1, we present the mathematical model of linearized elasticity for the description of the physical behavior of our shapes; then, in Section 3.2, we summarize Hadamard's method for describing variations of shape, and the inferred notion of shape derivative. Eventually, in Section 3.3, we recall the main features of the level set method, used in the numerical representation of shapes and their motions.

Let us emphasize once more that the selection of such 'geometric' (theoretical and numerical) representations of shapes -as opposed to the so-called density-based methods for topology optimization -is motivated by our intention to incorporate their geometric features (normal vector, etc...) into the optimization problem: see the constraint functionals proposed in Section 5.

Mechanical framework: compliance minimization in linear elasticity.

In this work, we limit ourselves with the model problem of compliance minimization in the stationary context of linearized elasticity, which already reflects a great deal of the complex physics involved.

From the mathematical viewpoint, the considered structures (or shapes) in the present work are bounded domains Ω ⊂ R d enjoying at least Lipschitz regularity (in our applications, d = 2, 3). These shapes are filled with an isotropic linear elastic material with Hooke's tensor A given by: ∀e ∈ S d (R), Ae = 2µe + λtr(e)I, where S d (R) denotes the set of symmetric d × d matrices, and λ, µ are the Lamé coefficients of the material, satisfying µ > 0, and λ + 2µ/d > 0.

The considered shapes Ω are clamped on a region Γ D of their boundary; they are subjected to traction loads g ∈ L 2 (Γ N ) d , exerted on a disjoint subset Γ N ⊂ ∂Ω, and body forces f ∈ L 2 (Ω) d (accounting for instance for gravity). In this situation, the displacement u Ω of Ω is the unique solution in the space

H 1 Γ D (Ω) d := u ∈ H 1 (Ω) d
, u = 0 on Γ D to the linear elasticity system:

(1)

       -div(Ae(u)) = f in Ω, u = 0 on Γ D , Ae(u)n = g on Γ N , Ae(u)n = 0 on Γ,
where e(u) = 1 2 (∇u + ∇u T ) is the strain tensor associated to u. In this situation, we are interested in minimizing a criterion J(Ω) of the domain. One such criterion is the compliance C(Ω) of shapes, which is equivalently defined as the elastic energy stored in Ω, or as the work of external loads:

(2)

C(Ω) = Ω Ae(u Ω ) : e(u Ω ) dx = Ω f • u Ω dx + Γ N g • u Ω ds.
The optimization is carried out on the set U ad of admissible shapes defined by

(3)

U ad = Ω ⊂ R d is a smooth bounded domain, Γ D ∪ Γ N ⊂ ∂Ω .
In most applications, some constraints, modelled by a functional P (Ω), should be imposed on shapes; there are classically two ways for incorporating them in the optimization problem:

• One solution consists in penalizing the performance criterion, in which case the optimization problem rewrites: min

Ω L(Ω), where L(Ω) := J(Ω) + mP (Ω),
where m > 0 is a given weighting coefficient. • They could alternatively be enforced as hard constraints, so that the optimization problem reads:

min Ω s.t. P (Ω)≤α J(Ω),
where α is a maximum threshold for the values of P (Ω). The perhaps most obvious example of a constraint is associated to the volume function P (Ω) = Vol(Ω) := Ω dx; the main purpose of this paper is to devise other penalty functionals associated to geometric constraints which could find applications in architectural design.

Remark 3.1. In the above framework, we have considered that boundary conditions (that is, the fixations and loaded regions) are predetermined, as is usually the case in applications of mechanical engineering (see the definition (3) of admissible shapes). However, observing the change in the optimal shape when various sets of fixations are used is of great interest for architects, since there is a certain latitude as for the choice of this set. Let us explain this point with an example extracted from [START_REF] Bando | Optimisation d'une structure et application architecturale[END_REF] 1 . In Figure 2, we see the optimized shapes for three optimization problems, where a uniform force is distributed on the upper side of the design domain, and which differ only by the nature and locations of the imposed fixations. For the upper case, the displacement is fixed at the two right corners of the lower side, while the corners on the left side are free to move in the x-direction. The optimized structural system has a clear physical interpretation: in the middle of the structure, straight bars in compression transfer the loads to an inverted arc in tension. Then, two inclined bars in compression, as well as some horizontal bars in tension, transfer the loads to the fixations. Setting the same type of fixations all along the lower left and right-hand sides, the material in the optimized shape is mainly concentrated in the centre of the structure. Finally, imposing null displacement on both sides, the tractor in the inferior part of the structure, accounting for the horizontal displacements is removed and the solution resembles an arc in compression, with some vertical bars attached.

Mathematical optimization framework: perturbations of shapes and differentiation with respect to the domain.

Most optimization algorithms rely on the knowledge of derivatives of the involved objective function. In the present context, we aim at optimizing functionals of the domain. In this situation, several notions of differentiation are available and we rely on Hadamard's boundary variation method [START_REF] Allaire | Conception optimale de structures[END_REF][START_REF] Henrot | Variation et optimisation de formes: une analyse géométrique[END_REF][START_REF] Murat | Etude de problèmes d'optimal design[END_REF] (see [START_REF] Novotny | Topological Derivatives in Shape Optimization[END_REF] about the notion of topological derivative): variations of a shape Ω are considered under the form (see Figure 3 for an illustration):

Ω θ := (Id + θ)(Ω), θ ∈ W 1,∞ (R d , R d ), ||θ|| W 1,∞ (R d ,R d ) < 1. Definition 3.1. A function F (Ω) is shape differentiable at a particular shape Ω if the underlying map- ping θ → F (Ω θ ), from a neighborhood of 0 in W 1,∞ (R d , R d ) into R, is Fréchet differentiable at 0. The corresponding derivative θ → F (Ω)(θ)
is the shape derivative of F at Ω, and the following expansion holds:

(4) F (Ω θ ) = F (Ω) + F (Ω)(θ) + o(θ) where lim θ→0 |o(θ)| ||θ|| W 1,∞ (R d ,R d ) = 0.
Notice that, in general, so that perturbations Ω θ of admissible shapes Ω ∈ U ad stay admissible, the vector fields θ have to be restricted to a set Θ ad of admissible perturbations defined by:

Θ ad = θ ∈ W 1,∞ (R d , R d ) smooth, θ = 0 on Γ D ∪ Γ N .
1 The unpublished work carried out in this reference was obtained with the collaboration of some of the authors, using the exact same algorithms as those detailed in this article Variations Ω θ of a shape Ω using Hadamard's method.

Recall that, according to the so-called Structure theorem (see e.g. [START_REF] Henrot | Variation et optimisation de formes: une analyse géométrique[END_REF], §5.9, or [START_REF] Delfour | Shapes and geometries[END_REF], Chap. 9, Th. 3.6), the shape derivative F (Ω)(θ) of a 'large' class of functions F (Ω) only depends on the values of the normal component θ • n of the considered perturbation on ∂Ω. Actually, as we shall see, the shape derivatives of the considered integral functions of the domain will enjoy the more precise structure ( 5)

F (Ω)(θ) = ∂Ω v Ω θ • n ds,
where v Ω is a scalar field depending on F . For instance, it is well-known (see e.g. [START_REF]Structural optimization using sensitivity analysis and a level-set method[END_REF]) that when F (Ω) is the compliance (2), one has:

(6) v Ω = -Ae(u Ω ) : e(u Ω ).
In the general case, from the expression (5), a descent direction for F (Ω) is revealed under the choice

(7) θ = -tv Ω n,
for a small, positive pseudo time-step t > 0.

3.3. Numerical framework: encoding shapes via the level set method.

Among the existing 'geometric' techniques for the numerical representation of shapes in the perspective of their optimization, we have selected the 'classical' level set method on a fixed mesh, due to its relative simplicity and robustness. To mention a few, alternative 'geometric' techniques for shape and topology optimization represent designs e.g. by means of a computational mesh, as in [START_REF] Allaire | Shape optimization with a level set based mesh evolution method[END_REF][START_REF] Christiansen | Topology optimization using an explicit interface representation[END_REF][START_REF] Persson | Mesh generation for implicit geometries[END_REF], or by a set of simple primitives [START_REF] Guo | Doing topology optimization explicitly and geometricallya new moving morphable components based framework[END_REF][START_REF] Zhang | A new topology optimization approach based on moving morphable components (mmc) and the ersatz material model[END_REF].

The level set method was pioneered in [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[END_REF], then introduced in the context of shape optimization in [START_REF] Allaire | A level-set method for shape optimization[END_REF][START_REF]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]; see the monographs [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] for extensive discussions around the features of this method. The main idea is to describe a shape Ω ⊂ R d via a scalar 'level set' function φ : R d → R with the properties (see Figure 4): During the optimization process the shape Ω ≡ Ω(t) evolves in (pseudo-)time with a motion driven by a normal velocity field V (t, x), which stems from the shape derivative of the objective criterion under consideration, in the line of Section 3.2 (see [START_REF] Allaire | A level-set method for shape optimization[END_REF]. The motion translates in terms of an associated level set function φ(t, x) as the well-known Hamilton-Jacobi equation: [START_REF] Amir | Topology optimization for conceptual design of reinforced concrete structures[END_REF] ∂φ ∂t + V |∇φ| = 0, which can be solved, e.g., on a Cartesian grid of a large computational box D, with an explicit second order upwind scheme [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 on a simplicial mesh of D (i.e. composed of triangles in 2d and tetrahedra in 3d) with a semi-Lagrangian scheme [START_REF] Strain | Semi-lagrangian methods for level set equations[END_REF].

(8)    φ(x) < 0 if x ∈ Ω, φ(x) = 0 if x ∈ ∂Ω, φ(x) > 0 if x ∈ c Ω.
Remark 3.2. In the shape and topology optimization context of interest in this article, the practical calculation of the velocity field V (t, x) cannot be performed exactly: its expression (6-7) involves the elastic displacement u Ω , the solution of (1). This system is posed on the shape Ω, and is typically solved with the Finite Element method, which requires a discretization of Ω by means of a mesh. In our context, Ω is not discretized exactly, and is only known via an associated level set function φ, defined at the vertices of a (fixed) mesh of the computational domain D. To get past this difficulty, we rely on an approximation which is quite common in the treatment of shape and topology optimization problems on a fixed mesh (see e.g. [START_REF]Structural optimization using sensitivity analysis and a level-set method[END_REF][START_REF] Bendsoe | Topology optimization: theory, methods and applications[END_REF]), namely the Ersatz material approximation: the void region D \ Ω is filled with a soft material, i.e. ( 1) may be accurately approximated by the same linear elasticity system posed on the whole domain D, where the Hooke's tensor A is replaced by that A ε defined by:

A ε (x) = A if x ∈ Ω, εA if x ∈ D \ Ω,
for a small value ε 1. This last system may be solved by a standard Finite Element method.

Preliminaries about the signed distance function

In this section, we collect some material about the notions of signed distance function and unit normal vector to a domain Ω ⊂ R d ; we also discuss their extensions to the context where the ambient space is equipped with a Riemannian structure, induced by an anisotropic metric tensor. This material plays a fundamental role in the expression of many geometric features of shapes (such as minimum and maximum length scales, see for instance [START_REF] Allaire | Mutli-phase structural optimization via a level-set method[END_REF][START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF][START_REF] Guo | Explicit feature control in structural topology optimization via level set method[END_REF][START_REF] Xia | Constraints of distance from boundary to skeleton: For the control of length scale in level set based structural topology optimization[END_REF]), in particular in the definitions of our geometric constraints in Section 5,

The present section purposely stays at an elementary level, and all the underlying technical details (including mathematical proofs) are postponed to Appendix A.

The signed distance function.

Let us start with the definition of the main notion of this section.

∀x ∈ R d , d Ω (x) =    -d(x, ∂Ω) if x ∈ Ω, 0 if x ∈ ∂Ω, d(x, ∂Ω) if x ∈ c Ω, where d(x, ∂Ω) = inf y∈∂Ω |x -y| is the usual Euclidean distance function to ∂Ω.
We now state several definitions and properties in connection with the signed distance function. Even though they will not be used as such in the sequel, they are very similar (yet much simpler to understand) to their anisotropic counterparts addressed in Section 4.3, and they are key in the proofs of the latter. Notice that we do not attempt to provide statements under minimal assumptions. Definition 4.2. Let Ω ⊂ R d be a bounded, Lipschitz domain.

• For x ∈ R d , the set of projections Π ∂Ω (x) of x onto ∂Ω is defined by: Π ∂Ω (x) = {y ∈ ∂Ω, |x -y|= d(x, ∂Ω)}; it is a non empty, closed subset of ∂Ω. When Π ∂Ω (x) is a singleton, its unique element is denoted as p ∂Ω (x) and is called the projection of x onto ∂Ω. • The skeleton Σ of Ω is the set of points x ∈ R d where d 2 Ω is not differentiable; note in passing that Σ ∩ ∂Ω = ∅, and that Σ has null Lebesgue measure as a consequence of Rademacher's theorem.

• The ray of ∂Ω emerging from a point y ∈ ∂Ω is the set:

ray ∂Ω (y) = p -1 ∂Ω (y) = x ∈ R d , Π ∂Ω (x) = {y} .
The following proposition is a pot-pourri of several results in [START_REF] Cannarsa | Representation of equilibrium solutions to the table problem of growing sandpiles[END_REF][START_REF] Delfour | Shapes and geometries[END_REF][START_REF] Henrot | Variation et optimisation de formes: une analyse géométrique[END_REF]] and [START_REF] Li | The distance function to the boundary, finsler geometry[END_REF], which are illustrated on Figure 5.

Proposition 4.1. Let Ω ⊂ R d be a bounded, Lipschitz domain. (i) A point x /
∈ ∂Ω has a unique projection over ∂Ω if and only if x / ∈ Σ. Then, d Ω is differentiable at x and its gradient reads: [START_REF] Amir | Reinforcement layout design for concrete structures based on continuum damage and truss topology optimization[END_REF] ∇d

Ω (x) = x -p ∂Ω (x) d Ω (x) .
In particular, |∇d Ω (x)|= 1 wherever it makes sense.

(ii) If Ω is of class C 1 , then for x / ∈ ∂Ω, x / ∈ Σ, (11)
∇d Ω (x) = n(p ∂Ω (x)), and

x = p ∂Ω (x) + d Ω (x)n(p ∂Ω (x)). (iii) If Ω is of class C 2
, the closure Σ also has null Lebesgue measure. Moreover, there exists a tubular neighborhood

U of ∂Ω such that Σ ∩ U = ∅. (iv) Still in the case where Ω is of class C 2 , d Ω is of class C 2 on R d \ Σ. Actually, d Ω is twice differentiable on R d \ Σ
, and ( 11) actually holds on R d \ Σ (in particular, it holds in the vicinity of ∂Ω).

⌃ • x • • ⇧ @⌦ (x) ⌦ • y • p @⌦ (y)
n(p @⌦ (y))

• z ray @⌦ (z) The foregoing properties are completed with the following proposition which describes the second-order derivatives of d Ω in terms of the second-order character of Ω (see again [START_REF] Cannarsa | Representation of equilibrium solutions to the table problem of growing sandpiles[END_REF] and [START_REF] Allaire | Mutli-phase structural optimization via a level-set method[END_REF]). Proposition 4.2. Let Ω ⊂ R d be a bounded domain of class C 2 ; for y ∈ ∂Ω, we denote by T y ∂Ω be the tangent plane to ∂Ω at y; II y : T y ∂Ω × T y ∂Ω → R stands for the second fundamental form of ∂Ω at y, and κ 1 (y), ..., κ d-1 (y) (resp. (e 1 (y), ..., e d-1 (y))) are the associated principal curvatures (res.p principal directions), obtained as the eigenvalues (resp. eigenvectors) of II y .

(i) For any point x ∈ R d , and p ∈ Π ∂Ω (x), one has:

(12) -κ i (p)d Ω (x) ≤ 1, i = 1, ..., d -1.
(ii) The closure Σ of the skeleton Σ reads as the disjoint union: Σ = Σ ∪ Γ, where Γ is the set of points x ∈ R d where one of the inequalities in ( 12) is an equality.

(iii) The projection mapping p ∂Ω is differentiable at any x ∈ R d \ Σ and its derivative reads:

(13) ∇p ∂Ω (x) =        1 -dΩ(x)κ1(p) 1+dΩ(x)κ1(p) 0 ... 0 0 . . . . . . . . . . . . . . . 1 -dΩ(x)κ d-1 (p) 1+dΩ(x)κ d-1 (p) 0 0 ... 0 0        , p ≡ p ∂Ω (x),
where the above matrix is written in the orthonormal basis (e 1 (p), ..., e d-1 (p), n(p)) of R d .

4.2.

The normal vector to a domain.

In the following, we equivalently denote by n Ω or n the unit normal vector to ∂Ω, pointing outward Ω. When Ω is of class C 2 , Proposition 4.1 indicates that n Ω : ∂Ω → S d-1 may be extended to an open neighborhood U of ∂Ω (and actually to R d \ Σ) via the formula:

n Ω (x) = ∇d Ω (x) = n(p ∂Ω (x)), x ∈ U.
This choice of an extension for the normal vector enjoys the interesting symmetry property: [START_REF] Bendsoe | Topology optimization: theory, methods and applications[END_REF] ∇n

Ω = ∇ 2 d Ω = ∇n T Ω .
In the following, we will need to evaluate the sensitivity of the normal vector n Ω with respect to the shape Ω; the main result in this direction is the following (see for instance [START_REF] Murat | Etude de problèmes d'optimal design[END_REF]):

Proposition 4.3. The normal vector Ω → n Ω has a material derivative in the sense that the mapping to Ω reads:

θ → n Ω θ • (Id + θ), from C 1,∞ (R d , R d ) into C(∂Ω) is Fréchet differentiable at 0. Its derivative ṅΩ (θ) reads: ṅΩ (θ) = ∇θ T n Ω , n Ω n Ω -∇θ T n Ω .
d M Ω (x) =    -d M (x, ∂Ω) if x ∈ Ω, 0 if x ∈ ∂Ω, d M (x, ∂Ω) if x ∈ c Ω, where d M (x, ∂Ω) := inf y∈∂Ω |x -y| M .
Before going any further, let us consider an intuitive example in the case of the two-dimensional shape depicted in Figure 6 (left), and of the metric tensor, written in the canonical basis:

(15) M = 0.1 0 0 1 .
Here, imposing a metric with a much lower coefficient in the horizontal direction implies that, from a point x ∈ R 2 , travelling until hitting ∂Ω 'costs' almost nothing (that is, does not increase much the travelled distance), whereas it does cost very much when the travel has some vertical component. Hence, in the vertical section of the cross, where any point can be connected to a point of the boundary by an horizontal travel, the anisotropic signed distance function is very low, whereas it is larger on the horizontal bar. In this section, we discuss the generalization of the properties of Section 4.1 in the anisotropic case. This material is new to the best of our knowledge, and we provide the proofs in Appendix A.

First, note that Definition 4.2 extends verbatim to the anisotropic context. Hence, we shall speak about:

• The anisotropic projection set Π M ∂Ω (x) (or point p M ∂Ω (x) when it is a singleton) of a point x ∈ R d , • The anisotropic skeleton Σ M of Ω, • The anistropic ray ray M ∂Ω (x) of a point x ∈ R d . Proposition 4.4.
Let Ω ⊂ R d be a bounded, Lipschitz domain, and M be a symmetric, positive definite matrix.

(i) For x / ∈ ∂Ω, the set Π M ∂Ω (x) is a singleton if and only if x / ∈ Σ M . Then, d M Ω is differentiable at x and its gradient reads:

(16) ∇d M Ω (x) = M (x -p M ∂Ω (x)) d M Ω (x) . In particular, M -1 ∇d M Ω , ∇d M Ω = 1 holds wherever it makes sense. (ii) If Ω is of class C 1 , then for x / ∈ ∂Ω, x / ∈ Σ M , ( 17
) ∇d M Ω (x) = n M (p), and x = p + d M Ω (x)n M (p), p ≡ p ∂Ω (x),
where the distorted normal n M (y) reads n M (y) :

= 1 M -1 n(y),n(y) M -1 n(y) (note that |n M (y)| M = 1.
We now assume Ω to be at least of class C 2 .

(iii) The set Σ M has zero Lebesgue measure. Moreover, there exists an open neighborhood U of ∂Ω such that U ∩ Σ M = ∅, and (17) actually holds on R d \ Σ M . (iv) For any point x ∈ R d , the matrix

I + 1 |M -1 2 n Ω (p)| d M Ω (x)M -1 2 II p M -1 2 , p ≡ p M ∂Ω (x),
is positive semi-definite. Moreover, if x / ∈ Σ M , this matrix is positive definite.

(v) The projection p M ∂Ω is differentiable on R d \ Σ M and its derivative reads:

(18) ∇p M ∂Ω (x) =          I + 1 |M -1 2 nΩ(p)| d M Ω (x)M -1 2 II p M -1 2 -1 0 0 0          , p ≡ p M ∂Ω (x),
the latter matrix being written in an orthonormal basis of R d of the form

M 1 2 τ 1 (p) |M 1 2 τ 1 (p)| , • • • , M 1 2 τ d-1 (p) |M 1 2 τ d-1 (p)| , M -1 2 n Ω |M -1 2 n Ω (p)| , for any set of tangent vectors τ i (p) ∈ T p ∂Ω, i = 1, ..., d -1 satisfying the identity M τ i (p), τ j (p) = δ ij .
4.3.2. Shape differentiation of the anisotropic, signed distance function.

We now proceed to the shape differentiation properties of the anisotropic signed distance function. Our first result is about the pointwise shape differentiability of d M Ω : Proposition 4.5. Let Ω ⊂ R d be a bounded and Lipschitz domain, and let

x ∈ R d \ ∂Ω. The mapping (0, ∞) t → d M Ω tθ (x) ∈ R is differentiable at t = 0 +
, and its derivative reads:

d dt (d M Ω tθ (x))(0 + ) = inf y∈Π M ∂Ω (x) y -x, θ(y) M d M Ω (x)
.

In particular, if Ω is of class C 1 , and x / ∈ Σ M , the above formula rewrites:

d dt (d M Ω tθ (x))(0 + ) = -n M (p), θ(p) M
, where p ≡ p M ∂Ω (x).

We now come to the differentiation of integral functionals involving d M Ω . Proposition 4.6. Let Ω ⊂ R d be a bounded domain of class C 1 , contained in a larger domain D, and let j : R → R be a function of class C 1 . Let J(Ω) be the shape functional defined by:

J(Ω) = D j(d M Ω (x)) dx.
Then J(Ω) is shape differentiable at Ω in the sense of Gâteaux, and its shape derivative reads:

∀θ ∈ W 1,∞ (R d , R d ), J (Ω)(θ) = - D j (d M Ω (x)) n M (p M ∂Ω (x)), θ(p M ∂Ω (x)) M dx

Three geometric criteria for architectural applications

In the shape and topology optimization literature, several geometric constraints have been proposed, which may find a great interest in architectural design (even though it seems that such applications have not yet been investigated):

• Constraints on the maximum and minimum feature size of shapes -hereafter referred to as 'maximum thickness' and 'minimum thickness' -have notably been dealt with in [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF][START_REF] Sigmund | Morphology-based black and white filters for topology optimization[END_REF]. Let us mention that whereas restricting the maximum thickness of shapes is of particular interest, mainly due to fabrication limitations, controlling the minimum feature size is rarely necessary for applications in architecture since the scale of civil structures is quite large. Penalizing small features could be used mainly to simplify the optimized design. • Another way to control the complexity of the shape is to add a penalization over the total perimeter (see for instance [START_REF]Structural optimization using sensitivity analysis and a level-set method[END_REF][START_REF] Yamada | A topology optimization method based on the level set method incorporating a fictitious interface energy[END_REF] on this topic), but the amount of personalized information offered by this ingredient is limited to the extent of the penalization.

• The possibility of imposing non-optimizable areas [START_REF] Clausen | Topology optimization with flexible void area[END_REF] leaves larger room for a personalized shape control. For instance, imposing fictitious holes in the working domain (that is, areas where the optimized structure cannot infringe) is quite simple and allows to explore different designs, using different strategies to bear the imposed loads (see the example of Figure 7, excerpted from [START_REF] Bando | Optimisation d'une structure et application architecturale[END_REF]). In the present work, we propose three novel geometric criteria that are very much suited for architectural applications:

• In Section 5.1, we propose a shape functional P m (Ω) which constrains shapes to stay 'close' to a user-specified design. • In Section 5.2, we introduce a shape functional P p (Ω) which constrains Ω to conform to a userspecified pattern (or general trend) rather than strictly a design. • In Section 5.3, a function P M a (Ω) is devised with the goal to reduce the anisotropic thickness of shapes with respect to a user-defined metric tensor M . This comes in handy when trying to improve the visibility of shapes, or in an attempt to prevent the emergence of too elongated members in some direction.

These criteria are intended to be inserted either as penalization in the minimized function, or as constraints in the optimization problem (see the discussion in the end of Section 3.1, and the examples of Section 6).

5.1.

Requiring that shapes stay close to a specific design.

The first functional of interest for our purpose aims at keeping the proposed design Ω close to a reference shape Ω T . In fact, if one tries to optimize the shape and topology of a structure starting from a reference design, nothing guarantees that the final optimized shape with ressemble the initial proposal. Thus, one shall try to formulate a functional penalizing some sort of deviation from the reference domain.

Suppose we want Ω to resemble some specific design Ω T , supplied via its signed distance function d Ω T . Morphing the design domain Ω to the target domain Ω T can be achieved by minimizing the functional:

(19) P m (Ω) = Ω d Ω T dx.
The intuition behind [START_REF] Cannarsa | Representation of equilibrium solutions to the table problem of growing sandpiles[END_REF] is that every region of Ω comprised inside Ω T (resp. inside c Ω T ) has a negative (resp. positive) contribution in P m (Ω). Hence, the minimization of P m (Ω) drives the regions of Ω lying inside Ω T to expand, and those lying outside Ω T to retract. In particular,

• If Ω T is connected, the unique local minimizer of P m (Ω) is Ω T itself, • In the general case where Ω T has several connected components, say O 1 , ..., O N , the local minimizers of P m (Ω) are exactly the sets i∈I O i , where I ⊂ {1, ..., N }, and the unique global minimizer of P m (Ω) is again Ω T .

The shape derivative of ( 19) is easily calculated as (see e.g. [START_REF] Henrot | Variation et optimisation de formes: une analyse géométrique[END_REF]):

Proposition 5.1. The functional P m (Ω) defined by ( 19) is shape differentiable at any shape Ω ∈ U ad , and its shape derivative reads:

(20) ∀θ ∈ Θ ad , P m (Ω)(θ) = ∂Ω d Ω T θ • n ds.
Remark 5.1.

• Functionals of the form [START_REF] Cannarsa | Representation of equilibrium solutions to the table problem of growing sandpiles[END_REF], based on the signed distance function have already been used in [START_REF] De Buhan | An optimization method for elastic shape matching[END_REF] (or prior in [START_REF] Bajcsy | Multiresolution elastic matching[END_REF]) in the contexts of shape morphing and shape registration. • Different functionals P (Ω) from ( 19) could be devised for the same shape matching purpose; for instance, one may imagine using a least-square functional of the form:

P (Ω) = D (d Ω -d Ω T ) 2 dx,
where D is a fixed computational domain. However, in our opinion, [START_REF] Cannarsa | Representation of equilibrium solutions to the table problem of growing sandpiles[END_REF] has the nice feature to enjoy a very simple shape derivative; besides, least-square functionals often cause conditionning issues in optimization problems (see for instance [START_REF] Dambrine | On the ersatz material approximation in level-set methods[END_REF]), and we prefer to avoid the latter insofar as possible.

5.2.

Imposing a resemblance with a user-defined pattern.

We are now interested in imposing that the considered shapes fit a user-defined pattern, defined throughout the ambient space R d , rather than sticking to a shape, strictly speaking. As an example, one may look for an optimized structure which enjoys a wavy boundary (see the numerical illustration in Section 6.2).

In our applications, the pattern is supplied as the 0 level set of a given smooth function g : R d → R. We assume that the gradient ∇g does not vanish, so that the unit normal vector n g to the level sets of g reads:

n g = ∇g |∇g| .
As a result of the optimization process, we are interested in that the normal vector n Ω of Ω be as close as possible to that of the pattern; this may be of interest either for aesthetic purposes, or to drive the light reflection properties of Ω, which are directly tied to the orientation of its boundary. The function P p (Ω) of interest is then:

(21) P p (Ω) = ∂Ω |n Ω -n g | 2 ds.
The shape differentiability of P p (Ω) is now studied in a slightly larger context -see Section A.2 for a proof of the forthcoming results. Proposition 5.2. Let P (Ω) be the shape functional defined by:

(22) P (Ω) = ∂Ω ϕ(n Ω ) ds,
where ϕ : R d → R is a function of class C 1 . Then P (Ω) is shape differentiable in the sense of Fréchet at any shape Ω ∈ U ad , and its derivative reads:

∀θ ∈ Θ ad , P (Ω)(θ) = Γ κ ϕ(n) θ • n ds - Γ ∇ϕ(n) • ∇ ∂Ω (θ • n) ds,
where ∇ ∂Ω f := ∇f -∇f, n Ω n Ω is the tangential gradient of a smooth enough function f : ∂Ω → R.

In particular, we infer the following corollary:

Corollary 5.1. The functional P p (Ω) defined by ( 21) is shape differentiable at any shape Ω ∈ U ad , and its shape derivative reads:

∀θ ∈ Θ ad , P p (Ω)(θ) = 2 ∂Ω (κ -div ∂Ω (n g ))θ • n ds,
where div ∂Ω V := div(V )-∇V n, n is the tangential divergenceof a smooth enough vector field V : ∂Ω → R d .

The result of Corollary 5.1 is quite intuitive: in order to minimize P p (Ω), the mean curvature of ∂Ω has to be adjusted to that of the pattern.

Penalization of the directional thickness of shapes.

In this section, we discuss a third class of constraints which have interesting applications in architectural design, namely constraints on the anisotropic maximum thickness of shapes.

Mathematical formulation of anisotropic maximum thickness constraints.

Let us start with a rigorous account of the notion of (isotropic or anisotropic) maximum thickness of shapes.

Definition 5.1. Let Ω ⊂ R d be a bounded domain, M be a symmetric, positive definite d × d matrix. Then,

• Ω has maximum thickness smaller than δ > 0 if:

∀x ∈ Ω, d Ω (x) ≥ -δ/2.
• Ω has anisotropic maximum thickness smaller than δ > 0 if:

∀x ∈ Ω, d M Ω (x) ≥ -δ/2.
Imposing constraints on the maximum thickness of shapes has countless applications in mechanical engineering. It is, for instance, a means to prevent the appearance of too thick regions, which could slow down the cooling process during manufacturing; see [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF] and references therein. However, to our knowledge, constraints on the anisotropic maximum thickness of shapes have not yet been considered, and could help in testing (then constraining) the thickness of shapes in only one (or several) selected direction -for instance the length of bars in certain directions, or the area of surfaces with a prescribed orientation.

To appraise the difference between the isotropic and anisotropic versions of the notion of maximum thickness, consider again the example of the cross Ω, depicted in Figure 6, (left), when the anisotropic metric tensor M is defined by [START_REF] Besserud | Structural emergence: architectural and structural design collaboration at som[END_REF]. In this case, the anisotropic signed distance function d M Ω (x) measures approximately the closest horizontal distance from a point x to a region of the boundary ∂Ω. Hence, if one aims at decreasing (in terms of absolute value) this value at any point x ∈ Ω, the solution (reflected by the shape derivative of Proposition 4.5) consists in reducing the 'horizontal thickness' of Ω. Hence, imposing a constraint on the anisotropic maximum thickness of Ω with respect to M is a way to constrain the length of its horizontal bars. Notice that relying on the notion of isotropic maximum thickness in such a situation is useless to this purpose, since the isotropic maximum thickness of Ω is already small in this situation.

Let us then construct a shape functional for imposing a bound δ > 0 on the anisotropic thickness of a shape Ω. Ideally, the desired constraint reads: ∀x ∈ Ω, d M Ω (x) ≥ -δ/2, which can be equivalently formulated using an integral penalty function: [START_REF] Clausen | Topology optimization with flexible void area[END_REF] P (Ω) = 0, where

P (Ω) = Ω d M Ω -δ/2 2 -dx,
where we denote by t -= max(0, -t) the negative part of a real number t ∈ R.

For reasons that are extensively detailed in [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF][START_REF] Michailidis | Manufacturing Constraints and Multi-Phase Shape and Topology Optimization via a Level-Set Method[END_REF], imposing maximum thickness constraints with this kind of functional turns out to be too strict, and produces undesirable numerical artifacts. The numerical experience reported in there suggests to impose this maximum anisotropic thickness as a constraint in the optimization algorithm (see the discussion in the end of Section 3.1) by using the following relaxed version of ( 23):

(24) P M a (Ω) ≤ δ 2 , where P M a (Ω) =     1 Ω h(d M Ω )dx Ω h(d M Ω ) (d M Ω ) 2 dx     1 2
.

In the above formula, h : R → R is a regularized characteristic function, defined as:

h(s) = 1 2 1 + tanh |s| -δ/2 α δ/2 ,
α > 0 being a parameter that controls the regularization of the constraint. The interpretation of h and α is the following: when α goes to 0, h converges to the function h * (s) = 1 if |s|> δ/2, 0 otherwise , a function that penalizes only regions where the maximum thickness constraint is violated. Of course, h * is way too steep to be used in practice, and we use a 'small', positive value for α; our experience suggests the heuristic choice α = 2 5 ∆x, where ∆x stands for the constant mesh size. Remark 5.2. This idea of using a Riemannian, anisotropic structure to measure lengths as a means of constraining the features of geometric objects in different ways depending on their orientation in space has already been used in other contexts, notably in meshing; see e.g. [START_REF] Frey | Mesh generation[END_REF][START_REF] Vallet | Anisotropic control of mesh generation based upon a voronoi type method[END_REF].

Shape derivative of the constraint functional P M a (Ω).

Let us now provide the main result concerning the shape differentiability of P M a (Ω); its proof is postponed to Section A.2.

Theorem 5.1. The constraint functional P M a (Ω) defined by ( 24) is shape differentiable at any admissible shape U ad , and its shape derivative in a direction θ ∈ Θ ad reads:

(25) P M a (Ω)(θ) = Γ D Ω θ • n ds, where D Ω (x) = ray M ∂Ω (x)∩Ω G(z) d-1 i=1 1 + d M Ω (z) |M -1 2 n(x)| λ i (x) d (z),
where d (z) stands for the one-dimensional Hausdorff measure on the rays ray M ∂Ω (x), the λ i (x), i = 1, ..., d-1, are the eigenvalues of the (d -1)

× (d -1) symmetric matrix M -1 2 II ∂Ω (x)M -1 2
, and the integrand factor G(z) is given by the expression:

(26) G(z) = - Ω h(d M Ω )dx h (d M Ω (z))d M Ω (z) 2 + h(d M Ω (z))2d M Ω (z) + Ω h(d M Ω )d M Ω 2 dx h (d M Ω (z)) 2P M a (Ω) Ω h(d M Ω )dx 2 .

Practical implementation of the maximum thickness constraint.

Unfortunately, the expression ( 25) is difficult to evaluate in numerical practice; taking inspiration from the 'approximate shape derivative' proposed in [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF] to handle the isotropic maximum thickness constraint, we rely on a quite heuristic simplification. This simplification, which proves very efficient in practice, builds upon the classical idea of changing inner products to infer a descent direction θ from the knowledge of the shape derivative of a functional (see e.g. [START_REF] Burger | A framework for the construction of level set methods for shape optimization and reconstruction[END_REF][START_REF] De Gournay | Velocity extension for the level-set method and multiple eigenvalues in shape optimization[END_REF]).

At first, an anisotropic diffusion equation is solved on the whole computational domain D; i.e. one calculates the solution Q ∈ H 1 (D) to the following variational problem: [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF] ∀v ∈ H 1 (D),

D α 2 r M -1 ∇Q • ∇v + Qv dx = Ω Gv dx,
where G is given in [START_REF] Dapogny | Optimisation de formes, méthode des lignes de niveaux sur maillages non structurés et évolution de maillages[END_REF], and the coefficient a r is set as

α 2 r =   2 max Ω |d Ω (x)| ∆x max i (M -1 ii )   (∆x) 2 .
A descent direction θ for the functional P M a (Ω) is then revealed as θ = -Qn. Let us comment shortly on the rationale behind this procedure (see [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF], §7.3 for further details). The function G carries all the needed information about the violation of the maximum thickness constraint [START_REF]Topology optimization of coated structures and material interface problems[END_REF]; in particular, it takes significant values in the regions of Ω where this constraint is violated, and almost zero values elsewhere. In this perspective, [START_REF] Dambrine | On the ersatz material approximation in level-set methods[END_REF] expresses a means to convey this information about constraint violation from internal regions of Ω to the boundary ∂Ω; in other terms, [START_REF] Dambrine | On the ersatz material approximation in level-set methods[END_REF] uses the volumetric information encoded in G to predict how the boundary ∂Ω should be deformed to guarantee a better fulfillment of the constraint [START_REF]Topology optimization of coated structures and material interface problems[END_REF].

Another means to achieve this information transfer is accounted for by [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF], which may be understood as a preconditionner for the shape optimization algorithm. The information contained in G undergoes a diffusion process from the regions where it takes large values to the boundary of Ω. This process is chosen to be anisotropic so as to favor diffusion in the direction where the maximum thickness constraint is prevailing. The coefficient α r is chosen so that the diffusion rate is large enough for the information to hit the boundary.

Numerical examples

In this last section, we present several numerical simulations to appraise the applicability of the shape optimization framework introduced in Section 3 and the geometric functionals of Section 5 in the field of architecture.

Our implementation differs slightly between the two-and three-dimensional cases. In both situations, a large computational box D is equipped with a fixed quadrilateral (Cartesian) mesh G, and each considered shape Ω is represented by means of a level set function φ (i.e. ( 8) holds) which is discretized at the vertices of G. The (constrained or unconstrained) optimization process is driven by an SLP-type algorithm similar to that presented in [START_REF] Dunning | Introducing the sequential linear programming level-set method for topology optimization[END_REF], implemented in Scilab [START_REF] Campbell | Modeling and simulation in Scilab/Scicos[END_REF].

When it comes to the Finite Element analyses involved in the resolution of the linearized elasticity system (1), our two-dimensional computations are performed using Scilab: the Cartesian mesh G of D also serves as a computational support for the execution of the Q 1 Finite Element method. In the three-dimensional case, we use a tetrahedral mesh T of D, together with P 1 Finite Elements, in the FreeFem++ environment [START_REF] Pironneau | [END_REF]; interpolation procedures are used to transfer data about the level set function φ from G to T , and conversely, to transfer data about the calculated elastic displacement u Ω from T to G.

In both situations, the anisotropic signed distance function d M Ω is calculated on G by solving an anisotropic counterpart to the so-called redistancing equation (see e.g. [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF] about the redistancing equation and [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF] and references therein about its anisotropic version), with a second-order numerical scheme [START_REF] Sethian | Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science[END_REF].

Let us emphasize that the fixed Cartesian mesh G (and tetrahedral mesh T in the 3d setting) could be easily replaced by a simplicial mesh of D, possibly changing in the course of the optimization process, so that it is consistently adapted to the shape; see for instance [START_REF] Frey | Mesh generation[END_REF] for related mesh adaptation techniques.

For all results, the considered elastic material is assumed to be isotropic, with normalized Young modulus E = 1 and Poisson's ratio ν = 0.33. The ersatz material (simulating the void) has Hooke's tensor εA, where ε = 1.10 -3 . 6.1. Imposing shapes to stick to a user-defined design.

Our first example deals with the geometric constraint introduced in Section 5.1, on the example of the optimal design of a two-dimensional bridge depicted on Figure 8. The working domain D has size 250 × 20; the considered structures are clamped on two small regions on their lower side, no body forces apply (i.e. f = 0 in (1)), and a uniform, unit vertical load g = (0, -1) is distributed on their upper side. Taking advantage of the symmetry of the test-case, the calculation is performed on only one half of D.

In a first experiment, we seek to minimize the compliance of the considered structures, under a volume constraint, namely: [START_REF] De Buhan | An optimization method for elastic shape matching[END_REF] min

Ω C(Ω) s.t. Vol(Ω) ≤ V T , where C(Ω) is the compliance (2).
In our applications, the maximum authorized volume is V T = 0.163|D|. Starting from the rather arbitrary initialization shown in Figure 9 (top), the optimized shape of Figure 9 (bottom) is obtained, whose compliance equals 9347.

Let us now assume that, out of aesthetic purposes, we would like that the optimal design resemble the target shape Ω T represented on Figure 10 (top), which uses the same amount V T of material. The compliance of Ω T is very high (C(Ω T ) = 48693), and we would like to improve its structural performance, that is to optimize the shape Ω in view of its compliance so that it has the same material volume |Ω|= V T , while staying 'close' to Ω T . Unfortunately, solving [START_REF] De Buhan | An optimization method for elastic shape matching[END_REF] with Ω T as an initial guess gives no guarantee that the resulting optimized shape will retain any resemblance with Ω T , as evidenced in Figure 10 (bottom). Although this final shape has a significantly improved performance over Ω T (its compliance equals 12817), it is not satisfactory from an aesthetics point of view. In order to incorporate the desired aesthetic information, we trade [START_REF] De Buhan | An optimization method for elastic shape matching[END_REF] for the new optimization problem (29), [START_REF] De Gournay | Velocity extension for the level-set method and multiple eigenvalues in shape optimization[END_REF] min

Ω L(Ω) s.t. Vol(Ω) ≤ V T , where L(Ω) := t C(Ω) C(Ω T ) + (1 -t) P m (Ω) P m (Ω T )
involving the penalization functional P m (Ω) defined by [START_REF] Cannarsa | Representation of equilibrium solutions to the table problem of growing sandpiles[END_REF], and a weighting coefficient t ∈ [0, 1]. Using the target shape Ω T as initialization for the problem [START_REF] De Gournay | Velocity extension for the level-set method and multiple eigenvalues in shape optimization[END_REF] and different values of t, the optimized shapes of Figure 11 are obtained. As expected, for low values of t, the topology of Ω T is retained by the optimized shape, while increasing t accounts for a smooth transition of optimal shapes from Ω T to the one of Figure 10 (bottom). The values of the compliance C(Ω) and of the penalization P m (Ω) for the resulting optimal designs are reported on Table 1, and the associated convergence histories are represented on Figure 12.

6.2.

Constraining shapes to fit a user-defined pattern. Table 1. Values of the compliance C(Ω) and of the penalization P m (Ω) for the optimal designs obtained in the 2d bridge example of Section 6.1.

Let us now apply the material of Section 5.2 to influence structures with a user-specified pattern. The pattern used in this example is given by the function g : D → R defined by (see Figure 13, left):

(30) ∀x = (x 1 , x 2 ) ∈ D, g(x) = x 1 - L 1 10 sin 2πx 2 3L 2 , L 1 = 1, L 2 = 2.
In a first example, we simply aim at evaluating the penalty functional P p (Ω) defined by ( 21): starting from the initialization of Figure 13 (middle), we solve the optimization problem [START_REF]Shape identification via metrics constructed from the oriented distance function[END_REF] min

Ω P p (Ω).
The algorithm converges to a local minimum, displayed on Figure 13 (right). As indicated by the shape derivative of Corollary 5.1, the mean curvature of the boundary at the optimal shape equals that of the motif at each point. This means that in regions where n • n g > 0, the algorithm tries to align the exterior normal vector with n g , while if n • n g < 0, a local minimum is achieved when n • n g = 0.

Let us now incorporate this aesthetic constraint in a structural optimization problem: the situation is that depicted on Figure 14. The shape Ω is enclosed in a computational box D of dimensions 1 × 2, and is clamped at some part of its lower side; a unit vertical load g = (0, -1) is distributed on its upper side and no body forces occur, i.e. f = 0.

We aim at minimizing the compliance under a volume constraint, and additionnally want to influence the shape with the sinusoidal pattern (30) of Figure 13 (left). To this end, we solve once more the optimization problem [START_REF] De Gournay | Velocity extension for the level-set method and multiple eigenvalues in shape optimization[END_REF], replacing P m (Ω) with P p (Ω). The results for different values of the weighting factor t are 30), (middle) initial shape, and (right) optimized shape for Problem [START_REF]Shape identification via metrics constructed from the oriented distance function[END_REF].

shown in Figure 15, where the influence of the pattern is obvious for decreasing values of t. The corresponding values for the compliance C(Ω) and the penalty P p (Ω) of the optimized shapes are reported in Values of the compliance C(Ω) and of the penalization P p (Ω) for the optimized designs obtained in the 2d pattern-fitting example of Section 6.2.

6.3. Illustrations of the constraint over the anisotropic maximum thickness of shapes. In this section, we now exemplify how the notion of anisotropic maximum thickness introduced in Section 5.3 may be used in architectural applications. In all the examples below, the considered optimization problem features an additional constraint corresponding to a maximum authorized value over the function P M a (Ω) defined in [START_REF]Topology optimization of coated structures and material interface problems[END_REF] (for a metric tensor M to be specified). Hence, the considered optimization problem reads: [START_REF] Dombernowsky | Three-dimensional topology optimisation in architectural and structural design of concrete structures[END_REF] min

Ω C(Ω) s.t. Vol(Ω) ≤ V T P M a (Ω) ≤ δ/2
, where δ > 0 is a threshold depending on the particular situation. 6.3.1. Simplification of the shape of a three-dimensional bridge.

We first show how structures may be simplified by enforcing a control on their anisotropic maximum thickness. The considered test case is that of the optimal bridge depicted on Figure 2, which is reproduced on Figure 16.

At first, a simple minimization of the compliance C(Ω) of the bridge under the volume constraint Vol(Ω) = V T is performed -i.e. ( 28) is solved -where the target volume is the fraction V T = 0.1|D| of that of the working domain D, and D acts as an initial shape. One can see that the resulting optimized shape Ω * contains regions with large thickness in the (yz)-plane, although its isotropic maximum thickness is low. The mechanical relevance of these parts is not clear and shall be further interpreted by the designer, to result in a more classical truss-like shape composed of bars in tension or compression.

To help this purpose, we now solve Problem (32) with a metric tensor M of the form

(33) M =   m 1 0 0 0 m 2 0 0 0 m 3   ,
where the coefficients m 1 , m 2 , m 3 are yet to be chosen. The threshold parameter δ is taken as δ = 0.3 min {m 1 , m 2 , m 3 }, and the optimal shapes associated to the values m 1 = 0.1, then m 1 = 0.01, and m 2 = m 3 = 1 are displayed on Figure [START_REF] Bradner | Parameters tell the design story: ideation and abstraction in design optimization[END_REF].

The associated values of the compliance of the optimized shapes are reported in Table 3, and the related convergence histories are those of Figure 17. 3. Values of the compliance of the optimized shapes of Figure 16.

One observes that the performance of all three results are quite similar. The thick part has disappeared and the material has been redistributed in the rest of the domain to retain rigidity. 6.3.2. Improvement of the visibility of a three-dimensional short cantilever.

Let us now show how a constraint on the anisotropic maximum thickness may improve the visibility of shapes.

Consider the 'short cantilever' test case of Figure 18 (left): in a working domain D of dimensions 1 × 2 × 1, the shape is clamped at its left side and two different loads g 1 , g 2 , defined by: g 1 = (-1, 0, 0), and g 2 = (0, -1, 0) are applied on the centre of its right-hand side. Again, no body forces apply: f = 0. The mechanical criterion whose minimization is sought is the total compliance S(Ω) and u i Ω is the elastic displacement of the shape when subjected to the surface loads g i (i.e. the solution of (1) with g replaced by g i ).

S(Ω) = C 1 (Ω) + C 2 (Ω), where C i (Ω) = Ω Ae(u i Ω ) : e(u i Ω ) dx,
First, we perform the minimization of S(Ω) under the volume constraint Vol(Ω) = V T , V T = 0.2|D|, without adding any constraint on the anisotropic maximum thickness of shapes (i.e. ( 28) is solved by replacing C(Ω) with S(Ω)). Doing so results in the optimized shape of Figure 18 (right).

One may argue that this shape is not particularly interesting, mainly due to the large bulk of material concentrated in the middle of the shape. Moreover, it is possible that the large surface in the (x, y) plane hinders the light entering from the z-direction. Therefore, we add a constraint on the anisotropic maximum thickness and solve the optimization problem [START_REF] Dombernowsky | Three-dimensional topology optimisation in architectural and structural design of concrete structures[END_REF] for δ = 0.024, by using different metric tensors of the form [START_REF] Dunning | Introducing the sequential linear programming level-set method for topology optimization[END_REF]. The optimized shapes appear in Figure 19 and the corresponding results in Table 4. Once more, the performance of all structures is quite similar, although the shapes are quite different. 3.

4.

5.

6.

7.

8.

9.

no thickness constraint m1=0.10 m1=0.01 In our last example, we show how a constraint over the anisotropic maximum thickness may serve as a manufacturing constraint for civil structures.

Iterations

To this end, let us consider the two-dimensional benchmark MBB beam example, of Figure 20, where only one half of the structure is considered, taking advantage of the symmetry of the problem (so that all In a first experiment, we minimize the volume Vol(Ω) of shapes, under a compliance constraint, that is:

min Ω Vol(Ω) s.t. C(Ω) ≤ C T ,
where the maximum authorized compliance C T equals 150. The result is displayed on Figure 21. The resulting, optimized shape shows a very long, horizontal bar in the upper side. Such members could be undesirable: because of their length, it is usually not possible to produce them in industry. In this kind of situation, engineers assemble such long bars by joining two (or more) smaller bars and constructing very strong connections between them (see Figure 22). Besides the increased cost of such connections, the final g = 0.5 beam may not behave mechanically in the predicted way, and so structural engineers generally wish to avoid them. , with a threshold value δ = 1.0, and where the anisotropy is dictated by the metric tensor:

M = m 1 0 0 m 2 .
The optimized shapes corresponding to the values m 1 = 0.1, m 1 = 0.01, and m 2 = 1 are presented on Figure 23: so that the distance in the horizontal direction is reduced, the upper feature is rotated and broken in two bars, while a third element is added at their joint. Recall that while the appearance of thinner, smaller bars could be a hindrance in other contexts, the typical scale involved in architectural is so large that even small members (on the reduced scale) are manufacturable. In this last section, we intend to deal with a more realistic example of how the proposed material in this article may help in the design of realistic civil structures.

The proposed test-case is inspired by the Pylon Design Competition organized by the Royal Institute of British Architects (RIBA) in 2014 [START_REF]of the Royal Institute of British Architects[END_REF]. This competition invited experts from various fields to search for innovating designs of electric pylons, with minimal impact on the surrounding environment, which are structurally coherent in the sense that they safely carry the loads transmitted by the cables.

Taking into account the specifications imposed by the organizers of the competition, and after several simplifications, the setting is that presented in Figure 24: the working domain has size 15 × 15 × 45; the structure is clamped on its lower side, and loads are applied at the seven points depicted in blue. Six load cases g i , i = 1, ..., 6 are considered, as represented in Figure 25: the first two correspond to extreme wind scenarii, while the last four correspond to situations where some wires have broken. We additionnally impose the symmetry of the structure with respect to the two planes containing the point (7.5, 7.5, 0) and parallel to the (xz)-and (yz)-planes respectively. As in the example of Section 6.3.2, the objective function of interest in this multiple load setting is the sum S(Ω) of the individual compliances C i (Ω) associated to each load case g i :

S(Ω) = 6 i=1 C i (Ω), where C i (Ω) = Ω Ae(u i Ω ) : e(u i Ω ) dx,
and u i Ω is the elastic displacement of the shape when the loads g i are applied (i.e. the solution of (1) with g replaced by g i ).

Several numerical experiments are conducted in this context, and the corresponding details are reported in Table 5.

•

• • • • • • x y z 0 @ 0 25 5 1 A 0 @ 0 0 50 1 A 0 @ 0 0 50 1 A 0 @ 0 0 50 1 A 0 @ 0 1 A 0 @ 0 1 A 0 @ 0 1 A (6)
Figure 25. Definition of the six load cases considered in the pylon design example of Section 6.4

6.4.1. Simplification of the optimal design and improvement of its visibility.

At first, we merely aim at minimizing S(Ω) under a volume constraint, i.e. we solve the optimization problem [START_REF] De Buhan | An optimization method for elastic shape matching[END_REF] with C(Ω) replaced by S(Ω), and with target volume V T = 0.1|D|. The resulting optimized shape is displayed in Figure 26. Although this shape looks structurally coherent, it contains several extended surfaces; as we have already mentioned, this is an undesirable feature, since it undermines the visibility of the shape, and makes it difficult to interprete the shape in terms of bars, which is a key step in the civil engineering design process.

One first idea towards circumventing this drawback simply consists in reducing the maximum authorized volume V T ; it is then expected that mechanically unnecessary extended surfaces will disappear. Acting this way on the target volume V T , which is given as a strict requirement in most industrial design situations, is perfectly legitimate in the present context: indeed, let us bear in mind that the designs resulting from the numerical optimization process are only conceptual outlines aimed at guiding the architectural design process. In this direction, we again solve Problem (28) (still with S(Ω) replacing C(Ω)) with target volumes V T = 0.075|D| and V T = 0.05|D|; the resulting shapes are displayed in Figure 27. Evidently, the structure gets simplified and the visual result is better.

An alternative way to achieve this simplification (thus allowing to explore different structural systems) consists in using a constraint on the anisotropic maximum thickness of shapes, in the spirit of Section 6.3, while keeping the initial target volume V T = 0.1|D|. More precisely, we now solve a version of Problem (32) The optimized shapes associated to the threshold values δ = 4 min {m 1 , m 2 , m 3 }∆x and δ = 3 min {m 1 , m 2 , m 3 }∆x (recall that ∆x is the mesh size) are represented in Figure 28. Notice that these shapes still present extended surfaces, which are not parallel to the two penalized planes. Interestingly, in the case where δ = 4 min {m 1 , m 2 , m 3 }∆x, a large quantity of material is concentrated at the bottom of the structure and the visual result is quite satisfying.

6.4.2. Imposing a resemblance between the optimal design and a target shape.

Let us now imagine that one desires a pylon design in which the material is concentrated around a central axis, as in the shape Ω T of Figure 29.

One possibility to achieve this purpose is to constrain shapes to stay close to Ω T , in the line of Section 5.1. More precisely, we solve Problem [START_REF] De Gournay | Velocity extension for the level-set method and multiple eigenvalues in shape optimization[END_REF] with target volume V T = 0.1|D|, and for different values of the weight t. The resulting optimized designs are displayed in Figure 30. Understandably enough, as the value of t increases, the optimized shape gets more and more influenced by the structural criterion S(Ω): its bending and torsional stiffness dramatically increases at its bottom side.

Exploring innovating designs by changing the boundary conditions.

Eventually, as exemplified in Remark 3.1, significantly different designs may be explored by acting on the imposed set of fixations. In this section, we slightly modify the definition of the test-case supplied in Figure 24 by shrinking the clamping surface Γ D into a 5 × 5 square centered at the middle of the bottom side of the working domain D. Grossly speaking, it is expected that doing so will entail optimized designs showing curved members, benefitting from the surrounding space of the working domain to connect Γ D to the supports of the loads.

In a first step, we simply minimize the mechanical criterion S(Ω) for several values of the target volume V T (that is, we solve Problem [START_REF] De Buhan | An optimization method for elastic shape matching[END_REF] with C(Ω) replaced by S(Ω)), with this new definition of the set Γ D of fixations for shapes. The associated optimized designs, represented in Figure 31, obviously show the aforementioned trend. Then, we incorporate two constraints on the anisotropic maximum thickness of shapes to the optimization problem: these constraints are identical to those detailed in Section 6.4.1, and are associated to the metric tensors defined in [START_REF] Feringa | An integral approach to structural optimization and fabrication[END_REF] (so that we solve exactly the same problem as in Section 6.4.1, up to the aforementioned change of boundary conditions). The resulting shapes are presented on Figure 32.

Concluding remarks

In this article, we have proposed a shape and topology optimization framework oriented towards applications in architecture. This framework emphasizes on a precise, geometric description of shapes: from the theoretical point of view, Hadamard's method is used to evaluate the sensitivity of optimization criteria with respect to variations of shapes, while their numerical representation is achieved via the level set method. Not only does this setting allow for an accurate account of the mechanical behavior of shapes, but it naturally lends itself to an efficient modelling of constraints of a geometric nature. In this perspective, we have devised three new shape functionals, which are either related to the mechanical performance of structures, or allow to add personal, aesthetic requirements into the shape optimization problem. In particular, our shape functionals make it possible to

• constrain shapes to resemble a specific design,

• constrain shapes to fit a user-defined pattern, • penalize the lengths of bars or the areas of extended surfaces of shapes, thus allowing to simplify their design, or to improve their visibility. Realistic numerical simulations in two and three space dimensions have been presented to confirm the efficiency of the presented framework. Ω. The following lemma will be used repeatedly in the following. Its proof is postponed to the end of this section.

Lemma A.1. Let Ω be a bounded domain of class C 2 and M be a symmetric, positive definite matrix. Then, (1) The unit normal vector n Ω M to ∂Ω M , pointing outward Ω M reads:

(36) ∀x ∈ ∂Ω, n Ω M (x M ) = com(M 1 2 )n Ω (x) |com(M 1 2 )n Ω (x)| = M -1 2 n Ω (x) |M -
∀x ∈ ∂Ω, M 1 2 II M x M M 1 2 = 1 |M -1 2 n Ω (x)| II x .
It follows immediately from [START_REF] Frattari | THE STRUCTURAL FORM: Topology Optimization in Architecture and Industrial Design[END_REF] that the skeleton Σ M is obtained from that Σ M of Ω M via Σ M = M -1 2 Σ M . Hence, for x / ∈ ∂Ω, d M Ω is differentiable at x if and only if x / ∈ Σ M ; using again [START_REF] Frattari | THE STRUCTURAL FORM: Topology Optimization in Architecture and Industrial Design[END_REF], one has then:

(37) p M ∂Ω (x) = M -1 2 p ∂Ω M (x M ).

Hence, for x / ∈ ∂Ω, x / ∈ Σ M , (16) follows and (i) is proved. Using again integration by parts on the boundary on the second term in the right-hand side of the above equality yields:

P p (Ω)(θ) = Γ κ|n Ω -n g | 2 θ • n ds - Γ (n Ω -n g ) • ∇ ∂Ω (θ • n) ds, = Γ κ|n Ω -n g | 2 -2κ(n Ω -n g ) • n Ω + 2(κ -div ∂Ω (n g )) θ • n ds,
which leads to the desired result.

Proof of Theorem 5.1. It follows from Proposition 4.6 that P M a (Ω) is shape differentiable at an arbitrary shape Ω ∈ U ad . An easy calculation then yields: We now proceed to achieve an expression of the above derivative with the convenient structure [START_REF] Allaire | Shape optimization with a level set based mesh evolution method[END_REF]. To this end, recalling that Σ M has zero Lebesgue measure (see Proposition 4.4), we apply the coarea formula to the mapping p M ∂Ω : Ω \ Σ M → ∂Ω, as a curvilinear version of the classical Fubini theorem (see [START_REF] Chavel | Riemannian geometry: a modern introduction[END_REF], or [START_REF] Allaire | Mutli-phase structural optimization via a level-set method[END_REF] in the same context): .

P M a (Ω)(θ) = ∂Ω M n M (x), θ (x)
This ends the proof.

Figure 1 .

 1 Figure 1. Front facade of the Qatar National Convention Center, in Doha (taken from [65]).

Figure 2 .Figure 3 .

 23 Figure2. Three optimized designs of a three-dimensional bridge associated to different sets of fixations (see[START_REF] Bando | Optimisation d'une structure et application architecturale[END_REF] for details).

Figure 4 .

 4 Figure 4. (Left) One shape Ω ⊂ R 2 , and (right) the graph of an associated level set function.

Definition 4 . 1 .

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

Figure 5 .

 5 Figure 5. Illustration of Definition 4.2 and Proposition 4.1: projection set of a point x ∈ Σ, projection point p ∂Ω (y) of y / ∈ Σ, and ray emerging from z ∈ ∂Ω.

4. 3 .

 3 The anisotropic, signed distance function.The material of Section 4.1 may be extended to the anisotropic context. Let M be a symmetric, positive definite d × d matrix (we do not discuss the case of space-dependent metric tensors, which is much more involved). In the anisotropic context, the space R d is equipped with the inner product •, • M and associated norm |•| M defined as: ∀x, y ∈ R d , x, y M = M x, y , and |x| M = x, x M , where •, • is the usual Euclidean inner product. Definition 4.3. Let Ω ⊂ R d be a bounded, Lipschitz domain. The signed, anisotropic distance function d M Ω

Figure 6 .

 6 Figure 6. (Left) One shape Ω (in black); (middle) isolines of the isotropic distance function d Ω ; (right) isolines of the anisotropic distance function d M Ω associated to the tensor (15).

Figure 7 .

 7 Figure 7. (Left) Reproduction of the definition of the test-case of Figure 2 modified by the addition of an 'empty' non-optimizable zone; (right) resulting optimal shape.

250 20 Figure 8 .

 208 Figure 8. Description of the two-dimensional bridge test case of Section 6.1.

Figure 9 .

 9 Figure 9. (Top) Initial shape, and (bottom) resulting optimized shape for the problem (28).

Figure 10 .

 10 Figure 10. (Top) 'target' domain Ω T provided by the designer; (bottom) optimized shape for the problem (28) starting from Ω T .

Figure 11 .

 11 Figure 11. Optimized shapes for the problem (29) setting: (a) t = 0.35, (b) t = 0.45, (c) t = 0.55, (d) t = 0.70. Weight t 0.7 0.55 0.45 0.35 Compliance C(Ω) 12927 13108 11558 12864 Value of P m (Ω) 81.48 60.76 17.96 -39.68

Figure 12 .

 12 Figure 12. Convergence histories for the optimization processes leading to the shapes in Figure 11: values of the objective function L(Ω) versus number of iterations.

Figure 13 .

 13 Figure 13. (Left) Isovalues of the 'pattern' function g defined in (30), (middle) initial shape, and (right) optimized shape for Problem (31).

Figure 14 .Figure 15 .

 1415 Figure 14. Boundary conditions in the mechanical test case of Section 6.2.

Figure 16 .

 16 Figure 16. (a) Description of the optimal bridge test case; (b) optimized shape without adding any thickness constraints; (bottom) optimized shapes for problem (32) setting: (c) m 1 = 0.1; (d) m 1 = 0.01.

6. 3 . 3 .

 33 Constraining the length of horizontal bars in a two-dimensional MBB beam.

Figure 17 .

 17 Figure 17. Evolution of (top,left) the compliance C(Ω), (top,right) the violation Vol(Ω) -V T of the volume constraint, and (bottom) the violation P M a (Ω)δ/2 of the thickness constraint during the optimization of the three-dimensional bridge examples of Section 6.3.1.

Figure 18 .

 18 Figure 18. (a) Description of the 3d short cantilever test case; (b) optimized shape without adding any thickness constraint.

Figure 19 .Table 4 .

 194 Figure 19. Optimized shapes for the short cantilever test case with different metric tensors M in the definition of the anisotropic maximum thickness (a) m 1 = 0.01; (b) m 1 = 0.01, m 2 = 0.01; (c) one constraint with m 1 = 0.01, m 2 = 1.0 and a second constraint with m 1 = 1.0, m 2 = 0.01; (d) m 3 = 0.01.

1 3Figure 20 .

 120 Figure 20. Setting of the MBB beam example of Section 6.3.3.

Figure 21 .

 21 Figure 21. (Left): Initial and (right) optimized shape obtained in the MBB beam example without adding any thickness constraint.

Figure 22 .

 22 Figure 22. Example of a junction between two HEA beams.

Figure 23 .

 23 Figure 23. Optimized shapes in the MBB beam example with a constraint on the anisotropic maximum thickness, and parameters (left) m 1 = 0.1; (right) m 1 = 0.01.

Figure 24 .

 24 Figure 24. Setting and boundary conditions of the pylon design example of Section 6.4.

Figure 26 .

 26 Figure 26. Resulting pylon from the minimization of S(Ω) with target volume V T = 0.1|D|.

Figure 27 .

 27 Figure 27. Optimized designs for the minimization of S(Ω) with target volumes (a) V T = 0.075|D| and (b) V T = 0.05|D|.

Figure 28 .

 28 Figure 28. Optimized design in the minimization of S(Ω) under a volume constraint V T = 0.1|D|, and anisotropic maximum thickness constraints associated to the parameter (a) δ = 4 min {m 1 , m 2 , m 3 }∆x; (b) δ = 3 min {m 1 , m 2 , m 3 }∆x (Section 6.4.1).

Figure 29 . 4 .M 1 2 1 2

 29411 Figure 29. Target shape Ω T in the example of Section 6.4.2.

Figure 30 .

 30 Figure 30. Optimized designs under the constraint that the shape resembles the design Ω T of Figure 29 (Section 6.4.2); values of the weighting coefficient (a) t = 0.2; (b) t = 0.3; (c) t = 0.4; (d) t = 0.6.

1 2 n

 2 Ω (x)| , where com(B) stands for the cofactor matrix of a d × d matrix B.(2) The second fundamental form II M of the boundary ∂Ω M satisfies the identity:

Figure 31 .

 31 Figure 31. Optimized designs in the minimization of S(Ω) performed in Section 6.4.3 with target volumes (a) V T = 0.1|D|; (b) V T = 0.075|D|; (c) V T = 0.05|D|.

  ∀θ ∈ Θ ad , P M a (Ω)(θ) = Ω G(d M Ω (x)) M n M (p M ∂Ω (x)), θ(p M ∂Ω (x)) dx.

d

  (z) ds(x),where we recall that ray M ∂Ω (x) = (p M ∂Ω ) -1 (x), and the Jacobian Jac(p M ∂Ω (x)) is defined as:Jac(p M ∂Ω (z)) = ∇p M ∂Ω (z)∇p M ∂Ω (z) T , z ∈ R d \ Σ M .The only remaining task is then to calculate this Jacobian, which we do by using the representation (18) of ∇p M ∂Ω in an orthonormal basis of R d :Jac(p M ∂Ω (z)) = p ∂Ω M (z))| λ i (p ∂Ω M (z)) -1

Table 2 .

 2

	Weight t	0	0.1	0.4	0.55	0.65	0.85
	Compliance C(Ω) 11.66 11.95 15.13 13.20 10.75 13.22
	Value of P p (Ω)	27.59 26.77 13.47 14.27 10.35 8.22

Table 2 .

 2

Acknowledgements. Part of this work has been supported by the OptArch project (H2020-MSCA-RISE-2015). The authors are indebted to A. Beghini, from Skidmore Owings & Merrill, for his help in the definition of the examples discussed in Remark 3.1 and Section 5. They also ackowledge fruitful discussions with F. Jouve (University Paris-Diderot), N. Aage (Technical University of Denmark) and O. Amir (Technion, Israel Institute of Technology). G. Allaire is a member of the DEFI project at INRIA Saclay Ile-de-France. A. Faure and R. Estevez also acknowledge support from the Labex CEMAM at Grenoble.

Appendix A. Mathematical details

In this appendix, we present the mathematical proofs of the results stated in Sections 4.3.1 and 5.

A.1. Proofs of the properties of the anisotropic, signed distance function. To obtain (ii), we simply apply [START_REF] Bajcsy | Multiresolution elastic matching[END_REF] to the domain Ω M , and rely on Lemma A.1; for x / ∈ ∂Ω, x / ∈ Σ M , the sequence of equalities

proves the first part of [START_REF] Burger | A framework for the construction of level set methods for shape optimization and reconstruction[END_REF]. As for the second part, one has:

Then, using [START_REF] Frattari | THE STRUCTURAL FORM: Topology Optimization in Architecture and Industrial Design[END_REF], [START_REF] Guo | Doing topology optimization explicitly and geometricallya new moving morphable components based framework[END_REF] and Lemma A.1, it follows:

which is the desired result.

Point (iii) is a straightforward consequence of the first two points (and notably the fact that Σ M = M -1 2 Σ M ), and Proposition 4.1. Likewise, point (iv) follows from Proposition 4.2, and the transformation rule of Lemma A.1 for the second fundamental form of Ω M .

Eventually, point (v) follows easily from a combination of (37), Lemma A.1, and the expression (13) of the derivative of the usual projection mapping.

We now provide the missing ingredient in the proof of Proposition 4.4.

Proof of Lemma A.1. The first point is classical (and actually holds in a much larger context than that of linear mappings), and we concentrate on that of the second one. To this end, we use the following elementary characterization: for a given tangent vector [START_REF] Guo | Explicit feature control in structural topology optimization via level set method[END_REF] together with [START_REF] Frey | Mesh generation[END_REF] yield:

)

and the desired result follows.

Proofs of Propositions 4.5 and 4.6. Again, both results may be equivalently inferred from their isotropic counterparts, by using the isomorphism x → M 1 2 x in the same way as in the proof of Proposition 4.4, or by adapting their original proofs (see [START_REF]Shape identification via metrics constructed from the oriented distance function[END_REF] and [START_REF] Dapogny | Optimisation de formes, méthode des lignes de niveaux sur maillages non structurés et évolution de maillages[END_REF]) to the present anisotropic context, and we omit further details about these points.

A.2. Details of the calculations of the shape derivatives of Section 5.

Proof of Proposition 5.2. The proof is carried out in the work [START_REF] Allaire | Structural optimization under overhang constraints imposed by additive manufacturing technologies[END_REF], and is reproduced in here for convenience.

For θ ∈ Θ ad , a change of variables in the boundary integral [START_REF] Christiansen | Topology optimization using an explicit interface representation[END_REF] yields (see e.g. [START_REF] Henrot | Variation et optimisation de formes: une analyse géométrique[END_REF] Prop. 5.4.3):

Using the well-known identity over matrices:

we obtain that θ → |com(I + ∇θ)| is Fréchet-differentiable at θ = 0, and

In this formula, we recall that div ∂Ω V := div(V )-∇V n, n stands for the tangential divergence of a (smooth enough) vector field V : ∂Ω → R d . Hence, using Lebesgue dominated convergence theorem and Proposition 4.3, it follows that θ → P (Ω θ ) is Fréchet-differentiable at 0, and that the corresponding derivative reads:

Using integration by parts on the boundary Γ (see [START_REF] Henrot | Variation et optimisation de formes: une analyse géométrique[END_REF], Prop. 5.4.9), together with the identity: