Increasing efficiency of iterative application of the force density method

The method for reducing total computational time in iterative application of the force density method, intended to attain prescribed force values in cables or prescribed lengths of cable segments, is described in the paper. In each step of the iterative procedure, linear systems are solved with the accuracy that takes into account differences between the calculated and required values. The rule that prevents excessively fast increase in accuracy, while maintaining it high enough not to compromise the convergence of the iterative process, is proposed.


Introduction
Cables have negligible flexural stiffness. Consequently, when subjected to forces perpendicular to their axes, they change their configuration to enable development of equilibrating internal tensile forces. The shape of a flexible cable structure is maintained by arranging cables into a net that "describes" a surface with a negative Gaussian curvature. Additionally, cables should be prestressed for tensile forces to be preserved at all load combinations. The shape of cable net structures is determined by the laws of statics: designing the structural form is a search for the system of forces in equilibrium [1,2]. The form finding process is therefore the first stage in the design of prestressed cable nets (Section 2.1). It is a non-standard, inverse task, since in a "standard" design the shape of a structure is already specified so that "only" a compatible distribution of internal forces needs to be determined. Furthermore, computational model equations are highly nonlinear (Section 2.2). In the force density method, developed in the early 1970s for the design of the Olympic Games facilities in Munich, the problem was linearized [3][4][5] (Section 3.1). However, the concept of force density, as a ratio of one static and one geometric magnitude, has arisen from formal manipulation of mathematical expressions. Hence it is often difficult to select a distribution of force densities that complies with aesthetic, functional and structural requirements. Using the nonlinear force density method, based on the least squares method and also introduced in article [3], some structural requirements can be satisfied in a systematic way: attainment of prescribed force values in selected cables or cable segments, attainment of prescribed cable segment lengths in the equilibrium net configuration, and attainment of prescribed unstrained cable segment lengths. An extension to the nonlinear force density method, enabling attainment of the prescribed support reactions, and indirectly assigning position of selected (free) nodes, is given in reference [6]. Another extension, described in [7], enables assignment of coordinates (all three, two, or just one) of some or all nodes. As indicated in references [1,2], an iterative application of the linear force density method enables either partial or complete attainment of prescribed force values in cables, while in [8] the iterative procedure is extended in order to attain prescribed lengths of the equilibrium configuration segments (Sections 2.3 and 3.2). An experienced designer can satisfy other, usually non-structural, requirements in a series of attempts through interactive work. An iterative application of the force density method sometimes demands timeconsuming calculations and can therefore be inappropriate for interactive work. The objective of our work is to reduce duration of iterative computations (Section 3.3, with programming implementation in Section 4, and examples in Section 5). Initial results of this research were presented at the 40 th Solid Mechanics Conference [9].

Form finding
The notion of form finding for a prestressed cable net implies determination of its initial equilibrium configuration, before application of service loads and without taking self weight into account. The notion of configuration involves both the geometric shape and distribution of prestressing forces in cables; in [10] the authors use the term "initial equlibrium problem" instead of the "form finding problem". A comprehensive review of various cable net form finding methods is given in reference [11]. An interesting idea involving application of graphical statics and reciprocal diagrams, which presupposes interactive work, is described in reference [12]. Taking into account some simplifying assumptions (explained in more detail in [8]), cable nets will be modeled as space trusses [1,3]: the points in which cables cross each other and where they are connected to the ground or to rigid structural elements are spherical joints, while cable segments between the crossing points are bars pinned at their ends ( Figure 1). Bar connecting joints i and j will be denoted {i, j}. Joints that are not supports will be called free nodes. Basic variables in the cable net form finding are: disposition of cables, positions of their crossing points in space, disposition and spatial position of supports, and values of prestressing forces in cables or their ratios (adapted from [10] for cable nets). The choice of the form finding procedure will depend on which of the variables are given and which are unknown. The disposition of cables, which determines the disposition and connectivity of their crossing points, i.e. the topology of the system of pinned bars in the computational model, is almost always given in advance. As the shape of the net is an equilibrium configuration of prestressing forces in the system of bars, the coordinates of free joints are the principal unknowns. Prestressing force values can be either given in advance, or obtained as unknowns by the conditions of nodal equilibrium or, possibly, by satisfying additional kinematic constrains. The shape of the net can be varied by changing the ratios of force values. The shape of the net can also be influenced by adding/removing supports or by altering their positions. GRAĐEVINAR 69 (2017) 12, 1075-1084 Increasing efficiency of iterative application of the force density method

