Procedural Metamaterials: A Unified Procedural Graph for Metamaterial Design

2.1 Cellular Architectures

The notion of metamaterials is very broad: even within graphics, HCI, and ML, recent research includes two-dimensional (2D) metamaterial sheets that are embedded in 3D [Konaković et al. 2016; Martínez et al. 2019; Signer et al. 2021]; interactive metamaterial mechanisms [Ion et al. 2016; , 2017; , 2019]; functionally graded structures for, e.g., spatially varying elasticity [Schumacher et al. 2015]; and multi-material composites with engineered properties [Gongora et al. 2021]. Although our approach may also apply to these domains, we restrict our focus to static 3D cellular metamaterials whose unit cells are regular (rather than random), tilable in \(\mathbb {R}^3\), and confined to a unit cube.

Trusses and Beams. Trusses and beams are often used for metamaterials, as they are easy to specify yet widely varied. The space of truss topologies alone is large enough to demand its own taxonomy [Zok et al. 2016], even before accounting for continuous parameters like vertex positions or thickness profiles over beams and their junctions. To explore this space, Jenett et al. [2020] and Frenzel et al. [2017] hand-designed several curved-beam structures exhibiting chiral, auxetic, rigid, and compliant behaviors. For diverse Poisson’s ratio and Young’s modulus values, Panetta et al. [2015] described 1,205 truss-based topologies using a tetrahedral decomposition of a cube and a graph over 15 possible nodes. Bastek et al. [2022] created 262 topologies by deforming and superimposing seven fundamental lattice units. Chen et al. [2018] used topology optimization to compute candidate structures, over which they fit graph templates to create parametrized truss families with extremal properties. Our representation is reminiscent of these approaches and trivially captures truss- or beam-based topologies that reside in a unit cube. However, we expand this powerful graph-based representation by porting it to architectures that are traditionally less amenable to exploration.

Solid Bulks. Solid bulks appear in many forms, including non-periodic spinodoid topologies [Kumar et al. 2020], foams [Ashby 2006], and open-cell porous structures [Tian et al. 2020]. However, we restrict ourselves to regular topologies, such as the topology-optimized, cubic-symmetric structures of Schumacher et al. [2015]. Chan et al. [2020] proposed another set of cubic-symmetric structures based on the level sets of 36 crystallographic space groups. Many CSG-like structures have also been hand- designed for specific properties like phononic bandgap [Muhammad and Lim 2021]. Our representation covers the vast majority of these regular structures.

Shellulars. Shellular units are often based on 2-manifold surfaces [Nguyen et al. 2016] or non-manifold surfaces like cubic lattices and honeycomb cells [Spadoni et al. 2014]. Several works also design functionally graded porous structures by spatially varying these base surfaces [Lu et al. 2014; Hu et al. 2022b]. Our representation captures a wide range of (non-)manifold surfaces, but the latter examples are out of scope, as we focus on individual cells.

TPMS Shellulars. TPMS are minimal surfaces that tile seamlessly along three mutually orthogonal directions in \(\mathbb {R}^3\). Non-self-intersecting TPMS are of particular interest, as they partition the surrounding volume into two or more distinct labyrinthine channels. Shellulars based on TPMS are smooth and uniquely suitable for additive manufacturing [Yan et al. 2021; , 2020; Hu et al. 2022b; Qin et al. 2017], making them ideal for tasks such as bone scaffolding [Ataee et al. 2018; Ambu and Morabito 2019], static mixing [Ouda et al. 2020], thermal energy management [Attarzadeh et al. 2022; Fan et al. 2022], and strength-preserving lightweighting [Qin et al. 2017]. TPMS are often created by Enneper–Weierstraß functions, specialized triangle meshes [Reitebuch et al. 2019], or level sets of implicit functions. The latter have enabled several TPMS explorer tools [Hsieh and Valdevit 2020; Al-Ketan and Abu Al-Rub 2021; Jones et al. 2021; Maskery et al. 2022] and methods for hybrid TPMS creation. Feng et al. [2021] and Khaleghi et al. [2021] construct hybrid TPMS by extracting isosurfaces from weighted sums of TPMS’ trigonometric approximations, while Chen et al. [2019] and Feng et al. [2021] use mesh Booleans. Zhang et al. [2021] propose hierarchical TPMS structures by assembling small-scale TPMS cells into larger ones. Although promising, these methods are limited by their reliance on structure-specific TPMS representations, which are difficult to construct and combine with other classes. Our approach addresses both limitations.

