## Abstract

Elastic gridshell is a class of net-like structure formed by an ensemble of elastically deforming rods coupled through joints, such that the structure can cover large areas with low self-weight and allow for a variety of aesthetic configurations. Gridshells, also known as X-shells or Cosserat Nets, are a planar grid of elastic rods in its undeformed configuration. The end points of the rods are constrained and positioned on a closed curve—the final boundary—to actuate the structure into a 3D shape. Here, we report a discrete differential geometry-based numerical framework to study the geometrically nonlinear deformation of gridshell structures, accounting for non-trivial bending-twisting coupling at the joints. The form-finding problem of obtaining the undeformed planar configuration given the target convex 3D topology is then investigated. For the forward (2D to 3D) physically based simulation, we decompose the gridshell structure into multiple one-dimensional elastic rods and simulate their deformation by the well-established discrete elastic rods (DER) algorithm. A simple penalty energy between rods and linkages is used to simulate the coupling between two rods at the joints. For the inverse problem associated with form-finding (3D to 2D), we introduce a contact-based algorithm between the elastic gridshell and a rigid 3D surface, where the rigid surface describes the target shape of the gridshell upon actuation. This technique removes the need of several forward simulations associated with conventional optimization algorithms and provides a direct solution to the inverse problem. Several examples—hemispherical cap, paraboloid, and hemi-ellipsoid—are used to show the effectiveness of the inverse design process.

## 1 Introduction

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^{2} generations with 5 × 10^{2} individuals in a population in genetic algorithm, corresponding to 5 × 10^{4} 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.

## 2 Numerical Method

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.

### 2.1 Discrete Elastic Rods Method.

In the discrete setting of DER, shown schematically in Fig. 2(a), the rod centerline is discretized into *N* nodes: [**x**_{1}, **x**_{2}, …, **x**_{N}], and (*N* − 1) edges: [**e**_{1}, **e**_{2}, …, **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 3*N*) and the twist angles (total *N* − 1) constitute the (4*N* − 1)-sized degrees-of-freedom (DOF) vector, **q** = [**x**_{1}, θ_{1}, **x**_{2}, …, **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.

*E*, shear modulus

*G*, and isotropic circular cross section, the elastic energies—stretching, bending, and twisting—are given by [18,19]

*a*)

*b*)

*c*)

*A*is the area of the cross section,

*I*is the area moment of inertia,

*J*is the polar moment of inertia,

*ε*

_{i}is the stretching strain associated with the

*i*th edge, $\Vert e\xafi\Vert $ is its undeformed length, $\kappa i=[\kappa i(1),\kappa i(2)]$ is the bending curvature at the

*i*th node ($\kappa \xafi$ is the curvature in the undeformed configuration), τ

_{i}= θ

_{i}− θ

_{i−1}+

*m*

_{i}is the twist at the

*i*th node (

*m*

_{i}is the twist of the reference frame, details can be found in Refs. [19,21]), and Δ

*l*

_{i}= (‖

**e**

_{i}‖ + ‖

**e**

_{i+1}‖)/2 is its Voronoi length. The strain measures, i.e.,

*ε*

_{i}, $\kappa i$, and τ

_{i}, can be expressed in terms of

**q**(specifically,

**x**

_{i−1}, θ

_{i−1},

**x**

_{i}, θ

_{i},

**x**

_{i+1}). The case of non-circular cross-section can be included in the above formulation with minor changes [14,18,19].

*q*

_{j}(

*j*th element of the vector

**q**), the elastic forces (associated with nodal positions) and elastic moments (associated with the twist angles) are given by

*j*is an integer between 1 to 4

*N*− 1 and $Fjint$ is the

*j*th element of the (4

*N*− 1)-size elastic force vector.

*N*− 1 equations of motion and update the DOF vector

**q**as well as its velocity (time derivative of DOF) $v=q\u02d9$ from time-step

*t*

_{k}to

*t*

_{k+1}=

*t*

_{k}+

*h*(

*h*is the time-step size):

*a*)

*b*)

*c*)

**F**

^{ext}is the external force vector (e.g., gravity, damping, or reaction force due to contact) and

