V V Rumyantsev
email: vvrum@ccas.ru.

Variational principles for systems with unilateral constraints

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.

1. Some definitions [START_REF] Suslov | Theoretical Mechanics[END_REF][START_REF] Chetayev | Theoretical Mechanics[END_REF] Let Oxyz be an inertial system of coordinates relative to which a system of point masses M with masses m ( = 1, . . ., N) is moving. If the points M of the system at an arbitrary instant of time t are unable to occupy an arbitrary position in space, or are unable to have arbitrary velocities, then such a system is said to be constrained. The conditions imposed on the positions or motions of a constrained system are called constraints. If the constraints do not allow the system to occupy arbitrary positions, they are called geometric or finite constraints. If the constraints only do not allow the points of the system at a given instant of time to have arbitrary velocities, the constraints are called kinematic or differential constraints. Finite constraints, and also differential integrable constraints, are called holonomic constraints, while differential non-integrable constraints are called non-holonomic constraints. If the constraints do not depend on time explicitly, they are called stationary or scleronomic constraints, and if they do, they are called non-stationary or rheonomic constraints.

Constraints are said to be bilateral if analytically they are expressed by equations of the form When the left-hand sides of these relations are positive, the constraints are said to be loosened, or that the system has left its constraints; when, however, the left-hand sides are equal to zero, the constraints are active, or strained. Constraints (1.1) and (1.2) will be assumed to be independent and smooth. We will confine ourselves below to linear differential constraints where B ␤ and D ␤ are certain functions of time and the coordinates. Here and below, summing over is carried out from = 1 to = N, summing over ␣ is carried out from ␣ = 1 to ␣ = a and summing over ␤ is carried out from ␤ = 1 to ␤ = b.

Elementary displacements ␦r , satisfying the conditions

(1.3)
at a fixed value of time (␦t = 0), are called virtual displacements of points of the system, and here the equality sign only occurs for non-loosening displacements which leave the system constrained.

Note that here there is a considerable defference between the virtual displacements in the cases of bilateral and unilateral constraints. In the former case, the virtual displacements are always reversible, whereas in the latter case they are irreversible, with the exception only of virtual displacements which leave the system constrained; this explains the double sign -equal to or greater than zero. Figuratively speaking, in the case of bilateral constraints, the system lies somewhere within configuration space, and its virtual displacements can occur in any direction submitted by the constraints, whereas in the case of unilateral constraints, the system lies on the boundary of the configuration space, as a result of which, when the virtual displacements are directed inwards, displacements in the opposite direction are impossible as they would lead beyond the configuration space. [START_REF] Lanczos | The Variational Principles of Mechanics[END_REF] For the motion of a mechanical system under the action of prescribed forces F applied, the constraints imposed on the system act at different points M of the system with certain forces R called constraint forces.

Axiom of constraints

The action of constraints imposed on a mechanical system can be replaced by the action of constraint forces R .

Corollary

If to the prescribed forces F acting on points of the system, constraint forces R are added, then the system can be thought of as being free of the constraints causing the constraint forces R .

The constraints imposed on the system depend on the physical nature of the mechanisms implementing the constraints, and therefore a constraint characteristic should be introduced into the mechanics in the form of an axiom that fixes realistic relations observed in experiments: In many cases the concept of ideal constraints is adopted as such an axiom. By definition, the total elementary work of ideal constraint forces R on an arbitrary virtual displacement of the system ␦r is equal to zero in the case of bilateral constraints or equal to or greater than zero in the case of unilateral constraints (1.4) and here the equality sign occurs only for non-loosening displacements.

We will now find expressions for the constraint forces corresponding to ideal constraints. We will first consider the non-loosening virtual displacements ␦r for which, in relations (1.3), the equality signs apply. We will multiply each of Eq. (1.3) by arbitrary multipliers -␣ and -␤ , sum them and then add them to Eq. (1.4); as a result, we obtain the equation from which we find expressions for the ideal constraint forces (1.5) where the number of constraint multipliers ␣ and ␤ is equal to the number a + b of constraint equations (1.1). In the sums on the right-hand side of Eq. (1.5), each term represents the action of one respective constraint on a given point M (1.6) In the case of unilateral constraints, we will retain expressions (1.6) for the constraint forces, on account of which the expression for the virtual work of the unilateral constraint forces coincides in form with the respective expression for bilateral constraints, i.e., whence, by virtue of relation (1.4), there follow additional conditions for the multipliers

(1.7)
It should be pointed out that, in the case of loosening virtual displacements, the left-hand side of relation (1.4) represents the elementary work of the constraint forces only in the nominal sense, namely if it is assumed that, during the displacement, the constraint forces retain their initial values. In this sense, inequality (1.4) can be regarded as indicating the relationship between the directions of the displacements and the constraint forces, rather than the work of the constraint forces. The work of an ideal unilateral constraint force on virtual displacement is always equal to zero: when the displacement leaves the system constrained, then but when the displacement causes system constraint violation, i.e. when the corresponding reactions vanish, and consequently ␣ = 0 and ␤ = 0 (␣ = 1, . . ., a; ␤ = 1, . . ., b). [START_REF] Suslov | Theoretical Mechanics[END_REF] In the literature, non-holonomic unilateral contraints are considered in two versions. In the first of these, the conditions for the coordinates and velocities are written in the form of inequalities, as in the present paper. In the second version, inequalities link the constraint forces, and the coordinates and velocities satisfy the constraint equations as long as the constraint forces belong to some open set in constraint force space. As soon as the constraint force reach the boundary of this set, the constraints disappear. This version is not considered in the present paper.

The virtual work principle, the fourier inequality

The equilibrium position of a system under the action of given forces is that position of the system in which it may for an indefinite time be at rest relative to a given reference system. The equilibrium positions of systems are usually considered when the prescribed forces F and constraints do not depend on time explicitly.

Let us complete the prescribed forces F with all constraint forces R ; then the mechanical system, according to the axiom of constraints, can be thought of as a system of free points under the action of forces F + R . For free points, the equilibrium equations have the form [START_REF] Newton | Mathematical Principles of Natural Philosophy (Philosophiae naturalis principia mathematica)[END_REF] We multiply these equations scalarly by the virtual displacements ␦r of points of the system, sum them and, taking into account inequality (1.4), obtain the relation (2.1) which holds for arbitrary virtual displacements ␦r from the equilibrium position. Relation (2.1) is not only necessary but also sufficient [START_REF] Suslov | Theoretical Mechanics[END_REF][START_REF] Chetayev | Theoretical Mechanics[END_REF] for equilibrium, i.e., it expresses the principle called the virtual work principle: for the equilibrium of the system the sum of elementary work of active forces on any virtual displacement of the system is either zero or negative; it is zero for non-loosening displacements and negative for loosening displacements.

This principle for the case of bilateral constraints was established by Bernoulli, and for the case of unilateral constraints by Fourier. [START_REF] Fourier | Mémoire sur la statique, contenant la démonstration du principe des vitesses virtuelles et la théorie des moments[END_REF] In the case of potential forces, when

where U(r 1 , . . ., r N ) is the force function, inequality (2.1) takes the form (2.2)

The d'Alembert-Lagrange principle

Let us consider the motion of a system subject to ideal constraints (1.1) and (1.2). Taking into account expressions (1.5) for the forces of the strained constraints, we can represent the equations of motion in the form The sum of the elementary work of active forces F and inertia forces -m r on any virtual displacement is equal to zero or is non-positive, depending on whether all constraints are bilateral or whether there are unilateral constraints among them; the equality sign occurs in the case of non-loosening virtual displacements, and the less than sign in the case of loosening displacements.

Principle (3.3) was derived as a consequence of the equations of motion (3.1). It is possible, conversely, to adopt principle (3.3) as the primary one and, from it, then to obtain the equations of motion (3.1).

Relation (3.3) defines the relationship between the virtual displacements ␦r that are allowed by constraints (1.1) and (1.2) imposed, the active forces F applied to the system and the accelerations r caused by the prescribed forces under the constraints imposed.

When all the accelerations r = 0 ( = 1, . . ., N) and, consequently, the system is in equilibrium, the d'Alembert-Lagrange principle (3.3) becomes the virtual work principle (2.1) (the Fourier inequality).

The quantity P = Fm r is sometimes called the lost force. In terms of this quantity, the d'Alembert-Lagrange principle can be represented by an inequality similar to that of (2.1)