Alternatively, Akbari et al. [2020]; , 2022] construct surface- and truss-based approximations of TPMS from physical principles by using 3D graphic statics to explore the dual graphs of geometric form and static forces. To design hybrid TPMS, they specify the target labyrinth(s), i.e., channels of negative space surrounding the surface, and then construct a surface that partitions the volume accordingly. This scheme is powerful, but even domain experts struggle to relate labyrinth–surface pairs [Schoen 2008]. Our approach takes input describing the surface itself rather than its voids.

As discussed in Section 1, a critical element of our TPMS approach is the CSCM, which was introduced by Lawson [1970] in \(\mathbb {S}^3\) and then adapted to \(\mathbb {R}^3\) by Karcher [1989] with a discrete approach by Pinkall and Polthier [1993]. This method (described in Section 4) leverages a surface’s associate family to construct complex TPMS indirectly. We introduce the first automated, easy-to-use pipeline for this method.

2.2 Relevant Modeling Methods

Procedural Modeling. Procedural models “encapsulate a large variety of shapes into a concise formal description that can be efficiently parametrized” [Krs et al. 2021], which lends them to a variety of tasks including 2D textures and shaders [Cook 1984; Perlin 1985; Shi et al. 2020; Hu et al. 2022a] and virtual world modeling [Prusinkiewicz and Lindenmayer 2004; Smelik et al. 2014; Whiting et al. 2009]. Graph-based models are of particular interest, as they are widely used in practice (e.g., SideFX Houdini, Blender, Adobe Substance Designer), and they are amenable to performance optimization [Boechat et al. 2016] and the intuitive specification of edits and constraints [Krs et al. 2021; Michel and Boubekeur 2021]. Our procedural graph builds on these ideas toward concise, intuitive metamaterial design.

Skeleton-based Modeling. Skeletons have long been explored for efficient shape representation, as many volumetric shapes are well approximated by lower-dimensional structures [Blum 1967; Bærentzen and Rotenberg 2021]. Tagliasacchi et al. [2016] survey the rich skeletonization literature, which generally tries to reduce a volume to a skeleton. By contrast, we specify skeletons to construct a volume. Although most approaches approximate shapes with curve networks, Tagliasacchi et al. [2012] observe that some shapes are best represented by meso-skeletons containing a mix of curves and surfaces. As such, we develop a concise meso-skeleton representation for shapes appearing in metamaterial design.

3.1 Skeletal Design Space

Driven by our five target classes, our skeletal design space must accommodate straight/curved beams, planar/curved shells, and basic volumetric primitives such as cuboids and spheres. We posit that line and surface skeletons are sufficient to capture these shapes, when paired with a few annotations. To motivate and codify this design space, we briefly examine the needs of each shape category.

Beams. Straight and curved beams are well represented by line skeletons, which follow a path given by an ordered list of vertices. Curves can be created from concise vertex lists via, e.g., natural cubic spline interpolation. Thus, we need only introduce a “smoothness” flag to determine whether individual segments should be straight or smoothed. The final cross section of the thickened beams can be controlled by a spatially varying thickness profile over the line.

