David Jourdan
email: david.jourdan@inria.fr

Shrink & Morph: 3D-printed self-shaping shells actuated by a shape memory effect

Fig. 1. Given a target surface (left), our method optimizes a dense set of deposition trajectories that form a printed plate. Once the plate is immersed in hot water, it deforms and closely approaches the target surface (right).

While 3D printing enables the customization and home fabrication of a wide range of shapes, fabricating freeform thin-shells remains challenging. As layers misalign with the curvature, they incur structural deficiencies, while the curved shells require large support structures, typically using more material than the part itself.

We present a computational framework for optimizing the internal structure of 3D printed plates such that they morph into a desired freeform shell when heated. This exploits the shrinkage effect of thermoplastics such as PLA, which store internal stresses along the deposition directions. These stresses get released when the material is heated again above its glass transition temperature, causing an anisotropic deformation that induces curvature.

Our inverse design method takes as input a freeform surface and finds an optimized set of deposition trajectories in each layer such that their anisotropic shrinkage deforms the plate into the prescribed surface geometry. We optimize for a continuous vector field that varies across the plate and within its thickness. The algorithm then extracts a set of deposition trajectories from the vector field in order to fabricate the flat plates on standard FFF printers. We validate our algorithm on freeform, doubly-curved surfaces.

CCS Concepts: • Computing methodologies → Computer graphics; • Applied computing → Computer-aided design.

Additional Key Words and Phrases: Additive Fabrication, Deployable Structures, Shape Memory Polymers.

INTRODUCTION

Freeform thin shell surfaces appeal to both designers and engineers for their aesthetics and performance. They are commonly used to form the outer shell of home appliances and handheld devices. However, they involve a significantly more complex and costly fabrication process than their flat or single-curved counterparts.

In particular, the additive manufacturing of thin shells usually requires large volumes of supports, while misalignment between shell curvature and printing layers leads to staircase artifacts and poor structural performance (see Fig. 2).

One approach to circumvent these limitations is to consider selfshaping surfaces: plates fabricated flat that deform into a specific geometry in response to different stimuli such as heat, pH, or moisture. These approaches avoid the aforementioned issues, show great potential for reducing fabrication costs, and allow postponing deployment until on-site placement.

Heat

In this work, we fabricate selfshaping surfaces with Fused Filament Fabrication (FFF), one of the least expensive and most widely used 3D printing techniques. We propose to use the shape memory effect of thermoplastics such as polylactic acid (PLA): polymer chains in the material tend to align along deposition trajectories when printed, their trajectory-aligned distribution turns to a more isotropic state when the material is reheated above its so-called glass transition temperature 𝑇 𝑔 ≈ 60 • C causing the thermoplastic material to shrink along the printing direction. Our experiments verify and investigate further this physical anisotropic property, which is crucial to our method. As detailed in Section 2, studies have shown that parameters such as speed, layer thickness, and printing temperature can influence the amount of length change induced. However, these parameters cannot be easily varied locally within the plate, and understanding their effect on the shrinking behavior requires precise calibration procedures.

Instead, we propose to optimize the geometry of the printing trajectories: by varying the local direction of the trajectories, we control the direction of the anisotropic deformation (see Figure 1). However, controlling the in-plane deformation is insufficient, as the same 2D deformation map may result in different 3D shapes. Conversely, curved shapes can exist without in-plane deformations, e.g., developable surfaces.

To gain more control over curvature, we propose to vary the deformation direction not only in-plane but also in the normal direction -the plate thickness -by having printing trajectories with different orientations from layer to layer.

We devise an inverse design method that takes as input a freeform surface and outputs an optimized set of deposition trajectories, such that when the fabricated plate is immersed in hot water, it deforms into the prescribed surface geometry. We optimize a continuous vector field that varies within the plate thickness to obtain freeform surfaces with non-trivial curvature. After optimization, the algorithm extracts a set of deposition trajectories from the vector field in order to fabricate the flat plates. Remarkably, by tailoring the path orientation along the plate layers, our method can control both Gaussian and normal curvature using a single base material.

Contributions. In summary, we introduce:

• An end-to-end framework for 3D printing self-shaping thinshells, capable of reproducing freeform surfaces. This is achieved through the sole control of the geometry of the printing trajectories within and across layers. • A novel formulation for morphing shape-memory polymers based on the theory of non-Euclidean plates that allows for variation in the normal direction of the plate to control both Gaussian and normal curvature. This takes into account specificities of 3D printing, in particular the gradation of shrinkage ratios between top and bottom layers. • A precise calibration procedure based on a study on the shapememory properties of PLA with respect to parameters such as heat, layer height, and printing speed.

An open-source implementation of our method is available at https://github.com/DavidJourdan/shrink-morph

RELATED WORK

We call self-shaping surfaces initially flat plates that deform into curved 3D surfaces without any external actuation forces. Such surfaces can result from natural processes such as growth or drying processes or can be manufactured to react to a given stimulus that will drive the deformation.

Heat

Fig. 2. Layer by layer fabrication (right) leads to staircase artifacts and requires support, while our method (left) enables to have layers aligned with the surface and does not need support.

Manufactured self-shaping surfaces come in a variety of forms and materials. Examples include hydrogels [START_REF] Sydney Gladman | Biomimetic 4D Printing[END_REF][START_REF] Kim | Designing Responsive Buckled Surfaces by Halftone Gel Lithography[END_REF][START_REF] Na | Grayscale gel lithography for programmed buckling of non-Euclidean hydrogel plates[END_REF]] and nematic elastomers [START_REF] Aharoni | Universal inverse design of surfaces with thin nematic elastomer sheets[END_REF]] which can expand when heated or exposed to water, hygroscopic wood-based materials [START_REF] Cheng | Programming material compliance and actuation: Hybrid additive fabrication of biocomposite structures for large-scale self-shaping[END_REF][START_REF] Grönquist | Analysis of hygroscopic self-shaping wood at large scale for curved mass timber structures[END_REF], pre-tensioned membranes combined with rigid elements [START_REF] Guseinov | Programming temporal morphing of self-actuated shells[END_REF][START_REF] Guseinov | CurveUps: Shaping Objects from Flat Plates with Tension-Actuated Curvature[END_REF][START_REF] Jourdan | Computational Design of Self-Actuated Surfaces by Printing Plastic Ribbons on Stretched Fabric[END_REF][START_REF] Pérez | Computational Design and Automated Fabrication of Kirchhoff-Plateau Surfaces[END_REF]] and more exotic materials such as pasta which self-forms when cooked [START_REF] Tao | Morphlour: Personalized[END_REF][START_REF] Tao | Morphing pasta and beyond[END_REF][START_REF] Wang | Transformative Appetite: Shape-Changing Food Transforms from 2D to 3D by Water Interaction through Cooking[END_REF].

Self-shaping PLA structures. Among these possible techniques, exploiting the shape-memory properties of PLA in fused filament fabrication presents the advantage of being easily accessible through the widely available FFF 3D printers. This leads to direct applications for educational, hobbyist, and repair projects, as well as ongoing research in 4D printing and robotics.

Many approaches rely on fabricating discrete elements that behave as actuators. These elements use a predefined trajectory pattern producing a known, calibrated behavior. The elements are then instantiated at specific positions and orientations to guide the deformation of the surrounding passive parts. Some techniques rely directly on shrinkage [START_REF] Chalissery | Highly Shrinkable Objects as Obtained from 4D Printing[END_REF][START_REF] Kongpyung | Shrink-Cells: Localized and Sequential Shape-Changing Actuation of 3D-Printed Objects via Selective Heating[END_REF][START_REF] Wang | 4DMesh: 4D Printing Morphing Non-Developable Mesh Surfaces[END_REF] while others rely on a bilayer effect to induce curvature [START_REF] Teunis Van Manen | Programming 2D/3D shape-shifting with hobbyist 3D printers[END_REF][START_REF] Wang | A-Line: 4D Printing Morphing Linear Composite Structures[END_REF][START_REF] Yu | Material characterization and precise finite element analysis of fiber reinforced thermoplastic composites for 4D printing[END_REF]: different predefined patterns with different shrinking behaviors are used for the top and bottom layers. For instance, An et al. [2018] uses two different materials -shrinking PLA and non-shrinking TPU -while [START_REF] Sun | ShrinCage: 4D Printing Accessories That Self-Adapt[END_REF] only heats specific parts using conductive PLA. As opposed to discrete elements, we explore a controlled, continuous deformation throughout the plate to obtain thin curved shells.

The Geodesy approach [START_REF] Gu | Geodesy: Self-Rising 2.5D Tiles by Printing along 2D Geodesic Closed Path[END_REF][START_REF] Gu | Inverse Design Tool for Asymmetrical Self-Rising Surfaces with Color Texture[END_REF] is closest to our objectives but is, however, limited to height-fields as the deposition trajectories always align with the isolines of the input field. In addition, the contraction rate is controlled by continuously varying layer thicknesses. This necessitates a material-specific calibration to precisely quantify the relationship between layer height and shrinking ability. This makes the method difficult to implement on commercial PLA filaments, which exhibit a wide range of behaviors. Although our method also needs some calibration, we only need to find the shrinking ratio for one set of printing parameters. [START_REF] Deshpande | Fab4D: An Accessible Hybrid Approach for Programmable Shaping and Shape Changing Artifacts[END_REF] experiment with controlling the orientation of the deposition trajectories. However, this is not incorporated into an inverse design method. This work demonstrates that heating elements can be locally embedded inside the surface to locally trigger shrinkage.

