Semyon Efremov

Jonàs Martínez

Sylvain Lefebvre

3D periodic cellular materials with tailored symmetry and implicit grading

Keywords: crystal symmetries, periodic and cellular structures, graded materials

Periodic cellular materials allow triggering complex elastic behaviors within the volume of a part. In this work, we study a novel type of 3D periodic cellular materials that emerge from a growth process in a lattice. The growth is parameterized by a 3D star-shaped set at each lattice point, defining the geometry that will appear around it. Individual tiles may be computed and used in a periodic lattice, or a global structure may be produced under spatial gradations, changing the parametric star-shaped set at each lattice location.

Beyond free spatial gradation, an important advantage of our approach is that elastic symmetries can easily be enforced. We show how shared symmetries between the lattice and the star-shaped set directly translate into symmetries of the periodic structures' elastic response. Thus, our approach enables restricting the symmetry of the elastic responses -monoclinic, orthorhombic, trigonal, and so on -while freely exploring a wide space of possible geometries and topologies. We make a comprehensive study of the space of symmetries and broad combinations that our method spans and demonstrate through numerical and experimental results the elastic responses triggered by our structures.

Introduction

3D periodic cellular materials are an essential ingredient in the manufacture of advanced functional parts. They allow triggering various elastic responses within a part by structuring its internal volume at a small scale. Ideally, a method to model cellular materials would allow for a wide gamut of elastic responses and geometries that could be freely chosen within the part volume.

A central line of research is the design of periodic cellular materials, defined from a representative volume element (or tile) repeated in space. The periodicity simplifies the design and analysis while limiting computational and storage requirements. We consider two types of cellular materials [START_REF] Gibson | Cellular solids: structure and properties[END_REF]: closed-cell structures which are formed by isolated void cavities (cells), and open-cell structures, formed by connected solid beams (possibly curved).

A key to achieving complex, dynamic behaviors [START_REF] Bertoldi | Flexible mechanical metamaterials[END_REF] is to allow for free spatial gradation of the structures within the part volume. Most techniques employ a catalog of tiles, pre-determined or optimized, placing different tiles with different elastic responses in different locations. Adjacent tiles will generally not connect properly, requiring applying constraints on the sets of tiles or optimizing after placement (see Section 2).

In this work, we take a distinct view on this problem. Rather than reasoning in terms of a fixed number of tiles, we start from a periodic lattice of nuclei (points) and associate each a parametric description of the local structure being produced -a 3D starshaped set. Each star-shaped set implicitly defines a distance function. Since the star-shaped set is possibly non-convex, we span a wide range of possible distance functions and a more diverse space of cellular solids compared to existing techniques. An efficient computational process grows the structures around each nucleus. Intuitively speaking, the cells centered around the lattice nuclei grow, obeying a law of growth governed by the star-shaped distance while being forbidden to overlap. Once the growth process is complete, we obtain a cellular structure, which can be extracted as either closed or open-cell. If the same starshaped distance is used at every nucleus, we obtain a periodic structure, from which a single period can act as a traditional periodic tile. If different star-shaped distances or lattices are used, we obtain spatial gradations.

Resorting to a classical Voronoi diagram under the star-shaped distance may lead to Voronoi cells with more than one connected component. Instead, the growth process ensures that each cell is connected by construction. More precisely, we consider a discrete growth process since it is challenging to formulate a continuous one: our growth process is atypical, guided by the growth of a non-convex set. To the best of our knowledge, continuous growth processes have been studied for convex distances (e.g. [START_REF] Aletti | Integration in a dynamical stochastic geometric framework[END_REF]) but not for star-shaped ones.

Our structures can be prescribed to enforce certain crystallographic geometric symmetries -and hence a given type of elastic response: monoclinic, orthorhombic, trigonal, and so on. The parametric star-shaped distances can also be easily interpolated, making it possible to transition from one geometry to another smoothly. Thus, our approach supports interpolation across different symmetries and topologies without any special treatment. It spans a broad space of possible geometries with controlled symmetries that can be freely combined and spatially graded. Finally, the compact parameterization of the distance and the lattice facilitates exploring the wide variety of possible cellular materials.

Definitions. We now define the terminology used throughout the paper. A tile, or representative volume element of a structure, is a periodic computational domain. The domain is covered by a regular grid, a set of discrete points each representing the center of a voxel (a cubic cell of the grid centered on the point). We define a nucleus as the starting point of growth of cavities or cells on the grid. We call a lattice the discrete set of nuclei distributed within the tile. The lattice has translational symmetries; therefore, it can be represented by periodic crystals [START_REF] Viterbo | Crystal lattices and crystal symmetry[END_REF]. Finally, a compact set S ⊂ R 3 is star-shaped with respect to the origin O if for all x ∈ S, the segment [O, x] is contained in S.

Related work

There is a large body of prior work studying periodic cellular materials, especially as 3D printing now enables the fabrication of small-scale internal structures [START_REF] Attene | Design, Representations, and Processing for Additive Manufacturing[END_REF]. A key challenge is to infill the interior of an object while grading physical properties. Earlier works focused on density alone, but recent approaches attempt to give finer control over the object's average elastic behavior.

Periodic trusses.

A subset of open-cell materials is represented by periodic trusses -structures composed of an organized set of struts. These structures are well-studied and find applications in multiple fields [START_REF] Ngo | Additive manufacturing (3D printing): A review of materials, methods, applications and challenges[END_REF][START_REF] Kadic | 3D metamaterials[END_REF][START_REF] Gu | Fracture of three-dimensional lattices manufactured by selective laser melting[END_REF] as they provide specific functionalities, i.e, higher strength and stiffness [START_REF] Hutchinson | The structural performance of the periodic truss[END_REF][START_REF] O'masta | The fracture toughness of octet-truss lattices[END_REF][START_REF] Song | Octettruss cellular materials for improved mechanical properties and specific energy absorption[END_REF]. The most well-known examples are, for instance, tetrahedral truss, pyramidal truss, octet-truss, and Kagome truss. Trusses are usually associated with corresponding crystal lattices. They are similarly divided by symmetries: cubic trusses correspond to cubic lattices, non-cubic trusses are obtained by an affine mapping of the cubic lattice [START_REF] Zok | Periodic truss structures[END_REF]. Although there is a range of techniques for periodic trusses, the resulting structures have a very regular nature and limit the achievable properties (geometrical or mechanical) of the material. Our method enables computation of the well-known trusses and structures with more sophisticated geometry (see Figure 1).

Single tile.

