Numerical Method for Direct Solution to Form-Finding Problem in Convex Gridshell

Traditional three-dimensional shell structures can resist external loads through their inherent shapes; however, if regular holes are made in the shell, with the removed material concentrated into the remaining strips, a structurally flexible gridshell can be achieved. Several spectacular architectures, e.g., Helsinki Zoo’s observatory tower and Centre Pompidou Metz, were manufactured with a network of one-dimensional beams, such structures serve both aesthetical and functional purposes in the civil engineering community [1]. Besides the construction of buildings in civil engineering [2–4], abundant applications in mechanical systems, e.g., micro/nanostructures [5–7], stretchable electronics [8,9], and bio-inspired patterns [10,11], employ gridshell as a major structural component in their design step to achieve specific functionalities. While the gridshells studied by Baek et al. [12,13] had joints that were free to rotate and twist, recent work by Panetta et al. [14] constrained the bending and twisting at the joints. This leads to non-trivial twisting and bending coupling between two rods at the joints, which can improve the robustness of the structure and increase the design space of the architectural shapes [3]. Computationally efficient numerical simulation tools for this class of structures can allow simulation-guided design and eliminate the need for painstaking trial-and-error prototyping.

In the computational mechanics community, modeling and simulation of thin elastic objects, e.g., rods and shells, are of sufficient general interest because of the preponderance of geometrically nonlinear deformation. Finite element method has been the most commonly used method in structural analysis over the past few decades [15,16]. Recently, discrete differential geometry (DDG)-based methods [17] are becoming popular in the computer graphics community to simulate the thin elastic structures, e.g., hair and clothes, due to the computational efficiency and the robustness in handling geometric nonlinearity, collision, and contact. Previous DDG-based methods have shown surprisingly successful performance in simulating slender structures, e.g., rods [18–22], beams [23], beam networks [24], ribbons [25,26], and plates/shells [26,27–30]. Gridshell, on the other hand, usually represents a curved surface comprised of multiple 1D elastic rods and differs from the traditional 1D rods or 2D shells. This leaves room for new numerical methods for accurate and efficient simulation of gridshells. Baek et al. first proposed a method based on discrete elastic rods (DER) to investigate the buckling instability and form-finding of gridshells [12] and found excellent agreement between experiments and simulations. A stiff spring is used in that framework to simulate the joint between two rods and the spring force is treated using an explicit approach. The joint between two intersecting rods is free to twist as well as rotate such that the twisting and bending coupling between two rods are not taken into account. This numerical framework was later used to study the elastic rigidity of hemispherical gridshells [13]. Numerical methods to capture the bending and twisting coupling at the joints either use a penalty energy between the neighboring material frames of two rods system [31,32], or a geometric constraint-based energy functional [14]. Finite element-based numerical methods to simulate this class of structures have also been introduced [33,34].

An even more intriguing feature of elastic gridshell is its form-finding process. Figure 1(a) shows three examples of convex 3D gridshells, whereas Fig. 1(b) describes the actuation process. In Fig. 1(b), the undeformed gridshell is planar and the extremities of the elastic rods fall on a closed curve, $G0$⁠. In order to actuate the gridshell, the end points of the rods are constrained to fall on a second closed curve, $G$⁠. The form-finding problem, i.e., the inverse problem in this case, calls for computation of $G0$ given the target 3D shape and the final boundary, $G$⁠. This transformation between the 2D planar structures and the complex 3D topologies by using the geometry and structural instability is of interest [5,6,11,34,35] and might lead many applications in mechanical systems [36–48]. Prior works on mechanically guided assembly of 3D structures range from macroscopic origami-inspired structures [49,50] to microscopic buckling of elastic ribbons attached to a pre-stretched substrate [6]. While a number of studies investigated the forward dynamics, we focus on a computationally efficient method to solve the inverse design problem of finding the initial planar shape with a given 3D target configuration. Reference [51] considered Chebyshev net theory to map a group of rods onto a given surface, e.g., human face, to design wire mesh. Prior works on the inverse problem include analytical solution to a pair of ordinary differential equations (ODEs) on the basis of Gauss equation [12], or numerical optimization coupled with physics-based simulations [11,14,30]. Recently, a genetic algorithm-based method [30] and an optimization-based simulation framework [14] have been introduced to study the form-finding problem in elastic gridshells; however, these methods require running the physics-based simulation numerous times in order to find the optimal solution, especially when a good initial guess is not available. As an example, the form-finding problem of a hemispherical gridshell in Fig. 1(a) may take approximately 10² generations with 5 × 10² individuals in a population in genetic algorithm, corresponding to 5 × 10⁴ forward simulations. The proposed method reduces this problem to a single forward simulation.