To the best of our knowledge, our method is the first to optimize for the orientation of trajectories across layers in order to achieve a prescribed, continuous deformation.

Thermomechanical studies of shape memory PLA. The shape memory effect of PLA is a complex phenomenon that depends on various factors. Several material science researchers studied the influence of different parameters on the deformation, such as nozzle temperature [START_REF] Chalissery | Highly Shrinkable Objects as Obtained from 4D Printing[END_REF], build plate temperature [START_REF] Kačergis | Influence of fused deposition modeling process parameters on the transformation of 4D printed morphing structures[END_REF], printing speed [START_REF] Kačergis | Influence of fused deposition modeling process parameters on the transformation of 4D printed morphing structures[END_REF][START_REF] Ravichandra | Additive manufacturing-enabled shape transformations via FFF 4D printing[END_REF][START_REF] Teunis Van Manen | Programming 2D/3D shape-shifting with hobbyist 3D printers[END_REF][START_REF] Zhang | Pattern Transformation of Heat-Shrinkable Polymer by Three-Dimensional (3D) Printing Technique[END_REF], thickness of the plate [START_REF] Ravichandra | Additive manufacturing-enabled shape transformations via FFF 4D printing[END_REF] or the flow of material [START_REF] Teunis Van Manen | Programming 2D/3D shape-shifting with hobbyist 3D printers[END_REF].

As the experiments depend on numerous factors, the results are not always consistent across the literature. Our own experiments showed that results may vary widely depending on the specific PLA filament used, with variations across vendors and colors. This is likely due to the influence of specific additives and pigments, whose exact composition is undisclosed. We provide the parameters used for our examples on specific, commonly available brands of PLA and detail in Section 4 our own experiments to find print parameters maximizing the shrinking performance.

Anisotropic parameterization for deployable structures. The flattening operation required for the inverse design problem to go from the curved target surface to the flat plate to be printed can be represented as a parameterization of a surface embedded in R 3 onto a flat, 2D domain. When materials are deployed from locally isotropic deformations -for instance, with auxetic materials -conformal maps have been shown to be the proper framework to geometrically capture their deformation [START_REF] Konaković | Beyond Developable: Computational Design and Fabrication with Auxetic Materials[END_REF]] and to build inverse design tools [START_REF] Chen | Bistable Auxetic Surface Structures[END_REF][START_REF] Konaković-Luković | Rapid Deployment of Curved Surfaces via Programmable Auxetics[END_REF].

When materials morph from anisotropic deformation fields, however, conformal maps cannot be used, and one has to resort to different formalisms. For instance, [START_REF] Griniasty | Curved Geometries from Planar Director Fields: Solving the Two-Dimensional Inverse Problem[END_REF] employed hyperbolic PDE solvers to solve the inverse design problem for nematic elastomers. This technique is not readily applicable to our case as it does not model a change of anisotropy across layers required to account for extrinsic curvature.

Another strategy used for the inverse design of nematic elastomers is to use a local-global solver to find a parameterization much in the same way as ARAP parameterization [START_REF] Aharoni | Universal inverse design of surfaces with thin nematic elastomer sheets[END_REF]. [START_REF] Panetta | Computational Inverse Design of Surface-Based Inflatables[END_REF] use this local-global solver to initialize their non-linear optimizer to design inflatable plates with oriented channels triggering curvature. We employ a similar optimization strategy with a very different simulation model: to account for the fact that the channels are much easier to bend across the air tubes than along, they add a custom term in their optimization, whereas our simulation model for anisotropic, non-euclidean plates allows us to find the proper scaling of the bending terms without introducing an additional free parameter.

Curved printing, conformal slicing. A different take on the fabrication of shells is to rely on curved printing [START_REF] Zhang | 𝑆 3 -Slicer: A General Slicing Framework for Multi-Axis 3D Printing[END_REF]]. In particular, [START_REF] Hong | Open5x: Accessible 5-Axis 3D Printing and Conformal Slicing[END_REF] propose a method to 3D print shells using conformal slicing: a 5-axis printer prints a curved surface on top of a PVA support structure that dissolves in water. This, of course, requires a custom 3D printer with 5-axis capabilities and does not remove the need for large volumes of support (compared to the fabricated shell volume).

OUR APPROACH

An overview of our inverse design approach is given in Figure 3. Our method takes as input a triangle mesh embedded in R 3 representing the desired target (deployed) surface shape and computes (1) a flattening of the surface into the plane, (2) a three-dimensional vector field covering the plane and varying across its thickness, and (3) outputs a set of deposition trajectories following the vector field for fabrication. We give below a brief overview of each step.

Flattening. We first obtain a planar parameterization of the input surface, which is suitable for optimizing the trajectory orientations. We use a local/global iterative algorithm akin to an As-Rigid-As-Possible parameterization [START_REF] Liu | A Local/Global Approach to Mesh Parameterization[END_REF]. A similar approach was used for nematic elastomers and inflatable structures [START_REF] Aharoni | Universal inverse design of surfaces with thin nematic elastomer sheets[END_REF][START_REF] Panetta | Computational Inverse Design of Surface-Based Inflatables[END_REF].

Vector field modeling. We solve for a vector field represented as plane angles 𝜃 throughout the plate's volume. However, the plate being thin -its thickness much smaller than its other dimensionswe use two angle fields 𝜃 1 and 𝜃 2 to model it with a linear variation across the thickness:

𝜃 (𝑥, 𝑦, 𝑧) = 𝜃 1 (𝑥, 𝑦) + 𝑧 ℎ 𝜃 2 (𝑥, 𝑦), 𝑧 ∈ [-ℎ 2 , ℎ 2 ]
where ℎ is the plate thickness.

Vector field optimization. The flattening procedure gives us an initial guess for the vector field, which is used to initialize our inverse optimization method determining optimized printing orientations. In particular, we obtain the value of 𝜃 1 from this flattening procedure.

The heart of our technique is an accurate forward simulation method that determines the shape that would be obtained given the vector field. We derive it from the theory of non-Euclidean plates [START_REF] Efrati | Elastic theory of unconstrained non-Euclidean plates[END_REF]] as detailed in Section 3.1.

We use the forward approach to find an equilibrium shape closest to the target, using the adjoint method from sensitivity analysis to optimize for the vector field. We detail the inverse optimization in Section 3.2. A key novelty of our technique is to model the anisotropic bending in the optimization process. This is essential to capture the material behavior adequately.

Trajectory extraction. In each fabrication layer, we generate equally spaced trajectories aligned with the optimized orientations. We rely on the Stripe Patterns algorithm [START_REF] Knöppel | Stripe Patterns on Surfaces[END_REF]]. The entire process is described in Section 3.2.3.

Forward problem

In order to tackle the inverse optimization, we first have to solve the forward problem: given trajectory orientations in a flat 3D printed plate, we determine its final configuration after re-heating above the glass transition temperature 𝑇 𝑔 .

More precisely, given an initial flat configuration Ŝ ∈ R 2 and a scalar field 𝜃 representing the orientations, we have to simulate how lengths on the surface change locally so as to obtain the final curved shape. This problem is reminiscent of similar phenomena where flat sheets deform under environmental changes (e.g., drying leaves, moist paper, etc. see [START_REF] Chen | Physical Simulation of Environmentally Induced Thin Shell Deformation[END_REF]). These phenomena can be approximated geometrically as a metric change within the rest configuration of the material. This initially flat plate therefore gets imbued with a reference configuration whose metric may represent a surface that is not embeddable in R 3 , which is why we refer to these structures as non-Euclidean plates [START_REF] Efrati | Elastic theory of unconstrained non-Euclidean plates[END_REF][START_REF] Pezzulla | Curvature-driven morphing of non-Euclidean shells[END_REF]]. We begin by describing in Section 3.1.1 a general elastic energy formulation (often attributed to W.T. Koiter) whose derivations can be found in more detail in [START_REF] Efrati | Elastic theory of unconstrained non-Euclidean plates[END_REF]]. We explain how to adapt this model for anisotropic materials and derive the exact quantities specific to our problem in Section 3.1.2. Finally, we solve the problem using a discretization similar to [START_REF] Sperl | Homogenized Yarn-Level Cloth[END_REF].

3.1.1 Theory of non-Euclidean plates. Let us consider a deformation from an initially flat plate of thickness ℎ to a curved surface, defined by the mapping Φ :

V = Ŝ × [-ℎ 2 , ℎ 2 ] → R 3 .
We define the curved surface as S = 𝜙 ( Ŝ) with 𝜙 being a parameterization of the midsurface (see Figure 4).

The metric g = dΦ 𝑇 dΦ measures the deformations in lengths and angles between the initial flat domain and the deformed domain. The elastic energy of the plate depends on its rest state and the corresponding metric ḡ measuring the deformation between the initial and rest configurations. Usually, the rest state and the initial state coincide, making ḡ the identity. However, plastic deformations like the shrinkage of PLA under heat modify the rest shape. The resulting rest state is often not embeddable in R 3 : there is likely no shape in R 3 that would not have residual stresses. However, even though the rest shape cannot exist in R 3 it remains possible to define ḡ for our case of plastic deformation, as detailed in Section 3.1.2.

Using standard plane-stress assumptions, we find that the metrics g and ḡ can be reduced to their planar 2 × 2 sub-parts C and C corresponding to the deformation of the mid-surface. Assuming the thickness ℎ is negligible compared to the dimensions of the surface Ŝ, C and C can be expanded into first and second fundamental forms:

C = a -2𝑧b + O (𝑧 2 ), C = ā -2𝑧 b + O (𝑧 2 ), (1)
Integrating the total elastic energy 𝑊 = ∫

