Lucas Nussbaum
email: lucas.nussbaum@loria.fr

Testbeds Support for Reproducible Research

Keywords: reference → Experimentation, • Networks → Network performance evaluation, experimentation, testbeds, reproducible research

come

INTRODUCTION

Reproducibility has been the focus of a lot of interest over the recent years, in science in general, and in computer science specifically. But most of this focus has been targeted at the reproducibility of data analysis, which is usually handled by a pipeline [START_REF] Roger | Reproducible Research[END_REF] of several steps involving various tools, starting from measured data and going up to figures and tables included in an article. The various steps of that pipeline involve code for pre-processing the measured data, analyzing data, and presenting data. And the ultimate goal for improving this pipeline is usually considered to be executable papers, where the whole process is fully automated.

However, this focus on reproducibility of data analysis ignores the important question of how that measured data is produced. Experimental computer science generally involves two main methods to acquire data about the systems under test: simulation, and experimentation on testbeds. Experimenting on a testbed is a challenging process, and usually involves many different tools: the testbed itself, usually composed of numerous services, but also experiment orchestration solutions [START_REF] Buchert | A survey of general-purpose experiment management tools for distributed systems[END_REF] ranging from shell scripts to complex frameworks, load or failure injectors, emulation solutions, measurement tools, etc. Each of those components has a huge impact on the experiment and the results that will be obtained from it.

In this paper, we provide an overview of computer science testbeds, and measure the status of three testbeds regarding reproducible research. We explore how design choices, services and features can contribute to reproducibility of experiments, and also highlight missing bits on testbeds and possible areas for improvement.

The paper is structured as follows. Section 2 provides an overview of testbeds, before Section 3 introduces in more details the three selected testbeds. Section 4 reviews how each testbed addressed various requirements for reproducible research. Some open questions are discussed in 5, before the paper is concluded in Section 6.

A SHORT PANORAMA OF TESTBEDS