In their seminal work, Sigmund et al. [START_REF] Sigmund | Tailoring materials with prescribed elastic properties[END_REF][START_REF] Sigmund | Design of smart composite materials using topology optimization[END_REF] showed how a single tile might be optimized to achieve a specific elastic response. Here, the elastic response is to be understood in terms of average or homogenized response [START_REF] Allaire | Shape optimization by the homogenization method[END_REF]. A novel direction of research was triggered: the inverse design of a periodic tile to achieve specific mechanical properties [START_REF] Sigmund | Materials with prescribed constitutive parameters: An inverse homogenization problem[END_REF][START_REF] Zhou | Design of graded two-phase microstructures for tailored elasticity gradients[END_REF][START_REF] Radman | Topology optimization of functionally graded cellular materials[END_REF][START_REF] Xia | Design of materials using topology optimization and energy-based homogenization approach in Matlab[END_REF][START_REF] Andreassen | Design of manufacturable 3D extremal elastic microstructure[END_REF]. While most of these approaches work in grids (densities or level sets), Stankovic and Shea [START_REF] Tino | Investigation of a Voronoi diagram representation for the computational design of additively manufactured discrete lattice structures[END_REF] optimize periodic beam lattices, obtained as the edges of the Voronoi diagram of a periodic point set. The point set is parameterized by only two translation vectors, making the search space approachable. We obtain similar properties by carefully designing compact parameterizations of star-shaped distances that extend the Euclidean distance. Other techniques consider using a parameterized tile [START_REF] Fryazinov | Multi-scale space-variant FRep cellular structures[END_REF]. In such a case, the tile definition exhibits parameters -e.g., the beam's radius -that can be spatially graded. For instance, Lei et al. [START_REF] Lei | Parametric design of Voronoibased lattice porous structures[END_REF] propose to parameterize the Voronoi diagram of regular structures, varying the beam radius and the number of Voronoi seeds in order to control the degree of porosity. The construction has to guarantee that there are no mismatches from one tile to the next, preventing disconnections from appearing. Thus, the properties can be spatially graded by varying the parameters; the tile remains the same throughout the volume. This can be driven from an optimizer (e.g., [START_REF] Panesara | Strategies for functionally graded lattice structures derived using topology optimisation for additive manufacturing[END_REF][START_REF] Konaković | Beyond developable: Computational design and fabrication with auxetic materials[END_REF][START_REF] Malomo | Flexmaps: Computational design of flat flexible shells for shaping 3D objects[END_REF]) or mechanical analysis (e.g., [START_REF] Li | Interior structural optimization based on the density-variable shape modeling of 3D printed objects[END_REF]).

Families of tiles.

A core issue with using a single tile is that the gamut of achievable elastic responses is inherently restricted. To circumvent this issue, different families of tiles are considered. For instance, Bickel et al. [START_REF] Bickel | Design and fabrication of materials with desired deformation behavior[END_REF] stack a selection of periodic tile layers from a catalog, optimizing for a target elastic response. Flexible sheet bonds the layers of tiles together.

A key issue with this approach is dealing with the transition between one tile and its neighbors. Generally, the structures will not properly match one another.

One approach is to constrain the tiles always to have compatible boundaries [START_REF] Cohen | Wang tiles for image and texture generation[END_REF]. Panetta et al. [START_REF] Panetta | Elastic textures for additive fabrication[END_REF][START_REF] Panetta | Worst-case stress relief for microstructures[END_REF] optimize such families of tiles with matching boundaries, then used in a coarser scale structure optimization. Ion et al. [START_REF] Ion | Metamaterial mechanisms[END_REF] and Coulais et al. [START_REF] Coulais | Combinatorial design of textured mechanical metamaterials[END_REF] assemble mechanisms from tiles with compatible boundaries. Tozoni et al. [START_REF] Tozoni | A low-parametric rhombic microstructure family for irregular lattices[END_REF] assemble 2D tiles in conforming (deformed) 2D quad-mesh grids, accounting for distortions.

A second direction [START_REF] Zhou | Design of graded two-phase microstructures for tailored elasticity gradients[END_REF] is to enforce a good match between neighboring tiles after or during the global structure optimization. Schumacher et al. [START_REF] Schumacher | Microstructures to control elasticity in 3D printing[END_REF] place tiles from a pre-built catalog in a grid and apply a post-process to ensure a smooth, connected result. Zhu et al. [START_REF] Zhu | Two-scale topology optimization with microstructures[END_REF] perform two-scale topology optimization using a set of multi-material tiles. Boundary compatibility is taken into account during the selection of tiles. Garner et al. [START_REF] Garner | Compati-bility in microstructural optimization for additive manufacturing[END_REF] formulate a two-scale topology optimization problem, optimizing the tiles and their connectivity in the final global assembly simultaneously.

Several approaches define interpolations between tile geometries to enrich the set of tiles and favor boundary agreement [START_REF] Cramer | Microstructure interpolation for macroscopic design[END_REF][START_REF] Wang | Concurrent design with connectable graded microstructures[END_REF][START_REF] Du | Connecting microstructures for multiscale topology optimization with connectivity index constraints[END_REF][START_REF] Schumacher | Microstructures to control elasticity in 3D printing[END_REF][START_REF] Garner | Compati-bility in microstructural optimization for additive manufacturing[END_REF]. Direct interpolation is only simple when the tiles share a common topology. Otherwise, an optimization is required to refine the result [START_REF] Du | Connecting microstructures for multiscale topology optimization with connectivity index constraints[END_REF][START_REF] Garner | Compati-bility in microstructural optimization for additive manufacturing[END_REF].

There are two fundamental differences between the approaches above and our work. First, our structures are parameterized and entirely described from a compact vector encoding the star-shaped distance: there is no fixed set of tiles. Second, the structures do not require any special treatment regarding spatial gradations: the growth process locally adapts the structure, introducing variations in geometry and topology.

We build from the work of Martínez et al. [START_REF] Martínez | Star-shaped metrics for mechanical metamaterial design[END_REF] that introduces the notion of combining a star-shaped distance with an underlying lattice in 2D. We explore the much wider space of symmetries offered by the three-dimensional case and design dedicated 3D star-shaped parameterizations offering a good expressivenesscompactness trade-off. Finally, we explore both open and closedcell cellular structures, a notion that does not exist in 2D.

Method

This section introduces the two fundamental ingredients of our method: the computation of the parametric cellular geometry and the analysis of the connection between parameters, geometric symmetries, and elastic symmetries. We start by introducing in Section 3.1 our computational method to generate parameterized cellular structures with a target symmetry. Then, in Section 3.2, we examine how geometric symmetries derive into symmetries of the elastic behavior.

Parameterized cellular geometry