V 1 2 (C -C) 2
M d𝑉 in the thickness direction gives the following expansion:

𝑊 = ℎ 4 ∫ Ŝ ∥a -ā∥ 2 M d𝐴 + ℎ 3 12 ∫ Ŝ b -b 2 M d𝐴 + O (ℎ 4 ). (2)
As the plate is non-Euclidean, there is not necessarily a shape in R 3 such that 𝑊 = 0. In the equation above, ∥•∥ M is a material norm.

Previous methods [START_REF] Chen | Physical Simulation of Environmentally Induced Thin Shell Deformation[END_REF][START_REF] Weischedel | A discrete geometric view on shear-deformable shell models[END_REF]] employ a so-called isotropic Saint-Venant-Kirchhoff norm, suitable for large displacements but small strains. However, since in our case the material exhibits an orthotropic response, we generalize this StVK norm so as to model orthotropic materials with Young's moduli 𝐸 1 in the orientation 𝜃 and 𝐸 2 in the orthogonal direction. The material norm of a matrix 𝝐 is then written as:

∥𝝐 ∥ M = (𝑅 𝑇 𝜃 ā-1 𝝐𝑅 𝜃 ) : C : (𝑅 𝑇 𝜃 ā-1 𝝐𝑅 𝜃 ) (3)
where C is the rank-4 linear elastic tensor of the orthotropic material.

We use the parameterization of [START_REF] Li | Stable Orthotropic Materials[END_REF] who write this tensor as a 3 × 3 symmetric matrix using Voigt notation:

C = 1 1 -𝜈 2 𝐸 1 𝜈 √ 𝐸 1 𝐸 2 0 𝜈 √ 𝐸 1 𝐸 2 𝐸 2 0 0 0 1-𝜈 2 √ 𝐸 1 𝐸 2 , ( 4
)
with a Poisson's ratio-like parameter 𝜈 which we set to 0.4 as it is a reasonable value for PLA. We experimentally determine an approximate ratio 𝐸 1 /𝐸 2 = 10 for the material we work with (see Section 4.3 for a discussion).

To find the forms a and b, the deformation map Φ is computed from the parameterization of its mid-surface 𝜙 and its normal map 𝑛:

Φ(𝑥, 𝑦, 𝑧) = 𝜙 (𝑥, 𝑦) + 𝑧𝜌𝑛(𝑥, 𝑦).
The growth factor 𝜌 accounts for the shell being allowed to scale in the normal direction. Most previous work on non-Euclidean surfaces [START_REF] Chen | Physical Simulation of Environmentally Induced Thin Shell Deformation[END_REF][START_REF] Efrati | Elastic theory of unconstrained non-Euclidean plates[END_REF]] assumes that the thickness of a plate stays constant during the deformation, making 𝜌 = 1. However, this is not the case for the shape memory effect of PLA. We instead follow [START_REF] Pezzulla | Curvature-driven morphing of non-Euclidean shells[END_REF] and consider that the material has a constant growth factor 𝜌 = 𝜆 3 along this direction, which we measure in Section 4. The thin-shell assumption therefore gives:

C = d𝜙 𝑇 d𝜙 + 2𝜌𝑧d𝑛 𝑇 d𝜙 + O (𝑧 2 ) = a -2𝑧b + O (𝑧 2 ).
(5)

This equation allows to identify a = d𝜙 𝑇 d𝜙 and b = -𝜌d𝑛 𝑇 d𝜙, with d𝜙 being the deformation gradient and d𝑛 is the so-called shape operator (see e.g. [START_REF] Grinspun | Computing discrete shape operators on general meshes[END_REF]). Note that the integration in Equation 2 is done over the initial shape. The fundamental forms a and b are computed on Ŝ as pullbacks of the fundamental forms used in, e.g., [START_REF] Chen | Physical Simulation of Environmentally Induced Thin Shell Deformation[END_REF]]. We choose to use these pullbacks and define everything on Ŝ as it makes the rest metric easier to define in this setting.

3.1.2 Simulation of shape memory PLA. As explained earlier, the printed PLA shrinks in the direction of the trajectories when heated. A secondary effect is that it also grows in the other directions. In the following paragraph, we estimate the rest fundamental forms ā and b by measuring this deformation.

𝜃 𝜆

1 𝜆 2
We call 𝜆 1 and 𝜆 2 the ratios measuring the in-plane distortion (i.e. the change of length) respectively in the path direction (represented by an angle 𝜃 ) and in the orthogonal direction, and 𝜆 3 the ratio measuring the change in length along the thickness in the normal direction. We explain in Section 4.1 how to measure experimentally these ratios, which are constant through the plate.

Noting that a vector oriented along 𝜃 will change its length by a factor 𝜆 1 and a vector in the perpendicular direction will change its length by a factor 𝜆 2 , we write:

C = 𝑅 𝜃 𝜆 2 1 0 0 𝜆 2 2 𝑅 𝑇 𝜃 . (6)
Assuming that 𝜃 is a linear function of 𝑧, we can write it as 𝜃 (𝑥, 𝑦, 𝑧) = 𝜃 1 (𝑥, 𝑦) + 𝑧 ℎ 𝜃 2 (𝑥, 𝑦) where 𝜃 1 is the mean angle value and 𝜃 2 represents the variation along the thickness, and make the following expansion

C = 𝑅 𝜃 1 𝜆 2 1 0 0 𝜆 2 2 𝑅 𝑇 𝜃 1 - 𝑧 ℎ 𝜃 2 (𝜆 2 2 -𝜆 2 1 )𝑅 𝜃 1 0 1 1 0 𝑅 𝑇 𝜃 1 + O (𝑧 2 ), (7)
which allows us to identify ā and b.

3.1.3 Discretization. We simulate the deformation induced by the shape memory effect by minimizing the energy 𝑊 defined in Equation 2. We represent the mid-surface S as a triangular mesh composed of a set of vertices 𝑉 and a set of triangles 𝑇 . We consider the fields 𝜃 1 , 𝜃 2 to be constant over each triangle. A vector of 2D coordinates describes the position of the vertices x in the initial flat position and by the 3D coordinates x in a deformed state. We use piecewise linear finite elements so that the strain is constant on each triangle. We can then obtain the following discretized energy:

𝑊 (x) = ℎ 4 ∑ 𝑡 ∈𝑇 ∥a -ā∥ 2 M + ℎ 3 12 ∑ 𝑡 ∈𝑇 b -b 2 M √ det ā. (8)
ā and b are computed using Equation 7 using the value of 𝜃 1 and 𝜃 2 of the triangle 𝑡. a and b are derived from Equation 5 as a = 𝐹 𝑇 𝐹 and b = 𝐹 𝑇 Λ𝐹 . F is the deformation gradient obtained by the initial positions { v𝑖 } 1,2,3 of the vertices of triangle 𝑡 and their deformed positions {𝑣 𝑖 } 1,2,3 :

𝐹 = 𝑣 1 -𝑣 0 , 𝑣 2 -𝑣 0 v1 -v0 , v2 -v0 -1 , ( 9
)
and Λ is the shape operator defined on the deformed surface. We use the triangle-averaged shape operator from [START_REF] Grinspun | Computing discrete shape operators on general meshes[END_REF]:

Λ = ∑ 𝑖=1,2,3 𝜃 𝑖 𝐴ℓ 𝑖 𝜏 𝑖 𝜏 𝑇 𝑖 , (10)
with 𝐴 being the area of the triangle, 𝜃 𝑖 the signed dihedral angle between 𝑡 and its 𝑖 th neighbor, ℓ 𝑖 the length of the common edge 𝑒 𝑖 and 𝜏 𝑖 a vector obtained by rotating 𝑒 𝑖 of 90°around the normal of the triangle. The final shape induced by the shrinking of the plate can then be computed as x 𝑓 = arg min x 𝑊 (x). We minimize𝑊 using a projected Newton strategy, computing derivatives and projecting the pertriangle Hessians on the PSD cone using TinyAD [START_REF] Schmidt | TinyAD: Automatic Differentiation in Geometry Processing Made Simple[END_REF]].

Inverse problem

The forward problem lets us predict the curved surface obtained by heating a flat printed object. We now tackle the inverse problem: given a curved surface S representing the mid-surface of the plate, we want to find a flat shape and the printing trajectories such that it will deform into S when heated above glass temperature.

Namely, we are looking for a parametric domain Ŝ and the 𝜃 1 and 𝜃 2 fields defined on it such that, when applying the forward simulation, they deform Ŝ into S.

We decompose this problem in two parts. First, we consider that the direction field is constant across the thickness, setting 𝜃 2 = 0. This lets us optimize Ŝ and 𝜃 1 , only considering the in-plane deformations. Second, to accurately control the bending behavior of the surface, we optimize for orientation variations in the thickness 𝜃 2 . We use sensitivity analysis [START_REF] Zehnder | SGN: Sparse Gauss-Newton for Accelerated Sensitivity Analysis[END_REF], minimizing the distance between the simulated and target shapes (Section 3.2.2). Finally, from the field 𝜃 = 𝜃 1 + 𝑧 ℎ 𝜃 2 , we compute the printing trajectories for each plate layer.

Parameterization.

As a first step we find a mapping 𝜓 : S → R 2 whose induced distortions correspond to the shape memory effect of the PLA as close as possible. So the differential d𝜓 ∈ R 2×2 -parameterized by local barycentric coordinates on tangent planes -has singular values equal to 𝜆 -1 1 and 𝜆 -1 2 . By optimizing 𝜓 we also obtain the parameterization domain Ŝ = 𝜓 (S) as well as the average orientation 𝜃 1 -the orientation on the mid-surface -as the direction of the maximal distortion of 𝜓 .