Many testbeds have been designed for Computer Science, to the point that drawing a comprehensive global picture of this world is extremely challenging. Even when limiting only to the networking and distributed systems research communities, testbeds still differ in their focus (object of study): most testbeds focus on a specific field, such as IoT (e.g. IoT-Lab [START_REF] Adjih | FIT IoT-LAB: A Large Scale Open Experimental IoT Testbed[END_REF]), Internet-scale protocols and services (e.g. Planet-Lab [START_REF] Chun | Planetlab: an overlay testbed for broadcoverage services[END_REF] first, then GENI [START_REF] Mcgeer | The GENI Book[END_REF], or OFELIA [START_REF] Suñé | Design and implementation of the OFELIA FP7 facility: The European OpenFlow testbed[END_REF]), wireless networking technologies such as WiFi (e.g. R2Lab, w-iLab.t), 4G, or Software Defined Radio (e.g. CorteXLab [START_REF] Massouri | CorteXlab: An Open FPGA-based Facility for Testing SDR & Cognitive Radio Networks in a Reproducible Environment[END_REF]), clouds, HPC, or Big Data. Some testbeds are able to address several of those fields, because the underlying requirements in terms of infrastructure are similar (for example, Chameleon, CloudLab and Grid'5000 all support cloud, HPC and Big Data experiments at least up to some level, using what could be described as reconfigurable data centers).

Different objects of study raise different requirements in terms of experimental support. Testbeds differ widely in terms of resources available (ranging from wireless sensors to physical servers), of levels of access and control (from the use of virtualization technologies to bare metal reconfiguration), of monitoring features, and of guarantees on the overall environment (level of multi-tenancy on servers and network links). As data points: Planet-Lab used a container technology to slice hosts between several concurrent users, resulting in varying load over time, and the network interconnecting hosts was the real Internet, with no guarantees on performance stability (available bandwidth). On the other hand, Emulab used a cluster of identical, physical hosts, allocated exclusively to users using bare metal deployment, and network conditions are emulated by the Emulab software stack, ensuring stability over time.

Additionally, testbeds are often organized in federations, which provide either identity-level federation (enabling users to access several testbeds using the same credentials), API-level federation (generally by supporting the GENI API), or data-plane federation (to aggregate various resources from various testbeds in a single experiment, usually to increase the scale of the experiment, or to combine different resources, e.g. IoT and clouds). For example, GENI can be seen as a tightly integrated federation of testbeds (racks) throughout the United States. Fed4FIRE (and its follow-up project Fed4FIRE+) is a less tightly integrated federation (not supporting data-plane federation between all testbeds) of a more diverse set of testbeds in Europe.

EVALUATED TESTBEDS

In order to compare testbeds that are actually comparable, this paper focuses on three testbeds with a similar scope and a similar architecture: Chameleon1 , CloudLab2 , and Grid'5000 [START_REF] Balouek | Adding Virtualization Capabilities to the Grid'5000 Testbed[END_REF]. Those testbeds all welcome experiments in the fields of Cloud Computing (e.g. development and evaluation of custom cloud stacks), Big Data and High Performance Computing (mainly thanks to the availability of HPC networks like Infiniband and accelerators such as GPUs and Xeon Phi). CloudLab, and to a lesser level Grid'5000, are also suitable for experiments in wired networking, thanks to the availability of network topology emulation services. All three testbeds are similar in terms of hardware: each site is composed of racks of servers located in a data centers, and sites are interconnected by a dedicated network. All three testbeds are actively used by many users on a daily basis. Each testbed is described in the following paragraphs.

Chameleon is a recent project (started in 2014) composed of two sites, for a total of 424 nodes. Its code base is based on OpenStack, with the addition of Grid'5000 tools, and custom developments.

The CloudLab project also started in 2014 (funded by the same NSF program as Chameleon), but actually has its roots in Emulab [START_REF] White | An Integrated Experimental Environment for Distributed Systems and Networks[END_REF] 3 . Emulab is both a testbed (located in University of Utah) and a software system, developed since 2001 mainly by University of Utah. That code base is also used on numerous other Emulab instances throughout the world 4 , such as imec Virtual Wall in Belgium. CloudLab is composed of three main sites, for a total of 1081 servers, and is also federated with other Emulab instances.

Grid'5000 is a testbed located mainly in France, composed of 8 sites, for a total of 830 nodes. Most of its software stack has been developed internally since the beginning of the project in 2003.

All three testbeds provide in-vitro experimentation: experiments performed on those testbeds receive very little or no influence from the outside world. This is very different from e.g. Planet-Lab [START_REF] Chun | Planetlab: an overlay testbed for broadcoverage services[END_REF], where conditions in the real Internet have a direct impact on experiment results, and require the use of many repetitions and statistical methods to ensure sound results that do not depend on conditions at a particular time.

However, supporting reproducible research still raises a number of challenges, and requires specific features described in the next section.

SERVICES AND FEATURES 4.1 Support for Reconfiguration

The first requirement for such testbeds is to support reconfiguration by experimenters. This is required both to enable experimenters to set up an experimental environment that meets their needs, and later, to recreate the same experimental environment in order to support repeatability. All three testbeds support installing a custom software environment on reserved nodes, at the bare metal, through Frisbee [START_REF] Hibler | Fast, Scalable Disk Imaging with Frisbee[END_REF] on CloudLab, OpenStack Ironic on Chameleon, and Kadeploy [START_REF] Jeanvoine | Kadeploy3: Efficient and Scalable Operating System Provisioning for Clusters[END_REF] on Grid'5000. Even if those solutions have technical differences, they provide similar functionality.

The approaches for generating the testbed-provided system images differ: while Chameleon and CloudLab don't provide information about the process, Grid'5000 uses Kameleon [START_REF] Ruiz | Reconstructable software appliances with kameleon[END_REF]. This tool supports creating images using a set of recipes similar to configuration management tools such as Puppet, Chef or Ansible, and also uses a cache to ensure that environments generated at different times, but from the same recipes, will not differ if data downloaded from the Internet (e.g. software packages) changed since the image was first generated.

At the networking level, testbeds differ significantly. Chameleon is limited to what is supported by OpenStack Neutron, which is suitable for ensuring network isolation, but not really for creating network topologies, e.g. splitting the network using many VLANs. Grid'5000 uses KaVLAN [START_REF] Balouek | Adding Virtualization Capabilities to the Grid'5000 Testbed[END_REF] to achieve that. CloudLab provides very advanced support for networking experiments, with the ability to describe topologies together with network emulation parameters (bandwidth, latency, etc.).

Support for Collecting Provenance

Another important requirement is the ability to collect data about the experimental environment, in order to provide context for the results obtained during the experiment: which nodes were used? What was their configuration and performance? Such information is required to repeat or replicate experiments 5 (in order to ensure that the type and configuration of resources are identical), and to reproduce experiments and confirm that results are consistent.

The most basic level of support for provenance collection is documentation of available resources. On CloudLab, this is done through a textual description on web pages 67 . Additional information about resources is also available through the AM API (see Section 4.4).

Grid'5000 and Chameleon share the same technical solution [START_REF] Margery | Resources Description, Selection, Reservation and Verification on a Largescale Testbed[END_REF]: a detailed description of all resources (nodes, network equipment) in a set of JSON documents made available through a REST API. Those descriptions are automatically verified using hardware inventory tools and regression tests [START_REF] Nussbaum | Towards Trustworthy Testbeds thanks to Throughout Testing[END_REF], to avoid errors that could be introduced due e.g. to broken hardware replacements. Those descriptions are also archived, to enable one to explore and reference the past state of the testbed. Finally, those descriptions can be explored through web interfaces 8 to enable experimenters to select adequate resources.

Support for Long-term Data Storage

It is also important to provide experimenters with a solution to store the large datasets used as part of the experiment, and the artifacts generated during the experiment. All three testbeds provide various services to that effect. ChameleonCloud provides a filebased object store based on OpenStack Swift. CloudLab provides file-and block-stores, with versioning and snapshotting thanks to the use of ZFS as the underlying storage system. Grid'5000 provides a NFS-based service, and two managed Ceph clusters (file, block, or object storage). However, none of the testbeds provide support for exposing that data on the web, e.g. as companion data for articles. External data repositories would have to be used, which is probably a reasonable compromise.

Support for Automation

A common way to achieve repeatability and replicability is to rely on automation. All three testbeds provide APIs for experimenters to programmatically discover, reserve and setup resources, either using scripts or more advanced experiment orchestration tools [START_REF] Buchert | A survey of general-purpose experiment management tools for distributed systems[END_REF].

CloudLab uses the SFA AM API, an XMLRPC API which is also used on many GENI and Fed4FIRE testbeds. ChameleonCloud exposes the API provided by the various OpenStack services it uses as building blocks. Grid'5000 provides a custom-made REST API for all its services, and is in the process of adding support for the SFA AM API as part of its integration in the Fed4FIRE+ federation.

Each testbed also provides a way to automate specific experiments, by providing a framework that will automate the steps to setup commonly needed configurations, like the instantiation of an OpenStack cloud. On Cloudlab, this is done through profiles; on Chameleon, through Appliances; and on Grid'5000, through specific tools (mainly for Ceph and OpenStack deployment).

OPEN QUESTIONS

There are open questions about testbeds and reproducible research.

Respective Responsibilities of Testbeds and Experimenters

A first question is the one of how roles should be split between testbeds and experimenters. How far should testbeds go with providing advanced services to experimenters? What should be left as a burden for experimenters? This is an open question for automation of experiments, for which some testbeds provide a framework, and for monitoring: while it is expected that testbeds help experiments retrieve metrics about the infrastructure (e.g. power usage and network traffic, on Grid'5000), how far should they go with helping experimenters manage their own metrics, for example through agents running on nodes and a monitoring service?

Load Generation

Another question that is of crucial importance on the testbeds surveyed in this paper is load generation. All three testbeds provide a pristine experimental environment, with no impact from the outside world. As a result, all background traffic, perturbations, variations or faults needed in order to build realistic experiments must be generated and controlled by experimenters, using generators based on traces or probability distributions. However, surprisingly none of the testbeds really provide services to that effect, probably because of the diversity of workloads that need to be injected (ranging from background network traffic to e.g. concurrent usage in a multi-tenant cloud), and of the difficulty to obtain traces.

Standardization and Federation of Efforts

An upcoming challenge in next years is to standardize interfaces and federate testbeds, to allow the portability of experiments, and to facilitate the reproduction or even the replication of experiments outside of their original testbed. This could be achieved thanks to works on testbeds federations in the GENI and Fed4FIRE+ projects.

CONCLUSIONS AND FUTURE WORK

In this paper, we reviewed the state of three testbeds regarding support for reproducible research and discussed some open questions. In the future, this work should be extended to cover other types of testbeds, such as wireless testbeds, and more distributed testbeds. A specific goal could be to define a set of good practices for testbed operators in other to cover basic common functionality and avoid frequent mistakes. Additionally, this work mainly focused on the features advertised by testbeds: it would also be interesting to actually try to perform the same experiment on all three testbeds and draw lessons from that.

https://www.chameleoncloud.org/

https://cloudlab.us/

http://www.emulab.net

https://wiki.emulab.net/Emulab/wiki/OtherEmulabs

terminology: https://www.acm.org/publications/policies/artifact-review-badging

https://www.cloudlab.us/hardware.php

http://docs.cloudlab.us/hardware.html

https://www.chameleoncloud.org/user/discovery/


Jonàs Martínez

Haichuan Song

Orthotropic k-nearest foams for additive manufacturing

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

Fig. 1. 3D printed objects. Le : A comfy chair inspired by the work of designer Lilian van Daal

INTRODUCTION

Additive manufacturing (AM) enables the fabrication of objects having internal microstructures, with feature sizes in the order of tens of microns. The microstructures modify the large scale behavior of the lled volumes, for instance making them lighter but still rigid enough for their intended purpose. Thus, microstructures can be used to fabricate metamaterials, the physical properties of which emerge from the local geometric arrangements.

The design of microstructures raises several challenges that directly relate to computer graphics. Their ne and detailed structures, intended to ll comparatively large volumes, quickly lead to tedious modeling, large data structures, and intractable simulation costs. In addition, a major interest in using microstructures with AM is the ability to grade material properties within a volume. This a ords for materials that can adapt to local conditions, such as being lighter and porous where a part is subject to lower stress. This signi cantly complicates the design task, as the internal structures have to be speci cally tailored for each object.

In this work we propose to focus on microstructures that behave like orthotropic materials; that is, their expected elastic behavior di ers along orthogonal axes, which can be freely oriented. Orthotropic materials are found in nature: most woods, as well as bones exhibit this behavior. Orthotropy allows the material to better adapt to uneven load scenarios: they can be oriented such that their most rigid axis aligns against the largest stress vectors [START_REF] Pedersen | On Optimal Orientation of Orthotropic Materials[END_REF].

In general, orthotropic materials o er more subtle trade-o s than isotropic ones. Beyond mechanical engineering, they also allow for greater design exibility as witnessed for instance by the work 121:2 • Martínez et al.

of designer Lilian van Daal1 . Unfortunately, such metamaterials are much more challenging to synthesize and analyze, and to the best of our knowledge no works have yet focused on automatically synthesizing orientable, orthotropic graded materials that can be fabricated.

We propose a novel technique to produce microstructures that result in a controlled, freely-orientable orthotropic elastic behavior (see Figure 1). The structures are e ciently generated by a procedural function, without the need for storing or simulating their entire geometry.

Our main contributions are:

• The de nition of an orientable, orthotropic foam well suited for modeling and additive manufacturing.

• An e cient procedural evaluation of the microstructures, which is remarkably simple to implement. • The analysis of the relationship between procedure parameters and elastic properties.

We demonstrate several applications of orthotropic graded foams on 3D printed objects.

RELATED WORK

We start by discussing techniques which optimize for material properties within objects in Section 2.1. The output of these methods can be used as a control eld for the parameters of a (meta)material de ned by a microstructure. We discuss in Section 2.2 microstructures for material design, and in Section 2.3 methods that consider modeling and computational challenges of in lling 3D shapes with microstructures. Finally, we discuss in Section 2.4 approaches that optimize internal structures within objects with similar objectives of rigidity and weight reduction.

Optimizing for material properties

The eld of topology optimization is dedicated to the optimization of material properties within a domain [START_REF] Eschenauer | Topology Optimization of Continuum Structures: A Review[END_REF][START_REF] Rozvany | Aims, Scope, Methods, History and Uni ed Terminology of Computer-Aided Topology Optimization in Structural Mechanics[END_REF]]. The methods optimizing for continuously spatiallyvarying material properties (e.g. density) are of direct interest to us [START_REF] Allaire | Shape Optimization by the Homogenization Method[END_REF][START_REF] Bendsøe | Generating Optimal Topologies in Structural Design Using a Homogenization Method[END_REF]. After optimization, the question of how to map the continuous, homogenized properties to an actual material arises. Microstructures are one possible approach to this problem [START_REF] Bendsøe | Generating Optimal Topologies in Structural Design Using a Homogenization Method[END_REF][START_REF] Suzuki | A Homogenization Method for Shape and Topology Optimization[END_REF]. It is worth noting that topology optimization via homogenization optimizes for a full elasticity tensor at every point (e.g. free material optimization [START_REF] Bendsøe | An Analytical Model to Predict Optimal Material Properties in the Context of Optimal Structural Design[END_REF][START_REF] Kočvara | Free Material Optimization for Stress Constraints[END_REF], laminate optimization [START_REF] Allaire | A Numerical Algorithm for Topology and Shape Optimization[END_REF]), and thus anisotropy information is available. In practice, it is dicult to physically produce materials with this type of control, and this information is often unexploited. One exception is the work of [START_REF] Pantz | A Post-Treatment of the Homogenization Method for Shape Optimization[END_REF], which maps non-regular lattices onto the homogenized solution of a compliance minimization problem. This however requires a global optimization over the domain. Our goal is to provide a metamaterial capable of producing an orthotropic elastic behavior while being simple to conform to a control elde.g. without having to solve for a global mapping problem.

In computer graphics there is also a strong interest in fabricating deformable, animated objects, as well as interactive editing of such designs. [START_REF] Bentkus | Design and Fabrication of Materials With Desired Deformation Behavior[END_REF] design materials with prescribed deformation from a set of prede ned base materials. [START_REF] Chen | Spec2Fab: A Reducer-Tuner Model for Translating Speci cations to 3D Prints[END_REF] propose a framework for the optimization of heterogenous multimaterial arrangements achieving speci c properties. [START_REF] Skouras | Computational Design of Actuated Deformable Characters[END_REF] optimize for a bi-material distribution to achieve a desired deformation under imposed displacements. [START_REF] Xu | Interactive Material Design Using Model Reduction[END_REF] optimize for distributions of di erent isotropic materials to meet displacement and internal force constraints. [START_REF] Ion | Metamaterial Mechanisms[END_REF] provide an advanced interface to allow users to paint material properties and create mechanisms embedding microstructures.

Microstructures for (meta)material design

The main approach to make microstructure design tractable is to rely on periodic structures [START_REF] Sigmund | Materials With Prescribed Constitutive Parameters: An Inverse Homogenization Problem[END_REF][START_REF] Sigmund | Tailoring Materials With Prescribed Elastic Properties[END_REF][START_REF] Sigmund | Design of Smart Composite Materials Using Topology Optimization[END_REF]. This o ers two signi cant advantages. First, the periodicity allows storing in memory a single tile, that is implicitly repeated in a regular grid covering the object. This enables compact storage, e cient display and processing by the AM device (see Section 2.3). Second, the theory of homogenization a ords for the computation of the parameters of an equivalent homogeneous material [START_REF] Allaire | Shape Optimization by the Homogenization Method[END_REF]]. Given these, the object can be simulated while abstracting away ne scale heterogeneities.

