Form-finding of shells containing both tension and compression using the Airy stress function

Pure-compression shells have been the central topic in the form-finding of shells. This study investigated tension-compression mixed-type shells by utilizing a NURBS-based isogeometric form-finding approach that analyses Airy stress functions to expand the possible plan geometry. A complete set of smooth version graphic statics tools is provided to support the analyses. The method is validated using examples with known solutions, and a further example demonstrates the possible forms of shells that the proposed method permits. Additionally, a guideline to configure a proper set of boundary conditions is presented through the lens of asymptotic lines of the stress functions.


Introduction
Shell structures, such as concrete, masonry, and metal shells, steel or timber gridshells, cable nets, and fabric structures, are elegant and light weight ( Figure 1). Their dominating load-carrying action is membrane action, which is a combination of tensile, compressive, and shear stresses or forces acting in a plane tangential to the surface of the structure. Shell structures containing only tensile or only compressive stresses cannot have unsupported boundaries with overhang but their boundary curves must draw curves shrinking inward ( Figure 2).
To allow for overhangs, a mix of tension and compression stresses is needed (Figure 3). Whenever there are compressive stresses, some bending stiffness is required to avoid buckling. Bending stiffness can also be added to allow for the shape of the shell to deviate from the pure membrane action shape ( Figure 4).
During the design of shell structures, membrane action is commonly secured using structural form-finding whereby a geometry is determined in such a way that no bending action is needed for the transfer of the dominating load. Most numerical methods for form-finding of shells simulate a physical model that might involve hanging chains or fabric that will be then inverted and form a compressive structure. The models are usually discretized using line elements connected to mutual nodes where external loading also is applied, and the geometry is iteratively updated based on the previously computed stress state. As a result, such methods are restricted to tensile only or compressive only shell structures.
In this study, we revisit, clarify, and develop the numerical form-finding method presented by Miki et al. 1 and the smooth version Graphic Statics tools discussed by Miki et al. 2 The method takes as input a pre-determined membrane stress state specified using the scalar-valued Airy 3 stress function φ defined on a plane. With φ given, the shape of the shell z is solved such that the vertical equilibrium equation is satisfied 4,5 given an appropriate set of boundary conditions. ρ is the upward vertical loading per unit plan area, which may be constant or may depend on the shape. Equation (1) is equivalent to equation (f) on page 462 in Theory of Plates and Shells, 6 whose authors attribute the first use of the Airy stress function in this manner to Pucher. 7 In Chapter 10.3 and the subsequent chapters of Theory and Practice of Membrane Shells, Csonka 5 discusses stress functions for shells on polygonal plans in pure compression and star-shaped plans with a mix of compression and tension. The stress function is reintroduced in the Appendix D. Some shell structures consist of several smooth surfaces joined via sharp kinks where concentrated forces arise, for example, several of Félix Candela's concrete shells ( Figure 3) and Jorn Utzon's Sydney Opera house (Figure 4, left). Such force concentrations can be represented in the Airy stress function by narrow regions with a high degree of curvature. If these narrow regions are taken to the limit, "folds" form on the stress function. 8 Then, the  piecewise smooth stress function can be seen as a hybrid between a smooth stress function and a stress polyhedral. 9 In graphic statics for bar frameworks such as planar trusses, stress polyhedrons can be used to establish a discrete force diagram, [10][11][12][13][14] which has a reciprocal relation to a discrete form diagram. Graphic statics originated in the 18th and 19th century [15][16][17][18][19][20] and have recently regained attention [21][22][23][24][25][26][27][28][29][30] with generalization into higher dimensions [31][32][33][34][35][36][37][38] and applications for the form-finding of she lls. 14,[39][40][41][42][43][44] Attempts have been made to establish a similar relation between a piecewise smooth stress function and the shape of the shell, 1 which require the computation of the Christoffel symbols of the second kind. To avoid the often tedious computation of the Christoffel symbol, a reciprocal relation between the piecewise smooth continuous stress function, a pre-computed continuous force diagram, and continuous form diagram (i.e. plan geometry) has been established. 2 If φ has a positive Gaussian curvature, it represents a state of pure compression or pure tension, whereas if the Gaussian curvature is negative, it represents a state of compression and tension. Thus, the method, in principle, allows for the form-finding of shells with any membrane stress state. The stress functions considered in this paper all have negative Gaussian curvatures, resulting in mixed tension-compression shell structures. For such cases, equation (1) becomes a hyperbolic partial differential equation, 4 which is often challenging to solve analytically.
In line with the methodology of Miki et al., 1 equation (1) is solved numerically using isogeometric analysis with NURBS surfaces as finite elements. [45][46][47][48] In isogeometric analysis, the computational model and the geometry model are the same, eliminating the need to discretize the smooth shell surface into flat panels or a network of straight-line elements. Without the need to mesh the surface, designers can concentrate on the form in early design stages, leaving discretization of the surface into structural elements such as panels and bars for later stages.
In this study, we only consider isogeometric models consisting of multiple NURBS surfaces, with their control points on the boundary edges shared at their intersections. This strategy makes the computational method simpler because there is no need to consider additional  boundary terms to glue surfaces. More general isogeometric models may contain B-rep representations, where each surface is defined by a combination of base surfaceusually a NURBS surface-and arbitrary boundary curves and holes. In this case, extra "weak" boundary terms need to be incorporated to stitch surfaces together. Nitsche's method is known as one of the standard methods to do this. [49][50][51][52] However, we leave it to the future work and only discuss simple cases in which each surface is rectangular in its UV-space and adjacent surfaces simply share control points.