We follow the local/global optimization scheme of [START_REF] Aharoni | Universal inverse design of surfaces with thin nematic elastomer sheets[END_REF]. The energy to minimize is:

𝐸 (𝜓, 𝑈 , 𝑉 ) = ∫ S d𝜓 -𝑈 𝜆 -1 1 0 0 𝜆 -1 2 𝑉 𝑇 2 𝐹 d𝐴. ( 11
)
The idea is to alternate between minimizing for𝜓 while keeping 𝑈 , 𝑉 constant on the one hand (the "global" problem) and minimizing 𝑈 and 𝑉 while treating d𝜓 as a constant on the other (the "local" problem). We follow the framework of [START_REF] Liu | A Local/Global Approach to Mesh Parameterization[END_REF] that we briefly remind below.

The global problem is solved by fixing 𝑈 and 𝑉 (from the initialization or a previous iteration) and minimizing for 𝜓 , it is solved by solving a linear system of 2𝑛 𝑉 linear equations -𝑛 𝑉 being the number of vertices-whose unknowns are the discretized positions x representing Ŝ in R 2 .

Then, by computing the singular value decomposition (SVD) of d𝜓 = 𝑅 𝛼 𝜎 1 0 0 𝜎 2 𝑅 𝑇 where 𝑅 𝛼 is a rotation of angle 𝛼 such that 𝑢 𝛼 = (cos 𝛼, sin 𝛼) is the direction of maximum shrink, the local problem admits the closed-form solution 𝑈 = 𝑅 𝛼 , 𝑉 = 𝑅.

Iterating between these two problems quickly converges to a discrete mapping 𝜓 matching the prescribed distortions in a least squares manner. Moreover we also get the associated flat domain Ŝ = 𝜓 (S) describing the flat plate to be printed and we can deduce the directions 𝜃 1 from the SVD of d𝜓 since 𝑢 𝛼 represents the direction of maximum shrinkage.

The local-global algorithm needs an initial parameterization. Unfortunately, different initializations can push the optimization to converge to vastly different local minima, some of which are not suitable for our fabrication purposes. In particular, we need a distribution of directions that are smooth enough to be able to trace continuous deposition trajectories. We found that first running asrigid-as-possible (ARAP) parameterization (i.e., running the same algorithm with 𝜆 1 = 𝜆 2 = 1) and then deforming the flat surface by a factor 𝜆 -1 1 along one axis and 𝜆 -1 2 provides a good initialization to this algorithm, which encourages smooth directions 𝑢 𝛼 in solutions.

In practice, this worked well across all our examples. In many cases, this initialization was enough to ensure smooth directions. If extra smoothness is required it is also possible to add to the optimization problem an explicit smoothness term as in [START_REF] Panetta | Computational Inverse Design of Surface-Based Inflatables[END_REF]. In their case, the smoothness term consists of the Dirichlet energy of the vector field of distortion directions 𝑢 2𝛼 (where 𝛼 is doubled to account for the 𝜋 symmetry of the 2-RoSy direction field):

min 𝜓 ∫ S 𝜎 1 -𝜆 -1 1 0 0 𝜎 2 -𝜆 -1 2 2 𝐹 d𝐴 + 𝑤 𝐷 ∫ S ∥∇𝑢 2𝛼 ∥ 2 d𝐴. ( 12
)
Since there is no obvious way to solve for this in the local-global framework, we propose to solve Equation 12 using a projected Newton strategy and computing the derivatives using TinyAD [START_REF] Schmidt | TinyAD: Automatic Differentiation in Geometry Processing Made Simple[END_REF]]. The difference between running this smoothing step or not is illustrated in Figure 5. As this smoothing step may not always be necessary, we show the cases where we use it in Table 1 as 𝑤 𝐷 = 0.1.

Directions optimization.

The previous optimization gives us both a domain Ŝ and average directions 𝜃 1 minimizing ∥a -ā∥ 2 M . We now seek to optimize orientation variations across the thickness 𝜃 2 . This parameter controls the shell mean curvature / bending. When optimizing for 𝜃 2 we keep 𝜃 1 and Ŝ fixed.

This problem cannot be solved by simply finding the fundamental forms ā and b which minimize Equation 2, i.e., directly minimizing b -b 2 M . Indeed, contrary to the ā case we cannot, in general, find 𝜃 2 such that b -b M ≈ 0. This means that even if we minimize the energy 𝑊 as a function of 𝜃 1 and 𝜃 2 in Equation 2, the target shape will not correspond to a stable equilibrium, and the simulation might significantly drift from the target shape.

We propose instead to minimize the distance between the target shape S and the result of the forward simulation applied on 𝜃 2 with 𝜃 1 , Ŝ fixed. We formulate the problem directly in the discretized domain. We consider vectors containing the position of all the vertices of the mesh. We call x 𝑇 the positions of the target shape, while x are the positions resulting from a forward simulation. We also use a vector 𝜽 2 containing the values of 𝜃 2 per vertex. The values of 𝜃 2 are averaged per triangle for the simulation.

The objective is to find the value of 𝜽 2 that minimizes the distance between x 𝑇 and x. The forward simulation results in a shape at equilibrium, giving us the constraint 𝜕𝑊 𝜕𝑥 = 0. We add two regularization terms to the minimization problem. First, 𝜽 2 𝑇 𝑀𝜽 2 (where 𝑀 is the vertex lumped mass matrix) ensures that 𝜃 2 is on the order of ℎ (and therefore small). Second, 𝜽 2 𝑇 𝐿𝜽 2 where 𝐿 is the cotan Laplacian. This ensures that the resulting field 𝜃 remains integrable. We obtain the following constrained optimization problem:

min 𝜽 2 (x -x 𝑇 ) 𝑇 𝑀 (x -x 𝑇 ) + 𝑤 𝑀 𝜽 2 𝑇 𝑀𝜽 2 + 𝑤 𝐿 𝜽 2 𝑇 𝐿𝜽 2 subject to 𝜕𝑊 𝜕𝑥 (x, 𝜽 2 ) = 0. ( 13
)
This problem can be converted to an unconstrained, nonlinear least squares program by expressing x as an implicit function of 𝜽 2 through the 𝜕𝑊 𝜕𝑥 = 0 constraint. x(𝜽 2 ) can be evaluated by solving the forward problem (Section 3.1). The energy to be minimized can then be expressed as a function of

𝜽 2 : 𝐸 (𝜽 2 ) = (x(𝜽 2 ) - x 𝑇 ) 𝑇 𝑀 (x(𝜽 2 ) -x 𝑇 ) + 𝑤 𝑀 𝜽 2 𝑇 𝑀𝜽 2 + 𝑤 𝐿 𝜽 2 𝑇 𝐿𝜽 2 .
𝐸 (𝜽 2 ) can then be differentiated using the adjoint method from sensitivity analysis as done by [START_REF] Zehnder | SGN: Sparse Gauss-Newton for Accelerated Sensitivity Analysis[END_REF]:

𝜕𝐸 𝜕𝜽 2 = -2 𝜕 2 𝑊 𝜕𝜽 2 𝜕x 𝜕 2 𝑊 𝜕x 2 -1 𝑀 (x -x 𝑇 ) + 2𝑤 𝑀 𝑀𝜽 2 + 2𝑤 𝐿 𝐿𝜽 2 . (14)
We then minimize this nonlinear least squares optimization problem using the method of [START_REF] Zehnder | SGN: Sparse Gauss-Newton for Accelerated Sensitivity Analysis[END_REF]. Optimizing the thicknesswise variations 𝜃 2 allows a significant decrease in the deviation between the target surface and the simulation (see Figure 6). Note that even though the difference between layer trajectories is small (see Figure 7), the induced deformation can be large.

7 mm 0 mm Fig. 6. Comparison with and without the optimization of 𝜃 2 . Distance ranges between 0 and 7 mm (for a 50mm-wide object).

Trajectory generation.

The final step is to extract oriented trajectories covering the plate Ŝ while following the field 𝜃 in each layer. The material will be deposited along these trajectories during fabrication. We fabricated the plate with 𝑛 different layers, so that each layer has orientations

𝜃 𝑘 = 𝜃 1 + ( 𝑘 -1 𝑛-1 -1 2 )𝜃 2 , 𝑘 ∈ [1, 𝑛].
We then use the Stripe Patterns algorithm [START_REF] Knöppel | Globally optimal direction fields[END_REF] to extract a set of equally spaced curves covering the domain and following the input direction field. This algorithm takes as input a 2-rotationally symmetric (2-RoSy) direction field in complex form, which in our case amounts to writing 𝑒 2𝑖𝜃 𝑘 . Finally, we order the trajectories by solving a Traveling Salesperson Problem (TSP), minimizing the traveling distance between curve endpoints. We use ORTools [Perron and Furnon 2023] using the Guided Local Search metaheuristic initialized by the Savings algorithm. As this NP-hard problem may take a long time to find a global minimum, we set a time limit of 5 seconds for each layer.

This set of ordered trajectories (displayed in Figure 3: Trajectory generation, with colors indicating the ordering from white to blue) can then be converted into instructions for the printer -typically G-Code -that is used to fabricate the plate.

MATERIALS STUDY

We discuss in this Section our experiments in finding the material parameters to feed to the inverse design optimization. To find the distortion ratios 𝜆 1 , 𝜆 2 and 𝜆 3 , we print 80 mm × 20 mm × 0.8 mm cuboids with parallel printing trajectories aligned with the longest sample dimension (see Figure 9). We heat each sample by immersing them in a 90 • C water bath and measure their final dimensions.1 Fig. 9. A sample before (top) and after shrinkage (bottom).