The base tile of a periodic microstructure can be optimized through an inverse homogenization problem to target a speci c material [START_REF] Andreassen | Design of Manufacturable 3D Extremal Elastic Microstructure[END_REF][START_REF] Radman | Topology Optimization of Functionally Graded Cellular Materials[END_REF][START_REF] Sigmund | Materials With Prescribed Constitutive Parameters: An Inverse Homogenization Problem[END_REF][START_REF] Xia | Design of Materials Using Topology Optimization and Energy-Based Homogenization Approach in Matlab[END_REF][START_REF] Zhou | Design of Graded Two-Phase Microstructures for Tailored Elasticity Gradients[END_REF]. [START_REF] Schumacher | Microstructures to Control Elasticity in 3D Printing[END_REF] extended this methodology to optimize families of periodic tiles with varying properties. [START_REF] Panetta | Elastic Textures for Additive Fabrication[END_REF] proposed a family of isotropic tiles based on periodic truss structures. In both these works, tiles with di erent properties can be spatially arranged to grade properties, such as obtaining varying degrees of elasticity. This requires a special treatment of the boundaries, either ensuring the tile borders are compatible across the entire tile set [START_REF] Panetta | Elastic Textures for Additive Fabrication[END_REF], or performing a global optimization step to choose tiles with best matching borders [START_REF] Schumacher | Microstructures to Control Elasticity in 3D Printing[END_REF]. Whenever using a periodic grid care must be taken where the grid intersects the object surface, as tiles are cut by the surface [START_REF] Robbins | An E cient and Scalable Approach for Generating Topologically Optimized Cellular Structures for Additive Manufacturing[END_REF].