Contributions
In this study, we breakdown the original method presented by Miki et al. 1 to make it more accessible and reproducible. While the earlier paper 1 restricted its attention to shells of pure compression, we concentrate on shells containing both compression and tension stresses. The contributions in this paper are as follows: for tension-compression mixed type problem is addressed (Section 4) through the lens of asymptotic lines. 5. Introduces an Airy stress function (Section 5) that can take an arbitrary plan geometry. An example problem is provided to demonstrate the possible forms enabled by the proposed stress functions.

Theoretical background and numerical method
In this section, starting from the traditional discrete graphic statics tools, their smooth versions are derived. Then, a NURBS-based isogeometric form-finding method is introduced. The method solves a linear system of equations when the surface area is evaluated on the projected plane. Its nonlinear version that accounts for the surface area accurately is also presented. Figure 5 shows an example of a force and form diagram pair. If each internal node in the form diagram is balanced by axial forces acting along the incoming lines, a reciprocal force diagram consisting of closed polygons can be constructed. If the form diagram is not in equilibrium, one or several polygons in the force diagram can not be closed. In the force diagram, the lengths of the edges of the polygons represent the magnitude of the corresponding force and are rotated by 90 compared to the corresponding line in the form diagram. 53 The reciprocal relation is such that a point, a polygon, and a line in a form diagram maps to a polygon, a point, and a line in the force diagram, respectively. The existence of a force diagram is equivalent to the existence of a planar-faced polyhedron (i.e. stress polyhedron 9 or Airy stress polyhedron) whose planar projection is the form diagram. [10][11][12][13][14] If the plane of a face in such a polyhedron is expressed as

Discrete force diagrams and Airy stress polyhedron
its corresponding point in the force diagram is located at a b , ( ) . At the same time, the "normalized" normal vector of the same face is given by Thus, a force diagram can be obtained by flattening all normal vectors of the polyhedron and connecting the points if their corresponding faces are adjacent ( Figure 6).