(3.4)
In the case of active potential forces, the d'Alembert-Lagrange principle takes the form (3.5) where U = U(t, r 1 , . . ., r N ).

Principle (3.3) was first obtained for bilateral constraints by Lagrange 6 by combining the d'Alembert and Bernoulli principles, and for unilateral constraints by Ostrogradskii. 7

The Gauss principle of least constraint 8

In the d'Alembert-Lagrange principle, the simultaneous (␦t = 0) positions r and r of the system in the actual and any infinitely close kinematically possible motions, differing by quantities ␦r = rr , satisfying conditions (1.3), are compared. Here, principle (3.3) does not have the nature of the principle of stationarity (extremality) of any function.

Gauss managed to convert the d'Alembert-Lagrange principle into the principle of least constraint by narrowing the set of motions compared: the actual motion is compared with kinematically possible motions in which, at a given instant of time t, the radius vectors r and velocities ṙ of the points M of the system have the same values as in the actual motion. Gauss called these motions conceivable -for them, the accelerations of points ṙ that are possible with the prescribed configuration and velocity of the system are cosidered.

At the instant of time t + , where is a small quantity, we have Neglecting terms of the third and higher orders in , with Gaussian variation, at the instant of time t + dt we obtain Gauss called the quantity Z [Eq. (4.

2)] a constraint and adopted it as a measure of the deviation of the actual motion from the motion of a system free from all constraints. Consequently, for the actual motion, constraint Z has the least value in the class of conceivable motions. This statement forms the Gauss principle of least constraint. For holonomic and linear non-holonomic constraints, the Gauss principle has the same generality as the d'Alembert-Lagrange principle. 2 [START_REF] Chetayev | A modification of Gauss principle[END_REF] Chetayev proposed an interesting modification of Gauss principle. Along with the actual motion of the system, we will examine the class of Gauss conceivable motions of the system () in the time interval from t to t + dt, determined at the instant of time t by the values of r and ṙ = v , and by constraints (1.1) and (1.2) imposed on the system.

Modification of Gauss principle: the Chetayev principle

The work of the prescribed forces F on an infinitesimal displacement of conceivable motion is equal to

where ␦ṙ is the variation of velocity of a point M within a time interval dt for the conceivable motion. From this expression we will subtract the expression for the work on the examined displacement of the forces m r that would be sufficient to create the actual motion if the points M were completely free:

As a result we obtain the expression (5.1) for the work on an elementary cycle consisting of forward conceivable motion in a field of prescribed forces and backward (reverse) motion in a field of forces that would be sufficient to create the actual motion if the points of the system were completely free.

The work on a similar cycle constructed for the actual motion is equal to We will denote by the operation of calculating the increment due to transition from the actual motion to a close conceivable motion. Since the forces F are independent of the accelerations, it follows that F = 0, and Eq. ( 5.3) can be represented in the form (5.4) According to the relation (4.1), we obtain Applying to this inequality the operation , we find (5.5) Consequently, the Chetayev principle of maximum work holds: the work A on an elementary cycle consisting of forward motion in the field of prescribed forces and reverse motion in the field of forces that would be sufficient to create the actual motion if the mechanical system were completely free has a maximum for the actual motion in the class of Gauss conceivable motions.

The Jourdain principle 10

Unlike the d'Alembert-Lagrange, Gauss and Chetayev principles, the actual motion in the Jourdain principle is compared with those kinematically possible motions for which at a given instant of time t the radius vectors r of the points of the system are the same as in the actual motion.

With Jourdain variation, at the instant in time t + dt we will have

The Hamilton-Ostrogradskii principle

Let us consider the motion of a holonomic mechanical system in the time interval from t 0 to t 1 , and also, for each instant t, the virtual displacements ␦r from position r occupied in the actual motion. We will that assume the virtual displacements ␦r are functions of t that belong to class C 2 and vanish at instants of time t 0 and t 1 corresponding to possible positions A 0 and A 1 of the system in the configuration space. The sequence of positions r + ␦r can be regarded as a varied or tortuous path that for a non-holonomic system does not, however, satisfy the constraint equations. [START_REF] Suslov | Theoretical Mechanics[END_REF] The variation of the kinetic energy of the system is given by Integrating with respect to t, we will have

The first term on the right-hand side of this equation is equal to zero, since, according to the condition, ␦r = 0 at instants of time t 0 and t 1 , which implies

The right-hand side of this equation is non-positive according to the d'Alembert-Lagrange principle (3.3), so that (7.1)

The Hamilton-Ostrogradskii principle states that the integral with respect to time of the sum of the variation of the kinetic energy and the work of the active forces on virtual displacements ␦r for the actual motion is equal to zero or is negative in comparison with motions along tortuous paths between the same positions A 0 and A 1 if the motion of the system along all paths is accomplished in the same time interval t 1t 0 , and here t 0 and t 1 are the initial and final instants of time that are common to all paths.

In the case of active potential forces, when F = grad r U ( = 1, . . ., N), we obtain the Hamilton principle

(7.2)
where L = T + V is Lagrangian function, i.e., for the actual motion (7.3) Principle (7.1) was originally stated by Hamilton 11 for the case of stationary holonomic systems with bilateral constraints under the action of potential forces, and was then extended by Ostrogradskii [START_REF] Ostrogradsky | Sur les integrales des équations générales de la dynamique[END_REF] to the non-stationary case, and also to the case of non-potential forces. For holonomic systems with unilateral constraints, principles (7.1) and (7.3) were investigated by Stavrakova. [START_REF] Nye | The Hamilton-Ostrogradskii principle for a system with unilateral constraints[END_REF]

The principle of least action in Lagrangian form

The first verbal formulation of the principle of least action was given by Maupertuis (1744) who then (1746) declared this principle as the universal law of motion and rest. [START_REF] Maupertuis | Les lois de mouvement et du repos déduites d'un principe métaphysique[END_REF] : "General principle. When in nature a certain change occurs, the amount of action necessary for this change is the least possible. The amount of action is the product of the mass of the bodies by their velocity and the distance that they travel".

The mathematical formulation of the principle of least action for the case of a point mass was given by Euler. [START_REF] Euler | Methodus inveniendi lineas curvas maximi minimive proprietate gaudentes sive solutio problematis isopermetrici latissmo sensu accepti[END_REF] Lagrange 6 derived the principle for a system of point masses with bilateral constraints.

Below we will assume that the holonomic constraints imposed on the system are stationary, impacts possible in the system as it comes against unilateral constraints or leave them being absolutely elastic (coefficient of restitution ␦ = 1) Ref. 1 and the applied active forces are potential with a force function U(r 1 , . . ., r N ) that is independent of time. Under these conditions the following energy integral exists (8.1) In fact, under these conditions, the actual displacements dr = ṙ dt ( = 1, . . ., N) are among the virtual displacements. Multiplying Eq. (3.1) scalarly by dr , and summing over all , we obtain, with the above assumptions, the equation from which integral (8.1) follows. If the relations between the Cartesian and generalized coordinates are timeindependent, the Lagrangian function L(q, q) = T + V will also be independent of time, and integral (8.1) will exist in this case also.

Left to itself, the mechanical system may choose its motions from motions with a given store of the total energy h, and therefore the set of motions being compared may be limited by condition (8.1). [START_REF] Chetayev | Theoretical Mechanics[END_REF] The actual motion will vary so that the initial position A 0 and the final position A 1 of the system remain unchanged. We will obtain the varied trajectory by producing at each instant of time t the virtual displacement ␦r (␦q i ) relative to the actual trajectory, where the instant of time t + ␦t corresponds to the position r + ␦r (q i + ␦q i ) on a trajectory being compared. Since the total energy H should be equal to h on all the trajectories considered, from this the velocity at the position r + ␦r is found. Thus, unlike earlier, the time is varied, and the duration of the motions being compared will differ from the duration of the actual motion between points A 0 and A 1 . Below we will consider ␦r (␦q i ) and ␦t to be functions of time t of class C 2 . [START_REF] Hölder | Über die Prinzipien von Hamilton und Maupertuis[END_REF] We will integrate expression (3.5) of the d'Alembert-Lagrange principle with respect to t in the range from t 0 to t -the instant the final position A 1 is reached. It expresses the principle of least action in Lagrangian form for ideal smooth unilateral constraints. The actual motion, being compared with other close motions occurring with the same constant energy H, satisfies condition (8.2); here, the variations of positions should be virtual displacements, and the initial and final positions of the system remain unvaried.