We compute the geometry of a periodic cellular structure with a periodic cell growth process. The growth process is compactly parameterized by a periodic lattice L ⊂ R 3 and a star-shaped set S ⊂ R 3 . Each nucleus p ∈ L is the origin of a cell. Cells evolve according to a growth law and are forbidden to overlap. Cells do not approach each other closer than a radius R > 0. The interspace left between the cells defines the solid phase of the structure. The cells are connected sets, and the inter-space left between them has at least a local diameter of R.

The growth process operates over a 3D grid G. It is algorithmically defined and guided by the star-shaped set S. A distance function d S (p, q), defined between two points p and q relying on S, governs the growth process. At the end of the growth, we obtain a periodic closed-cell geometry. If desired, we can extract an open-cell structure as a subset of the closed-cell structure.

Figure2 shows an overview of the cellular geometry computation. A significant strength of our parameterization is that combining different symmetries of L and S leads to predictable symmetries of the periodic cell geometry. As we will discuss later in Section 3.2, these geometrical symmetries directly translate to symmetries in the elastic response.

Crystal symmetries

In this work, we consider the standard classification of symmetries given by the theory of crystallography. A crystal's defining property is its symmetry: performing certain symmetry operations on the crystal leaves it unchanged. Let us recall that the following isometric transformations, denoted throughout the paper by T , are the ones exhibited by crystals [START_REF] Viterbo | Crystal lattices and crystal symmetry[END_REF]:

• Rotation over an angle of 2π n around a given axis (n-fold symmetry).

• Reflection where each point is reflected with respect to a plane. • Rotoinversion combination of a rotation and an inversion in a point on an axis. • Rotoreflection rotation about an axis and reflection in a plane perpendicular to an axis.

Thirty-two crystallographic point groups cover crystal symmetries, and only (1,2,3,4,6)-fold rotational symmetries are possible [START_REF] Monaco | Fundamentals of crystallography[END_REF][START_REF] Bradley | The mathematical theory of symmetry in solids: representation theory for point groups and space groups[END_REF]. This article always follows either the Schoenflies [START_REF] Schoenflies | Krystallsysteme und Krystallstructur[END_REF] or the Hermann-Mauguin [START_REF] Sands | Introduction to crystallography[END_REF] notation to specify point groups in three dimensions (see the Appendix B).

Lattice parameterization

We use the so-called Bravais lattices to represent L since they can be easily parameterized and represent all the crystal systems [START_REF] Viterbo | Crystal lattices and crystal symmetry[END_REF]. Moreover, Bravais lattices provide us a comprehensive and well-studied framework to analyze the link between the distance, lattice, and cellular structure symmetries. A Bravais lattice is a periodic set of points defined as

L = {i⃗ a + j ⃗ b + k⃗ c} (1)
where (i, j, k) are any integers and (⃗ a, ⃗ b, ⃗ c) are primitive vectors that lie in different directions and span the lattice. There can be additional lattice points on the center of the cell (body-centered), on the edges (a, b, c-centered), and its faces (face-centered). See Appendix C for more details on the Bravais lattices featured in this article. In addition, we used in some results a diamond cubic crystal structure.

Star-shaped set parameterization

Let us recall that a compact set S ⊂ R 3 is star-shaped with respect to the origin O if for all x ∈ S, the segment [O, x] is contained in S. We choose to parameterize the boundary ∂S with a function f S in a spherical coordinate system r, φ, θ centered in O and given by the radial distance to the boundary r > 0, the polar angle φ ∈ [0, 2π], and the azimuthal angle

θ ∈ -π 2 , π 2 0 < r min ≤ f S (φ, θ) = r ≤ r max < ∞ (2)
Please note that the function f S is continuous and single-valued.

We propose two complementary methods to parameterize Equation 2, each having specific advantages for certain symmetries. Together, these two parameterizations cover the whole space of possible crystal symmetries (see Appendix B).

1) Spherical harmonics. Spherical harmonics provide a natural basis for radial signals. They are akin to a Fourier series defined on the sphere leading to a compact form to parameterize a starshaped set analytically (e.g., see Figure 2a).

2) Spherical polyhedra. A spherical polyhedron is a partition of the sphere in which the surface is divided by arcs into bounded and convex regions called spherical polygons. We derive spherical polyhedra by projecting regular or semi-regular polyhedra vertices onto the unit sphere (see Figure 3a). Thus, we can strictly impose the symmetries T of S by selecting a polyhedral symmetry group, including the tetrahedral, octahedral, and icosahedral symmetries.

In this paper, we restrict ourselves to Platonic and Archimedean solids. We parameterize each spherical polygon as follows. Let c be the center of a convex spherical polygon P. We define a parameterized local coordinate system centered on c, that for each point of P gives the radius of S. We can impose symmetries (e.g., rotational) to the local coordinate system. Distance function. Once ∂S is compactly parameterized, we can define the distance function to guide the growth process. A distance function d S : R 3 × R 3 → R + is a mapping which represents a distance from one point to another in R 3 that satisfies the identity axiom ∀p, q ∈ R 3 d S p, q = 0 ⇔ p = q Let p, q ∈ R 3 be two points and let φ pq and θ pq be the polar and azimuthal angle given by the vector ⃗ pq = qp. The star-shaped distance from p to q is

d S (p, q) = ⃗ pq f S (φ pq , θ pq ) (3)
where ∥.∥ is the Euclidean norm. Intuitively speaking, this distance corresponds to the scaling factor, which makes ∂S centered at p to touch q (see Figure 4). Since f S is non-negative an bounded, we can verify the following set of axioms

d S p, q ≥ 0 ∀p, q (4)
d S p, q = 0 ⇔ p = q (5)

d S αp, αq = αd S p, q α ≥ 0, ∀p, q (6)
Although the distance d S does not enforce the symmetry and triangle inequality axioms of metrics, it falls into a class of distances called q 1 , q 2 -quasimetrics since it satisfies a weaker form of the triangle inequality [START_REF] Arutyunov | Theory of (q 1 , q 2 )-quasimetric spaces and coincidence points[END_REF]

d S p, q ≤ q 1 d S p, n + q 2 d S n, q (7)
where q 1 , q 2 > 0 are constant values and n is a different point from p and q. As the Euclidean metric enforces the triangle inequality axiom, following the definition of the star-shaped distance

q -p f S φ pq , θ pq ≤ q 1 n -p f S φ pn , θ pn + q 2 q -n f S φ nq , θ nq (8)
and the imposed properties of our functions (see Equation 2) a bounded upper bound of the coefficients q 1 , q 2 can be derived as follows ) such that it touches q. In that particular case d S (p, q) = 2.

q 1 = max p,n,q∈X f S φ pn , θ pn f S φ pq , θ pq = max p,n,q∈X f S φ nq , θ nq f S φ pq , θ pq = q 2 (9) (a) p, q ∈ R 3 (b) p + ∂S (c) p + d S (p, q)∂S