Shells. Shells are best represented by surface skeletons annotated with a spatially varying thickness profile, as above. Surface skeletons are also amenable to an ordered-vertex parametrization, as surfaces are frequently generated over a target boundary loop. For example, the widely studied Plateau problem spans a fixed boundary with a minimal surface, which locally minimizes surface area and has zero mean curvature everywhere [Harrison and Pugh 2016; Wang and Chern 2021]. Minimal surfaces can also be given by the free boundary problem, in which the input boundary is not fixed: each boundary portion is restricted to lie in its given plane, but its specific shape is inferred automatically as it “slides” along the plane to improve the metric [Karcher 1989]. For general TPMS construction, both boundary types are required. To capture them, we devise a pair of sliding solvers that take smoothness-annotated boundaries as input: smooth boundary portions are permitted to slide, while non-smooth portions remain fixed. The mixed minimal sliding solver is used when at least one fixed edge is present, as in Figure 2 (top). For fully sliding boundaries—which may otherwise degenerate—we introduce the conjugate solver, based on our extended CSCM. Finally, we introduce direct solvers to generate (not necessarily minimal) surfaces over fully fixed boundaries; here, non-smooth boundary portions remain straight while smooth portions are interpolated. All of our solvers assume disc topology (genus 0), as higher-genus surfaces can generally be decomposed. Critically, this also holds for TPMS: despite having genus \(\ge\)3 [Meeks 1975; Garbuz 2010], TPMS can be decomposed via their symmetry lines (see Section 4.1).

Fig. 2.

View Figure

Volumetric Primitives. When combined with simple thickening methods, lines and surfaces yield many basic primitives. For example, by offsetting the thickness along the skeleton’s normal direction(s), we can transform a line into a cylinder or a square- bounded planar surface into a cuboid. Similarly, by sweeping a sphere of the desired radius across our line/surface skeleton, we can create cylinders/cuboids with rounded ends/edges. This yields a sphere if the underlying line has length zero, as in Figure 2 (bottom).

Unified Design Space. In summary, line and surface skeletons capture our full set of target structures. Each skeletal element can be given by (1) a smoothness-annotated path over 3D vertices, (2) the skeletal type/solver, and (3) a spatially varying thickness profile over the element’s domain. For sliding solvers, we also provide (4) a set of bounding planes. Sections 4 and 5 explore the technical aspects of this design space.

3.2 Unified Procedural Graph

To facilitate the four-stage modeling process posed above, we create a procedural graph that unifies our skeletal design space with other pertinent operations. As shown in Figures 1 and 2, each graph node performs an operation such as vertex creation, line/surface inference, mirroring, or skeleton thickening. Each node also has properties that control its behavior. For example, an edge chain has a smoothness flag, and each surface node has a solver type and a thickness profile. By chaining and sequentially evaluating these nodes, we can form a variety of structures. We can also directly integrate nodes for performance evaluation, including, e.g., voxelization and simulations for material property prediction. Section 6 provides the full set of operations, along with their properties and implementation details.

3.3 User Design Process

To conceptualize a structure in our representation, users work backward through the stages of our procedural graph. First, users identify symmetries (e.g., mirrors, rotations) to reduce the structure to its smallest representative unit(s). After discounting thickness, they arrive at the structure’s fundamental skeleton (FS), which resides in the fundamental bounding volume (FBV). We support three scalable FBV primitives (cuboid, triangular prism, and tetrahedron) and custom FBVs. An FBV typically occupies a small part of the unit cube, though this need not be the case (Figure 2, bottom). Structures may also contain multiple FS, each residing in a unique FBV. Once each FS is identified, users can begin building the graph. The constituent vertices and edge chains for each FS are given by tracing its guiding/bounding path and classifying each portion as (non-)smooth. If a boundary contains both smooth and non-smooth portions, then multiple edge chains are required (Figure 2, top). The edge chains are then fed into a skeleton node, which can be transformed as needed to fill the unit cube and then thickened into an object. As the graph takes shape, users must verify that all nodes’ required conditions are met, e.g., for sliding solvers, all smooth boundary portions must lie on an FBV face. All such requirements are detailed in Section 6.

3.4 Outline

To build up a full understanding of our representation, we first explain our conjugate surface solver (Section 4), which is critical for many popular TPMS structures, and also our primary technical challenge. In Section 5, we discuss the remaining surface and line solvers, together with our spatially varying thickness annotations. The procedural graph is detailed in Section 6 and then evaluated for expressiveness, compactness, and intuitiveness in Section 7, which includes a user study.

4.1 Background and Overview