Nodal equilibrium equations
Mathematical formulation of the cable net form finding is very simple: it all comes to equilibrium equations of free nodes that are acted upon only by prestressing forces in connected bars. The set of bars connected to the node i will be denoted b i . For each free node i, three scalar equilibrium equations, representing the vanishing of force projection sums on three coordinate axes, can be written as follows: (1) where S i,j denotes the value of the force in connected bar {i, j} and is the length of that bar. Since (2) these equations are nonlinear.
If the net contains n free nodes, the obtained system will have 3n equations of the form shown in (1). If b is the number of bar elements, the number of potential unknowns will be 3n+b, since every set consisting of n coordinate triples (x i ,y i ,z i ) and b force values S i,j , and satisfying these equations, forms an equilibrium configuration. Since equations (1) do not contain coefficients that express constitutive relations between cable extensions and force values, the form finding problem is a pure static problem.

Generalized minimal nets and kinematic constrains
Equilibrium equations (1) can be interpreted as a condition for the minimum of the function e (3) where N and b denote set of free nodes and set of bars, respectively. The nets with specified prestressing force values, which satisfy that condition, can be called generalized minimal nets. Namely, if the values of all forces are equal, S i,j = S, the solution of equations (1) is the shape for which the sum of cable lengths is smaller than in any other shape of the net with the same topology [1,2,8,13].
The possibility of assigning different force values in different cables increases the number of possible cable net forms. For example, if forces in other cables stay fixed, the selected cable will tighten as its force increases and its length will decrease. At the same time, the spatial polygonal line of the cable axis straightens and approaches the rectilinear connection of its ends (example in Sections 5.1 and 5.2).
To enable the (generalized) minimal configuration with specified force values, the sliding of cables one over another during the prestressing procedure should not be prevented ( Figure 2) [1, 2].

Figure 2. Generalized minimal net
However, it may happen that two or more nodes slide along a cable into a single point, regardless of the values of forces. The unconfined approaching and merging of nodes cannot be prevented by increasing the force value in that cable, compared to forces in sliding cables. In our computational model, the sliding along a cable can be prevented by specifying lengths of segments between the crossing points, i.e. the lengths of bars into which a cable is divided. But, if we specify the bar lengths, the force values required to accomplish them become the unknowns [8].
Using Lagrange multipliers l i,j , the specified bar lengths can be expressed as kinematic constraints, and so the function e becomes (4) where b c is the set of bars with specified length . Since the number of Lagrange multipliers l i,j is equal to the number of unknown force values, the condition for stationary points of the function e c is a system with an appropriate number of equations.
The application of Newton-Krylov methods in solving the system of equations (1) (function e minimum condition), and the system of stationarity conditions of function e c , is analysed in reference [8]. In the first case, the analysis confirmed the well-known fact already described in [2,14,15]: the duration of calculation depends upon the assumed initial form, and the domain of convergence ("basin of attraction") exhibits a highly irregular, fractal shape. On the other hand, the problem of

Initial version of the linear force density
By specifying ratios (5) for each bar nonlinear equations (1) are linearized: The force-to-length ratios q i,j are called force densities [3][4][5].
If instead of force values S i,j force densities q i,j are specified, unknowns are force values and free node coordinates, which gives b+3n unknowns. Unknown force values and unknown nodal coordinates are connected by additional b equations (5). But the specified values q i,j are constant coefficients in equations (6), and so the system of 3n equilibrium equations is decomposed into three independent systems with n unknowns. The unknowns of the first, second and third systems are {x k } k N , {y k } k N and {z k } k N , respectively. The systems have the same system matrix, but they differ in vectors on the righthand side, since their components are q i,s x s , q i,s y s and q i,s z s for s S, where S is a set of supports. By obtaining free node coordinates, we can calculate bar lengths i,j and thereafter force values, since (5) gives S i,j = q i,j i,j .