Here, we develop a numerical method for the inverse form-finding problem of gridshells. Different from above analytic and optimization methods that typically require numerous “forward” simulations to predict the deformation of the gridshell under various boundary conditions imposed by $G$ and $G0$⁠, this method implements a mechanics-based forward simulation of a gridshell draping around the target rigid shape under gravity. This single forward simulation can offer an excellent solution to the inverse problem. The simulation relies on a DER-based numerical framework, where both the rods and the joints are represented by the discrete elastic rod model. Discrete equations of motions, based on the balance of elastic and external forces, are solved to update the structural configuration with time. The main contribution of this paper is a numerical method for form-finding of convex gridshells based on contact [52]. In Fig. 1, we show several 3D configurations of convex gridshell structures as well as their corresponding initial planar boundaries constructed by the contact-based method described in this paper. The boundary of the 2D undeformed shape, $G0$⁠, can be almost exactly obtained by draping the elastic gridshell under gravity over the rigid 3D target surface. This calls for simulation of contact between the gridshell and the rigid surface and is handled via the modified mass method [27,53]. Discrete simulations are naturally suited to handle contact, which underlines the need for DDG-based methods in the study of form-finding of gridshells. The initial planar pattern of grid can be easily obtained by only running the physically based simulation once, which can significantly reduce the computational time when solving the form-finding problem.

Our paper is organized as follows. In Sec. 2, we discuss the DER-based numerical framework of gridshell simulation, with a focus on geometric decomposition and bending-twisting coupling at rotational joints. Then, in Sec. 3, we introduce a modified mass-based contact algorithm for the form-finding problem associated with gridshells. Finally, concluding remarks and avenues for future research are presented in Sec. 4.

In this section, we discuss the forward physically based simulation of gridshells. Gridshell is a type of structure that comprises multiple one-dimensional rods connected through joints. These joints may twist and rotate [14]. Here, in Sec. 2.1, we first briefly review the DER method for a single elastic rod, then extend this method to a numerical framework for the simulation of gridshells in Sec. 2.2. Finally, Sec. 2.3 presents two simple cases to demonstrate the twisting and bending coupling between two rods at the joint.

In the discrete setting of DER, shown schematically in Fig. 2(a), the rod centerline is discretized into N nodes: [x₁, x₂, …, x_N], and (N − 1) edges: [e₁, e₂, …, e_N−1], with e_i = x_i+1 − x_i, where i = 1, …, N − 1. Each edge, e_i, has an orthonormal adapted reference frame ${di(1),di(2),ti}$ and a material frame ${mi(1),mi(2),ti}$⁠; both the frames share the tangent t_i = e_i/‖e_i‖ as one of the directors. The reference frame is updated at each time-step through parallel transport in time, and, referring to Fig. 2(b), the material frame can be obtained from a scalar twist angle θ_i, see Ref. [21] for a detailed exposition of the DER algorithm. Nodal positions (total 3N) and the twist angles (total N − 1) constitute the (4N − 1)-sized degrees-of-freedom (DOF) vector, q = [x₁, θ₁, x₂, …, x_N−1, θ_N−1, x_N], of the discrete rod. Based on this kinematic representation, in the remainder of this section, we discuss the formulation of elastic energies, elastic forces, and the time stepping procedure of the rod solver.

In addition to the coincidence of two nodes, there is a non-trivial coupling between the twisting and bending modes at the joints, e.g., twisting rod 1 in Fig. 3(a) can cause rod 2 to rotate. Here, we consider the pin-joints with a specific constraint for rotations at the contact area. To account for this coupling at the joints, we decompose the basic gridshell element into four elastic rods in Fig. 3(b): the first two are the physical rods denoted as rod 1 and 2; the other two are linker rods with 3 nodes to model the joints. Hereafter, we use subscripts to denote quantities associated with the physical rods, e.g.,x_1,1 is the first node on the first rod, and superscripts when associated with linker rods, e.g., x^1,1 is the first node on Linker 1. Each rod can be simulated by the conventional DER method. A penalty energy can be used to account for the coupling between twisting and rotating at the joints. For the first linker and the first physical rod, the penalty energy is