An associate family \(F^\phi (u,v)\) is a set of minimal surfaces that can be continuously transformed into one another by varying the scalar \(\phi\). Some well-known families transform a catenoid into a helicoid or a Schwarz P surface into a Schwarz D (Figure 4). For special pairs of surfaces \(S_1, S_2 \in F^\phi (u,v)\), the family can be parametrized as follows: (1) \(\begin{align} F^{\phi }(u,v) &= \cos \phi \cdot S_1(u,v) + \sin \phi \cdot S_2(u,v), \end{align}\) (2) \(\begin{align} &= \Re (e^{-i\phi } \cdot [ S_1(u+iv) + iS_2(u+iv) ])\ , \end{align}\) where \(\Re (z)\) returns the real part of a complex number z. As suggested by the complex formulation, the special surfaces \(S_1\) (at \(\phi =0\)) and \(S_2\) (at \(\phi =\frac{\pi }{2}\)) are said to be conjugate to one another. More generally, any two members \(F^{\theta }\) and \(F^{\theta +\frac{\pi }{2}}\) are conjugate to one another and every minimal surface is part of such a pair. These surface pairs have several pertinent properties:

Fig. 4.

View Figure

The symmetry lines noted in Proposition 4.1 3 are critical for minimal surfaces (particularly TPMS), as they allow the surface to be extended while preserving smoothness:

As hinted in Section 3, symmetry lines can also be used to decompose a complex surface into a fundamental patch P bounded by symmetries. Then, by Proposition 4.13, P has a conjugate \(\bar{P}\) bounded by the opposite symmetries. Moreover, it is possible to construct one surface from the other via a “conjugation” process [Karcher 1989; Pinkall and Polthier 1993]. Thus, the problems of solving for P and \(\bar{P}\) are equivalent. Since the term “conjugate” is overloaded—as the relationship between P and \(\bar{P}\) as well as the procedure that maps between them—we call P the primary surface and \(\bar{P}\) the dual.

The above equivalence is the foundation of the CSCM, which is most powerful when the dual \(\bar{P}\) can be solved for more easily than our primary target, P. For example, when P is a free boundary problem, its boundary \(\partial P\) contains only planar (sliding) symmetries; thus, \(\partial \bar{P}\) is a simple polygonal contour that admits a Plateau solution [Gray and Micallef 2007]. The CSCM easily recovers P by solving and conjugating \(\bar{P}\). The CSCM has fewer advantages when \(\partial P\) has mixed symmetries, as \(\partial \bar{P}\) is similarly complex. As such, we only apply the CSCM (via our conjugate solver) when P is a free boundary problem; other cases are deferred to our mixed minimal solver (Section 5.1.1).

The CSCM (Figure 3) contains several steps that are well established and tangential to our extension; for a full description, we refer readers to our Supplementary Materials. Here, we need only discuss the construction process for \(\partial \bar{P}\) (from Figures 3(b) to (c)), as this is the major limitation of existing CSCM algorithms.

By Proposition 4.1 1, the surface normal at each dual vertex \(\bar{v} \in \partial \bar{P}\) must match that of the corresponding primary vertex \(v \in \partial P\). Moreover, since P must be orthogonal to both bounding planes incident on v (Proposition 4.21), the normal at v must align with the intersection line on which v sits. By Proposition 4.1 2, the angle¹ between adjacent edges incident on \(\bar{v}\) equals the dihedral angle spanned by the planes incident on v. These facts prescribe everything about \(\partial \bar{P}\) except the length of each edge. Boundaries with four vertices are determined up to scaling due to the loop closure condition, but for N vertices, there are generally \((N-4)\) edge lengths that can be set arbitrarily.

Previous works require users to manually compute these edge lengths using an iterative process that requires an understanding of the relationship between \(\partial \bar{P}\) and P [Karcher and Polthier 1996]. This is especially taxing when aligning P with a given FBV, as the edges’ relative lengths determine the alignment—thus, they cannot be set arbitrarily. Improper edge lengths may also preclude a valid solution entirely: for the period problem, in which multiple portions of \(\partial P\) reside on the same bounding plane (see Table 2, “Wei’s genus 4”), most length assignments fail [Karcher and Polthier 1996]. Existing methods rely on a manual intermediate value argument that cannot readily expose the permissible values. Our approach identifies a valid solution while obscuring the CSCM details and reducing the required user expertise (as evidenced by the user study in Section 7.3).

