Fig. 1. The cells of a quadrilateral mesh are optimized to become quasi-rhombic; then material properties (variable Young's moduli and Poisson's ratio) are assigned to the cells. The assigned material properties are used to evaluate the geometric parameters of a tileable microstructure, encoded with a smooth spline map.

New fabrication technologies have significantly decreased the cost of fabrication of shapes with highly complex geometric structure. One important application of complex fine-scale geometric structures is to create variable effective elastic material properties in shapes manufactured from a single material. Modification of material properties has a variety of uses, from aerospace applications to soft robotics and prosthetic devices. Due to its scalability and effectiveness, an increasingly common approach to creating spatially varying materials is to partition a shape into cells and use a parametric family of small-scale geometric structures with known effective properties to fill the cells.

We propose a new approach to solving this problem for extruded, planar microstructures. Differently from existing methods for two-scale optimization based on regular grids with square periodic cells, which cannot conform to an arbitrary boundary, we introduce cell decompositions consisting of

Advances in fabrication of highly complex geometry using additive manufacturing and other technologies resulted in new opportunities for shape design. In particular, small-scale, topologically and geometrically complex structures make it possible to achieve variable effective material properties using a single material for fabrication, including material properties not easily obtained by other means (such as negative Poisson's ratio) or material properties needed for precise control of deformation behavior of a shape. Realizing the potential of microstructures requires automatic generation and high-level control of the geometry, absent from commonly used geometric modeling tools.

A variety of approaches were developed for generating small-scale structures. Global methods, like topology optimization, are flexible but require very expensive computations to obtain high-quality results for fine-scale structures (cf. [Aage et al. 2017]). Other approaches partition the problem into different scales, using a variety of methods to generate a small-scale structure locally from a coarsescale assignment of spatially-varying target properties. In this paper we present a method for design of small-scale structure families supporting an approach of the second type; specifically,

• Partition an input shape into quadrilateral cells, possibly with irregular connectivity; each cell is assigned target elastic properties;

• Assign to each cell a geometric microstructure, chosen from a family of such structures directly parameterized by their effective elastic properties and cell shapes.

The main goals for the choice of the microstructure family include:

(1) cover a broad range of material properties, (2) be simple enough for fabrication at a small scale (i.e. avoiding thin features and small holes), (3) handle a range of cell shapes, (4) be easily tileable without significant modifications, (5) depend smoothly on the target effective elasticity tensor to avoid discontinuities in transitioning between varying material properties, and (6) be efficient to compute to enable tiling of large lattices.

The question of designing such microstructures is extensively studied in the literature, although advances in design of practical families covering large ranges of material properties are more recent. Existing techniques fall somewhat short of meeting the requirements above. First, these methods are exclusively based on periodic cell tilings with squares, regular triangles, or hexagons: only a single cell shape is used. This restricts the shapes that can be tiled to those constructed of this single cell type. Other shapes require either cutting cells, or deforming cells; both changes to the cell shape lead to a significant deviation from the intended deformation behavior (Section 8).

Most closely related previous work constructs large libraries of cell geometries [Panetta et al. 2015;Schumacher et al. 2015], with a separate shape or topology optimization performed for a dense set of values of material properties (Young's modulus and Poisson's ratio combinations for isotropic materials), and a specific base material. While this is generally adequate, if variable cell shapes are used, the size of the needed library will become prohibitively large (and still be limited to a discrete set of shapes). Furthermore, a different set of microstructures needs to be computed if some aspects of the properties of the base material change. These approaches do not provide the guarantee of smooth dependence of the resulting geometry on the material properties, and, as a consequence, a large number of samples are needed and interpolation between samples cannot be used.

In this paper, we address these problems for extruded, planar microstructures, aiming to meet the requirements listed above. The main aspects of our approach include the following.

• We choose to use rhombic cells, among different classes of tileable cells, since they have the flexibility to approximate arbitrary boundaries and can be connected in arbitrary orientations while being sufficiently simple to define microstructures for all members of the family in a compact way.

• We develop a parametric family of structures (Figure 2) completely described by: (1) Eight geometric parameters, defined as smooth spline functions P (𝐸, 𝜈; 𝛼) of material properties (Young's modulus 𝐸 and Poisson's ratio 𝜈 in the isotropic case) and the cell shape parameter (rhombus angle 𝛼).

(2) The domain of this parameterization in (𝐸, 𝜈) space is defined by three simple linear constraints ℓ (𝛼) dependent on the angle. This enables efficient optimization of material properties.

• Our family is tileable, i.e., the structures for adjacent cells connect with little or no modifications. Our tiles have matching topology and positions in the boundary. In addition, due to the smoothness of our mapping functions, if the transition of materials is smooth, the change in the connection region will also be smooth and almost no difference exists in the shape between neighboring cells.

• We demonstrate that the set of geometries is universal, i.e., can be used for any base material (although the ranges of realizable material properties do change).

• We present an algorithm for optimizing a quadrangulation of a planar domain to minimize the deviation of quads from the rhombic shape and demonstrate that this approach yields nearlyrhombic cells so that our microstructure family can be used.

We will make the data for computing microstructure geometric parameters, and constructing microstructure geometry from these parameters, publicly available.

We build on the foundation of previous work on construction of microstructure families and mesh quality optimization.

Periodic homogenization. Homogenization is a central tool in our construction. We use an extension of the FEM-based formulation used in [Panetta et al. 2017[Panetta et al. , 2015]], which, in turn, goes back to [Allaire 2002], and is widely used in the literature. [Schumacher et al. 2018] presents an efficient homogenization approach tailored to rod structures. In these papers, regular lattices are used. We instead

2. An array of microstructures with Young's modulus 𝐸 varying from 0.005 to 0.1 and 𝜈 varying from -0.15 to 0.5. Top: obtained by pointwise inverse homogenization of [Panetta et al. 2017]; Bottom: our family.

perform homogenization on a range of rhombic cell shapes, by transforming the problem to a domain where the cell is square, and replacing the target tensor with a transformed one. Transformations of the homogenized material properties when the base material is changed are considered in [Cherkaev et al. 1992].

Microstructure design and optimization. Many papers considered different aspects of microstructure design, see e.g., books [Allaire 2002;Cherkaev 2000;Cioranescu and Donato 1999;Milton 2002;Torquato 2002], and references in [Schumacher et al. 2015] and [Panetta et al. 2015], which represent topology and shape optimization approaches to microstructure optimization respectively. We refer to these methods as pointwise inverse homogenization, as they construct microstructures separately for each combination of effective material parameters.

Other examples of topology-optimization based work can be found in [Bendsùe 1989;Bendsùe and Sigmund 2003;Chen et al. 2018;Nakasone and Silva 2010]. Initially, mostly the problem of identifying extremal microstructures are considered (i.e., microstructures with properties at the boundary of the ranges that can be achieved). Recent works, e.g., [Ostanin et al. 2018;Panetta et al. 2015;Schumacher et al. 2015;Zhu et al. 2017] consider the question of constructing families spanning a broad range of elastic properties. In particular, [Ostanin et al. 2018] shows that near-optimal ranges of isotropic material behaviors can be achieved in 2D with hexagonal and triangular cells and a small number of microstructure parameters. We discuss the differences from this work in Section 4. [Milton et al. 2017] provides a characterization of achievable elastic tensors in terms of energies in 2D and 3D.

Most of the works designing microstructure families use an ad hoc approach to connecting structures corresponding to adjacent cells. [Garner et al. 2019] takes a more systematic approach for structures obtained using topology optimization by adding additional terms to the functional. A recent concurrent work [Martínez et al. 2019] introduces a metric that can be used to interpolate between a variety of microstructures on 2D regular grids, allowing to create smooth variation of material properties similar to our construction (but limited to regular lattices). [Konaković-Luković et al. 2018] uses a special type of 2D triangular auxetic structure to effect conformal surface deformations. This method requires domain meshing with triangles close to regular. Similarly, a recent paper [Malomo et al. 2018] uses 2D spiral microstructures for controlling deformation of sheets into a target shape. Varying geometric properties of spirals allows for a restricted control over material properties.

As an alternative to periodic microstructures, [Martínez et al. 2016;Martínez et al. 2017] construct randomized printable structures with control over Young's moduli both for isotropic and anisotropic target properties, but cannot independently control the Poisson's ratio. [Ion et al. 2016[Ion et al. , 2019] ] describe a simple set of small-scale twodimensional structures that can be assembled into mechanisms, and a computational design tool creating this type of mechanisms. While a range of properties can be obtained by changing the basic structure parameters their properties are difficult to control precisely.

A number of works apply the two-scale approach in a different way in 2D, using a simple rectangular cell structure and obtaining a directional field and scalar fields for geometric parameters defined on a regular grid. A field-aligned coarse mesh is extracted from the directional field and then filled with microstructures with parameters determined by the scalar fields [Gil-Ureta et al. 2019;Groen and Sigmund 2017;Groen et al. 2019]. Compared to these works, we do not require the mesh cell orientations to be physically meaningful: in our case the effective material distribution and meshing are independent.

Global topology optimization. In [Aage et al. 2017;Liu et al. 2018;Wu et al. 2016], topology optimization was scaled up to highresolution uniform and adaptive 3D grids. [Wu et al. 2018] performs high-resolution topology optimization with additional constraints to create an evenly distributed porous small-scale structure minimizing compliance for specific loading scenarios. The two-scale methods based on microstructures can always be combined with topology optimization to improve efficiency or resolution achievable in a given time.

Fabrication. [Hollister 2005;Kang 2010;Lin et al. 2004a,b] have demonstrated fabrication of optimized microstructures in the context of bone scaffold and fusion cage design. [Andreassen et al. 2014;Bückmann et al. 2012;Greaves et al. 2011;Schwerdtfeger et al. 2011] have shown the possibility of manufacturing auxetic materials. The idea of manufacturing objects with spatially varying properties using tileable structures also appears in [Hiller and Lipson 2009]. [Bickel et al. 2010] designs and fabricates objects satisfying an input deformation by optimizing for the best combination of stacked layers of their multi-material 3D printer's base materials. [Skouras et al. 2013] applies discrete material optimization to achieve desired deformations of complex characters with actuation, fabricating the results with multi-material printing. et al. 2015] enforces approximate printability via constraints on the skeleton defining the structure. In our planar case, there is no need to impose fabrication constraints other than requiring that the geometry should stay connected, and the minimal feature size is large enough for the chosen fabrication process.

Mesh optimization. The optimization of the element shape of discrete meshes has been studied in many disciplines. In the context of finite element simulations, the shape is optimized to reduce the distortion introduced by the geometric map [Brewer et al. 2003;Livesu et al. 2015]. Similarly, for texture mapping and quadrilateral meshing applications the distortion of a map from a surface to a plane is minimized by either evolving a Tutte's parameterization [Aigerman et al. 2014;Degener et al. 2003;Fu et al. 2015;Hormann and Greiner 2000;Kovalsky et al. 2016;Poranne and Lipman 2014;Sander et al. 2001;Schüller et al. 2013;Shtengel et al. 2017;Smith and Schaefer 2015;Sorkine et al. 2002;Zhu et al. 2017] or recovering from a possibly inverted initial guess [Aigerman and Lipman 2013;Fu and Liu 2016;Kovalsky et al. 2015;Lipman 2012]. These maps are also commonly used for mesh deformation applications [Bouaziz et al. 2012;Sorkine and Alexa 2007], and special constraints are often used in architectural geometry to generate meshes with planar faces [Bouaziz et al. 2014;Deuss et al. 2015;Jiang et al. 2015;Poranne et al. 2013;Tang et al. 2014].

We introduce an algorithm to optimize a quadrilateral mesh to have rhombic elements, a requirement that, to the best of our knowledge, has never been studied before. Our algorithm follows the paradigm introduced in ShapeUp [Bouaziz et al. 2012], alternating the projection to the space of rhombic shapes and a continuous optimization, and it relies on the solver introduced in [Rabinovich et al. 2017] to speed up convergence.

Problem. Our goal is to approximate a distribution of (potentially variable) material properties given by a spatially variable symmetric elasticity 4-tensor s 𝐶 (x) on a two-dimensional polygonal domain Ω, by partitioning it into constant-material cells, and assigning a microstructure to each cell. We assume that a single base material with known elastic properties given by the tensor 𝐶 base is used for all small-scale geometry.

The geometry of each cell 𝑄 is chosen so that its homogenized or effective elasticity tensor matches that of the assigned material s 𝐶 (𝑄).

Informally, the effective elasticity tensor can be understood as follows. If the cell 𝑄 were repeated periodically, a block of material consisting of sufficiently large number of such cells would behave as if it were made of a homogeneous material with elasticity tensor s 𝐶 (𝑄). In this work, we focus on isotropic target material distributions, for which the elasticity tensor is defined by a pair of parameters 𝑀 = (𝐸, 𝜈), Young's modulus and Poisson's ratio.

We focus on isotropic materials since the effective orthotropic elasticity tensors that rhombic cells can produce have principal axes aligned with the diagonals of the rhombi; while for some purposes this broader space can still be used, the space of possible materials becomes mesh-dependent, complicating and restricting material optimization. In contrast, considering isotropic materials makes the result much more mesh-independent.

As mentioned in Section 1, our goal is to design the geometry in individual cells to have the following properties: (1) cover a broad range of material properties, (2) be simple enough for fabrication at small scale (i.e. avoiding thin features and small holes), (3) handle a range of cell shapes, (4) be easily tilable without significant modifications, (5) depend smoothly on the target effective elasticity tensor to avoid discontinuities in transitioning between varying material properties, and (6) be efficient to compute to enable tiling of large lattices. Notation. We use rhombic cells 𝑅 and a specific parametric cell geometry structure, which respects the two reflection symmetries of the rhombus. The shape of the rhombic cell is described by a single parameter, 𝛼 ⩽ 𝜋/2, the smaller angle of the rhombus.

The geometry of each cell is defined by a vector of parameters p ∈ R 𝑛 , 𝑛 = 8, (Figure 3); the geometry is generated from parameter values using a 2D version of the method of [Panetta et al. 2017], described in Appendix A. The reasons for this choice are summarized in Section 4. We assume that for each parameter a range [𝑝 𝑚𝑖𝑛 𝑖 , 𝑝 𝑚𝑎𝑥 𝑖 ] is provided. The 𝑛-dimensional box of admissible parameter values is denoted 𝐵 𝑝 . In particular, for any choice of parameters, the resulting homogenized elasticity tensor is orthotropic, and completely determined by four components, e.g., Young's moduli 𝐸 1 , 𝐸 2 along two directions, shear modulus 𝐺, and one of two Poisson's ratios 𝜈 12 .

The function 𝐻 (p; 𝛼) : 𝐵 𝑝 → R 4 , the homogenization function, computes the effective elasticity tensor components of a cell with angle 𝛼 from its geometry parameters. We describe a way to compute this function on arbitrary parallelogram cells, including rhombic, in Section 5. We use 𝐻 𝐸 (p) and 𝐻 𝜈 (p) to denote the Young's modulus and Poisson's ratio corresponding to 𝐻 (p) when it is isotropic. Solution overview. Our approach is composed of two main components:

• We construct a parametric family of cell structures solving the inverse homogenization problem for a range 𝐷 (𝛼) ⊂ 𝐻 (𝐵 𝑝 ; 𝛼) of material properties (𝐸, 𝜈) ∈ 𝐷 (𝛼), where 𝐸 is the Young's modulus and 𝜈 is the Poisson's ratio, for a range of cell angles 𝛼, [𝛼 𝑚𝑖𝑛 , 𝜋/2]. This family is described by the material-togeometry map P (𝐸, 𝜈; 𝛼) : 𝐷 (𝛼) → 𝐵 𝑝 .

• We optimize a quadrangulation to obtain quasi-rhombic cells, so that the map P (𝐸, 𝜈, 𝛼) can be used to fill in the small-scale geometry.

The difficulty of the problem of inverse homogenization is due to several factors. The homogenization map 𝐻 (p) is straightforward, but expensive, to compute, because it requires several finite-element solves as a part of its computation. More fundamentally, for a fixed 𝛼, we need to invert the map 𝐻 : 𝐵 𝑝 ⊂ R 𝑛 → R 4 on a two-dimensional subspace Iso ⊂ R 4 of isotropic materials. This is challenging, since the inverse is far from unique. Last but not least, the map depends on the base material properties 𝐶 base , in addition to the homogenized parameters and cell geometry.

We introduce a novel, smooth, closed form material-to-geometry map P (𝐸, 𝜈; 𝛼, 𝐶 base ) covering a broad range of material properties and rhombus angles 𝛼, that can be used for an arbitrary isotropic base material with tensor 𝐶 base . This map is uniquely determined by selecting a suitable four-dimensional subspace of the geometric parameter space R 𝑛 , and can be represented in a very compact form by a set of 3D spline.

In the next sections, we discuss our choice of cell shape (Section 4), review homogenization (Section 5) and show how it can be extended to rhombic cells. Then we explain how we solve the inverse homogenization problem (Section 6), and how to compute a tiling with quasi-rhombic cells (Section 7).

The choice of cell structure (i.e., topology of small-scale geometry, and parameters defining the geometry) is not unique: many choices may have similar behavior. We outline the heuristics we have used to select our structure for quadrilateral cells and briefly compare to alternatives.

Our choice of structure follows the general procedure outlined in [Panetta et al. 2015], applied to planar square cells. On the one hand, we want to minimize the number of parameters and topological complexity of the structure, as complex cells are difficult or impossible to manufacture. On the other hand, we want to maximize the coverage that can be achieved by the structure, that is, the range of material properties (𝐸, 𝜈). To obtain a rough prediction of coverage, we run a coarse sampling of geometric parameter space, evaluating the material parameters for each sample. In these sweeps, we enforce square symmetry (Figure 5), so that the resulting elastic tensor has only three free parameters (𝐸, 𝜈, 𝐺), and we look at the coverage we get in the projection to the (𝐸, 𝜈) plane as an indicator of what one can expect for inverse homogenization.

Among simplest topology patterns, the specific one we have chosen has the largest area covered. Expanding this area requires increasingly complex topology which negatively affects manufacturability. Choice of cell shape. As we would like to partition arbitrary shapes into cells, using square cells is not possible: in general, we cannot conform to an arbitrary boundary without introducing cells of other shapes. At the same time, the cells need to be close to periodically tileable, i.e., if we use quadrilateral cells, close to parallelograms. Parallelograms have only central symmetry, which does not restrict the possible homogenized (averaged) elasticity tensors much; on the other hand, rhombic cells have two reflectional symmetries with respect to diagonals, and, as a result (Section 6), its elasticity tensor is orthotropic, as long as the small-scale geometry satisfies the same constraint. This considerably simplifies the construction of the material-to-geometry properties map, by decreasing the dimension of the space of possible materials.

As an alternative, regular triangular and hexagonal cells as shown in [Ostanin et al. 2018] can span a large range of material properties, and have an important advantage of being isotropic by construction, i.e., any geometry can be used on a cell as long as it has the symmetries of the cell (regular triangle or regular hexagon). While this significantly simplifies the problem of constructing geometries with target properties for regular cells, using distorted triangular or hexagonal cells as we need to tile an arbitrary shape, nullifies this advantage.

For triangular cells, we can enforce an additional symmetry by requiring the cells to be isosceles; in a periodic tiling, this is equivalent to using rhombic cells. At the same time, triangular cells are more restrictive in terms of choosing structure topology, if it needs to have the symmetry of the regular triangle. Finally, and perhaps most importantly, due to non-existence of regular tetrahedral tilings in 3D, methods based on quadrilateral tilings generalize better to 3D. For these, reasons, we choose rhombic cells.

To define the actual microstructure geometry, we use an implicit function parametrized by an edge skeleton vertices and radii, using the method of [Panetta et al. 2017], This method, on the one hand, allows for explicit control of topology through the skeleton connectivity, on the other hand, allows for merges of parts of the structure when these run into each other due to parameter choices, which is critical for robustness of the optimization process.

We briefly review the standard homogenized tensor computation on square cells and show how it extends to rhombic cells. We also discuss how the inverse homogenization problem can be solved using optimization.

Our formulation follows [Panetta et al. 2015]. Suppose that a planar domain Ω is tiled by identical square cells 𝑄, each with the same microstructure 𝜔. At a large scale, we can consider averaged deformations: in the limit of infinitesimal cells, these deformations can be viewed as linear on each cell, and the actual deformation as a sum of the averaged one and a local fluctuation on the cell. The averaged deformation is the elastic deformation when the solid Ω is viewed as a completely filled volume with variable material properties defined by cell small-scale geometry, in the limit of zero cell size. The averaged deformation s u satisfies the macroscopic elasticity equation

where the elasticity tensor s 𝐶 is the effective elasticity 4-tensor with entries s 𝐶 𝑖 𝑗𝑘𝑙 and the 𝜀 (u) ≔ On a single cell, in the limit of zero cell size in a periodic tiling, we can view elastic deformations as having a constant averaged strain 𝜀 (s u) and a periodic łmicroscopic fluctuationž component w, having zero average strain by periodicity. We consider fluctuations w for three constant average strains forming a basis for all possible constant strains 𝑒 11 = e 1 ⊗e 1 , 𝑒 22 = e 2 ⊗e 2 , and 𝑒 12 = e 1 ⊗e 2 +e 2 ⊗e 1 .

Then the expressions for the components of s 𝐶, are given by [Panetta et al. 2017]:

To obtain the microscopic fluctuation term, the force balance equation in the cell 𝑄 is solved for each of the three basis strains:

Inverse homogenization and shape derivatives. The conceptually straightforward, but technically difficult, way to solve the inverse homogenization problem for a given target tensor s 𝐶 * is by shape optimization, i.e. by varying the geometry 𝜔 using shape parameters p as variables to minimize a functional. In the case of inverse homogenization, the functional penalizes the difference of the homogenized elasticity tensor and the target s 𝐶 * :

) with the norm taken to be e.g., the Frobenius norm. A combination of regularizing terms can be used as 𝐸 𝑟𝑒𝑔 , the most common being the volume of the structure defined by p and the proximity to the initial solution [Panetta et al. 2017[Panetta et al. , 2015]].