While extremely e cient and widely adopted, periodic microstructures have a number of disadvantages. The underlying regular grid makes it di cult to smoothly grade the structure properties along arbitrary elds. In particular, when the material is orthotropic, the grid makes it di cult to follow an orientation eld seamlessly across grid cells. [START_REF] Martínez | Procedural Voronoi Foams for Additive Manufacturing[END_REF] deviate from periodic microstructures and propose to generate aperiodic, stochastic open-cell foams with prescribed isotropic elasticity, through a procedure akin to solid procedural textures. This builds upon well known properties of open-cell foams [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], the Young's modulus of which is strongly correlated to their volume (e.g. denser foams produce sti er materials following a quasi-linear relationship). We discuss in Figure 3 possible ways to obtain an orthotropic behavior from such foams, however none proved ecient. Another possibility, not in the Figure, is to rely on anisotropic Voronoi diagrams (e.g. [START_REF] Lévy | Variational Anisotropic Surface Meshing with Voronoi Parallel Linear Enumeration[END_REF]). For instance, Nervous System uses such diagrams on surfaces to model 3D printed jewelry1 . However, extraction of curved cell boundaries (edges in 3D) of anisotropic diagrams remains computationally expensive. Disconnected islands may also appear under varying metrics. In comparison, our approach is simple to implement, procedural, and produces structures made of straight beams that exhibit a controlled orthotropic behavior.

E icient infilling with microstructures

Given a speci c microstructure geometry the question of how to ll a target shape with this detailed geometry arises. [START_REF] Wang | A Hybrid Geometric Modeling Method for Large Scale Conformal Cellular Structures[END_REF] conform truss latices in a thick shell below an object's surface through an e cient procedure that directly outputs a STL model. [START_REF] Rosen | Computer-Aided Design for Additive Manufacturing of Cellular Structures[END_REF] considers the challenges of designing parts with internal lattices and proposes a slicing procedure which works directly from the truss skeleton. [START_REF] Chen | 3D Texture Mapping for Rapid Manufacturing[END_REF] extends texture mapping and signal specialized parameterization to in ll volumes with densityvarying truss latices. [START_REF] Brennan-Craddock | The Investigation of a Method to Generate Conformal Lattice Structures for Additive Manufacturing[END_REF] studies several approaches to apply microstructures within objects, and proposes a slicing algorithm exploiting the periodicity of the structures for e ciency. [START_REF] Pasko | Procedural Function-Based Modelling of Volumetric Microstructures[END_REF] and [START_REF] Fryazinov | Multi-Scale Space-Variant FRep Cellular Structures[END_REF] explore procedural formulations of periodic microstructures. This enables e cient visualization through ray-tracing, and e cient fabrication by streaming slice images to the printer. Li et al. [2015] rely on such microstructures to optimize the internal cross-sections of an object. [START_REF] Vidimče | Foundry: Hierarchical Material Design for Multi-Material Fabrication[END_REF]2013] explore a voxel approach to de ne complex procedural structures within objects, with an emphasis on multi-materials, usability, and slicing e ciency.

Our approach answers the requirements that provide e cient procedural synthesis and direct slicing for scalability of manufacturing. Compared to lattices, we propose a signi cantly di erent viewpoint by considering stochastic structures, that do not require solving for a mapping to conform to a vector eld. Despite its randomized nature, our approach de nes a metamaterial in a principled way, and a ords for a precise control of the orthotropic behavior.

Globally optimized internal structures

A number of approaches optimize internal structures globally. [START_REF] Lu | Build-To-Last: Strength to Weight 3D Printed Objects[END_REF] optimize closed-cell foams formed by the faces of a Voronoi diagram. Wu et al. [2016a] solve for a global, high-resolution topology optimization problem under local material density constraints, which produces porous structures. Wu et al. [2016b] subdivide rhombic structures in a global optimization process, using smaller and denser cells in regions of high compliance. Other approaches optimize sparse truss structures within the volume of objects, such as to obtain rigid but lightweight 3D prints [START_REF] Wang | Cost-E ective Printing of 3D Objects With Skin-Frame Structures[END_REF][START_REF] Zhang | Medial Axis Tree-an Internal Supporting Structure for 3D Printing[END_REF]. Recent software for AM, such as Autodesk Netfabb2 , Autodesk Within3 or nTopology Element4 , provide (a) Stress in x direction,

E x = σ L ∆ x x , x = - ∆ x ∆ x x . (b) Stress in direction. E = σ L ∆ , x = - ∆ x ∆
.

the ability to model regular and randomized truss latices that can be globally optimized.

The focus of these global approaches is di erent from ours: they optimize structures answering a speci c scenario, while we seek to de ne a parameterized metamaterial. The metamaterial is only later applied within the object, following a coarse resolution control eld, without further optimization -thus metamaterials scale to arbitrarily large volumes.

BACKGROUND ON ORTHOTROPIC MATERIALS

This section recalls important notions regarding orthotropic materials. For the sake of clarity we limit the exposition to the 2D case. We provide details for the 3D case in supplemental material. All the formulas consider axis-aligned orthotropic materials. The interested reader can also refer to the text book by Jones [1975] (Chapter 2), and to the work of Li and Barbic [2015] who study orthotropic materials in the context of animation and simulation.

For linear elasticity, the stress-strain relation is σ = Cϵ, where C is the elasticity tensor, σ is the stress, and ϵ is the strain. Inversely ϵ = Sσ , where S is the compliance tensor.

Stress-strain relations. Ideal 2D orthotropic materials are de ned by ve physical properties: E x , E , x , x , G x , which are respectively the Young's moduli along x/y, the Poisson's ratios (xy/yx), and the shear modulus. Young's moduli and Poisson's ratio are illustrated in Figure 2. We recall below the expression of the compliance tensor S or tho , the expression of C or tho is found by inverting it.

S or tho = s 11 s 12 0 s 12 s 22 0 0 0 s 66 (1)

s 11 = 1 E x , s 22 = 1 E , s 12 = - x E x = - x E , s 66 = 1 G x
Note that the tensor symmetry (s 12 = s 21 ) leads to:

x E x = x E (2)
Thus, there are only four independent material properties.

PROCEDURAL ORTHOTROPIC FOAM GENERATION

The di culty to modify Voronoi foams to obtain orthotropy motivated our search for a di erent type of structure (see Figure 3 and Section 2.2).

We retain an open cell-foam design: our structures are made of connected beams. This has advantages regarding fabrication, in particular on powder fusion AM processes (e.g. Selective Laser Sintering, common in industrial applications). The open structure allows non fused material to exit the part. We generate the foams by rst constructing a graph embedded in 2D or 3D (Section 4.1), which is then thickened to obtain the beam structure. The graph does not have to be explicitly constructed and stored. Instead it is generated locally whenever needed, on-the-y (Section 4.2). The structure is in nite, aperiodic and stochastic. Most importantly, its parameters can be changed locally without globally impacting the structure. This a ords for e cient spatial grading and orientation of properties.