Table 1.

View Table

Table 1. Node Types and Their Inputs. Brackets Around an Input Indicate that there can be Multiple Nodes of the Same Type Fed into the Given Node

Table 2.

View Table

Table 2. TPMS Boundary Loops. We Represent Several Well-Known TPMS using a Small Set of Vertices and a Single Boundary Loop. Column “S” Indicates the Name and Type of Each Structure: A Blue Label Indicates a Direct Surface, Green Indicates a Mixed-Minimal Surface, and White Indicates a Conjugate Surface. For Each Structure, We show the Boundary Loop and Bounding Planes (“vtx/edge”), the Generated Surface Patch (“patch”), and the Final Volumetric Structure (“obj”)

4.2 Our Edge Length Optimization

To address the limitation of existing CSCMs, we devise a fully automated solver (Algorithm 1) for free boundary problems given by a set of boundary planes B and an ordered edge loop L over these planes. We propose an optimization loop that finds edge lengths l for the dual contour \(\partial \bar{P}\), such that the final surface P conforms to B.

5.1 Surfaces

To complete our surface design space, we describe the mixed minimal and direct surface solvers. Figure 5 highlights the variation within our surface space by showing the result of each solver over identical annotated boundary loops and planes.

View Figure

Fig. 6.

View Figure

5.2 Lines

Line skeletons are represented by a sequence of 3D points, which can form any open, non-branching path or simple closed loop (if the endpoints are identical). Along non-smooth edge chains, neighboring vertices are connected by straight lines. Each smooth edge chain is interpolated to form a natural cubic spline that passes through the input points with C2 continuity everywhere. We use natural cubic splines, because they are simple, intuitive, and familiar, thanks to their widespread use in other modeling tools. To address concerns that are specific to cellular metamaterial design, we adjust the standard spline solver to permit (1) C2 continuity along closed loops and (2) curves that tile smoothly across a periodic boundary. We address ill-defined curves (with fewer than four vertices) by computing a quadratic spline (three vertices) or a straight line (two vertices).

5.3 Spatially Varying Thickness

To control how each skeletal element should be instantiated as a volumetric object, we include a spatially varying thickness function over the element’s domain. This function is given by sparse input: the user specifies the thickness at a few sample points, and then these values are interpolated over the full domain. The values are computed and stored w.r.t. a simple parametrization of the domain rather than actual vertex positions along the skeletal element. For lines, we use a 1D parametrization \(\gamma (t)\) for \(t\ \in [0, 1]\). The full thickness profile is linearly interpolated from the provided sample points. For surfaces, which have

disc-topology, we define thickness w.r.t. a UV-domain computed via the as-rigid-as-possible parameterization in libigl [Liu et al. 2008; Jacobson et al. 2018]. Then, as shown in the inset, we interpolate the thickness by treating each sample point (white dot) as a handle and computing bounded biharmonic weights over the domain [Jacobson et al. 2011].

6.1 Vertex and Edge Chain Nodes

A vertex node is the simplest node in our graph: it places a vertex at the 3D location given by its “position” property. Vertex nodes do not accept any input node connections. The next simplest node is the edge chain, which accepts a set of vertex nodes as input, and then instantiates a path over these vertices in the traversal order that is stored as a property of the edge chain node. As discussed, each edge chain also has a “smoothness” flag that is used by the subsequent line or surface node(s) to facilitate the variations described in Section 5.

6.2 Line Nodes

Our line node accepts a set of edge chains that are continuously traversable, i.e., the end of one edge chain is the start of another, such that they form a simple, non-branching path (open or closed). As discussed in Section 5, the shape of a line is determined by the smoothness of each constituent edge chain. A line node can accept any combination of smooth and fixed edge chains, as long as they are continuously traversable. For more complicated paths, multiple line nodes must be defined. The spatially varying thickness function for each line node is given as in Section 5.3.

6.3 Surface Nodes