Continuous force diagram
For a stress function φ φ = , x y ( ) , the normalized normal vector is given by When partial differentiation is performed to a function multiple times, the order of the differentiation is insignificant. This often-overlooked condition, that is allows any choice of φ automatically satisfies the horizontal equilibrium equations, Similarly, we choose any pair of functions X and Y that satisfy For such a pair, the stress tensor is given by Thus, from equations (5) and (9), which are the components of the normal vector in equation (4). Hence, in analogy with the relationship between the normal vector of a stress polyhedron and the discrete force diagram, X Y , ( ) is the continuous force diagram corresponding to a continuous stress function. 2 Note that the existence of a stress function φ is not necessary to construct a force diagram ( , ) X Y . In fact, it is sufficient to require the conditions given in equation (8). This condition is known as "curl-free" and is essentially an integrability condition of gradient vectors. It is also explained in a structured manner in exterior calculus. 55 Curl-free means that the deformation gradient is symmetric and consequently has no components of rotation. The force diagram can be considered as a deformed shape of the form diagram, and the local deformation at each point can be described by a deformation gradient. The deformation gradient maps a unit circle in the form diagram (Figure 7(a)) to an ellipsoid in the force diagram ( Figure 7(b)). Without any rotational deformations, there are two orthogonal directions where the local deformation has stretches only (the red lines in the figure) and these coincide with the principal directions of the stress tensor; but the first and the second directions are swapped (Figure 7(b)).

Introduction of curvilinear coordinate system
So far, an orthonormal coordinate system x y , ( ) has been used. However, for shells with arbitrary boundary shape, a curvilinear coordinate system u v , ( ) is often preferred. The mapping between the two systems can be described well by the terms of differential geometry. 56 Books such as Ciarlet. 57-59 are a good reference because they discuss differential geometry from the viewpoint of the theory of elasticity.
be two symmetrical maps. Then, the maps ( , ) ( , ) u define a pair of parametric surfaces representing the form and force diagrams that can be expressed as Because the inverse of ( , ) ( , ) u v x y → exists (or we assume singularity/overlapping does not occur in our surfaces), exits. If it is given by equation (10), it can be computed by where The deformation tensor F on a u v , where g u and g v are the dual base vectors of x( , ) u v and G u and G v are the base vectors of X ( , ) u v . For details on the computation of the base vectors and the extraction of the components F xx , F F xy yx = , and F yy , see the Appendix A. The deformation gradient tensor F given by equation (16) is equivalent to the one given by equation (11), and both are symmetric. Therefore, which is an equivalent condition to the curl-free condition. A special case of the curl-free condition is when Figure 8), which is an interpretation of the perpendicular condition of the discrete force diagrams to the continuous ones.
In general, a computation of the second derivatives of the stress function on a curvilinear coordinate system requires the computation of the Christoffel symbols of the second kind, Γ ij k , describing the curvature of the isocurves of a parametric surface (see e.g. Green and Zerna, 4 Eisenhart, 56 and Ciarlet [57][58][59] ). However, through the above formulation, the first differentiation is computed on ( , ) u v , yielding a vector on the ( , ) u v system, which is transformed to the ( , ) x y system using equation (14). Then, a differentiation of ( , ) X Y with respect to ( , ) u v is computed using equation (16), but ( , ) X Y is already a vector in the ( , ) x y system. Therefore, there is no need to compute the Christoffel symbols. Moreover, the precomputed force diagram, ( , ) X Y , assists the evaluation of the crease curves in a piecewise smooth stress function (Section 2.7). This is where the proposed method takes advantage of the original method by Miki et al., 1 substantially simplifying the computations.