For clarity, in the following sections we illustrate the 2D version of our foams. Nevertheless, the description of the procedural generation extends analogously to 3D.

Graph structure

Our approach for generating orthotropic foams involves two algorithms. The rst is a point distribution algorithm, which produces the nodes of the structure. The second algorithm connects the nodes together, producing the skeleton of the beams. Each algorithm has parameters which impact the structure geometry, and ultimately its elastic behavior.

In this work we always use beams of minimal printable thickness. For the same volume this lets us maximize foam density.

Point distribution. Our structures rely on an isotropic point distribution similarly to procedural Voronoi foams [START_REF] Martínez | Procedural Voronoi Foams for Additive Manufacturing[END_REF]]. We generate one random point in each cell of a virtual grid covering the space, which provides a crude but e cient approximation of a Poisson disc distribution [START_REF] Worley | A Cellular Texture Basis Function[END_REF]]. The subdivision scheme of [START_REF] Martínez | Procedural Voronoi Foams for Additive Manufacturing[END_REF] is used to locally increase or decrease the point density.

Edge connections. In absence of anisotropy, we generate the edges of the graph by connecting each point to its k closest neighbors. For su ciently large values of k (as discussed later), this produces fully connected isotropic structures.

The edges of the graph are the essential component in inducing anisotropy in the elastic behavior of the structures. As seen in Figure 3 (top), a global stretch of an isotropic graph introduces an excellent orthotropy in its elastic response. The e ect of the stretch is to bias the angular distribution of edge lengths, making the equivalent material sti er along the direction with longest edges. Our key insight is to produce a similar bias, and hence a similar elastic response, without having to rely on a global stretch. Instead we introduce local stretches that can be easily controlled. In the limit our approach converges towards an ideal laminate (like a global stretch), which reach orthotropic elasticity bounds [START_REF] Lipton | Optimal Bounds on E ective Elastic Tensors for Orthotropic Composites[END_REF]].

We produce the angular edge length bias by stretching the distance computation when selecting the k-nearest neighbors. This encourages longer edges to appear in the direction of largest stretch, as k-nearest neighbors are selected further away (see the inset below and Figure 4). Section 5 provides an in-depth analysis of the impact of the stretch of metric on the orthotropic elastic behavior.

More precisely, let us consider two points p i ,p j in the embedded graph in dimension d. We denote the frame orientation at point p as Θ(p) : R d → [0,π ] a , with a = {1, 3} for respectively the 2D and 3D cases (Euler angles). We denote the stretch at p as H (p) : R d → R d + , which, without loss of generality is assumed to be an ordered length vector. In 2D (h u ,h ) with h u h and h = 1. In 3D (h u ,h ,h w ) with h u h h w and h w = 1. We de ne the metric tensor as [START_REF] Du | Anisotropic Centroidal Voronoi Tessellations and Their Applications[END_REF]:

M (p) = E T U E, E = R(Θ(p)), U = diag(H (p) -2 ) (3)
where R(Θ) is the n dimensional clockwise rotation matrix.

The asymmetric anisotropic distance from p i to p j is:

d p i (p i ,p j ) = (p i -p j ) T M (p i )(p i -p j ) (4)
In order to simplify the synthesis of the structure, we de ne a symmetric distance d (p i ,p j ) as:

d (p i ,p j ) = d p i (p i ,p j ) + d p j (p j ,p i ) 2 (5)
This distance is used to select the k-nearest neighbors, producing an anisotropic distribution of edges. We here conjecture that this translates to an orthotropy in the elasticity tensor. This, however, is not obvious. We verify it is the case in Section 5. We also quantify the relationship between the metric anisotropy and the impact on the Young's moduli and Poison's ratios in each direction.

Connectivity. A major concern for fabrication is to obtain fully connected structures. Studies in graph theory indicate that k-nearest random graphs are connected with very high probability for even small values of k [START_REF] Balister | Connectivity of Random k -nearest-neighbour Graphs[END_REF]]. Our case is, in fact, more favorable: when using a point distribution similar to [START_REF] Worley | A Cellular Texture Basis Function[END_REF] one can easily guarantee that at least one sample exists in each virtual grid cell. In practice we use k = 6, which with uniform density always achieves connectivity even with strong anisotropy.

We experimentally evaluate the impact of anisotropy on connectivity by producing random instances of the foam in a nite, large domain. We assume the structure to be attached to the boundary (i.e. the surface of the object), and track for any disconnected component inside. We use a square domain of unit size covered by a grid of 32 2 points, and a maximum tested stretch ratio of h u = 80. The maximum edge length was measured at 0.7. Running 7k tests on stretch ratios from 1 to 80 did not produce any disconnected component.

While we did not observe any detrimental e ect of anisotropy on connectivity, disconnections may occur when varying density. In fact, following the rst steps of the proof in [START_REF] Xue | The Number of Neighbors Needed for Connectivity of Wireless Networks[END_REF], we can easily setup a "trap" as illustrated in the inset. Such a trap consists in sandwiching a narrow band of low density in between two bands of much higher density. As a result, the nodes from the coarse band will take their k-nearest neighbors in the highdensity bands only, creating a gap in connectivity. Figure 5 reveals this behavior on our foams. Such cases seldom occur in practice as the density control elds are typically smooth ; if needed the connectivity of the structure can be checked e ciently through a The dense part is 64 times denser. As long as the low-density band is not too narrow, the "trap" situation (right most) does not occur. sweeping scheme (line in 2D, slice in 3D) as the structure synthesis is procedural. Larger values of k also alleviate this issue.

Locality. Locality is an important property for the e cient evaluation of the foam geometry: we seek for an algorithm where a sub-set of the foam can be generated without having to globally synthesize the result. To achieve this, the search radius (with the usual L 2 distance) to nd the k-nearest neighbors of a point has to be bounded.

Since we guarantee one sample per jittered grid cell, such an upper bound exists. Finding a tight upper bound is however di cult. We need to nd the smallest radius that will contain at least k points -this relates to the Gauss circle problem generalized for ellipses [START_REF] Bentkus | Design and Fabrication of Materials With Desired Deformation Behavior[END_REF][START_REF] Hardy | On the Expression of a Number as the Sum of Two Squares[END_REF]], but is not strictly the same as we rely on a jittered grid. Nevertheless, from these works we can derive a lower bound (the radius that must at least be searched) that is given in 2D by c

• h u k/(π • h u • h ), and in 3D by c • h u 3 k/(π • h u • h • h w )
, where c is the grid cell size length. In practice, we observe that the lower bound follows a trend similar to the experimental upper bound, see Figure 6. We therefore set the upper bound proportional to the lower bound, in practice using a factor of 2. As shown Figure 6 the 2D upper bound remains well above the measured edge lengths for the stretch ratios we use. During evaluation we use the bound obtained from the maximum stretch in the control eld.

Procedural evaluation of orthotropic foams