Our surface node accepts a set of edge chains that form a simple closed edge loop, over which we instantiate a surface patch. Users first select a surface type from among mixed minimal, conjugate, and direct. This selection determines the algorithms used to interpret the boundary and solve for the final surface, as detailed in Section 5. It also dictates the requirements and properties for the surface node. For direct surfaces, the user need only provide the energy weight \(\alpha\) (see Section 5.1.2) and a non-degenerate, simple, closed boundary loop over vertices located anywhere in the FBV.

For conjugate and mixed-minimal surfaces, each sliding segment must lie on an FBV face. Mixed-minimal surfaces permit vertices anywhere on an FBV face (see Table 2, “Deformed H”). Conjugate boundaries are more restricted to satisfy the properties of Section 4.1. Specifically, all vertices must lie on FBV edges and form property-respecting configurations, as shown in the inset. The bottom inset is invalid, because the surface normal (blue) cannot align with the FBV edge as required.

We also have dual surface and associate family nodes to explore the families given by Equation (1). This lets us capture and explore intermediate surfaces like the gyroid (see Section 7.1).

6.4 Mirror, Transform, and Group Nodes

To build full translational units from a structure’s FS, we provide standard geometric transformations such as the mirror node. We also support translation, rotation, and scaling via the transform node, as well as the ability to combine a set of skeletal elements into one unit via the group node. Each of these nodes takes one or more “skeletons” as input, which may be a line, a surface, or the output of one of the transformation nodes above. The user can select whether the transformation is applied to a copy of the input skeleton (“doCopy”=true) or to the input skeleton itself (“doCopy”=false).

6.5 Object, CSG Boolean, and Voxel Nodes

To generate volumetric objects, our object node thickens each skeletal element based on its interpolated thickness profile. First, we instantiate a regular grid over a unit cube, which will store an indicator function that tracks whether each grid point is inside or outside of the object. The grid resolution is a property of the object node, along with the desired extrusion method. As suggested in Section 3, we support two extrusion approaches: spherical and normal. For spherical extrusion, we densely sample the skeleton and—for each sample point—we splat a sphere onto the indicator function grid, using a sphere diameter equal to the sample point’s prescribed thickness. For normal extrusion, the offsets are applied along the normal direction at each sample point. For surfaces, we offset each vertex by one-half of the desired thickness along each orientation (\(\pm\)) of its normal. For lines, we sweep a scaled circular cross section along the skeleton to create a tube. These methods are illustrated in Supplementary Materials. After extrusion, we compute the indicator function for each resulting geometry and union them together. To preserve tileability, we require identical function values for each pair of corresponding grid points on opposite boundary faces. Finally, we perform marching cubes on the indicator function [Lorensen and Cline 1987] to extract the object mesh.

We also provide a CSG Boolean node that supports union, intersection, and difference operations on a set of volumetric objects, as shown in Supplementary Materials. For efficiency, we compute Boolean operations on the underlying indicator functions for each object; then we enforce tileability and perform marching cubes as before.

Finally, we use the indicator function to implement our voxel node, which creates voxelized meshes suitable for property prediction.

6.6 Metamaterial Property Nodes

As a proof of concept for fully integrated design, we implement two property prediction nodes: material matrix and phononic bandgap.

The material matrix node computes the stiffness tensor of the volumetric structure given by our graph. We implement the homogenization approach of Panetta et al. [2015] to compute the equivalent stiffness tensor matrix C, , use hex not tetsassuming periodic boundary conditions and a linear elastic material over our voxel mesh. All linear systems were solved using Intel MKL Pardiso [Schenk and Gärtner 2004]. We also use C to compute and display the “material sphere” that illustrates the (an)isotropy of each structure (Figure 11(c)).

View Figure

View Figure

Fig. 9.

View Figure

Fig. 10.

View Figure

Fig. 11.

View Figure

Fig. 12.

View Figure

Fig. 13.

View Figure

Fig. 14.

View Figure