The principle of least action in Jacobian form

Jacobi [START_REF] Jacobi | Vorlesungen über Dynamik[END_REF] iminated the time t from the expression of action (8.2) using the energy integral and gave the principle of least action a geometric form.

Consider the configuration space defined by the line element [START_REF] Chetayev | Theoretical Mechanics[END_REF] and let us represent the kinetic energy of the system in the form

(9.1)
where is a parameter. Let the family of surfaces (q 1 , . . ., q n ) = c intersect the actual trajectory of the system, and also infinitesimally close trajectories drawn through the points A 0 and A 1 of configuration space. Each of these curves can be considered to be specified by its coordinates represented as functions of . Let the virtual displacements of the system from some point of the actual trajectory to a point of a being compared trajectory correspond to the same value of , and let the positions A 0 and A 1 correspond to the values of 0 and 1 respectively. We will express from Eq. (3.1) and replace T with U + h; we then obtain Thus, Lagrange variational principle with fixed lower and free upper bounds is transformed to Jacobi variational principle with fixed bounds 0 and 1 and fixed ends A 0 and A 1 .

The Euler-Lagrange equations for problem (9.3) have the form 2 (9.4) (q i = dq i /dτ). These are the differential equations of the actual trajectory in the configuration space. The motion in time is found by integrating Eq. (9.2), which gives t as a function of the parameter .

If dt = 0 on the trajectory considered, then, as the parameter , it is possible to select the time t measured on the actual trajectory, for which integral (8.1) occurs. Here, the Eq. (9.4) take the form of Lagrange's equations of motion of the second kind which also proves the Jacobi principle.

Note 2 that in mechanical problems it is not possible to fix the values of t 0 , t 1 , A 0 , A 1 and h in advance.

The Suslov-Voronets principle for non-holonomic systems

Consider a mechanical system with the generalized coordinates q ( = 1, . . ., r + b) which identically satisfy all holonomic constraints. Suppose, furthermore, that the system is subjected to b non-holonomic constraints 1 (10.1) so that the virtual displacements satisfy the conditions The symbol ␦B ␤ denotes expressions that are linear and homogeneous in ␦q ␤ . If Eq. (10.1) are integrable, then all ␦B ␤ ≡ 0.

We integrate the expression of the d'Alembert-Lagrange principle in generalized coordinates from t 0 to t 1 (10.4) assuming that, for t = t 0 and t = t 1 , all independent variations ␦q ( = 1, . . ., r) and consequently also dependent variations ␦q r+b are equal to zero (␤ = 1, . . ., b).

We will have and also, taking relations (10.3) into account Substituting these results into relation (10.4), we obtain the Suslov-Voronets principle [START_REF] Suslov | Theoretical Mechanics[END_REF][START_REF] Voronets | On the equations of motion for non-holonomic systems[END_REF] (10.5)

(1. 1 )

 1 and unilateral if they are expressed by inequalities of the form(1.2)

(3. 1 )

 1 Multiplying each of the equations in system (3.1) scalarly by ␦r , and summing over all = 1, . . ., N, we obtain (3.2) from which, taking into account inequality (1.4), we obtain an analytical expression for the d'Alembert-Lagrange principle(3.3)

Substituting

  ␦r into expression (3.3) of the d'Alembert-Lagrange principle, we will have (4.1) or, taking into account that the prescribed active forces F = F (t, r, ṙ) are independent of the accelerations r , and introducing into cosideration the function (4.2) we will rewrite relation (4.1) in the form (4.3) from which ␦ 2 Z = m (␦r ) 2 > 0.

(6. 1 )

 1 Substituting this expression into expression (3.3) of the d'Alembert-Lagrange principle, we obtain the Jourdain principle (6.2)

  We obtain Since the previous relation takes the form Varying the energy integral (8.1), we obtain ␦T = ␦U, and consequently the last relation can be represented in the form (8.2)

(10. 2 )

 2 The operations of differentiation and variation are permutable only for = 1, . . ., r, i.e. and for the remaining b dependent velocities, from inequalities (10.1) and (10.2) we derive the relations(10.3)

Acknowledgement

I wish to thank V. F. Zhuravlev for drawing my attention to the monograph Ref. 19.


Jonàs Martínez

Jérémie Dumas

Sylvain Lefebvre

Procedural Voronoi Foams for Additive Manufacturing

Keywords: 3D printing, additive manufacturing, procedural modeling, material design Concepts:, Computing methodologies → Shape modeling;

rigid flexible Figure 1: Our method procedurally generates structures with graded material elasticities, which can be directly fabricated. Here the user paints elasticity on a 3D model to create a flexible figurine. Model: Moomin (thing:1173447) by Jeroentjj.

Introduction

Additive manufacturing enables the fabrication of objects having unprecedented complexity. This capability is often understood in terms of the overall shape of objects. However, it is also possible to fabricate parts filled with microstructures -having intricate internal details in the order of tens of microns. The macro-scale mechanical properties of the object are then directly influenced by the geometry of the microstructures. In particular, careful design of the structures affords for parts that are lighter while remaining rigid enough for their intended use. This reduces material usage, shipping and transportation costs. In addition, inner structures can progressively vary within the object and adapt to varying rigidity requirements between regions subject to different stresses.

There are several challenges to achieve these goals. First, the fine scale geometry of the structures has to produce the desired large scale elastic behavior. This often implies formulating challenging global optimizations, either to directly synthesize the fine scale geometry or to fill the shape with precomputed microstructures (see Section 2). Second, the structures are very small compared to the size of the objects (e.g. tens of microns versus tens of centimeters). Therefore, the meshes describing the models become quickly prohibitively large, posing important computational challenges for simulation, visualization, and fabrication. Third, the structures have to enforce fabricability constraints. The main industrial processes have different requirements depending on whether they locally deposit material (e.g. fused filament fabrication, resin droplets) or whether they locally solidify a bed of material itself acting as a support (e.g. selective laser sintering). In the first case the structures should not present any disconnected parts during fabrication, while in the second case they should not enclose non-solidified material in pockets.

To answer these challenges we draw inspiration from proce-dural noises in computer graphics, where an infinite amount of content is produced at low, constant memory cost while precisely controlling the statistical properties of the produced noises [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF]. This hints at the possibility of generating procedural, stochastic microstructures that directly exhibit the desired elastic behavior, without further optimization. Our approach explores this idea and defines a procedure to synthesize open-cell foams that enforce fabrication requirements while having precisely controlled elastic properties. The foam parameters can vary spatially to follow gradients of elasticity.

As the microstructures are procedural we only generate their details when needed for fabrication. Typically, evaluation happens at the slice level, just before sintering or curing a layer of material. The microstructures are stochastic and aperiodic in nature. Stochasticity results in an exceptionally good isotropic behavior, and lets us grade the properties without introducing discontinuities along pre-defined boundaries. Aperiodicity removes the need for a global optimization when conforming the structures to a surface.

Our contributions are:

• The definition of procedural Voronoi foams that can be evaluated very efficiently, have precisely controlled isotropic elastic behavior, and can be spatially graded to produce gradients of elasticity.

• A methodology to derive an inverse mapping, from a target elasticity to the parameters driving the microstructure generation.

• A complete implementation that maps well to stackless, massively parallel architectures.

• Applications to the 3D fabrication of objects filled with procedural Voronoi foams, with proper handling of the outer object frame.

The result is an algorithm that generates microstructures with a prescribed elastic behavior on the fly, during fabrication. No optimization process is required to adapt to a different object or to match a graded elasticity field. There is no limit to the size of the printable objects as the microstructures never fully reside in memory. Therefore, our approach will naturally scale with future technology as the resolution and size of printable object increases.

Scope of the Paper. Isotropic elastic materials are described by two parameters: Young's modulus and Poisson's ratio. Intuitively the Young's modulus captures how rigid or soft an object is, while the Poisson's ratio captures how one dimension stretches with another. In this work we focus on varying the Young's modulus while preserving the Poisson's ratio of the base material.

For spatially varying elasticity we assume the target scalar field is given as input and do not consider its computation.

Related Work

