An overview of the natural force density method and its implementation on an eflcient parametric computational framework**

The new paradigms of parametric modelling have been proving promising on the advance of systems for analysis and design of taut (or tensile) structures. With this premise, the presented work consist on the development with a form-finding tool for Computer Aided Design(CAE) and Computer Aided Engineering (CAE) integration using VPL (Visual Programming Language), in the context of parametricmodelling. Themethods used in the implementation are the Force Density Method (FDM) and the Natural Force Density Method (NFDM), taking advantage of the linear solution approach provided, suitable for fast form-finding computational procedures. The program is implemented as a Grasshopper plug-in and it is named BATS (Basic Analysis of Taut Structures), which enables parametric definition of boundary conditions for the form-finding. The program structure and benchmarks with other available Grasshopper plug-ins for taut structures form-finding are presented, showing considerably superior performance using BATS.


Introduction
Taut structures are defined by their characteristic state, in which the structural members work only in tension. They are usually flexible systems and are subjected to large displacements upon changes of the intensity and distribution of the loads.
Since Frei Otto's pioneering works in the late 1950s [1], they became an important research field for many architects and engineers. They provide lightweight structures, composed mainly of cables and membranes, which have no bending stiffness, and thus neither compression stiffness. Generally, it is not possible to define the shape of a taut structure a priori, but an appropriate shape is sought, compatible to an initial stress field imposed to the system. Once a viable configuration is found, the structural behavior under external loads can be determined considering geometrically non-linear analysis.
The Force Density Method (FDM) [2] provides the most convenient alternative for shape-finding, approximating a continuous surface by a network of linear elements. By its turn, the Natural Force Density Method (NFDM) [3] is an extension of the FDM that preserves the linearity of the original method, and overcomes some of its limitations to cope with irregular meshes which may arise from non mapped tessellations of surface geometries.
The imposition of natural force densities and the imposition of 2nd Piola-Kirchoff stresses on a reference configuration are equivalent [4], a property previously described for the original force densities by [5]. [4] also recognized that if the solution of the NFDM is recursively used as new reference configurations, the process might converge in few iterations to a configuration in which the resultant Cauchy stresses equal to the imposed 2nd Piola-Kirchhoff stresses. That means that for an isotropic homogeneous stress state, the final converged viable configuration will be a minimal surface. That approach is computationally more efficient for finding minimal shapes than other methods, such as pseudo-dynamic systems and geometrical optimizations. This makes NFDM well suited for structural design soft-wares, specially when the use of a parametric workflow is featured, as the user can quickly search for shapes, for instance changing the initial stress and then observing the final shape in real-time.A comparison between different methods on Ix-Cube 4.0 software, which has the NFDM implemented, was performed by [6] showing faster and more accurate results provided by NFDM. Although it is not studied in detail in this paper, the FDM and NFDM can also be applied for the form-finding of compression and mixed tension and compression structures (see [7]).
The new paradigms of parametric modelling have been showing promising on the advance of systems for analysis and design of tensile structures. The possibility to integrate Computer Aided Design (CAD) and Computer Aided Engineering (CAE) systems inside a parametric workflow allows better rationalization of the design process and automation of tedious tasks by linking CAD/CAE data for conception, analysis and detailing of structures. Finite Element Analysis (FEA) packages inside the visual programming language (VPL) Grasshopper have been showing outstanding results for structural design, as can be seen at [8], [9] and [10].
In this work, BATS (Basic Analysis of Taut Structures) is presented. BATS is a Grasshopper add-on which applies procedures for form-finding of taut structures and funicular shells inside parametric environment, and it was firstly prototyped for studies regarding the use of NFDM for funicular shell form-finding by [7]. The main code of BATS was updated from a pure C# code to a multi-language code with C# code for pre and post-processing using Grasshopper and Rhino API, which calls the solution from C++ functions, optimized for computational speed. This computational approach overcomes some issues presented on the previous code, as numerical inefficiency, mainly due to the use of the Open Source C# Library Math.NET, and non-optimized assembly of the FEA model from raw CAD geometry. The new code permits updates between many viable shapes in terms of milliseconds with an acceptable mesh refinement.
Benchmarks considering computational speed is presented considering saddle and catenoid minimal surface form-finding. BATS is compared to SATS (System for Analysis of Taut Structures, the first implementation of NFDM in MATLAB) [11], and Kangaroo solver [12], which is a multiphysics simulation Grasshopper add-on that uses pseudodynamic methods for many shape finding procedures, including minimal surfaces, with good computational performance.