The phononic bandgap node predicts a structure’s ability to prevent the propagation of waves in certain frequency ranges, toward applications such as frequency filters, beam splitters, waveguides, and sound/vibration protection devices. The blocked ranges are known as bandgaps, and they are often the result of structural frequency-filtering mechanisms such as Bragg scattering and local resonances. We predict the bandgap using the approach of Åberg and Gudmundson [1997], which generates a set of dispersion curves showing the structure’s eigenmodes over varying wave vectors (Figure 12, light blue). We process these curves in search of horizontal bands through which no curves pass (Figure 12, grey rectangles). Each such area is a bandgap, as it indicates that the given frequency range has no viable transmission path through the structure.

Fig. 15.

View Figure

7.1 Conjugate Surface Construction

To evaluate our CSCM, we reproduce a variety of TPMS whose FS are given by free boundary problems, including the popular Schwarz P, Neovius, and Schoen I-WP structures. As shown in Figure 6, our surfaces exhibit strong agreement with the ground-truth TPMS given by the Enneper–Weierstraß functions integrated over a unit cell. We obtained each ground-truth surface via a publicly available Mathematica notebook [Weber 2018d; , 2018a; , 2018b; , 2018c], from which we extracted and uniformly rescaled exactly one translational cell of unit size.

Figure 6 also shows our accurate reproduction of the gyroid, which is one of very few TPMS that do not contain any straight or planar symmetry lines. This seems to render the gyroid incompatible with our method. However, the gyroid is a known member of the Schwarz P/D associate family: assuming that D occurs at \(\phi =0\), the gyroid occurs at \(\phi \approx 38\) [Karcher 1989]. Thus, we can construct a patch of the gyroid by creating the P and D structures via our CSCM, then interpolating via the associate family node.

We ensure that our optimization converges reliably for all of our examples, including those in Table 2 and hundreds of randomly generated structures (Figure 9). Even the most intensive structures in Table 2 , quantifyconverge within 80 iterations, as shown in Figure 7. Moreover, Figure 8 demonstrates that our algorithm’s output is stable w.r.t. the input vertex positions and initial edge lengths \(l^{\textrm {init}}\). The only difference is that more accurate initial guesses permit faster convergence, as evidenced by Tables 2 and 4.

Table 3.

View Table

Table 3. Varied Cellular Metamaterials. We Represent a Variety of Metamaterial Architectures, Including Some from Literature and Others of Our Own Creation

Table 4.

View Table