Modeling for Fabrication. In the recent years there has been a strong interest in techniques helping the user to design physical objects with controlled behaviors, for instance to avoid fragile regions and reinforce them [START_REF] Umetani | Cross-sectional structural analysis for 3D printing optimization[END_REF][START_REF] Zhou | Worst-case structural analysis[END_REF][START_REF] Stava | Stress relief: Improving structural strength of 3D printable objects[END_REF], to allow the user to balance shapes [Prévost et al. 2013;[START_REF] Bächer | Spin-it: Optimizing moment of inertia for spinnable objects[END_REF], or to produce strong yet lightweight inner structures [START_REF] Wang | Costeffective printing of 3D objects with skin-frame structures[END_REF][START_REF] Lu | Build-to-last: Strength to weight 3D printed objects[END_REF]]. These latter approaches are related to our work in that they seek to produce structures of maximum rigidity under a prescribed volume to strengthen 3D printed objects. However, these approaches rely on global optimization of as rigid as possible macro-scale structures, while we focus on the procedural generation of micro-structures with controlled elastic properties.

Graded Elasticity. In the field of mechanical engineering [START_REF] Sigmund | Tailoring materials with prescribed elastic properties[END_REF][START_REF] Sigmund | Design of smart composite materials using topology optimization[END_REF] introduced the optimization of micro-scale structures to achieve specific macro-scale behaviors, such as controlled elasticity. A key idea is to consider the behavior of a composite material, made of infinitely many repetitions of a base material tile. This limit behavior is captured by the theory of homogenization for linear elasticity [START_REF] Allaire | Shape optimization by the homogenization method[END_REF]] which relates the unit material tile to the elastic properties of its (infinite) periodic tiling. We provide more details in the supplemental material. Therefore, most techniques pose an inverse homogenization problem, optimizing for a tile producing a target elastic behavior [START_REF] Sigmund | Materials with prescribed constitutive parameters: An inverse homogenization problem[END_REF][START_REF] Zhou | Design of graded two-phase microstructures for tailored elasticity gradients[END_REF][START_REF] Radman | Topology optimization of functionally graded cellular materials[END_REF][START_REF] Xia | Design of materials using topology optimization and energy-based homogenization approach in Matlab[END_REF]. [START_REF] Andreassen | Design of manufacturable 3D extremal elastic microstructure[END_REF] include fabrication constraints in the structure optimization.

In the field of computer graphics the design of materials for fabrication has become an important direction of research, as we strive to enable artists and designers to physically realize virtual models of deformable objects [START_REF] Skouras | Computational design of actuated deformable characters[END_REF][START_REF] Xu | Interactive material design using model reduction[END_REF][START_REF] Pérez | Design and fabrication of flexible rod meshes[END_REF]. [START_REF] Bickel | Design and fabrication of materials with desired deformation behavior[END_REF] proposed a data-driven approach to design materials with prescribed deformation. Base materials are stacked by an optimizer to obtain the target properties. [START_REF] Schumacher | Microstructures to control elasticity in 3D printing[END_REF] extend this idea in two ways. First, elementary material tiles are optimized by homogenization to cover a large spectrum of elastic behaviors. Second, a process globally optimizes for a choice of tiles in a grid covering the object, to achieve the desired, spatially varying elastic behavior. The process considers fabrication constraints and connectivity between adjacent tiles. [START_REF] Panetta | Elastic textures for additive fabrication[END_REF] take a different approach, by optimizing for an optimal family of elementary tiles among a large -but restricted -set of possibilities. The tiles take into account fabrication and connectivity constraints. Once the best family is determined, it is possible to arrange the tiles in a grid such as to obtain a spatially varying elasticity.

These works achieve impressive results at relatively low computational costs, thanks to the regularity of the grid and the periodicity assumptions: homogenization is performed on a single base tile instead of having to rely on a global optimization, and the periodicity affords for compact descriptions [START_REF] Pasko | Procedural function-based modelling of volumetric microstructures[END_REF]]. However, the tile-based approach has a number of drawbacks. First, grading the material by changing the tiles within the grid requires careful handling of the interface between neighboring tiles. Either the tiles are pre-optimized to have matching boundaries [START_REF] Panetta | Elastic textures for additive fabrication[END_REF], or a global selection process is required to strike a compromise between continuity across neighbors, grid discretization and elasticity objectives [START_REF] Schumacher | Microstructures to control elasticity in 3D printing[END_REF].

The global optimization does not scale well with the size of the structures and objects -especially as it is desirable to produce the smallest possible tiles to converge towards the limit behavior computed by homogenization. Our approach solves this thanks to a direct and simple relationship between elastic properties and geometric parameters. Second, it is often desirable to conform the elasticity field to the object surface, for instance having a more rigid crust within a dis-tance of the surface. The global periodic nature of the grid makes this difficult. One possible way is to optimize for a 3D parameterization of the grid within the object. While this is a topic of intense research (related to hexahedral meshing), this remains difficult [START_REF] Staten | Why is hex meshing so hard? Presentation at Sandia National Laboratories[END_REF]], and the effect on the final elastic properties is hard to precisely establish. Instead, we seek to produce stochastic, aperiodic structures which are by nature simpler to conform to the gradients of a field as they do not require specific spatial alignments. Finally, the isotropy of periodic tilings is equal to that of the optimized base tile, while the isotropy of stochastic foams further improves with larger extents of foams.

Efficient Microstructure Generation. [START_REF] Pasko | Procedural function-based modelling of volumetric microstructures[END_REF] considered procedural definitions of periodic microstructures.

The parameters of the microstructures can vary spatially to produce graded materials [START_REF] Fryazinov | Multi-scale space-variant FRep cellular structures[END_REF], for instance to reinforce an object following a cross-sectional stress analysis [START_REF] Li | Interior structural optimization based on the density-variable shape modeling of 3D printed objects[END_REF]. [START_REF] Brennan-Craddock | The investigation of a method to generate conformal lattice structures for additive manufacturing[END_REF] computes the intersection of the object and periodic microstructures on a per-slice basis. A frame structure is built on the object surface by subtracting the microstructure cells from a thick surface shell. Our frame structure follow a similar intuition for the case of Voronoi foams. OpenFAB [START_REF] Vidimče | OpenFab: A programmable pipeline for multi-material fabrication[END_REF]] provides a specialized language to describe procedural microstructures. The geometric details are efficiently evaluated at slicing time, streaming voxels to the printer. Similar to these approaches, we evaluate the microstructures procedurally during slicing. However, the microstructures we generate are aperiodic graded foams whose Young's modulus is directly and precisely controlled.

Recent works consider the problem of generating nonperiodic microstructures with varying density. The adaptive voids approach [Medeiros e [START_REF] Medeiros E Sá | Adaptive voids[END_REF]] relies on a subdivision scheme to produce denser structures near the surface of an object. [START_REF] Brackett | An error diffusion based method to generate functionally graded cellular structures[END_REF] perform a sequential dithering of a density field to keep a subset of points, which are then used to define an open-cell foam. While closely related to our work, these approaches do not explicitly control the Young's modulus of the produced structures, nor afford for an efficient parallel evaluation.

Recent software for additive manufacturing propose microstructure generation packages. In particular, Within-Labs [Autodesk 2016] proposes trabecular structures resembling Voronoi foams. While the parameters can be varied, to the best of our knowledge there is no direct control of the Young's modulus and the structures are not defined by a procedure akin to procedural solid textures.

Open-Cell Foams. Our procedural microstructures belong to a specific class of microstructures known as open-cell foams. These structures occur naturally and can be obtained from physical processes in a variety of materials including metals. Therefore, in the field of mechanical engineering there has been a strong interest in modeling and analyzing the properties of these structures.

Interestingly, naturally occurring open-cell foams are often idealized as edges of Voronoi cells [START_REF] Gibson | Cellular solids: structure and properties[END_REF]].

As we propose to generate foams from procedurally generated Voronoi diagrams, these works are highly relevant for the study of the mechanical behavior of our structures. In particular, from these studies we can expect the following properties: 1) the elastic behavior of open-cell irregular foams is highly isotropic [START_REF] Luxner | Numerical simulations of 3D open cell structures -influence of structural irregularities on elasto-plasticity and deformation localization[END_REF], and 2) their