Given the properties of connectivity and locality, we devise a procedural evaluation scheme. A procedural solid texture function, which returns either empty of solid for any coordinate is possible. However Fig. 7. Foam procedural per-block evaluation. The foam is oriented according to a sine wave pa ern. Le : Two blocks of foam generated independently for di erent overlapping regions of space. The foam outside the colored rectangles is ignored a er synthesis; it only serves as padding to ensure correct computations inside. Right: Overlapping the two blocks produces a seamless result: their content is deterministic and exactly matches in the overlap region. The space spawned by the procedure is infinite.

it would be computationally wasteful due to the repeated k-nearest searches. Instead, we propose a per-block synthesis scheme that groups computations in a local region of space.

Let us assume our point distribution algorithm is pseudo-random and can generate point-sets deterministically in any box in space. This is the case of the jittered grid approach we favor. Let us consider an axis aligned box with min corner c min and max corner c max . To produce the foam geometry within the box [c min ,c max ] d , we generate the point distribution in an enlarged box [c min ,c max ] d ⊕ L max where ⊕ is the morphological dilation and L max is the bound on maximum edge length, computed from the maximal stretch in the eld. This guarantees that all the points required for the knearest search within the smaller box [c min ,c max ] d are available within the enlarged box. Thus, the set of edges matches the set that would have been generated globally. In other words, this evaluation is deterministic: generating the foam in overlapping boxes produces the same geometry within the overlapping regions. This process and its results are shown in Figure 7.

Using this approach, we can e ciently generate the foam in any spatial box, in a time that only depends on the box extent (the foam domain being in nite). In particular, this can be used for frontto-back rendering during visualization, or during slicing for AM where only a thin slab of foam has to be computed [START_REF] Pasko | Procedural Function-Based Modelling of Volumetric Microstructures[END_REF]. Multiple boxes of foam can be generated in parallel.

ANALYSIS OF ORTHOTROPIC FOAMS

For the foams to be useful as a metamaterial it is important that the equivalent material properties correlate to the synthesis parameters. If this relationship exists and can be characterized, one can easily ll a volume with a foam that produces a target elastic behavior. We use numerical homogenization (Section 5.1) to verify the link between the anisotropy of the foam geometry and the orthotropy of the elasticity tensor (Section 5.2). We then analyze the material space spawned by the foams (Section 5.3) for varying synthesis parameters. We conclude by mechanical tests on printed samples (Section 5.4). We provide additional analysis in supplemental material for homogenization versus detailed simulation and for the isotropic case.

In the remainder, we consider an isotropic base material of Young's modulus E = 1 (without loss of generality since the Young's modulus values are linearly proportional to that of the base material). We consider a base Poisson's ratio of = 0.3 which is representative of plastic materials.

Homogenization

We analyze the macroscale behavior using numerical homogenization in grids of hexahedral (3D) or quad (2D) linear elements. This computes the tensor of an equivalent homogeneous material for a given sample of our foam. We provide additional details on numerical homogenization in the supplemental material.

Similarly to [START_REF] Martínez | Procedural Voronoi Foams for Additive Manufacturing[END_REF], we generate periodic versions of the foams in a grid of varying spatial extent, and apply numerical periodic homogenization to each. We t the homogenized elasticity tensor C hom to an ideal orthotropic tensor C or tho of orientation O (see Section 3), which retrieves both the orthotropic parameters and a tting error. We increase the spatial extent until the tting error is negligible, ensuring enough foam is observed to obtain a reliable result (grid of 300 2 elements in 2D, 80 3 in 3D). For now, let us assume the orthotropic orientation O is known. We discuss its selection in Section 5.2.

To perform the tting we consider the logarithmic Euclidean distance. This has the desirable property of evaluating to the same value for both the elasticity and compliance tensors, and provides more accurate results for elasticity compared to the Frobenius norm [START_REF] Moakher | The Closest Elastic Tensor of Arbitrary Symmetry to an Elasticity Tensor of Lower Symmetry[END_REF]. The logarithmic Euclidean distance for the elasticity tensor is:

d L (C or tho ,C hom ) = log(C or tho ) -log(C hom ) F (6)
The optimization is constrained such that C or tho is positive denite [START_REF] Ting | Positive De niteness of Anisotropic Elastic Constants[END_REF]]. We use a gradient-based optimization with nonlinear constraints [START_REF] Johnson | The NLopt Nonlinear-Optimization Package[END_REF]]. Since the objective function involves matrix logarithms, its gradient is computed with numerical di erentiation.

Orientation conjecture

We generate anisotropic foam geometries by locally stretching the metric used to connect a point to its k-nearest neighbors (Section 4.1). In this section we verify that the geometric stretch translates to an orthotropy of same angle in the elasticity tensor. We perform the analysis in 2D.

To observe this, in particular on low stretches, we setup the following numerical experiment. We generate 6800 foam samples for stretch values h u in {1.5, 2.5}, where h u aligns with the horizontal axis, h with the vertical axis. We t an orthotropic tensor to each, rotating the expected angle of orthotropy O (Section 5.1), and considering the tting error. If the material is indeed orthotropic with direction O, then the tting error is expected to be small. Otherwise, the tting error is expected to increase as the directions of orthotropy do not match.

Rather than using absolute values of tting error, which are difcult to interpret, we consider the normalized measure proposed in [START_REF] Moakher | The Closest Elastic Tensor of Arbitrary Symmetry to an Elasticity Tensor of Lower Symmetry[END_REF]:

d N = d L (C or t ho ,C iso ) 2 d L (C hom ,C iso ) 2 ∈ [0 . . . 1] (7)
where C iso is a tted isotropic tensor (obtained by the tting process of Section 5.1, using an idealized isotropic tensor, see supplemental material), and C or t ho is the tted orthotropic tensor. A value of one indicates that the tted tensor is indeed close to ideally orthotropic, while zero indicates the opposite. A value above 0.9 is considered to be a good agreement [START_REF] Moakher | The Closest Elastic Tensor of Arbitrary Symmetry to an Elasticity Tensor of Lower Symmetry[END_REF]].

We perform this analysis for three di erent uniform point distributions: uniform random coordinates, Poisson disc, and the jittered grid approach we use in practice. Results are reported in Figure 8. As can be seen, for angles where O does not align with the metric stretch the error increases (lower values of d N ), while the agreement is excellent when axes align. This veri es that the stretch in the metric translates to a similarly aligned orthotropy in the material.

Another interesting observation is to consider how randomness impacts the properties (i.e. the spread of points around the average behavior). The spread is wider for a uniform random coordinate distribution, and narrower for Poisson disc, while the jittered grid lies in between. The spread of the jittered grid actually narrows on larger stretch values, and becomes similar to that of the Poisson disc distribution. This indicates that there is little di erence in quality between using a jittered grid or a Poisson disc distribution, especially with high orthotropy.

Material space

We now analyze the space of materials covered by our technique. We generate a large number of orthotropic 2D samples, and evaluate each with numerical homogenization. We vary the stretch h u (up to 15 to cover a wide range) and the density. The orientation is xed to zero since we veri ed numerically in Section 5.2 the orientation conjecture. We exploit symmetries in orthotropy when exploring the material space. The samples all have the same spatial extent but are randomized: the dataset may contain samples with similar parameters, but they are di erent random realizations.

