Fast Stiffness Matrix Calculation for Nonlinear Finite Element Method

We propose a fast stiffness matrix calculation technique for nonlinear finite element method (FEM). Nonlinear stiffness matrices are constructed usingGreen-Lagrange strains, which are derived from infinitesimal strains by adding the nonlinear terms discarded from small deformations. We implemented a linear and a nonlinear finite element method with the same material properties to examine the differences between them.We verified our nonlinear formulationwith different applications and achieved considerable speedups in solving the system of equations using our nonlinear FEM compared to a state-of-the-art nonlinear FEM.


Introduction
Mesh deformations have widespread usage areas, such as computer games, computer animations, fluid flow, heat transfer, surgical simulations, cloth simulations, and crash test simulations.The major goal in mesh deformations is to establish a good balance between the accuracy of the simulation and the computational cost; achieving this balance depends on the application.The speed of the simulation is far more important than the accuracy in computer games.The simulation needs to be real-time in order to be used in games so free-form deformation or fast linear FEM solvers can be used.However, high computation cost gives much more accurate results when we are working with life-critical applications such as car crash tests, surgical simulators, and concrete analysis of buildings; even linear FEM solvers are not adequate enough for these types of applications in terms of the accuracy.
For realistic and highly accurate deformations, one can use the finite element method (FEM), a numerical technique to find approximate solutions to engineering and mathematical physics problems.FEM could be used to solve problems in areas such as structural analysis, heat transfer, fluid flow, mass transport, and electromagnetics [1,2].
We propose a fast stiffness matrix calculation technique for nonlinear FEM.We derive nonlinear stiffness matrices using Green-Lagrange strains, themselves derived from infinitesimal strains by adding the nonlinear terms discarded from infinitesimal strain theory.
We mainly focus on the construction of the stiffness matrices because change in material parameters and change in boundary conditions can be directly represented and applied without choosing a proper FEM [3].Joldes et al. [4] and Taylor et al. [5] achieve real-time computations of soft tissue deformations for nonlinear FEM using GPUs; however, they do not describe how they compute stiffness matrices; thus, we cannot implement their methods and compare them with our proposed method.Cerrolaza and Osorio describe a simple and efficient method to reduce the integration time of nonlinear FEM for dynamic problems using hexahedral 8-noded finite elements [6].We compare our stiffness matrix calculations with Pedersen's method [3] to measure performance and verify correctness.We achieve a 142% speedup in calculating the stiffness matrices and a 15%