Linear method
The horizontal equilibrium of the shell is ensured by equation (7). To complete the equilibrium, the shape of the shell z must be determined so that the vertical equilibrium given by equation (1) is fulfilled. Using the principle of   virtual work, a finite number of simultaneous equations can be obtained as an approximation of the vertical equilibrium equation. This approach is equivalent to the Galerkin method, 61 and the details are provided in the Appendix B. Assume the vertical loading ρ is constant, so that the solution z of equation (1) becomes a linear problem. Later on, in Section 2.5, ρ will be taken as the selfweight of the shell. Then, ρ is no longer constant but depends on the actual area of the shell, turning equation (1) into a non-linear problem.
We introduce a curvilinear coordinate system ( , ) ξ η on the shell that is different from that of the force and form diagrams u v , ( ). Let h ξ and h η be the base vectors and h ij be the first fundamental form of the ξ η , ( ) system, and we denote the reference states of the properties of the shell that are not affected by differential operators by a bar − . Then, from the principle of virtual work, we get: is a small area element, S ij the components of the second PK stress tensor re-evaluated on ( , ) ξ η , and ρ a constant loading per unit area, all measured in the "plan" of the form diagram. The Einstein summation convention applies, and the indices i and j are either 1 or 2, where 1 corresponds to the ξ -direction and 2 corresponds to the η -direction. Details on the extraction of the stress components on the ( , ) ξ η system are given in the Appendix A. If the shell is represented by multiple finite elements or parametric surfaces (e.g. triangles, quadrilaterals, or NURBS-surfaces) defined by n independent control points, all functions defined on the shell, such as h ij and z , become functions of the control point coordinates. By packing the coordinates of these control points in a 3n-dimensional column vector r and replacing the variation of functions δ f by ∇f δ r , where δ r is a 3n-dimensional arbitrary column vector and the gradient operator is r r  n (20) equation (18) becomes Since the horizontal equilibrium is automatically ensured, solving equation (21) in terms of z is enough to obtain a complete equilibrium. In case equation (21) is solved in terms of x and y , the plan of the solution should match the form diagram.
Note that equation (21) is linear in r . By distributing integration points and performing a numerical integration (see Section 2.8), a system of linear equations is obtained. With a sufficient number of fixed points given, this system of equations can be solved easily by computing an inverse matrix. Although presented quite differently, Pauletti and Pimenta 62 have already pointed out that having a constant second PK stress tensor results in a linear system of equations.

Nonlinear iterative method
If ρ is the self-weight of the shell, equation (18) only gives an approximate solution. The exact solution is given by is the actual element area measured on the shell. Since da depends on h ij which is a property of the sought solution z, equation (22) is a non-linear problem requiring an iterative solution process.
In each iterative step, the residual forces of the system are computed by where ∇ ∇ ∇ 2 = T and ∇ 2 h ij is a constant matrix. Thus, the matrix S h ij ij ∇ 2 does not change during the iterations. Then, with h ij , the inverse of h ij , and using the relation where D is a generic differential operator that can be either a δ , , d or ∇ operator, the system of equations to solve in each step becomes Note that in the first step, solving equation (26) is the same as solving equation (21). The convergence of the solution is often improved when the second term in the parentheses of equation (26) is omitted, that is, it is equivalent to repeating the linear method but with the updated loading. Another good strategy is to scale down δ r by 0.5, which contributes to the stability of the method.

Point and area supports
A control point can be easily fixed to a specific height by eliminating the control point parameters completely from the system of equations. However, if NURBS surfaces are used as finite elements, the control points do not pass through the surface except at the corner points. Thus, to constrain an arbitrary point of such a surface to a specific height, a linear spring can be used.
The elastic energy of such a spring is where z is the prescribed height of the point, z( ) r is a linear function that converts the control points to the height of the point, and k is a penalty factor.
When such springs are incorporated into the system, the equilibrium of the entire system is springs membrane term where the spring term is linear in r and the membrane term is discussed in Section 2.4 or Section 2.5 depending on the loading condition. Similarly, area supports defined by closed curves can be imposed on the shell by the use of a set of linear springs distributed along the curve. This approach is used later on in Section 5 when discussing some examples.

Piecewise smooth stress functions
Whenever the shell is made up of several smooth surfaces joined by sharp kinks along their intersection curves, the stress function must be a piecewise smooth function. For such shells, the equilibrium is sustained by concentrated forces acting along the kinks on the shell. The axial force along a kink in the shell surface can be calculated by measuring the "jump" of the normal vectors on either side of the corresponding kink in the stress function. The same entity can be found by measuring the "jump" in the force diagram. Denoting the "jump" by L , the second PK stress tensor along a kink in the surface 1 becomes and the principle of virtual work adding up to equation (18) or equation (22) is is a length element on the intersection curve as seen in the plane of the form diagram.