Iterated force density method
Each solution of the system (6), for any one of ∞ b possible distributions of force densities in bars of the net of a given topology, is an equilibrium configuration. However, it is difficult to predict the force density distribution that will give the imagined net shape or the required distribution of prestressing forces. An iterative application of the force density method will be called the iterated force density method (IFDM). In that procedure, the linear force density method is applied in each step, while force densities are determined in a given step based on the requirements and results from the preceding step.
To attain the force value in the bar {i,j}, the force density in the bar is calculated in the k-th iteration step according to the expression (7) where is the force value calculated in the preceding step. Force values can be different in different bars [1,2]. A special case are minimal or geodetic nets where force values are distributed uniformly [13]: .
Analogously, the required length of the bar {i,j} and the required distance of nodes i and j can be attained in the equilibrium configuration of a net by calculating the force density in the k-th step according to one of the following expressions (8) where is the bar length/node distance calculated in the preceding step [8].
Both the force value in the bar and the bar length cannot be assigned at the same time, but it is possible to assign the approximate force value and then, after few iteration steps, switch over to obtaining the length. The iteration is terminated when i (9) where t S and t l are the prescribed tolerances, which can differ in different bars.
Since every solution of equations (6) is an equilibrium configuration, the described iterative procedure converges towards the target solution through a sequence of equilibrium configurations. If the iterative procedure is terminated before the required tolerance is achieved, or even if it cannot be achieved, the net obtained is still in equilibrium. Another advantage of iterative procedure is that in addition to the required force values, the required lengths of cable segments can also be obtained without introduction of Lagrange multipliers. However, numerical experiments have shown that the examples in which force values and lengths of cable segments are given demand more iteration steps than the examples with an approximately equal number of segments in which only force values are prescribed. The differences are however not so large compared to Newton-Krylov methods. It is stated in reference [8] that the described procedure is very rapid. This claim needs to be rewritten to be more accurate. The procedure is faster than the Newton-Krylov methods, but in the initial version, for form finding of cable nets with a large number of segments in real time, it takes quite a long time to complete if a higher accuracy is needed. A possible solution for speeding up this process is shown in the next section.