Young's modulus relate almost linearly to their geometric parameters (thickness, density) [START_REF] Gibson | Cellular solids: structure and properties[END_REF][START_REF] Roberts | Elastic properties of model random three-dimensional open-cell solids[END_REF]. This is an exceptionally good property for our purpose of producing graded elastic materials, as the simple relationship affords for a direct derivation of structure parameters given an elasticity target. We verify that our structures meet these observations in Section 4.3.

Prior studies also indicate that the Poisson's ratio of opencell foams does not significantly vary [START_REF] Gibson | Cellular solids: structure and properties[END_REF], and our structures indeed share this limitation.

Procedural Voronoi Foam Generation

We now introduce our approach for the procedural generation of open-cell foams. We describe the procedural synthesis of the structures in Section 3.1 and discuss implementation in Section 3.2. We explain how to derive the parameters of the structures from a desired target elastic behavior in Section 4 and present results and applications in Section 5.

Procedural Generation

We seek to define aperiodic procedural structures akin to procedural textures in computer graphics. The structure is defined as a function F : R 3 → {0, 1} which is evaluated during display and slicing at every point in space, at the desired resolution.

To have the computational advantages of procedural textures, F has to follow a number of requirements [START_REF] Lagae | A survey of procedural noise functions[END_REF]], that we summarize as follows: 1) F has to evaluate in constant time and constant memory regardless of the point of evaluation.

2) The size of F -its program and built-in data -has to be independent from the size of the generated content. This is the case of our technique which produces arbitrary large aperiodic content from a constant, small memory footprint (a few hundreds of bytes).

In addition, the structure has to enforce geometric requirements to be printable. First, there should be a minimal number of pockets (holes) enclosing printing material. This advocates for an open-cell structure made of beams along the edges of a cellular structure. Second, there should be no disconnected parts appearing during fabrication. It is easy to see that convex cells enforce this property everywhere but at the boundary (which we discuss in Section 5.3). The convex cells of Voronoi diagrams are therefore well suited.

Our procedural Voronoi foams are defined by two parameters: the density ρ of Voronoi seeds per unit volume (seeds/mm 3 ), and the radius τ of the beams along the edges (mm). Density may vary spatially and is given as a function, i.e. at a given point x the desired density is ρ(x). We assume the variations to be smooth compared to the size of the structure cells.

Open-Cell Voronoi Foams. We seek to design a procedural function Fρ,τ that produces beams of thickness 2τ along the edges of a Voronoi cell structure having a density ρ. In addition, we would like to allow for ρ to be spatially varying.

Our procedural generation is inspired by the seminal work on cellular solid textures by [START_REF] Worley | A cellular texturing basis function[END_REF], revised to produce an open-cell structure. We further extend the procedural scheme to afford for spatially varying densities.

Worley defines procedural cellular textures by using pseudorandom sequences to generate seeds in a virtual grid, fol- lowing a Poisson distribution (not to be confused with a Poisson disc distribution). Given an evaluation point q the seed closest to q is determined. The seed id is used to derive a color value at q, for instance coloring each Voronoi cell differently. The set of grid cells that can have an influence is limited to cells neighboring q. This stems from the fact that each cell contains at least one seed -a deviation from a pure Poisson distribution to enable efficiency [START_REF] Worley | A cellular texturing basis function[END_REF]] -and, therefore, the closest seed is necessarily within the 2-ring of neighboring cells. This leads to the constant time evaluation property, as the number of considered seed points remains below a constant everywhere in space.

Our algorithm achieves similar properties while generating the edges of a Voronoi open-cell foam. Given an evaluation point q our goal is to write a function returning 0 (empty) if q is not inside a beam of the structure and 1 (solid) otherwise. The pseudo-code of our structure generation is in Algorithm 1, and a graphical overview is in Figure 2. It starts by generating all the seeds that could contribute to the definition of the Voronoi cell by calling algorithm GatherSeeds (line 1) -we detail this algorithm later. We next enumerate all the line equations of the edges (line 4) and for each compute pl, the point closest to q on the line (line 5). If the distance between q and the edge is greater than the beam radius τ (line 6), q is not influenced by this edge. Otherwise, q might be inside a beam. To be certain that it is the case, we have to verify that pl is indeed on an edge of the Voronoi cell containing q. A counter example is illustrated in Figure 3. We verify this by considering whether pl belongs to the Voronoi cell of a neighboring seed (lines 8-11). If that is the case, then we ignore this line equation (line 10). Otherwise, q is in a solid region and we return 1 (solid) line 13.

We next discuss how the seeds are gathered and generated.

Gathering Seeds. Algorithm 1 requires all seeds that can influence the result at q. If some required seeds were to be missed, the produced structure could fail to print or break. However, we only need to be conservative: as long as we have a superset of the required seeds, the algorithm will produce the correct result.

We generate seeds in a grid, and we guarantee that all grid cells receive at least one seed. This bounds the number of grid cells that we have to consider. As explained in Figure 4, the Voronoi cell of a seed cannot be influenced beyond a 2-ring of neighbors.

Algorithm 1: EvalStructure Fρ,τ (q) Input: Density field ρ, beam radius τ , query point q.

Output: Algorithm 2 gives the pseudo-code for gathering the seeds around q. The evaluation point q might belong to the Voronoi cell of any of the seeds within (at most) a 2-ring radius. We therefore first search for the seed closest to q, and then gather seeds with a 2-ring around the closest seed. SubdivideCell produces at least one seed per grid cell -possibly more with spatially varying density. We detail this algorithm next.

Voxel state ∈ {0, 1} 1 N ← GatherSeeds(q) ; // seeds influencing q 2 for i ← 1 to |N | do 3 for j ← i + 1 to |N | do 4 bl ← BisectorLineEquation(N [0], N [i], N [j]); 5 pl ← ClosestP ointOnLine(q, bl); 6 if q -pl τ then 7 accept ← true; 8 for k ← 1 to |N | do 9 if N [0] -pl > N [k] -
Algorithm 2: GatherSeeds(ρ, q) Input: Density field ρ, query point q.

Output: Set of seeds possibly influencing q Seed Generation and Density Control. We now describe SubdivideCell. It generates seeds in each cell, at least one and possibly more by subdividing to locally adapt to density. Conceptually our technique is based on a primal subdivision of the coarsest density grid, where each parent cell is split regularly in eight children. It is important to recall, however, that the grids are never stored: all computations happen implicitly and on the fly. We refer to the first level of grid cells as the coarse grid cells (between 2mm and 5mm in our implementation). Another important design goal of our subdivision process is to avoid bias: The statistics of the point distributions remain constant -up to a scaling factor -at all density levels.

1 N ← ∅ ; 2 V isited ← ∅ ; 3 cq ← GridCellEnclosing(q) ; 4 closest ← ∞ ; 5 for cell ∈ T woRingN eighborhood(cq) do 6 V isited ← V isited ∪ {cell};
The pseudo-code is given in Algorithm 3 for the general case of a spatially varying density field. Given a grid cell of size l and center c, we compute the number of seeds it has to contain as l 3 × ρ(c) (line 2, cell volume times density). We evaluate the density at the cell center, but more elaborate schemes could be used (e.g. multi-point evaluation, or preintegration of ρ in a summed area table). The field ρ is clamped to a minimum value to ensure that the coarsest cells always receive at least one seed. When the target number of seeds in a cell is above 2 3 -there is more than exactly one seed per subdivision child of the current cell -we recurse and subdivide the cell (line 13). Otherwise, we randomly select n = l 3 × ρ distinct children and draw exactly one seed in each (lines 4-7). We take into account the remaining fraction f = t -n by drawing an additional sample in a next child cell, with probability f (lines 8-10). Note that all random number generators are pseudo-random sequences seeded by the grid cell coordinate. Therefore, for a same initial grid cell, the exact same set of seeds is produced.

The pseudo-code in Algorithm 3 cannot be implemented directly on massively parallel architectures (GPUs) due to the recursive calls. We present in the supplemental material a stackless iterative version that maps well to massively parallel processors. ). While this value may also vary spatially (τ (q)) we only vary density in our approach (see Section 4.2).

We next discuss how to select the parameters to reach a target elastic behavior.

Implementation

We implement our technique in an image-based slicer, which either directly sends images to the printer (SLA) or extracts contours (SLS). We slice at a resolution between 10µm and 50µm per pixel.