Numerical integration
In this work, a standard Gauss quadrature was used to numerically compute the integrals over the surfaces and curves. 45 Because the stress components are calculated from a stress function that is not a polynomial in general, a sufficient number of integration points is hard to estimate. Hence, the number of integration points was controlled by an integration multiplier, denoted by NI.
We denote the number of control points in the u and v directions by NU and NV , respectively, and the degree of the NURBS surface by Dim . Then, a NURBS surface is typically divided into (NU-Dim)×(NV -Dim) smaller patches, and each of them is dependent on (Dim +1)×(Dim +1) control points. 63 In this study, (Dim ×Dim ×(NI)) 2 Gauss integration points were distributed on each patch (i.e. n = Dim ×NI in the commonly used Gaussian quadrature rule tabulation is used). NURBS surfaces with Dim = 3 and NI = 1 2 or work for most cases, resulting in n or = 3 6 for the Gauss quadrature rule. In the worst case n = 4 is enough.

Convergence study
A variational operator δ has infinite degrees of freedom (DoF) by definition. The Galerkin method restricts the variational operator to a finite number of DoF equal to the DoF of the finite element model. If a solution exists in the original PDE, by increasing the DoF in the system, the numerical solution should converge to the exact solution.
In this study, the convergence towards the exact solution is evaluated by studying the total energy of the system given by The elastic energy of the springs is omitted because they will have high energy if used as hard constraints.
It is known that there are incompatible boundary conditions (BCs) that make hyperbolic PDEs ill-posed (i.e. solutions do not exist). In such cases, numerical solutions hardly ever converge. Since a mix of tension and compression stresses results in a hyperbolic PDE, attention must be paid to such situations since the proposed numerical method returns a solution regardless of the problem is well-posed or not. A guideline to choose proper BCs is provided in Section 4.

Benchmark problems
In this section, we present two simple benchmark problems. We only consider Dirichlet type BCs and BCs involving higher order derivatives will not be discussed. With more complex problems, by taking the advantage of the fact that the NURBS surfaces have higher order non-zero derivatives, methods that require computations of higher order derivatives, such as the least squares method, may be considered. However, we leave it to future work.

Hypar
Consider a simple stress function that is defined on a square that spans ( 1, 1) − − , (1, 1) − , (1,1) , and ( 1,1) − . In Figure 9, the purple surface corresponds to φ , and the gray flat planes are the stress-free region outside the structure apart from the reaction forces seen as the kinks at the top.
The force diagram of φ is given by ( , ) = ( , ) X Y x y − , which is curl-free. Figure 10 depicts the force diagram and the "jumps" of the normal vectors between the smooth pink surface of φ and the surrounding flat gray planes. By adding edge beams, the shell can be supported by only two points, and the magnitudes of the "jumps" represent the axial forces in the beams.
Substituting equation (33) into equation (1) gives for which there exists a known solution, The solution satisfies the boundary equilibrium between the edge beams and the shell. Note that equation (35) is only a particular solution. With more complex BCs, the solution is a sum of the particular solution and a general solution (a solution of equation (34) where the right-hand side is set to zero), which satisfies the BCs. While a general solution always exists in an elliptical problem, in a hyperbolic problem, there exist incompatible BCs, where general solutions do not exist. In general, a hyperbolic PDE only accepts compatible BCs, otherwise the problem becomes ill-posed. This issue is further addressed in Section 4.
This problem can be used as a useful benchmark problem, and the proposed method was tested. The green surface in Figure 9 represents the solution obtained using the proposed method with ρ = 4.0 as load. The height of the obtained green shell is 1/ 4 of ρ , which matches equation (35).
The blue surface in Figure 11 is a plot of the left-hand side of equation (1) multiplied by 1/ ρ obtained numerically with a NURBS surface with 32 32 × control points, Dim = 3, and an integration point multiplier NI = 1. It clearly shows a drop of −1 from the orange surface representing zero loading. This indicates that the method provides the sought solution. Figure 12 shows a convergence analysis of this example. As shown in Figure 12, the total energy remains flat even if the number of control points was increased. This is because the hypar can be accurately represented with even a NURBS surface of 2 2 × control points and with Dim = 3.