Unfortunately, due to non-uniqueness of the solution, resulting geometry parameters typically do not depend smoothly on the input material properties which has a number of negative consequences. We compare direct pointwise optimization to our approach in Section 8. To minimize the functional (4) efficiently, we need derivatives of 𝐽 (p) with respect to the geometry parameters p. Note that this requires differentiating s 𝐶 (p); the computation of s 𝐶 requires solving elasticity PDEs on a domain depending on p. This type of derivatives are called shape derivatives; these are computed by solving an adjoint PDE, which in our case is identical to elasticity, but with redefined volume forces. The details of the computation can be found in [Panetta et al. 2017].

Discretization. To solve the PDEs (3), as well as similar equations needed for shape derivative computation, we use a standard FEM discretization. The domain is remeshed at each iteration using marching squares to extract the boundary and Triangle [Shewchuk 1996] to mesh the domain.

We discretize the cell problems with quadratic triangle elements, which we found essential for accurate stress and homogenized tensor evaluation. We use straight-edged elements (subparametric FEM) for representing the geometry to simplify meshing and the shape derivative formulas (edge nodes are placed at the edge midpoints).

In our experiments, running 2D homogenization usually took under 1s for meshes with an order of 4000 vertices. The mesh resolution (and hence the number of vertices) is determined by the maximal triangle area constraint (flag -a in Triangle) and the marching squares grid size (used to compute the boundary polygon), which we chose to be 5 × 10 -4 and 256 respectively for most of our experiments. Figure 7 shows an example of the resolution we used.

