Algebraic 3D Graphic Statics: Constrained Areas
AAlgebraic 3D Graphic Statics: Constrained Areas
Masoud Akbarzadeh a and M´arton Hablicsek a, ba Polyhedral Structures Laboratory, School of Design, University of Pennsylvania, Philadelphia, USA b Department of Mathematics, Leiden University, Leiden, NetherlandsJuly 31, 2020
Abstract
This research is a continuation of the Algebraic 3D Graphic Statics Methods. It provides algorithms and (numer-ical) methods to geometrically control the magnitude of the internal and external forces in the reciprocal diagrams of3D/Polyhedral Graphic statics. 3D graphic statics (3DGS) is a recently rediscovered method of structural form-findingbased on a 150-year old proposition by Rankine and Maxwell in Philosophical Magazine. In 3DGS, the form of the structureand its equilibrium of forces is represented by two polyhedral diagrams that are geometrically and topologically related.The areas of the faces of the force diagram represent the magnitude of the internal and external forces in the system.For the first time, the methods of this research allow the user to control and constrain the areas and edge lengths of thefaces of general polyhedrons that can be convex, self-intersecting, or concave. As a result, a designer can explicitly controlthe force magnitudes in the force diagram and explore the equilibrium of a variety of compression and tension-combinedfunicular structural forms. In this method, a quadratic formulation is used to compute the area of a single face based onits edge lengths. The approach is applied to manipulating the face geometry with a predefined area and the edge lengths.Subsequently, the geometry of the polyhedron is updated with newly changed faces. This approach is a multi-step algorithmwhere each step includes computing the geometry of a single face and updating the polyhedral geometry. One of the uniqueresults of this framework is the construction of the zero-area, self-intersecting faces, where the sum of the signed areas of aself-intersecting face is zero, representing a member with zero force in the form diagram. The methodology of this researchcan clarify the equilibrium of some systems that could not be previously justified using reciprocal polyhedral diagrams.Therefore, it generalizes the principle of the equilibrium of polyhedral frames and opens a completely new horizon in thedesign of highly-sophisticated funicular polyhedral structures that are beyond compression-only systems.
Keywords:
Algebraic three-dimensional graphic statics, polyhedral reciprocal diagrams, geometric degrees of freedom,static degrees of indeterminacies, tension and compression combined polyhedra, constraint manipulation of polyhedral dia-grams.
Recently, the geometry-based structural design methods,known as Graphic Statics, has been extended to 3D dimen-sions based on a historic proposition by Rankine and Maxwellin Philosophical Magazine [23, 4, 1, 8, 16, 15, 11, 13, 16, 17,18].In this method which is called
3D Graphical Statics usingReciprocal Polyhedral Diagrams , the equilibrium of the forcesin a single node is represented by a closed polyhedron or apolyhedral cell with planar faces (Figure 1a). Each face ofthe force polyhedron is perpendicular to an edge in the formdiagram, and the magnitude of the force in the correspond-ing edge is equal to the area of the face in the force polyhe-dron. The sum of all area-weighted normals of the cell mustequal zero that can be proved using the divergence theorem[24, 6, 7, 5, 1, 16]. In some cases, a cell can have complex c)b)a) n k ̭ f k f k f i n i ̭ f i A = | f k | (c)(b)(a) n k f k f k f i n i f i A = | f k | Figure 1: (a) A 3D structural joint with an applied force andinternal forces in its members; (b) the form diagram/bar-noderepresentation of the same joint in the context of 3DGS; and(c) the force diagram/polyhedron representing the equilib-rium of the same node in 3DGS.1 a r X i v : . [ c s . C G ] J u l i,j,2 n i,j,1 f i,j f i,j,1 f i,j,1 (a) (b) (c) n i,j n i,j,1 f j,i = A i,j,1 ∙ n i,j,1 + A i,j,2 ∙ n i,j,2 > 0 f j,i = A i,j ∙ n i,j > 0 f j,i = A i,j,1 n ∙ i,j,1 + A i,j,2 ∙ n i,j,2 < 0 ++ - n i,j,2 i,j,2 f i,j,2 f Figure 2: From left to right: (a) a convex face with a positiveforce direction (out of the page); (b) a complex face with twoenclosed regions and a positive net force direction; (c) anda complex face with two enclosed regions and a negative netforce direction. faces (self-intersecting) and which has multiple enclosed re-gions (Figure 2a). The direction and the magnitude of theforce corresponding to a complex face can be determined bysumming the area-weighted normals of all of the enclosed re-gions (Fig. 2b, c). As a result the direction of the internalforce in the members of the structure might flip based on thedirection of the face of a single force cell.Exploiting the potentials of 3DGS in design and engineer-ing requires the ability to manipulate the geometric diagramswithout breaking the reciprocity between the two and in-stantly observe the effect of the change in the other diagram.The existing computational tools for the design and manipu-lation of the reciprocal polyhedrons of 3DGS are quite limited.Moreover, controlling and optimizing the magnitude of forcesby changing the areas of the faces needs efficient algorithms.
In 2016, [1] showed that the reciprocal diagrams of 3DGS canbe constructed in a procedural (step-by-step) approach in aparametric software by assigning constraints between recip-rocal components of each diagram that allows simultaneouscontrol over the geometry of both diagrams. This method isextremely time-consuming and tedious for structures with alarge number of nodes and members [5, 1]. [6] developed acomputational algorithm that could receive convex polyhedralcells as a primal and construct its reciprocal diagram itera-tively within a certain tolerance defined by the user [6, 1].This method cannot deal with (non-convex) self-intersectingpolyhedrons or explore tension and compression equilibrium.Moreover, controlling the areas of the faces was computation-ally quite expensive. In 2018, [13] proposed a method called
Disjointed Force Polyhedra where the equilibrium of the sys-tem was computed by constructing a single convex polyhe-dron for each node using Extended Gaussian Image algorithm[14, 10, 19] and matching the areas of the shared faces [12, 13]. This method allows controlling the areas of the convex cells,but it breaks the reciprocity between the two diagrams. More-over, it cannot control the areas of the self-intersecting faces.Recently, [9] developed an algebraic formulation relating thegeometry of the reciprocal polyhedral diagrams using a lin-ear system of equations. This method can directly constructthe dual from a given primal in one step [3, 9, 2]. Althoughthe previous formulation could immediately construct the re-ciprocal polyhedral diagrams, it did not provide any insighton how to control the areas of the faces corresponding tothe magnitude of the forces in the form diagram. Moreover,the geometrically constrained constructions were also not ad-dressed.
This paper provides a robust algebraic method to constructpolyhedrons with assigned areas and edge lengths of theirfaces from which its reciprocal dual can be constructed as thestructural form. The formulation introduced in this paperrelates the areas of the faces of the polyhedral system to itsedge lengths allowing the combination of this method withthe previous algebraic formulation to control the areas of thefaces. Moreover, the methods of this research can computethe areas of self-intersecting faces with constraint edges whichhas never been addressed in the literate previously. Specifi-cally, this approach can construct zero-area, self-intersectingfaces in the system, where the sum of the signed areas of aself-intersecting face is zero. The existence of such faces in theforce diagram allows to either remove forces in the boundaryconditions or internal forces and therefore, describe internalforce equilibrium that previously was not possible using re-ciprocal polyhedral diagrams.The paper is organized as follows. A quadratic formulationto compute the area is introduced for a single face based on itsedge lengths (Section 2.3). Then, a methodology is describedto manipulate the geometry of the face with a predefined areaand edge lengths (Section 2.4). Subsequently, the geometryof the polyhedron is updated with the newly changed faces(Section 2.5). This approach is a multi-step algorithm whereeach step includes computing the geometry of a single faceand updating the polyhedral geometry. In the end, the dualstructural form is updated with the new magnitude of the in-ternal or external forces (Section 2.6). Section 4 shows theapplication of this method in the design of funicular struc-tures with zero force members or reactions in the boundaryconditions.
We denote the algebra objects of this paper as follows; ma-trices are denoted by bold capital letters (e.g. A ); vectorsare denoted by lowercase, bold letters (e.g., v ), except theuser input vectors which are represented by the Greek letters2 ν and ξ ). Table 1 encompasses all the notation used in thepaper. The first step in the methodology is to link the areas ofthe faces of the polyhedral system to their edge lengths.The definition of the area based on the edge lengths willresult in a quadratic function per face of the polyhedron.As a result, controlling the areas of the faces of the poly-hedral system requires solving a complex system of non-homogeneous quadratic equations simultaneously. This com-plex system of quadratic equations is usually solved us-ing optimization methods. However, to the knowledgeof the authors, these methods usually fail in computingself-intersecting/compression- and-tension combined systems.Moreover, the objective of our research is to provide amethodology to control the areas of the faces without per-turbing the system drastically. Therefore, in this section, weprovide a simple methodology to solve these quadratic equa-tions sequentially by preserving the main geometric featuresof both the form and force diagrams.In this method, there are two types of equations; (i) thequadratic equations that compute the areas of the faces basedon the edge lengths; and (ii) the linear equations that providethe geometry of the faces of the polyhedrons with user-definededge lengths as constraints. The quadratic equations of thefaces are solved using the linear equations around the edgesof each face with constrained edge lengths. Each quadraticequation for a face area has as many variables as the numberof edges of the face which results in a variety of significantlydifferent solutions. However, we can control the solution spaceby reducing the number of variables to one. This allows usto find a solution for the quadratic equation with a limitedgeometric perturbation in the system.In the following sections, we introduce the steps to developthe quadratic equation to compute the area of a face of apolyhedral system based on the edge lengths, and then wedevelop the non-homogeneous linear equation system describ-ing the equilibrium equation for the face with predefined edgelengths. In the end, we show how to solve the equation systemand how to recompute the geometry of the form.
In the previous paper [9], we showed how to write the equi-librium equations for a system of polyhedral cells with planarfaces. For each face f i , we can write an equation based on itsedge lengths that shows the closeness of the face. By choosinga normal vector for each face f i , we can obtain a consistentorientation of the edges. We denote the unit direction vector Topology Description
Γ primal diagramΓ † dual, reciprocal diagram v e f c v † † e † † f † † Matrices M f Area matrix of the face f E f Equilibrium matrix of the face f L f Matrix of predefined edge lengths of the face f E p Equilibrium matrix of the polyhedral system p E † Equilibrium matrix of the dual L p Matrix of predefined edge lengths of the polyhedron p B f Constraint matrix of the face f B p Constraint matrix of the polyhedron p B + p Moore-Penrose inverse of B p RREF f RREF of ( textbfB f | b f ) E † equilibrium matrix of the dual diagram( E † ) + Moore-Penrose inverse of E † Vectors n consistent unit normal vector q vector of edge lengths u j direction vector of edge vector e j b f constraint vector for the face f l f vector of predefined edge lengths of the face f b p constraint vector for the polyhedron p l p vector of predefined edge lengths of the polyhedron p q nci vector of nci edge lengths of a face f i q fix vector of fixed edge lengths of a face f i q nfd vector of nfd edge lengths of a face f i q ci vector corresponding to the edge length of the ci edge e † j edge vector of e † j in Γ † u † j direction vector of edge vector e † j q † vector of edge lengths of Γ † Parameters ν parameter for the MPI method to solve Eq. 30 ξ parameter for the MPI method to solve Eq. 33 Other O centroid of a face h i,j (signed) distance of v j from e i H i average (signed) distance of the v j from e i A f area of face fµ i,j the scalar ratio between e i × e j and n η i,j the scalar ratio between u i × u j and n GDoF f Geometric Degrees of Freedom of face fCGDoF f Constrained Geometric Degrees of Freedom of face fe fix (user)-selected fixed edges of a face e ind list of independent edges of a face e nfd list of nfd edges of a face e nci list of nci edges of a face e pfix (user-)selected fixed edges of the polyhedron q ci length of the ci edge r rank of RREF f Table 1: Nomenclature for the symbols used in this paper andtheir corresponding descriptions.3 j corresponding to the edge vector e j of Figure 3a. Sinceeach face provides a closed loop of edges, the sum of the edgevectors has to be the zero vector. Thus, we obtain a vectorequation for the edge lengths q j of e j as (cid:88) e j u j q j = (1)where the sum runs over the edges e j of the face f i .Thus, each face f i of the polyhedron p provides three equa-tions for the edge lengths, which can be described by a [3 × e ]matrix, E f i E f i q = (2)where q denotes the vector of the edge lengths of the polyhe-dron. This equation describes the geometry of the face.Similarly, we can obtain a [3 f × e ] matrix, E p describing thegeometry of the entire network. Here f denotes the numberof faces and e denotes the number of edges in the network. Inother words, we have a linear equation system E p q = (3)where q denotes (again) the vector of edge lengths of the poly-hedron p . Each solution of the linear Eq. system 3 representsa network, whose edges are parallel to their associated edgesof the original network with different edge lengths. In this section, we explain, how to develop a quadratic systemof equation for a face f i of a polyhedral network based on theedge vectors of the face after [21].Consider a face f i with k vertices: v , v , ..., v k − . Wedenote the edge e i by the vertices v i and v i +1 . Let n be achosen unit normal vector of the face f i . Using the right-hand rule, the normal n provides the direction of the edges.We denote the directional edges by e , e , ..., e k − vectors.For the sake of simplicity, in the following explanation, weuse a cyclic order of the edges, meaning e k and e k +1 will alsodenote the edges e , and e , etc.We can compute the area of the face f i by : • dividing the face f i into k triangles given by the e i andthe geometric center of vertices, O ; and • compute the area of each triangle by computing its height H i which is the distance from O to the edge e i (Figure3b). Average height
The height H i , the perpendicular distance of O from the edge e i , is the average of the signed projected distances h i,j of thevertices v , v , ..., v k − from e i (Figure 4d). For instance, h , , h , , and h , are the signed distances ofthe vertices v , v and v from the edge e . The h , and h , are zero since the vertices v and v lay on edge e . Therefore,the area A i of the triangle O , v i and v i +1 can be written as A i = 12 | e i | H i = 12 k | e i | k − (cid:88) j =0 h i,j , and the total area of the face equals A f = k − (cid:88) i =0 A i = 12 k (cid:88) ≤ i,j ≤ k − | e i | h i,j . (a) v v v v v e e e e e o n f (a) v v v v v e e e e e o n f H H H H H (a) v v v v v e e e e e o n f (b) v v v v v e e e e e o n f H H H H H Figure 3: (a) Vertices and edge vectors of the face f with anormal direction n f ; and (b) dividing the face into triangleswith a base e i and the height H i from the centroid of thevertices.We are looking for a formula for the area A i based on theedge vectors. Thus, let us compute the h i,j based on theedge lengths of the face f i . Recall, that the vertices v i and v i +1 are on the edge e i , hence, h i,i and h i,i +1 are zero (seeFigure 4a, Figure 4d). The first vertex that can contributenon-trivially is v i +2 , and the height, h i,i +2 , can be computedby constructing the triangle of the vertices v i , v i +1 and v i +2 (Figure 4a). We denote the signed area of this triangle by A i,i +2 that can be computed using the two following methods: • the area can be found by the height h i,i +2 (Figure 4 a): A i,i +2 = 12 | e i | h i,i +2 . (4) • also, the cross product of e i and e i +1 provides the signedarea: A i,i +2 n = 12 ( e i × e i +1 ) (5)Note that the sign of the area in Eq. 4 is defined by the h i,i +2 where it can only be negative in a concave or aself-intersecting polygon.From Eqs. 4 and 5, we get the following12 ( e i × e i +1 ) = 12 | e i | h i,i +2 n . (6)4 a) v v v e e o n f (a) v v v v e e e o n f (a) v v v v v e e e e e o n f (a) v v v v v e o n f (a) v v v v v e o n f (a) v v v v v e o n f H H H h h h h h h e + e + e e + e h h h h h h (a) v v v e e o n f (b) v v v v e e e o n f (c) v v v v v e e e e e o n f (d) v v v v v e o n f (e) v v v v v e o n f (f) v v v v v e o n f H H H h h h h h h e + e + e e + e h h h h h h Figure 4: Finding the average height H i for each edge e i byconstructing triangles starting from (a) e i and e i +1 , and (b) e i and e i +1 + e i +2 until we find all heights of the vertices in(c) and (d); and repeating the same process for other edgesto find all H i s in (e) and (f).We have two vectors on the sides of the above equation. Onone hand, it does not make sense to divide vectors by vectors.On the other hand, the vectors e i and e i +1 are perpendicularto n . Hence, the cross product of vectors e i and e i +1 is eitherparallel or opposite to n . Therefore, there exists a scalar µ i,i +1 so that µ i,i +1 n = e i × e i +1 . Thus, we can think of µ i,i +1 as µ i,i +1 = e i × e i +1 n . (7)Going back to Eq. 6, using the notations introduced above,we obtain a formula for h i,i +2 : h i,i +2 = e i × e i +1 | e i | n = µ i,i +1 | e i | . (8) The Eq. 8 is not of the form of a linear equation, becausethe scalar µ i,i +1 depends on the lengths of the edges e i and e i +1 . By dividing both sides of the Eq. 7 by the edge lengthsof e i and e i +1 , the scalar µ i,i +1 will become η i,i +1 = 1 | e i | · | e i +1 | · µ i,i +1 . (9)The scalar η i,i +1 will only depend on the directions of e i and e i +1 and not on the edge lengths. Using this notation, we canwrite h i,i +2 as h i,i +2 = η i,i +1 | e i +1 | . (10)In Eq. 10, the height is expressed as a product of a scalar η i,i +1 that depends only on the direction of the edges, and theedge length of e i +1 . Therefore, we obtained the height h i,i +2 based on the linear function of the edge length.In the next step, we will compute h i,i +3 . Consider thetriangle with vertices v i , v i +1 , v i +3 (see Figure 4b). Notethat the vector from v i +1 to v i +3 is the vector e i +1 + e i +2 . We compute the area A i,i +3 of the triangle using two meth-ods. The area can be computed by using the height h i,i +3 (Figure 4b): A i,i +3 = 12 | e i | h i,i +3 . Also, the cross product of e i and e i +1 + e i +2 provides thesigned area: A i,i +3 n = 12 ( e i × ( e i +1 + e i +2 )) . As a consequence, we obtain the following relationship.12 ( e i × ( e i +1 + e i +2 )) = 12 | e i | h i,i +3 n We combine this equation with Eq. 6 to obtain e i × e i +2 + | e i | h i,i +2 n = | e i | h i,i +3 n (11)Again, we would like to divide Eq. 11 by the vector n .Similarly, the vectors e i and e i +2 are perpendicular to n ,hence e i × e i +2 is either parallel or opposite to the directionof n . We, again, introduce the scalar µ i,i +2 satisfying µ i,i +2 n = e i × e i +2 . With this notation, Eq. 11 becomes µ i,i +2 + | e i | h i,i +2 = | e i | h i,i +3 , so we obtain a recursive formula h i,i +3 = h i,i +2 + 1 | e i | µ i,i +2 . (12)5he scalar µ i,i +2 depends on the lengths of the edge vectors e i and e i + , however η i,i +2 defined below does not: η i,i +2 = 1 | e i | · | e i +2 | µ i,i +2 . The scalar η i,i +2 only depends on the directions of the edges.Using this notation and Eq. 10, Eq. 12 becomes a linearexpression for h i,i +3 h i,i +3 = h i,i +2 + η i,i +2 | e i +2 | = η i,i +1 | e i +1 | + η i,i +2 | e i +2 | (13)Repeating this process for the rest of the edges (Figure 4c)results in a formula for h i,i + l h i,i + l = (cid:88) ≤ m ≤ l − η i,i + m | e i + m | . Finally, we compute the average of these heights, H i byrepeating the same process for all the edges of the face (Figure4e,f): H i = 1 k k − (cid:88) j =2 h i,i + j = 1 k k − (cid:88) j =2 (cid:88) ≤ m ≤ j − η i,i + m | e i + m | = (14)= 1 k k − (cid:88) j =1 ( k − j − η i,i + j | e i + j | (15)As a consequence, we obtain a quadratic formula for thearea of the face f in the edge lengths of the edges of fA f = 12 k − (cid:88) i =0 | e i | H i = 12 k k − (cid:88) i =0 k − (cid:88) j =0 ( k − j − η i,i + j | e i || e i + j | . (16) Quadratic form
The next step is to turn the quadratic equation 16 into aquadratic form with a matrix. Note, that we can computethe coefficients ( k − j − η i,i + j in Eq. 16 without knowing the edge lengths. As a result, wecan formulate the right-hand side of the Eq. 16 in a quadraticform given by a matrix M f , whose entries are given by thecoefficients: M f,ij = ( k − j + i − η i,j if j > i ( i − j − η i,j if j < i i = j . (17)Thus, we can rewrite Eq. 16 as a quadratic form A f = 12 k q T M f q (18)where q is the column vector of the edge lengths q = ( | e | , | e | , ..., | e k | ) T . Symmetric matrix
Usually, the matrix M f is not a symmetric matrix. How-ever, the computations may become simpler if the matrix issymmetric. Indeed, the matrix M f can be turned into a sym-metric matrix. Since (cid:0) q T M f q (cid:1) T = q T M Tf q the quadratic form given as12 k q T M Tf q also computes the area of the face.As a consequence, A f = 12 k q T M f + M Tf q , and in this case, the corresponding matrix M f + M Tf M f appearing in Eq. 16 is alwayssymmetric. In this section, we develop a method to reconstruct a givenface f i by constraining particular user-defined edge lengthsand the target area for the face. To give a general overviewof our approach, consider the face f i of Figure 3a. Initially,without any constraint, we have five unknowns which are theedge lengths of the five edges e − . There are three equilib-rium equations based on Eq.2, in − x , − y , − z around the face f i , and one of them is redundant [3]. As a consequence, thedimension of the possible solutions, i.e. the possible faces,satisfying the equilibrium equations is e − . Instead of solving the quadratic area equation 18 for threeunknowns, we constrain two of them and solve the quadraticequation for only one unknown. Using this technique, wecan significantly reduce the complexity of finding a solutionfor this quadratic equation. This provides additional designpossibilities for the user, as we either allow the user to defineup to two edge lengths out of three or we use the existing edgelengths for two unknown edges and compute the area basedon the last unknown edge.In general, our goal is to simplify solving the quadraticequation of the face by solving it for only one edge length.6 omputing
GDoF f using RREF The dimension of the possible solutions for the geometry ofthe face is called the Geometric Degrees of Freedom (
GDoF f ).In fact, GDoF f describes the dimension of the family of poly-gons with edges parallel to the edges of the initial face whichis always equal to: e − . The GDoF is also equal to the number of independent edgesin each face. In fact, the lengths of the independent edges candefine the lengths of the rest of the edges and the geometryof the face [2]. The independent edges can be found using theReduced Row Echelon form method (RREF).Specifying the edge lengths for the e − e − Constrained Geometric Degrees of Freedom(
CGDoF f ) The user-defined edge lengths provide linear equations for theedge lengths that are in general non-homogeneous. As a con-sequence, the dimension of the solution space for possible ge-ometries of the face may decrease significantly.The user may over-determine the system, for instance,by assigning too many edge lengths. To avoid this prob-lem, we compute the dimension of the constrained solutionspace, called the Constrained Geometric Degrees of Freedom(CGDoF), using RREF.The result of this method classifies the edges into the fol-lowing classes: • the fixed edges, e fix : the edges chosen by the user withpredefined edge lengths (these edges are always depen-dent edges of the equation system); • the non-fixed dependent edges, e nfd : the dependentedges which are not predefined by the user, and • the independent edges, e ind .To solve the quadratic equation for the area, we reduce thenumber of independent edges e ind to one, by assigning the ex-isting edge length for all independent edges except one. Thelast remaining independent edge is called the critical indepen-dent edge, e ci , the length of which we find using the quadraticequation. This method will be described in detail in the nextsections. Defining the Constrained equations for a face
In Section 2.3, we expressed the area of a face polygon basedon a quadratic form of the edge lengths. Now, we developlinear, non-homogeneous constrained equations describing thegeometry of the face with preassigned lengths for certain edgesof the face.We can write the edge e i with a predefined length q i as aconstraint vector equation in the following way: l Ti q = q i (20)where l i is the [ e ×
1] column vector whose entries are all zero(0) except at the index of e i where it is one (1).Similarly, multiple constraints, i.e other fixed edge lengths,can be written as a matrix equation L f q = l f (21)where the rows of L f are the row vectors l Ti and l f is thevector whose entries are the q i , the predefined edge lengths.Together with the equilibrium equations of 2, we obtain allthe linear equations describing the linear constraints whichresults in the constraint equation system B f q = b f (22)where the matrix B f is obtained by stacking the matrices E f and L f B f = (cid:18) E f L f (cid:19) , and the vector b f is obtained as stacking the zero vectorand the vector l f together b f = (cid:18) f (cid:19) We call the matrix B f the constraint matrix and the vector b f the constraint vector. Analyzing the constraint equation system (RREF)
The constraint equation system, Eq. 22, is, in general, a non-homogeneous, linear equation system. The solution space ofthis equation system is often not a linear subspace but ratherthe empty set or an affine subspace of the possible solutionspace R e . Here, e denotes the number of edges. The di-mension of this affine subspace is the Constrained GeometricDegrees of Freedom of the face ( CGDoF f ). The CGDoF f is the geometric degrees of the face after applying the edgeconstraints by the user.It is also possible to have no solutions for Eq. 22. In thiscase, we say that the CGDoF f is −∞ . If there exists a solu-tion to Eq. 22, then the CGDoF f is a non-negative integer,which is the dimension of the affine subspace formed by thesolutions. When the CGDoF f is zero, the constraint equa-tions have a unique solution. If the CGDoF f is positive, then7here are many significantly different solutions to the con-straint equations.In order to compute the CGDoF f , we use the reduced rowechelon form ( RREF f ) of the matrix obtained from the con-straint matrix, B f , and the constraint vector b f : (cid:0) B f b f (cid:1) . The
CGDoF f can be easily computed from this reduced- row-echelon form, but the following two possibilities might occur: • If there exists a row of
RREF f , so that the last entry isone, but all other entries are zero: (cid:0) ... (cid:1) then, the constraint equation, Eq. 22 have no commonsolution. In this case, the CGDoF f is −∞ , and the userneeds to modify the constraints and/or release some ofthe constrained edges from his/her input. • Otherwise, we have at least one solution. In this case,the
CGDoF f equals e − r where e is the number of edgesof the face and the number of columns of the constraintmatrix B f , and r is the rank of the RREF f . This rankequals the number of pivots that also equals the rank of B f . Solving the area equation
The main idea of manipulating the edge lengths of the face toobtain a required area is to solve the area equation Eq. 18 byreducing the number of unknown edges into a single unknownedge length.We reduce the unknowns by finding all the independentedges of the constraint equation system 22 and assigning ei-ther the current values or a user-defined values to them.From now on, we assume that there exists a solution to Eq.22. The columns in
RREF f corresponding to the pivots arecalled the pivotal columns. The non-pivotal columns corre-spond to the so-called independent edges whose lengths canbe manipulated freely. Once the values for the independentedges are set (possibly by the user), there is a unique solutionto Eq. 22.The edges corresponding to the pivotal columns are theedges depending on the independent ones, these edge lengthswill be updated so that the Eq. 22 is satisfied.Our method for solving the area Eq. 18 is to set as manyvalues of the lengths of the independent edges as possible. Inthis case, it is one less than the number of independent edges: e − r −
1. The length of the last independent edge, q ci , is thelength of the critical independent edge, e ci . This length willbe treated as a variable for which we will solve the equation.To solve Eq. 18, we organize the edges according to theform of RREF f into vectors: e fix = e e ci = e e nci = e e nfd = e e nfd = e v v v v v GDoF f = 5 - 2 = 3CGDoF f = 3 - 1 = 2 e fix = e e ci = e e nci = e e nfd = e , e A f = 1357.06 (e.g.) n f n f e fix = e e ci = e e nci = e e nfd = e e nfd = e v v v v v GDoF f = 5 - 2 = 3CGDoF f = 3 - 1 = 2 e fix = e e ci = e e nci = e e nfd = e , e A f = 1357.06 (e.g.) n f Figure 5: A sample face showing the edge vectors, its normaldirection and the choices of e fix , e cr , e nci , and e nfd . • the vector corresponding to the critical edge is defined asan [ e ×
1] column vector q ci : q ci,i = (cid:40) q ci if i is the index of the e ci q ci , the edge length of the e ci edge, is the unknownand we will solve the area equation for q ci ; • the vector of the edge lengths of the non-critical inde-pendent edges, e nci , is defined as q nci which is an [ e × q nci = (cid:40) q i if i is the index of an e nci edge0 otherwisewhere q i are the current edge lengths of the e nci edges; • Similarly, the vector of the edge lengths of thefixed/predefined edges, e fix is defined by q fix as q fix = (cid:40) q i if i is the index of a user-selected edge0 otherwiseThe indices of fixed edges are indices of some of the piv-otal columns. These q i are fixed in the beginning of theproblem and won’t be updated; • Finally, the vector of edge lengths of the non-fixed depen-dent edges, e nfd is defined by vector q nfd as q nfd = (cid:40) q i if i is the index of an e ndf q i of the e nfd edges can be computedfrom the lengths of the edges corresponding to e ci , e nci ,and e fix . The edge lengths will be updated in order tosatisfy the area equation.After setting up the notations, we begin to solve the areaEq. 18.8ince any edge is either e ci , e nci , e fix or e nfd , we have anequality q = q ci + q nci + q fix + q nfd . (23)Moreover, the lengths of the e nfd depend linearly on thelengths of the e ci , e nci and e fix , hence there exist an [ e × e ]square matrix, D and vectors d , g so that: q nfd = Dq nci + d q ci + g (24)The matrix D and the vectors d and g can be computed fromthe RREF form easily as follows.The matrix D is a matrix whose entries are mostly 0 exceptat the entries corresponding to the columns of e nci and tothe rows of the e nfd , where the value is the opposite of thecorresponding value in the RREF f matrix of (cid:0) B f b f (cid:1) . The vector d is a vector whose entries are mostly zero (0)except at the entries corresponding to the column of the e ci edge and to rows of the dependent edges, where the value isthe opposite of the corresponding value in the RREF f matrixof (cid:0) B f b f (cid:1) . The vector g is the contribution coming from the fixededges. This is a vector whose entries are mostly zero (0),except for entries corresponding to the indices of the e nfd edges, when the entry is the last entry of the correspondingrow of the RREF f matrix (cid:0) B f b f (cid:1) . We simplify Eq 23 slightly. Consider the [ e × e ] squarematrix Id nci , which is the identity restricted to the e nci edges,and 0 elsewhere. Since, Id nci q nci = q nci we have that q nfd + q nci = D (cid:48) q nci + d q ci + g where D (cid:48) is the matrix D + Id nci .Thus, by Eq. 23, we have q = q ci + D (cid:48) q nci + q fix + d q ci + g . Let us denote d (cid:48) by the vector obtained by adding a 1 to thevector d at the index of the critical edge, i.e d (cid:48) q ci = q ci + d q ci .Then, we have q = D (cid:48) q nci + q fix + d (cid:48) q ci + g (25)Now, we can solve the area Eq. 18 by plugging in the right-hand side of Eq. 25 into q : the quantity k ( D (cid:48) q nci + q fix + d (cid:48) q ci + g ) T M f ( D (cid:48) q nci + q fix + d (cid:48) q ci + g ) (26)computes the area of the face, A f . Rearranging the terms, weobtain a quadratic equation for q ci : aq ci + bq ci + c = 0 (27)where a = d (cid:48) T M f d (cid:48) (a) A f = 1357.06 (e.g.) v v v v v (b) A = -1557.31A = = -377.63 v v v v v A f = A + A + A = 1357.06 n f n f n f n f n f n f Figure 6: There are at least two significantly polygon differentgeometries that represent the same area of a polygon.and b = 2 d (cid:48) T M f ( D (cid:48) q nci + q fix + g )and c = ( D (cid:48) q nci + q fix + g ) T M f ( D (cid:48) q nci + q fix + g ) − kA f We can solve this quadratic equation (using the quadraticformula) to obtain possibly two solutions for q ci : q ci = − b ± √ b − ac a Number of solutions
Depending on the target area, A f , we might have differentnumber of solutions for Eq. 27. It is possible to have no solu-tion, a unique solution, or two significantly different solutions(see Figure 6).Depending on the signature of A , a large positive or a smallnegative prescribed area ensures that we have multiple (two)solutions. Once we computed the length q ci of the e ci edge, we can up-date the lengths of the dependent edges using Eq. 24. Now,all the lengths of the edges are computed. The face corre-sponding to the edge lengths has the required area and theedges of the face satisfy the constraint equation system, Eq.22 while only the lengths of the e nfd edges and the length ofthe e ci were manipulated (Figure 8).The above discussion is summarized in Algorithm 2 in Sec-tion 3.2.9 v v v v v v v v v (b)(a) n f n f A f = ∑ A i = A = - = n f n f A = - 1.94 v v v v v v v v v v v v v (b)(a) A f = ∑ A i = A = - = n f n f A = - 1.94 v v v Figure 7: (a) Multiple steps to compute the new face geometrywith the preassigned area (for visualization purposes only);and (b) the computed face geometry with zero area.
Consider the pentagon of Figure 3. The matrix M f for thisexample is the following: − . − .
85 0 .
506 00 0 − . − .
358 0 . .
925 0 0 − . − . − .
012 0 .
679 0 0 − . − . − .
761 0 .
288 0 0
Here, the order of the edges is given as e (0 , , e (1 , , e (2 , , e (3 , and e (4 , . The user selects e (4 , as the fixed edge andassigned area of the face to be zero.The matrix B f is given as .
289 1 0 . − . − . − .
957 0 0 .
631 0 . − .
380 0 0 0 1 and the vector b f is . . Here the first two rows describe the equilibrium equations,Eq. 2, –for the x − and y − coordinates. The third equationis the constraint equation, Eq. 21, for the fixed edge whoselength is 41 .
78. In our case, the fixed edge corresponds to thelast entry.The
RREF f matrix is − . − .
709 0 − . . − .
529 0 43 . . . (28)As a consequence, the CDGoF f can be computed as e − r =5 − (c) A f = e FIX = e , e (d) A f = e FIX = e , e (b) A f = e FIX = e , e (a) A f = ∑ A i = e FIX = e , e v v v v v v v v v v v v v v v v v v v v v v v v v v v v v e e e e e e e e n f n f n f n f (c) A f = e fix = e , e (d) A f = e fix = e , e (b) A f = e fix = e , e (a) A f = ∑ A i = e fix = e , e v v v v v v v v v v v v v v v v v v v v v v v v v v v v v e e e e e e e e n f n f Figure 8: (a) to (d) various zero-area computation for a start-ing face with the area A f and different chosen fixed edges.Moreover, we can identify the e ci , e nci , e fix and e nfd edgesfrom the RREF matrix, 28 as follows.The fixed edge e (4 , corresponds to the fifth entry. Since,the CDGoF f equals to 2, we have two more (5 − − e nfd edges, e (0 , and e (1 , ,given by the other pivotal columns. The e ci edge was chosento be the edge corresponding the fourth entry, e (3 , . Finally,the e nci edge is the remaining edge, e (2 , for which we solvethe quadratic area equation (see Figure 5).Now, we compute the coefficients of Eq. 27 to solve for theedge length, q ci of the e ci edge. Here, the target area, A f iszero. a = − . , b = − . , c = − . q ci = − .
974 and 212 . . As a consequence, we get two significantly different solutionsfor the geometry of this face. The updated (self-intersecting)face corresponding to q ci = − .
974 is shown in Figure 7.
In this process, a user would select multiple internal/externalfaces and edges of a polyhedral system and would assign target10reas for each face and edge lengths for each edge to computethe new geometry of the polyhedron. Computing the geome-try of a system of polyhedral cells with pre-assigned areas andfixed edge lengths in one step is a complex task. We proposea multi-step, inductive process to tackle this problem.
In each step, we compute the geometry of a single face f i withthe assigned area as described in Section 2.4. Then we up-date the polyhedron with the new edge lengths using Moore-Penrose Inverse (MPI) Method and move to the next face andrepeat this process until there is no face left to change.First, we identify the fixed edges from the list e fix whichlie on the face f i with prescribed edge lengths, and use RREFmethod to solve the quadratic area Equation 18 described inSection 2.4.In the next step, to preserve the new geometry of the face f i , we consider all the newly generated edges of the face asfixed edges for the entire polyhedral system. I.e., we updatethe list of fixed edges e fix for the entire polyhedron with thenewly computed edge lengths of f i . Similarly to Section 2.4, Eq. 22, the linear constraint equa-tions for a polyhedral system can be described by two differentkinds of linear equations. First, we have the equilibrium equa-tion system, Eq. 3 describing the topology of the polyhedralsystem. Second, we have the linear constraints coming fromthe prescribed edge lengths.As a result, we obtain an equation system that describesthe equilibrium of the polyhedron with the prescribed edgelengths of the fixed edges: B p q = b p (29)where the constraint matrix B p = (cid:18) E p L p (cid:19) is built from the equilibrium matrix of the polyhedron E p (seeEq. 3) and the constraint equations L p q = l p coming fromthe fixed edges (see Eq. 21). Similarly, the vector b p = (cid:18) p (cid:19) is obtained from the edge vector of the fixed edges l p . Now, we propose to solve Eq. 29 using the Moore-PenroseInverse (MPI) method. The MPI method is a technique to solve a general non-homogeneous equation system of the ma-trix form B p q = b p (30)where the matrix B p and the vector b p are given.We represent MPI of the matrix B p by B + p that satisfiesthe following matrix equations B p B + p B p = B p and B + p B p B + p = B + p . Assume that the vector b p is of the form B p q for somevector q . Multiplying the first equality by q , we have B p B + p B p q = B p q . As a consequence, we obtain B p B + p b p = b p . (31)Therefore, if a solution to Eq. 30 exists, then the Eq. 31 hasto be satisfied. Similarly, if Eq. 31 is satisfied, then the vector B + p b p is a solution to Eq. 30 This provides an effective toolto check whether Eq. 30 has a solution or not.From now on, we assume that Eq. 31 holds, in other words,we assume that Eq. 30 has a solution. In this case, any vector q of the form q = B + p b p + ( Id − B + p B p ) ν (32)solves the linear equation system Eq. 30 where Id is the iden-tity matrix and ν is any column vector of the right dimension.In fact, these are all the solutions to Eq. 30. Summarizingthe above discussion, we have at least one solution to Eq. 30if and only if Eq. 31 holds for b p . Moreover, if there is asolution to Eq. 30, then all solutions have the form of Eq. 32[20, 22].In Eq. 32, the parameter ν is freely chosen by the user tocontrol the solution. In our examples, we take ν to be theinitial edge lengths of the polyhedron, resulting in a solutionto Eq. 32 which is the new geometry for the initial polyhe-dron with the prescribed area for face f i . In this case, thenew edge lengths are the best fit (least squares) to the initialedge lengths. Also, in many cases, only certain parts of thepolyhedron change significantly (see Figure 9 and 12).Another approach could be to take ν to be the vector whoseentries are all 1, in this case, we get a solution with well-distributed edge lengths. The previously discussed method can compute the geometryof a polyhedral system with multiple faces with prescribedareas in an inductive process. In each step, we update thegeometry of the polyhedron using Eq. 32 with the newlycomputed face whose edge lengths are added to the list offixed edges. The new edge lengths change the constraintsequations L p q = l p to compute the polyhedral geometry.We summarize the process in Algorithm 3.11 m) (n) (o) A f0 = A f1 = e FIX = e , e (l) A f0 =
0, A f1 >
0 , e FIX = e , e (k)(j) A f0 , A f1 > f0 =
0 , e FIX = e (h)(g) (f) A f0 = f1 < e FIX = e (e)(d) A f0 , A f1 > f =
0 , e FIX = e (b)(a) A f > v v v v v v v v v v v v e e e f n f n f - n f n f - n f v v v v v
0, 1 v
5, 4 e e v
0, 1 v
5, 4 e e v v v v e e - n f1 - n f1 n f1 n f0 n f1 n f0 v
0, 3 v
1, 2 v
0, 3 v
1, 2 v v v v v v v v v v v v v v v v n f1 n f1 n f0 v v v v n f1 v
4, 5 v
6, 7 e e e e e e e e e e n f (m) (n) (o) A f0 = A f1 = e fix = e , e (l) A f0 =
0, A f1 >
0 , e fix = e , e (k)(j) A f0 , A f1 > f0 =
0 , e fix = e (h)(g) (f) A f0 = f1 < e fix = e (e)(d) A f0 , A f1 > f =
0 , e fix = e (b)(a) A f > v v v v v v v v v v v v e e e f n f n f - n f n f - n f v v v v v
0, 1 v
5, 4 e e v
0, 1 v
5, 4 e e v v v v e e - n f1 - n f1 n f1 n f0 n f1 n f0 v
0, 3 v
1, 2 v
0, 3 v
1, 2 v v v v v v v v v v v v v v v v n f1 n f1 n f0 v v v v n f1 v
4, 5 v
6, 7 e e e e e e e e e e n f Figure 9: (a), (d), (g), (j), (m) Multiple polyhedral geometries with selected faces (orange) and user-assigned fixed edges;(b), (e), (h), (k), (n) the face area computation and visualization in multiple steps (for visualization purposes only); and (c),(f), (i), (l), and (o) the resulting polyhedral geometries with zero areas for the selected faces.12 k e i (c) (d)(b)(a) e i n k e j † = n j . q j † , ∑ e i † = j f i † Є Γ † f i † f i † e i † m e i Є Γ e i Figure 10: (a) An input force polyhedron as primal and itscorresponding (b) funicular polyhedron as the dual; (c) goingaround each edge of the primal with its attached faces (c)provides the direction of the edge vectors of the correspondingface (e) in the reciprocal diagram where the sum of the edgevectors must be zero.
In [3, 9] three different methods were described to generatethe dual diagram from a given primal. Let us call the startingdiagram the primal , Γ, and the reciprocal perpendicular poly-hedron dual , Γ † (Fig. 10a, b). The vertices, edges, faces, andcells of the primal are denoted by v , e , f , and c respectively,and the ones of the dual are super-scripted with a dagger ( † )symbol (Figure 10a,b).Since the face f i † is a closed polygon, the sum of the edgevectors e j † should be zero. Hence, we obtain a vector equationsimilar to Equation 1 (cid:88) f j u j † q † j = where the sum runs over the attached faces f j of the edge e i of the primal Γ; u † j denotes the unit directional vector corresponding to the edge vector e † j ; and q † j denotes the edgelength of e † j in the dual Γ † .Similarly, as before, each vector equation for the face ofthe dual diagram yields three linear equations for the edgelengths, and we obtain a linear equation system for the edgelength vector q † which can be described by a [3 e × f ] matrixthat we call the equilibrium matrix of the dual diagram, E † , E † q † = . (33)In [3, 9] three different methods were described to generatethe dual diagram from Eq. 33. In this paper, we chooseto use the Moore-Penrose Inverse (MPI) method to initiallyconstruct the dual diagram before we apply any changes tothe force and its face areas. The MPI method of this sectionis as same as the method described in Section 2.5.3 with b p being the zero vector. As a result, the solutions to Eq. 33 canbe described as q † = (cid:16) Id − ( E † ) + E † (cid:17) ξ. Here ( E † ) + denotes the MPI of the equilibrium matrix E † .For the parameter ξ , we choose the vector whose entries are all1 to obtain a dual diagram with well-distributed edge lengths. The direction of the internal forces
The initial direction of the internal force as compression ortension is stored and altered after the computation of theforce with prescribed areas. The tensile force members areupdated in the form if the normal of a face in the force di-agram flips after the computation. As shown, the geometryof a face can become self-intersecting in some cases. On suchoccasions, if the area of the region with the initial normal di-rection is bigger than the area with the flipped normal, thenthe direction of the initial internal force does not change; oth-erwise, the direction of the internal force will flip. If the faceis a zero area face, the member will carry no force and candisappear in the form diagram (Figure 12.)
In this section, we explain the computational setup for themethods we explained in Section 2. The input for this frame-work is a polyhedral system with planar faces and is con-sidered as the force diagram for the methods of 3D graphicstatics. we compute the dual geometry of the updated primaldiagram according to [3, 9].The user can initially select certain edges in the systemand assign a target length per selected edge. Similarly, s/hecan select multiple faces and assign an area per face. Ourmethod is a sequential computational setup that updates thegeometry of the polyhedral system for each user-selected faceat each step. For instance, let’s assume that the user selects13 tartInput Primal geometry v i , e i , f i area matrix M f i Equlibrium matrix E f i Solve area equation for q ci Compute the new geometry for face f i Update Primal geometry E p † q † = StopU/ISelect faces f i-k & define target area A f i Select edges e fix (fixed or predefined lengths) f ≠ Ø Yes No i < k CGDoF f CGDoF f > = 0 Yes No Yes
RREF( B f i | b f i ) If there is a solution for q ci U/Iupdate A f i update e nfd of the face f i e ci , e nci , e fix and e nfd obtain { L f i q = l f i B f i q = b f i Update B p q = b p solve B p q = b p using MPI method compute CGDoF p CGDoF p ≥ YesYes No
Updating the Dual geometry retrieve e fix of f i from e fix Yes No > = i = i +1 Equilibrium matrix E p and GDoF p (E p ) e ci , e nci , e fix and e nfd e fix e fix ≠ Ø NoStop NoStop
Figure 11: The flowchart expanding multiple steps in com-puting the primal geometry/force diagram with preassignededge lengths and face areas, and the updated dual geometryas the form diagram. three faces with a target area per face. We start from a face,compute the new geometry of the face with the target area,update the geometry of the polyhedral system based on thenew geometry of the face, and move on to the next face andcontinue the computation until there is no face left.Accordingly, in each step of the computation, we constructthe area matrix and the equilibrium matrix for the user-selected face. We, then, add the user-defined edge constraintsas non-homogeneous linear equations and compute the Con-strained Geometric Degrees of Freedom,
CDGoF using RREFmethod explained in Section 2.4. In the next step, we computethe geometry of the face, and then we update the geometry ofthe polyhedral system using MPI method described in Section2.5.After the multi-step computation process is completed, weupdate the direction and the magnitude of the forces in themembers of the dual that was initially constructed.The above description can be summarized into three mainsections as shown in the flowchart of Figure 11. These sectionsare as follows: • computing the new geometry for a face with constrainededges and areas; • updating the new geometry of the polyhedral systembased on the newly-computed face geometry and thefixed edges; and, • updating the internal forces in the members of the dualbased on the new force magnitudes.Sections 3.1, 3.2, and 3.3 provide additional details of thealgorithms used in this process. These algorithms include:computing the area matrix for a face; face reconstruction withconstrained edges and target area; and updating the geometryof the polyhedron with constrained areas of faces and edgelengths. Algorithm 1 receives a face f i with an ordered list of vertices[ v , ..., v k − ] as an input and outputs a symmetric matrix M f i which is used in the quadratic form, Eq. 18, to compute thearea of the face.First, we choose an (arbitrary) normal vector for the face,by taking the unit cross product of the first two consecutiveedge vectors. Then, we construct the matrix M f i row byrow as follows: starting from each vertex v n , we create anordered list of directed edges e i , and compute the scalars η i,j as explained in Eq. 9.Once the whole matrix M f i is constructed, the algorithmoutputs a symmetric matrix (see Eq. 19) to be used in thequadratic form for computing the area of the face f i .14 lgorithm 1: Computing the area matrix M f i Input: f i : [ v , v n , ..., v k − ] the ordered list of verticesaround f i Output: M f i the area matrix of the face f i . Function M row ( v row , n f ) : v row : [ v j , v j +1 , ..., v j − ] v j for v n ∈ v row do e ←− e i [ v n , v n +1 ] for e i ∈ e doe i ←− (cid:104) x n +1 − x n , y n +1 − y n , z n +1 − z n (cid:105) v n to v n +1 | e i | ←− l i e i for e i ∈ E don ,i = e × e i e and e i if x n f (cid:54) = 0 x coordinate of n f then η i = x n ,i / ( x n f ∗ | e i | ∗ | e | ) elseif y n f (cid:54) = 0 y coordinate of n f then η i = y n ,i / ( y n f ∗ | e i | ∗ | e | ) if z n f (cid:54) = 0 z coordinate of n f then η i = z n ,i / ( z n f ∗ | e i | ∗ | e | ) M e i = ( k − ∗ i − ∗ η i e i where k is the length of V row M row ←− M e i return M row begine ←− v to v n f i ←− e × e f i for v i ∈ { v , ..., v k − } do v v i ←− [ v i , v i +1 , ..., v i − ] v i M v i = M row ( v v i , n f i ) M f i ←− M v i M v i M f i = (cid:16) M f i + M Tf i (cid:17) The Algorithm 2 updates the geometry of a face with con-strained edges and a target area. The input of this algorithmis a user-selected face f i with a target area A f i and a (user-selected) edges of e fix with prescribed edge lengths.First, we compute the linear constraint equations, Eq. 22,where the equilibrium equations describe the geometry ofthe face and the constraint equations come from the (user-selected) edge constraints.Once, the constraint equation system, Eq. 22, is cre-ated, we compute the Constrained Geometric Degrees of Free-dom of f i , CDGoF f i , of the face using RREF method. If CDGoF f i = −∞ , the algorithm stops, i.e., the constrainedequation system, Eq. 22, cannot be solved. In this case, theuser may modify the input by selecting less constrained edgesor by selecting different edges.If CDGoF f i is at least zero, the algorithm classifies theedges of f i into ci , nci , fix and nfd edges. Next, the weconstruct the area matrix M f i using Algorithm 1. Using thearea equation, Eq. 18, and Eq. 27, we compute the edgelength q ci of the edge e ci . If Eq. 27 has no solution, thealgorithm may ask the user to modify the target area. TheEq. 27 often has two solutions, the user may choose fromthose particular solutions (see Figure 6).Once a solution is chosen for q ci , the algorithm updatesthe lengths of the e nfd (see Section 2.4.1) and outputs theupdated lengths of the edges of f i . The last algorithm updates the geometry of a polyhedron af-ter the new geometry of each face is computed. In fact, thisis a multi-step process, in each step, we update the geome-try of the polyhedron after updating the geometry of a face.Note that, in each step, the list of fixed edges is updated af-ter computing the geometry of each face, in other words, wesolve for the geometry of the polyhedron with more and moreconstraints.The Algorithm 3 describes the computation of the new ge-ometry of a polyhedron with a given list of constrained edges, e fix . It computes the linear constraint equation system, Eq.29. This is a non-homogeneous linear equation system whichcan be solved using MPI method. The parameter ν in theMPI method can be chosen by the user or can be the vectorof the initial values of the edge lengths of the polyhedron. Figures 9, 12, and 13 illustrate the potentials of using this ap-proach in polyhedral transformation with target areas of faces15 a) (a) (a) Џ † (a) Џ † (a) e † = n . q , e † = n . q (a) e † = 0 , e † = - n . q (a) (a) (a) Џ (a) Џ (a) A , A >0 (a) A = 0 , A < 0 n e † n f f = 0 n e † n -n -n e † e † n e † n e † e † e † -n -n (d) (e) e † = - n . q † , e † = - n . q † (f) Γ † (a) Γ † (b) e † = n . q † , e † = n . q † (c) e † = 0 , e † = - n . q † (a) A , A , A >0 A = 0, A , A < 0 Γ Γ A , A >0 A = 0 , A < 0 n e † n f f = 0 n e † n -n -n e † e † n e † n e † e † e † -n -n Figure 12: (a) Initial force polyhedron of the Figure 9j (top) and its reciprocal form (bottom) as a compression-only system;(b) highlighted internal faces of the polyhedron before transformation (top) and the corresponding members in the form(bottom) ; (c) zero area faces and their updated normal directions and the resulting zero-force members in the form; (d), (e)and (f) highlighted faces in the second step of the transformation and the updated form with new internal force distribution.16 lgorithm 2:
Updating the geometry of the face (UF)
Input: f i : [ v , ..., v k − ] face with ordered list of vertices A f i target area for the face e fix : [ e m , ..., e q ] list of constrained edges Output: Q : [ q , ..., q k − ] list of edge lengths of f i witharea A f i . begin e ←− [ e , ..., e k − ] for e i ∈ e doe i ←− (cid:104) x n +1 − x n , y n +1 − y n , z n +1 − z n (cid:105) v n to v n +1 | e i | ←− l i e i u i = e i / | e i | e i E x ←− x u i x -coordinates of u i E y ←− y u i y -coordinates of u i E f i ←− E x , E y × e ] equilibrium matrix ofthe face f i for e i ∈ e fix dol i , q i B f i , b f i RREF (( B f i | b f )) CGDoF f i if CGDoF f i = −∞ then no solution = ⇒ end program or ask user tomodify the input. elseq nci D , d , g D (cid:48) , d (cid:48) ←− D , d M f i : a, b, c M f i q ci ←− a, b, c q nfd ←− q ci , q nci , q fix Q = [ q , ..., q k − ] Algorithm 3:
Updating the geometry of the polyhedron
Input: e fix : [ e m , e n , ..., e q ] : list of fixed edges Output: Q p : [ q , ..., q e ] the updated list of edge lengthsof the polyhedron beginfor e i ∈ e fix dol i , q i B p , b p B + p ←− MPI( B p ) Q p ν (Eq.29)and constrained (selected) edge lengths, and their correspond-ing structural forms. In the examples shown in Figure 9, theintention is to highlight certain properties of the constrainedpolyhedral computation. In all the examples the chosen faceis highlighted by an orange shade and the constrained edgesare highlighted by blue color.For instance, in Figure 9a, the target area for the chosenface is zero and an edge has been chosen that does not belongto the selected face. Figure 9b shows an animated drawingfor clarification purposes. As it is shown in Figure 9c, theselected face turns into a self-intersecting face. In the nextexample, in Figure 9d, the top face is set to be zero withthe same constrained edge. As a result of this computation,the face collapses to a line due to its rectangular geometry.Also, the normal of the side faces, n f , flips which means thatthe magnitude of the internal force in the corresponding formwill change (9h, i). In the Example of Figure 9g, one of thevertical side faces with a constrained top edge is chosen andthe target area is also set to zero. As illustrated in Figures9h and 9i, the rest of the vertical faces will disappear afterthe polyhedral computation. This method can certainly beapplied to more than one face in the polyhedral system. InFigure 9j, one external and one internal face are chosen witha zero area target and as illustrated in Figures 9k-o, the com-putation process proceeds sequentially by first, computing thegeometry of the face f in Figure 9l and then recomputing thepolyhedral geometry to solve for the new face f . In this pro-cess, multiple other faces will also turn into a self-intersectingface as shown in Figure 9l. In the final geometry, the face f will collapse to an edge e (Figure 9o).In most of these examples, the target areas were intention-ally set to be zero to highlight its results in the reciprocalstructural form. Figure 12a shows the same polyhedron ofthe Figure 9j. In this figure, certain faces in the force di-agram and their corresponding edges in the form diagramare highlighted to show the effect of changing the areas onthe magnitude of the internal forces. Starting from the force17olyhedron of Figure 12a (top) and its compression-only form(bottom), faces f and f are emphasized in Figure 12b (top)and their corresponding compression-only edges in the form.The result of the zero area computation results in face f toturn into a self-intersecting face together with other similarfaces attached to the edges of the face f . Besides, face f is flipped as well as the direction of its normal vector. Notethat as a result of this transformation, the internal force inthe corresponding member of the face f is decreased to zero.This is a fascinating effect in the equilibrium of polyhedralframes as it describes the internal equilibrium of forces in apolyhedral system where the edges or the members can beremoved from the system without disturbing the internal andexternal equilibrium. The zero-force edges have been pre-viously observed in some polyhedral reciprocal diagrams by[16], but there was no method to compute them particularlyin self-intersecting faces as described in the methodology sec-tion. Also, the change in the direction of the normal of theface f results in reversing the internal force in the edge e ofthe form diagram.Figure 12d emphasizes the faces f , f , and f and facesthat will be affected as a result of the second step of thecomputation. As it is shown in Figure 12e, the normal of thefaces f and f will invert the direction of the internal forcesin the form diagram. This transformation also removes theapplied load in the system as the area of the face f is zero. Inanother example, the zero area faces are used to completelyremove the external horizontal forces in the system. Figure 13shows a force diagram and its corresponding form in anothertransformation. In this example, one of the vertical facesis chosen and the target area is set to become zero. Notethat as animated in the drawings of Figure 13b, all the sidefaces collapse into a line and the top and bottom faces ofthe polyhedron will become coplanar. This transformationresults in the disappearance of all the horizontal applied loadsin the system as shown in Figure 13f. The most interestinggeometric outcome of this process is the transformation of theinternal face f in this process. Face f changes its directionand so does its corresponding edge at the boundary of theform diagram. The resulting structure shows a funnel shapecompression-only structure with tensile members on the top. This paper provides an algebraic formulation alongside withalgorithms, and numerical methods to geometrically controlthe areas of the faces of general polyhedrons of the recip-rocal diagrams of 3D/polyhedral graphic statics. The pre-sented methods bridge the gap between the previously devel-oped algebraic methods for the construction of the reciprocalpolyhedral diagrams and controlling the magnitude of inter-nal and external forces by changing the areas of the faces.This method for the first time allows the user to manipulate both convex and complex faces and explore the compressionand tension combined features in structural form-finding us-ing 3DGS. Controlling the areas of complex faces has neverbeen addressed in the literature prior to this research as theprevious approaches mainly dealt with convex polyhedrons.Thus, this research opens a new horizon understanding theequilibrium of both tension and compression forces beyondthe existing compression-only polyhedral funicular forms.The paper explains the process of turning geometric con-straints such as edge lengths and target face areas of the re-ciprocal polyhedral diagrams into algebraic formulations com-patible with the previously developed method by [9] and [3].This research describes a quadratic formulation to computethe geometry of a face with a target area and provides a linearformulation to consider the edge lengths as constraints. Solv-ing an equation system including both linear and quadraticequations is a highly complex task. The key idea in our pro-posed method is to reduce the number of unknowns in the(quadratic) equation system of a face using Reduced RowEchelon (RREF) method. Computing the updated geome-try of the polyhedral diagrams is achieved by Moore-PenroseInverse (MPI) method. In this approach, multiple faces andedges can be selected as constraints, and the new geometryof the polyhedral system is computed in a sequential process.The paper also describes the Constrained Geometric De-grees of Freedom (CGDoF) of the linearly constrained poly-hedral systems and opens up a door to explore a wide vari-ety of interesting geometries satisfying the initial equilibriumequations with selected edge lengths. The algorithms andnumerical methods provide an interactive tool for the user tostudy and manipulate large-scale general polyhedral diagramsby assigning face areas and edge lengths.
Future work
The existing algorithm deals with each face at each step andadds the newly computed edge lengths to the initial con-straints of the polyhedral system. As a consequence thereis no control over the number of constraints generated in eachstep and the eventual number of constraints to compute theentire polyhedral group. Therefore, in certain cases, depend-ing on the chosen edge by the user, or the geometric degreesof freedom of the entire system, the polyhedral computationbecomes over-constrained. This property proposes an inter-esting problem for future research.Another interesting future direction is the study of thedifferent solutions of the same initial, constrained problem.As mentioned in Section 2.4, for a single face there can betwo, significantly different polygons satisfying the linear andquadratic constraint equations. As a consequence, for a poly-hedral system with n assigned face areas, there can be 2 n significantly different updated polyhedral systems that canalso be explored in future research.18 f) e † = - n . q , | f | = A f (e) e † = n . q (d) Џ † | f | = A f (d) Џ A f > f > f < f > f = e v e v v v v
1, 2 v
0, 3 v v e v
1, 2 v
0, 3 v v v v v
2, 1 v v e † e † n f n -n f = 0 n (f) e † = - n . q † , | f | = e † = n . q † (d) Џ † | f | = A f Џ A f >
0 A f > f < f > f = e v e v v v v
1, 2 v
0, 3 v v e v
1, 2 v
0, 3 v v v v v
2, 1 v v e † e † n f n -n f = 0 n Figure 13: (a) A force polyhedron and a selected face with zero area target, and two fixed edges on the top; (b) thetransformation animation (for visualization purposes only); (c) the resulting force polyhedron where the selected face and allother side faces collapse to a line; (d), (e), and (f) the force and form diagrams and their transformation after changing theareas. 19
Acknowledgment
This research was partially funded by National Science Foun-dation Award (NSF CAREER-1944691).
References [1] M Akbarzadeh.
3D Graphic Statics Using PolyhedralReciprocal Diagrams . PhD thesis, ETH Z¨urich, Z¨urich,Switzerland, 2016.[2] M. Akbarzadeh and M. Hablicsek. Geometric degrees offreedom and non-conventional spatial structural forms.In
Impact: Design With All Senses, Design ModellingSymposium , Berlin, Germany, September 23-25 2020.[3] M. Akbarzadeh, M. Hablicsek, and Y. Guo. Developingalgebraic constraints for reciprocal polyhedral diagramsof 3d graphic statics. In
Proceedings of the IASS Sympo-sium 2018, Creativity in Structural Design , MIT, Boston,USA, July 16-20 2018.[4] M Akbarzadeh, T Van Mele, and P Block. Equilibriumof spatial networks using 3d reciprocal diagrams. In J.B.Obr ebski and R. Tarczewski, editors,
Proceedings of theInternational Association for Shell and Spatial Structures(IASS) Symposium 2013 , Wroclaw, Poland, September2013.[5] M Akbarzadeh, T Van Mele, and P Block. 3d GraphicStatics: Geometric construction of global equilibrium.In
Proceedings of the International Association for Shelland Spatial Structures (IASS) Symposium , Amsterdam,August 2015.[6] M Akbarzadeh, T Van Mele, and P Block. On the equi-librium of funicular polyhedral frames and convex poly-hedral force diagrams.
Computer-Aided Design , 63:118–128, 2015.[7] M Akbarzadeh, T Van Mele, and P Block. Spatialcompression-only form finding through subdivision ofexternal force polyhedron. In
Proceedings of the In-ternational Association for Shell and Spatial Structures(IASS) Symposium , Amsterdam, August 2015.[8] A Beghini, L L Beghini, J A Schultz, J Carrion, andW F Baker. Rankine’s theorem for the design of cablestructures.
Structural and Multidisciplinary Optimiza-tion , 2013.[9] Mrton Hablicsek, Masoud Akbarzadeh, and Yi Guo.Algebraic 3d graphic statics: Reciprocal constructions.
Computer-Aided Design , 108:30 – 41, 2019.[10] B. K. P. Horn. Extended gaussian images.
Proceedingsof the IEEE , 72(12):1671–1686, Dec 1984. [11] J. Lee.
Computational Design Framework for 3D GraphicStatics . PhD thesis, Zurich, 2018.[12] J. Lee, T. Van Mele, and P. Block. Area-controlled con-struction of global force polyhedra. In
Proceedings of theIASS Symposium 2017 , Hamburg, September 2017.[13] J. Lee, T. Van Mele, and P. Block. Disjointed force poly-hedra.
Computer-Aided Design , 99:11–28, June 2018.[14] James J. Little. An iterative method for reconstruct-ing convex polyhedra from extended gaussian images. In
Proceedings of the Third AAAI Conference on ArtificialIntelligence , AAAI83, page 247250. AAAI Press, 1983.[15] J C Maxwell. On reciprocal figures and diagramsof forces.
Philosophical Magazine and Journal Series ,4(27):250–261, 1864.[16] A McRobie. Maxwell and Rankine reciprocal diagramsvia Minkowski sums for 2D and 3D trusses under load.
International Journal of Space Structures , 31:115–134,2016.[17] A McRobie. The geometry of structural equilibrium.
R.Soc. open sci. , 4(160759), 2017.[18] A McRobie, M Konstantatou, G Athanasopoulos, andL Hannigan. Graphic kinematics, visual virtual workand elastographics.
R. Soc. open sci. , 4(170202), 2017.[19] Shankar Moni. A closed-form solution for the reconstruc-tion of a convex polyhedron from its extended gaussianimage. [1990] Proceedings. 10th International Conferenceon Pattern Recognition , i:223–226 vol.1, 1990.[20] E H Moore. On the reciprocal of the general algebraicmatrix.
Bulletin of the American Mathematical Society ,26(9):385–396, 1920.[21] Joseph O’Rourke. Convex polyhedra realizing given faceareas.
CoRR , abs/1101.0823, 2011.[22] R Penrose. A generalized inverse for matrices.
Mathemat-ical Proceedings of the Cambridge Philosophical Society ,51(3):406–413, 1955.[23] M Rankine. Principle of the equilibrium of polyhedralframes.
Philosophical Magazine , 27(180):92, 1864.[24] G G Stokes.