The results for the 2D foams are summarized in Figure 9. The variations in normalized Young's moduli along each axis are shown in the top part. An isotropic material would only cover a diagonal, as the orthotropic axes would have the same Young's modulus. Our structures cover the full spectrum of Young's moduli and ratios between them. The main factor in controlling the ratio between the Young's moduli is, as expected, the metric stretch. The main factor in controlling the absolute Young's modulus values is the density.

The plot at the bottom of Figure 9 reveals the behavior of the Poisson's ratio together with the Young's moduli. Figure 9 plots only x as x follows from Equation (2). The Poisson's ratio remains positive, and sharply increases on small volumes and high stretches. Figure 10 presents the analysis of Young's moduli on 3D foam samples (cubes). As can be seen, the foams provide a good coverage of possible ratios between E x , E , E z . We limited the exploration to a stretch of at most 14 in all directions, but the space could be explored further. The Poisson's ratio behaves as with 2D foams. Please refer to the supplemental material for the plots of the 2D and 3D shear modulus, and 3D Poisson's ratio.

In summary, our foams cover the Young's modulus space densely. However, for a xed orientation of orthotropy, our foams have only two (2D) or three (3D) independently controlled elasticity parameters, so they do not o er a complete control. For instance, in 2D, choosing both Young's moduli xes the Poisson's ratio and shear modulus. Varying k changes the sti ness for a same density, but does not signi cantly change the covered material space.

The user selects a foam by specifying a desired pair (2D) or triplet (3D) of Young's modulus. The foam parameters are obtained by interpolating from the three closest samples in our dense datasets (distance weighted interpolation). If any of the three is further than a threshold derived from sampling density, the point is rejected as not covered by the foam pre-computed material space. Projecting to a closest match is left for future work. In 2D, for 13000 randomly tests we obtain an average relative Young's modulus error of 0.84%, and in 3D for 350 tests an error of 2.06%. The higher error in 3D is explained by the lower sampling of the material space.

Measurements on printed samples

We verify our numerical analysis on actual 3D printed samples. We focus this analysis on 2D versions of the foam, which are vertically extruded to obtain a 3D solid object. We print the structures using SemiFlex lament on an Ultimaker 2+ with a 0.25 mm nozzle, with a custom slicer. Each beam is exactly one lament wide. Please refer to the supplemental material for details.

We printed a family of samples with varying degrees of metric stretches. We additionally print a full sample, which acts as a reference for the (3D printed) full base material (cross-hatched, 100% in ll). The test speci cations and results are summarized in Figure 11. Figure 11a shows the raw data measured from the sensor. For tting the linear material parameters we only use the rst half millimeter of displacement (from 0 to 0.5 mm), as the material enters non-linear elasticity beyond. Note that the non-linear behavior is only a very gently bent, quasi linear curve.

Figure 11b compares the normalized predicted Young's modulus to the measured one, relative to the reference. That is, we divide each measured sample Young's modulus by the Young's modulus measured on the reference specimen. The overall behavior matches remarkably well, as the shapes of the curves closely correspond. The vertical o set between the measure and the prediction indicates that our samples are overall sti er than expected. This is most likely explained by the fabrication process, and in particular by the imperfect in lling of the fabricated reference (changes to the reference sti ness moves the curve up/down).

We provide in the graph both numerical predictions (homogenized and detailed simulations). It is expected that the homogenized material is slightly di erent since it considers an in nite periodic medium. The close agreement between them shows that the samples are su ciently large to capture the macroscale foam behavior.

We report the result of repeated compression tests in Figure 12. Compared to tensile tests, under compression the foams enter a non-linear regime after less deformation, which is due to the buckling of elongated edges. Nevertheless, all specimens exhibit a similar non-linear behavior that follows the same proportions as their Young's modulus ratios. Repeated testing of the same specimens reveal that after the rst few tests the behavior changes in the non-linear regime -while it remains very similar in the linear one. Interestingly, subsequent tests behave more stably, indicating a burn-in behavior that converges towards a stable response. This could be explained by failures in places of high local stresses, which no longer in uence the behavior after the rst few compressions. We report in Figure 13 compression tests on 3D foam specimens.

APPLICATION TO ADDITIVE MANUFACTURING

In this section we discuss the use of our foams for the design of 3D objects with orthotropic elasticity. Figure 1 (left) illustrates an orthotropic chair design inspired from the work of designer Lilian van Daal. The initial foam design is two dimensional and generated within the chair pro le. It is then extruded into the nal object, producing orthotropy in the seat and the back. Figure 1 (right) is an illustrative model of a gear made lighter by an orthotropic foam oriented against the surface, keeping it resistant to pressure. Figure 14 is a tree model with a foam oriented vertically, along the inner distance eld, showing how longitudinal strength can be preserved. Figure 10 (right) illustrates a 3D printed cube of orthotropic foam, printed on a B9 Creator resin printer. The cube is 35 × 35 × 35 millimeters with Young's moduli ratios of E z /E = 20.5,E /E x = 8.4, for stretch h x = 1.0,h = 2.0,h z = 6.0. Figure 15 show two cylinders lled with our foams and having di erent orthotropic behaviors. In each cylinder, the local stretch rigidi es one direction while making the other two more exible.

Regarding results with deformations, please keep in mind our analysis is performed on linear elasticity (small deformations), and (a) To minimize errors due to the fabrication process all specimens were printed together and measured in one session. The stretch along the compression direction is di erent for each specimen (h u = 1, 1.5, 3, 6), all other parameters are equal (h = h w = 1). Tests are performed on an Instron 3345 testing machine. The compression speed is 0.5 mm/ min, the maximum compressive strain is 1.5 % using a high precision 50 N sensor. (a) Force measured as a function of relative compressive strain. (b) Pictures of the foams (front view, initial position). (c) The first two columns report the Young's modulus (E z ) for each specimen, for simulation (relative) and experiment (MPa). For comparison the two last columns show the ratios over the average value, for each column (i.e. how much more rigid/flexible each is with respect to the average). As can been seen the ratios agree well overall. therefore we cannot predict large deformation behaviors accurately. Nevertheless, the designed structures trigger the expected global deformation. Under small deformations the analysis is accurate (see Figure 11 and 13).

h u = 1 h u = 1.5 h u = 3 h u = 6 (b)
Figure 16 illustrates how our foams can be used in conjunction with the results of topology optimization. We run the SIMP method [START_REF] Bendsøe | Optimal Shape Design as a Material Distribution Problem[END_REF]] using a penalty of p = 1 (continuous material distribution) and a volume constraint of 40%. We then extract the stress eld from the FEM simulation. Material density, principal stress direction and magnitude are mapped to the foam parameters. The extremities are isotropic, the center is rigid along the main axis but flexible radially.

While this case is essentially illustrative, we believe future integration of our foams within free material optimization methods will allow to further exploit the material space.

LIMITATIONS AND FUTURE WORK