The Nonlinear FEM with Green-Lagrange Strains
We use tetrahedral elements for modeling meshes in the experiments.Overall, there are 12 unknown nodal displacements in a tetrahedral element.They are given by [2] {} = { (, , In global coordinates, we represent displacements by linear function by For all 4 vertices, ( 2) is extended as Constants   can be found as where V −1 is given by det(V) is 6, where  is the volume of the tetrahedron.If we substitute (4) into (2), we obtain Because of the differentials in strain calculation,  is not used in the following stages.If we expand (6), we obtain For tetrahedral elements, to express displacements in simpler form, shape functions are introduced ( 1 ,  2 ,  3 ,  4 ).They are given by In our method, nonlinear stiffness matrices are derived using Green-Lagrange strains (large deformations), which themselves are derived directly from infinitesimal strains (small deformations), by adding the nonlinear terms discarded in infinitesimal strain theory.The proposed nonlinear FEM uses the linear FEM framework but it does not require the explicit use of weight functions and differential equations.Hence, numerical integration is not needed for the solution of the proposed nonlinear FEM.Instead of using weight functions and integrals, we use displacement gradients and strains to make the elemental stiffness matrices space-independent in order to discard the integral.We extend the linear FEM to the nonlinear FEM by extending the linear strains to the Green-Lagrange strains.We constructed our linear FEM by extending Logan's 2D linear FEM to 3D [2].To understand Green-Lagrange strains, we must first see how they differ from the infinitesimal strains used to calculate the global stiffness matrices in a linear FEM. Figure 1 shows a 2D element before and after deformation, where the element edge  with initial length  becomes     .The engineering normal strain is calculated as the change in the length of the line The final length of the elemental edge can be calculated using By neglecting the higher-order terms in (11), 2D infinitesimal strains are defined by By the definition, the nonlinear FEM differs from the linear FEM because of the nonlinearity that arises from the higher-order term neglected in calculation of strains.The strain vector used in the linear FEM relies on the assumption that the displacements at the -axis, -axis, and -axis are very small.The initial and final positions of a given particle are practically the same; thus, the higher-order terms are neglected [7].When the displacements are large, however, this is no longer the case and one must distinguish between the initial and final coordinates of the particles; thus the higher-order terms are added into the strain equations.By adding these high-order terms, 3D strains are defined as [8] which leads to The Green-Lagrange strain tensor is represented in matrix notation as where {} is the nodal displacement, [  ] is the linear, and [ NL ] is the nonlinear part of the [ 0 ] matrix [3].For a specific element, [  ] and [ NL ] are constant, as with the [] matrix in the linear FEM.With the variation of {} [9], (15) becomes Gathering the strain components together, we can rewrite ( 15) and ( 16) as The linear part of the ) is the same as the [] matrix in the linear FEM.Calculating [ 0 ] becomes more complex with the introduction of the nonlinear terms.After finding the nonlinear strains, these equations are combined with the shape functions to find matrix [ 0 ]: The most frequently used terms, which are the nine displacement gradients for calculating the nonlinear strains, are   /,   /,   /,   /,   /,   /,   /,   /, and   /.Using ( 8) for displacements, we can construct the displacement gradients using the partial derivatives of the shape functions.They are represented by where   represents   /.
We can evaluate the partial derivatives of the shape functions as follows (for the 1st node of [ NL ]): Using ( 17) and (20), we obtain [ 0 ] for the 1st node (21).Similarly, using ( 17) and (20), we obtain [ 0 ] for the 1st node (22). Consider The FEM is derived from conservation of the potential energy, which is defined by where E strain is the strain energy of the linear element and  is the work potential.They are given by where the engineering strain vector {} is From (24), the engineering stress vector  is related to the strain vector by The secant relations are described by the matrix [].We substitute (15) and ( 26) into (24), obtaining the element stiffness matrix We can discard the integrals as we did for the linear FEM.
[ 0 ], [], and [ 0 ] are constant for the four-node tetrahedral element, so (27) is rewritten as The secant stiffness matrix which is

Introducing nodal forces, we obtain
With the equilibrium equation and cancelling {}  , the whole system for one element reduces to By substituting {} with , we obtain Finally, only nonlinear displacement functions remain, which are solved with Newton-Raphson to find the unknown displacements  [10].Element residuals are necessary for the iterative Newton-Raphson method.The element residual is a 12 × 1 vector for a specific element.The residual for a specific element is defined as Having determined   , we can now express (32) in expanded vector form as () (3,1) +   () (3,2) +   () (3,3) + ⋅ ⋅ ⋅ +   () (3,12) . . .
The tangent stiffness matrix []   (  ) is also necessary for the iterative Newton-Raphson method.The tangent stiffness matrix is also 12 × 12 matrix, like the elemental stiffness matrix.However, the tangent stiffness matrix depends on residuals, unlike the elemental stiffness matrix.Elemental stiffness matrices are used to construct residuals and the derivatives of the residuals are used to construct the elemental tangent stiffness matrices.We can express the elemental tangent stiffness matrix for a specific element as Newton-Raphson method is a fast and popular numerical method for solving nonlinear equations [10], as compared to the other methods, such as direct iteration.In principle, the method works by applying the following two steps (cf.Algorithm 1): (i) check if the equilibrium is reached within the desired accuracy; (ii) if not, make a suitable adjustment to the state of the deformation [11].An initial guess for displacements is needed to start the iterations.The displacements are updated according to In the proposed nonlinear FEM,  is the vector that keeps the information of the nodal displacements.Instead of making only one assumption, we make whole  vector initial guess in order to start the iteration. Consider where  is residual of the global stiffness matrix [] calculated in (33) and   is the tangent stiffness matrix calculated in (34).
At every step, the vector  and the matrix   are updated for every element with the new   values.Then,  and   are assembled as we did with for the global stiffness matrix  and the global force vector  in linear FEM.Boundary conditions are applied to the global  vector and the global   matrix.Using the global  vector and the global   matrix, we have is updated with the solution of (37).Consider Then, we check if the equilibrium is reached within the desired accuracy defined by  as      (  )     ≤ . (39) After the desired accuracy is reached, the unknown nodal displacements are found.

Experimental Results
The proposed nonlinear FEM and Pedersen's nonlinear FEM were implemented using MATLAB programming language.The visualizer was implemented with C++ language and connected to the solver using the MATLAB engine [12], which allows users to call the MATLAB solver from C/C++ or Fortran programs.The simulation results, interaction with the 3D model, and the 3D models themselves were visualized using OpenGL, and the nose experiment was visualized using 3ds Max [13].To speed up the nonlinear FEM, we used MAPLE's symbolic solver [14], which is integrated into MATLAB.We conducted all the experiments on a desktop computer with a Core i7 3930 K processor overclocked at 4.2 GHz with 32 GB of RAM.We used linear material properties for the models in the experiments.We used 1 GPa for Young's modulus () because polypropylene has Young's modulus between 1.5 and 2 GPa and polyethylene HDPE has Young's modulus 0.8 GPa, which shows plastic properties and they are close to 1 GPa.Because most steels and plastic materials undergo plastic deformation near the value of 0.3, we used 0.25 for Poisson's ratio (]).Our simulation is static so we used a single load step in all experiments.Multiple load steps are used when the load forces are time-dependent or the simulation is dynamic [15].We conducted four experiments, each having different number of elements to observe the speedup for both stiffness matrix calculation and for the solution of the system.As expected, the proposed method and Pedersen's method produced same amount of nodal displacements in all experiments.The first experiment was conducted for a cube mesh with eight nodes and six tetrahedral elements.Figure 2 shows that the cube is constrained at the upper four nodes and pulled downwards with a small amount of force (one unit force for each of the upper four nodes).This experiment was conducted with a small mesh in order to carefully examine the nodal displacements and strains for each element.Figure 3 shows the initial and final positions of the nodes for the linear and nonlinear FEMs, respectively.As seen in Figure 3, the   linear and nonlinear methods produce similar displacements when the force magnitude is small.Tables 1 and 2 show the displacements at force applied nodes (first, second, third, and fourth) using the linear and nonlinear FEMs, respectively.Figure 4 shows that displacement increases linearly with force magnitude.However, as expected, the nonlinear FEM behaves quadratically due to the nonlinear strain definitions.Figure 5 depicts the convergence of the Newton-Raphson method for the nonlinear FEM.
The second experiment was conducted on a beam with 90 nodes and 216 tetrahedral elements.Figures 6(a   twisted at both ends.Figure 7 shows the final shape of the beam mesh for both the proposed and Pedersen's methods.Tables 3 and 4 show the displacements at force applied nodes (green nodes) for the second experiment using the linear and nonlinear FEMs, respectively.The third experiment was conducted with a cross mesh of 159 nodes and 244 tetrahedral elements.We aimed to observe if there is a root jump occurring when solving the system for a high amount of force (50 N units), and its effect Journal of Applied Mathematics  on the computation times for both methods.Figure 8 shows that the cross shape is constrained at the blue nodes and pushed towards the green nodes.Figure 9 shows the final shape of the beam mesh for both the proposed and Pedersen's methods.Tables 5 and 6 show the displacements at force applied nodes (green nodes) using the linear and nonlinear FEMs, respectively.We conducted fourth experiment with a liver mesh of 465 nodes and 1560 tetrahedral elements.Figure 10 shows that the mesh is constrained at the blue nodes and pulled from the green node (30 N units) in the direction of the arrow.We aimed to observe the similar amount of speedup like previous experiments for a high density mesh.Figure 11 shows the final shape of the beam mesh for both the proposed and Pedersen's methods.Tables 7 and 8 show the displacements at force applied node (node number 271) using the linear and nonlinear FEMs, respectively.
Computation times of the finite element experiments are required to compare how much faster our proposed method is than Pedersen's.When comparing nonlinear FEMs, we calculated the computation times to construct the stiffness matrices as well as the computation times of the nonlinear FEM solutions to determine how different calculations affect them.Table 9 depicts the computation times for the stiffness matrix calculation and Table 10 depicts the computation times for the system solution.Table 11 shows the iteration counts to solve the system using the Newton-Raphson procedure.
The speed-up columns of Tables 9 and 10 depict the speedups of the proposed method compared to Pedersen's method for the stiffness matrix calculation and the system solution using a single thread, respectively.The speedup is calculated as follows: Speedup = Runtime (Pedersen's method) Runtime (The proposed method) .Our proposed method outperforms Pedersen's method.On the average, it is 142% faster at computing stiffness matrices because Pedersen's method uses more symbolic terms.However, both methods use Newton-Raphson to solve nonlinear equations, which takes approximately 90% of the computation time.Thus, the overall speedup decreases to 15% on average.In experiment 3, because of more iterations due to root jumps, there was a performance loss against Pedersen's method.
We also applied our method for corrective operation on the misshapen nose of a head mesh.The head mesh is composed of 6709 nodes and 25722 tetrahedral elements (see Figure 12(a)); all the operations were performed in the nose area of only 1458 tetrahedral elements.Figure 12(b) shows that the head mesh is constrained at the blue nodes and pushed upwards at the green nodes and Figure 12(c) shows the result of the nonlinear FEM.

Conclusions and Future Work
We propose a new stiffness matrix calculation method for nonlinear FEM that is easier to analyze in terms of constructing elemental stiffness matrices and is faster than Pedersen's method.The proposed method is approximately 2.4 times faster, on average, at computing stiffness matrices and 15% faster at computing the whole system than Pedersen's method.

Figure 2 :
Figure 2: A 10 × 10 × 10 m cube mesh with eight nodes and six tetrahedra is constrained at the blue nodes and pulled downwards from the green nodes.

Figure 3 :
Figure 3: (a) The initial and final positions of the nodes for the linear FEM.(b) The initial and final positions of the nodes for the nonlinear FEM.

Figure 4 :
Figure 4: Force displacements (in m) at node 4 for the linear and nonlinear FEMs.

Figure 5 :
Figure 5: Newton-Raphson convergence graphics for the nonlinear FEM.The graph is plotted using the logarithmic scale.
) and6(b)   show that the beam is constrained at the blue nodes and