Table 4. Procedural Graph Statistics. For Each Structure, We List the Total Number of Nodes (#all) As Well As a Breakdown by Node Type: Vertex (#vtx); Edge (#e); Line (#l); Surface/Dual Surface (#s), Transformation/Mirror/Group (#T), Object (#o), and Voxel (#vox). We Also show Structure Type (“type”): “c” Refers to Conjugate Surface, “d” Denotes Direct Surface, “m” is Mixed Minimal Surface, and “l” is for Lines. Finally, We List the Time Cost to Evaluate Each Graph (t [s])

7.2 Representing Established Cellular Structures

We use our procedural graph to recreate structures from literature that span all of our target classes. As a simple test, we first re-implement the exhaustive truss-based exploration strategy of Panetta et al. [2015] and arrive at the same collection of 1,205 valid topologies in our representation. We also recreate a number of other structures found in literature. Tables 2 and 3 show the final geometry for a selection of our structures, which were created from scratch in our GUI, following the design process of Section 3.3. Table 4 gives a detailed list of nodes used for each structure, along with the computation time required to evaluate our full procedural graph. The minimal, median, average, and maximal number of nodes used is 12, 16, 19.1, and 34, respectively. This demonstrates that our graph is both lightweight and versatile. Furthermore, the computation time shown in column t demonstrates that our system is able to quickly produce the final structure from the sparse graph input, with median, average, and max time costs being 0.72, 2.31, and 21.91 s, respectively. Thus, our method can be used to create and modify metamaterials at interactive rates.

7.3 User Study

To evaluate the expressiveness, compactness, and ease-of-use of our approach, we invited 10 participants to model several metamaterials using our procedural graph. We sought participants with varying degrees of expertise in metamaterials, 3D modeling, and minimal surface theory. Many participants self-identified as a novice in one, two, or all three domains, but nobody identified as an expert in all domains.

Procedure. The study generally took 3 to 4 hours per user and contained five stages: (1) a pre-survey; (2) a brief presentation introducing the project goal and our representation; (3) a guided modeling session to practice conceptualizing/building structures using our representation; (4) an independent modeling session, in which the participant constructed six target structures; and (5) a post-survey. Supplementary Materials provides a detailed account of each stage, as well as the participants’ responses/results. As noted to our users, the study is primarily concerned with the intuitiveness and flexibility of the procedural graph representation, not the interactive tool.

Our primary experiment occurred in stage (4), as participants independently modeled six structures given by target 3D meshes. The target structures spanned all of our major classes, so we could examine our method’s overall expressivity and ease-of-use. To reduce undue burden on the participants, we asked them to focus on reproducing each target’s main structure rather than precisely inferring continuous parameter values for, e.g., thickness or vertex positions.

Main Results. All 10 users successfully reproduced all six structures, independent of prior experience. A subset of the structures are shown in Table 5, along with statistics about the time and number of nodes required to represent them. Of the 60 total modeling tasks, 56 (\(93\%\)) were completed in \(\le\)30 minutes and 45 (75%) were completed in \(\le\)20 minutes. In the post-survey, users also indicated high levels of confidence that they could implement unseen structures of the various classes in the future. Moreover, the overwhelming majority of users (90%) agreed that the process of modeling a diverse set of metamaterials would be easier/more intuitive in terms of our proposed procedural graph than it would be in terms any single other representation.

Table 5.

View Table

Table 5. User Study Results. For Each of the Six Modeling Tasks, We Show Two Randomly Selected User-Created Structures (All Structures are Shown in the Supplement). The Construction Time and Number of Nodes Used for Each Structure are Reported As the “avg (min, max)” Measured Across All 10 Participants

Curved shells presented the largest challenge, as users uniformly reported the lowest degree of confidence in—and highest degree of difficulty with—these structures. Nevertheless, all users succeeded in the curved shell tasks, and 90% of them expressed that it would have been more difficult or impossible to represent these structures using any other approach. This is particularly true of Hao Chen’s and Wei’s TPMS, as neither structure currently has a trigonometric approximation; this typically renders them inaccessible to designers, but our representation provided novice access to both. Moreover, the presence of a traditionally challenging period problem (see Section 4.1) in Wei’s genus 4 was a non-issue for our users; in fact, this structure required the lowest average modeling time of all. Overall, this user study confirms the expressiveness, compactness, and ease-of-use of our proposed representation, as even novice modelers can rapidly and faithfully realize a wide range of design intents.

8.1 Automated Structure Generation

Due to its compact form, our representation is conducive to automatic exploration. As a proof of concept, we devise simple exploration strategies for structures in two classes: straight trusses and shellulars. For a detailed explanation of these strategies, see the Supplementary Materials. Our truss exploration can generate a virtually unlimited number of structures due to the relatively unconstrained space for trusses. Valid shellular graphs are considerably more restricted, but our methods still generate hundreds of orthotropic and general asymmetric shellular structures. In particular, we obtained 1,000 direct structures in approximately 1 hour, 498 conjugate structures in 4 hours, and 500 general asymmetric structures in 3 hours. This exploration returned an enormous collection of unstudied direct shellulars and a mix of established and novel TPMS shellulars, as shown in Figure 9. The presence of established structures confirms that our representation encompasses critical regions of the cellular metamaterial design space, while the presence of novel structures indicates its potential for innovation. Our structures are also tilable by construction and (in all observed cases) physically realizable via inkjet-deposition additive manufacturing.² Figure 10 shows photographs for four of our fabricated structures, which each feature a \(3\times 3\times 3\) tiling of our unit cell, with a total size of 9 cm in each dimension and a minimum wall thickness of 1 mm.

8.2 Material Properties

Our work also presents opportunities for performance-driven design and shape optimization, as even our randomly generated structures exhibit a wide range of interesting material properties. We consider predictions for the stiffness tensor and phononic bandgap, using a base material with Young’s modulus \(E=1\) Pa, mass density \(\rho =1\) kg/m\(^3\), and Poisson’s ratio \(\nu =0.45\).