Algorithmic growth process

The algorithmic growth process is analogous to the 2D version proposed in [START_REF] Martínez | Star-shaped metrics for mechanical metamaterial design[END_REF]. Let us recall that the growth process operates on a 3D grid G of integer labels. First, we label each voxel containing a nucleus of the lattice L as its origin. Intuitively, the growth process operates similarly to a flood fill algorithm, using a priority queue Q ordered by the distance d S from a voxel to its emanating cell. Each voxel of G is incrementally labeled with the cell's identifier to which it belongs. We consider that the cell identifiers are greater than zero. The cells are forbidden to overlap during the growth and become closer to other cells by a radius R. Importantly, both the grid G and the growth process are periodic. At the end of the process, if a voxel is non-labeled (equal to zero), it belongs to the closed-cell structure. See Algorithm 1 for precise details.

Algorithm 1 Growth(L, S, R)

1: G = 0 ▷ Initialization 2: for l ∈ L do 3:
Let v be the voxel containing the lattice point l 4:

Insert (l, v) in Q 5: while Q is not empty do ▷ Growth 6:
Let (l, v) ∈ Q with smallest distance d S (l, center(v))

7:

Remove (l, v) from Q 8:

Let i be the index of v in G 9:

Let k > 0 be the identifier of l 10:

if

G[i] = 0 then 11:
Let v R be the voxels at distance at least R from v 12:

if v R only contains the labels 0 or k then 13:

G[i] = k 14:
for v n in 6-voxel neighbourhood of v do 15:

Insert (l, v n ) in Q return G
Open-cell structure. We can extract an open-cell structure from the closed-cell structure computed by Algorithm 1. Let D > 0 be the side length of a grid voxel and ϵ ≥ 1 be a tolerance value used to control the extraction process. We extract the voxels that belong to the closed-cell structure whose voxel neighborhood of radius R + ϵD contains at least three different cells: the voxels that lie approximatively close to the different cells' ridges. In our implementation, we choose a tolerance factor of ϵ = 1.

Geometrical symmetry parameterization

We now analyze the connection between the symmetries of the lattice L, the star-shaped set S, and the resulting symmetries of the cellular structure computed by Algorithm 1 in terms of crystallographic point groups. More point groups exist, but we only consider periodic structures associated with periodic crystals. Please note that seven types of elasticity tensor of periodic materials can be associated with the corresponding seven existing crystal systems, apart from the isotropic one [START_REF] Ting | Generalized Cowin-Mehrabadi theorems and a direct proof that the number of linear elastic symmetries is eight[END_REF].

Let T L and T S be the point groups corresponding to the lattice L and the star-shaped set S, respectively. Let us denote T L∩S = T L T S a subgroup shared by both T L and T S . All the point groups consist of only linear transformations as they are subgroups of the orthogonal group. By the definition of T L and T S , we have that ∀g ∈ T L∩S : g(S) = S and g(L) = L. Therefore, for any lattice point p ∈ L and transformation g ∈ T L∩S there is a point n ∈ L such that

g(p + S) = g(p) + g(S) = n + S ( 10
)
which hints that the discrete growth process is invariant to any transformation from T L∩S (see Figure 5) up to the voxelization of the computational domain G. Thus, to compute an approximate cellular structure with a user-given symmetry T , it suffices to ensure that both L and S share T . Growth instabilities. The growth process works very reliably for most distances. However, in some singular cases, we witness degenerate cases of the structure's geometry. This is mostly explained by growth instabilities between points close to the boundary of two or more cells. Let p, q ∈ L be two lattice points. The following implicit function characterizes the bisector points x ∈ R 3 with respect to p and q

p -x q -x = f S (φ px , θ px ) f S (φ qx , θ qx ) (11)
Note that small perturbations on p and q may significantly change x, mainly when f S varies abruptly. Thus, a discrete front of growing voxels may quickly become unstable and break the structure's expected symmetry (see Figure 6a) -we denominate these occurrences as symmetry breaking issues. We have observed, as [START_REF] Martínez | Star-shaped metrics for mechanical metamaterial design[END_REF] did, that increasing the resolution of the grid G or increasing the value of R (see Figure 6b) alleviates these symmetry-breaking events. Note also that the voxel grid captures rotational symmetries approximatively. In general, increasing the grid resolution or slightly rotating the distance attenuates this problem (see Figure 7). It also helps to avoid S, L, and G sharing additional symmetries apart from the target ones. Coincident cases. For some specific lattice and distance choices, the resulting cellular structure may have more symmetries than the ones provided by the point group T L∩S (see Figure 8 for a simple two-dimensional example). We denominate such occurrences as coincident cases. Our intuition is that coincident cases are mainly explained by a high number of shared symmetries between S, L, and G. Nevertheless, the cellular structure consistently exhibits the symmetries that correspond to T L∩S , i.e., T L∩S provides a lower bound of symmetries.

Cellular geometry gradation

The implicit parameterization offered by our method enables, in turn, the implicit gradation between different cellular structures. Let us recall that unlike most previous work (see Section 2), our method does not require any additional technique to properly transition between different cell geometries while spanning an ample space of symmetries. We distinguish two types of spatial gradation: varying the star-shaped set and varying the lattice type. Both types of geometry gradations can be used together or separately.

Star-shaped set gradation. We propose two different options. 1) Constant: S is given for each lattice point in L (see Figure 10a). 2) Linear interpolation: S is given only for a subset of lattice points in L. Let us define a set of star-shaped sets {S i } N i=1 with origin points centered at O and a set of corresponding functions { f S i x, y } N i=1 representing parametrised boundaries of these sets. Let us define the following functions

f Σ x, y = N i=1 α i f S i x, y , d Σ x, y = y -x f Σ x, y , (12)
where

0 ≤ α i ≤ 1, N i=1 α i = 1. ( 13
)
As { f S i x, y } N i=1 is the set of continuous single-valued functions, therefore, f Σ will be continuous and single-valued too and, as a result, represent the boundary of a star-shaped set with the same point of origin. Also it is non-negative and bounded

0 < r min ≤ f Σ x, y ≤ r max < ∞ ( 14
)
It can be verified that d Σ x, y is a star-shaped distance and satisfies the metric axioms stated in Equations 4 to 7. Therefore, d Σ x, y is a q 1 , q 2 -quasimetric. Let all the sets {S i } N i=1 be invariant to an operator g ∈ T S . Since g is a linear operator then

g S i N i=1 = S i N i=1 → f g(S i) x, y N i=1 = f S i x, y N i=1 (15)
Thus, a linear interpolation between an arbitrary number of parametrized star-shaped sets spawns a parametrized star-shaped set with the shared symmetry and, therefore, imposes that symmetry to the resulting structure. This property can be useful to explore the space of interpolated star-shaped sets given by S i N i=1

