Filip Křikava
email: filip.krikava@inria.fr

Martin Monperrus
email: martin.monperrus@univ-lille1.fr

TTC'15 Live Contest Case Study Transformation of Java Annotations

Java 5 introduced annotations as a systematic mean to attach syntactic meta-data to various elements of Java source code. Since then, annotations have been extensively used by a number of libraries, frameworks and tools to conveniently extend behavior of Java programs that would otherwise have to be done manually or synthesized from external resources. The annotations are usually processed through reflection and the extended behavior is injected into Java classes using aspect-oriented techniques or a direct byte code modification. However, in some cases, class-level instrumentation might not always be available neither desirable and therefore the transformation is done at the source code level.

In this case study we focus on such source-level transformation. Concretely, the task is to inject behavior specified by an annotation library that encapsulates common programming concerns such as logging, caching and a failure retry. The objective is to explore how are the contemporary transformation tools suitable for programming language transformations.

Introduction

Java 5 has been extended with annotations that provide a convenient way to supply meta-data to various element of Java source code such as classes and methods. Since then, annotations have been extensively used by a number of Java libraries, frameworks and tools. One of the advantage of using annotations is that they attach syntactic meta-data directly to the relevant element of Java source code rather than decoupling them to external resources or setting them up manually through some API.

Among other things, annotations can be used to express crosscutting concerns such as logging or caching that would otherwise cluttered the source code as they are difficult to encapsulate using regular programming techniques [3]. Essentially, an annotation acts as a marker for a transformation tool which applies the expressed concern into the corresponding Java element (e.g., a class, a method). Majority of these tools rely on Java reflection API and some sort of class instrumentation using aspect-oriented programming (AOP) techniques or direct byte code manipulation to process the these annotations. While this is usually the preferred way, in some cases, the class level instrumentation might not available neither desirable. The problem is that it might (i) lead to leaky abstraction since the transformation is transparent to a developer, (ii) introduce some performance penalties, or (iii) bring some extra non-trivial dependencies.

TTC'15 Live Contest

In this case study we focus on the situation when class level instrumentation is not possible and develop a solution for annotation processing based on source level transformation. Concretely, the transformation task is to inject behavior specified by an annotation library that covers common programming concerns such as logging, caching and retry certain type of failures. The objective is to explore how the current model transformation tools are suitable for programming language transformations; particularly in comparison to some dedicated program transformation libraries such as Spoon [4]. The idea to bring together different communities on a common transformation problem.

In the next section we present the sample annotation library and describe the tasks of the case study in more details. After that, Section 3 overviews criteria that will be used to evaluate submitted solutions.

All supporting material, as well as this document, is available on the case study github page [START_REF]ASF: Simple Logging Facade for Java[END_REF] . The contestants are invited to submit an issue shall there be any ambiguity about the transformation tasks that are carried in this case study.

Case Description

As the annotation library, we use a simplified version of jcabi-aspects [2], a collection of Java annotations together with AOP-based processing tool that allows developers to conveniently apply some common crosscutting concerns into Java applications. In order to balance the development effort in this case study, we have selected the following annotations from the jcabi-aspects library: -@RetryOnFailure for an automated retry of failed method execution, -@Cacheable for a simple data caching of parameterless methods return values, and -@Loggable for an automated logging of method calls.

The core transformation involves traversing Java source code and extending the annotated methods with the behavior specified by these annotations. The details of each annotation including the specification of its arguments are provided in its associated javadoc.

In the following text we provide an overview of the individual transformation tasks that are part of this case study. The transformations are illustrated on an example rather than formally described which we believe should be sufficient [START_REF] Bugayenko | Useful AOP aspects[END_REF] given the nature of the transformation. In the following we discuss the details of the individual annotation transformations. The code for these transformation tasks is provided on the github page in the ttc15.tranj.examples package.

Running Example

To better illustrate the transformation tasks of this case study, we use the following example. Let us consider a class that is used to download a content of an URL. A possible [START_REF] Kiczales | Aspect-oriented programming[END_REF] When the basic functionality is implemented, we would like to extend it with the following features:

-Failure handling. The method should become more robust and accommodate for some of the inevitable network delays by retrying the download in the case an exception occurs. It should, however, only retry the call in the case when it makes sense-i.e. when it is possible to recover. For example, in the case of SocketTimeoutException exception, the execution should be retried (preferably after a delay), while UnknownHostException it should fail immediately.

-Caching. The URL should not be consulted every single time the method is called, but instead it should be keep the result of the last invocation in memory and reuse it for certain amount of time.

-Logging. Each call to the method should be logged. The logging message should present the state of the current instance at the method entrance as well as how long the invocation took at the method exit. All exceptions handled within the method body should be also logged.

Instead of coding this manually, we would like to use annotations to declaratively specify the above concerns, and have a transformation tool automatically synthesize code similar to the listing shown in Section A. Concretely, the only change to the original code should be the following annotations: The __getCacheLastAccessed and __getCacheContent defined on line 3 and 4 are cache bookkeeping fields used for storing the information about the last access time and method result respectively. The condition on line 10 checks the validity of the cache based on the annotation lifetime. Finally, the method exit point is replaced with the stored result on lines 24-25. Clearly, this have to be done for all method exit points in the case the annotated method have multiple ones.

Task 3: Logging

The last transformation task introduces logging. Any method annotated with @Loggable should have conditionally (depending on the annotation options) logged (i) its entry point, (ii) all its exit points, and (iii) all exceptions that are caught within the method body.