**M**is the diagonal mass matrix comprised of the lumped masses [54]. The Jacobian associated with Eq. (3a) necessary for Newton’s iteration can be expressed as

*i*and

*j*are integers between 1 to 4

*N*− 1, the mass (or angular mass) associated with the

*i*th DOF is

*m*

_{i}, and both the elastic internal force, $Fiint$, and external force, $Fiext$, are evaluated at

*t*=

*t*

_{k+1}. If the gradient of the external forces cannot be evaluated, $\u2202Fiext/\u2202qj$ can be ignored so that these forces are handled explicitly. Importantly, the Jacobian $J$ is a banded matrix and the time complexity of this algorithm is

*O*(

*N*) (the computational time linearly scales with the number of nodes [19]). This computational efficiency has motivated its application in the animation industry (e.g., hair simulation for movies) as well as its adoption in mechanical engineering [55–57].

### 2.2 Discrete Elastic Gridshells.

*j*th node on the

*i*th rod within the gridshell system is denoted as

**x**

_{i,j}. The twist angle of the

*j*th edge on the same rod is θ

_{i,j}. In Fig. 3(a), the two nodes,

**x**

_{1,3}and

**x**

_{2,3}, from two different rods overlap at the joint. A straightforward method to enforce the coincidence of two nodes at the joint is a linear spring-like energy of the form

*C*is the Lagrange multiplier. Its negative gradient, −(∂

*E*

_{c}/∂

**q**), is included as an external force in Eq. (3a). The Hessian of this energy can be trivially computed to aid the Newton’s method in the solution of the equations of motion. An alternative approach to enforce this condition will be discussed next.

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

*C*

_{1}and

*C*

_{2}represent the stiffness of the joint against the rotation and twist coupling. A similar penalty energy exists between the first linker and the second physical rod. At sufficiently high values of

*C*

_{1}and

*C*

_{2}, the rods at the joint cannot twist or rotate with respect to one another. We use

*C*

_{1}=

*C*

_{2}= 10

^{6}

*EI*in the current numerical investigation after a convergent study [30]. The external force and Jacobian associated with these energies can again be trivially computed. We should keep in mind that, when we decompose the basic element of the gridshell structure into two rods and two linkers, the mass and stiffness of the rods at the joint should not be double counted, e.g. the lumped mass at the joint node should be divided by four and then used as the mass associated with

**x**

_{1,3},

**x**

_{2,3},

**x**

^{1,2}, and

**x**

^{2,2}.

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].

### 2.3 Demonstration of Bending and Twisting Coupling.

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.

**x**

_{1,1}≡ [

*x*

_{1,1}(

*t*),

*y*

_{1,1}(

*t*),

*z*

_{1,1}(

*t*)], highlighted in Fig. 4(b), where

*a*)

*b*)

*c*)

*L*is the total length of the undeformed rod, and ω = 10 rpm. The location of the middle node on the first rod that falls on the joint is kept fixed with time to avoid rigid body motion. Due to these fixed degrees-of-freedom, the first rod rotates about the

*z*-axis with a prescribed angular velocity ω. Also, because of the coupling of bending between two rods at the joint, the second rod also rotates about the

*z*-axis at the same angular velocity.

## 3 Form-Finding

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.

### 3.1 Modified Mass Method.

*z*=

*f*(

*x*,

*y*). The position of a node,

**x**

_{i}(

*t*

_{k}), in a discrete gridshell structure at time

*t*=

*t*

_{k}, approaches the target rigid surface,

*z*=

*f*(

*x*,

*y*). If the rigid surface is not accounted for, the position of this node at the next time-step is, say,

**x**

_{i}′(

*t*

_{k+1}) ≡ [

*x*

_{0},

*y*

_{0},

*z*

_{0}]. In the time marching scheme of the simulation, if this node falls under the target surface so that

*correction*is required to move

**x**

_{i}′(

*t*

_{k+1}) onto the target surface along the the surface normal vector

**p**

_{n}‖ is the magnitude of the vector

**p**

_{n}, and

*ψ*is the angle between

**n**

_{xy}(the normal vector to the

*x*−

*y*plane), and

**p**

_{n}(the surface normal vector), as shown in Fig. 5(b).

*i*th node,

**x**