with user-defined symmetries. On Figure 9 we show an example of the linear interpolation between two given star-shaped sets S 1 and S 2 (see Figure 10b): for α ∈ [0, 1] we have αS 1 + (1 -α)S 2 .

Lattice gradation. We have considered the simple case of constant transition between different lattice types (see Figure 11). We leave as future work a smoother interpolation between lattices since it is non-trivial for lattices whose representative tiles have a different amount of lattice points.

Surface frame structure

The open-cell structures generated by our method can be easily conformed to a surface by considering the closed-cell structure within a distance from the surface, in a similar manner as [START_REF] Martínez | Procedural Voronoi foams for additive manufacturing[END_REF] This generates an enveloping frame structure (see Figure 12). Note that some small disconnected parts may appear in the frame.

Elastic behavior of cellular solids

This section details how we calculate the cell structures' elastic behavior and discuss the links between geometric symmetries and symmetries of the elastic behavior.

Periodic homogenization method

We perform a periodic homogenization of the cellular geometry, theoretically assuming that the cellular structure's size is infinite and that all identical periodic cells deform in the same way (see Figure 13). A commonly made assumption to homogenize the elastic behavior of cellular structures. We restrict ourselves to linear elasticity (infinitesimal strains) to give a numerical understanding of the elastic symmetries. In that case, the stress tensor σ is linearly proportional to the strain ϵ with respect to the elasticity tensor C (fourth-order symmetric tensor)

σ ⃗ x = C ⃗ x : ϵ ⃗ x (16)
Periodic homogenization seeks to determine the homogenized elasticity tensor C of a periodic structure. All the homogenized elasticity tensors shown in this article follow the Voigt notation [50]

                σ 11 σ 22 σ 33 σ 23 σ 13 σ 12                 =                  C 11 C 12 C 13 C 14 C 15 C 16 C 22 C 23 C 24 C 25 C 26 C 33 C 34 C 35 C 36 C 44 C 45 C 46 C 55 C 56 C 66                                     ϵ 11 ϵ 22 ϵ 33 2ϵ 23 2ϵ 13 2ϵ 12                    (17)
We compute the homogenized elasticity tensor C with an FFTbased method because they are well-suited for 3D grids and computationally efficient [START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF]. We have implemented the method of Monchiet et al. [START_REF] Monchiet | A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF] and Willot et al. [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF]. The main difference between both methods is that [START_REF] Monchiet | A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF] considers a continuous Green operator that is less dependent on the grid step while [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF] considers a discrete one. For most of the results, we use [START_REF] Monchiet | A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF]. When convergence becomes challenging, we resort to [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF].

For all the numerical results provided in this paper, we consider a solid phase with Young's modulus E = 100, a very soft phase (void) with E = 0.1, a homogeneous Poisson's ratio of 0.3, and a convergence threshold of ψ = 10 -4 . We provide in Appendix A all the implementation details on how we compute the homogenized elasticity tensor.

Geometric and physical symmetries

It is known that there exists a connection between the number and types of symmetries for crystal structures and the resulting symmetries of the effective elasticity tensor [START_REF] Viterbo | Crystal lattices and crystal symmetry[END_REF][START_REF] Chadwick | A new proof that the number of linear elastic symmetries is eight[END_REF][START_REF] Dvorak | Micromechanics of Composite Materials[END_REF]. In particular, the following possible types of materials differ by the amount and type of invariant transformations (see Appendix B for additional details):

• Triclinic: no particular symmetry (point group -1 or 1).

These solids are completely anisotropic. • Monoclinic: one plane of reflectional symmetry.

• Orthorhombic (or rhombic, orthotropic): three perpendicular planes of reflectional symmetry. • Trigonal and tetragonal: one axis three-fold and four-fold rotational symmetry, respectively. • Hexagonal (or transversely isotropic): one axis of 6-fold rotational symmetry and a plane of reflectional symmetry perpendicular to this axis (plane of isotropy). • Cubic: nine planes of elastic symmetry. Three of these nine planes are perpendicular to the coordinate axes, the other six have normals producing π 4 -angle with the latter. • Isotropic structures have any plane as a plane of symmetry.

Results

We now explore the variety of structures and elastic behaviors spanned by our method, focusing on symmetries. Please keep in mind that our method can seamlessly transition between all the shown structures of the same class (open or closed-cell) throughout the result sections.

Numerical results

Our method gives a tight control of the resulting periodic cellular material's symmetries. It produces all the elastic symmetries apart from strictly isotropic, a particular case in 3D (contrary to the 2D setting) that requires further optimization [START_REF] Rossi | Numerical technique for the 3D microarchitecture design of elastic composites inspired by crystal symmetries[END_REF]. Figure 21 and Figure 22 showcase galleries of close and open-cell structures pertaining to different elastic symmetry.

The galleries reveal the wide variety of symmetries, geometries, and topologies our method can produce. Nevertheless, the symmetries are precisely controlled and predicted while increasingly complex elasticity tensors are obtained. This can be seen by looking at the homogenized tensors in the last column. We recall the form of symmetric elasticity tensors in Table B.1 -note how some coefficients appear multiple times or are zero in these. Comparing the tensors with numerical results in Figure 21 and Figure 22 reveals a consistent agreement between the shared symmetries of L, S, and the elastic behavior. Occasional slight deviations are due to the voxel grid discretization. Ultimately, all the results of this article would benefit from a higher resolution of G.

It is worth noting the interplay of the distance and lattice symmetries. In some cases, the lattice has a low symmetry and the distance a high symmetry (e.g., open cell monoclinic structure in Figure 22, rhombic-dipyramidal distance, and prismatic lattice), while conversely others use a lattice with high symmetry and a distance with low symmetry (e.g., closed-cell triclinic structure in Figure 21, hexagonal-pyramidal distance, ditrigonal scalenohedral lattice). Only the shared symmetries exist in the result (see Section 3.1.5).

Also it is important to mention that changing the walls and beams' thickness would impact the overall structure's stiffness (but leave the symmetries unchanged).

We also report the timing results of Figures 21 and22 in Tables 14 and15, done with a Intel Xeon E-2186G (3.80Ghz) equipped with an NVidia GTX 1080. We consistently observe that the most time-consuming part is homogenization, and in particular, for open-cell structures.

Experimental results

