1 Introduction
Origami has emerged as a technology to construct threedimensional (3D) folded structures by folding a piece of flat material. The folded configurations that can be attained depend on the crease pattern which consists of mountain/valley creases and vertexes. Among the various crease patterns, the standard Miuraori pattern is probably the most wellstudied and widely used pattern in engineering. It was proposed for the packaging and deployment of large membranes in space
miura1970proposition; miura1985method; the pattern and its variants have found applications as the foldcores of sandwich structures miura1972zeta; klett2011designing; heimbs2013foldcore; ma2018origami. Recently, unique properties of origami structures, such as negative Poisson’s ratio and multistability, have been revealed by researchers and the Miuraori tessellation has inspired the design of mechanical and acoustic metamaterials, see wei2013geometric; schenk2013geometry; silverberg2014using; pratapa2018bloch among others.The standard Miuraori tessellation is made up of a number of identical unit cells, or Miura cells, and each cell consists of four congruent parallelograms, see Fig.1(a), (b) and (c). Each interior vertex is of degree4 with one pair of creases symmetric with respect to the other collinear pair. For stiff facets connected by soft creases, it would be much harder to bend the facets than to fold along the creases, thus the morphing of the origami structures can be described by the rigid origami model in which all the facets are assumed rigid and only folding along the creases are allowed. Treated as the rigid origami, the Miuraori tessellation constitutes a mechanism which folds with a single degree of freedom (DOF) before all the facets collapse into the same plane. In all the partially folded states, the Miuraori tessellation can only exhibit inplane motion and tiles the plane. To approximate 3D curved surfaces, variations in the edge length and sector angles of the cells should be involved, see Fig.1(d) for an illustration of a generalized Miura cell. However, for a general quadrilateral mesh, rigidfoldability is nontrivial since the compatibility between all the rigid facets during the folding process can lead to an overconstrained mechanism tachi2009generalization; dieleman2019jigsaw. It has been shown by Tachi that standard Miuraori pattern can be perturbed into general ones which satisfy the developability, flatfoldability and other constraints tachi2009generalization; tachi2010freeform_quad. The resulting generalized patterns can approximate some curved surfaces when partially folded. Most importantly, it has been proven that if an intermediate folded state exists for a quadrilateral mesh consisting of developable and flatfoldable vertexes, the pattern is rigidfoldable limited only by the avoidance of selfpenetrations tachi2009generalization; lang2017twists. Based on the findings, Tachi has developed the “Freeform origami” software which can modify the standard Miuraori tessellation (and other tessellations) into some curved structure by interactively moving the vertexes and projecting the perturbation into the constraint space tachi2010freeform; tachi2010freeform_quad. However, it is difficult for the software to tackle inverse design problems since constraints regarding approximating target surfaces are not considered yet. By altering a single characteristic of the Miuraori cell, Gattas et al. have investigated the parameterizations of firstlevel Miuraori derivative patterns gattas2013miura. Lang and Howell proposed a method to construct rigidfoldable quadrilateral meshes based on the inherent relations between the sector and fold angles of flatfoldable vertexes lang2018rigidly. Recently, Dieleman et al. developed a systematic approach by first identifying a number of rigidfoldable motifs as jigsaw puzzle pieces and then fitting them together to obtain large rigidfoldable crease patterns dieleman2019jigsaw.
Several approaches have been developed for the design of cylindrical and axisymmetric origami structures zhou2015design; wang2016folding; song2017design; hu2019design. Due to the symmetry of the target cylindrical/axisymmetric surfaces, the crease patterns observe collinear creases which are parallelly (radially) spaced for the cylindrical (axisymmetric) cases. Though restrictive in the designed geometry, the approaches in the references zhou2015design; wang2016folding; song2017design; hu2019design can be used to design folded structures which fit two given target surfaces. Dudte et al. have proposed an optimizationbased algorithm to find the intermediate folded form which approximates general curved surfaces dudte2016programming. In this algorithm, the target surface is first tiled with a Miuralike tessellation. To ensure the approximation, the four corner vertexes of each cell (vertexes 1, 3 ,7 and 9 of each cell shown in Fig.1(d)) are fixed as the initial guess whilst the positions of the other vertexes are searched by the optimization method to meet the constraints, including the planarity of all the quadrilaterals and the developability at all the interior vertexes. However, it was reported that the algorithm fails to find the folded form that satisfies the flatfoldability condition for general curved surfaces. Due to the aforementioned theorem by Tachi tachi2009generalization, the patterns constructed by Dudte et al. are not rigidfoldable, i.e., the quadrilateral facets will be bent along the diagonals when it is folded from the flat state to the target state. To approach the flatfoldability (and thus the rigidfoldability), the flatfoldability constraint is replaced by inequalities with auxiliary tolerance dudte2016programming. However, this remedy can only reduce the bending energy of the quadrilateral whilst the crease patterns consist of triangular facets during the folding process. While introducing diagonal creases enlarges the space of rigidfoldable configurations, extra efforts are required when folding the flat pattern into the target form due to the introduced DOFs and the complex interaction between the energy potential and the rigid origami waitukaitis2015origami. The failure in satisfying the flatfoldability condition may be caused by the fact that all the four corner vertexes of each cell are fixed during the optimization which reduces the number of variables in the algorithm. In fact, the position of corner vertexes should also be rearranged to enlarge the searching space such that flat and rigidfoldable crease pattern may be found for a curved target surface.
In this paper, based on the constrained optimization, we present the construction of generalized Miuraori tessellation for 3D curved parametric surfaces which can be folded rigidly from the flat state to the target form with a single degree of freedom. The construction utilizes general quadrilateral mesh with developable and flatfoldable interior vertexes whilst the mountain and valley crease assignment is kept the same as that of the Miuraori pattern. The target surface is first tiled with triangulated Miuralike cells and then the optimization method is used to enforce the constraints, including the planarity of all the quadrilaterals, developability and flatfoldability at all the interior vertexes. To ensure the fitting, a portion of vertexes on the one side is attached to the target surface; another portion of vertexes on the other side is also attached to the second target surface if a second target surface is to be approximated. For the attached vertexes, the parametric coordinates are adopted as the unknown variables whilst the Cartesian coordinates are the unknowns for the other vertexes. Several examples of fitting a single target surface and two target surfaces are considered to demonstrate the feasibility of the algorithm.
2 Approximate a single target surface
2.1 Initial tessellation construction and the constraints
Fig.1(c) shows the standard Miuraori cell made up of four congruent parallelograms. The interior vertex is of valence four with one valley crease and three mountain creases or vice versa. Besides, one pair of creases are symmetric with respect to the other collinear pair. The Miura tessellation is the repeating of the unit cells and tiles the plane in its partially folded states. To approximate curved surfaces, the edge length and sector angles of the cells should vary across the pattern while the basic rules of Miura pattern, including developability, flatfoldability and rigidfoldability, have to be fulfilled. A generalized Miura cell is illustrated in Fig.1(d) which consists of four quadrilaterals and the cell is developable and flatfoldable, i.e., the sector angles around the central vertex satisfy
(1) 
for developability and flatfoldability, respectively tachi2009generalization; lang2017twists.
As pointed out by Dudte et.al dudte2016programming, the basic steps for constructing tessellations for curved surfaces based on optimization are: first, generate an initial Miuralike tessellation as an initial guess for the optimization algorithm which does not have to satisfy the constraints; then, enforce the constraints, including developability, flatfoldability and etc, by the constrained optimization with appropriate objective function. Similar ideas have also been adopted by Bhooshan to generate curved folding based on the dynamic relaxation framework bhooshan2016interactive.
For expository purposes, we consider the target surface with as an example. The target surface can be expressed in the parametric form as . Fig.2 shows the construction of the initial triangulated Miuralike tessellation:

A base mesh is first generated by the and coordinate lines with and for and . In other words, the  and coordinate lines are equispaced in the  and coordinates, respectively.

Then, the vertexes on the coordinate lines with and 8 are moved along the coordinate by , see Fig.2(b).

Finally, the vertexes on the coordinate line with and 8 are moved along the normal of the target surface by ; the quadrilateral mesh is triangulated by connecting vertex() with vertex() for and vertex() with vertex() for .
The vertexes of each quadrilateral are noncoplanar in general and a representative triangulated cell is shown in Fig.3(a). It should be remarked that the initial triangulated Miuralike tessellation is not unique. By changing and (comparable to in general), different initial tessellations can be constructed. Compared with the construction process in the reference dudte2016programming, the current construction suits for a large class of parametric surfaces and avoids the process of merging vertexes shared by adjacent cells.
The triangulated Miura cell in Fig.3
(a) can be specified by the coordinates of the 9 vertexes. The planarity of the quadrilateral facets are to be restored in the optimization process. For the quadrilateral 1254, the planarity of the quadrilateral facet requires that the volume of the hexagon formed by the vectors
, and equals 0, i.e.,(2) 
where is the vector from vertex to vertex and is the Cartesian coordinate for the vertex. The planarity for the other quadrilaterals can be expressed analogously. Due to the triangulation, all the inner vertexes are of degree6, see Fig.2(c). In Fig.3(a), the sector angles are denoted anticlockwise as through around the central vertex5. The developability condition at vertex5 is
(3) 
When the coplanarity of the quadrilaterals 2365 and 5698 are restored, the diagonal lines 35 and 59 becomes redundant and the vertex5 reduces to degree4 vertex. In the optimization process, the flatfoldability condition at vertex5 is
(4) 
It is clear that the identity, , holds if Eq.(3) and Eq.(4) are satisfied. Considering the orientation of the quadrilateral facets and triangulation, the numbering of the sector angles in Fig.3(a) can be used for vertexes on the coordinate lines while the numbering shown in Fig.3(b) should be adopted for the vertexes on the the coordinate lines such that Eq.(4) is consistent for all interior vertexes. To facilitate the optimization algorithm, the derivatives of the constraints in Eq.(2), (3) and (4) with respect to the variables are supplied in the A. The triangulated Miuralike cell will turn into a generalized Miura cell if the constraints on quadrilateral facet planarity, developability and flatfoldability are satisfied simultaneously.
2.2 The variables and objective function
To enforce a close approximation to the target surface, a portion of the vertexes is attached to the target surface. For the cell in Fig.3(a), it is most likely that the vertexes in the set {1,4,7,3,6,9} reside on the one side while vertexes in {2,5,8} on the other. Here, the corner vertexes of each cell, i.e., vertexes {1,7,3,9} in Fig.3(a), are chosen to be attached to the target surface. The attached vertexes are allowed to slide on the target surface. For the attached vertex, it is natural to use their Cartesian coordinates as the unknown variables while enforcing the attachment by satisfying the equation of the target surface, for example in the form of . This approach is still feasible when the target surface has explicit expression regarding , and coordinates. However, for many parametric surfaces, it would be difficult to derive such a form and numerical root finding and derivative calculations are required. Instead, we adopt the parametric coordinates {} as the unknown variables for the attached vertex, i.e., , thus the vertex naturally resides on the target surface. Aside from the attached vertexes, the Cartesian coordinates are adopted as the unknown variables for the other vertexes.
With the triangulated Miuralike tessellation as an initial guess, the optimization algorithm is utilized to adjust the vertex positions such that the constraints are fulfilled while minimizing appropriate objective function. Let be the set of all the vertexes; is set of auxiliary cell contour edges for all the triangulated Miura cells, see the dashed lines in Fig.3; be the set of all the edges and the total number of edges is . Assuming that the initial triangulated Miuralike tessellation is a reasonable guess, the following objective function is suggested
(5) 
where is the edge length of the current configuration whilst is the counterpart of the initial guess; is the characteristic length; is the coordinate for the vertex of the current configuration whilst is the counterpart of the initial guess. Furthermore, the parametric coordinates {} are the unknowns for the attached vertexes, i.e., , whilst the Cartesian coordinates are the variables for the other vertexes. The first term on the left hand of Eq.(5) is to preserve the relative distance between the vertexes such that no crease or boundary edge degenerates. The second term is to ensure that the vertexes are as close as possible to the initial guess, which also prevents the entire tessellation from sliding on the target surface. Both terms are dimensionless such that the objective function is scale independant.
For the tessellation with triangulated Miuralike cells, there are quadrilaterals and each quadrilateral brings in one planarity constraint, see Eq.(2). Besides, there are interior vertexes and, at each inner vertex, Eq.(3) and Eq.(4) should be satisfied for the local constraints on developability and flatfoldability. Thus, the total number of constraints is
(6) 
It is clear that the total number of vertexes is . As the four corner vertexes of each cell are chosen for the attachment, there are vertexes with their parametric coordinates as the unknowns and each parametric vertex owns 2 variables. The total number of unknown variables and the excess DOFs are, respectively,
(7) 
When , i.e., there are more constraints than variables, the optimization algorithm generally fails to find a generalized Miuraori tessellation if there were no sufficient redundancy among the constraints. For the case of , it is required that to ensure that .
The objective function in Eq.(5) together with the constraints in Eq.(2), (3) and (4) can be solved by standard optimization routines. Here, we adopt the interiorpoint method implemented in the fmincon routine of MATLAB. To facilitate the solution, the derivatives of the objective and constraint functions with respective to (w.r.t) the variables should be provided for fmincon and are discussed in the A.
For the construction of generalized Miuraori tessellation with 44 cells for , the fmincon converges in 93 iterations and the constraints are satisfied to the machine accuracy. The converged configuration is shown in Fig.4. It has been proved by Tachi that if there exist a partially folded form of the quadrilateral mesh with interior developable and flatfoldable vertexes, then the quadrilateral mesh is rigidfoldable limited only by the penetration avoidance (see Theorem 2 in the reference tachi2009generalization). Thus, the origami is rigidfoldable with a single DOF. The crease pattern is given in Fig.4(b) in which the vertexes indicated by the blue dot are exactly on the target surface in the folded form shown in Fig.4(a). Since the origami constitutes a single DOF mechanism, the folded state can be characterized by any one of the dihedral angles between two adjacent facets. Here, as indicated in Fig.4(a), the dihedral angle at the red crease is chosen which equals approximately 87.4° for the target state. Following the procedures in Section 4 and Appendix A of the reference lang2018rigidly, the dihedral angles at all the interior creases and thus the 3D folded form can be obtained for given . It is checked that no penetration occurs between the facets, see B for detailed discussions. Several snapshots during the folding process are shown in Fig.5.
2.3 Further design examples
This subsection considers several examples to further illustrate the design algorithm which is implemented in MATLAB. The code is executed on a desktop with Intel(R) Core(TM) i76700 (8 cores, 3.41Gz). The details of the design cases are summarized in Table 1. The most time consuming part of the whole construction is to solve the optimization problem by the interiorpoint method of fmincon. As a rough reference of computational efficiency, the number of iterations and computational time (only the solution time of the fmincon routine) are listed. Besides, it is checked that the constraints including quadrilateral facet planarity, developability and flatfoldability at the interior vertexes are fulfilled to machine accuracy (of order less than ). From the Table 1, it takes more iterations to search the solution for the sphere case. The converged folded forms, the crease pattern and the fully folded state are, respectively, shown in the first, second and third columns of Fig.6. It can be seen that the crease patterns for the cases (a), (c) and (d) are roughly symmetric with respect to the middle horizontal line.
3 Fit two target surfaces
3.1 Construction process
It is of practical interest that the origami, when partially folded, the vertexes on each side approximate different given surfaces. Without losing generality, the two surfaces are denoted as lower surface and upper surface . The construction of the initial triangulated Miuralike tessellation follows similar procedure as that for the single target surface. For clarity, the case with and is considered as an example. Fig.7(a) shows the base mesh generated on the lower surface by the  and coordinate lines. The coordinate lines are equispaced in the parametric coordinate and so are the coordinate lines. The vertexes on the coordinate lines ( and 4) are then moved along the coordinate by to form the mesh shown in Fig.7(b). After that, the vertexes on the coordinate lines with and 4 are moved to the upper target surface by changing the coordinate of each vertex to . The vertexes on the blue lines in Fig.7(c) are then on the upper target surface. The initial construction does not constitute a generalized Miuraori tessellation as the constraints, i.e., quadrilateral facet planarity, developability and flatfoldability, are not fulfilled in general.
For the cell shown in Fig.3(a), the vertexes {1, 3, 7, 9} are attached to the lower surface whilst the vertex5 is attached to the upper surface. For the attached vertexes, the associated parametric coordinates are the variables whilst the other vertexes own the Cartesian coordinates as the variables. The construction then resorts to the optimization algorithm (the fmincon routine in the current implementation) to enforce the constraints on quadrilateral facet planarity, developability and flatfoldability, i.e., Eq.(2), (3) and (4), respectively. For the case in Fig.7, the fmincon convergences in 85 iterations and the final configuration is shown in Fig.7(d) and (e). The crease pattern is shown in Fig.7(f). It is checked that the constraints are fulfilled to machine accuracy. Due to the Theorem 2 of the reference tachi2009generalization, the pattern is rigidfoldable and no penetration is observed before it is fully folded.
For a tessellation with cells, the number of constraints are the same as Eq.(6). Compared to the single target surface case, the central vertex (vertex5 in Fig.) of each cell is also attached to the upper surface, thus the total number of variables and the number of excess DOFs are
(8) 
respectively. For the case with , it is required that to ensure that , i.e., the number of variables should be lager than that of the constraints. When , the optimization algorithm generally fails to find a generalized Miuraori tessellation.
3.2 Further design examples
This subsection considers several examples to further illustrate the design algorithm for approximating two target surfaces. The details of the implementation are the same as those in section 2.3. Table 2 summarizes the parameters for each cases. The resulting folded forms are shown in the first and second columns of Fig.8 whilst the crease patterns and the fully folded states are listed in the third and fourth columns, respectively.
3.3 Airfoil example
In this section, the lower and upper target surfaces are adopted from the NACA2412 airfoil, see Fig.9. The parametric equations are given in the C. It is difficult to find explicit expressions regarding } for the target surfaces. Thus, if the Cartesian coordinates of all the vertexes were chosen as the variables, it would be inconvenient to enforce the attachment of vertexes to the target surfaces through additional constraints regarding the variables }. Using the parametric coordinates for the pertinent vertexes as the variables avoids this difficulty. Besides, the vertexes at the bottom and tip ends remain respectively in the and planes during the optimization which is enforced by adding linear constraints on the pertinent coordinate ( for the bottom end and for the top end) or coordinate for vertexes with parametric variables ( for the bottom end and for the top end).
First, a tessellation with 312 cells is considered. Similar to the previous examples, the vertexes {1,3,7,9} and {5} of each cell are, respectively, attached to the lower and upper target surfaces. The resulting folded form is shown in Fig.10(a) and (b) with different view angles. Fig.10(c) shows the crease pattern in which the blue dots and purple squares are on the lower and upper surfaces when the pattern is folded to the target state. From Fig.10(a) and (b), it can be seen that a large number of vertexes (51 out of the 175) extrude out of the target surface. By observing the relative position of the vertexes on each zigzag line in Fig.10(b) and (c), the vertexes attached to the lower and upper target surfaces are reselected which is shown by the blue dots and purple squares in Fig.10(f). From the same initial triangulated Miuraori tessellation, the optimization algorithm yields the improved folded form, see Fig.10(d) and (e), and few vertexes (10 out of the 175) are outside of the region. Next, the tessellation with 416 cells is considered. When the vertexes {1,3,7,9} and {5} of each cell are, respectively, attached to the lower and upper target surfaces, the algorithm fails to find a valid folded form. This failure may indicate that too many vertexes are attached to the target surfaces. To enlarge the solution space, fewer interior vertexes are attached to the target surfaces, see Fig.11(c). The folded form is shown in Fig.11(a) and (b).
4 Conclusion
This paper has presented the construction of rigidfoldable generalized Miuraori tessellations which can approximate curved parametric surfaces based on constrained optimization. The algorithm is capable to approximate a single or double target surfaces by attaching portions of vertexes to the pertinent target surfaces. All the inherent properties of standard Miuraori fold, including quadrilateral facet planarity, developability and flatfoldability, are enforced through optimization. To enlarge the searching space, the attached vertexes are allowed to slide on the target surface through using their parametric coordinates as the variables. Tessellations for target surfaces with varying curvatures in both directions have been constructed which are rigidfoldable with a single degree of freedom.
A limitation of the current algorithm is that the number of attached vertexes is restricted such that the optimization problem is not overconstrained, see Eq.(7) and (8) for the approximation of single and double target surfaces, respectively. While a reselecting strategy is suggested and illustrated in Section 3.3, vertexes to be attached to the target surfaces are still chosen empirically in both layout and number. Neverthless, the present algorithm can be an efficient tool to generate small, rigidfoldable, curved patches which may be used as building blocks to piece up larger rigidfolable crease patterns lang2017twists; lang2018rigidly; dieleman2019jigsaw. It is noted that the level of approximation is not exactly controlled in the present algorithm which may be incorporated into the objective function in future works.
Acknowledgment
The authors would like to thank Dr. Shuo Feng for insightful discussions. This work is supported by the National Natural Science Foundation of China (Nos. 11272303, 11072230, and 11802107), the Fundamental Research Funds for the Central Universities of China (Nos. WK2480000001 and WK2090050040) and the support of Entrepreneurship and Innovation Doctor Program in Jiangsu Province.
Appendix A Derivatives of the objective and constraint functions
This appendix discusses the derivatives of the objective and constraint functions which should be provided to enhance the efficiency of the optimization solver. The derivatives with respective to (w.r.t) the Cartesian coordinates are considered first. The derivatives of the objective function in Eq.(5) and the quadrilateral facet planarity constraint in Eq.(2) w.r.t to the Cartesian coordinates can be ready obtained and are not detailed here. For an interior vertex, the developability and flatfoldability constraints are the sum of the pertinent sector angles around the vertex. Each sector angle is determined by the three vertexes of the triangle where it resides. For a typical triangle shown in Fig.12, the sector angle is given by
(9) 
The derivatives of w.r.t the vertex coordinates are (also see Eqs.(12), (13) and (14) of reference tachi2010freeform)
(10a) 
(10b) 
(10c) 
where is the unit normal of the triangle. Referring to the unit cell in Fig.3
, the derivatives of the developability and flatfoldability constraints can then be obtained by adding the contributions from the sector angles around the central vertex. The derivatives w.r.t the parametric coordinates are obtained by the chain rule. For instance, assume that the vertex1 of the triangle in Fig.
12 is attached to the target surface and its unknown variables are , then the derivatives of the sector angle w.r.t the parametric coordinates are(11) 
Appendix B Check of selfintersection
Thanks to the Theorem 2 of the reference tachi2009generalization, the crease pattern is rigidfoldable if an intermediate folded form exists for the quadrilateral mesh consisting of developable and flatfoldable vertexes. Nevertheless, the physically valid folded form of the crease pattern is limited by selfintersection. For the generalized Miuraori pattern considered here, it is obvious that no selfintersection will be observed when the crease pattern is folded only by a small amount. As the folded form is of single degreeoffreedom (specified by the characteristic fold angle ) and the shape morphing is smooth, contacts between the facets will happen when reaches a certain critical value. After the converged design has been obtained for each case, we can view the threedimensional folded form for given and visually check whether global penetration occurs. For instance, for the case of in Section 2.1 and 2.2, the folded form with is plotted in Fig.13. There is no penetration even for this almost fully folded state.
Appendix C Lower and upper surfaces of the airfoil
Fig.9(a) shows the cross section of the NACA fourdigit series airfoil in which is the airfoil chord. Introducing the parametric coordinate , the camber line and the thickness distribution are, respectively, moran2003introduction
(12a) 
(12b) 
in which , and are parameters. For the NACA2412 airfoil considered in the example, the parameters are and . The lower and upper curves of the cross section are respectively
(13) 
and
(14) 
where
The cross sections at the bottom () and tip () ends are prescribed with the chord and , respectively, see Fig.9(b). The lower curves of the cross sections at the bottom and tip ends are respectively
(15) 
The lower surface of the airfoil is a ruled surface which is given by
(16) 
The upper surface of the airfoil is defined in a similar way, i.e.,
(17) 
where
(18) 
are the upper crosssectional curves at the bottom and top ends, respectively.
Comments
There are no comments yet.