Force Density Method
One of the first alternatives for form-finding was the definition of force densities, proposed by Linkwitz [2] and Sheck [13], in the context of cable nets.
The FDM is based on the equilibrium of each node in a cable net. With reference of the forces at the system in Figure 1, the resultant of internal forces acting on node i is Where N ij is the interaction force among the nodes i and ⃦ is the unit vector oriented from i to j. Applying equilibrium conditions gives equation 2 Relating nodes with nodal displacements results in a nonlinear equation system. However, defining at each element a force density n ij (equation 3), a linear system with 3n equations is obtained (4), and with boundary conditions imposed, the system can be easily solved.
Still, it is convenient for large meshes to adopt a matrix notation instead of the presented one. The Cartesian coordinates of the system, as the external and internal forces, can be expressed by 3 global vectors X, F and P, respectively  [3] defined by equations 5 are local coordinate and force vectors of each node. It is also convenient to substitute the sum inside the definition of P from the n nodes to b elements connecting those nodes. The vector of nodal coordinates and the vector of internal forces of the the element, as show in Figure 2, are related by equation 6.
x e = A e x ; P e = A eT p e Where as the director cosine of the element inside the global Cartesian coordinate system. A e is defined as the Boolean incidence matrix of the element.  [3] Once these definitions are given, the equation 4 can be rewritten as equation 7, where K d = ∑︀ e b=1 A eT k e d A e is the force density stiffness matrix of the system, with k e d defined in equation 8 is the element local stiffness matrix, with n e as the element force density and I 3 as an order 3 identity matrix.