Reduction of computation time
Three linear equation systems with equivalent system matrices are solved in each step. The LU decomposition, which enables simultaneous back-substitution of several right-hand side vectors, is therefore a particularly suitable solution method [16]. Since these matrices are sparse, a sparse variant of the LU decomposition should be used for larger systems [17]. Nevertheless, the limited fill-in of the matrices will still occur, and our numerical experiments have shown that it can be significant. Increasing efficiency of iterative application of the force density method Iterative procedures for solving linear systems do not require any fill-in of system matrices. The group of Krylov procedures is very efficient [18] although three systems must be solved separately when these procedures are applied. The entire procedure of iterative application of the force density method will be called the outer loop, while the iterative procedure for solving the linear system will be called the inner loop. The tolerance t eq needs to be specified to enable iterative solution of linear systems. This tolerance depends on tolerances for required bar lengths and bar forces, t ℓ and t S , respectively, contained in the criteria for the outer loop termination. Lengths and force values are functions of nodal coordinates. The error expected in the determination of these coordinates is greater then t eq . Namely, in the condition for termination of an iterative procedure, the tolerance t eq is compared with the norm of the residual, ‖r (k) ‖< t eq , where k is a step of the inner loop. Since the error of the solution and the norm of the residual are related by ‖e (k) ‖= ‖A -1 r (k) ‖ ≤ ‖A -1 ‖ ‖r (k) ‖, the error of the solution can be up to ‖A -1 ‖ times larger than the residual. In addition, the length calculation error can be up to 2√3 times larger than the error with which nodal coordinates are determined. Systems (6) should therefore be solved at least with the tolerance of (10) where a is an estimate of A -1 (as such estimates are rather crude, there is no need for the factor 2√3). The application of Krylov methods reduces duration of the outer loop significantly. First of all, in an outer loop step, the solutions obtained in the preceding step are taken as the first approximations of solutions of linear systems. And, more importantly, the idea of the Newton-Krylov or, more evocatively called, inexact Newton methods, is borrowed. This provides a balance between accuracy of the linear system solutions and the amount of computation done in a single step of the outer loop [19]. Namely, if the computed force values or bar lengths are far from those required, it makes sense to solve linear systems (6) only approximately, by reducing accuracy with an increase of the deviation from the required value. It should be noted that the aforementioned advantage of the iterative application of the force density method is thereby sacrificed: since linear systems are in almost all steps solved only approximately, the obtained configurations, except for those in the last steps, are not in an equilibrium. Therefore the tolerance with which the systems (6) are solved is gradually decrease, from the relatively loose initial tolerance t max to the specified tight tolerance t eq that provides the final equilibrium configuration. If we wish to solve the system in the k-th step of the outer loop with the tolerance that reflects "distances" of calculated force values and lengths from the required ones, t (k) needs to depend on errors and . Since these errors are not known before the systems are solved, we will use values from the preceding step, and . The dependences can be expressed as and (11) Let's assume, to begin with, that functions f S and f ℓ are linear: and (12) Also, we will take that the relationship between and is equal to the relationship between t eq and in the final step of the outer loop and, similarly, that the relationship between and is equal to the relationship between t eq and : and (13) Since and , we can write and (14) At the end of the outer loop, when and reach t S and t l , it should be and . But, it will never be or , errors will be either slightly higher or slightly lower than the specified tolerances. If errors are smaller than the specified tolerances, the condition (9) for the termination of the outer loop is satisfied and the system (6) will never be solved with the tolerance t eq . Therefore, slightly smaller values will be used: and (15) The analysis of examples has shown that the tolerance is increasing too rapidly. Therefore, quadratic functions are taken: and (16) The reduction rate of errors and will also be taken into account by considering the ratio of these values in two consecutive steps [20]: and (17) where h is the "dumping" that provides additional control. The tighter of the tolerances and should be taken and therefore (18) In the "inexact" procedure, the tolerance of the step is determined by the greater of the errors and . Besides, GRAĐEVINAR 69 (2017) 12, 1075-1084 Elizabeta Šamec, Krešimir Fresl, Maja Baniček the tolerance should not be smaller than t eq . Therefore, the constraint is introduced as follows: (19) We have noticed that in some examples the convergence of the outer loop is not uniform: it can happen that or . That is why it should be confirmed that the tolerance in k-th step is equal to or smaller than the tolerance of the previous step: (20)

About the computer programme
The genuine programme code in which the calculations were conducted is written in programming package SageMath [21]. The flow chart of the programme is shown in Figure 3. Linear equation systems are solved by functions splu() and cg() contained within the programming library SciPy [22] that is included in SageMath. Sparse matrices are represented using the csc_matrix class (Compressed Sparse Column Matrix) from the SciPy library.

Minimal net with rigid supports
The first example is the net with "rigid" supports: all cables have fixed nodes at both ends. The net is spread over a ground-plan area [-a,a] 2 . The heights of supports are given by the expression (21) The net has two families of 23 cables, which gives 23·23=529 inner nodes, 4·23=92 boundary nodes, and 46·24=1104 bars. Inner nodes are free, and so three systems of 529 equations with 529 unknowns need to be solved in each iteration step. Unit force densities are assigned to all bars in the first step. The resulting net is shown in Figures 4.a and 4.c. Force values in bars range from 1,668 to 2,903.
The minimal net, determined by iterative application of the force density method, is shown in Figures 4.b and 4.d. Force densities in each step of iterative procedures are calculated according to Increasing efficiency of iterative application of the force density method expression (7), where the unit value of i,j is taken for all {i, j} (since the shape of the net does not depend on force values, but only on their ratios, the same shape is obtained for any positive ). Final force densities range from 0,090 to 1,197. When solving the linear system using the sparse LU decomposition, the required tolerance t S = 10 -4 is achieved after 576 steps. The same number of steps of the outer loops is required to solve the system using the conjugate gradient method (CG). The tolerance with which linear systems have to be solved is specified as t eq = 5×10 -7 . If zero vectors ( ) are taken as initial approximations of solutions, the total number of inner steps for all three systems is 78258. However, if solutions from preceding steps ( ) are taken as initial approximations, the total number of inner steps is reduced to 34505. In the proposed "inexact" procedure, linear systems are also solved using conjugate gradients (ICG). The tolerance is achieved after 557 steps of outer loop with the total of 16201 steps of inner loops. In this and all other examples, h = 0.025 in (8). Results for this example (Ex. 1), as well as for the other three examples, are concisely shown in Table 1 at the end of this section.