The procedural foam is implemented as an OpenCL kernel processing all pixels of a slice in parallel on a GPU. In addition, we implement a supersampling procedure to obtain gray-scale exposure levels on SLA processes. This is done by performing supersampling around the evaluation point q in a small loop around lines 6-13 of Algorithm 1, counting the number of times q + lies within the beam.

Exact beam contours could be extracted by locally constructing the explicit geometry of the Voronoi cells around q in the manner of Algorithm 1. We found this approach less convenient than our implicit description that fits very well existing massively parallel architectures, and avoids having to union beam geometries explicitly.

Foams With Controlled Elasticity

We now consider the problem of selecting the microstructure parameters (ρ, τ ) to achieve a target elastic behavior. This is done by statistical analysis of the homogenized elasticity tensor of Fρ,τ , for different parameters.

Analysis Using Homogenization

For a choice of (ρ, τ ) we produce a periodic version of the structure in a base volume, and apply homogenization to compute the elastic tensor of the corresponding periodic media. By doing this for many choices of (ρ, τ ) we reconstruct the underlying relationship between ρ, τ and the elastic properties (Young's modulus, Poisson's ratio). The graph resulting from this analysis is shown in Figure 5. It is obtained for a material with a unit Young's modulus (= 1) and a Poisson's ratio v = 0.3 which matches most plastics. Only changing the Poisson's ratio requires re-running the simulation as results scale linearly with Young's modulus.

We perform homogenization on a grid of hexahedral elements (see supplemental material for details). The size of the elements is chosen to properly capture the beams; our experiments showed that using half the diameter gives stable results. The stiffness of each element is computed by calling our implicit procedural function with super-sampling, which returns the volume of structure intersected by the element.

Homogenization of Aperiodic Foams. With this approach we are making an important assumption. Our structures are not periodic: they form an aperiodic foam of constant density in space. We therefore assume that the overall behavior of the aperiodic foam is similar to the periodic behavior of a sufficiently large base volume. To determine the size of the base volume we rely on the expected isotropic behavior of irregular open-cell foams [START_REF] Luxner | Numerical simulations of 3D open cell structures -influence of structural irregularities on elasto-plasticity and deformation localization[END_REF]]. We perform homogenization for volumes of increasing spatial extent and consider the deviation of the computed tensor from a perfectly isotropic tensor.

The homogenized tensor is a full symmetric tensor composed of 21 independent variables. We approximate the full tensor with the tensor of an isotropic material and consider the residual error. The tensor of an isotropic material is expressed from only two independent variables -Young's modulus and Poisson's ratio -and we denote it as C I (E, v). Its expression is given in the supplemental material.

Similarly to [START_REF] Schumacher | Microstructures to control elasticity in 3D printing[END_REF] we approximate the homogenized elasticity tensor C H by minimizing: Each dot is the result of homogenization for a given (ρ, τ ). For each (ρ, τ ) we generate multiple instances of foams using different random seeds. As can be seen the dots tightly cluster, indicating that randomness has little impact on the elastic properties. The maximum deviation after the fitting is 3.3% of the Young's modulus. The data points are used to fit a polynomial function (see Section 4). For clarity, the graph shows a subset of our data that extends to higher densities for lower beam radii.

ξ(C H ) = min E,v C I (E, v) -C H 2 F -1 v 0.5 0 E EM (1)

divergence from isotropy). As the spatial extent of the simulated volume grows the tensor tends towards the ideal isotropic tensor.

, where EM is the upper bound Young's modulus of the solid base material. The expression of C I is nonlinear with respect to E, v and we therefore solve this minimization with a nonlinear optimizer [START_REF] Johnson | The NLopt nonlinear-optimization package[END_REF]].

The approximation error measures whether homogenization was able to properly capture the isotropic behavior of the structure. Figure 6 reports the error in isotropy for various grid sizes. In practice, we find that a volume having at least 60 × 60 × 60 hexahedral elements provides a stable result across a wide range of beam radii.

Collecting Data

Points. In a one-time precomputation we generate a large number of pairs (ρ, τ ). We uniformly sample ρ and τ to cover the range that can be fabricated (see Section 4.2). For each selection of (ρ, τ ), we generate a number of different structure realizations by varying the global random seed. We compute the Young's modulus and Poisson's ratio obtained by applying homogenization on each unit volume. This gives us a large number of data points characterizing the elastic behavior of the structures, as shown in Figure 5.

We then fit a polynomial function P on the experimental results to correlate the density, beam radius, and Young's modulus. We optimize for the coefficients of a degree 4 polynomial by least square fitting, minimizing i (P(ρi, τi) -Ei) 2 . The resulting polynomial fit is shown in Figure 5.

Deriving Parameters for a Target Elasticity

We now describe how to select ρ and τ to achieve a target Young's modulus. As can be seen in Figure 5, for a given target elasticity multiple choices of ρ and τ are possible: the full isocurve where P(ρ, τ ) = Etarget. Our preferred strategy to select the values of ρ and τ is to obtain structures that are as dense as possible, to remain close to the limit behavior computed by homogenization. The limiting factor is the minimal printable beam radius, which we denote as τmin.

To maximize density, we fix τ = τmin. As densities increases with a fixed radius, we reach a point where the block of matter is full. There is no need to further increase density, which leads to a density upper bound ρmax = 1 (2τ min ) 3 . The resulting curve is extracted from the polynomial shown in Figure 5 by intersecting it with the plane τ = τmin -in practice we use τmin ∈ [0.2, 0.3] mm depending on the printer.

For high densities, just before reaching a full block of matter, small pockets start to appear. This is due to the edges of the Voronoi cells merging together. As measured in Figure 7 the volume captured by pockets never exceeds 4%, and is negligible for relative volumes below 0.7. Pockets have a negligible impact under a relative volume of 0.7 and go up to 4% of the full volume around 0.9 relative volume. Thanks to the quasi-linear relationship between relative volume and Young's modulus, we can expect at most a similar error on the elastic behavior.

Other strategies could be developed. For instance, for aesthetic purposes a user could use thicker beams and a lower density to reveal the structures at the same elasticity.

Elastic Behavior: Properties and Analysis

The collection of data points we produced for the parameter fitting lets us fully characterize our procedural foams. We now compare our results to the open-cell foams literature and verify our microstructures behave as expected.

Elasticity

The linear elastic response of low density (0.04 < ω < 0.5) open-cell foams is described by the following model [START_REF] Gibson | Cellular solids: structure and properties[END_REF]]:

E f Es = C ω f ωs n (2)
, where the subscript f denotes the foam and s denotes a solid block of base material.

ω f ωs defines the relative volume ω. The constants C and n change depending on the foam family [START_REF] Roberts | Elastic properties of model random three-dimensional open-cell solids[END_REF].

Figure 8a shows the relationship between the relative volume and the Young's modulus E of our procedural foams. We fit the model for low density foams on the range 0.04 < ω < 0.5, using C = 0.85 and n = 1.95. As can be observed, our data fits this model very well.

Figure 8b shows the Poisson's ratio of our structures, and as can be seen it remains stable. While different materials will give different values of Poisson's ratio, it always remains stable around the Poisson's ratio of the base material as observed by [START_REF] Gibson | Cellular solids: structure and properties[END_REF] on open-cell foams. The base material has a Poisson's ratio of 0.3. The Poisson's ratio tends to 0.5 for low densities, and to the base material Poisson's ratio for high densities.

Randomization

Our procedural foams are based on a (pseudo) random process. Therefore, it is important to consider whether different realization using the same base parameters ρ, τ produce a consistent Young's modulus. Figure 5 reveals this as for each pair (ρ, τ ) we produced three different realizations. Close inspection reveals that the different realizations do not perfectly match. The largest deviation to the fitted polynomial is 3.3% of the Young's modulus and therefore is negligible. This is in agreement with the literature where it was observed that randomization has little influence on the linear elastic properties of open-cell foams [START_REF] Van Der Burg | On the linear elastic properties of regular and random open-cell foam models[END_REF]].

Regularity

Studies of naturally occurring foams often discuss the regularity of the foam, which in our case measures the minimum distance between Voronoi seeds in a bounded domain using the same number of seeds. A Poisson disc distribution provides the most regular cases, while a random point process is the less regular. Our procedural foams lie in-between.