In our numerical implementation, at every time-step, the equations of motions for the physical and linker rods are independently solved. This allows us to take advantage of the banded nature of the Jacobian matrix. The penalty forces in Eqs. (5)–(6) are then calculated and included as external force in the next time-step, i.e., the penalty forces are treated explicitly. An alternative to this approach of solving a number of smaller systems and subsequently bringing them together is to solve a large system, consisting of all the physical and linker rods, with an implicit treatment of the penalty forces. The large system would no longer have a banded Jacobian matrix since the Hessian matrix of the penalty energies would occupy non-banded entries within the Jacobian. A second alternative is to forego the use of penalty energies and treat the overlapping nodes (e.g. x_1,3, x_2,3, x^1,2, and x^2,2) and edges (e.g. θ_1,2 and θ^1,1) with the same degrees-of-freedom. For example, instead of using 3 × 4 degrees-of-freedom for the overlapping nodes in Fig. 3(b), we can introduce three degrees-of-freedom, x_joint, for the joint node and apply the sum of forces from all the four nodes onto the newly introduced single node. A simulation code developed for one method can be easily re-purposed to employ a different method. While solving extremely large systems, correct choice of the time integration scheme may depend on the computer memory as well as the degree of parallelism. A detailed comparison among the explicit method (used in the current study), implicit method, and the mapping method can be found in Ref. [30].

We use two simple demonstrations to show the coupling between two rods at the joint. In Fig. 4(a), we show the response of the basic element of a gridshell when one rod is twisted. The first twist angle of one rod, θ_1,1, referring to Fig. 4(a), is rotated with a prescribed angular velocity, ω = 10 rpm, such that θ_1,1(t) = ωt. Due to the two linkers at the joint, the centerline of the other rod rotates about the former one with the same angular velocity ω. This demonstrates the twisting coupling between two rods.

We now focus on the inverse problem of obtaining the undeformed 2D shape, given the target 3D shape. This method relies on draping the gridshell around a rigid body with the target geometry. In this section, first, the discrete gridshell algorithm is coupled with the modified mass method [27,53] to simulate the draping process; then, the procedure to obtain the initial boundary, $G0$ in Fig. 1(a), is described, accompanied by a number of examples.

Every time-step in simulation accounting for the contact with a rigid surface may require integration of the equations of motion twice. The first solve is the predictor step that determines if any node fell under the target surface. The optional second solve is the corrector step that is only necessary if any node was detected to fall through the rigid surface. A pseudo code of contact-based form-finding simulation is provided in Algorithm 1.

In Figs. 6(a1)–6(c1), the undeformed planar gridshells are located above the 3D rigid surfaces described by Eqs. (13). The elastic rods are symmetrically distributed about the x and y-axes in case of the hemisphere (17 × 17 grid) and the paraboloid (15 × 15 grid); for the hemi-ellipsoid, on the other hand, there are 11 rods along the x-axis and 15 along the y-axis. Note that the rod number for each case is determined by the size of the desired shapes, i.e. we want at least one node on each rod to contact the target surface. The planar gridshells are dropped under a gravity-type load that is large enough to drape the structure around the target rigid surface. In Fig. 6, gravitational acceleration of g = 9.81 m/s² was sufficient. Figures 6(a2)–6(c2) shows the deformed shapes of the gridshells. Parts of the gridshell are in contact with the rigid surface (located above the x − y plane) and the other parts remain suspended under gravity below the x − y plane. The suspended parts (i.e. nodes that fall below the minimum z-coordinate of the target rigid shape) are trimmed in Figs. 6(a3)–6(c3) to obtain the new extremities (first and last nodes) on each elastic rod. This describes the final boundary $G$ of the form-finding problem (also see Fig. 1). In Figs. 6(a4)–6(d4), the extremities upon trimming are mapped back to the initial planar gridshell, i.e., the planar shape is also trimmed to get rid of the suspended portions. This gives the initial boundary, $G0$⁠, of the gridshell. Then, the target 3D pattern described in Eqs. (13) can be obtained by moving the nodes on the extremities of the rods from the initial footprint, $G0$⁠, to the final boundary, $G$⁠.

The analytical solution to the initial boundary in case of a hemisphere [12] is also shown in Fig. 6(a4). For the cases of paraboloid and hemi-ellipsoid, the analytical solutions are not easy to derive and, therefore, we compare the planar boundaries obtained from the draping process and the ones found by genetic algorithm-based optimization [30] in Fig. 6(b4) and 6(c4). The good match indicates the correctness and the validity of the proposed method. Even when the solution from the process outlined in Fig. 6 is not accurate enough, it provides an excellent initial guess for conventional optimization algorithms.