Annotating the URLDownload.get() method with @Loggable should produce the following code:

}Listing 1 :

 1 Java implementation is shown in Listing 1. public class URLDownload { private final URL url; public URLDownload(String url) throws MalformedURLException { this.url = new URL(url); } public byte[] get() throws IOException { try (InputStream input = url.openStream()) { ByteArrayOutputStream buffer = new ByteArrayOutputStream(); byte[] chunk = new byte[4 * 1024]; int n; while ((n = input.read(chunk)) > 0 ) { buffer.write(chunk, 0, n); } return buffer.toByteArray(); } } Basic version of URLDownload class.

@

  RetryOnFailure(attempts = 3, delay = 1000, retry = { SocketTimeoutException.class }, escalate = { UnknownHostException.class }) → @Cacheable(lifetime = 1000) @Loggable public byte[] get() throws IOException { ... } given period of time. Annotating the URLDownload.get() method with @Cacheable(lifetime = 1000) should produce the following code 4 : public class URLDownload { // added bookkeeping fields private long __getCacheLastAccessed = 0; private byte[] __getCacheContent = null; ... public byte[] get() throws IOException { // added condition if (System.currentTimeMillis() -__getCacheLastAccessed < 1000 && __getCacheContent != null) { → return __getCacheContent; } try (InputStream input = url.openStream()) { ByteArrayOutputStream buffer = new ByteArrayOutputStream(); byte[] chunk = new byte[4 * 1024]; int n; while ((n = input.read(chunk)) > 0) { buffer.write(chunk, 0, n); } // added __getCacheContent = buffer.toByteArray(); __getCacheLastAccessed = System.currentTimeMillis();

  logger private final org.slf4j.Logger __logger = org.slf4j.LoggerFactory.getLogger(URLDownload.class); → ...

https://github.com/fikovnik/ttc15-live-contest

Shall there be any ambiguity, the readers are invited to raise an issue on the case study github page.

In order to keep the case study concise and clear we omit the use of any external library.

Too keep things simple, we do not concern ourselves with invalidating the cache and thus freeing occupied memory.

Task 1: Retrying on a Failure

The first transformation tasks should handle the @RetryOnFailure annotation. The objective is to extend a given method with a simple failure handling strategy that retries failed invocations according to given options. Annotating the URLDownload.get() method with The original method code-i.e. what is between the lines 9-20 in Listing 1, is wrapped by another try-catch block (line 7) which is itself in a while loop (line 6). The number of attempts (3) is projected in the condition if (__retryCount > 3) (line 17) and the delay is translated in a current thread sleep on lines 21-25. Finally, we make a distinction in the exceptions that on which the call should be retried (the retry parameter -SocketTimeoutException) and the ones that should be escalated (the escalate parameter -UnknownHostException) in the generated catch clauses in line 10 and 13 respectively.

Task 2: Caching

The second transformation task concerns the @Cacheable annotation. The objective is to extend a given method with a simple caching strategy that keeps a result of a method invocation for a TTC'15 Live Contest Line 3 defines an instance of SLF4J logger [1]. Lines 9-12 and 25-27 insert the entry and exit point logging respectively. The entry point logging message additionally includes the state of the object-i.e. the value of all its declared (non-inherited) fields.

Because, there is no exception handling code-i.e. no catch clause, only entry and exit points are traced. If we would, however, apply the logging annotation on the method which is the output of @RetryOnFailure transformation, log messages should be also added right after every catch clause. The result should therefore look like: The code transformation should be performed in the declared order-i.e., first applying @RetryOnFailure, then @Cacheable and finally adding @Loggable. The expected result is shown in Section A.

Evaluation

The evaluation will be done in two parts. First, the submitted solutions will be evaluated by the case study authors on their correctness (the percentage of successful transformations) and performance (how long each transformation takes). Second, all the TTC'15 participants will judge the solutions on their conciseness, readability and the overall usability and suitability of the transformation tool for the given task.

The solution is expected to be packed in the way that it is runnable from a command line shell. The program should take two arguments, a path to a valid Java 5 source file input and a path to where the result transformation should be stored.

$ ./my-solution.sh URLDownload.java SynthesizedURLDownload.java

The main part of the evaluation is the assessment of the usability of a given transformation tool to the problem being addressed. Usability of a programming language, a library or a tool is difficult to assess as it tends to be subjective since it largely depends on the preferences and background of its users. We therefore leave it for the TTC attendees as they shall represent a mixture of preferences and backgrounds.


Franc ¸ois Jouve Grégoire Allaire
email: gregoire.allaire@polytechnique.fr

Georgios Michailidis
email: michailidis@cmap.polytechnique.fr

Molding direction constraints in structural optimization via a level-set method

In the framework of structural optimization via a level-set method, we develop an approach to handle the directional molding constraint for cast parts. A novel molding condition is formulated and a penalization method is used to enforce the constraint. A first advantage of our new approach is that it does not require to start from a feasible initialization, but it guarantees the convergence to a castable shape. A second advantage is that our approach can incorporate thickness constraints too. We do not adress the optimization of the casting system, which is considered a priori defined. We show several 3d examples of compliance minimization in linearized elasticity under molding and minimal or maximal thickness constraints. We also compare our results with formulations already existing in the literature.

Introduction

The increasing number of publications on industrial applications of shape and topology optimization reflects the interest of engineers to introduce these techniques in the design process of mechanical structures. Especially in case of complicated problems, where mechanical intuition is very limited, shape and topology optimization can serve as a valuable tool both in order to obtain an optimized structure and to accelerate the design process.

Among the several methods that appeared in the literature, such as SIMP (Solid Isotropic Material with Penalization) method [START_REF] Bendsøe | Optimal shape design as a material distribution problem[END_REF][START_REF] Zhou | The coc algorithm, part ii: topological, geometrical and generalized shape optimization[END_REF][START_REF] Bendsoe | Topology optimization: theory, methods and applications[END_REF], the homogenization method [START_REF] Allaire | Shape optimization by the homogenization method[END_REF][START_REF] Allaire | Optimal design for minimum weight and compliance in plane stress using extremal microstructures[END_REF][START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF], the phase field method [START_REF] Bourdin | Design-dependent loads in topology optimization[END_REF][START_REF] Takezawa | Shape and topology optimization based on the phase field method and sensitivity analysis[END_REF][START_REF] Zhou | Multimaterial structural topology optimization with a generalized cahn-hilliard model of multiphase transition[END_REF][START_REF] Blank | Relating phase field and sharp interface approaches to structural topology optimization[END_REF] or the Soft Kill Option [START_REF] Mattheck | Design and growth rules for biological structures and their application to engineering[END_REF][START_REF] Harzheim | Optimization of engineering components with the sko method[END_REF], the level-set method for shape and topology optimization [START_REF] Allaire | A level-set method for shape optimization[END_REF][START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF][START_REF] Osher | Level set methods for optimization problems involving geometry and constraints: I. frequencies of a two-density inhomogeneous drum[END_REF][START_REF] Sethian | Structural boundary design via level set and immersed interface methods[END_REF][START_REF] Wang | Structural shape and topology optimization in a levelset-based framework of region representation[END_REF] seems to fulfill industrial requirements in a satisfying way. Using a level-set function to describe implicitly the boundary of a shape [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[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] allows topological changes to appear in an easy way, while the geometric nature of the method is a benefit for the study of problems where the position of the interface plays a significant role (stress constraints, thermal problems with flux across the boundary, etc.). The method is independent of the objective function under study [START_REF] Allaire | A level-set method for vibration and multiple loads structural optimization[END_REF][START_REF] Allaire | Damage and fracture evolution in brittle materials by shape optimization methods[END_REF][START_REF] De Gournay | Velocity extension for the level-set method and multiple eigenvalues in shape optimization[END_REF][START_REF] De Gournay | Shape and topology optimization of the robust compliance via the level set method[END_REF] and the ability to adapt the mesh on the boundary [START_REF] Allaire | Topology and geometry optimization of elastic structures by exact deformation of simplicial mesh[END_REF][START_REF] Allaire | Shape optimization with a level set based mesh evolution method[END_REF][START_REF] Xia | A level set solution to the stress-based structural shape and topology optimization[END_REF] alleviates possible numerical difficulties due to the "ersatz" material or to the discontinuity of the material properties. Moreover, industrial design introduces significant constraints according to the fabrication method, the tooling limitations and the total cost that can be afforded. Some of them are essentially geometric constraints, related to a notion of local thickness. We have shown in our previous work [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF] that thickness can be explicitly controlled using a level-set method, which constitutes a great advantage for the industrialization of the method. Such constraints are of great significance for cast parts, i.e. structures that are intended to be constructed by casting. Casting [START_REF] Campbell | Castings[END_REF] is the fabrication process where molten liquid is poured into a cavity formed by molds. The final structure is obtained after solidification of the liquid and removal of the molds. Thus, the structure should have such a shape, so that the construction and the removal of the molds is possible without destroying either the structure or the mold. This is called the "molding constraint". The casting process imposes further specifications of mechanical nature on the shape of the structure, mainly related to the solidification and filling process. In [START_REF] Michailidis | Manufacturing Constraints and Multi-Phase Shape and Topology Optimization via a Level-Set Method[END_REF], we argued that such constraints can be translated, in the context of topology optimization, into geometrical constraints on the maximum and minimum allowable feature size, since the complete casting system (molding, solidification and filling system) is usually designed after the structural form definition. According to the choice of shape and topology optimization algorithm, different ways have been proposed to handle the molding constraint. In the framework of the SIMP method, which is a density method, Zhou et al. [START_REF] Zhou | Topology optimization with manufacturing constraints[END_REF] implemented a penalization scheme that favors higher densities at the lower part of the structure. Leiva et al. [START_REF] Leiva | An analytical directional growth topology parameterization to enforce manufacturing requirements[END_REF] have chosen to introduce directly the growth direction in the parametrization of the problem, while methods of topology control, such as connectivity and growth direction control, have been applied for the Soft Kill Option [START_REF] Harzheim | Topshape: An attempt to create design proposals including manufacturing constraints[END_REF]. A complete review of these methods and a comparison of results of topology optimization with and without manufacturing constraints can be found in [START_REF] Harzheim | A review of optimization of cast parts using topology optimization: Ii -topology optimization without manufacturing constraints[END_REF][START_REF] Harzheim | A review of optimization of cast parts using topology optimization: I -topology optimization with manufacturing constraints[END_REF]. In the framework of the level-set method, the only works on the topic -to our knowledge-are those of Xia et al. [START_REF] Xia | A level set based method for the optimization of cast part[END_REF][START_REF] Xia | Simultaneous optimization of cast part and parting direction using level set method[END_REF]. In [START_REF] Xia | A level set based method for the optimization of cast part[END_REF] the authors have proposed a molding condition on the design velocity, i.e. a modification of the descent direction that ensures the castability of the shape at each iteration, provided that the initial shape is also castable. In this work, the molding system is a priori defined. In [START_REF] Xia | Simultaneous optimization of cast part and parting direction using level set method[END_REF] the authors have added the optimization of the draw direction in the optimization problem. The same choice for the design velocity is done. Although the method allows those topological changes that do not come in conflict with castability, it is mentioned in [START_REF] Xia | A level set based method for the optimization of cast part[END_REF] that the shape cannot expand orthogonally to the casting direction. This is a great disadvantage in case one wants to impose a minimum thickness contraint.

In the present paper, we introduce a new approach to handle the molding constraint in the framework of the level-set method for shape and topology optimization. A pointwise constraint is formulated using the signed distance function and a penalty functional is then constructed to turn the constraint into a global one. A shape derivative [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] Mohammadi | Applied shape optimization for fluids[END_REF], [START_REF] Murat | Etude de problèmes d'optimal design[END_REF], [START_REF] Sokołowski | Introduction to shape optimization[END_REF] is calculated for this new functional and a simple penalization method is applied which guarantees that the optimal shape is castable at convergence. A first advantage of our new approach is that it does not require to start from a feasible initialization. This is of course a key feature since, in many industrial problems, it is very hard to find out a feasible design to start with. A second advantage is that our approach can incorporate thickness constraints, contrary to the previous method in [START_REF] Xia | A level set based method for the optimization of cast part[END_REF][START_REF] Xia | Simultaneous optimization of cast part and parting direction using level set method[END_REF]. The contents of our paper is as follows. Section 2 describes our model shape optimization problem. For simplicity we focus on compliance minimization with volume constraint: the main difficulty on which we shall focus is the addition of further molding and thickness constraints. Section 3 is a short review of the level-set method. Section 4 discusses the casting process while Section 5 introduces our new molding direction constraint. We also recall the approach of Xia et al. [START_REF] Xia | A level set based method for the optimization of cast part[END_REF][START_REF] Xia | Simultaneous optimization of cast part and parting direction using level set method[END_REF], as well as the "uniform cross-section surface constraint" of Yamada et al. [START_REF] Yamada | A topology optimization method based on the level set method incorporating a fictitious interface energy[END_REF], which simplifies a lot the shape of the desired molds. Section 6 is devoted to the computation of the shape derivatives of these molding constraints. Finally Section 7 features our 3-d numerical results which are obtained in the finite element software SYSTUS of ESI-Group [START_REF] Systus | User's manual[END_REF], which is well adapted to an industrial context. Our results were partially announced in [START_REF] Allaire | Casting constraints in structural optimization via a level-set method[END_REF].

Setting of the problem

Our goal is to optimize a shape Ω ⊂ R N (N = 2 or 3), a bounded domain occupied by a linear isotropic elastic material with Hooke's law A (a positive definite fourth-order tensor). Typically, the boundary of Ω is comprised of three disjoint parts, such that ∂ Ω = Γ D ∪ Γ N ∪ Γ 0 , with Dirichlet boundary conditions on Γ D , nonhomogeneous Neumann boundary conditions on Γ N and homogeneous Neumann boundary conditions on Γ 0 . We introduce a working domain D (a bounded domain of R N ) which contains all admissible shapes, that is Ω ⊂ D. The volume and surface loads are given as two vector-valued functions defined on D, f ∈ L 2 (D) N and g ∈ H 1 (D) N . The displacement field u is the unique solution in H 1 (Ω ) N of the linearized elasticity system

       -div (A e(u)) = f in Ω , u = 0 on Γ D , A e(u) n = g on Γ N , A e(u) n = 0 on Γ 0 , (1)
where e(u) is the strain tensor, equal to the symmetrized gradient of u. A classical choice for the objective function J (Ω ) to be minimized is the compliance (the work done by the loads). It reads

J(Ω ) = Ω f • u dx + Γ N g • u ds = Ω A e(u) • e(u) dx. (2)
A typical shape optimization problem is inf

Ω ∈U ad J (Ω ) , (3)
where U ad is the set of admissible shapes. Imposing that all shapes belong to the working domain D and that they satisfy a volume constraint 0 < V < |D|, a possible choice of admissible set is

U ad = {Ω ⊂ D such that |Ω | = V } . ( 4
)
As it is well-known [START_REF] Allaire | Conception optimale de structures[END_REF][START_REF] Bendsoe | Topology optimization: theory, methods and applications[END_REF], problem (3) may lack an optimal solution. Numerically, the non-existence of a minimizer of (3) is reflected by the fact that approximate numerical solutions are mesh dependent (the finer the mesh the more details or finer members in the solution). Classically, to obtain existence of optimal shapes, one needs to restrict further the admissible set U ad by imposing additional geometrical, topological or smoothness constraints [START_REF] Allaire | Conception optimale de structures[END_REF][START_REF] Murat | Etude de problèmes d'optimal design[END_REF][START_REF] Sokołowski | Introduction to shape optimization[END_REF].

In order to find a descent direction for advecting the shape, we rely on the Hadamard method of shape differentiation, following the approach of Murat and Simon [START_REF] Murat | Etude de problèmes d'optimal design[END_REF]. Starting from a smooth reference open set Ω , we consider domains of the type

Ω θ = Id + θ (Ω ), with θ ∈ W 1,∞ (R N , R N ).
It is well known that, for sufficiently small θ , (Id + θ ) is a diffeomorphism in R N .

Definition 1. The shape derivative of J(Ω ) at Ω is defined as the Fréchet derivative in W 1,∞ (R N , R N ) at 0 of the application θ → J (Id + θ )(Ω ) , i.e.

J (Id + θ )(Ω ) = J(Ω ) + J (Ω )(θ ) + o(θ ) with lim θ →0 |o(θ )| θ = 0 , where J (Ω ) is a continuous linear form on W 1,∞ (R N , R N ).
A classical result states that the shape derivative J (Ω )(θ ) depends only on the normal trace θ • n on the boundary ∂ Ω [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], [START_REF] Sokołowski | Introduction to shape optimization[END_REF]. We refer to [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF] for various examples of shapes derivatives in the elasticity setting, including that for compliance.

3 Level-set framework

Shape representation

We favor an Eulerian approach and use the level-set method [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[END_REF] to capture the shape Ω on a fixed mesh. Then, the boundary of Ω is defined by means of a level set function ψ (see Figure 1) such that

   ψ(x) = 0 ⇔ x ∈ ∂ Ω ∩ D, ψ(x) < 0 ⇔ x ∈ Ω , ψ(x) > 0 ⇔ x ∈ D \ Ω .
During the optimization process the shape is being advected with a scalar (normal) velocity V (x) derived from shape differentiation, as we will see in the sequel. The advection is described in the level-set framework by introducing a pseudo-time t ∈ R + and solving the well-known Hamilton-Jacobi equation

∂ ψ ∂t +V |∇ψ| = 0. ( 5
)
using 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].

Fig. 1 Level-set representation of a shape.

Signed distance function

We recall that if Ω ⊂ R N is a bounded domain, then the signed distance function to the boundary ∂ Ω is the function R N x → d Ω (x) defined by :

d Ω (x) =    -d(x, ∂ Ω ) if x ∈ Ω 0 if x ∈ ∂ Ω d(x, ∂ Ω ) if x ∈ R N \ Ω , where d(x, ∂ Ω ) is the usual Euclidean distance from x to ∂ Ω .
Very often, the Hamilton-Jacobi equation ( 5) is initialized, or re-initialized, with the signed distance function. However, at later times t, the level set function ψ(t, x), solution of (5), is not a signed distance function. Furthermore, the functions ψ and d Ω do not share the same boundary conditions (see [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF] for details). Therefore one cannot retrieve geometrical informations on the shape Ω (t) from ψ(t, x). However, at every time t it is not hard to compute the signed distance function of Ω (t). As in the case of thickness constraints [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF], we shall use this signed distance function to get all necessary information for the formulation of our molding constraints.

Ersatz material

Using the so-called "ersatz material" approach, we extend the state equations to the whole domain D. To do this, we fill the holes D \ Ω by a weak phase that mimics the void, but at the same time avoids the singularity of the rigidity matrix. More precisely, we define an elasticity tensor A * (x) which is a mixture of A in Ω and of the weak material mimicking holes in D \ Ω

A * (x) = ρ(x)A with ρ = 1 in Ω , 10 -3 in D \ Ω . (6)
Decomposing the boundary ∂ D of the working domain in three parts

∂ D = ∂ D D ∪ ∂ D N ∪ ∂ D 0 ,
and demanding that the shape boundary

∂ Ω = Γ D ∪ Γ N ∪ Γ 0 must further satisfy Γ D ⊂ ∂ D D , Γ N ⊂ ∂ D N ,
where ∂ D 0 supports homogeneous Neumann boundary conditions, the displacement u is finally computed as the solution of

       -div (A * e(u)) = f in D, u = 0 on ∂ D D , A * e(u) n = g on ∂ D N , A * e(u) n = 0 on ∂ D 0 . (7)

Casting process

We give in this section a short description of the casting process. A simplified sequence of steps for the construction of a cast part is the following:

1. Molds are used in order to create a cavity, having the shape of the structure that we intend to construct. 2. The cavity is filled with molten liquid metal. 3. The liquid solidifies. 4. The molds are removed and the cast part is revealed.

There are many different types of casting (metal casting, sand casting, investment casting, etc.) and the choice among them depends on the type of cast part. Each type inserts different constraints on the casting process. We address the interested reader to [START_REF] Campbell | Castings[END_REF] for a complete presentation of the casting process. Here we confine ourselves to permanent mold casting, in which the molds are removed without being destroyed. We call parting direction the direction along which one mold is removed and parting surface, the surface on which different molds come in contact [START_REF] Xia | A level set based method for the optimization of cast part[END_REF]. Note that several molds can be used in the casting system and each one has its own parting direction (see Figure 3). The parting surface between two molds can be predefined or it can be constructed after the optimization using suitable methods [START_REF] Ahn | Separating an object from its cast[END_REF][START_REF] Fu | The application of surface visibility and moldability to parting line generation[END_REF]. In most of the industrial applications, planar parting surface are preferred because of reasons of cost and simplicity [START_REF] Xia | A level set based method for the optimization of cast part[END_REF]. Each of the above steps introduces different constraints in the shape of the cast part. In this work we are mainly interested to ensure the feasibility of the last step, i.e. the removal of the molds. Thus, we need to impose that the cast part has such a shape, so that the molds can actually be removed after the end of the solidification process. Let us give a 2d example of the above mentioned. Suppose that for an optimization problem like the one described by equations ( 1) to ( 4) we obtain the optimized shape Ω , shown in Figure 2. In Figure 3 we see that depending on the molding system considered, this shape can be moldable or not. In the right figure of Figure 3, some parts of the shape oppose to the removal of the molds in their corresponding parting direction. The construction of the casting system is usually based on the intuition of the caster. Changes on the number and on the position of the molds can turn a non-moldable shape to a moldable one. The design of the whole casting system is very difficult (if possible) to be formulated mathematically and be subjected to optimization. Works in this direction are mostly concerned with parametric or shape optimization of parts of the molds [START_REF] Lewis | Thermal optimisation in the sand casting process[END_REF][START_REF] Tortorelli | Optimal design of nonlinear parabolic systems. part i: Fixed spatial domain with applications to process optimization[END_REF] or of the riser [START_REF] Tavakoli | Optimal riser design in sand casting process by topology optimization with simp method i: Poisson approximation of nonlinear heat transfer equation[END_REF]. In the present paper, we do not consider the optimization of the casting process, but the molding system is considered a priori defined.

5 Formulation of the molding direction constraint

Molding direction condition on design velocity

A molding direction condition on the design velocity was proposed by Xia et al. in [START_REF] Xia | A level set based method for the optimization of cast part[END_REF], which is inspired by Fu et al. [START_REF] Fu | The application of surface visibility and moldability to parting line generation[END_REF]. According to these authors, if a shape is feasible with respect to the molding direction specification for its corresponding molding system, then the boundary of the structure ∂ Ω can be divided into m disjoint parts Γ i , i = {1, ..., m}, such that Γ i ∩Γ j = / 0, j = {1, ..., m} , ∪ m i=1 Γ i = ∂ Ω and Γ i can be parted in the direction d i . Thus, a molding direction condition for this shape is

d i • n(x) ≥ 0, ∀x ∈ Γ i , (8)
where n(x) is the exterior unit normal at x ∈ Γ i . The shape on the left in Figure 4 satisfies the condition (8), while the shape on the right does not. In fact, as it is mentioned in [START_REF] Xia | A level set based method for the optimization of cast part[END_REF], undercuts (slots that hint the removal on the mold in its parting direction) and interior voids are not allowed.

Based on the the molding condition (8), Xia et al. [START_REF] Xia | A level set based method for the optimization of cast part[END_REF] proposed the following method : starting from a shape that satisfies the constraint [START_REF] Allaire | Casting constraints in structural optimization via a level-set method[END_REF], consider an advection velocity of the form

θ i (x) = λ (x)d i , ∀x ∈ Γ i . (9)
In this way, the shape remains always moldable, since no undercut can be created during the advection of the shape with this type of velocity and no interior void can be nucleated. The topological changes that can occur using this advection velocity cannot turn the shape from moldable to non-moldable [START_REF] Xia | A level set based method for the optimization of cast part[END_REF]. This method, despite its simplicity and effectiveness, presents two major drawbacks. First, the shape must be initialized as being castable so that it can satisfy the mold- ing constraint during the entire optimization process ; this is a severe limitation on the choice of admissible initial guess shapes. Nevertheless, if such an initialization can be found, then it turns out that the method is flexible enough, especially in 3d, in order for complicated topologies to appear from very trivial initializations. Second, and more important from our point of view, the very form (9) of the advection velocity does not allow all possible deformations of the shape, including those which are required for some other constraints. As it is stated in [START_REF] Xia | A level set based method for the optimization of cast part[END_REF], there is no component of the advection velocity normal to the parting direction. Therefore, the shape can shrink by extinction of some part, but it cannot expand normal to its corresponding parting direction. As an example, consider the case where a minimum thickness constraint is also applied [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF]. Then, if the measured thickness is in a direction orthogonal to the parting direction, the shape cannot expand in this orthogonal direction (in order to meet the constraint of minimal thickness) because it can move only parallel to its parting direction. Therefore, in such a situation, the thickness constraint will not be respected. Therefore, it is necessary to formulate a more general molding constraint, free of the above limitations.

Generalised molding constraint

A first idea for a generalized way to treat the molding direction constraint consists simply in regarding (8) as a pointwise constraint in our optimization problem. Then, it can be exactly penalized as we shall do in ( 16) to compute its shape derivative.

A second idea is to use the signed distance function to the boundary of the domain to derive all necessary information, as we have done for thickness constraints in [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF]. Denoting Ω the actual shape and D the design domain, a generalized molding direction constraint can be formulated as:

d Ω (x + ξ d i ) ≥ 0 ∀x ∈ Γ i , ∀ξ ∈ [0, dist(x, ∂ D)] , (10)
or equivalently

d Ω (x + ξ d i ) ≥ 0 ∀x ∈ Γ i , ∀ξ ∈ [0, diam(D)] , (11)
where we denote diam(D) = sup x,y {dist(x, y), x, y ∈ D} the diameter of the fixed domain D. We prefer to use formulation [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF] instead of [START_REF] Allaire | A level-set method for shape optimization[END_REF], in order to avoid the dependence of the term dist(x, ∂ D) on the shape Ω . Intuitively, this formulation says that, starting from a point on the boundary, which will be casted in the direction d i and travelling along this direction, we should not meet again some part of the structure (see Figure 5). In case that the parting surface is not defined a priori, but is revealed at a second step after the design has been completed and for a system of two molds (see Figure 3, right image), the constraint (11) becomes:

d Ω (x + ξ sign(n • d)d) ≥ 0 ∀x ∈ ∂ Ω , ∀ξ ∈ [0, diam(D)] . (12)

Uniform cross-section surface constraint

Another useful constraint for cast parts is the so-called "uniform cross-section surface constraint" [START_REF] Yamada | A topology optimization method based on the level set method incorporating a fictitious interface energy[END_REF], since it simplifies a lot the shape of the desired molds. To our knowledge, Yamada et al. [START_REF] Yamada | A topology optimization method based on the level set method incorporating a fictitious interface energy[END_REF] were the first to study this type of constraint in shape and topology optimization using a combination of a phase-field and a levelset method. The constraint states that the cast part should have a uniform constant thickness along some direction d. An example of a uniform cross-section cantilever of thickness h is given in Figure 6. The boundary conditions may not be uniform along this direction and therefore the problem cannot be reduced to a 2d problem. We can formulate this type of constraint at least in two ways. The first formulation states that the normal to the boundary cannot have a non-zero component in this direction d:

d • n(x) = 0, ∀x ∈ ∂ Ω \ ∂ D. (13)
A second way to enforce the constraint is to limit the admissible advection fields θ . Starting from an initial guess shape that has a uniform cross-section along the desired direction d and constraining the advection fields to be zero along this direction, the thickness along d will not change. In fact, this is the easiest way to follow, since no mathematical constraint is imposed in the optimization process and the calculation of the velocity field is reduced to a 2d problem, as we will see in the next section.

By enforcing the constraint (13), the feasibility of the shape is guaranteed for casting along the direction d, i.e. this constraint is a sufficient but not a necessary condition. 6 Shape derivative

Derivative of the condition on design velocity

Xia et al. proposed in [START_REF] Xia | A level set based method for the optimization of cast part[END_REF] a modification of the advection velocity according to [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF], that guarantees a descent direction. Starting from the general form of the shape derivative for a functional J(Ω )

J (Ω )(θ ) = ∂ Ω θ (s) • n(s)V (s)ds = m ∑ i=1 Γ i θ i (s) • n(s)V i (s)ds
and considering admissible advection fields of the type (9), we get

J (Ω )(θ ) = m ∑ i=1 Γ i λ i (s)d i • n(s)V i (s)ds,
and choosing

λ i (s) = -V i (s)d i • n(s), ∀i = 1, ..., m
for each part Γ i of the boundary ∂ Ω , the shape derivative becomes

J (Ω )(θ ) = - m ∑ i=1 Γ i (d i • n(s)) 2 (V i (s)) 2 ds ≤ 0,
which shows that the chosen advection velocity is indeed a descent direction. We replace the Hamilton-Jacobi equation ( 5) by the linear transport equation

θ i (s) = -V i (s)(d i • n(s))d i , ∀i = 1, ..., m (14)
∂ ψ ∂t + θ • ∇ψ = 0, (15)
where the vectorial velocity θ is an extension of the advection velocity [START_REF] Bendsøe | Optimal shape design as a material distribution problem[END_REF] and the normal n is the normal associated with the initial shape.

Derivative of the generalized molding constraint

We start with the derivation of constraint [START_REF] Allaire | Casting constraints in structural optimization via a level-set method[END_REF]. One advantage of this constraint is that it is of local nature, i.e. it contains information only for points on the boundary without searching along rays emerging from them. On the other hand, it contains the exterior normal vector, whose derivation is more complicated than the one of the signed distance function. In a first step, a global penalty functional can be formulated as

P GMC (Ω ) = ∂ Ω [(d • n(s)) -] 2 ds, (16)
with the usual notations ( f ) + = max ( f , 0) and ( f ) -= min ( f , 0).

Proposition 1. For a smooth shape Ω , the shape derivative of ( 16) reads

P GMC (Ω )(θ ) = ∂ Ω θ (s) • n(s) 2d • ∇ s (d • n(s)) --H(s)[(d • n(s)) -] 2 ds, ( 17
)
where H is the mean curvature and ∇ s the tangential gradient.

Proof. Using a classical result about shape derivation of integrals with shapedependent integrands (see Proposition 6.28 in [START_REF] Allaire | Conception optimale de structures[END_REF]), the shape derivative of ( 16) reads

P GMC (Ω )(θ ) = ∂ Ω θ (s) • n(s) H(s)[(d • n(s)) -] 2 + ∂ ([(d • n(s)) -] 2 ) ∂ n ds + ∂ Ω ∂ ([(d • n(s)) -] 2 ) ∂ Ω (θ )ds = ∂ Ω θ (s) • n(s) H(s)[(d • n(s)) -] 2 + 2(d • n(s)) -∂ (d • n(s)) ∂ n ds + ∂ Ω 2(d • n(s)) -d • n (s)(θ )ds = ∂ Ω θ (s) • n(s) H(s)[(d • n(s)) -] 2 + 2(d • n(s)) -d • ((∇n)n) ds + ∂ Ω 2(d • n(s)) -d • n (s)(θ )ds, (18)
where n (s)(θ ) is the shape derivative of the normal. Under the smoothness assumption on the shape, there exists an extension of the unit normal in a tubular area around the boundary by n(x) = ∇d Ω (x). Now, the unit normal satisfies the equa-tion |n(x)| 2 = 1 from which differentiating both sides, we obtain (∇n)n = 0. Thus, equation ( 18) reduces to

P GMC (Ω )(θ ) = ∂ Ω θ (s) • n(s)H(s)[(d • n(s)) -] 2 ds + ∂ Ω 2(d • n(s)) -d • n (s)(θ )ds. ( 19
)
What remains is the calculation of the shape derivative of the unit normal to the boundary. From Lemma 4.8 in [START_REF] Murat | Etude de problèmes d'optimal design[END_REF], we have that the transported of the unit normal n(Ω , x) is

n((Id + θ )(Ω ), x + θ (x)) = ((I + ∇θ ) -1 ) T n |((I + ∇θ ) -1 ) T n| = n -(∇θ ) T n + o(θ ) 1 -(∇θ ) T n • n + o(θ ) = (n -(∇θ ) T n + o(θ ))(1 + (∇θ ) T n • n + o(θ )) = n(Ω , x) -(∇θ ) T n + ((∇θ ) T n • n)n + o(θ ),
and so the Lagrangian shape derivative of the unit normal is

Y (θ , x) = -(∇θ ) T n + ((∇θ ) T n • n)n.
Since by the Hadamard structure theorem [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], [START_REF] Sokołowski | Introduction to shape optimization[END_REF], the shape derivative in the direction θ depends only on the normal component θ •n on the boundary ∂ Ω , we can restrict our attention to a vector field θ (x) of the form θ (x) = w(x)n(x), where w is any scalar function. In such a case, we find that (∇n

) T θ = w(x)(∇n) T n = 0 and thus Y (θ , x) = -(∇θ ) T n -(∇n) T θ + ((∇θ ) T n • n)n + ((∇n) T θ • n)n = -∇(θ • n) + n • ∇(θ • n) n = -∇ s (θ • n) = -∇ s (w(x)).
The Eulerian shape derivative of the unit normal reads

n (x)(θ ) = Y (θ , x) -∇nθ (x) = Y (θ , x) = -∇ s (w(x)).
The same result was found in [START_REF] Laadhari | On the equilibrium equation for a generalized biological membrane energy by using a shape optimization approach[END_REF], using similar variational principles. In view of the above results, equation [START_REF] Campbell | Castings[END_REF] becomes

P GMC (Ω )(θ ) = ∂ Ω w(s)H(s)[(d • n(s)) -] 2 ds - ∂ Ω 2(d • n(s)) -d • ∇ s w(s)ds.
On the other hand, using the identity (see [START_REF] Laadhari | On the equilibrium equation for a generalized biological membrane energy by using a shape optimization approach[END_REF])

∂ Ω a • ∇ s bds + ∂ Ω (∇ s • a)b ds = ∂ Ω a • n H b ds,
where a is a vector field and b is a scalar field, we deduce

P GMC (Ω )(θ ) = ∂ Ω w(s) ∇ s • (2(d • n(s)) -d) -H(s)[(d • n(s)) -] 2 ds = ∂ Ω w(s) 2d • ∇ s (d • n(s)) --H(s)[(d • n(s)) -] 2 ds,
which completes the proof.

Lemma 1. The shape derivative ( 17) can be also written in the form

P GMC (Ω )(θ ) = ∂ Ω θ (s) • n(s) N-1 ∑ i=1 κ i (s)(d • e i (s)) 2 (-sign((d • n(s)) -)) -H(s)[(d • n(s)) -] 2 ds,
where κ i are the principal curvatures of ∂ Ω at a point s ∈ ∂ Ω and e i the associated principal curvature directions (i = 1, ..., N -1).

Proof. For a point s ∈ ∂ Ω we can write ∇ s n(s) in the form (see [START_REF] Allaire | Mutli-phase structural optimization via a level-set method[END_REF]):

∇ s n(s) = N-1 ∑ i=1 κ i (s)e i (s) ⊗ e i (s). ( 20
)
Substituting ( 20) in ( 17) yields the desired result.

We now switch to the derivation of the other constraints ( 11) and ( 12), which are pointwise constraints of the same type as the minimum thickness constraint in [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF]. Therefore, the same steps need to be followed for their shape derivation and the final extraction of a descent direction. For the sake of completeness, let us mention once more the basic steps of this procedure.

For constraint [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF] we formulate a penalty functional of the form

P GMC (Ω ) = m ∑ i=1 Γ i diam(D) 0 (d Ω (s + ξ d i )) -2 dξ ds,
while for constraint [START_REF] Allaire | Damage and fracture evolution in brittle materials by shape optimization methods[END_REF], it reads

P GMC (Ω ) = ∂ Ω diam(D) 0 (d Ω (s + ξ sign(n(s) • d)d)) -2 dξ ds.
The two functionals are of the same type and can be written in compact notation (see Figure 7)

P GMC (Ω ) = ∂ Ω diam(D) 0 (d Ω (x m )) -2 dξ dx, (21)
where x m denotes an offset point of the boundary which is either

x m = x + ξ d i or x m = x + ξ sign(n(x) • d)d.
The following result was obtained in [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF] and we recall it without proof.

Proposition 2. The shape derivative of ( 21) reads

P GMC (Ω )(θ ) = ∂ Ω diam(D) 0 θ (x) • n (x) H (d Ω (x m )) -2 + 2 (d Ω (x m )) -∇d Ω (x m ) • ∇d Ω (x) dξ dx - ∂ Ω diam(D) 0 θ x m|Ω • n x m|Ω 2 (d Ω (x m )) -dξ dx,
where x m|Ω is the orthogonal projection of x m on ∂ Ω .

The advantage of the formula of Proposition 2, compared to that of Proposition 1, is that it does not contain any tangential derivative of the normal, or equivalently principal curvatures, which are notably hard to compute with great accuracy.

Remark 1. As already noticed in [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF], a descent direction can be found in a second step, after identifying the linear form of the shape derivative with another scalar product. The idea is similar to that of regularization, as described in [START_REF] De Gournay | Velocity extension for the level-set method and multiple eigenvalues in shape optimization[END_REF]. More precisely, solving the variational formulation

D α 2 reg ∇Q • ∇v + Qv dx = P (Ω )(v) ∀v ∈ H 1 (D), ( 22
)
where α reg > 0 is a positive number (of the order of the mesh size) which controls the regularization width, yields a solution Q ∈ H 1 (D). Then, choosing a vector field θ = -Q n, we obtain a guaranteed descent direction for P GMC since, taking v = -Q in ( 22), we get

P GMC (Ω )(-Q n) = - D α 2 reg |∇Q| 2 + Q 2 dx.
Fig. 7 Offset point x m and its projection x m|Ω on the boundary.

Derivative of the uniform cross-section constraint

For the constraint (13), a quadratic penalty functional reads

P UCS (Ω ) = ∂ Ω \∂ D (d • n(s)) 2 ds, (23)
which highly ressembles ( 16) and thus its shape derivation is omitted here. If we work in the feasible set of shapes which are constant in the direction d, it is even simpler to take into account this constraint. We consider shapes Ω = S × [0, h] where S is a surface perpendicular to d (see Figure 6). In this case, we force the advection velocity to be zero along the direction d of uniform thickness. Starting from the general formula of a shape derivative

J (Ω )(θ ) = ∂ Ω V (s)θ (s) • n(s)ds, where ∂ Ω = ∂ S × [0, h],
we use Fubini's theorem for the shape derivative

J (Ω )(θ ) = ∂ S h 0 V (ξ )θ (s) • n(s)dξ ds = ∂ S θ (s) • n(s) h 0 V (ξ )dξ ds.
From this, a descent direction is revealed for the uniform cross-section optimizable boundary ∂ S with the choice

   θ (s) = -n(s) h 0 V (ξ )dξ , ∀ s on ∂ S, θ (s) • d = 0, ∀ s on ∂ S, (24)
where n(s) is the normal to ∂ S which satisfies n • d = 0. Another, simple way to treat this constraint is through the regularization of the velocity field via equation [START_REF] Fu | The application of surface visibility and moldability to parting line generation[END_REF]. Choosing α reg to be a tensor, instead of a positive scalar, we can smooth the advection field in an anisotropic way. Then, equation ( 22) is rewritten as

D N ∑ i=1 a 2 i ∂ Q ∂ x i ∂ v ∂ x i + Qv dx = J (Ω )(v) ∀v ∈ H 1 (D), (25)
where α reg = ∑ N i=1 a i e i ⊗e i is the regularization tensor in the canonical basis (e i ) i=1,...,N of R N . For example, if we want a vector field θ of the type of ( 24) with d = e 2 , we can set a 2 >> a i , i = 2. Then, the solution Q of (25) will be constant in the x 2 direction and the descent direction θ = -Q n(s), where n(s) is the normal to ∂ S, satisfies θ • d = 0. In other words, starting from an initial shape that respects the constraint and regularizing the advection field in the way just described, we obtain a final optimized shape with a uniform cross-section.

Numerical examples

We have coded all numerical examples herein in the finite element software SYS-TUS of ESI-Group [START_REF] Systus | User's manual[END_REF]. A quadrangular mesh has been used both for the solution of the elasticity system and for the level set function. For the elasticity analysis, Q1 finite elements have been used, the Young modulus E is normalized to 1 and the Poisson ratio ν is set to 0.3. The "ersatz material" is considered to have the same Poisson ratio, while its Young modulus is set to 10 -3 .

Molding direction

The three-dimensional box-like structure of Figure 8 is our test case to apply several molding direction constraints and compare the corresponding optimized shapes. The entire domain is used for the analysis and is discretized using 40 × 40 × 20 Q1 elements. We minimize the compliance under an equality constraint for the volume. The optimization problem reads min

Ω ∈U ad ∂ Ω g • uds s.t. Ω dx = a V |D|, (26)
where u is the solution of ( 7) and a V ∈ [0, 1] determines the final volume of the structure as percentage of the volume of the working domain D. An augmented Lagrangian method is used here to enforce the constraints, as in our previous work [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF]. We refer to [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF] for the formula of the shape derivative for the compliance. At a first step, we impose no molding constraint and solve the optimization problem [START_REF] Harzheim | A review of optimization of cast parts using topology optimization: I -topology optimization with manufacturing constraints[END_REF] for a V = 0.2 using the arbitrary initialization of Figure 9(a). The optimized shape after 250 iterations is shown in Figure 9(b). Let us now solve the same optimization problem for a cast part that must comply with a predefined molding system. For example, if we want to use one mold in the design domain D, remove it in the direction d = (0, 0, 1) and impose the plane z = 0 to be a possible parting surface, then obviously the shape in Figure 9(b) is no more feasible. Of course, we cannot hope that starting with a different, even much simpler, initialization, we can obtain a castable optimized shape without enforcing a molding direction constraint. As we have mentioned before, in the absence of thickness constraints, we believe that the method of [START_REF] Xia | A level set based method for the optimization of cast part[END_REF], as described in section 5.1, gives quite satisfying results with a very simple implementation. The only restriction is that the initial design must satisfy the constraint. For later iterations it is enough to impose the molding direction condition on the design velocity. Starting with a full-domain initialization (see Figure 10(a)) and taking the initial level-set function equal to the signed distance function to the upper part of the domain, we choose an advection velocity of the type [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF], where d = (0, 0, 1), and we obtain the optimized shape of Figure 10(b). A comparison of the performance with the previous test case is proposed on Figure 11.

More flexibility in shape variations is given if the casting direction is set as d = (0, 0, 1) and no parting surface is imposed. In this case, the design domain D can contain two molds, one removed in the direction d and the second in the opposite direction (-d). The same full-domain initialization, with the signed distance function to the upper part of the domain, could be chosen, but this would unfortunately result in a final system with just one mold. Instead, it is more efficient to take the initial level-set function equal to the signed distance function both to the upper and lower part of the domain. The optimized shape is shown in Figure 12(b).

As expected, a completely different optimized shape is obtained if we change the casting direction. Separating the molds horizontally, in the direction d = (1, 0, 0) and imposing no specific parting surface, yields the optimized shape of Figure 14. In both Figures 12 and14 we see that topological changes can take place by "pinching a thin wall" [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF], even though we started from a full-domain initialization.

Molding direction & Maximum Thickness

A constraint on the maximum local thickness can be combined with the molding condition on the design velocity without any difficulty a priori. The reason is that the maximum thickness constraint gradient will be of uniform sign, tending always to reduce the thickness (and the volume) of the shape. As we have mentioned in section 5.1, when an advection velocity of the type ( 9) is chosen, the shape can shrink, but not expand normal to the casting direction. Adding a maximum thickness constraint to the test case of Figure 10, where the shape is casted along the direction d = (0, 0, 1) and the plane z = 0 is chosen as a possible parting surface, we solve the optimization problem min

Ω ∈U ad ∂ Ω g • uds s.t. Ω dx = a V |D|, P MaxT (Ω ) =    Ω f (d Ω (x))d Ω (x) 2 dx Ω f (d Ω (x))dx    1 2 ≤ d max /2, (27)
where d max = 0.2 and f is a regularization function that reads (see [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF]):

f (d Ω (x)) = 0.5 1 + tanh |d Ω (x)| -d max /2 α f d max /2 ,
α f > 0 being a parameter that controls the regularization of the constraint. Using the same initialization as in Figure 10(a) and enforcing z = 0 as the only possible parting surface, the optimized shape after 250 iterations and the convergence diagrams are shown in Figures 16 and17.

Molding direction & Minimum Thickness

Suppose now that we want to add a minimum thickness constraint with d min = 0.4 in the shape of Figure 10(b). The molding condition ( 9) is no more a suitable method to follow (see section 5.1) and we shall instead combine a minimum thickness constraint with the generalized molding constraint [START_REF] Allaire | A level-set method for shape optimization[END_REF]. The previously optimized shape is taken as an initial guess to solve the problem min

Ω ∈U ad ∂ Ω g • uds s.t. Ω dx = a V |D|, P 1 (Ω ) = P MinT (Ω ) = ∂ Ω d min 0 (d Ω (s -ξ n (s))) + 2 dξ ds = 0, P 2 (Ω ) = P GMC (Ω ) = ∂ Ω diam(D) 0 (d Ω (s + ξ d)) -2 dξ ds = 0, ( 28
)
without any condition on the advection velocity. An optimized shape for the optimization problem ( 28) is shown in Figure 18(b). The convergence diagrams for the penalty functionals P 1 and P 2 are shown in Fig. 19.

We now switch to a minimum thickness constraint of d min = 0.3 and to the case of two mold in the z-direction, as applied to the shape of Figure 12(b). In this case, the optimization problem reads Fig. 12 Plots of the optimized shape for the optimization problem [START_REF] Harzheim | A review of optimization of cast parts using topology optimization: I -topology optimization with manufacturing constraints[END_REF], setting d = (0, 0, 1) as casting direction, no a priori defined parting surface and using the molding direction condition [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF].

min

Ω ∈U ad ∂ Ω g • uds s.t. Ω dx = a V |D|, P 1 (Ω ) = P MinT (Ω ) = ∂ Ω d min 0 (d Ω (s -ξ n (s))) + 2 dξ ds = 0, P 2 (Ω ) = P GMC (Ω ) = ∂ Ω diam(D) 0 (d Ω (s + ξ sign(n • d)d)) -2 dξ ds = 0 (29)
Note that the constraint P 2 (Ω ) is different from that in [START_REF] Laadhari | On the equilibrium equation for a generalized biological membrane energy by using a shape optimization approach[END_REF]. We obtain the optimized shape of Figure 20. The convergence diagrams for the penalty functionals P 1 and P 2 are shown in Figure 21. 102.07 and no parting surface (Figure 12). With casting direction d = (0, 0, 1), no parting surface and minimum thickness constraint (Figure 20(bottom)).

105.87 With casting direction d = (1, 0, 0) 114.13 and no parting surface (Figure 14). With casting direction d = (0, 0, 1)

123.68 and parting surface at z = 0 (Figure 10(bottom)). With casting direction d = (0, 0, 1), parting surface at z = 0 134.68 and minimum thickness constraint (Figure 18(bottom)). With casting direction d = (0, 0, 1), parting surface at z = 0 143.65 and maximum thickness constraint (Figure 16).

Uniform cross-section

The 2 × 0.5 × 1 three dimensional cantilever of Figure 22, discretized by 40 × 10 × 20 Q1 elements, is chosen as test case to apply the uniform cross-section surface We now look for an optimized shape with a uniform cross-section along the yaxis. Starting from the initial shape of Figure 24(a), which has five uniform holes along this direction, we regularize at each iteration the velocity field for the advection of the shape in an anisotropic way, based on equation [START_REF] Harzheim | A review of optimization of cast parts using topology optimization: Ii -topology optimization without manufacturing constraints[END_REF], with a much higher regularization coefficient in the y-direction (a y >> a x , a z ). In our example, a x = a z = 2∆ x, ∆ x being the uniform mesh size, has been used to regularize the advection velocity in a small region around the shape boundary in the direction of the xand z-axis, while a y = √ 10 has been set to create a uniform velocity along this direction. The optimized shape is shown in Figure 24(b). The convergence diagrams for the compliance and the volume are shown in Figure 25. Acknowledgements The authors acknowledge fruitful discussions and helpful remarks from Marc Albertelli (Renault) and Charles Dapogny (Laboratoire Jean Kuntzmann, CNRS, Université de Grenoble). This work has been supported by the RODIN project (FUI AAP 13). G. Allaire is a member of the DEFI project at INRIA Saclay Ile-de-France. Fig. 16 Optimized shape for the optimization problem [START_REF] Henrot | Variation et optimisation de formes: une analyse géométrique[END_REF], with a maximum thickness constraint, setting d = (0, 0, 1) as casting direction and using the molding direction condition [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF].

Fig. 2

 2 Fig. 2 Possible optimized shape of a cast part.

Fig. 3

 3 Fig. 3 Left: moldable shape; right: non-moldable shape.

Fig. 4

 4 Fig. 4 Top: moldable shape; bottom: non-moldable shape.

Fig. 5

 5 Fig. 5 Checking castability along the parting direction d at the point x ∈ ∂ Ω .

Fig. 6 (

 6 Fig. 6 (a): uniform cross-section cantilever of thickness h; (b): cross-section S.

Fig. 8

 8 Fig. 8 Boundary and loading conditions for a 3d box.

Fig. 9

 9 Fig. 9 Initialization (top) and optimized shape (bottom) for the optimization problem (26) without a molding constraint.

Fig. 10

 10 Fig. 10 Initialization (top) and optimized shape (bottom) for the optimization problem (26), setting d = (0, 0, 1) as casting direction, z = 0 as a possible parting surface and using the molding direction condition[START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF].

Fig. 11

 11 Fig. 11 Compliance (top) and volume (bottom) convergence histories for the results in Figures 9(bottom) and 10(bottom).

Fig. 13

 13 Fig. 13 Compliance (top) and volume (bottom) convergence diagrams for the results in Figures 9(bottom) and 12.

  Fig. 13 Compliance (top) and volume (bottom) convergence diagrams for the results in Figures 9(bottom) and 12.

Fig. 14

 14 Fig.14Plots of the optimized shape for the optimization problem[START_REF] Harzheim | A review of optimization of cast parts using topology optimization: I -topology optimization with manufacturing constraints[END_REF], setting d = (1, 0, 0) as casting direction, no a priori defined parting surface and using the molding direction condition[START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF].

Fig. 15

 15 Fig. 15 Compliance (top) and volume (bottom) convergence diagrams for the results in Figures 9(bottom) and 14.

  Fig. 15 Compliance (top) and volume (bottom) convergence diagrams for the results in Figures 9(bottom) and 14.

Fig. 17

 17 Fig. 17Compliance (top) and volume (middle) convergence diagrams for the results in figures 10(bottom) and 16; convergence diagram for the maximum thickness functional (bottom) P MaxT (Ω ) for the optimized shape in Figure 16.

  Fig. 17Compliance (top) and volume (middle) convergence diagrams for the results in figures 10(bottom) and 16; convergence diagram for the maximum thickness functional (bottom) P MaxT (Ω ) for the optimized shape in Figure 16.

Fig. 18

 18 Fig. 18 Optimized shapes under: a molding constraint (top) and: a molding and minimum thickness constraint (bottom), with a predefined parting surface at z = 0.

Fig. 19

 19 Fig. 19 Convergence diagrams for the penalty functionals: P 1 (top) and: P 2 (bottom), for the results of Figure (18)(bottom).

Fig. 20

 20 Fig. 20 Optimized shapes under: a molding constraint (top) and: a molding and minimum thickness constraint (bottom), without a predefined parting surface.

Fig. 21

 21 Fig. 21 Convergence diagrams for the penalty functionals: P 1 (top) and: P 2 (bottom), for the results of Figure (20)(bottom).

Fig. 22

 22 Fig.22Boundary conditions for the "uniform cross-section" test case.

Fig. 23

 23 Fig. 23 Initialization (top) and optimized shape (bottom), without the "uniform cross-section" constraint.

Fig. 24

 24 Fig.24Initialization (top) and optimized shape (bottom), with a "uniform cross-section" constraint.

Fig. 25

 25 Fig. 25 Compliance (top) and volume (bottom) convergence diagrams for the results of figures 23 and 24.

Table 1

 1 Compliance of the optimized structures.

	Compliance


G Allaire

A review of adjoint methods for sensitivity analysis, uncertainty quantification and optimization in numerical codes

Keywords: adjoint method, sensitivity, uncertainty, optimization

The goal of this paper is to briefly recall the importance of the adjoint method in many problems of sensitivity analysis, uncertainty quantification and optimization when the model is a differential equation. We illustrate this notion with some recent examples. As is well known, from a computational point of view the adjoint method is intrusive, meaning that it requires some changes in the numerical codes. Therefore we advocate that any new software development must take into account this issue, right from its inception.

Introduction

The adjoint method originates in the theory of Lagrange multipliers in optimization. As an example, consider the minimization of a function J from R n into R, under the equality constraint F, a function from R n into R m , min J(x) such that F(x) = 0 (1) Introducing the Lagrange multilplier p in R m and the Lagrangian L(x,p) = J(x) + p.F(x) (2) the optimality condition for (1) (under some qualification conditions) is the stationarity of the Lagragian, namely J'(x*) + p*.F'(x*) = 0 (3) The adjoint method is an extension of this approach in the framework of optimal control theory. In this context, the variable x is the union of a state variable y and a control variable u, while the constraint F(y,u)=0 is the state equation which gives y in terms of u. In such a case the Lagrange multiplier p is called the adjoint state. This approach was developed by Pontryagin [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF] for ordinary differential equations and by Lions [17] for distributed systems, i.e., partial differential equations.

Adjoint method

Setting of the problem

For simplicity in the exposition we consider a finite dimensional model problem. Denote by u in R k a control or a parameter which is the ultimate optimization variable. The state of the system is denoted by y in R n and is defined as the solution of the following state equation A(u) y = b (4) where b is a given right hand side in R n and A(u) is an invertible n by n matrix. Since the matrix depends on u, so does the solution y. The goal is to minimize, over all admissible controls u, an objective function J(u,y) under the constraint (4). The difficulty of the problem is that y depends nonlinearly on u. We assume here that u is a continuous variable (the case of discrete variables is much more subtle). To numerically minimize the objective function J, the most efficient algorithms are by far those based on derivative informations. Therefore, a key issue is to compute the gradient of J(u,y(u)). This is the purpose of the next section.

Computation of a gradient

By the chain rule lemma the gradient of j(u)=J(u,y(u)) is j'(u) = ∂uJ(u,y(u)) + y'(u) ∂yJ(u,y(u)) ( 5) Differentiating (4) with respect to u yields A(u) y'(u) = -A'(u) y(u) (6) From ( 6) and [START_REF] Allaire | Structural optimization using sensitivity analysis and a levelset method[END_REF] we deduce the so-called direct gradient j'(u) = ∂uJ(u,y(u)) -A(u) -1 A'(u) y(u) ∂yJ(u,y(u)) (7) If the dimension k of the vector variable u is small, formula [START_REF] Allaire | Flow Control of Curved Air Ducts using Topological Derivatives[END_REF] is viable in the sense that it is not too costly for computing j'(u). Indeed, using [START_REF] Allaire | Flow Control of Curved Air Ducts using Topological Derivatives[END_REF] requires k linear solving of (6) (recall that the notation ' is the gradient with respect to u, a vector in R k ). Solving (6) is a costly operation since A(u) is an n by n matrix and n is usually very large (it is the number of degrees of freedom in the state equation, ranging from 10 4 up to 10 7 in the most recent applications). However, if k is as large as n (like in structural optimization, parameter identification, data assimilation, etc.), then formula [START_REF] Allaire | Flow Control of Curved Air Ducts using Topological Derivatives[END_REF] is useless ! Rather one should use the adjoint method which goes as follows. Introduce the so-called adjoint equation A T (u) p = -∂yJ(u,y(u)) (8) where A T (u) is the adjoint or transposed matrix. Multiplying (8) by y'(u) leads to A T (u)p.y'(u) = -y'(u).∂yJ(u,y(u)) (9) while multiplying (6) by p gives A(u)y'(u).p = -A'(u)y(u).p [START_REF] Bendsoe | Topology Optimization. Theory, Methods, and Applications[END_REF] Comparing ( 9) and ( 10) which have equal left hand side yields to a simplification of (5) as j'(u) = ∂uJ(u,y(u)) + A'(u)y(u).p (11) On the contrary of ( 7), formula ( 11) is cheap to evaluate, whatever the size of k: it requires only to solve the additional adjoint equation [START_REF] Ammari | An introduction to mathematics of emerging biomedical imaging[END_REF]. It largely outperforms [START_REF] Allaire | Flow Control of Curved Air Ducts using Topological Derivatives[END_REF] and it is thus the method of choice for computing the gradient of the objective function j(u)=J(u,y(u)). There are at least two ways for introducing an adjoint in a computer code. Either, a so-called analytic adjoint, i.e. ( 8), is implemented. Or a program for the adjoint is obtained by automatic differentiation [14], [START_REF] Hascoët | The Tapenade automatic differentiation tool: principles, model, and specification[END_REF] of the code solving the state equation ( 4). The only drawback of (11) (and of the whole adjoint approach) is that it is intrusive from a computational software point of view, namely it requires the development of an adjoint solver. This is quite easy if you know well your software. It is impossible if your software is a black box with no access to the source code. Even, if the source files are available, it could be a nightmare if they are not documented and/or not written in a clear and systematic manner... Figure 1: Optimal mast obtained with the level set method for shape and topology optimization

Applications

There are many applications of the adjoint approach.

The most classical ones are in control theory [17], [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF], optimization [START_REF] Allaire | Conception optimale de structures[END_REF], sensitivity or perturbation analysis, uncertainty quantification, inverse problems [START_REF] Tikhonov | Solutions of ill-posed problems[END_REF] and data assimilation [15]. I will review some of these applications, in connection with my own experience, therefore not trying to be exhaustive. The goal of this list is to convince the reader that, whatever the field of application, any numerical simulation code will sooner or later need an adjoint. Therefore, it is crucial to design any new software with this view in mind.

Optimal design

A classical application of the adjoint approach is the optimal design of structures or systems. The number of optimization variables describing the system is so large that it is the only viable approach. For example in fluid mechanics, the number of variables is related to the number of parameters, necessary to describe the boundary of an airplane [13], [START_REF] Mohammadi | Applied shape optimization for fluids[END_REF]. In the context of topology optimization for mechanical structures, the number of design variables is even larger since any cell of a "hold-all" computational domain is potentially a design variable: either it is full of material or empty. For any topology optimization method (be it homogenization [START_REF] Allaire | Shape optimization by the homogenization method[END_REF], SIMP [START_REF] Bendsoe | Topology Optimization. Theory, Methods, and Applications[END_REF], phase field [START_REF] Chambolle | Design-dependent loads in topology optimization[END_REF] or level set [START_REF] Allaire | Structural optimization using sensitivity analysis and a levelset method[END_REF]) the number of design variables is so large that again the adjoint method is the only way for computing a gradient and perform efficient optimization. Topology optimization (using an adjoint approach) has been very successful those last years and we present two results obtained with the level set method [START_REF] Allaire | Structural optimization using sensitivity analysis and a levelset method[END_REF], [START_REF] Allaire | Shape optimization with a level set based mesh evolution method[END_REF] in Figures 1 and2.

Figure 2: Optimal chair obtained with the level set method and an exact mesh of the structure Figure 3: Topological derivative or sensitivity for adding a vortex generator in an air duct

Perturbation analysis

The adjoint approach is also very helpful for sensitivity analysis, i.e., in the determination of the effect of small perturbations of the data on the solution. This is a very old and classical topic in nuclear reactor physics [START_REF] Planchard | Méthodes mathématiques en neutronique[END_REF] where the adjoint state is called the importance function. It is a classical tool in meteorology too [START_REF] Errico | What is an adjoint model ?[END_REF].

Figure 4: Boundary conditions for the mast test case

Compared to the setting of section 2, the variable u is now interpreted as a coefficient or a data in the model (4) and j(u)=J(u,y(u)) is a quantity of interest, the sensitivity of which has to be computed. Formula ( 11) is precisely the prediction of the sensitivity of the output j(u) with respect to perturbations of the data u.

Figure 5: Optimal mast without uncertainties A variant of this perturbation analysis appeared recently under the name of topological derivative [START_REF] Garreau | The topological asymptotic for PDE systems: the[END_REF], [START_REF] Sokolowski | On the topological derivative in shape optimization[END_REF]. It amounts to compute a criterion which predicts if the nucleation of a small hole or inclusion is going to decrease or not the objective function. It has been very successful, as an additional topology optimization tool, in structural mechanics [START_REF] Allaire | Structural optimization using topological and shape sensitivity via a level set method[END_REF], inverse problems [START_REF] Ammari | An introduction to mathematics of emerging biomedical imaging[END_REF] and fluid mechanics. In figure 3 an example from [START_REF] Allaire | Flow Control of Curved Air Ducts using Topological Derivatives[END_REF] is displayed: the blue (negative) region indicates where it is advantageous to put a vortex generator (a small winglet) on the boundary of an air duct in order to obtain a homogeneous flow distribution at the outlet.

Uncertainty quantification

Estimating the variations of a computational results in terms of possible uncertainties of the data or parameters is a very important topic which has attracted a lot of attention those last years. Very often, this task is associated to some probabilistic representation of the parameters and of the results. Besides the direct approach of Monte-Carlo simulation (which is very CPU time-consuming) one can use spectral methods or polynomial chaos [16]. Still these methods are quite intensive in terms of CPU. Another possibility is to rely on a worst-case approach which makes the problem deterministic if the worst-case scenario can be identified (which is not always possible). Recently a linearization approach for small uncertainties was suggested in [START_REF] Allaire | A linearized approach to worst-case design in parametric and geometric shape optimization[END_REF] and [START_REF] Babuška | Worst case scenario analysis for elliptic problems with uncertainty[END_REF] which makes the problem fully deterministic and explicit. Here is the main idea. Suppose that small perturbations δ occur in the right hand side of (4), namely A(u) y(u,δ) = b + δ (12) The perturbed solution is y(u,δ) = y(u) + z with z=A(u) -1 δ. Let us assume that the quantity of interest is j(u,δ) = J(u,y(u,δ)) (13) Assuming that the norm of the perturbation is bounded by a small number m>0, the worst-case function is w(u) = max|| δ|| < m j(u,δ) (14) Evaluating ( 14) is a difficult task in general. However, if m is small we can linearize (13) and compute (14) as w(u) ≈ j(u,0) + max|| δ|| < m ∂y J(u,y(u,0)) . z (15) which yields, since z=A(u) -1 δ, w(u) ≈ j(u) + m || A(u) -T ∂y J(u,y(u)) || (16) As usual inverting the matrix A T (u) is too costly, so we resort to the adjoint approach again. Multiplying the adjoint equation ( 8) by z and comparing with the equation A(u)z=δ multiplied by p leads to -∂y J(u,y(u,0)) . z = δ . p (17) which in turn gives the worst-case function w(u) ≈ j(u) + m || p || (18) In other words the worst case deviates from the standard ouput by the norm of the adjoint p multiplied by the magnitude m of the uncertainties. This is a very simple (albeit linearized) way of assessing and quantifying uncertainties which, again, relies on the notion of adjoint.

In [START_REF] Allaire | A linearized approach to worst-case design in parametric and geometric shape optimization[END_REF] we used formula [START_REF] Mohammadi | Applied shape optimization for fluids[END_REF] to perform shape and topology optimization of mechanical structures under various type of perturbations (on loads, material properties or geometry). In Figures 4, 5 and 6 are displayed boundary conditions and optimal shapes for minimal compliance under a volume constraint when uncertain vertical loads are applied in the blue zone. As can be seen here, these perturbations can have a drastic effect on the optimal topology.

Conclusion

In this paper we briefly recalled the adjoint theory for computing derivatives or sensitivities of an objective function depending on the solution of a state equation. In an increasingly large variety of cases (as illustrated in the figures above), the number of variables is so large that there is no viable alternative to the adjoint method. We also collected various problems where the adjoint is a key technique: of course optimal design, control and inverse problems, but also sensitivity analysis and uncertainty quantification. These types of problems occur in virtually all fields of science and engineering. Despite its efficiency the adjoint approach suffers from one major drawback which is the difficulty of implementing an adjoint solver (be it analytic or by automatic differentiation) in a preexisting code which was not designed for accommodating an adjoint solver. Note however that programming an adjoint along the early development of a numerical software is a rather easy task. Therefore we advocate that any new software planning should incorporate, in its early stages, all the necessary tools for later incorporating an adjoint (when the need will be felt and, surely, it will happen at some times). These ingredients include various data structures, the possibility of defining "dummy" loads, similar to the right hand sides of the adjoint equation ( 8), the definition of the linearized adjoint operator (if an analytic adjoint is going to be chosen) or the "clean" programming of the routines that will be submitted to an automatic differentiation tool like Tapenade [START_REF] Hascoët | The Tapenade automatic differentiation tool: principles, model, and specification[END_REF]. These minor precautions will later save a huge amount of development time, not to speak of the even larger savings in computational efficiency compared to any other method, either not relying on gradients or computing mere approximations (like finite differences).

Figure 6 :

 6 Figure 6: Optimal mast with uncertainties

  elasticity case", SIAM J. Control Optim., 39(6), 1756-1778, 2001. [13] M. Giles, N. Pierce: "An Introduction to the Adjoint Approach to Design", Flow, Turbulence and Combustion 65: 393-415, 2000. [14] A. Griewank, A. Walther: "Evaluating derivatives. Scientific computation, Springer, New York, 2010. [17] J.L. Lions: "Optimal control of systems governed by partial differential equations", Die Grundlehren der mathematischen Wissenschaften, 170, Springer-Verlag, 1971.

	Principles	and	techniques	of	algorithmic
	differentiation", Second edition, Society for
	Industrial and Applied Mathematics (SIAM), 2008.
	[15] F.X. Le Dimet, O. Talagrand: "Variational
	algorithms for analysis and assimilation of
	meteorological observations: theoretical aspects",
	Tellus A, 38A(2), 97-110, 1986.
	[16] O. Le Maître, O. Knio: "Spectral methods for
	uncertainty quantification. With applications to
	computational	fluid	dynamics",

Acknowledgement

Part of this work was performed in the framework of the RODIN project (sponsored by FUI AAP 13). This work was done in collaboration with many colleagues, among them M. Albertelli, Ch. Dapogny, G. Delgado, P. Frey, F. Jouve, G. Michailidis. The author is also a member of the DEFI project at INRIA Saclay Ile-de-France.