Our seed generation draws at least one sample per grid cell. This is known to be a crude -yet very efficientapproximation of a Poisson disc distribution [START_REF] Cook | Stochastic Sampling in Computer Graphics[END_REF]]. The approximation is less evenly spaced than a high-quality Poisson disc distribution. However, less regular foams exhibit a better isotropy [START_REF] Luxner | Numerical simulations of 3D open cell structures -influence of structural irregularities on elasto-plasticity and deformation localization[END_REF], which is also desirable. Therefore, we believe that our jittered grid approximation is a good compromise between computational efficiency and elastic properties, while allowing for graded densities.

Applications and Results

We now experimentally challenge the behavior of our foams. In Section 5.1 we compare their isotropy to recent tile-based techniques, in Section 5.2 we measure actual printed samples, and in Section 5.3 we describe a complete application to fill 3D models with graded elasticity.

Isotropic Behavior Versus Periodic Tiles

Figure 9 compares the isotropy of the tensors computed by homogenization using tile-based methods. The method of [START_REF] Panetta | Elastic textures for additive fabrication[END_REF] explicitly constructs isotropic tiles. When discretized in our homogenization process they provide close to perfect isotropy; the measured error is due to the limited numerical precision. At a same scale and resolution our structures exhibit lower isotropy; however, performing homogenization with increasingly larger volumes reveals that the residual quickly decreases (see Figure 6). This shows that filling larger volumes with procedural foams improves isotropy -a property that stems from the aperiodic stochastic nature of Voronoi foams. Figure 9 also shows a tile with a less isotropic behavior (manually designed to match [START_REF] Schumacher | Microstructures to control elasticity in 3D printing[END_REF]). As larger volumes are simulated, the isotropy of our structures quickly improves (see also Figure 6).

Experimental Results on Printed Samples

We verify the predicted elastic behavior of our foams on printed samples. We print two families of four samples with varying Young's modulus. Each family uses a different random seed. All samples are printed on a B9Creator V1.2 using red-cherry material. The two families are visible in Figure 10.

We put each sample on a high precision scale and impose to each a same small displacement. We then measure the weight applied to the sample onto the scale (canceling out the structure weight). Due to the limited maximum measurable weight of the scale we consider samples with relatively close Young's modulus. The results are reported in Figure 10. We obtain a good correspondence between the measured data and the prediction; the average deviation of the Young's modulus from linear correlation (dashed line in Figure) is 0.001434, that is less than 0.2% absolute error, and 25% relative error.

For reproducibility please note that the SLA print process, the subsequent cleanup and the material curing all have significant variabilities. We noticed on a few prints (2 out of 16) an incorrect final weight -we discarded these samples. We also noticed that the elastic behavior significantly changes in the first few hours after the print completed. Similarly to [START_REF] Panetta | Elastic textures for additive fabrication[END_REF] we wait for 24 hours before taking measurements.

Procedural Foam with Elasticity Gradients

The process we have described in Section 4 is capable of generating graded structures, with a spatially varying Young's modulus. We exploit this property to generate objects with controlled elastic properties. The user inputs the spatially varying Young's modulus as a scalar field in space, which we denote E(x), with x a point in space. Our technique does not put any requirement on how the field is encoded (e.g. implicit function or interpolated from a grid), but it is expected that the field varies smoothly compared to the size of the Voronoi cells. Violating this expectation will not result in an incorrect structure, but the produced elasticity gradient will not be a good match to the input field.

The target Young's modulus field is converted upon lookup into a target density following the approach proposed in Section 4.2 -this requires a simple tabulation computed from Figure 5. This directly drives ρ in Algorithm 1. Unless otherwise specified, the beam radius remains fixed at τmin.

There are two main usage scenarios for our technique. A first scenario is to fill the inside of objects that have to remain rigid -simply adapting the inner density to varying stresses. In such a case the outer hull of the object remains solid and there is no additional challenge. The second scenario is to produce objects that can deform -flexible prosthesis and robot parts, toys. In such a case, printing the outer hull of the object would be detrimental to its flexibility. We therefore propose an object frame generation well suited to our approach. This frame is visible on all our 3D printed objects. It is fully procedural and only assumes that we have access to a (narrow band) distance field from the surface.

Object Frame. The key idea of our frame generation is to intersect the faces of the Voronoi diagram with a thick shell just below the surface, while taking care of cases where faces are almost parallel to the surface. Figure 11a explains our frame generation process. q is the query point and we want to know if it belongs to the solid structure or not. Here we consider two seeds s1, s2 among the closest seeds from q as obtained from Algorithm 1. If the point q is too close from the surface, we need to keep the intersection of the Voronoi faces with the frame (in addition to the regular Voronoi edges). Simply using the distance from q to the bisector of {s1, s2} would yield the region colored in light orange, whereas we would like to select the point q iff it belongs to the beam of center c and thickness t (in light purple). However, we don't know the exact position of c, but only the distance to the border d, the distance to the bisector x, and an estimate of the angle formed by the surface and the bisector at point q, that we call α. Assuming the surface is locally planar the rest follows from basic trigonometry: knowing l1 + l2 = d sin(α) , and l2 = x tan(α) , we can compute y and test whether x 2 + y 2 t 2 . Printed Results. We now apply our approach to produce a variety of 3D printed results. We use two different printers: a B9Creator with red-cherry resin and an Autodesk Ember with standard clear resin. All these prints are prepared using our image based in-house slicer using the implicit procedural foam generation.

Our SLA printers require support structures. The objects we printed do not themselves require support -a property preserved by our microstructures but for a few cases along the object frame. These cases are rare, for illustration we selected one in Figure 11b (top of rightmost case). Table 1 shows the amount of volume filtered out due to those constraints. We filter them out during out-of-core slice generation, keeping track of connected components from one slice to the next. Adapting support techniques to our microstructures is left as future work -a possibility would be to connect a standard support to the closest microstructure beam. Figure 12 shows the location of isolated voxels filtered out before printing on SLA printers. SLS printers do not require support and can directly print our objects.

Figure 12: Illustration of the isolated voxels filtered out for resin printers. Removed percentage for each example is reported in Table 1. Greyscaled values in voxels (section 3.2) have been thresholded at isovalue 0.5 to show a binary design.

Figure 13 shows a very simple case of graded material applied to a 3D model. Figure 14 is a detailed 3D model filled with a spatially varying structure, denser along the surface and coarser inside. Note how the frame preserves the surface details. Figure 17 is a similar case on a sliced sphere, revealing the deformation behavior. Figure 18 is a model with varying elasticity. The model deforms as expected when pushing its head sideways. The object frame has a limited impact on the deformation behavior, while producing a much more visually pleasing object. It is however difficult to quantify precisely the mechanical influence of the frame. Figure 15 illustrated an android finger with built-in flexible articulations. Figure 1 shows a more complex case of a painted object, where we added a support manually. Figure 16 shows how the elasticity field can be controlled, e.g. by the SIGGRAPH logo. Figure 19 illustrates how the elasticity field can be used to produce additional effects, such as anisotropic behaviors.

Table 1 gives detailed statistics on all these prints, including average time per-slice. The throughput of our implementation averages to 7.8M pixel/sec on the slice images on an Intel R Core TM i7-4770K @ 3.50GHz, 16 GB RAM with a Titan Black NVidia GPU.

Large Objects. Our approach scales trivially with object size. While we are currently limited in the size of objects we can print, we illustrate this by producing a microstructure in a model that is about one meter in size. The sliced equivalent has 5 tera-voxels -but of course only a single slice (15402 × 19814) would have to fit in the printer memory.

Figure 20 provides closeups of the depth map obtained by raycasting the implicit structure generated by our algorithm.

Discussion and Limitations

Our technique has a number of limitations. The foams only exhibit their target properties when a sufficiently large volume is printed -a limitation shared by all approaches relying on homogenization. This is however aligned with our goal of producing dense microstructures in large objects.

Compared to regular structures the stochastic nature of the foams produces more localized stresses. We have observed that a few beams fail under large deformations, perhaps earlier than on regular structures. This is specially the case along object frames with the Ember standard resin, which is more brittle. This would require further studies, noting that crushing behaviors of naturally occurring open-cell foams have received some interest [START_REF] Gaitanaros | On the crushing response of random open-cell foams[END_REF]].

Contrary to tile-based approaches we currently cannot provide spatially varying Poisson's ratios. Studying stochastic structures that can vary both Young's modulus and Poisson's ratio is an interesting direction of future work. Generating structures with anisotropic behaviors is another natural venue for further studies.