Bowl
Consider a stress function φ and a shell z that both have a radial symmetry defined in a polar coordinate system. Then, the vertical equilibrium becomes where r C ≥ ≥ 0 is the parameter in the radial direction in the polar coordinate system. By assuming the same shape for the shell and stress function, so that φ( ) = ( ) r Az r where A is a constant, the solution to equation (36) is given by with r B ″ ″ 1 . Thus, at the circumference where r B = , there are concentrated stresses parallel to the circumference direction. If vertical support is provided at r C = , the edge at r B = can "float" without any support and the circumferential concentrated stress is tensile.
Let A = 1/ 2 , B = 1, C = 0.2 , and ρ = 1. Then, the total height of the shell is given from the analytical solution as z B z C z z ( ) ( ) = (1.0) (0.2) = 0.5867 − − , and the cross section of the shell is plotted in Figure 13. Solving z numerically with the proposed method gives the shape shown in Figure 14 whose profile and total height, 0.58 , match the analytical solution.
The blue surface in Figure 15 is a plot of the left-hand side of equation (36) multiplied by 1/ ρ obtained numerically with a NURBS surface with 32 32 × control points, Dim = 3 , and integration point multiplier NI = 1. Similarly as in the hypar example, the blue surface sits at a level −1 below the red surface representing zero loading, again indicating that the method provides the sought solution.
A convergence analysis was also conducted. Figure 16 shows plots of the energy as a function of the number of control points, proving the solution converges to a single solution.

Reaction forces at supports
Axial force in beams (0,1) (1,0) Figure 10. Force diagram and the "jumps" of the hypar problem. Figure 11. Left-hand side of vertical equilibrium obtained using 32 32 × control points (blue surface) and zero-loading reference (orange surface) for the hypar problem.
Unlike pure-compression or pure-tension zones, there is a restriction on the BCs in tension-compression mixed zones. In this example, because we assumed a radial symmetry when deriving the stress function, fixing the rim of the oculus along a flat circle should work, that is, it is compatible to the problem. This compatible boundary condition may not be the unique one and there may be more, but, in general, an arbitrarily determined boundary condition is incompatible, making the problem ill-posed.
As such, in general, a hyperbolic partial differential equation is well-posed only when a proper set of compatible boundary conditions is given. This makes solving the hyperbolic problems substantially different from solving elliptic problems, that is, pure-compression or pure-tension ones. The next section further elaborates the topic about compatible and incompatible BCs of hyperbolic PDEs.

Asymptotic lines and boundary conditions
Within areas of pure compression or tension, support BCs cause singularities that dissipate rapidly and only affect the solution locally, allowing free placement of the supports in such areas. However, in areas with a mix of compression    and tension, BCs give rise to singularities that are transferred through characteristic lines with no dissipation, often making the hyperbolic PDE ill-posed. 5,64 Such incompatible BCs should be avoided. The four green surfaces in Figure 17(a) are the numerical solutions obtained by imposing two compatible sets of BCs and two incompatible ones on the hypar example introduced in Section 3.1. As in Figures 11 and 15, the blue surfaces of Figure 17(a) are supposed to drop by 1.0 below the reference surface if the PDE is well-posed. However, the incompatible BCs cause disturbances seen as bumpy ridges that run across the solution, and the blue surfaces do not drop by 1.0 but form spikes.
A second-order hyperbolic PDE is often called a wave equation, and waves propagate along so-called characteristic lines. Thus, a perturbance transfers from one end of a characteristic line to the other end with no dissipation, and if no compatible BC exists at the other end, the problem becomes ill-posed (i.e. no solution exists). Even if solutions do not exist, the proposed method returns a numerical solution with no errors. Hence, extra attention must be paid to this issue.
In the problems discussed in this paper, the characteristic lines are the projections of the asymptotic lines of the stress function to the ground plane (see Chapter 4 in Csonka, 5 Chapters 1 and 2 in Sanchez-Palencia et al., 65 and Appendix C). An asymptotic line is a line on a surface whose normal curvature is zero, and it only exists in areas where the Gaussian curvature is negative. Typically, the asymptotic lines are a group of diagonal lines that run across the negative Gaussian curvature area, intersecting the principle curvature lines at roughly 45 degrees.
The general solution of the hypar problem is z x y F x y G x y ( , ) = ( ) ( ) + + − with F and G arbitrary functions that are determined such that boundary conditions are satisfied, and the asymptotic lines are such that they result in characteristic lines given by x y − = const. and x y + = const. , which run parallel with the edges of the shell. As seen in Figure 17(a), the BCs that contradict the general solution are incompatible or, in other words, BCs that do not have compatible BCs at the end of the characteristic lines are incompatible. Figure 17(b) and (c) show more examples of compatible and incompatible BCs imposed on the bowl example of Section 3.2 respectively on a problem reproducing Félix Candela's Xotimilco restaurant (see Appendix D). The BCs for these problems were guessed based on the knowledge of radial symmetry. However, in general, identifying compatible BCs for hyperbolic PDEs is a challenge.
To avoid this challenge, "safe" areas with positive Gaussian curvature may be added to the stress function in which the supports are placed. This strategy is illustrated in the examples in Figure 18 where the purple surfaces are the used Airy stress functions with their asymptotic lines drawn in white. Though not shown, the "blue surfaces" dropped properly for all cases. The safe-area-strategy is used in the examples discussed in Section 5.