We 3D printed different cellular structures obtained by our method. We use an SLS printer Lisa Sinterit and a TPU Flexa Black powder. Figure 16 shows some 3D printed samples of open-cell materials generated by our method. These are highly flexible and exhibit the expected symmetries in their directional behaviors. Figure 12b shows a 3D printed sphere filled with an open-cell material that joins gracefully with the outer frame (see Section 3.1.7). Figure 17 shows a progressive, linear interpolation between two closed-cell structures. Note how our method smoothly interpolates between these complex geometries. Figure 18 shows a concentric gradation, where a different cellular geometry is used in the center. Our method correctly interpolates through this constant implicit gradation, which is usually challenging since the change of geometry is not aligned with the underlying lattice of periodic tiles. Figure 19 shows an elongated box with different open-cell materials and constant change in the control field. The different elastic properties give rise to interesting non-uniform deformations. As long as the structures use a same lattice, progressive transitions are achievable (see Figure 17).

Figure 20 shows experimental results of compression tests. We obtain a good agreement between the predicted and measured values in terms of the relative Young's moduli of different directions.

              (d) Elasticity tensor

Conclusions

We have introduced a method for generating cellular materials with a compact and straightforward parameterization covering a broad spectrum of elastic symmetries. We provide a first study that attempts to shed some light on the materials given by our parameterization. Notably, our method can generate periodic tiles or spatial gradations implicitly without requiring to deal with spatially-varying topological or geometrical changes. We have numerically verified the expected symmetry of the results and provided a comprehensive enumeration of various means of achieving a given symmetry. Figure 20: Comparison of the experimental and numerical relative Young's moduli of three different specimens. The experimental Young's modulus in each orthogonal direction x, y, and z was measured with an Instron 3345 machine. The direction z is the printing the direction. First and second column: for each specimen (triclinic, monoclinic, cubic), we divide the three values of Young's moduli (experimental, numerical) with its average to obtain the relative value. The third column shows the division between the corresponding values of the first and second columns ; values close to one indicate that the relative moduli correspond.

Limitations and future work. Our method can have wide applications in the design and optimization of cellular materials and opens a number of topics for future work. One of the most relevant is the parametric optimization of periodic materials through our compact parameterization, for example, to obtain 3D isotropic auxetic materials [START_REF] Andreassen | Design of manufacturable 3D extremal elastic microstructure[END_REF]. Large-scale computational exploration of our parameterized materials could shed light on the achievable gamut of elastic behaviors: our current implementation of the growth process is relatively slow and could be optimized with parallelization. Bringing our method to large-scale grids is not trivial and requires further study. We could improve computational scalability by precomputing different periodic tiles and restricting the growth process computation at the interface between them.

Besides, it would be interesting to study more complex lattices (e.g., quasicrystalline). Finally, it would be useful to investigate further how to improve the robustness of Algorithm 1 to deal with symmetry-breaking events.

Aknowledgements

This work was partly supported by the ANR MuFFin (ANR-17-CE10-0002) and "Région Lorraine and FEDER". We thank Pierre-Alexandre Hugron and Pierre Bedell for providing their help with 3D printing and making the photos of the results.

Star-shaped set S Lattice L T L∩S Periodic cell 3 × 3 Tensor              
              6/mmm 1 1 Triclinic              
21.742 6.144 6.525 -0.001 -2.000 -0.001 17.346 6.020 -0.002 -0.312 -0.001 23.814 -0.003 -3.912 0.000 5.066 0.000 -1.421 6.838 0.000 6.012

              m 4/mmm m Monoclinic               16
              1 m 3m 1 Triclinic              
3.393 2.276 1.318 0.000 -0.821 0.000 6.828 2.913 0.000 -1.388 0.000 4.128 0.000 -1.181 0.000 3.041 0.000 -1.385 1.277 0.000 1.569

              mmm 2/m 2/m Monoclinic              
0.447 0.291 0.050 0.000 0.000 0.000 0.450 0.107 0.000 0.000 0.000 1.758 0.000 0.000 0.000 0.127 0.000 0.000 0.115 0.000 0.101

              222 6/mmm 222 Orthorhombic              
1.362 0.500 0.680 0.040 0.002 -0.004 1.367 0.656 -0.038 0.004 -0.015 3.079 -0.003 0.001 -0.003 0.136 0.000 0.003 0.134 0.039 0.453

              6 3m 3 Trigonal              
6.710 4.235 2.925 0.000 0.000 -0.003 6.710 2.925 -0.001 0.000 -0.003 4.115 0.000 0.000 0.000 1.913 0.123 0.000 1.913 0.000 3.0314

              422 4/mmm 422 Tetragonal              
0.859 0.612 0.105 0.000 0.000 0.000 0.858 0.100 0.000 0.000 0.000 1.474 0.000 0.000 0.000 0.070 0.000 0.000 0.068 0.000 0.119 It is known that methods with continuous Green's operator may have convergence problems when dealing with porous materials with a high contrast ratio between the solid and soft phases (i.e., inclusion and matrix) [START_REF] Schneider | Computational homogenization of elasticity on a staggered grid[END_REF][START_REF] Herve | Convergence of iterative methods based on neumann series for composite materials: Theory and practice[END_REF]. In contrast, the discrete Green's operator tends to provide better convergence even for arbitrarily large contrast between the phases. In this paper, we prioritized Monchiet et al. [START_REF] Monchiet | A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF] method since the Green's operator's continuous representation does not depend on the grid step and is suitable for grids with low resolution. We used Willot et al. [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF] for high contrast scenarios or if there was a need to compute a smoother displacement field (e.g., Figure 13). Both implementations were compared with the publicly available software CraFT (version 1.0) 1 that implements the method of Moulinec et al. [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]. For a set of sample structures, the comparisons gave similar results compared to CraFT. Below we expose a brief description of [START_REF] Monchiet | A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF][START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF]. For more information, the reader is welcomed to refer to the papers directly.

              622 6/mmm 622 Hexagonal               14
Hooke's law represents the elastic deformation problem in combination with equilibrium and strain-displacement equations

σ ⃗ x = C ⃗ x : ϵ ⃗ x , divσ ⃗ x = 0, ϵ ⃗ x = 1 2 ∇u ⃗ x + ∇ T u ⃗ x .
As it is more closely described in [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF][START_REF] Monchiet | A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF], the linear Hooke's law boils down to the Lippmann-Schwinger equation with the average media elasticity tensor C 0 , leading to Green's operator corresponding to the latter

ϵ ⃗ x + Γ 0 ⃗ x ⋆ δC ⃗ x : ϵ ⃗ x = 0.
Appendix A.1. Continuous Green's operator iterative scheme The homogenization algorithm of [START_REF] Monchiet | A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF] is based on a spectral decomposition of the Green's tensor on the following basis:

E 1 ⃗ ξ = 1 2 ⃗ k ⊥ ⊗ ⃗ k ⊥ , E 4 ⃗ ξ = ⃗ k ⊥ ⊗ ⃗ k + ⃗ k⊗ ⃗ k ⊥ , E 2 ⃗ ξ = ⃗ k ⊗ ⃗ k, E 5 ⃗ ξ = ⃗ k ⊗ ⃗ k ⊥ , E 3 ⃗ ξ = ⃗ k ⊥ ⊗ ⃗ k ⊥ -E 1 ⃗ ξ , E 6 ⃗ ξ = ⃗ k ⊥ ⊗ ⃗ k.
Here

⃗ k ⃗ ξ = 1 |ξ| 2 ξ i ⊗ ξ j , ⃗ k ⊥ ⃗ ξ = ⃗ i -⃗ k ⃗ ξ , ⃗ i is the identity tensor of
the second order and ⃗ ξ = ξ 1 , ξ 2 , ξ 3 is a frequency vector. Then the tensor projectors

Q ⃗ ξ = E 1 ⃗ ξ + E 3 ⃗ ξ , P ⃗ ξ = E 2 ⃗ ξ + E 4 ⃗ ξ , Q ⃗ ξ : ϵ ⃗ ξ = 0, P ⃗ ξ : σ ξ = 0
can be used for a convergence test:

P ⃗ ξ : σ ξ 2 σ ⃗ ξ 2 < ψ.
Then Green's tensor associated with the reference medium stiffness tensor C 0 = 3k 0 J + 2µ 0 K with J = 1 3 ⃗ i ⊗ ⃗ i and K = I -J ( ⃗ i and I are the second-order and the fourth-order identity tensors, respectively) can be defined by the following expression:

Γ 0 ⃗ ξ = 1 λ 0 + 2µ 0 E 2 ⃗ ξ + 1 2µ 0 E 4 ⃗ ξ ,
Here parameters λ 0 and µ 0 represent the Lamé coefficients of the reference media. In this work, they were computed as

λ 0 = 1 2 (λ 1 + λ 2 ) , µ 0 = 1 2 µ 1 + µ 2 ,
1 https://lma-software-craft.cnrs.fr/

where the indexes 1, 2 correspond to the matrix and the inclusion materials. Finally, the iterative method is presented on the following scheme

ε i ⃗ x = F -1 ε i ⃗ ξ , σ i ⃗ x = C ⃗ x : ε i ⃗ x , σ i ⃗ ξ = F σ i ⃗ x ,
convergence test,