Computing the homogenization map 𝐻 (p) is the essential step in our process for obtaining a parametric family of structures. We use a modification of the standard homogenization method on square cells to handle arbitrary rhombi. Homogenization for an arbitrary rhombic cell (more generally parallelogram) can be transformed to homogenization on a square cell by a change of basis for elasticity tensors 𝐶 base and s 𝐶.

Suppose we have a cell 𝑅 (rhombus with unit sides) and we want to write all our equations on a domain 𝑄 (unit square). Let 𝐹 be the affine map 𝑄 → 𝑅 and let its inverse be 𝐺.

Proposition 1. Let (e 1 , e 2 ) be a nonorthogonal basis aligned with the sides of 𝑅; in this coordinate system, the domain is a unit square, i.e., the map 𝐹 maps (e 1 , e 2 ) to the unit coordinate vectors. If displacements u satisfy the elasticity equation on domain 𝑅, -∇ • [𝐶 : 𝜀 (u)] = f, then if ũ and f are the displacements and forces expressed in the coordinate system (e 1 , e 2 ), they satisfy the elasticity equation

where the transformed tensor 𝐶 ′ components are given by

The proposition is verified directly by the change of variables in ( 1) and (3).

Note that the components of 𝐶 can be expressed in terms of the components of 𝐶 ′ using the inverse relation:

This observation applies both to the microscopic equation on the domain in 𝑅 filled with the base material, as well as to the homogenized equation. This leads to the following 3-step procedure for determining the homogenized properties of a rhombic cell:

• Transform the base material tensor 𝐶 base to the coordinate system (e 1 , e 2 ) using ( 6).

• Compute the homogenized tensor s 𝐶 ′ on 𝑄, as described in Section 5.1.

𝐶 back to the original coordinate system, using (7).

This allows us to solve the forward problem: given the geometry in a rhombic cell, compute the homogenized elasticity tensor. As the problem is reduced to square-cell problems, shape derivatives are computed exactly in the same way, and the problem is discretized using the same approach.

We describe our solution to the inverse problem in the next section.

In this section, we explain the steps for constructing a map P from isotropic material parameters to geometry parameters. We assume in most of the exposition that 𝐶 base is fixed, and show how this assumption can be removed in Section 6.3. We use 𝐸 = 1 and 𝜈 = 0 for the base material, but the specific choice does not matter.

The fact that the rhombic cells always have orthotropic elasticity tensor is of critical importance in our construction, since it determines the dimension of the image of 𝐻 (p, 𝛼). We consider this fact in more detail and define the measures of deviation from isotropy tailored for orthotropic materials.

The rhombus has two reflection symmetries with respect to its diagonals, and we consider geometries that by construction have the same symmetries. This means that the periodic structure obtained by repeating the elementary cell 𝑅 is invariant with respect to these two transformations. As a consequence, the elasticity tensor s 𝐶 is also invariant with respect to these transformations. This implies orthotropy:

Proposition 2. [Love 1944] If an elasticity tensor is invariant with respect to reflections about two orthogonal axes, then it is orthotropic with respect to these axes and in the coordinate system aligned with these axes has the form (in Voigt notation)

This means, in particular, that only two additional constraints need to be satisfied to obtain an isotropic elastic tensor:

We use two anisotropy measures, corresponding to these two constraints:

We note that principal directions of the tensor are exactly the diagonals of the rhombus, so cannot be set independently of the mesh directions.

6.2 Map construction for 𝛼 = 𝜋/2 (squares)

Recall that the material-to-geometry map P (𝐸, 𝜈; 𝛼) for a given 𝛼, is not uniquely defined; we describe a simple way to obtain a unique initial map on a smaller initial domain 𝐷 (𝛼), for 𝛼 = 𝜋/2, which we then expand. For clarity, we drop the dependence on 𝛼 in this section.

We explain the general form of construction, to clarify how it can be used for other types of structures (e.g., 3D), and explain how it applies in the specific case of planar rhombic structures. Recall that 𝑝 = 8 is the number of geometric parameters and we denote 𝑚 = 4 the number of material parameters for 2D orthotropic materials.

We start with the homogenization map 𝐻 : 𝐵 𝑝 → R 𝑚 , 𝑝 ⩾ 𝑚. In general, a 𝑝 -𝑚-dimensional submanifold of 𝐵 𝑝 is mapped to each orthotropic elasticity tensor in 𝐻 (𝐵 𝑝 ). 𝐼𝑠𝑜 corresponds to the set of possible isotropic materials in the space of material parameters R 𝑚 , e.g., it is defined by equations 𝑎 1 = 0 and 𝑎 2 = 0 for orthotropic materials.

Our key idea is to restrict 𝐻 to a carefully chosen 𝑚-dimensional affine subspace 𝑉 ⊂ R 𝑝 , for which 𝐻 has large coverage, i.e. 𝐻 (𝑉 ∩ 𝐵 𝑝 ) ∩ Iso, the set of isotropic material properties covered by geometries in 𝑉 , is not too far from 𝐻 (𝐵 𝑝 ) ∩ Iso. We refer to 𝑉 as the transversal subspace.

𝐻 ′ : 𝑉 → R 𝑚 , the restriction of 𝐻 to 𝑉 , is locally injective near all points where it is non-degenerate.

The map p ortho : 𝐷 = 𝐻 (𝐵 𝑝 ) → R 𝑝 , can be constructed by inverting 𝐻 ′ as described below, and its restriction to isotropic materials Iso ∩ 𝐷 yields the desired material-to-geometry map p. Note that if

Fig. 8. Transversal subspace 𝑉 illustrated for the simplified case of a 3parametric structure and 2-dimensional property space (𝐸, 𝜈) (In our construction, we consider spaces that are 8 and 4-dimensional respectively).

The curved surfaces are isosurfaces of 𝐸 (p). The map 𝐻 , restricted to 𝑉 , is bijective.

non-isotropic orthotropic materials are desired, this map also can be used to compute these.

As the 𝐻 (p) depends on p smoothly, almost everywhere we expect P (𝐸, 𝜈) to be smooth as well: the only places where a jump is possible is in the areas where the the differential of 𝐻 (p) is degenerate, or at the boundary of 𝐻 (𝐵 𝑝 ), which naturally defines the domain 𝐷 for P (𝐸, 𝜈). In practice, we have not observed the degeneracies of P (𝐸, 𝜈) for the shape and material parameter ranges that we have considered, and the smooth approximation we obtain matches the samples closely.

Defining the transversal subspace 𝑉 . We start with identifying a 𝑚-dimensional affine subspace containing a large range of geometric configurations with isotropic elasticity tensor, i.e., has large isotropic coverage. A direct approach to this would be to sample 𝐻 sufficiently densely on 𝐵 𝑝 , select points for which the anisotropy measures 𝑎 1 and 𝑎 2 are small, and fit an 𝑚-dimensional hyperplane to this set, e.g., via PCA. Unfortunately, densely sampling the 𝑝-dimensional set (even for 𝑝 = 8) is prohibitively expensive, as each evaluation of 𝐻 requires meshing and solving three elasticity problems at sufficiently high resolution.