For a physical understanding of this method, we consider the balance of forces. Each node in the simulation is balanced by three forces: (1) gravity, (2) contact force from the target rigid surface, and (3) elastic forces (primarily bending). This competition of forces yields a deformed shape that conforms to the target surface. On the other hand, in the “pop-up” fabrication process [12] of gridshell where the nodes on $G0$ are moved to $G$⁠, gravity and contact forces are replaced by forces acting on the extreme nodes (located on the boundary) by an external agent. Our results show that, surprisingly, the deformed shape remains almost the same despite substitution of gravity and contact with boundary conditions on a handful of nodes. Next, to demonstrate that the numerical method is robust against initial grid spacing, in Figs. 7(a) and 7(b), we show the hemispherical gridshell with different grid spacings. Here, the distance between two parallel rods are Δs = 4 cm (for Fig. 7(a)) and Δs = 5 cm (for Fig. 7(b)); the rod number is changed to 13 × 13 and 11 × 11, respectively, to ensure that each rod comes in contact with the target surface. As shown in Figs. 7(a2) and 7(b2), the initial planar grids match well with the analytical solution in both of these two cases.

Next, we highlight the computational efficiency of the presented contact-based numerical simulation of elastic gridshells. In Fig. 8(a), we plot the computational time as a function of time-step size, h, for three different cases. Here, the number of DOFs in all scenarios is fixed at ∼8000. The total simulation times are approximately 50 s (for hemisphere), 30 s (for paraboloid), and 20 s (for hemi-ellipsoid), separately. Then, referring to Fig. 8(b), we show the reliance of the computational time on the number of DOFs. The time-step size in this figure is set to be h = 0.1 s. Unsurprisingly, computational time dramatically increases as the number of DOFs increases. The simulations are performed on a single thread of Intel Core i7-6600U Processor @ 3.4 GHz. Overall, reasonable predictions can be obtained within one minute.

Via comparison with Refs. [12,30], the effectiveness of our proposed direct-contact method can be indirectly validated. The proposed method relies on the approximation that replacing the reaction forces at the boundary of the gridshell with a gravity-like load on all the nodes and constraints from the rigid target surface results in the same deformed shape of the gridshell. In the future, limitations of this approximation can be analyzed using theoretical mechanics. We also limited ourselves to convex surfaces with analytical solutions. Future directions of work can include surfaces with negative curvature and concave shapes; our contact-based form-finding method fails to handle negative curvature, e.g. hyperbolic surface. In this case, the final footprints are no longer planar. In addition, the proposed frictionless contact-based method would fail to fully drape around non-convex geometry. One possible approach in these cases is to use the proposed method to obtain an initial guess and then use this guess to start a more comprehensive inverse design process, e.g., genetic algorithm [30] and generative adversarial networks [58]. Another potential direction is to design multiple simple gridshells and then sew the solutions together to achieve more complex (potentially concave) gridshells. This concept was alluded to in Ref. [12].

The amount of twist along the rods in the examples explored in this manuscript is negligible. We did not find any significant effect of the coupling between bending and twisting on multiple rods at the joints. In the future, the role of the mechanics of the joints on the overall shape of the gridshell can be explored.

We introduced a numerical framework for the simulation of gridshells and solved the form-finding problem directly, without any numerical optimization. For the forward physical simulation, we first decomposed the gridshell as well as its joints into multiple elastic rods, such that each component can be treated using the well-established DER method. For the inverse problem of form-finding, we formulated a modified version of the discrete gridshell simulation algorithm by coupling it with the modified mass method to account for the contact between an elastic gridshell and the target rigid 3D surface. We showed that the gridshell, upon draping around the target shape, can be simply trimmed to directly get the initial planar boundary. A good match between the analytical solution and the contact-based result in case of a hemispherical target shape indicates the potential use of our method in form finding problems. Here, we limited ourselves within the convex surfaces with analytical solutions. The shape construction for arbitrary surfaces may need to introduce the frictional contact between the stretchable gridshells and target surfaces. We hope that our results and methodology will instigate future work on buckling induced mechanically guided assembly in physical systems (e.g., pop-up actuation of a planar grid to a target shape) from macro scale (e.g., domes in architecture) to micro-scale (e.g., controlled buckling of slender rods for stretchable electronics).

We acknowledge support from the National Science Foundation (Award # IIS-1925360) and the Henry Samueli School of Engineering and Applied Science, University of California, Los Angeles.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The authors attest that all data for this study are included in the paper.