More complex example
A stress function that can take a general plan geometry where v is bounded by v min > 0 and v max , and Q v ( ) defines the generatrix profile curve of the stress function. Then, from equation (10), the force diagram follows as where a prime is used to denote differentiation.  If ′ Q v ( ) = 0 at an edge of the shell, the stress function tangentially touches the horizontal plane intersecting the edge of the stress function. Such points on the form diagram x are mapped to the center of the force diagram X . Therefore, the distances between any proximity points of the edge as measured in the force diagram are always zero. Hence, there is no horizontal thrust acting on the edge. In contrast, if ′ ≠ Q v ( ) 0 at an edge, the stress function forms a kink with the horizontal plane intersecting the edge.
As pointed out before, for example, when discussing the kinks between the hypar stress function and the surrounding planes in Figure 10, a kink represents concentrated stresses acting along the edge, which are balanced by stresses acting in the perpendicular direction. The magnitude of the axial force can be obtained by measuring the "jump" from the corresponding points of the edge in the force diagram to its center, which represents the axial force of the reinforcing beam placed along the perimeter of the edge.
The proposed stress function, form diagram, and force diagram of equations (38) and (39) are valid not only for smooth guide curves but also for piecewise smooth guide curves. In those cases, the stress function also becomes a piecewise smooth function, with creases along lines intersecting the discontinuities in the guide curve. As discussed in Section 2.7, such stress functions result in "jumps" in the force diagram representing concentrated axial forces running along the corresponding creases in the shell, and their magnitudes are given by equation (29). Figure 19 shows a form diagram generated using a piecewise smooth guide curve with v min = 0.2 and v max = 1.0 . Vertical supports are provided along three short stretches of the exterior perimeter of the form diagram. Along these stretches, the form diagram is modified with small semi-circles that ensure a local area of pure compression close to the support avoiding problems of incompatible BCs.
Let the stress function φ = ( ) Q v be such that where l is a parameter embedded to the equation to control shape of the guide curve (generatrix). Then, a stress function that comprises three smooth areas and three crease curves is obtained, as shown on the left in Figure 20, resulting in a force diagram shown to the right. As illustrated in Figure 20, if one travels from point a to b to c on the stress function, corresponding travel from A to B to C takes place on the force diagram where point B exists on both sides of the "jump," representing the crease running through point b on the stress function. Figure 21 shows two numerical form-finding results obtained using the discussed functions φ , x , and X , one using the linear method (Section 2.4) and the other the nonlinear method (Section 2.5). Each leaf is represented by 2 NUBRS surfaces with 32 16 × control points, Dim = 3 and NI = 2 , and l of equation (40) was set to 0 .
The resulting shells are rather similar, with the nonlinear solution a little higher than the linear solution. The difference is due to the accurate account of the surface area in the nonlinear method. The shells were obtained with loading ρ = 0.05 − , and the nonlinear method converged, as proved by the history of the residual norm shown in Figure  22. However, when the load was increased to ρ = 0.15 − , the nonlinear method did not converge. As such, the nonlinear method does not converge when the loading is large. However, fortunately, in the cases we tested, although it depends on the problems, the height of the shell at the maximum loading ρ that gave a convergence gives forms that are usable as architectural spaces.