Conclusion

The main advantages of our approach stem from the implicit formulation of stochastic microstructures. This has significant computational advantages by allowing for evaluation at slicing time and by avoiding to resort on global optimizations for each new object. The aperiodic and stochastic nature of the foams provide a simple and efficient way to grade the structures and to conform to target elasticity fields in space, without introducing artificial boundaries.

While we took a strict interpretation on the procedural nature of our structure generation, it is clear that other computational schemes could be envisioned. The important fundamental properties are 1) that the evaluation of the structure remains local and independent from the overall size of the domain and 2) that the elastic properties relate to the structure parameters through a simple relationship avoiding complex parameter fitting during evaluation.

Our approach deviates significantly from both the periodic tiling of microstructures and the optimization of macrostructures, by making a link between microstructures and procedural solid textures with controlled statistics in computer graphics. We believe there are many other such structures to be discovered, and hope our work will spark further interest in procedurally generated, stochastic microstructures.

  Voronoi edge beams.

Figure 2 :

 2 Figure 2: 2D overview of Algorithm 1. (a) Input query point q and density field ρ. (b) The seeds that could contribute to the Voronoi cell of N [0](which is the seed closest to q) are gathered (Algorithm 2). (c) The bisectors of the seed pairs influencing q are computed. (d) Finally, the algorithm checks whether q lies inside a beam of radius τ along the Voronoi edge.

Figure 3 :

 3 Figure 3: In this counterexample, all the seeds lie on the same 3D plane (cross section shown). The closest point pl to the bisector line of N [0], N [i], and N [j] is at a distance τ -, where > 0. However, pl does not belong to any edge of the Voronoi cell of N [0] (shaded in blue).

7Figure 4 :

 4 Figure 4: Left: We consider the Voronoi cell of a seed located anywhere in the green square, and how it is influenced by another seed in the red square. The red grid square is within the second ring of neighbors of the green grid square.Each pair of seeds in the green/red squares defines a possible bisector, which might be a face of the Voronoi cell of the green seed. The opposite half-space cannot belong to the Voronoi cell. This is illustrated for a single pair of seeds in the figure (green/red dots). We define the shadow of the red square as the intersection of all the half-spaces from all possible pairs of green/red seeds. The shadow is shown in gray ; it represents the region of space that cannot possibly be part of the Voronoi cell of the green seed: regardless of the position of the seeds, the shadow is always cutout by the bisectors. Right: The shadows of all grid squares in the 2-ring completely cover the space beyond. Therefore, no seed outside of the 2-ring can have an influence on the Voronoi cell of the green seed.

Algorithm 3 :

 3 SubdivideCell Input: Density field ρ, cell center c, cell size l. Output: A set of seed, with a density driven by ρ 1 N ← ∅ ; 2 t ← l 3 × ρ(c) ; // target number of seeds in current cell 3 if t 2 3 then 4 I = RandomP ermutation(subcells); 5 nmin = t ; // minimum number of samples to draw 6 for i ← 0 to nmin -1 do 7 N ← N ∪ {RandomSampleInSubcell(I[i])}; 8 p ← random(0, 1); 9 if p (t -nmin) then 10 N ← N ∪ {RandomSampleInSubcell(I[nmin])}; 11 else 12 for subcell ∈ currentCell do 13 N ← N ∪ SubdivideCell(ρ, subcell.center, l/2); 14 return N Beam Radii. The radii of the beams are directly controlled by τ during evaluation (Algorithm 1, line 6

Figure 5 :

 5 Figure 5: Young's modulus as a function of radius and density.Each dot is the result of homogenization for a given (ρ, τ ). For each (ρ, τ ) we generate multiple instances of foams using different random seeds. As can be seen the dots tightly cluster, indicating that randomness has little impact on the elastic properties. The maximum deviation after the fitting is 3.3% of the Young's modulus. The data points are used to fit a polynomial function (see Section 4). For clarity, the graph shows a subset of our data that extends to higher densities for lower beam radii.

Figure 6 :

 6 Figure 6: Relative volume -of the structure over simulated volume -and error ξ(C H ) of the approximation of the ideal isotropic tensor (i.e. divergence from isotropy).As the spatial extent of the simulated volume grows the tensor tends towards the ideal isotropic tensor.

Figure 7 :

 7 Figure 7: Relative volume ω and volume of pockets (%).Pockets have a negligible impact under a relative volume of 0.7 and go up to 4% of the full volume around 0.9 relative volume. Thanks to the quasi-linear relationship between relative volume and Young's modulus, we can expect at most a similar error on the elastic behavior.

Figure 8 :

 8 Figure 8: Relative Young's modulus and Poisson's ratio. The base material has a Poisson's ratio of 0.3. The Poisson's ratio tends to 0.5 for low densities, and to the base material Poisson's ratio for high densities.

[Figure 9 :

 9 Figure 9: Relative distance ξ to an ideal isotropic tensor, for a fixed voxel size and a resolution of 80 × 80 × 80. Both pairs have a similar Young's modulus. The small isotropy error of the tile from [Panetta et al. 2015] (left) is due to limited numerical precision and decreases with finer discretization.As larger volumes are simulated, the isotropy of our structures quickly improves (see also Figure6).

Figure 10 :

 10 Figure 10: Compression tests on printed samples. The measured weight is expected to be linearly correlated with the target Young's modulus.

  Scheme of the query. (b) A cylinder with 3 possible frames.

Figure 11 :

 11 Figure 11: Top: 2D planar cut orthogonal to the surface (top black line). The goal is the produce the beam which cross section is drawn as a circle. The purple line is a Voronoi face intersected by the view plane. Bottom: Extreme case of a frame on a small cylinder. From left to right, no frame, frame obtained by intersecting the Voronoi faces with the surface naively, frame obtained with our approach.

Figure 14 :

 14 Figure 14: Knight model with a dense crust, and lowdensity interior. Model: Knight (thing:33804) by andreas.

Figure 15 :

 15 Figure 15: Finger with articulations. Right: deformed object. Note that the inferior articulation is made wider and thus bends further. Model: Olivier hand (AIM@SHAPE Project).

Figure 16 :

 16 Figure 16: SIGGRAPH logo driving the density field.

Figure 17 :

 17 Figure 17: A dome with variable density. No frame is generated on the cut revealing the interior.

Figure 18 :

 18 Figure 18: Cute octopus. From left to right: 3D model (original), printed model, printed with deformation (finger pressure). Model: Cute Octopus (thing:27053) by MakerBot.

Figure 19 :

 19 Figure 19: Varying the density in one direction only can also produce parts with anisotropic behavior. The same weight of 8 kg has been applied to the cube placed in different orientations, leading to a visible deformation in the direction of the anisotropy.

Table 1 :

 1 Statistics on the examples shown in the paper. Beam thickness varies between 0.2mm and 0.4mm.

	Example		Extent (mm)	# Voxels	Volume % Filtered	Time per slice (ms)
	Moomin	fig. 1	26.7 × 40.8 × 51.9	534 × 815 × 1038	6.44%	0.005%	68.34
	Ellipsoid	fig. 13	30.9 × 30.9 × 41.1	617 × 617 × 822	6.30%	0.001%	37.28
	Knight	fig. 14	26.1 × 30.0 × 50.55	521 × 600 × 1011	12.50%	0.023%	20.25
	Finger	fig. 15	25.0 × 23.25 × 70.5	500 × 465 × 1410	23.35%	0.006%	28.03
	SIGGRAPH logo fig. 16	20.0 × 40.0 × 80.0	400 × 800 × 1600	5.73%	0.003%	69.18
	Half-dome	fig. 17	25.0 × 50.0 × 25.0	500 × 1000 × 500	19.49%	0.025%	71.22
	Octopus	fig. 18	41.7 × 41.1 × 28.8	833 × 822 × 576	17.27%	0.009%	150.22
	Anisotropic cube fig. 19	40.0 × 40.0 × 40.0	800 × 800 × 800	26.86%	0.005%	113.52
	Forest dragon	fig. 20 770.1 × 990.7 × 961.7 15402 × 19814 × 19234	N/A	N/A	1666.91

flexible rigid Figure 13: A simple graded ellipsoid.

Acknowledgments

We would like to thank anonymous reviewers for their valuable comments, and the members of our research team for help with printing. This work was supported by ERC grant ShapeForge (StG-2012-307877).