Different cable forces
The shape of the net can be influenced by changing the ratio of force values in individual cables. The net with the same layout disposition and the same topology is shown in Figure 5. In the net in Figures 5.a and 5.c, unit force densities are assigned to all bars that belong to the convex family of cables, while the force densities in bars of the concave family of cables are equal to one half of the height of their supports. In the generalized minimal net in Figures 5 b) and d), force values in all bars of the concave family of cables are equal to one half of the height of their supports, while bars of the convex family have unit forces. Force densities are calculated according to expression (7) in the iterative procedure of the form finding of the generalized minimal net. Since the topology of the net is equal to the one given in the previous example, the number of equations is the same. The steps required to achieve tolerance t S = 10 -4 in various computational procedures, with t eq = 5×10 -7 in CG and in final steps of "inexact" procedure, are listed in Table 1, rows Ex. 2.

Net with edge cables
Cable net structures are often made with edge cables. Inner cables of the net are connected to edge cables, and edge cables are, at their end points, connected to "rigid" structural elements.
The only known coordinates are the coordinates of these few supports and so the form finding of edge cables becomes a part of the form finding of the entire net. The net shown in Figure 6 has four edge cables stretched between four corner points. Three corners are in the same horizontal plane and the fourth one is raised above it.  (7), while those in edge cable bars are calculated using expression (8).
To terminate the outer loop, both conditions (9) need to be satisfied, where t s = t ℓ = 10 -4 . The inner loop tolerance is t eq = 5×10 -7 . The number of steps required to achieve given tolerances in different computational procedures is listed in Table 1, rows Ex. 3. The resulting net is shown in Figures  6.b and 6.d.

Net over octagon
In the final example, the net has nine supports: eight on the edges and one inner "high point". Edge cables are stretched between the edge supports. The "ridge" and "valley" cables are stretched between the "high point" and edge supports. In addition to the lengths of edge cable bars, the lengths of bars of "ridge" and "valley" cables are also specified. The net contains 90 cables in total. As in the previous example, specified bar lengths are arithmetic means of the lengths obtained in the first step, in which force densities in the corresponding bars were several times larger than in other bars (Figures 7.a and 7.b). Lengths are specified for Tolerances that control outer loop are t s = 10 -4 and t ℓ = 10 -3 . The net has 832 free nodes and, therefore, three systems, each containing 832 equations, need to be solved. The tolerance in iterative procedures is once again t eq = 5×10 -7 . The resulting net is shown in Figures 7 b) and d), while required steps are listed in Table 1, rows Ex. 4.

Conclusion
The form finding of prestressed cable structures is usually conducted in a series of attempts in which the designer tries to meet a set of architectural and structural requirements. The described iterative application of the force density method facilitates attainment of specified force values in bars or specified bar lengths, but in some cases the computation can be time-consuming and therefore unsuitable for interactive work. We believe that the proposed method for reducing the computation time will enable integration of the force density method and its iterative application in an interactive environment for form finding, where the designer can change boundary conditions and force values in cables, as well as mark the cables whose segment lengths must not be changed.
Extensive numerical experiments show that the proposed method is almost always efficient and robust, although there are cases in which the efficiency strongly depends on constants in the proposed termination rule that keeps the accuracy of linear solutions high enough not to compromise the convergence of the iterative process, but prevents it from to quickly becoming unnecessarily high. This rule is not the only possibility, which leaves ample room for further research.