By considering a coarse sampling of 𝐵 𝑝 , we experimentally observed that, for our setting, 𝐵 𝑝 contains an axis-aligned hyperplane subset of dimension 𝑝 ′ , 𝑝 ⩾ 𝑝 ′ ⩾ 𝑚, which has good coverage. Intuitively, this corresponds to fixing the geometry parameters that can be fixed without decreasing the range of material properties; we refer to these as redundant. This reduction leads to a practical two-stage algorithm: first, remove the redundant parameters that have the least effect on coverage. Once no dimensions can be dropped, use denser sampling and PCA for the final reduction to 𝑚 dimensions. In our case, 𝑝 ′ = 5, so only one extra dimension needs to be eliminated.

Next, we discuss these two steps in greater detail:

(1) Finding redundant parameters. We compute a coarse regular grid sampling 𝑞 𝑗 = 𝐻 (p 𝑗 ) of the full homogenization map 𝐻 : 𝐵 𝑝 → R 𝑚 . Each geometric parameter 𝑝 𝑖 is uniformly sampled at values 𝑐 𝑖 𝑘 , 𝑘 = 1 . . . 𝑛 𝑖 . In our case, a total of 540𝑘 data points were collected (see more information in Appendix B). Then we extract a subset of p 𝑗 with the deviation from anisotropy sufficiently small (in our case 0.005 for both 𝑎 1 and 𝑎 2 ). Consider slices 𝜋 𝑖𝑘 consisting of all samples p 𝑗 with a fixed value of 𝑝 𝑗 𝑖 = 𝑐 𝑖 𝑘 : by direct check of the coverage area of each slice, we identify the geometric parameters 𝑖 for which there is a 𝑘 such that coverage of 𝜋 𝑖𝑘 is close to complete coverage, i.e., the coverage area obtained with the initial sweep, where we use all parameters as variables. Figure 9 shows the observed coverage for a sequence of expanding sets of fixed parameters.

For the rhombic structures, we observe that the offset parameters 𝑝 1 , 𝑝 2 , and the thickness parameter 𝑝 5 all have such values 𝑐 𝑖 𝑘 (0.5, 0.25, and 0.3 respectively). This is a particularly convenient set of parameters to fix, since this allows us to merge cells perfectly as the radius on the boundary is fixed and thus the same for all cells. Fixing 𝑝 1 , 𝑝 2 , and 𝑝 5 reduces the dimension of the geometry space to just 5, leading to a 5D affine subspace 𝑉 5𝐷 . It remains to reduce the dimension by 1.

(2) Reduce to 𝑚 dimensions by PCA. The restriction of the map 𝐻 to 𝑉 ′ ∩ 𝐵 𝑝 , where 𝑉 ′ is the 𝑝 ′ -dimensional subspace obtained by dropping dimensions, can be sampled much more densely (we use 15 samples per dimension for 5-dimensional 𝑉 ′ ). To obtain the 𝑚-dimensional transversal space 𝑉 , we approximate 𝐻 (p) with a linear map on 𝐵 𝑝 : 𝐻 (p) ≈ 𝑄 p + 𝐻 0 , where 𝑄 is a 𝑚 × 𝑝 ′ matrix, and p is obtained from p by discarding the fixed coordinates. We compute 𝑄 and 𝐻 0 using PCA. Let 𝑊 be the matrix spanning the nullspace of 𝑄, i.e., 𝑄𝑊 = 0, and 𝑊 has dimension 𝑝 ′ × (𝑝 ′ -𝑚) and maximal rank. For all points p + 𝑊 𝑢, where 𝑢 is (𝑝 ′ -𝑚)-dimensional, the value of the map is constant. Then the condition 𝑊 𝑇 ( p -p0 ) = 𝑊 𝑇 p + 𝑐 0 = 0, where p0 is a point in 𝑉 ′ , defines an 𝑚-dimensional subspace 𝑉 ⊂ 𝑉 ′ , perpendicular to the constant-value affine spaces, on which 𝑄 p+𝐻 0 is one-to-one. For rhombic structures, we obtain 𝑝 ′ -𝑚 = 1, and 𝑊 has dimensions 𝑝 ′ × 1, i.e., it is a vector in this case. 𝑊 is approximated well by a vector with nonzero components 𝑣 3 = 1 and 𝑣 7 = 1, and with 𝑐 0 = -0.82.

This gives us a complete description of 𝑉 = {p | 𝑝 1 = 0.5, 𝑝 2 = 0.25, 𝑝 5 = 0.3, 𝑝 8 = 0.82 -𝑝 4 }. Geometrically, 𝑉 is described as the space of structures where node 2 is fixed at the default position at barycentric coordinates (0.5, 0.5, 0), the radius of node 1 connecting the cell to other cells is fixed at 0.3, and the radius and displacement of node 4 are related linearly.

Inverse of 𝐻 ′ on the isotropic subspace. First, we construct a piecewise-linear interpolant of the samples of 𝐻 ′ . The map 𝐻 ′ is defined as the restriction of 𝐻 to 𝑉 ∩ 𝐵 𝑝 . We sample 𝐻 ′ on a regular grid in 𝑉 . For the rhombic structure, we use 𝑝 3 , 𝑝 4 , 𝑝 6 , 𝑝 8 as coordinates for sampling. We sample at the points of the grid 𝑞 𝑗 = (𝐸 𝑗 , 𝜈 𝑗 ) contained inside 𝐵 𝑝 .

As samples 𝑞 𝑗 are on a regular 𝑚-dimensional grid, they form mD cubes, which we then subdivide into 𝑚-simplices. We discard all simplices with one of the two possible orientations (as 𝐻 ′ is not guaranteed to be injective everywhere, simplices may have both orientations; we choose the predominant orientation to be the łcorrectž one). The map is defined by linear interpolation on each simplex.

For each sample 𝑞 on a regular 2D grid of samples in 𝐼𝑠𝑜, we search for m-simplex 𝑆 ℓ in 𝑉 , such that 𝐻 (𝑆 ℓ ) contains 𝑞. While the simplices found in this way may be not unique, in practice we do not observe this problem. The value of P (𝑞) is obtained by interpolating the values of p at the corners of each simplex 𝑆 ℓ . Finally, we fit a B-spline to the sampled values of 𝑞 with Laplacian regularization. We use 20 × 20 control points for rhombic structures.

Increasing coverage of the material-to-geometry map. Constraining the map 𝐻 to 𝑉 makes the map invertible, and we observe the inverse to be smooth. However, this also restricts the coverage, as the inverse image 𝐻 -1 (𝐼𝑠𝑜 ∩ 𝐻 (𝐵 𝑝 )) in 𝐵 𝑝 is not necessarily contained in 𝑉 . We can expand the range by the following procedure using optimization-based inverse homogenization, similar to [Zhu et al. 2017].

For each regular-grid sample point 𝑞 𝑗 in (𝐼𝑠𝑜), for which P (𝑞 𝑗 ) is not defined, but which is adjacent to a point 𝑞 ℓ where P is defined, we initialize the shape optimization for the functional (4) using the spline approximation to P obtained above evaluated at 𝑞 𝑗 . Then we optimize using all 8 geometry parameters as variables (and not only the 4 independent variables used in the sampling procedure), thus taking the value of P (𝑞 𝑗 ) out of 𝑉 . If the optimization converges to a value sufficiently close, using a tolerance of 0.005 for Poisson's ratio and 0.001 for Young's modulus, we include the additional point in the set, and refit the splines to the new set of points.

We note that, in principle, this process may suffer from the same flaws as the original inverse homogenization of [Panetta et al. 2015] (cf. Figure 10) as the map P (𝐸, 𝜈; 𝛼) may become non-smooth for lower values of 𝛼. However, the procedure defined above prevents this from happening, and the map P (𝐸, 𝜈; 𝛼 𝑚𝑖𝑛 ) remains smooth (Figure 18). Bounds for the domains 𝐷 (𝛼). For applications using optimization of material properties on a partitioned domain (Section 8), it is useful to ensure that the values used in the optimization stay within the coverage zone of the material-to-geometry function p.