Figure 6 :Figure 7 :
Figure 6: The beam mesh is constrained at the blue nodes and twisted at the green nodes.(a) Front view; (b) side view, which shows the force directions applied on each green node.

Figure 8 :
Figure 8: The cross mesh is constrained at the blue nodes and pushed towards the green nodes.

Figure 9 :Figure 10 :
Figure 9: Nonlinear FEM results for the both proposed and Pedersen's methods: (a) initial and final wireframe meshes are overlaid; (b) initial and final shaded meshes are overlaid.

Figure 11 :Figure 12 :
Figure 11: Nonlinear FEM results for both the proposed and Pedersen's methods: (a) left: wireframe surface mesh; right: wireframe surface mesh with nodes; (b) left: shaded mesh; right: shaded mesh with nodes.

Table 1 :
The displacements (in m) at nodes 1, 2, 3, and 4 using the linear FEM for the first experiment.The displacements of nodes 5 to 8 for all axes are zero.

Table 2 :
The displacements (in m) at nodes 1, 2, 3, and 4 using the nonlinear FEM for the first experiment.The displacements of nodes 5 to 8 for all axes are zero.

Table 3 :
The displacements (in m) at green nodes using the linear FEM for the second experiment.The displacements at blue nodes are zero.

Table 4 :
The displacements (in m) at green nodes using the nonlinear FEM for the second experiment.The displacements at blue nodes are zero.

Table 5 :
The displacements (in m) at green nodes using the linear FEM for the third experiment.The displacements at blue nodes are zero.

Table 6 :
The displacements (in m) at green nodes using the nonlinear FEM for the third experiment.The displacements at blue nodes are zero.

Table 7 :
The displacements (in m) at green nodes using the linear FEM for the fourth experiment.The displacements at blue nodes are zero.

Table 8 :
The displacements (in m) at green nodes using the nonlinear FEM for the fourth experiment.The displacements at blue nodes are zero.

Table 9 :
Computation times (in seconds) of stiffness matrices (Pedersen: Pedersen's nonlinear FEM; Proposed: proposed nonlinear FEM; and Speed-up: relative performance comparison of the stiffness matrix calculation of the proposed nonlinear FEM method with Pedersen's nonlinear FEM method using single thread).