1 Introduction
Metamaterials are structures of long-standing interest, as they induce material properties that differ from those of their constituent base material(s). Metamaterials often exhibit behaviors that are not found in nature, such as tuneable compliant, chiral, auxetic and non-reciprocal behaviors [Jenett et al.
2020; Panetta et al.
2015; Ou et al.
2018], and impressive strength-to-weight ratios [Qin et al.
2017].
The behavior of a given metamaterial is primarily governed by its
cellular architecture, which is the regular or random spatial arrangement of solid regions and voids used to fill a designated volume [Schaedler and Carter
2016]. The set of possible cellular architectures is uncountably large, even when we restrict our attention to, e.g., the space of regular cellular solids that fill a unit cube and tile periodically in
\(\mathbb {R}^3\) , as shown in Figure
1. Figure
1 also displays the many architectural elements that occur within our subspace, such as straight/curved beams, thin shells, and solid bulks. This breadth is powerful, as each architectural class offers a unique set of strengths [Bertoldi et al.
2017; Surjadi et al.
2019; Gibson et al.
2010].
However, these distinct classes complicate metamaterial design, because no existing representation is well suited for all structures. Generic representations like voxel grids can express any structure, but the specification is cumbersome and difficult to edit, since even simple changes (e.g., thickening a beam) trigger many independent voxel updates. Class-specific representations are often more practical, as each structure’s description is compact and editable. However, the underlying specifications are as varied as the structures themselves: trusses and beams are often given by graphs; solid bulks stem from constructive solid geometry (CSG) operations; and thin-shell cellular (shellular) structures use surface meshes or implicit functions. These specifications are not immediately compatible with one another, so it is difficult to explore over more than one class.
In some cases, it is difficult to explore variations even within a given class, as each structure requires a unique derivation. For example, many shellular metamaterials are based on TPMS, which are precisely defined as the integral of an Enneper–Weierstraß function. For ease of use, TPMS are commonly approximated by, e.g., level sets of an implicit trigonometric function. However, both function types are structure specific and difficult (if not impossible) to derive. This limits the set of structures accessible to engineers and, thus, available for exploration.
To alleviate these challenges, we propose a procedural graph representation that streamlines the process of metamaterial design for a wide range of common classes, including generic TPMS. Our representation is
specific to cellular metamaterial design, which allows us to capitalize on characteristics like symmetries while ensuring that our representation is equally suitable for all target classes. Moreover, our representation is
compact,
intuitive, and
easily editable, such that it is amenable to manual or automated design space exploration. Finally, it is
expressive enough to represent not only
known structures, as shown in Figure
1, but also
novel structures containing elements from one or more class(es).
The most critical aspect of our representation is the unified skeletal design space, which uses simple elements like lines and surfaces to capture the wide range of shapes found in metamaterials. Each skeletal element is constructed from a small set of simple, high-level specifications comprised of vertices and edges, a bounding volume, and the type of skeletal element (e.g., line, surface) to be instantiated. Each element is also annotated with a spatially varying thickness function that determines how it should eventually be thickened into a physically realizable volume (e.g., beam, shell, or solid bulk).
Although this design space naturally accommodates most of our target classes, the TPMS shellulars present a considerable obstacle. Our TPMS approach is grounded in mathematical principles, yet compatible with the design space described above and intuitive enough to be used by novices. A critical element of our approach is the
conjugate surface construction method (CSCM), which is “one of the most powerful techniques to construct minimal surfaces with a proposed shape in mind” [Karcher and Polthier
1996]. However, existing CSCM algorithms are largely inaccessible, as they require extensive domain expertise and human intervention. We embed the CSCM in an optimization loop to realize the first fully automatic version of this pipeline, making it accessible to a wider audience.
In summary, our contributions include the following:
•
a practical algorithm for TPMS via our extended CSCM,
•
a unified skeletal design space that compactly expresses the thickness-annotated skeletons for a wide range of metamaterials, including the five major classes in Figure
1, and
•
an intuitive procedural graph representation that facilitates the exploration and evaluation of novel structures.
We validate our approach by constructing hundreds of structures with diverse architectures. We also conduct a user study to verify our representation’s intuitiveness and ease-of-use. Finally, although we defer guided search strategies and physical property validation to future work, we show the potential of our approach by defining simple, random exploration schemes that automatically generate truss and shellular structures with a wide range of material properties.
3 Overview
As suggested by Figure
1, regular cellular architectures are well suited for skeleton-based design, as they are often highly symmetric structures derived from lines, surfaces, and easily reducible solid primitives. We could imagine modeling any such metamaterial with four simple steps: (1) build the skeleton for a small fundamental piece of the structure, (2) assign a spatially varying thickness profile
\(T(p)\) for each point
p of the skeleton, (3) apply any transformations (e.g., mirroring, rotation) required to fill the tiling unit, and (4) realize the final volumetric object according to
\(T(p)\) . To achieve this simple approach, we must address three challenges. First, we characterize the required skeletal elements and parametrize them in a concise, consistent manner. Then, we embed our boundary portion is restricted skeletal design space in a representation that captures the four-step approach above. Finally, we envision a user-centric modeling process that is easy and intuitive.
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).
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.
6 Procedural Graph for Metamaterials
Building atop our skeletal design space, we introduce a procedural graph that facilitates the full metamaterial design process, from initial shape specification to material property prediction. We summarize our node types in Table
1, including their inputs, properties, and requirements. In this section, we expand on the design logic and implementation for each node.
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)).
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.
7 Results
We implement our procedural graph in C++, with an OpenGL-based GUI for interactive design. We performed all experiments on Ubuntu 20.04 using an AMD Ryzan 5950X CPU (16 cores) with 64 GB RAM.
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.
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.
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.
9 Limitations and Future Work
Although our approach covers a wide variety of metamaterial architectures, it has several limitations. For example, our thickening operations do not necessarily preserve separation between features, which may cause issues for, e.g., interpenetrating lattices [White et al.
2021] or kirigami structures with cuts. Our thickening operations are also presently limited to simple sphere- or normal-extrusion with spatially varying thickness. We could provide more flexibility via other solidification routines, such as user-defined cross-section profiles for generating triangular prisms or I-beams. To mitigate areas of high stress concentration, we could also introduce blending methods that control the shape of junctions between abutting skeletal elements. The convex-hull-restricted blending approach [Panetta et al.
2017; Colli Tozoni et al.
2020] offers one such solution for beam structures, which could be integrated into our system by adding blending parameters on the vertices of line skeletons. However, this approach would need to be generalized to accommodate surfaces, solid bulks, and subtractive Boolean operations, as well as skeleton intersections that occur away from user-defined vertices.
Our work is also presently limited to a single cell of a regular metamaterial residing in a unit cube. We could expand this scope by exploring additional tilable cells such as hexagonal prisms or rhombic dodecahedral honeycombs. We could also integrate routines for stochastic porous structures like foams or for partitioning our graph into explicit, reusable subgraphs that could facilitate the design of complex hybrid structures or multi-scale designs [Zhang et al.
2021; Rao et al.
2019]. Finally, future work could consider ways to smoothly connect a set of distinct cells to facilitate the construction of larger volumes with, e.g., functional grading [Hu et al.
2022b; Al-Ketan and Abu Al-Rub
2021] or smooth transitions between structures of different classes (e.g., trusses to shells).
In a mathematical vein, we could also explore the limits of our CSCM approach w.r.t. the period problem discussed in Section
4.1. If robust, then our approach could facilitate the discovery of new examples or perhaps even a general solution for the problem of handle insertion. Alternatively, we could incorporate, e.g.,
minimal twin surfaces [Chen
2019] or
Willmore surfaces [Willmore
1965], which are a superset of minimal surfaces with
constant mean curvature. In conjunction with volume-preserving metrics, this could lead to a host of interesting new structures.
We also look forward to devising optimization schemes over our representation to permit the automatic discovery and interactive user-in-the-loop design of metamaterial structures with extremal properties. Toward this goal, our system would benefit from the development of a robust, principled GUI and expanded simulation capabilities, including nodes for, e.g., non-linear simulators, stress–strain curves, and tetrahedralization. It would also be interesting to explore whether we could perform analysis and/or property prediction directly on our graph representation to permit faster exploration and optimization techniques free of meshing- or simulation-related bottlenecks and sensitivities.