Natural Force Density Method
The NFDM preserves the linearity of the original FDM meanwhile issues related to irregular triangular meshes are over-  [3] come. The natural force density derives from the natural forces defined at the natural membrane finite element, first proposed by Argyris [14]. Pauletti [3] redefines the formulation of previous developments into a more concise notation. The nodes and sides at Figure 3 are numbered anticlockwise, where each node index is coincident with the face index in front of it. Node coordinates are referenced at both global and local Cartesian system, with the local representation expressed by " ∧ ".
To define the internal vector forces, it is convenient to define natural forces N i , which are parallel to its faces. The relationship between the internal force vector of the element p e and the internal forces are defined in equation 9. After some algebra (see [15]) it is possible to rewrite the natural loads vector N = [N1N2N3] T as in equation 10, where V = At is the element volume, L = diag {l 1 , l 2 , l 3 } a diagonal "lenght matrix" and T is a transformation matrix, given by equation 11.
Comparing equations 9 with 6 suggest the definition of natural force densities n i = N i /l i , and using the relation given at equation 10 gives us a vector of natural force densities (equation 12).
With the definition of n it is possible to assign the natural force density stiffness matrix k nd as in equation 13, and the solution of a system with n linear elements and m triangular elements is given by equation 7 with K d = K (equation Although this linear procedure produces a viable configuration of a bending-free shapes, the imposition of natural force densities to achieve a uniform stress field is nonviable. However, the NFDM process can be evaluated recursively by assigning a constant stress fieldσ 0 and repeatedly determining K, usingσ =σ 0 on Eq. 12, where the new reference configuration is updated every step. [4] shows that this process converges the solution to a minimal surface, as the Cauchy stress resultants for this converged solution are equal to the uniform second Piola-Kirchhoff stress resultants, which is an analogous solution to a mathematical minimal surface by the soap bubble analogy. This process finds minimal surface shapes in quite few linear steps, as will be discussed further, and is an interesting approach as requires low computational cost and produces a viable configuration at every step, feature that have clear advantages compared to non-linear methods that converges to a minimal solution but usually through a series of nonequilibrium configurations.
The NFDM can find also non-minimal shapes by imposition of non-isotropic and non-uniform stress fields. A obvious example is application of a single step of the NFDM, which gives a viable shape, yet have a non-uniform stress state. The imposition of non-isotropic solutions can be done Figure 4: Definition of an orthotropic initial stress field by using director places [15]. by finding natural force densities by an initial orthotropic stress field using director planes on the elements, as detailed by [15]. [15] also proposed a method to extend the NFDM to quadrangular elements, by coupling multiple triangular elements on both directions to reduce biased stress results.

CAD/CAE integration by parametric workflow
The steps of design and analysis of taut structures presents some computational challenges, specially in the communication between CAD (Computer Aided Design) and CAE (Computer Aided Engineering) systems. CAD generally have resources for geometrical definition of complex shapes and details, where CAE is used for computational physical analysis. [16] defines the main issues on CAD/CAE integration, and address the main issues to loss of data, compatibility during the process and lack of automation. Those issues are due the different characteristics between the processes. A CAD model is mainly a computational representation of the geometry, not necessarily having attributes and properties. These features are crucial inside CAE environment, as it needs data such as material properties, physical interfaces and element types. Many CAE softwares offers CAD modelling features inside its interfaces, but still very limited compared to CAD specialized softwares. Also, some FEA packages are only code implementation, which need a external pre and post processor to work properly as a design tool.
On the other hand, specialized CAD software as Rhinoceros3D [17] have been delivering scripting features with API's (Application Programming Interface) giving the developer access to the CAD system geometrical objects inside an object-oriented-programming environment. With this features it is possible to develop CAE structural applications with full use of CAD capabilities, using geometrical objects as parameters for structural classes.
In CAD scripting context, mainly two types of programming are featured: Visual Programming Languages (VPL) and Textual Programming Languages (TPL). A VPL consists on visual block elements which contains algorithms and can be manipulated in a logical sequence of inputs and outputs. On other hand, TPL systems relies on a sequence of linear characters which describes the commands the program should execute [18]. VPL's are advantageous as do not require a broad knowledge on programming syntax, and are easier to use for general CAD users. The possibility of association of inputs and outputs in VPL's are interesting for CAE, as association between geometrical and analysis data can be made. However, VPL's are strictly dependent on TPL's, as each block component on the VPL is produced from a TPL code.
In terms of TPLs, applications that depends on numerical methods requires, besides code with low algorithmic complexity, a programming interface that offers nice memory allocation tools. Classical examples for high performance languages for linear algebra procedures are FOR-TRAN( with LAPACK and BIAS libraries) and C++ (with Eigen numerical library). On other hand, CAD software normally features APIs from more modern programming languages as C# and Phyton, as they are more versatile to work and delivers better features for object-orientedprogramming and dynamic types, which fits well for API purposes as the developer can access and correlate all geometric classes and its methods given from the CAD software.
Working with multiples languages and types of programming can then lead towards a nice computer performance with user-friendly interface, that can be established in 3 levels: (1) User level, or the VPL itself; (2) Application level, for communication between the user and the numerical level; (3) Numerical processing level, which solves the problem desired with optimized code.
Rhinoceros is a CAD software and provides it's VPL Grasshopper, which have components for manipulation of many Rhino geometries and so forth enables a parametric workflow inside it's environment. The user can define a logical sequence of events describing the project in function of predefined parameters, and instantly retrieve geometry and analysis feedback within changes on them. Both Rhinoceros and Grasshopper have their own C# API's (Rhinocommon (RC) and GrasshopperSDK (GS)), Where GS depends on RC, and allow custom components and parameters development for many purposes. Many CAE extensions were developed as Karamba3D [8], Kiwi3D [9] and Beaver [19] for structural analysis, Ladybug [20] for Thermar/Solar Analysis, and Butterfly [20] and Eddy3d [21] for CFD analysis.
The use of a parametric workflow in FEA is then a promising candidate for solving CAD/CAE intercommunication problems, as it can automate processes converting geometric information and other relevant data for analysis. The conversion of data is automated by assembling series of algorithms, where any change in geometry returns a suitable new structural model, avoiding rework in mesh generation and other possible tedious tasks.
Thus, structural analysis software within a parametric environment enables a user-friendly yet powerful design and analysis process. The environment allows the user to implement routines for specific and interconnected tasks based on the geometry and properties of the membrane. In addition, The VPL routines pattern allows the user to perform the analysis and extract and reuse the results for other routines such cutting patterns, connection details and verification of structural components by technical standards [22].

BATS implementation
In order to implement an efficient and robust form-finding method inside parametric workflow environment, the Grasshopper plugin BATS was developed. It was firstly prototyped in [7], and was not concerned yet about computer efficiency. However, it's resulting code appointed for two main issues regarding computational performance.
The first issue is related to the linear system solution and other linear algebra processes required. Matrices can be highly memory consuming, considering that given n nodes, corresponding each one to 3 degrees of freedom, the number of variables inside the matrix is on the order of (3n) 2 . This means that the complexity solving linear systems are in exponential order. Also, modern languages as C# and Python works with automatic memory management and dynamic types, which can be undesired for numerical purposes evolving matrices. In addition, it is well know in numerical analysis area that languages like FORTRAN and C++ provides tools for memory management which can optimize computational cost, and this naturally leads to more efficient numerical libraries available. In C#, a wideused numerical library is the Open-Source library Math.NET. Although quite user-friendly and versatile, it shows high computational cost solving large linear system as also other linear algebra procedures.
On other hand, the C++ open-source library Eigen have proven highly efficient, leading to fast results even for large systems. Benchmarks shows that Eigen performs almost the same as Math Kernel Library (MKL), a commercial nu-merical library well know by it's good performance [23]. That makes use of these resources highly desirable inside parametric workflow, as this can lead to practically realtime visualization of form-finding results. The advantage of C++ in Grasshopper development relies on great interoperability between C# and C++, where C++ functions in numerical level can be called by C# code.
The second issue is due the assembly of a FEA model based on geometrical information, where the program should index all geometrical nodes and assign the connectives between the elements. This is not a trivial task, as requires searching algorithms for equal points in 3D space. Linear search, which is the most obvious approach and was used in prototyping, provides a O(n) complexity, and can be highly consuming when dealing with thousands of nodes. Other search algorithms can be performed in order to improve the complexity of the algorithm. Binary search and also R-Tree search can deal with the problem respectively with O(log n) and O(log M n), where M denotes to the maximum entries in R-Tree nodes [24]. This task can be done in application level, as do not have matrices or linear algebra applications involved. Yet, it can still be a tedious process, taking more computational time than the solver, as evolves a search problem on a large list of 3D points.
The VPL components were developed in a C# application that uses C++ unmanaged functions from the developed numerical solver. Figure 7 shows the general workflow of the developed tool. First, raw Rhino/Grasshopper data is converted into structural analysis data, which is then used to assemble the structural model. All this data is structured in object-oriented-programming in C#.  The next step is to solve the model, where the C# application calls the solver function. In this step the structured data is transformed in a series of arrays, as the C# and C++ interoperability is straightforward only for native C types. In the C++ code, the data is then restructured in objectoriented style, but it have different purposes than the C# data. The C# application data aims at general structural data, that any FEA user would be familiar with, and the C++ solver data focus on attributes and methods for assembly of the external force vector and stiffness matrix. Although firstly required for simplicity in data transfer between two programming languages, the scheme of have independent data structures for the application and the solver showed to be of benefit, as each implementation can deal with their own needs separately.
Latter, the C++ code outputs the results also in arrays, and the C# application retrieves the results and prepare them for visualization inside the Rhino/Grasshopper interface. The details of each part of the developed tool, respectively the VPL components, the C# application code and the C++ numerical code, are discussed in detail the following subsections.  Figure 7. A series of components retrieves geometrical data of meshes, lines and points and it's respective supplementary data, for force density linear elements, natural force densities mesh elements, supports and loads.

VPL components
Each component outputs custom parameters which contain all data provided, and are used to assemble the model, subsequently sending the model parameter to the solver component, which updates the model object with the shape found and its respective stress field. Results can be retrieved by special components, where membrane view, cable view and support reaction view are available.
With that workflow, the user can modify the formfinding parameters (i.e. initial stress fields) as well as the initial geometry and its constraints. That makes the process of form-finding very versatile, as shape finding definition can be defined by geometrical and structural design parameters inside a CAD environment and with real-time updates within parameter changes. Figure 9 shows the class diagram of the application code in C#, that is the core for the development of the VPL components. It is responsible to intercommunicate Rhino/Grasshopper with BATS numerical engine, and contains a collection of classes with geometric variables defined as Rhino geometrical objects in RhinoCommon, with supplementary information relevant for the analysis.

Application level
Node class collects variables regarding it's 3D geometrical point in RhinoCommon and have it's FEA index associated, which is determined on a further step inside Model class. For post-processing and colour plot purposes, it also have methods and fields for resultant nodal stress definition. Support Class contains a Node class field, which also have only geometrical information before the preprocessing. More fields are needed for these objects, as support plane, prescribed displacements and constrained DOFs. Reactions field and methods are needed in order to post-process the model after solution.  Class contain the geometrical line defined by two non-indexed nodes (i.e. geometrically informed only) which represents the element and other information described in Section 2. Also, a material density field is available for funicular shell purposes.
NFD classes contains a single triangular or quadrangular element defined by it's non-indexed nodes and other relevant information described in Section 3. Similar to FD class, it also contains fields regarding material density. For simplicity and assembly process, a auxiliary NFDMesh class is used, as this permits the definition of properties in a whole mesh as also speeds up model assembly as connectivity data of the mesh can be directly passed through the model. Elements have also fields for post-processing results, as the resultant stress fields.
Model class contains list field of all classes listed above, and has the relevant methods for pre and post-processing. It also contains the C++ calling method for the solve process. The pre-process is done by calling NodeIndex method, which retrieves all tridimensional geometrical information and provides index assignment of the global model nodes as also the reference global supports and element nodes indexes. BuildtoC and SolveonC methods grabs all structured data and assign series of data arrays to be sent to the C++ function, which will handle the numerical processing and retrieves the NFDM results. At last, post-process on each object type is made by PostProcessElements, PostPro-cessNodes and PostProcessSupports methods.
Every application class have it's corresponding parameter wrapper inside Grasshopper, which makes possible connection of objects data between components. That makes possible the assembly structure proposed in the VPL components, as the VPL do not recognize automatically custom classes that are not on RhinoCommon or GrasshopperSDK.

Numerical level
Numerical level refers to series of C++ functions related to a main function, which is called from application level. The main function must call it's variables considering interoperability between both C# and C++ variable types. C# provides compatibility of main C native data types when calling C++ functions using DllImport resources [25]. Arrays from the native types, that are needed due to size of data, can be passed considering them as pointers in C++, that will access directly the memory data called in C# code.
Sparse matrices are used as they storage only non-null values, considerably reducing memory allocation. However, sparse matrix build-up can be slow when applied changes on global indexes, as its structure is not directly correlated to dense matrix indexes [26]. Hopefully, sparce matrix assembly can be optimized by directly providing a set of indices and corresponding values, which are called triplets, defined by three values (i, j, value).
For matrix assembly, the global stiffness matrix K is partitioned between free and constrained DOFs, resulting on four sub-matrices K , K fc , K cf , Kcc, as the solution process requires only the unconstrained part, which can cause economy in highly constrained systems, and reaction forces can be directly obtained using the results with the other submatrices. Therefore, a list of triplets is needed for each submatrix, which is done by the Partition method, which stores each triplet list for further matrix assembly. Figure 10 shows auxiliary classes for numerical model assembly. Arrays data are rebuilt in a scheme of classes DOFs → Nodes → Elements. This helps to organize matrix assembly of the linear system in a object-oriented scheme. Elements have information about force densities and also of its nodes, which defines its type (linear, triangular or quadrangular). Node objects contains geometrical position, its nodal index and correponding DOFs. DOF object then collects unpartitioned and partitioned indices. That approach turns matrix assembly simple yet effective, as it avoid further searching algorithms for constrained DOFs during assembly, as they are directly considered before this step. The same partitioning process is done in the force vector F in F f and Fc. LL T factorization with Cholesky Decomposition is used to solve linear systems, as K preserves the symmetry and positive definitive properties required. At last, outputs delivers the new geometry with the resultant stress field and reactions.
Addressing these main issues retrieves good performance even for large systems,as is shown in Section 6.1.

Examples and benchmarks
Examples and benchmarks are shown in order to evaluate BATS reliability and computational performance. All tests were made with a personal computer with the following specifications: Intel i7 9750h processor (2.4Ghz), 16GB RAM. BATS is compared for minimal surface form-finding with two other solvers: Kangaroo [12] and SATS [11].
Kangaroo is a multi-physics solver with a wide-range of physical simulations and features coupling between many physical objectives. In terms of membrane and cable design, it offer both minimal edge lengths and minimal surface solutions, where the latter is obtained by an iterative mathematical procedure imposing zero Gaussian's curvature on the mesh.
SATS is a MATLAB program which contains routines for form-finding and analysis of taut-structures and was the first program which NFDM was implemented. As BATS is a further development of SATS with better user-interface, only computer performance is analysed.
For all benchmarks it is assumed that the 200 th iteration of NFDM is an acceptable minimal surface. [6] and [15] shows that for the shapes analysed the NFDM converges to these solutions in few steps. Then the geometrical error for each procedure can be evaluated and analysed at each iteration by the following relation: Figure 11 shows an initial mesh made from a simple NURBS patch. The model is defined by a square mesh of NFDM elements with FD elements in the mesh boundary, supported by the 4 corner edges. The addition of FD elements in the boundary adds a restriction to the minimal shape finding procedure which gives the best-fit solution for a minimal shape considering the tension in the FD elements. For the models assembled in Kangaroo, both considers an additional boundary restriction with linear spring elements at the border. This restriction is important as the surface degenerates if only the corner points are set fixed. More details on this matter can be seen at [27]. Performance analysis was made in three different mesh densities: 10x10, 25x25 and 50x50 grid. On all tests, the shape found for the saddle surface matches for the NFDM and minimal surface procedure with boundary restrictions on Kangaroo. However, the minimal edges procedure on Kangaroo presented a deviation from the exact solution. This is explained by the fact that the minimal edges procedure considers linear elements that are dependent on line length and mesh topology, and the correct definition for force densities is non intuitive. The definition for minimal edges on each line addresses the direction of principal stresses to the lines axis, and this constraint leads to the error observed. Therefore, the minimal edges length method, which is analogous to the original force density method, does not provides a correct solution.

Minimal saddle surface with boundary restriction
The following graphs ( Figure 12, 13 and 14) shows the evolution of the geometrical error as a function of both iteration number n and the time elapsed t.
With exception of the 10x10 grid solution, the mathematical minimal surface procedure on Kangaroo presented initial increase of the error, which starts to converge when the other methods are already or almost on convergence domain. It can be seen that in both terms of n and t BATS presented direct convergence to the solution, requiring a maximum of 20 iterations to converge (in the 50x50 grid case).
Both methods in Kangaroo presented an oscillatory behaviour for the error, and required more iterations to con- verge. Even that the minimal edge surface method required fewer time per iteration, the NFDM property of directly find viable shapes at each iteration, makes it more efficient as quite few steps are required to achieve convergence.
Interactive pseudo-dynamic methods presented on Kangaroo are explicit methods and so forth requires more steps to reach a solution, as they are dependent on the time step at each iteration. The procedure tries to find a equilibrium state in the pseudo-dynamic scheme by imposing constraints on the motion, which error function in  terms of time/iteration presents the behaviour of damped oscillations. However, these methods can be helpfull when dealing with non-linear structural analysis of tensile nonlinear materials, as the problem can be derived only from the internal and external force vectors and avoids the use of the tangent stiffness matrix, that might be singular due to the local instabilities on the membrane. This method for non-linear systems solution is the so called Dynamic Relaxation method, that, although lacks computational efficiency compared to direct methods, it is more robust and achieve solutions which are quite difficult to using a Newton-Raphson solver. Further reading in this topic is addressed to [28] and [29]. Figure 15 gives the time elapsed on each procedure for convergence. BATS performs considerably quickly than other implementations. SATS elapsed time is presented and requires the same iterations than BATS, as also uses the NFDM procedure, however presented lower performance. The main differences on the implementation of BATS and SATS is the language (C++ and MATLAB), numerical library (Eigen and MATLAB), sparse matrix assembly (optimized with triplets and no optimization) and linear system solving (factorization with Cholesky decomposition and generalized linear system solver).

Minimal catenoid surface
This example is defined from an initial mesh of a cylinder of radius r and height h, for the form-finding of a minimal catenoid shape. For this class of shapes, a minimal surface can be found analytically with the Goldschmidt solution, which states that a minimal surface of a catenoid, with two equal radius rings as boundary, only have a solution for h < 1.3254868r [30].
As explained before, the minimal edges procedure is not equivalent to a minimal surface form-finding and will be neglected in this example. Three different heights are defined in order to check the Goldschmidt limit: h = 1.30r (Figure 16), h = 1.32r ( Figure 17) and h = 1.34r (Figure 18). Besides a apparent gradient of stresses showed on the stress colour plot in the catenoid, the maximum stress variation is on order of 10 − 4, and then it can be clearly assumed that the stress distribution is uniform for all converged BATS results.
Expected results arises for h bellow the Goldschmidt limit, where both procedures converges to the same shape.Above the limit, both solutions diverge into a degenerative mesh as no solution is viable. However, it can be seen that for h = 1.32r, which is quite close to the limit, BATS converges to a solution whereas Kangaroo don't, showing better sensitivity of results using the NDFM.
The graph of geometrical error in function of time elapsed is presented on Figures 19 and 20. It can be observed that for h = 1.30r both solutions behaves in the     same manner as the saddle surface example. On other hand, h = 1.32r shows Kangaroo almost reaching a solution but diverging the minimum error is achieved, which happens at both BATS and Kangaroo solution for h = 1.34r.
However, as NFDM produces at each iteration an viable solution, it's possible to retrieve the last iteration in which the stress field deviance decreases. For h = 1.34r case, Figure 21 shows the maximum and minimum stress acting on the surface, where when the values are equal a uniform stress field is achieved. It was observed that the shape with  the stress field closest possible to a uniform state is found in the 20 th iteration, and it can be seen on Figure 22.
The error between each iteration and this solution is plotted in Figure 23, where the convergence domain and divergence domain of the error function can be observed.

Conclusion
This work showed a implementation of the Natural Force Density Method in a novel framework composed by Grasshopper components, application code in C# for component programming, and a numerical code in C++.
Comparison with the first implementation showed the importance of efficient memory management and the choice of numerical libraries and methods. Comparison with other form-finding methods in Kangaroo shows how an efficient implemented implicit method can perform highly better than a nicely implemented explicit method. Besides it, the developed tool delivers viable configurations at all iterations as also deliver the resultant stress field associated with the shape.
The combination of a robust form-finding method as NFDM with adequate numerical programming techniques and efficient numerical libraries inside a parametric environment produced a reliable, fast and versatile formfinding tool for the shape finding of taut structures and funicular shells. Future research is on applying the computational framework for non-linear analysis and cutting pattern routines for taut structures, taking advantage of the numerical efficiency presented.

Funding information:
The authors state no funding involved.
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.

Conflict of interest:
The authors state no conflict of interest.