Our approach has a number of limitations and opportunities for future work. The range of Poisson's ratio spawned by the method, while quite large, remains positive. Whether or not stochastic structures with negative Poisson's ratios can be devised remains an open question. We do not provide any method to generate beams covering the modeled objects. This is more di cult than with cellular structures, where the cells have a geometric footprint that can be subtracted from the surface [START_REF] Brennan-Craddock | The Investigation of a Method to Generate Conformal Lattice Structures for Additive Manufacturing[END_REF]. Currently, we trim or cancel the beams crossing the surface, and optionally rely on a solid shell skin. Our foams print reliably on lament printers (extruded 2D foams) and powder based systems (SLS). While we did successfully print on SLA (resin) printers two over three prints, the foams do not fully comply with the requirements. In particular, the foam graphs might exhibit unsupported local minima with respect to the build direction. The number of local minima is very small : 0.57% of the graph nodes on average (does not depend on stretch, unless reaching extreme values, e.g. h u > 100). We therefore trim these out prior to printing. The main di culty are however the elongated near horizontal beams, which tend to move with resin tanks motion.

As future work, we would like to consider other properties such as large scale deformations, fatigue and crushing/rupture behaviors. Finally, we only consider single material prints. It would be interesting to investigate a similar methodology for the design of porous structures in a multi-material setting, with beams in di erent materials.

CONCLUSIONS

We presented a stochastic, procedural, orthotropic metamaterial which can produce elastic behaviors with a wide range of Young's moduli ratio along orthogonal directions. Thanks to its stochastic nature and aperiodicity, the material properties can be graded freely in space, both in terms of local frame orientation and Young's moduli along orthogonal directions. This is achieved by mapping physical properties back to the procedural generation parameters: angles, metric stretch, and density. The foam scales trivially through procedural generation.

While our foams can be used for a wide range of design tasks and to interpret the results of topology optimization, we believe a tighter integration can be achieved with free material optimization methods [START_REF] Bendsøe | An Analytical Model to Predict Optimal Material Properties in the Context of Optimal Structural Design[END_REF][START_REF] Kočvara | Free Material Optimization for Stress Constraints[END_REF], constraining the elasticity tensor to remain in the covered material space.

Fig. 2 .

 2 Fig. 2. Uniaxial stress cases for 2D orthotropic materials.

ACMFig. 3 .

 3 Fig. 3. Inducing orthotropy from isotropic Voronoi foams. Anisotropy is tested numerically. First row: A Poisson disc distribution is used to compute a Voronoi diagram that is then stretched (×3). This structure exhibits a good orthotropic elastic behavior, but generalizing the approach would require solving for a global mapping. Second row: The Poisson disc distribution is stretched before computing the Voronoi diagram. The structure exhibits li le orthotropy. Third row: An anisotropic Poisson disc distribution is used to compute a Voronoi diagram. The structure exhibits li le orthotropy. Fourth row: Voronoi diagram of a Poisson disc distribution, where the beam radii are biased depending on their angle. The final structure exhibits only weak orthotropy (Young's moduli ratio of 1.38 for a radius factor -max over min -of 9.42).

  Fig. 4. A varying metric field (top), and the resulting k -nearest microstructure (bo om), k = 6.

Fig. 5 .

 5 Fig. 5. Le to right: Abrupt density change in a k-nearest foam (k = 6).The dense part is 64 times denser. As long as the low-density band is not too narrow, the "trap" situation (right most) does not occur.

Fig. 6 .

 6 Fig. 6. Edge lengths in a k-nearest foam (k = 6), within a ji ered grid of size 32 × 32 covering the unit square. The metric stretch goes up to an (extreme) value of h u = 80. The curves outline the lower and upper bounds.

Fig. 8 .

 8 Fig. 8. Numerical verification of the material orthotropy for uniform random coordinates, ji ered grid, and Poisson disc point distributions. The normalized measure d N is plo ed against stretch angle O (polar axis in degrees), and used as color code (bo om colorbar).

Fig. 9 .

 9 Fig. 9. Material space, please refer to the text for discussion. Top: Young's moduli along the main axes of orthotropy, color coded by the metric stretch parameter. The uneven spacing is due to the non-linearity between the sampled parameters and the elastic behavior. Middle: Foam samples. Bo om: Poisson's ratio x and Young's moduli.

Fig. 10 .

 10 Fig. 10. Covered material space for 3D foams. Le : Young's moduli along the three orthogonal axes, color coded by volume. Top right: Projection of the Young's moduli along the x and axes. Bo om right: printed sample.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Uniaxial tensile tests on a MTS ® 4/ML testing machine, with a force sensor of 100 N, displacement speed of 1 mm/ min. The specimens have a dimension of 18 × 108 mm and thickness 3.2 mm. Nine of them correspond to varying ratios of stretch h u /h , one is the reference base material. (a) Measured force/displacement curves (raw data). (b) Predicted (homogenized and detailed simulation) vs measured normalized Young's modulus. (c) Specimen with h u /h = 2 being tested.

  Fig. 13. Compression tests of 3D k-nearest foams with a beam radius of 0.25 mm within cubes of 20 × 20 × 20 mm, printed on a B9Creator (DLP resin printer).To minimize errors due to the fabrication process all specimens were printed together and measured in one session. The stretch along the compression direction is di erent for each specimen (h u = 1, 1.5, 3, 6), all other parameters are equal (h = h w = 1). Tests are performed on an Instron 3345 testing machine. The compression speed is 0.5 mm/ min, the maximum compressive strain is 1.5 % using a high precision 50 N sensor. (a) Force measured as a function of relative compressive strain. (b) Pictures of the foams (front view, initial position). (c) The first two columns report the Young's modulus (E z ) for each specimen, for simulation (relative) and experiment (MPa). For comparison the two last columns show the ratios over the average value, for each column (i.e. how much more rigid/flexible each is with respect to the average). As can been seen the ratios agree well overall.

Fig. 14 .

 14 Fig. 14. 3D printed tree with internal orthotropic foam. The foam is oriented vertically along the surface.

Fig. 15 .

 15 Fig. 15. 3D printed cylinders with orthotropic material. Top: The center segment is made more rigid horizontally and flexible vertically. Bo om:The extremities are isotropic, the center is rigid along the main axis but flexible radially.

Fig. 16 .

 16 Fig. 16. Filling a design generated by topology optimization (SIMP [Bendsøe 1989]) with our foam. Orthotropy aligns with the main direction of stress. Density is controlled from the output of topology optimization.

ACM Transactions on Graphics, Vol. 36, No. 4, Article 121. Publication date: July 2017.

Chair by designer Lilian van Daal: url. The chair was manually modeled in Rhinoceros, which took hundreds of hours, see here and here.

http://n-e-r-v-o-u-s.com/blog/?p=7465

https://www.netfabb.com/

http://www.autodesk.com/products/within

http://www.ntopology.com/

ACKNOWLEDGMENTS

We would like to thank Olivier Godard (Institut Jean Lamour, UMR 7198 CNRS -Université de Lorraine) for his help during mechanical experiments, Fabrice Neyret for valuable comments, and the anonymous reviewers for their useful suggestions. This work received support from Région Lorraine and FEDER.

This work is supported by ERC grant ShapeForge (StG-2012-307877).