To keep these constraints efficient, we approximate the coverage area by a convex polygon bounded by a set of a maximum of 6 lines: 2 horizontal lines defining minimum and maximum Young's modulus 𝐸 𝑚𝑖𝑛 and 𝐸 𝑚𝑎𝑥 , 2 vertical lines for minimum and maximum Poisson's ratio 𝜈 𝑚𝑖𝑛 and 𝜈 𝑚𝑎𝑥 , and 2 additional slanted lines, 𝑠 𝑙 and 𝑠 𝑟 , that allow for larger non-rectangular domains (see Figure 13).

Fig. 13. Lines defining the shape of our convex domain.

The lines are defined to maximize a weighted sum of the area and the height (maximum Young's modulus) of the convex polygon while keeping it completely inscribed in the domain coverage region. Below, we explain this process in more detail.

In our coverage data, we first identify two sets of points on the boundary of our covered area. For each value of Young's modulus, the minimum achieved Poisson's ratio is added to 𝐵 𝑙 , while the maximum should be in 𝐵 𝑟 . Then, we iteratively run the following two steps, trying different values for lines 𝜈 𝑚𝑖𝑛 , 𝜈 𝑚𝑎𝑥 , 𝐸 𝑚𝑖𝑛 and 𝐸 𝑚𝑎𝑥 at each time:

(1) Filtering. We filter out boundary points in 𝐵 𝑙 and 𝐵 𝑟 that are outside the square defined by 𝜈 𝑚𝑖𝑛 , 𝜈 𝑚𝑎𝑥 , 𝐸 𝑚𝑖𝑛 and 𝐸 𝑚𝑎𝑥 , creating new boundary sets 𝐵 ′ 𝑙 and 𝐵 ′ 𝑟 .

(2) Optimization. We optimize the slanted lines as follows. We create one linear constraint for each boundary point 𝑏 in 𝐵 ′ 𝑙 (resp. 𝐵 ′ 𝑟 ), making sure that slanted line 𝑠 𝑙 (resp. 𝑠 𝑟 ) is to the right (resp. left) of 𝑏. After building our constraint matrix, we run the fmincon() solver in Matlab to optimize our objective function (composed of area and height of the convex domain) given our new linear constraints.

At the end, the objective function values obtained for all iterations are compared and we select the 6 lines corresponding to the best value achieved. In our own case, we decided to fix minimum and maximum Young's modulus (𝐸 𝑚𝑖𝑛 and 𝐸 𝑚𝑎𝑥 ) to 0.005 and 0.32, due to our observed coverage. However, we tried multiple values for maximum and minimum Poisson's ratio (𝜈 𝑚𝑖𝑛 , 𝜈 𝑚𝑎𝑥 ). Notice that these values can also be chosen according to the application.

We extend the map to all angles in the range [𝛼 𝑚𝑖𝑛 , 𝜋/2]. We use 𝛼 𝑚𝑖𝑛 = 𝜋/4 since the coverage becomes very small below this value. In our experiments, the range [𝜋/4, 𝜋/2] proved to be sufficient, as, after quad optimization, the minimum angle for quads was larger than 𝜋/4. It is however possible that for complex geometries smaller angles are needed. We perform the extension incrementally for a sequence of angles 𝛼 𝑖 = 𝜋/2 -𝑖 𝑑, where 𝑑 = 1.25 • . Using the geometric parameters for 𝛼 𝑖-1 as a starting point, we use shape optimization to obtain a point for the same target value (𝐸, 𝜈) but for 𝛼 𝑖 , with a regularization term ∥𝒑 -𝒑 0 ∥ 2 penalizing deviation from the initial value of each parameter, where 𝒑 and 𝒑 0 represent the current and the initial parameter values during the optimization, respectively.

Once the whole range of angles is covered, we fit a 3D spline in variables (𝐸, 𝜈, 𝛼) to the whole set of points. The resulting change of geometry parameters as a function of 𝛼 is smooth (Figure 15). For spline fitting, we use least-squares on the whole cube. An extra regularization term, the Laplacian of the function at grid points, is added in the cost, which promotes smoothness and eliminates the null space corresponding to the regions with a few or no data points.

Complete map description. The final map P is defined as:

• a 3D uniform spline approximation on a 12 × 12 × 12 grid of 8 geometric parameters, for a total of 8 × 12 3 coefficients;

• a 1D splines in 𝛼 for each of the coefficients defining the lines of the convex polygon approximating the domain coverage region 𝐷 (𝛼).

The complete description of the map is relatively compact, requiring less than 14k coefficients. The plots for the eight components are shown in Figure 17 for 𝛼 = 𝜋/2, while Figure 18 shows plots for two of the parameter when 𝛼 = 𝜋/4. In addition, Figure 11 illustrates (non-)smoothness of the map obtained by pointwise shape optimization, compared to our map, for small changes of 𝐸 and 𝜈 for three choices of (𝐸, 𝜈).

Change of the base material. In principle, we would need to compute a new map P each time the base material changes, as the homogenization map 𝐻 depends on the choice of material. However, two observations show that s 𝐶 for any base material (𝐸 ′ , 𝜈 ′ ) can be obtained from a base material (𝐸, 𝜈) by a simple transformation. The first simple observation is that the expressions for the components of s 𝐶 and the solution for (3) depend linearly on the base elastic tensor 𝐶 base , i.e., if the base Young's modulus 𝐸 is replaced with 𝑠𝐸, the homogenized one 𝐸 𝐻 becomes 𝑠𝐸 𝐻 .

The less trivial part are the changes due to changes in the Poisson's ratio. Fortunately, if 𝜈 is replaced with 𝜈 ′ , the new homogenized elasticity tensor s 𝐶 ′ is characterized by the CLM theorem, a corollary of which for isotropic materials can be formulated as follows.

Theorem 6.1. [Cherkaev et al. 1992] Suppose a periodic 2D medium with holes is filled with material with Young's modulus 𝐸 and Poisson's ratio 𝜈, and the effective (homogenized) properties of this medium are 𝐸 𝐻 and 𝜈 𝐻 . Suppose 𝜈 is replaced with 𝜈 ′ . Then the effective Young's modulus is preserved, and the Poisson's ratio changes as

This immediately leads to the change of variables in the function 𝐻 (𝐸, 𝜈) that corresponds to the change of base material properties (𝐸, 𝜈) → (𝐸 ′ , 𝜈 ′ ): Discussion of alternative approaches. The closest works on 2D microstructures are [Ostanin et al. 2018] and [Martínez et al. 2019], which introduce low-parametric families of structures. However, both approaches are limited to regular grids in a profound way: they rely on symmetries of triangular/hexagonal grids to produce isotropic effective properties, a requirement that is lost on deformed lattices with fewer symmetries.

Regular quad/hex grids are used by topology optimization-based methods [Schumacher et al. 2015;Zhu et al. 2017], and in the family of structured proposed in [Panetta et al. 2017[Panetta et al. , 2015]]. While these methods could be adapted to rhombic tiles by applying a bilinear map to their structures, this construction introduces large errors (see Section 8 for more details).

Our construction is based on the optimization techniques introduced in [Panetta et al. 2017[Panetta et al. , 2015]], by addressing two of their limitations:

(1) we lift the regular lattice cell assumption (for the 2D planar case), directly targeting rhombic tiles and (2) we address the lack of smoothness of P (cf. Figures 17 and11).

We tessellate the polygonal domain Ω into (nearly-)rhombic elements by starting with a quadrangulation of Ω and optimizing the shape of its cells.

Quad mesh generation. Our algorithm can be applied to any planar quadrilateral mesh: for the examples in the paper, we used a combination of meshes automatically generated using the Integer Grid Maps algorithm [Bommes et al. 2013] (Figure 34, disk), Pattern-Based Quadrangulation [Takayama et al. 2014] (Figure 34, bar), and manually designed using Blender (Figure 34, pliers, ghost). Generally, denser meshes are easier to optimize for our algorithm since they have more degrees of freedom.

Quad to rhombi projection. Our algorithm uses an alternating optimization inspired by ShapeUp [Bouaziz et al. 2012]: we compute a target rhombic shape for each quad (see inset), and then optimize the vertex positions to match the target shape of each element in least-squares sense. To obtain the ideal rhombic shape 𝑅 for a quad 𝐹 , we fit a linear transformation 𝑇 of a canonical unit square 𝑄 (a parallelogram) to 𝐹 , after translating its barycenter to the origin:

where 𝑓 𝑖 and 𝑞 𝑖 are the vertices of 𝐹 and 𝑄 respectively. Letting 𝜃 𝑇 be the smaller angle of the parallelogram, we define the unit target rhombic shape 𝑅 𝑢 as the rhombus with unit length edges and minimal angle 𝑚𝑎𝑥 (𝜃 𝑇 , 𝜃 𝑚 ), where 𝜃 𝑚 is the smallest angle by our parameter sweep (45 degrees). Capping the minimal angle pushes the optimization toward rhombic shapes within the reachable range of our parametric family. The scaled target rhombic shape 𝑅 is computed by translating (by a vector 𝑡) and scaling (by a diagonal matrix 𝐴) 𝑅 𝑢 to best fit the quad 𝐹 : min

where 𝑟 𝑖 and 𝑓 𝑖 are the vectors of vertices of 𝑅 𝑢 and 𝐹 respectively.

Exponential rhombic distortion energy. Equipped with a target reference for each quad, we define a distortion energy by penalizing the deviation of each element from its ideal rhombic shape using the symmetric Dirichlet energy [Smith and Schaefer 2015]:

where 𝐽 𝐹 is the Jacobian of the bilinear map mapping 𝐹 to its target rhombic shape 𝑅. Ideally, we would like to minimize the 𝐿 ∞ norm of the distortion integrated over the entire domain Ω, which is however a challenging problem to solve numerically. We opt for an approximation, using a continuation approach [Rabinovich et al. 2017] minimizing the following exponential energy for increasing values of 𝑠:

where 𝐴 𝑄 is the area of the quad 𝑄.

Boundary projection. Without additional constraints, minimizing (14) would change the shape of the boundary of Ω (Figure 19). We add a penalty term to avoid changes to the boundary, while allowing nodes to slide on it:

where 𝑥 𝑖 are the coordinates of a boundary vertex of a point on the boundary of the quad mesh, and 𝑦 𝑖 = Φ(𝑥 𝑖 ) is the closest point on the boundary of Ω.

Optimization. We create a rhombic mesh by minimizing the combined energy:

with 𝜆 = 10000 a weight balancing the two terms. We use the solver proposed in [Rabinovich et al. 2017] to find vertex positions minimizing eq. ( 16). This solver uses a local/global optimization strategy in order to minimize eq. ( 16). In the local step, the target rhombi are rotated to minimize 𝐸 D (𝑥) (eq. ( 4) in [Rabinovich et al. 2017]), and we compute 𝑦 𝑘 𝑖 = Φ(𝑥 𝑘 𝑖 ), the projection on the boundary 𝜕Ω for border vertices after iteration 𝑘. In the global step, the Jacobian 𝐽 𝑄 in 𝐸 D is replaced by 𝐽 𝑘 𝑄 = 𝐽 𝑄 (𝑥 𝑘 ), 𝐸 𝜕 (𝑥, Φ(𝑥)) is replaced by 𝐸 𝜕 (𝑥, 𝑦 𝑘 ), and 𝐸 is minimized wrt 𝑥.

In this section, we explore the advantages of our approach, compared to "naive" application of existing techniques, and its limitations. Specifically, we explore the following features of our method:

• the use of boundary-aligned, irregular quad meshes instead of regular grids with cut cells;

• the optimization of the cell geometry structure taking cell shape into account, compared to remapping geometry optimized for square cells;

• continuous dependence of geometric parameters on material parameters, compared to pointwise inverse homogenization used in previous work;

• mesh independence and the effects of deviation from rhombicity on the accuracy of the method.

Finally, we show how our method handles a number of test problems. Many of the comparisons are performed on models with variable material properties generated using material optimization, as in [Panetta et al. 2015]. Specifically, the deformation on one part of the boundary is given, and the material properties for interior cells are optimized to obtain a user-defined target deformation on an other part of the boundary. For example, we use a plier shape for several examples; the prescribed deformation moves the handles together, and the target deformation moves the jaws together. In the optimization, the cells of the mesh are used as finite elements with constant material properties.

Comparison with regular grids with cut cells. The simplest approach that allows one to use existing techniques with microstructures optimized for squares is to overlay a regular grid on the geometry of the shape and remove parts of the cells outside the input shape boundary. We used a regular grid with square side length equal to the average edge length of our rhombic geometry to obtain approximately the same number of cells. The material optimization was performed with the same boundary conditions and target deformations but for each mesh separately (cut regular grid vs rhombic cells), rather than, e.g., interpolating target material properties from one mesh to the other. For the purposes of FEM simulation used in the material optimization, cut cells are triangulated. In both cases we assign the geometry for each cell based on the computed material. As an additional step for the regular cell grid, we crop the final fine-scale geometry using the original object's boundary.

We run the deformation simulation on full fine-scale geometry for comparison. Figure 20 shows an example with significant deformations compared to the ground truth. We note that one feature of our method, the continuous dependence of geometric parameters on material parameters remains useful even for objects that can be meshed with regular grids, and improves accuracy. To experiment with higher resolution quad meshes, we applied one level refinement on each of the two configurations shown in Figure 20, keeping the same material distribution. The result is presented in Figure 21. For both methods, the error is reduced, although it is still considerably higher when using regular grids with cut cells. Notice, however, that (although being a valid experiment) refinement is not a very practical solution in real world scenarios since most fabrication technologies have resolution limitations. Comparison with remapping square cell structures. In principle, the structures obtained for square cells can be directly applied to rhombic cells, by simply applying a bilinear transformation to the geometry. Figure 22 shows the effect of not using structures specialized to set with specific angles: the increase in the error is substantial.

Comparison with pointwise inverse homogenization. In this experiment, the same mesh and the same variable material properties, resulting from optimization, are used with two different ways of assigning a structure to each cell (Figure 23). The first approach performs the nearest-neighbor lookup in the database of structures obtained by pointwise shape optimization [Panetta et al. 2017], the other using our smooth parametric family P (𝐸, 𝜈; 𝛼). In both cases, the structures are mapped to the arbitrary quads of the mesh by a bilinear map. We observe that the result is closer to the reference simulation using a variable-material solid. Figure 24 shows that the 𝐸 -𝜈 space coverage of two methods is close.

Mesh independence. If all cells are assigned the same effective material properties, the resulting deformation should not change significantly as we change the domain discretization. Fig. 23. Left: deformation of an object constructed using a table of structures obtained using pointwise inverse homogenization as in [Panetta et al. 2017]. Right: deformation of the object constructed using our family, with the same boundary conditions. The gray shape in the background shows the ground truth obtained by simulation on a solid quad mesh with materials assigned per quad. . Top: deformation of a shape with constant assigned materials and regular connectivity. Bottom: deformation of the same shape, with irregular mesh connectivity. The deformed pattern (simulated using FEM on a fine mesh) is compared to the simulated homogeneous mesh with equivalent material properties (𝐸 = 0.05, 𝜈 = 0.4). In both cases, a uniform vertical load is applied on the top, while we keep a 0 Dirichlet condition at the bottom.

We performed two experiments to evaluate the effect of changing the mesh. The first experiment (Figure 25) confirms that the deformations obtained for an irregular mesh on a bar, which requires different rhombic structures for different cells, are nearly the same as for a regular mesh with identical structure, optimized for the same square cell everywhere.

In the second experiment, we study how deviation from rhombic shape affects the result, when the cell shape is intentionally distorted to be far from rhombic. We create a twisted structure ś a compound cell subdivided into 4 non-square subcells, and compute the properties of the material composed of such compound cells, using homogenization on the whole cell. Specifically, starting from a 4 × 4 square quad mesh on the cell, with each quad assigned 𝐸 = 0.05, 𝜈 = 0.0, we gradually rotate the cross formed by the edges incident at the central vertex of the mesh (5 degrees at a time), and optimize the shape of the elements of the twisted mesh (keeping the cross fixed, Figure 26). Figure 27 shows how the Frobenius norm of the distance between the target elasticity tensor used to obtain the structure in each subcell and the homogenized elasticity tensor of the compound cell changes with rotation angle. In the plot, we also show the error we observe if we use the structure optimized for 𝛼 = 𝜋/2 everywhere.

Notice how the error stays stable when we increase the rotation angle, which is not the case when using only the structures for 𝛼 = 𝜋/2 in all quads. We emphasize that this experiment was intentionally conducted with a very coarse grid of subcells. The error of our method further decreases for a finer grid: if we subdivide each quad and placing the same initial structure in all sub-quads (see Figure 28). For example, for a rotation angle of 10 degrees, we had an initial error of 0.0123 with our method, while, with one level of refinement, the error is almost halved, decreasing to 0.0069, and two levels decrease the error to 0.0040.

Non-rhombicity error. In our approach, all cells are approximated by rhombi for the purposes of computing the cell structure. To study the effect of deviation of cell shapes from rhombic, we experimented with a simple example shown in Figure 29. Using a 2 × 2 quad grid in a compound cell, and starting from 4 perfect squares, we move the center vertex in the direction of the top-right corner by an increasing distance 𝑑 in each axis, which increases the nonrhombicity of every cell. For each 𝑑, the homogenized properties of the material consisting of compound cells is computed.

We compute then, for each value of 𝑑, the rhombic errors of each of the 4 quads as the sum of the distances from the quad vertices to the vertices of the closest rhombus (which can be obtained as explained in Section 7). We also compute the target error as the Frobenius norm of the difference between the homogenized elasticity tensor for the compound cell and the target one. Results are shown in Figure 30. The distance to target material changes smoothly with the total rhombic error. Moreover, to analyse quantitatively the quality of Scalability. Using the same example shown in Figure 29 with 𝑑 fixed at 0.25, we also experimented with the scalability of our solution by running the entire pipeline (quad optimization, material optimization, material to geometry mapping and final mesh construction) on different levels of refinement on the initial quad mesh.

The result is shown in Figure 31. Notice that our running time scales linearly with the number of cells of the quad mesh.

Merge error. Once the geometry is defined per cell, it needs to be merged into a single geometry for the whole object. For our choice of pattern (but without using our smooth spline fit), adjacent cells share a single node and the structures of two cells may have different geometric parameters (radii) assigned to the shared node: these radii need to be averaged, affecting the properties of cells. In contrast, with our microstructures family, if the transition of materials is smooth, the change in radii will also be smooth and almost no averaging will be required. In fact, for square cells (𝛼 = 𝜋/2), the radius at the boundary will always be very close to 0.3 (see Figure 17), because of the way we fix 𝑝 5 in our map construction, as described in Section 6.2. Figure 32 shows an example where there is a significant mismatch between radii assigned to the shared node by two cell structures that have same material properties (𝐸 = 0.0315, 𝜈 = 0.75 and anisotropy measurements 𝑎 1 and 𝑎 2 in the order of 10 -5 and 10 -4 , respectively). After the merge, the new homogenized properties can be obtained by building a new square cell as shown on the right side of Figure 32.

The new material has significantly different properties, not being isotropic anymore (𝑎 1 = 0.028 and 𝑎 2 = 0.037). Interpolation accuracy. Finally, we show that our spline approximation is quite accurate. Figure 33 shows the absolute differences |𝐻 𝜈 (P (𝐸 𝑗 , 𝜈 𝑗 )) -𝜈 𝑗 | and |𝐻 𝐸 (P (𝐸 𝑗 , 𝜈 𝑗 )) -𝐸 𝑗 |, at points (𝐸 𝑗 , 𝜈 𝑗 ) not present in the originally sampled set: these are differences between actual homogenized material properties of the geometric structures obtained using our map P, and the values at which P was evaluated. Smoothness of the map P is essential for the interpolation of the relatively sparse set of values to be accurate.

Examples of manufactured structures. Finally, we applied our complete pipeline for a set of 2D examples similar to [Panetta et al. 2015], with a crucial difference being that we use irregular, boundaryconforming meshes optimized for rhombic shape of cells. The percell material property distribution is obtained using the optimization described in [Panetta et al. 2015]: we prescribe fixed Dirichlet boundary conditions on a part of the boundary, and target positions of another (target) part, optimizing the material distribution to obtain the desired deformation.

As we worked with 2D elasticity, we used 0.5ž sheets of closed-cell foam with small (< 0.2 mm) foam cell size to fabricate our examples using a laser cutter. We also experimented with other materials, as 1/8ž acrylic sheets, which allow for higher resolution when compared to foam (see pliers example in Figure 34). These examples have primarily illustrative purpose. We expect more practical applications to be accomplished by integration into new CAD design pipelines e.g, nTop Platform [nTopology 2020], which already supports a range of small-scale parametric structures.

We have demonstrated that it is possible to construct a family of geometric structures for rhombic shapes with a large range of angles, parameterized directly by material parameters of these shapes. The map from the material parameters to the geometric parameters is smooth and represented in a compact form as a set of 3D splines. Within explicitly defined bounds, any structure in the family is isotropic. While the material-to-geometry map depends on the base material properties, the transformation between different materials can be achieved with a simple analytic formula; for this reason, the proposed family can be viewed as independent of the base material properties.

We demonstrate that using this family improves the accuracy with which target material properties can be approximated using cellular structures manufactured from a single material.

This work is a first step towards explicitly defined material-togeometry maps of this type. We note that as an intermediate step in our construction we have obtained a map from arbitrary orthotropic, not just isotropic, properties to geometric parameters. Orthotropic materials may be useful in many contexts, and our work can be extended to this. Previous work [Milton et al. 2017;Ostanin et al. 2018] demonstrates that much broader, nearly optimal, range of parameters can be obtained by using more complex structure topology;

expanding the range further is another possible future direction of work.

This work forms a foundation for developing similar parametric families in 3D, since, although requiring more parameters and computational time, we do not foresee any fundamental problems in extending the construction of cell geometry families to 3D structures. As a proof-of-concept, we produced a smooth family of microstructures for cubic cells in 3D, following exactly the process described in Section 6.2, with 2D homogenization replaced with 3D. We used a cell geometry with cubic symmetry, guaranteeing orthotropy with same Young's modulus in all directions, which means having a three dimensional material space (𝑚 = 3). Initially, 9 parameters defining the geometry were considered (𝑝 = 9). This number was reduced to 7 and then to 5 with parameter elimination step (𝑝 ′ = 5). Then, we used PCA for the final reduction from 5 to 3. Figure 35 shows the initial and also the final coverage with our 3-parameter space, from which a nonlinear map for isotropic structures can be extracted. We note that the second component of the method, decomposing shapes into hexahedral cells of suitable shape, e.g., close to rhombohedra, requires much greater adaptation, although hex meshing methods can be used as a starting point.