Verification
A convergence study was conducted by increasing the number of control points from 20 10 × to an extreme of 400 200 × . Dim = 3 NURBS surfaces were used as finite elements, and the integration point multiplier was set to NI = 2 . As shown in Figure 23, the solutions were very stable even if the number of control points was increased. Figure 24 is a plot of the energy versus the number of control points, and the solutions converged.
So far, in the form-finding process, only an ideal membrane stress state is considered, and thus, structural properties such as thickness, stiffness, and deformation are not studied.
Using the shell obtained from the nonlinear method, a series of deformation analyses was performed by varying the thickness of the shell from a realistic one to an extremely small one. The three crease curves where modeled as solid beams with cross-section width × height = 0.2 m × 0.3 m, and standard steel was used for both the shell and beams. The analysis was performed using NURBS-based shell and beam elements provided by Kiwi3D, 66 which consider both membrane and bending stiffness. The NURBS surfaces and curves obtained in the form-finding stage was used as finite elements without meshing. Table 1 shows the maximum deflection with varied thicknesses. Even for a small thicknesses, the deformations are kept rather low. This indicates that the shell is working primarily in membrane action, just as intended.
Another way to check if the shell works in pure membrane action is to study the ratio between the membrane energy and the bending energy of the shell. As can be seen in Figure 25, the membrane energy becomes more dominant as the thickness of the shell is decreased.   × control points, Dim = 3 and NI = 2 . As expected, the blue surface drops by 1.0, except at two areas which correspond to the locations of the supports.
From these observations, it can be concluded that the obtained solution is an approximation of the original equilibrium equation and that the shell works by membrane action.
Remarkably, few asymptotic lines of the stress function ( Figure 27, left) emerge as waves running on the formfinding result (Figure 27, right). Compared to the bumpy ridges in the numerical solutions obtained with incompatible BCs discussed in section 4, these waves look clean, and there are no spiky noises in the blue surface except the areas near the supports. Thus, these waves can be distinguished from the numerical errors caused by the incompatible BCs; rather, it can be considered that a solution of a wave equation is obtained. Figure 28 shows an architectural rendering that illustrates an implementation scenario of the example.

Conclusion
The numerical form-finding of membrane action shells containing a mix of tensile and compressive stresses involves solving a hyperbolic second-order partial differential equation. In this study, an existing NURBS-based isogeometric approach originally designed for the formfinding of pure compression or pure tension membrane shells was developed further. It is demonstrated that the developed method can be used to solve tension-compression mixed type hyperbolic equilibrium equations with a satisfactory accuracy. The method takes as input a form diagram and a stress function from which a force diagram is computed, and their reciprocal relations are discussed.
In areas where the stress function describes a state of both tension and compression, it has negative Gaussian curvature and thus asymptotic lines, and it was recommended to put supports away from such areas. Otherwise, the boundary conditions must be compatible with the equilibrium equation to make the problem well-posed.
The proposed method was tested using simple problems whose analytical solutions are known, and the numerical and analytical results are shown to have a good agreement. A solution to a more complex problem was also demonstrated and verified, and the formulation of the problem can easily be modified to suit a wide range of architectural applications.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was supported by Skidmore, Owings, & Merrill from 2015 to 2020.  where g ij is the inverse of g ij . The same relations can be used for a force diagram. The deformation gradient F is given by equation (16) and its components on a ( , )

ORCID iDs
x y coordinate system: where e x and e y the orthonormal base vectors of the ( , )