We measured the distortion ratios (𝜆 1 , 𝜆 2 , 𝜆 3 ) for several sets of printing parameters with varying extrusion temperatures and layer heights. For each set, we printed 4 samples, measured their lengths before and after shrinking, respectively 𝑙 0 and 𝑙 and computed the distortion ratio 𝜆 = 𝑙 𝑙 0 . We can equivalently (as was done in previous work) compute a contraction rate 𝑐 = 100 𝑙 -𝑙 0 𝑙 0 (in %, see Figure 8). Despite small variations, we observed that the measurements of 𝜆 1 , 𝜆 2 , 𝜆 3 were consistent across each set, standard deviations are displayed as error bars in Figure 8. We choose to use the average of a set for each ratio.

During our experiments, we noticed that the shrinking behavior depends heavily on the PLA used, with filaments from different vendors and colors showing varying shrinkage performance. We believe these differences are due to additives used in the composition of the PLA-based filament. For a given choice of filament, we observed that factors such as layer height, deposition widths, and nozzle temperature impact the shrinkage magnitude. We, however, found no dependency on the printing speeds contrary to prior works [START_REF] Kačergis | Influence of fused deposition modeling process parameters on the transformation of 4D printed morphing structures[END_REF][START_REF] Ravichandra | Additive manufacturing-enabled shape transformations via FFF 4D printing[END_REF]. This, again, may be due to differences in exact filament formulations.

For the material we used in our examples -Polymaker Polymax PLA (grey)-the best contraction rate we obtain is by using an extrusion temperature of 190 • C with a layer height of 0.08mm and a nozzle of 0.4mm diameter. We obtain a contraction rate of 41% which is higher than previous work on heat-shrinkable PLA [START_REF] Gu | Geodesy: Self-Rising 2.5D Tiles by Printing along 2D Geodesic Closed Path[END_REF][START_REF] Ravichandra | Additive manufacturing-enabled shape transformations via FFF 4D printing[END_REF][START_REF] Teunis Van Manen | Programming 2D/3D shape-shifting with hobbyist 3D printers[END_REF], even though better performance can be obtained with more exotic materials (see e.g. [START_REF] Chalissery | Highly Shrinkable Objects as Obtained from 4D Printing[END_REF]).

We use these parameters in our examples and the corresponding ratios 𝜆 1 = 0.58, 𝜆 2 = 1.08, and 𝜆 3 = 1.52.

Taking the non-zero rest curvature into account

Besides shrinking and expanding of ratios 𝜆 1 , 𝜆 2 , 𝜆 3 , the test cuboid samples also exhibited a constant, non-zero curvature after being heated. We hypothesize that this is due to re-heating during the layerby-layer fabrication. As illustrated in Figure 10, the lower layers are partially annealed by the hot nozzle re-heating them when depositing material above. This diminishes their shrinkage performance and, in turn, causes a bilayer effect: the bottom layers shrink less than the top layers, biasing the surface to curve downwards.

We, therefore, incorporated this bilayer effect into the material model, considering that 𝜆 1 varies linearly across the thickness. We modify the rest metric Equation 6as follows

C = 𝑅 𝜃 (𝜆 1 + 𝑧 ℎ Δ𝜆 1 ) 2 0 0 𝜆 2 2 𝑅 𝑇 𝜃 (15)
where Δ𝜆 1 measures the variation of 𝜆 1 according to 𝑧. We assume here that Δ𝜆 2 = Δ𝜆 3 = 0, even though 𝜆 2 and 𝜆 3 may also be subject to variations due to reheating during printing, we found their variations to be negligible. We measure the natural curvature 𝜅 taken by rectangular samples using the Kappa plugin from Fiji [START_REF] Schindelin | Fiji: an open-source platform for biological-image analysis[END_REF]] (Figure 11) and deduce the value of Δ𝜆 1 from these measurements. The samples have parallel horizontal trajectories, so for them, we have 𝜃 1 = 𝜃 2 = 0, thus

C = 𝑅 𝜃 𝜆 2 1 0 0 𝜆 2 2 𝑅 𝑇 𝜃 + 2 𝑧 ℎ Δ𝜆 1 𝜆 1 𝑅 𝜃 1 0 0 0 𝑅 𝑇 𝜃 + O (𝑧 2 ). ( 16
)
We identify the second fundamental form as

b = - Δ𝜆 1 𝜆 1 ℎ 𝑅 𝜃 1 0 0 0 𝑅 𝑇 𝜃 . ( 17
)
Making the simplifying assumption that the samples form ribbons with constant curvature 𝜅 when deforming, we find the second fundamental form is constant on the surface We only need to modify the forward method using Equation 15 to take this effect into account. The modification is naturally integrated into the inverse problem through the constraint 𝜕𝑊 𝜕𝑥 = 0 that uses the forward simulation. Note that other phenomena that could occur for different materials may be added to the method similarly by only correcting the forward simulation.

b = 𝜆 3 𝜆 2 1 𝜅 0 0 0 , (18)

Estimating the directional moduli

As described in Section 3.1.1, we observed that our material has an orthotropic behavior. We could however not measure the exact directional Young's moduli of the material. Doing proper measurements to capture the value of 𝐸 1 and 𝐸 2 is indeed challenging: one would have to set up cantilever tests inside a thermal chamber at 90 • C and then find a way to recover these coefficients from the experimental results, which is far from trivial since the test samples have non-zero rest curvature. Therefore we leave such precise measurements as future work.

Instead we experimentally set the ratio 𝐸 1 /𝐸 2 such that the simulation matches a test example printed for this purpose (Figure 12). The example, initially shaped like an ellipse squished on the poles, deploys into what looks like a pinched hemisphere because of the orthotropic behavior of the material (Figure 12 right): it is easier to bend across the trajectories than along them. The simulation using 𝐸 1 = 𝐸 2 = 1 gives a shape with a near constant curvature that does not match the fabricated object (Figure 12 left) since it does not account for the orthotropy. We obtain a much closer simulation result using 𝐸 1 = 10 and 𝐸 2 = 1 (Figure 12 middle).

RESULTS

We demonstrate our method on a variety of target surfaces, visible in Figures 1,3 and 13. We printed our examples on a Creality CR10S Pro using Polymaker Polymax PLA (grey) and immersed them in a 90 • C water contained heated using an Anova Precision sous-vide cooker. The simulation and inverse problems are computed on a 6-core Intel Xeon E-2186G processor.

The first two rows of Figure 13 show a hemisphere and a saddle having positive and negative Gaussian curvature, respectively. The Beetle example in Figure 1 reproduces the hood of a VW Beetle, which, together with the Hat and Mask model in Figure 13 showcase more complex surfaces with spatially varying curvatures. The shells we create are typical of covers designed for enclosures and appliances (e.g., a computer mouse, a hand rest, a game console, a curved phone stand) and would be very difficult to fabricate otherwise. As they start as flat plates, they all have a smooth surface finish. They are easily flexed, showing no fragility in the bending direction.

Please refer to the accompanying video for heating and deformation sequences and turntable inspection of the results.

We provide optimization performance timings and printing times for these examples in Table 1. The optimization typically takes less than one hour. The field optimization is usually the most expensive step, followed by trajectory extraction. Printing times range from one to four hours, depending on the plate size. Note that both processes require no user intervention (apart from removing the fabricated plate from the printer): the optimization runs autonomously, and the printing is made trivial by the fact that we fabricate a simple plate. Both could be optimized further to reduce the time taken.

The "Dist. to sim. " and "Dist. to optim. " columns in Table 1 measure the distance between the simulated model before/after optimizing 𝜃 2 and the target shape. As can be seen, the optimized surface is a close fit to the target with only a few percents max deviation.

Of course, the fabrication process introduces numerous approximations. First and foremost, the discretization in trajectories and layers, followed by imprecisions in the printer setup and uncertainties in the material parameters. Thus, while the results are visually similar to the target shapes, we measured the actual error between the printed result and the target shape for the beetle example. We scanned and reconstructed the fabricated shell (visible in Figure 1) using the Meshroom photogrammetry software [START_REF] Griwodz | AliceVision Meshroom: An open-source 3D reconstruction pipeline[END_REF]]. We then measured the distance between each point of the reconstructed model and the target mesh. The error is visualized in Figure 14, as well as the superposition between the target and reconstructed surfaces. We obtain an average error distance of 2.2mm and a maximal distance of 9.7mm for a model of size 129mm. While clearly imperfect, it would be impossible to design such a plate by hand that would come even remotely close to the target. For cases where a tight fit to another object is required, the natural flexibility and strength of the shell would easily allow to snap-fit it in place.

Finally, Figure 6 demonstrates the importance of optimizing the trajectory orientations across the plate thickness (𝜃 2 ), on the left the example was simulated with 𝜃 2 = 0 while the example on the right is the output of the trajectory optimization of Section 3.2.2.

LIMITATIONS AND FUTURE WORK

Developable surfaces. Our approach sometimes fails at reproducing developable surfaces. For example, we printed a near-developable fashion surface model with one piece wrapped around the body. Even though the optimization succeeded at satisfying the constraints (see how the deviations were greatly minimized, model Top in Table 1), the plate did not wrap around the correct geometry, as shown in Figure 15. In such a case, the user may manually guide the plate towards the target when it is still hot, as shown in Figure 15, right. 2This issue appears as developable surfaces are a very specific case for our approach: they undergo no change in length when deforming, while our choice of optimization strategy is based on changes in length. There is also no secondary curvature to guide Fig. 13. Our results for the models saddle, hemisphere, hat and face: (from left to right) the target mesh, our optimized printing trajectories, the resulting simulated shape and the fabricated object.