_{i}= [

*x*

_{i},

*y*

_{i},

*z*

_{i}], the updated form of Eq. (3a) is

*m*

_{i}is the lumped mass; $\Delta vipre$ is the prescribed change in velocity that can be obtained from the prescribed displacement, $\delta i$; $Fiint$ is the three-element elastic force vector on the

*i*th node; $Fiext$ is the external force vector; and the modified mass matrix is

**p**

_{n}is the constrained direction when free DOF=2; and

**p**

_{n}and

**q**

_{n}are the constrained directions when free DOF=1. Note that when a node is free, $\Delta vipre=0$, and Eq. (12) can thus be simplified to Eq. (3). If the node is fully constrained ($Si=0$), Eq. (12) reduces to $\Delta x\u02d9i(tk+1)\u2212\Delta vipre(tk+1)=0$ and the change in position (as well as the velocity) is enforced to take the prescribed value. In our case, we only constrain the node along surface normal,

**p**

_{n}, such that the number of free DOFs of the

*i*th node is 2. In our numerical implementation, we employ inelastic collision between the

*i*-th node and the target 3D surface, i.e. once the node is in contact with the target surface, we manually set its velocity to zero at the end of the current time-step.

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.

#### Modified mass-based contact simulation of elastic gridshell

**Require:***m* — Number of rods

**Require:***n* — Number of linkers

**Require:**$q(i)(0)$, where $i\u2208[1,m]$ — Initial DOFs of rods

**Require:**$q(i)(0)$, where $i\u2208[1,n]$ — Initial DOFs of linkers

**Require:***h* — time step size

**Require:***T* — total simulation time

**Ensure:**$q(i)(T)$, where $i\u2208[1,m]$ — Final DOFs of rods

**Ensure:**$q(i)(T)$, where $i\u2208[1,n]$ — Final DOFs of linkers

$k\u21900$

$tk\u21900.0$

**while**$tk\u2264T$**do**

$tk+1=tk+h$

$k\u2190k+1$

**for**$i=1$ to $i=m$**do**

$solved\u21900$

**while** solved = 0 **do**

Update DOFs in *i*-th rod ($q(i)$)

$solved\u21901$

**for** each node on *i*-th rod **do**

$[x0,y0,z0]\u2190$$x,y,z$ coordinates of the node

**if**$z0<f(x0,y0)$**then**

Constrain the normal direction of $xi,j$

$solved\u21900$

**end if**

**end for**

**end while**

**end for**

**for***i* = 1 to *i* = *n***do**

Update DOFs in *i*-th linker (**q**^{(i)}) based on Eq. (3*b*)

**end for**

**end while**

**return q**_{(i)} (*T*), where $i\u2208[1,m]$.

**return q**^{(i)} (*T*), where $i\u2208[1,n]$.

### 3.2 Initial Boundary From the Draping Method.

*s*= 1.2 m, rod radius

*r*

_{0}= 1 mm (and, therefore, second moment of inertia, $I=\pi r04/4$, polar moment of inertia, $J=\pi r04/2$, and cross-sectional area $A=\pi r02$), Young’s modulus

*E*= 1 MPa, shear modulus

*G*=

*E*/3 (assuming incompressible material), material density

*ρ*= 1.0 g/cm

^{3}, distance between two parallel rods Δ

*s*= 3 cm, discrete edge length $\Vert e\xaf\Vert =5mm$, and distance between planar X-shell and 3D surface top (choosing arbitrarily) is

*H*= 5 cm. As long as the rod can be assumed to be soft enough and inextensible, these parameters do not significantly influence the actuated shape of the gridshell [12]. The geometries of the target surfaces are given by

*a*)

*b*)

*c*)

*R*

_{h}= 0.2 m (for hemisphere);

*R*

_{p}= 0.2 m and

*H*

_{p}= 0.12 m (for paraboloid);

*a*

_{e}= 0.2 m,

*b*

_{e}= 0.15 m, and

*c*

_{e}= 0.12 m (for hemi-ellipsoid).

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^{2} 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.

### 3.3 Computational Time.

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.

### 3.4 Limitations.

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.

## 4 Conclusions

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).

## Acknowledgment

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.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

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.