ε i+1 ⃗ ξ = ε i ⃗ ξ -Γ 0 ⃗ ξ : σ i ⃗ ξ ,
with the initial strain tensor value In this table we consider two standard nomenclatures for crystal point groups (Hermann-Mauguin [START_REF] Sands | Introduction to crystallography[END_REF] and Schoenflies [START_REF] Schoenflies | Krystallsysteme und Krystallstructur[END_REF]. For a comprehensive introduction to 3D elasticity tensor symmetries see [START_REF] Dvorak | Micromechanics of Composite Materials[END_REF]. Please note that for monoclinic, tetragonal, trigonal, and hexagonal crystal systems the elasiticy tensor may vary depending on the axis of symmetry (in the table we only show one of such possibilities). Let m be the number of independent constants of the elasticity tensor. Let Y k l be the spherical harmonic function of degree l and order k.

ε 0 ⃗ x = 0 ∀ ⃗ ξ 0, E ⃗ ξ = 0.

Figure 1 :Figure 2 :

 12 Figure 1: Our approach is able to compute some well-known periodic trusses and a vast space of open-cell structures with enhanced geometrical and mechanical diversity. Here we show some classic examples of periodic open-cell structures computed with our growth method using the Euclidean distance (in gray color) and other open-cell structures using the same underlying lattice with tailored symmetries (in green color)

Figure 3 :

 3 Figure 3: Parameterization with spherical polyhedra. (a) is regular spherical polyhedron, that has six identical faces and (b) is the corresponding star-shaped set with six identical spikes, each corresponding to a face of the spherical polyhedron.

Figure 4 :

 4 Figure 4: Illustration of 3D the star-shaped distance. (a) Two points p and q are linked by a line segment. (b) Boundary set ∂S with origin at p. (c) Scaling of ∂S by a factor of d S (p, q) such that it touches q. In that particular case d S (p, q) = 2.

Figure 5 :

 5 Figure 5: Some star-shaped sets are centered in a lattice, both sharing a hexagonal symmetry (6/mmm (D 6h ) for the lattice and 6 (C 6 ) for the distance). Here we can visually appreciate how a rotation of angle factor of π 6 around the origin of any star-shaped set keeps the periodic arrangement invariant (i.e., it leads to a 6 (C 6 ) point group).

Figure 6 :

 6 Figure 6: A simple example of growth instabilities. In this case, the periodic lattice L is a single point in the cube center. The star-shaped distance has orthorhombicdipyramidal symmetries (3 orthogonal reflection planes). (a) Due to the discrete size of G, we obtain a result with symmetry-breaking artifacts (red regions). (b) Increasing the radius helps to overcome these events, obtaining the expected periodic open-cell (c).

Figure 7 :

 7 Figure 7: An example of slight rotational symmetry-breaking event. (a) A hexagonal distance leads to a closed-cell structure with a slight deviation from the expected 6-fold symmetry. (b) A slight rotation of the distance leads to a structure with better 6-fold symmetry.

Figure 8 :

 8 Figure 8: Illustrative coincident case in two-dimensions. The star-shaped set S has two reflectional planes. The lattice has four-fold symmetry. The resulting cellular structure has four-fold symmetry.

S 1 S 1 +S 2 2 S 2 Figure 9 :

 11229 Figure 9: Linear interpolation between S 1 and S 2 . S 1 is a set with monoclinic domatic symmetry, while S 2 has hexagonal dipyramidal symmetry. Both sets share a monoclinic symmetry exhibited by the interpolating set S 1 +S 2 2 .

Figure 10 :Figure 11 :

 1011 Figure 10: Closed-cell geometry gradation: central part with S 2 and the rest with S 1 .

Figure 12 :

 12 Figure 12: Different surface frame structures, geometrically smoothed. A rendered structure bounded by a sphere (a) and a printed one (b). A 8x8x8 structure (c), which is fitted in a cylinder (d) and a torus (e).

  0.188 0.145 -0.002 -0.003 -0.022 0.359 0.155 -0.001 0 -0.015 0.301 0.007 0.009 0.030 0.136 0.006 -0.003 0.137 0.004 0.184

Figure 13 :Figure 14 :

 1314 Figure 13: Given a periodic structure (a) we compute its homogenized elasticity tensor (d) from the prescribed strain and stress fields over G. From these, fields we can also visualize, for instance, the overall magnitude of displacement fields (b), or the Von mises stress the analyze the yielding of the structure (c).

Figure 15 :

 15 Figure 15: Computational timings for the open-cell structures in Figure 22

Figure 16 :

 16 Figure 16: 3D printed periodic cellular materials.

Figure 17 :

 17 Figure 17: 3D printed, progressive gradation between two closed-cell materials.

Figure 18 :

 18 Figure 18: 3D printed concentric grading, with a sharp transition.

Figure 19 :

 19 Figure 19: 3D printed elongated with spatially varying structures.

Figure 21 :

 21 Figure 21: Gallery closed-cell

Figure 22 :

 22 Figure 22: Gallery open-cell

  21.159 7.115 7.052 -0.067 1.7572 -1.371 21.911 6.730 1.070 0.477 -1.579

	19.245 0.110 0.762 1.262 7.361 0.535 0.726 7.382 0.052 8.275

  C 11 C 12 C 13 C 14 C 15 C 16 C 22 C 23 C 24 C 25 C 26 C 33 C 34 C 35 C 36 C 44 C 45 C 46 C 55 C 56 C 11 C 12 C 13 0 C 15 0 C 22 C 23 0 C 25 0 C 33 0 C 35 0 C 44 0 C 46 Spherical hexahedron / octahedron / cuboctahedron Table B.1: Crystal systems, associated point groups, and its the corresponding definition of the star-shaped set S (either using spherical harmonics or spherical polyhedra).

						
	Triclinic	Pedial	1	C 1	1 + 0.8 Y 1 2 φ + 0.5π, θ + Y 1 3 φ, θ + Y 2 2 φ + 0.35π, θ	               	C 66
		Domatic Prismatic	m 2/m	C 1h C 2h	1 + Y 1 2 φ, θ + Y 2 3 φ -0.25π, θ + Y 4 3 φ, θ 1 + Y 1 4 φ, θ + Y 4 6 φ, θ	    	C 55 0 C 66	                	13
						Symmetry plane 2
	Orthorhombic	Rhombic-disphenoidal	222	D 2	1 + Y 2 3 φ -0.25π, θ + Y 4 5 φ -0.375π, θ	                 C 11 C 12 C 13 0 C 22 C 23 0 C 33 0 C 44 0 0 0 0 C 55 0 0 0 0 0 C 66
		Rhombohedral Trigonal-trapezohedral Ditrigonal-pyramidal Ditrigonal-scalenohedral 3m 3 32 3m	S 6 D 3 C 3v D 3d	1 + Y 3 4 φ, θ + Y 6 6 φ, θ + Y 6 8 φ -0.08π, θ 1 + Y 3 4 φ -π 6 , θ + Y 3 5 φ, θ + Y 3 6 φ -π 6 , θ 1.3 + Y 3 5 φ, θ + Y 3 6 φ, θ + Y 6 7 φ, θ 1 + Y 3 4 φ, θ + Y 6 6 φ, θ	     C 11 C 12 C 13 C 14 0 C 11 C 13 -C 14 0 C 33 0 0 C 44 0 C 44	0 0 0 0 C 14 (C 11 -C 12 )/2	                	6
						Axis rotational symmetry 3
	Tetragonal	Tetragonal-pyramidal Tetragonal-disphenoidal 4 4 Tetragonal-dipyramidal 4/m Tetragonal-trapezohedral 422 Ditetragonal-pyramidal 4mm	C 4 S 4 C 4h D 4 C 4v	1 + Y 4 4 φ -0.12π, θ + Y 4 5 φ, θ + Y 4 6 φ, θ 1 + Y 2 3 φ, θ + Y 4 4 φ, θ + Y 4 6 φ -0.08π, θ 1 + Y 4 4 φ -0.12π, θ + Y 4 6 φ, θ + Y 4 8 φ, θ 1 + Y 4 5 φ -0.125π, θ + Y 4 6 φ -0.250π, θ 1 + Y 4 5 φ, θ + Y 4 6 φ, θ	                 C 11 C 12 C 13 0 C 11 C 13 0 C 33 0 C 44 0 0 0 0 C 44 0 0 0 0 0 C 66	                	6
	Hexagonal	Tetragonal-scalenohedral 42m Ditetragonal-dipyramidal 4/mmm Hexagonal-pyramidal 6 Trigonal-dipyramidal 6 Hexagonal-dipyramidal 6/m Hexagonal-trapezohedral 622 Dihexagonal-pyramidal 6mm	D 2d D 4h C 6 C 3h C 6h D 6 C 6v	1 + Y 2 3 φ, θ + Y 4 4 φ, θ 1 + Y 4 6 φ, θ + Y 4 8 φ, θ 1 + Y 6 7 φ, θ + Y 6 8 φ, θ + Y 6 9 φ -0.1π, θ 1 + Y 3 5 φ, θ + Y 3 7 φ, θ + Y 6 8 φ -0.1π, θ 1.1 + Y 6 6 φ -0.08π, θ + Y 6 8 φ, θ + Y 6 10 φ -0.15π, θ 1.2 + Y 6 6 φ, θ + Y 6 7 φ -π 12 , θ + Y 6 8 φ, θ + Y 6 11 φ -π 12 , θ 1.3 + Y 6 6 φ, θ + Y 6 7 φ, θ + Y 6 8 φ, θ + Y 6 10 φ, θ	Axis rotational symmetry 3 C 11 C 12 C 13 0                  0 0 C 11 C 13 0 0 0 C 33 0 0 0 C 44 0 0 C 44 0 (C 11 -C 12 )/2	                	5
	Cubic	Ditrigonal-dipyramidal Dihexagonal-dipyramidal 6/mmm 6m2 Tetartoidal 23 Diploidal m 3 Gyroidal 432 Hextetrahedral 43m Hexoctahedral m 3m	D 3h D 6h T T h O T d O h	1 + Y 3 5 φ, θ + Y 6 10 φ, θ 1 + Y 6 6 φ, θ + Y 6 8 φ, θ + Y 6 10 φ, θ Spherical snub tetrahedral Spherical cube Spherical snub cube Spherical tetrahedron / truncated tetrahedron	Axis rotational symmetry 3 C 11 C 12 C 12 0                  0 0 C 11 C 12 0 0 0 C 11 0 0 0 C 33 0 0 C 33 0 C 33	                	3

Appendix A.2. Discrete Green's operator iterative scheme

The method is described through displacements u ⃗ x and corresponding discrete Green's function

where

Here the reference media parameter was also computed as

The convergence criterion was set as

Finally, the scheme of the method looks as follows

The initial condition is

Appendix B. Symmetry notation and definitions

We describe the notation used for the different symmetries and the corresponding definition in Table B