the deformation and rigidify the shell. As a result, our optimizer does not perform well in these cases.

Nevertheless, since 𝜃 2 allows us to control normal curvature, it is possible to find 𝜃 1 , 𝜃 2 that would allow to create a developable surface such as a cylinder for example. However, our optimization strategy is not well suited for this case and does not always yield the best value for 𝜃 1 . How to robustly achieve this within a more general optimization framework is a clear path for future work.

High frequency details. The thickness of the plates imposes a natural limit to the achievable normal curvature, which in turn limits the achievable details. This is properly captured in simulation (see the face model in Figure 13) such that its impact on the result is predictable. Scaling up the design for a same thickness can help capture more high-frequency details relative to the scale.

Multistability in shells.

A shell structure can exhibit multiple stable configurations but typically tends to deform towards the configuration with the lowest energy. As we run a static simulation at each iteration of our optimization, the lowest-energy state should be the one closest to the input shape -it was the case in all our results.

Table 1. For each result shown in the paper, we report the computation times separated between parameterization t param , directions optimization t optim and trajectory generation t traj , the print time, the distance between the simulated shape and the target before/after optimizing 𝜃 2 (given as percentages of the bounding box diagonal), the distance between the optimized shape and the target, the Dirichlet smoothing parameter 𝑤 𝐷 = 0.1, and the height of the initial plate in mm.

Model

Nb Two-step directions optimization algorithm. The first part of our optimizer consists in finding a flat parameterization domain Ŝ in which we can later optimize the directions (Section 3.2.1 ). A side effect of this parameterization step is we can already get an initial value for 𝜃 1 by computing the singular vector corresponding to the maximum shrinking direction. In our experiments, we found that this initial value of 𝜃 1 is sufficient so that we only optimize for 𝜃 2 in the second step (Section 3.2.2), which effectively translates into a greedy algorithm with two distinct phases. Nevertheless, it could be possible to add 𝜃 1 as optimization variables in Equation 13and therefore consider both 𝜃 1 and 𝜃 2 in the optimization, but it would come at the cost of a more challenging optimization landscape and longer computation times.

Mesh quality and refinement. Both the parameterization and the simulation objectives are sensitive to mesh quality. Therefore, we make sure to use evenly-tessellated meshes and apply Loop subdivision on them until the number of vertices is around 5000, which in our experiments strikes a good balance between computation time and mesh quality.

CONCLUDING REMARKS

We propose a method that optimizes a plate and the orientation of its fabrication trajectories such that when dipped in hot water, it deforms into a target surface. Compared to prior methods using a similar shrinkage technique, our method is the first to rely on continuous variations in trajectory orientation across the plate and its layers.

Our technique works with commercially available PLA filament with minimal calibration on standard FFF 3D printers. Thus, we expect it to be widely adopted in the maker, hobbyist, design, and education communities, where 3D printing is already widespread.

The inverse design technique provides a tight match to the target in simulation, but there are deviations between the target shapes and fabricated shell after deformation. Still, considering the many sources of approximations and uncertainty, we hope to have demonstrated the great potential of the approach, which is already useful for producing non-trivial shell surfaces.

Fabricating a flat object presents many advantages: printing is robust, there are no supports, no layering (staircase) defects. The plate can be shipped flat and deformed when needed. The forming process, performed in hot water, is simple and fun. It would be trivial to attach a layer atop the plate before deformation, such as conductive traces or colored patterns. We believe there are many open venues for future work, both in improving the technique and finding creative applications of the fabrication method.

Fig. 3 .

 3 Fig.3. Overview of our inverse design method. From an input surface, our system generates a 2D parameterization and a set of direction fields for each layer of a corresponding flat plate. Trajectories are extracted along the direction field, allowing for fabrication in PLA material on an extrusion 3D printer. Once dipped into hot water, the plate deforms into the desired surface.

Fig. 4 .

 4 Fig. 4. Ŝ represents the midsurface of the plate in the initial parametric space, S is the curved deformed midsurface. 𝜙 and 𝜓 = 𝜙 -1 are the mapping between them.

Fig. 5 .

 5 Fig. 5. Comparison of the parameterization with and without the additional smoothness term. Left: result of the local/global algorithm without any smoothness enforced, right: result with the additional smoothing step. The maximum shrinkage angle 𝛼 = 𝜃 1 is displayed here in the range [ -𝜋 /2, 𝜋 /2].

Fig. 7 .

 7 Fig. 7. Comparison between the first (left) and last (right) generated layers, with closeups on the bottom row.

4. 1 Fig. 8 .

 18 Fig.8. Experiments with varying the extrusion temperature and layer height: lower temperatures and smaller layer heights result in higher contraction rates.

Fig. 10 .

 10 Fig.10. Bilayer effect. When the nozzle extrudes material for a top layer, the layers below are re-heated locally (a). As a consequence, the bottom layers accumulate more heat than the top ones. We hypothesize that this relaxes the stresses within the bottom layers more than for the top layers (b). This translates to a linear variation of shrinkage which induces curvature in the final sample (c).

Fig. 11 .

 11 Fig. 11. Image-based measurement of the curvature 𝜅.

  therefore b = b gives us Δ𝜆 1 = -ℎ𝜆 3 𝜆 1 𝜅. ACM Trans. Graph., Vol. 42, No. 6, Article 187. Publication date: December 2023.

Fig. 12 .

 12 Fig. 12. Using the same initial state, we compare the simulation with 𝐸 1 /𝐸 2 = 1 (left), the simulation with 𝐸 1 /𝐸 2 = 10 (middle), and the fabricated example (right).

Fig. 14 .

 14 Fig. 14. Comparison between the fabricated beetle and the target mesh: (top) distance between fabricated and target visualized on the fabricated object, (bottom) superposition of the fabricated object (blue) and target surface (grey).

Fig. 15 .

 15 Fig. 15. A failure case: (left) target surface, (middle) deployed surface after heating, (right) deployed surface manually shaped when still hot.

  . faces t param // t optim // t traj (min'sec")

				Print time	Dist. to sim.	Dist. to optim.
				(h:mm)	(avg // max)	(avg // max)	w D height (mm)
	Hemisphere	17152	10.7" // 2'16" // 1'07"	2:26	3.36 // 10.28	0.25 // 3.02	0.1	150
	Saddle	2048	0.1" // 0'04" // 0'58"	1:11	2.20 // 6.05	0.47 // 2.15	0	100
	Beetle	8448	1.7" // 2'01" // 2'21"	2:49	5.41 // 16.23	0.34 // 1.41	0	200
	Face	3996	1.1" // 0'49" // 3'50"	4:18	2.81 // 12.56	0.36 // 2.04	0.1	200
	Top	3584	0.7" // 1'01" // 1'28"	2:53	5.85 // 45.62	0.24 // 2.95	0	250
	Hat	11008	16.5" // 0'57" // 1'03	1:08	1.34 // 9.31	0.30 // 1.18	0.1	100
			0	10mm

ACM Trans. Graph., Vol. 42, No. 6, Article 187. Publication date: December 2023.

ACM Trans. Graph., Vol. 42, No. 6, Article 187. Publication date: December 2023.Shrink & Morph: 3D-printed self-shaping shells actuated by a shape memory effect • 187:3

ACM Trans. Graph., Vol. 42, No. 6, Article 187. Publication date: December 2023.Shrink & Morph: 3D-printed self-shaping shells actuated by a shape memory effect • 187:5

Note that the presence of water is not necessary for the shape memory effect to occur but merely provides a suitable environment in which the temperature stays stable.

We did not guide any other result into shape, see accompanying video for the dipping process

ACKNOWLEDGMENTS

We thank the anonymous reviewers for their helpful feedback. We would also like to thank the following individuals: Etienne Vouga for insightful discussions about our shell material model, Emilie Yu for

modeling the Beetle, Hat, and Top test surfaces, and the developers of the open-source libraries which we used in our work: Geometry Central [Sharp et al. 2019a], Polyscope [Sharp et al. 2019b], TinyAD [START_REF] Schmidt | TinyAD: Automatic Differentiation in Geometry Processing Made Simple[END_REF], and libigl [START_REF] Jacobson | libigl: A simple C++ geometry processing library[END_REF]. This work was funded by the Inria Challenge "DORNELL".


Farah Abdoune

Vincent Cheutet

Maroua Nouiri

Olivier Cardin

Digital Twin for decision-support: An insight into the integration of simulation models into Digital Twin architecture

Keywords: Digital Twin, simulation, integration, decision-support, reconfiguration, architecture

The fourth industrial revolution aims to achieve greater flexibility and adaptability in manufacturing systems through the use of information and communication technologies. The Digital Twin technology has emerged as a promising solution to support human-centered decision-making in this context. Despite growing interest in this area, there is still a lack of applications that integrate decision-support functionality and emphasize the relationship between real-time Digital Twin models and what-if simulation models. Hence, this paper discusses the integration of simulation models into a Digital Twin architecture to assist operators in making appropriate decisions. A proof of concept is presented to demonstrate this approach's feasibility and open up perspectives for further research in this area.

Introduction

Seen as an important pillar of the Industry 4.0 paradigm, the Digital Twin (DT) has gained a strong following in both the scientific and industrial communities during the last few decades [START_REF] Grieves | Digital twin: manufacturing excellence through virtual factory replication[END_REF]. As it has been expanding to new sectors and applications, the DT paradigm has known a plethora of definitions, which can demonstrate a lack of maturity on this subject [START_REF] Semeraro | Digital twin paradigm: A systematic literature review[END_REF]. Although there is no one-size-fits-all definition of DT, it can be thought of as one or more simulation models of a real system that are always connected to the physical counterpart or are only briefly connected, to achieve accurate emu-lation and mirroring the operational conditions of the real systems, employing the best available representations, phys-ical and virtual models, and decisionmaking tools [START_REF] Villalonga | A decision-making framework for dynamic scheduling of cyber-physical production systems based on digital twins[END_REF]. According to NASA's literature, three key functions of the DT can be summarized [START_REF] Glaessgen | The digital twin paradigm for future nasa and us air force vehicles[END_REF]:

• Prediction is the conduct of research before the system runs.

• Monitoring and control of the system state while the system is running.

• Diagnosis is the investigation of unexpected problems that occur during the system's operation.

From a manufacturing perspective, the real-time connection allows for permanent updating, monitoring, optimizing, and enhancing industrial production system functions [START_REF] Cimino | Review of digital twin applications in manufacturing[END_REF]. These systems are often characterized by a dynamic environment where different types of events can occur. Hence, the need for adaptation and reconfiguration is predominant. To this end, the DT can be a promising technology to support the decision-making process. For instance, DT may generate "what-if" scenarios based on the state of the real system and, if necessary, send feedback to the human operators to modify the system behavior [START_REF] Agrawal | Digital twin: Where do humans fit in?[END_REF]. This implies combining DT models known generally as emulation models with offline simulation models. Following that, ideas from Boschert and Rosen [START_REF] Hehenberger | Digital twin-the simulation aspect[END_REF] and/or Schluse and Rossmann [START_REF] Schluse | From simulation to experimentable digital twins: Simulation-based development and operation of complex technical systems[END_REF] can be adopted to choose an appropriate simulation model from a model data-base coupled with collected real-time data. Alternatively, automatic simulation model generation could be used by combining semantic data models and operational data from the real system [START_REF] Kaiser | Model-based automatic generation of digital twin models for the simulation of reconfigurable manufacturing systems for timber construction[END_REF]. Although the use of DT to steer decisionmaking provides advantages, there are limitations to its implementation. [START_REF] Tao | Digital twin shop-floor: a new shop-floor paradigm towards smart manufacturing[END_REF] assert that the primary obstacle is ensuring integration, communication, and synchronicity between the physical and virtual domains, which often necessitates a technological infrastructure comprising sensors, intelligent systems, databases, and other components.

In this context, a major challenge is having a representation that is as close as possible to the real system by the simulation model to test different scenarios. This entails the integration of the relevant data in the simulation model to initialize it with the current state of the system but also provide enough information to estimate the future state of the system through the connection to the information system. While it is acknowledged that there are numerous additional challenges due to the complexity of decision-support systems, this research will solely concentrate on the challenge of integrating the simulation models into a DT architecture to test scenarios, focusing on the initialization mechanism of these models to best represent the real system.

The remainder of the paper is structured as follows: section 2 provides the state of the art, section 3 presents the methodology, and section 4 details the proof of concept. Finally, section 5 gives concluding remarks and perspectives.

State of the art

To achieve a holistic perspective about DT and decision-support, and contextualize the scope of this study, different existing DT architectures are examined in this section for starters. Among the most referenced architecture, one can find the Grieves architecture [START_REF] Grieves | Digital twin: Mitigating unpredictable, undesirable emergent behavior in complex systems[END_REF] based on 3-dimension architecture including the physical entity, the virtual model, and the connection between them. Tao et al. [START_REF] Tao | Digital twin driven prognostics and health management for complex equipment[END_REF] proposed a five-dimensional DT architecture that contains, in addition to the three original components, the DT service and data exchange between entities. According to Kritzinger et al. [START_REF] Kritzinger | Digital twin in manufacturing: A categorical literature review and classification[END_REF], another perspective based on integration level allows for a thorough grasp of how a DT should be defined. However, this positioning is not consistent with our point of view on manufacturing systems, with a clear integration of the human in the decision loop.

In a domain close to the DT paradigm, a Cyber-Physical Production Systems (CPPS) have been introduced as connected systems enhancing the decision-making process in real-time. They are systems of systems of autonomous and co-operative elements connecting in situation-dependent ways, on and across all levels of production, from processes through machines up to production and logistics networks, enhancing decision-making processes in real-time, response to unforeseen conditions, and evolution over time [START_REF] Cardin | Classification of cyber-physical production systems applications: Proposition of an analysis framework[END_REF]. From all existing CPPS architectures, the one from Cardin et al. [START_REF] Cardin | Digitalization and Control of Industrial Cyber-Physical Systems: Concepts, Technologies and Applications[END_REF] illustrates the positioning between the physical part and the cyber (i.e. digital) part, with a clear integration of the human in the decision loop for reconfiguration (Fig. 1a). In the same perspective, Bouleux et al. [START_REF] Bouleux | Requirements for a digital twin for an emergency department[END_REF] give a clear definition of what the DT represents in manufacturing (Fig. 1b) by illustrating the connection between the physical system, so-called physical twin (PT) and the DT that can integrate virtual models, simulation, optimization, etc. where the human plays a central role in decision-making process. In this context, there has been a growing interest in the development of decision-support systems based on DTs in the manufacturing sector. The following is a selection of relevant literature on the subject at hand. One approach discussed in the literature is the use of a tier-based context model with a generic ontology and model, with a focus on external factors such as the case study presented in [START_REF] Sahlab | Extending the intelligent digital twin with a context modeling service: A decision support use case[END_REF]. Another paper presents a simulation-based decision support tool that analyzes the activities occurring in distribution and production facilities by developing two simulation models [START_REF] Coelho | Simulation-based decision support tool for in-house logistics: the basis for a digital twin[END_REF]. Some papers propose the use of DTs to optimize production scheduling, rescheduling, and order management through the aggregation of several DTs representing different physical assets [START_REF] Kunath | Integrating the digital twin of the manufacturing system into a decision support system for improving the order management process[END_REF]. In [START_REF] Neto | Digital twindriven decision support system for opportunistic preventive maintenance scheduling in manufacturing[END_REF], a DT architecture is used to schedule preventive maintenance interventions, exploiting real-time opportunities and periodically accessing an updated digital model of the shop floor. Other studies propose the integration of discrete event simulation and fore-casting methods for operational planning [START_REF] Dos Santos | A decision support tool for operational planning: a digital twin using simulation and forecasting methods[END_REF]. Decision-making based on the interlinkage of DTs is achieved through the semantic ontology model in [START_REF] Meierhofer | Digital twin-enabled decision support services in industrial ecosystems[END_REF]. An architecture combining a real-time DT of logistics systems with simulation logic is discussed starting with the current state in [START_REF] Korth | Simulation-ready digital twin for realtime management of logistics systems[END_REF]. To better manage patient routes in Emergency Departments (EDs), an agent-based architecture of a DT is suggested in [START_REF] Moyaux | An agent-based architecture of the digital twin for an emergency department[END_REF]. This DT is regularly synchronized with the IS of the physical system. A comprehensive overview of DTs in the manufacturing domain is proposed in [START_REF] Pires | Decision support based on digital twin simulation: A case study[END_REF] which defines a conceptual architecture that considers simulation capabilities to optimize production processes. Another study presents a systematic literature re-view of the use of simulation as a DT to support decision-making, highlighting various uncertainties and discussions around the simulation model [START_REF] Dos Santos | Decision support in productive processes through des and abs in the digital twin era: a systematic literature review[END_REF]. While the use of simulation in DT is a relatively recent phenomenon in production processes [START_REF] Murphy | Representing financial data streams in digital simulations to support data flow design for a future digital twin[END_REF][START_REF] Karakra | Pervasive computing integrated discrete event simulation for a hospital digital twin[END_REF], other comparable approaches suggest that this concept has been investigated for several decades. For instance, Katz and Manivannan [START_REF] Katz | Exception management on a shop floor using online simulation[END_REF] and Cardin [START_REF] Cardin | Contribution of online simulation to production activity control decision support-application to a flexible manufacture system[END_REF] proposed linking a simulation model to manufacturing equipment for monitoring purposes, an approach known as "Online Simulation". It is noted that the objectives of this method are similar to those employed in DT approaches, and similar designations have been used over the years to refer to this technique [START_REF] Dos Santos | Decision support in productive processes through des and abs in the digital twin era: a systematic literature review[END_REF].

To conclude, it is observed from these papers that offline simulation models for decision-support are frequently employed to evaluate static and dynamic scenarios coupled with classic control approaches based on operational research and/or artificial intelligence for instance, or as a DT by often retrieving the necessary information for initializing the simulation model directly from the actual system. The Manufacturing Execution System (MES) as part of the information system (IS) can partially provide this information. This solution has several advantages. The collected data is reliable and represents the real state of the production system, The simulation initializes relatively quickly and it requires a minimal hardware architecture to be effective. However, a significant challenge arises because the MES only possesses partial knowledge of the system's complete state at any given time. This limitation is due to the placement and characteristics of sensors and detectors within the production system, leading to uncertainties in the state information provided by the MES. Thus, the objective of this paper is to integrate the simulation models with real-time emulation models in the DT architecture for decision support. More precisely, the study will focus mainly on the link to initialize the simulation models based on the information captured from the emulation models and the information systems.

Integration of simulation into DT architecture

This section presents first the definition of system state and model state to lay the foundation, then the integration strategy is discussed and finally, the synchronization levels are detailed.

Real System vs. virtual model

To be able to save the current state of the physical system with the emulation model and push it to the simulation model, it is important to understand some sensitive notions such as the state of the system, the state of the model, and the model's properties. A definition of these notions is suggested as follows:

• System state: The state of a system can be defined as a set of data that is used to describe the system's behavior at a specific moment and estimate the future behavior. • Model properties: Static properties of a model refer to the characteristics of the model that remain constant over time. These properties include features such as the position and speed of conveyors, the capacity of resources, and the number and types of entities involved in the process. On the other hand, the dynamic properties of a model are those that change over time. These properties include variables such as the quantity and type of products on the conveyor, the utilization of resources, and the status of entities as they move through the system.

In both emulation and simulation models, static properties remain constant, while dynamic properties are captured from the emulation model and sent to initialize the simulation model with the current state of the system. The emulation model can store this data as intermediate data models, which can be read by the simulation model to reconstruct the current state of the real system. The emulation model is designed to serve the simulation model for decision-making purposes, rather than fitting precisely to the real system for monitoring or anomaly detection purposes for instance. Therefore, the design of the emulation model is influenced by the simulation model since it is the one performing the decision-support task. Consequently, careful attention is given to capturing the necessary information that represents the system state to be transferred and read by the simulation model during the design process. For this reason, it can be seen in Fig. 2 that the information flow comes from the top to the bottom represented in blue arrows: the operator requests the simulation, the simulation requests the emulation model to capture the current state; then the data flow represented in red arrow goes from the emulation to the simulation to deliver the necessary data. Fig. 2 will be further detailed in section 3.2. Fig. 2: Integration of simulation with DT real-time models

Integration strategy

It is important to initialize the simulation model to closely match the initial state of the real system. An alternative approach is to employ an emulation model that is connected to the actual system for continuous monitoring. This model maintains a constant and seamless connection with the system, enabling prompt adjustment and synchronization of its estimations as soon as new information becomes avail-able. This approach minimizes uncertainty in the data by ensuring a continuous stream of reliable data directly from the real system.

Real-time data captured from the actual system is used to feed the emulation model, which subsequently supplies data to the simulation model. The controllers of the real system provide sensory information in discrete forms, typically represented as binary values (0s and 1s). Nevertheless, the emulation model can maintain continuous synchronization with the real system by utilizing the collected data to estimate nominal behavior. For example, when tracking the position of a pallet on a conveyor-a continuous data stream over time-the controllers report the pallet's position at each sensor point. This allows the emulation model to synchronize with the sensor data, bridging the gap between consecutive sensors. Even in instances where the real system does not provide information between two sensors, the emulation model estimates the behavior between them and synchronizes this estimation with each sensor input.

To minimize uncertainty and achieve a high level of resemblance between the simulation model and the real system, the emulation model conveys a continuous estimated value for the pallet's position on the conveyor. Consequently, the information required for initializing the simulation model may exceed that necessary for control purposes in this scenario.

The emulation model provides the current system state obtained from the controllers. Additionally, complementary data are required from the information system (MES, e.g.) to supply production data to the simulation model. This enables the simulation model to estimate the system's future state. Fig. 2 illustrates the integration of simulation models with Digital Twin (DT). The emulation models operate on a near-real-time horizon, maintaining continuous synchronization with the physical system. In contrast, the simulation model is triggered only upon operator demand, resulting in asynchronous requests

The simulation model's time horizon is referred to as "look-ahead," as it rapidly tests future scenarios. Further details regarding these two levels of synchronicity are provided in section 3.3. System data from sensors, actuators, or PLCs are consistently fed to the emulation model to maintain an up-to-date system state. Meanwhile, production data from the MES are directly supplied to the simulation model as needed. It's important to note that this process can vary depending on the application's purpose, although an attempt has been made to consider a generic scenario.

Synchronization levels

Tao et al. [START_REF] Tao | Digital twin-driven product design framework[END_REF] highlight that a crucial component of applying DT is its capacity for coevolution and synchronization between the physical and virtual realms. DT facilitates more efficient decision-making in such scenarios and is envisioned to adapt to the physical system with each improvement step, leading to a mutually beneficial development of both environments.

Two levels of synchronization can be distinguished in Fig. 2: the first one is be-tween the physical system and the emulation model and will be referred to as a "Synchronous update". In a previous study, the synchronization mechanism be-tween the physical system and the DT emulation model has been detailed and demonstrated in a case study by the authors [START_REF] Abdoune | Real-time field synchronization mechanism for digital twin manufacturing systems[END_REF]. This synchronization is done on events as follows: When an actual process occurs, it prompts the DT emulation model to generate an event. There are two scenarios: firstly, if the physical system lags behind the DT then the DT needs to wait for the physical system, and secondly, if the event occurs in the physical system before it is in the DT, then the DT should move forward to be aligned with the physical system. The second level of synchronization will be denominated as "Asynchronous request" and it represents the link between the emulation and simulation models. When triggered on demand, the synchronization is conducted by pushing all the current events happen-ing in the emulation model representing the current state of the physical system to the simulation model at once.

Proof of concept

For this work, a proof of concept is provided in this section to help the validation of the discussed integration.

Experimental setting

To demonstrate the suggested mechanism, a preliminary proof of concept has been carried out on a conveyor system that is responsible for transporting pallets. The conveyor system is equipped with two position sensors and RFID tag readers that can read the RFID tags on the pallets and is controlled by a PLC. The conveyor feeds two workstations; the operator would need to choose where to send the pal-lets depending on parameters like processing time. For the sake of this study, the physical system is represented by a model. We acknowledge that both the case study and the system used are straightforward. However, our objective is to im-plement the suggested technique, establish its viability, and illustrate its intricacy.

The emulation model is designed using the discrete simulation software Flexsim and connected to the real system via an OPC-UA server to collect data from sensors and RFIDs. Although we tend to believe that the two models are the same, this is not the case. As explained before in section 3.3, the static properties are the same; however, the dynamic properties are different. Additionally, the simulation time horizon is different between the two models. Whereas the emulation model acquires real-time data from the system and is synchronized each time it gets the data with the real system and represents only the current state of the system, the simulation model operates at a much higher speed. Fig. 3 illustrates the implementation framework that will be detailed next.

Implementation

Upon the occurrence of an external event, the operator may request a reconfiguration by defining the scenarios to be tested and sending a request to the simulation model. The simulation model requests the emulation model to deliver the necessary information. The emulation model is continuously updated and synchronized with the physical system. At this stage, the emu-lation model records all current events and generates a table to capture the current state, which is subsequently saved as a CSV file for instance. In this example, the properties of interest include the pallet number, its position on the conveyor, and its ID, which is defined as its access path in the model.

The data file containing the current system state is imported into the simulation model, which restores the system to its previous state. The simulation is then executed for a defined period, during which performance metrics, such as average time in the process and work in progress, are measured and recorded. These metrics are presented to the operator as a log of results, and he can make appropriate decisions or request further testing of alternative scenarios if deemed necessary. The focus of this case study did not involve the integration of the MES due to its complexity. Fig. 3 illustrates the implementation framework with the detailed steps highlighting five layers namely (a) the physical system (physical device along with the sensors, actuators, and PLCs), (b) the monitoring made up from the real-time connected emulation model under Flexsim software along with the synchronization mechanism, (c) the initialization which represents the detailed link between the emulation and simulation models and is made up of the intermediate data model to capture the current state, (d) the exploration which represents the test of different scenarios under simulation models and finally, (e ) the decision layer where the operator requests for reconfiguration or take the appropriate decisions.

Conclusion and perspectives

The current phase of Industry 4.0 envisions achieving greater flexibility and adaptability in manufacturing systems by progressively leveraging information and communication technologies. In this regard, the DT technology shows promise to support the human-centered decision-making process. Despite the growing interest of the scientific community in this field, there is still a dearth of applications that include decision-support functionality emphasizing the link between real-time models and what-if simulation models.

In this paper, we focused on integrating simulation models into a DT architecture to supplement real-time models and test what-if scenarios, thereby assisting operators in making appropriate decisions. We shed light on the relationship be-tween real-time and what-if simulation models by discussing various implications, such as possible integration strategies, synchronization levels, and system states. We proposed a small proof of concept to demonstrate these concepts. The paper presents new research opportunities. First, the case study can be expanded to a more complex case, such as an industrial production system, to aid operators in their daily decision-making. This will entail integrating more complex infrastructures, such as MES or ERP, to enhance the functionality of the simulation model. Future work could also concentrate on proposing a decision-support-centered DT architecture by studying decision perimeters and frequencies.

  (a) A CPPS architecture[START_REF] Cardin | Digitalization and Control of Industrial Cyber-Physical Systems: Concepts, Technologies and Applications[END_REF] (b) A DT simplified architecture[START_REF] Bouleux | Requirements for a digital twin for an emergency department[END_REF]

Fig. 1 :

 1 Fig. 1: Comparison between two architectures

Fig. 3 :

 3 Fig. 3: Implementation framework

Nantes University, École Centrale Nantes, CNRS, LS2N, UMR 6004, F-44000 Nantes, France

University Lyon, INSA Lyon, Univ. Jean Monnet Saint-Étienne, Univ. Claude Bernard Lyon 1, Univ. Lyon 2, DISP UR4570 e-mail: vincent.cheutet@insa-lyon.fr

Acknowledgement

The authors would like to thank the GDR MACS for funding the short-term mobil-ity that allowed this study to take place.