A Moore-Penrose continuation method based on a Schur complement approach for nonlinear finite element bifurcation problems

Continuation methods have proved to be very powerful tools when solving large ﬁnite element problems. However, implementation of these methods often require modiﬁcations to the standard ﬁnite element method. As a ﬁnite element code is already very complex, we would like to implement the continuation method as efﬁciently as possible. In this paper, we present a new implementation technique based on a Schur complement approach for the Moore–Penrose continuation method. This method facilitates the detection of bifurcation points and also enables branch following. Numerical examples will be presented and analyzed using the proposed approach.


Introduction
The finite element simulation of very large deformations of hyperelastic materials is still a challenging problem. These problems are generally driven by a loading parameter and it is often observed that for some values of this parameter, the solution varies extremely rapidly due to geometric and/or material non linearities, often leading to the break down of the solution process.
The typical solution strategy is often based on Newton like methods. The loading parameter is set to a given value, the Jacobian matrix stemming from the finite element discretization is constructed, the corresponding linear system is solved to correct the solution and this is repeated until convergence. The loading parameter is then somehow increased (''by hand'') and the process is repeated until the total load has been imposed on the structure. One crucial step in terms of computational cost is the solution of the linear systems involving the Jacobian matrix. For moderate size problems, an LU decomposition can be used but for very large problems, problem specific preconditioned iterative solvers must be designed and implemented.
In some situations, following the solution branch by increasing heuristically the loading parameter becomes extremely difficult and results in the divergence of the process. This is clearly the case in the neighborhood of limit points or bifurcation points where the Jacobian matrix becomes singular. Numerical continuation methods have proved to be a very powerful tool when dealing with these kinds of problems. They can be roughly divided into two main branches, predictor-corrector methods and piecewise linear methods (see [2]). Both types of methods share many common features and can be numerically implemented in similar ways. However, predictor-corrector methods generally perform best when high accuracy is needed and are thus frequently used in practice. Many different continuation methods exist (see for example [14,11,24,25]), but in this paper, only the Moore-Penrose (also known as Gauss-Newton) continuation method (see [2,13,12]), which is a predictor-corrector method, will be considered.
Nonlinear structural analysis leads invariably to collapse and buckling analysis [7,31]. Buckling is a mathematical instability which is characterized by a sudden failure of the structure when critical loads are reached. The value of these critical loads, which are associated with bifurcation points, are quantities of interests [29,21,30] and correspond to eigenvalues of the system. Detecting such points and following the different solution branches passed these critical loads are thus important and useful to better understand the physical properties of the problem we are solving. The analysis of structural instabilities, which includes the detection of bifurcation points, is generally based on a continuation method (see [11,6,14,18,27,1,22,28]). The implementation of the continuation method presented in this paper will therefore offer an efficient and simpler approach for these analyses.
In the Moore-Penrose method, the loading parameter becomes an additional unknown of the problem allowing a more precise control of the deformation, increasing its speed whenever possible and slowing it in some critical regions. This is done by modifying the linear system in the Newton method to take care of the additional unknown. This therefore requires a potentially important modification of an existing finite element code. In particular, new iterative solvers must be implemented to solve these increased linear systems.
As many of these problems have thousands or even millions of degrees of freedom, the continuation method should be implemented as efficiently as possible in order to reduce the already high computational cost. In order to achieve this goal, information that has already been calculated for the finite element method should be used whenever possible. This is precisely the main objective of this work. We will present a generalization of the Moore-Penrose continuation method that uses only information already available in a standard finite element code. In particular, all linear systems that have to be solved are based on the classical Jacobian matrix of the problem. Therefore, if an LU decomposition of this matrix is available, it can be used directly and if a specific iterative method was used, it can still be used without any modification to follow the solution curve, detect bifurcations and possibly follow various solution branches.
This paper thus presents a new implementation technique for the Moore-Penrose continuation method when applied in a finite element context. This new approach, which takes advantage of information already available, also facilitates the detection of bifurcation points. Section 2 briefly recalls the Moore-Penrose continuation method as well as the classical numerical implementation of this method. Section 3 is devoted to the details of the new implementation technique. Section 4 explains how the approach can be used to easily detect bifurcation points and finally, Section 5 is devoted to the validation of the proposed algorithm. In particular, a study of the classical elastic beam buckling problem will be presented and analyzed.

Moore-Penrose continuation method
The finite element discretization of nonlinear elasticity problems leads to nonlinear problems of the form: where u consists of the degrees of freedom. In general, u is the displacement field, but in the case of mixed formulations, u consists of the degrees of freedom for the displacement and the pressure. In most situations, a loading parameter k, corresponding to either external forces or prescribed displacements or both, implicitly drives the deformation. For highly nonlinear problems, applying the desired loading will frequently cause a breakdown of the numerical method used to solve this problem. Therefore it is standard practice to gradually increase the value of the parameter until the desired loading, k max , is attained. Typically, a Newton method is used to solve the nonlinear system. The following steps are therefore repeated until convergence is achieved: 1. Given an initial estimate u 0 and a fixed load k. 2. Solve the linear system: 3. Update the solution: where F 0 u is the Jacobian matrix of F and d k u is the correction vector. The loading parameter k is then increased by a certain quantity and the process is repeated until the total load k max has been imposed. The increment in k is performed more or less heuristically, depending on the convergence of the Newton method. This is highly inefficient and may result in very small increments or even divergence of the algorithm in the neighborhood of turning points or bifurcation points on the solution curve.
A more efficient approach is to use continuation methods where the loading parameter k is explicitly introduced in the system of nonlinear equations which is now expressed as: with F a smooth function of R Nþ1 into R N . The major difference is that the vector of unknowns, x, now consists of the displacement plus the loading parameter. Starting from a point x ðiÞ 2 R Nþ1 satisfying Fðx ðiÞ Þ ¼ 0, and given a vector v ðiÞ tangent to the solution curve at this point, the goal is to follow the curve up to the point where k ¼ k max . It is then natural to try to follow the solution curve by first making a prediction step of length h i in the tangential direction: As X 0 is not likely to be a solution of Eq. (1), the Newton method can be used to obtain the next point x ðiþ1Þ 2 R Nþ1 on the solution curve.
Assuming that the rows of AðX 0 Þ are linearly independent (i.e. AðX 0 Þ has full rank), which is generally the case since F 0 u ðX 0 Þ is the Jacobian matrix of the standard finite element method and is therefore invertible, solution of system (3) is given by: This therefore leads to the following generalization of the Newton method, which is also called the Moore-Penrose continuation method: For k ¼ 0; 1; 2; . . . ; k max 1. Calculate the Moore-Penrose correction: 2. Update the solution vector: 3. If kFðX k Þk 6 e F and kX kþ1 À X k k 6 e x , convergence attained: x ðiþ1Þ ¼ X kþ1 We note that k max represents the maximum number of iterations allowed while e F and e x are the desired tolerances on F and x respectively.
To obtain the correction d k x in Eq. (4), the Moore-Penrose pseu- The explicit computation of this matrix is extremely costly in realistic situations where the matrix A is large. It is however easily seen that d k x is also a solution of the following system: where L k is any vector in the kernel of the rectangular matrix AðX k Þ.
This can easily be verified by replacing d k x (i.e. Eq. (4)) in the system. The major difficulty that arises with this approach however is the construction of vector L k as it must be in the kernel of AðX k Þ. In an iterative process, a natural way to overcome this difficulty is to approximate vector L k by a vector V k in the kernel of AðX kÀ1 Þ instead of AðX k Þ. With this approach, an initial approximation V 0 is needed and will be chosen as the tangent vector v ðiÞ . At each iteration k, the following system will therefore be solved: As proposed in [13], we can replace the preceding algorithm by the following one which splits the solution of system (6) (including the construction of the vector V k ) into two linear systems sharing the same matrix. We note that in the following, T k represents a correction on the tangent vector.
For k ¼ 0; 1; 2; . . . ; k max . 1. Solve the linear system: . Update the solution vector: . Update the tangent vector: Normalization of the tangent vector: If kFðX k Þk 6 e F and kX kþ1 À X k k 6 e x , convergence attained: Note that system in step 2 ensures that AðX k ÞðV k À T k Þ ¼ 0. The vector Z ¼ V k À T k is thus in the kernel of AðX k Þ and therefore, the normalized vector V kþ1 ¼ Z kZk that will be used at the next iteration will be in the kernel of AðX k Þ as expected. A geometric interpretation of the Moore-Penrose algorithm is illustrated in Fig. 1. Starting from a point x ðiÞ on the solution curve, we follow the tangential vector V 0 (tangential to Aðx ðiÞ Þ) to get the initial estimate X 0 .
From there, we follow the direction perpendicular to V 0 to reach X 1 . Defining V 1 as a vector tangential to AðX 0 Þ, we then follow the direction orthogonal to V 1 , and so forth.
To obtain the next point on the solution curve, two systems need to be solved at each iteration. In a finite element context, especially for problems with a high number of degrees of freedom, constructing and solving these two systems can be costly and should be avoided. In the next section, we will show that the explicit construction of these systems, which in matrix form are given by: is not necessary. We note that in systems (10) and (11), vector V k was decomposed into its components in the u and k directions The correction vectors were also decomposed in the same way (i.e. d k

Schur complement
The Schur complement of the block M 11 of a matrix where x and a are column vectors of dimension p and y and b are column vectors of dimension q. As for the blocks M 11 ; M 12 ; M 21 and M 22 , they are matrices of dimension p Â p; p Â q; q Â p and q Â q respectively. Assuming M 11 is nonsingular the Schur complement is invertible (see [23]) and we get The Schur complement approach presented here is a mathematical tool used in various circumstances in structural analysis or in fluid dynamics (for instance, the static condensation as proposed by Bathe [7] uses the Schur complement approach on slave nodes). Since the Schur complement is a matrix of dimension q Â q, the Schur complement approach thus reduces the problem of inverting a ðp þ qÞ Â ðp þ qÞ matrix (matrix M) to that of inverting two matrices, a p Â p matrix (matrix M 11 ) and a q Â q matrix (the Schur complement).
In the case where q ¼ 1, this approach is very advantageous since the Schur complement, which we will denote by b, is a scalar. The solution of system (12) is therefore easily obtained.
The solution of system (12) is therefore obtained by solving two systems, both of which involve the same matrix M 11 .
For the Moore-Penrose algorithm, we have two systems of the form (12) to solve (see systems (10) and (11)), where the matrix M is the same for both systems and is given by: k is a block of dimension 1 Â 1 and matrix M 11 is given by F 0 u ðX k Þ which is nothing but the classical Jacobian matrix in a classical Newton method. Using a Schur complement approach to solve these two systems would therefore necessitate solving 4 linear systems involving F 0 u ðX k Þ. Fortunately, the specific form of systems (10) and (11) allows some simplifications as we will see in the next section.

Schur complement approach applied to the Moore-Penrose method
Applying the Schur complement approach to system (10) leads to the following steps: at each iteration k.
Concerning system (11), we note that it is not necessary to redo all the above steps. Indeed, the computation of W k u and b k is exactly the same for both systems. Moreover, step 1 gives: and from step 4 and 5 we get The following algorithm summarizes the implementation technique of the Moore-Penrose continuation method when using a Schur complement approach.
For k ¼ 0; 1; 2; . . . ; k max . 1. (a) Solve to obtain the standard Newton correction (for a fixed value of k): Calculate the Moore-Penrose correction for the load parameter k: Apply the Moore-Penrose correction on the standard Newton correction: (a) Calculate the correction for the k component of the tangent vector: Calculate the correction for the u component of the tangent vector: Update the solution vector: . Update the tangent vector: Normalization of the tangent vector: If kFðX k Þk 6 e F and kX kþ1 À X k k 6 e x , convergence attained:

.1. Important remarks
Both systems needed for the Moore-Penrose continuation method can thus be solved easily using the Schur complement approach. Furthermore, one interesting thing to notice is that d k u , the Moore-Penrose correction for the N degrees of freedom at iteration k, is calculated using the formula d k u;k is the correction obtained from the standard Newton method. The Moore-Penrose corrections can thus simply be obtained by modifying the corrections obtained without continuation. This also means that the Moore-Penrose continuation algorithm can be incorporated numerically in the finite element code by simply adding a postcondition to the Newton method solver. This postcondition is added to modify the correction obtained by the Newton method in order to obtain the desired Moore-Penrose correction. This postcondition can be activated when necessary and does not require the construction of matrices that are not already constructed for the finite element method. When using a direct solver to obtain the standard Newton correction D k u;k , solving the additional system (to obtain W k u ) is not very costly since the LU decomposition of the Jacobian matrix is already known. When an iterative method is used (e.g. in the case of very large finite element 3D problems) to obtain the correction D u;k , this approach is also advantageous as the same iterative method can be used to solve the additional system and obtain W k u . Iterative methods are very sensitive, and there is no guarantee that the iterative method developed to obtain D k u;k would enable us to solve the augmented systems (10) and (11). As the systems needed for most predictor-corrector continuation methods have the same form as the ones needed for the Moore-Penrose continuation method, we note that a similar approach could be used for other continuation methods as well.

Numerical algorithm
The following algorithm describes how the Moore-Penrose continuation method (based on a Schur complement approach) can be incorporated numerically in a finite element code. We will assume that the load is applied in terms of k, the continuation parameter which varies from 0 to 1. At k ¼ 1, the desired load will have been applied completely and the solution obtained will be the desired one. If the algorithm converges for 0 < k < 1, the solution will be considered as an intermediate solution and the algorithm will continue until convergence is attained for k ¼ 1 (if possible).
The Moore-Penrose continuation is added in the nonlinear loop that solves the system given by Fðu; kÞ ¼ 0 and can be described by Algorithm 1. We note that the algorithm assumes that the variational formulation is solved using the Newton method and that X k ¼ ½u k k k > .
Convergence attained END end if end if end for In Algorithm 1, a first order forward finite difference scheme (in regards to k) is used to approximate the initial tangent V 0 . The step size used for this finite difference scheme is denoted by e k while Tol u and Tol k are the desired tolerances for u and k respectively. We also note that the derivative of F with respect to k is also calculated (approximated) with a first order forward finite difference scheme.
Using this approach, we therefore only exit the nonlinear solver loop once the desired load has been attained at k ¼ 1. In an updated Lagrangian context (see [7,20]), since the geometry is only updated once the convergence is attained, this means that the desired load should not be imposed for the complete simulation, but only for one updated Lagrangian step. Obviously, if the imposed load is too high, the algorithm will not be able to converge to k ¼ 1 without remeshing. In this situation, a mechanism is therefore needed to exit the loop and update the mesh at the largest load attained. We note that this mechanism, which enables a user to exit the loop with either the last converged solution or the second last converged solution (if we want a state of the solution that enables convergence for at least one more step), was added to Algorithm 1.
Nonconvergence before k ¼ 1 could also mean that the step size h is too large. At each converged intermediate solution, the algorithm always tries to find a next point on the solution curve. To find this next point, the first step always consists in completing a prediction step of step size h in the tangent direction. If h is too large, the prediction point will be too far from the solution curve and convergence problems will be encountered. A mechanism therefore also has to be added to the algorithm in order to modify the value of h if needed. This mechanism must verify, at each intermediate solution, if the value of h should be increased or reduced to obtain the next point. In general, the number of iterations needed for convergence is the criterion used to determine if the value of h should be modified or not, and if so, in which way. This mechanism should also enable a restart from the last converged solution (with a reduced value of h) when the algorithm does not converge (and this before exiting the loop completely). We note that these two mechanisms were implemented numerically in our finite element code, but were not mentioned in Algorithm 1 in order to simplify the reading.

Detection of bifurcations
The approach proposed in Section 3 for the implementation of the Moore-Penrose continuation method facilitates the detection of bifurcation points. In fact, as we will see, most of the informations needed for the detection of bifurcation points are already known when using a Schur complement approach to solve the nonlinear system of equations, which means that important information about the problem can be obtained at a low computational cost. In this work, we will concentrate on the detection of simple bifurcation points only.
Let x 2 R Nþ1 be a simple bifurcation point of equation FðxÞ ¼ 0 (see [18] for a formal definition). We thus have that Fð xÞ ¼ 0 and that two distinct continuous solutions passes through x (as illustrated in Fig. 2). To detect a simple bifurcation point, we can use the following result.  The detection of bifurcations can therefore be made by solving an additional system (i.e. system (14)), which is simply an augmented system of the system used for the Moore-Penrose continuation. Clearly, for ss, the solution of this system is As for the vectors B and C, we note that these can be chosen randomly due to the high probability that they will satisfy the conditions stated above. In a finite element context (with continuation), the system to solve can be written as: where the vectors B; C and Y have been decomposed into their components in the u and k directions respectively (i.e. B ¼ ½B u B k > , . Using the Schur complement approach to find G ¼ ½G u G k > in (16), we get Finally, by replacing expressions (18) in (16), we obtain: At a given iteration k, the solution of system (17) is thus given by: F 0 k ðX k Þ. We note that b k and W k u are already known and do not need to be recalculated. To obtain the solution, we therefore only need to solve one more linear system with the matrix F 0 u ðX k Þ (to calculate P k u ) as well as compute a few scalar products. Furthermore, if the goal is only to detect bifurcation points, vector Y (i.e. Y u and Y k ) does not need to be calculated (because the bifurcation can be detected with the sign of s only). We also note that s does not need It is important to comment a little further on the computation of the critical value s. The expression (15) (or simply the determinant of the Jacobian matrix) can be used for small matrices like those encountered in the study of small dynamical systems to detect bifurcation points. For large systems as those obtained in finite element simulations (when an LU decomposition is used), the computation of these determinants would end up with machine overflows. This is why we use the expression given in Eq. (16). However, while following the solution curve, vector G can become orthogonal to vector C thus producing a singularity in the computation of s. This is also an indication that the matrix MðxðsÞÞ, defined by is singular for some value of s. The algorithm, as it is now, may detect a change of sign even though s is not defined there. To avoid this situation, a supplementary precaution is added to the algorithm when a change of sign is detected. We simply make another choice for vectors B and C in system (14) in order to verify if these new vectors also generate a change in the sign of s. To avoid unnecessary computations, vector B is left unchanged in our algorithm for the second verification test, but vector C is chosen as J À1 ðxðs 0 ÞÞB where s 0 represents the last converged point before the detected point.
Making this choice only requires to store in memory the last computed value of vector G and to compute the scalar product of this vector with vector G at the current converged point to see if the sign is positive. If so, it implies that s does not change sign with the use of these new vectors and that consequently the detected point is not a real bifurcation point.
As for vector Y, it can be used to follow a specific bifurcation branch. In fact, in the case of a simple bifurcation point where only two distinct solution branches intersect, vector Y can be used as a prediction direction (instead of the usual tangent vector) to find a point on the new branch. More specifically, by denotingx an approximation of x (obtained by looking at the sign of s) and YðxÞ the first part of the solution of the augmented system (14) used to compute sðxÞ, one can take YðxÞ as a prediction direction and complete one continuation step starting from the point ðx; YðxÞÞ to obtain a first point on the new branch. Once this first continuation step has been completed, the standard continuation method can be used to follow the new branch.

Elastic beam problem
To verify that the bifurcation detection algorithm works well, we can first consider the classical elastic beam buckling problem. As the bifurcation modes for this problem are well known, it is    Table 1 Point detected at v ¼ 0:032893 (second point of the pair): comparison of the six smallest eigenvalues. Comparison is made just before (v1 ¼ 0:032425) and just after (v2 ¼ 0:033150) the detected point.

Antisymmetric modes
Many bifurcation modes exist for the elastic beam problem (Hill and Hutchinson [16], Bardet [5], Triantafyllidis et al. [26]). On the main solution branch, the beam should remain straight and continue to widen with increasing load. This is precisely what is observed in Fig. 3. Using the procedure described in Section 4, we have been able to detect different bifurcations and this is illustrated in Fig. 4.
As an additional precaution, we present in Table 1 the six smallest eigenvalues for the different matrices used during the computation in the region of v ¼ 0:0329 (the second point in Fig. 4). We see clearly that the Jacobian matrix F 0 u has an eigenvalue that passes through zero (which confirms the presence of a bifurcation point). As expected, the same thing is also true for the Moore-Penrose matrix JðxÞ. The matrix M is in turn invertible. Now, by using the approach described in Section 4 to change branches, we are able to bifurcate towards the different bifurcation modes (which are all antisymmetric). Fig. 5 illustrates the different solutions obtained while Fig. 6 shows the bifurcation diagram for these same solutions. We note that the curve that passes through the detected bifurcation points (indicated by the symbol () in the figure) represents the fundamental branch.

Symmetric modes
Symmetric bifurcation modes also exist for the elastic beam problem. As described in Brun et al. [9], these bifurcations become possible after the accumulation point of the antisymmetric modes (which defines a surface instability). However, it is important to mention that these modes are generally difficult to obtain numerically with asymmetrical meshes. We can however hope to obtain these modes by using a mesh completely symmetric from the start. We thus used a regular symmetric mesh (not illustrated) composed of 4320 quadrilaterals (20 Â 216) and 4557 nodes and we have been able to detect the different antisymmetric modes (at approximately the same values of v than for the asymmetric mesh of Fig. 3(a)). We have also obtained the surface instability at an approximate load of v ¼ 0:451 (which is characterized by a very wavy deformation mode when we approach it, as is illustrated in Fig. 7) as well as the symmetric modes afterwards. However, as the symmetric bifurcation modes are often very close to one another, they become very difficult to detect (because if we move too quickly along the solution curve, the algorithm will not be able to detect the two consecutive changes in the sign of s). To detect these modes, it is thus preferable to first look at the appearance of the vector Y at the different intermediate solutions in order to see the different mode changes appear. Afterwards, the different bifurcations can be detected more precisely by moving slowly on the solution curve in the regions where the different mode changes occur. Fig. 8 illustrates some of the symmetric modes for the elastic beam problem and indicates the approximative load at which these different bifurcations occur. As we can see, the symmetric bifurcation modes decrease with increasing load, which is different from the antisymmetric modes where the opposite is true. For the antisymmetric modes, the greater the load, the wavier the bifurcation modes are. For the symmetric modes, the wavier modes occur at smaller loads. We note that this phenomenon was also observed in the numerical results of [9]. We also note that Fig. 8 illustrates the appearance of vector Y at the different bifurcation points and not the bifurcated solutions (i.e. solutions that we would obtain if we were to follow the different branches). However, to show that we can also follow the different symmetric branches, we have used the approach described in Section 4 to bifurcate towards the first symmetric mode which takes place at a load of v ¼ 0:747. Fig. 9 shows the solution obtained.

Computational time
To show the efficiency of the proposed algorithm, we have also solved the buckling beam problem with the classical approach described in Section 2. Table 2 compares the computational times when solved using the asymmetrical mesh of Fig. 3(a). For this test, a total load of v ¼ 0:2 was imposed (while staying on the main solution branch) and the same level of accuracy was attained with both algorithms. We also note that the detection of bifurcation points algorithm was not activated during the simulation as it had never been implemented in the classical algorithm (due to the computational cost). However, the cost associated with the detection of bifurcation points when using the new algorithm is very low (approximately 0.1 of a second).

Three-dimensional case
We can also validate the approach in the tridimensional case. Let us consider the beam illustrated in Fig. 10(a), which is a 3D generalization of the previous example, and solve using Q 2 À Q 1 elements. We note that the regular mesh is composed of 1250 hexahedra (5 Â 5 Â 50) and 6696 nodes.
In the same way as for the 2D case, we illustrate the evolution of the deformation of the beam (on the fundamental branch) in Fig. 10. As for the bifurcations detected by the algorithm, we have illustrated them in Fig. 11. As we can see, the same bifurcations are detected in the 3D case, which is normal since we are solving the same problem but with an extra dimension. We were also able to bifurcate towards the different bifurcation modes, as illustrated in Fig. 12.

Finite strain incompressible elasticity problem
To test the robustness of the proposed approach, let us now consider the finite strain incompressible elasticity problem recently presented in [3,4] where a square of initial configuration X 0 ¼ ½À1; 1 Â ½À1; 1, with coordinates x and y and meshed as in Fig. 14(a), is considered. For this problem, the top edge of the configuration will be traction free, while all the other edges will be fixed. The body will be subjected to a vertical uniform body force F and a Mooney-Rivlin model with parameters k ¼ 40; 000; 000; c 1 ¼ 20 and c 2 ¼ 0 will be used to describe the incompressible elastic material. We note that these parameters correspond to the parameters used in [4] with the value of k chosen to represent a purely incompressible elastic material. P 2 À P 1 elements will also be used for the computation.
Without continuation, the standard Newton approach converges to the trivial solution of this problem, which is given by u ¼ 0 and p ¼ ÀFð1 À yÞ, until convergence problems are encountered at a load of approximately F ¼ 265. This load corresponds to the critical value detected in [4]. Fig. 13 illustrates the deformation associated with this load. As we can see, the solution seems to have diverged from the fundamental branch (which corresponds to the trivial solution) at this breaking point and seems to indicate the presence of a bifurcation point. Our proposed Moore-Penrose continuation algorithm can then be activated to see if a bifurcation point was encountered and if better results can be attained.
By activating our proposed algorithm, the solution initially converges on the fundamental branch and a bifurcation point is detected at an approximate value of F ¼ 264 (which explains the convergence problems of the standard Newton approach). By following the bifurcation branch associated with this point, much larger deformations than with the standard Newton approach can be attained, as can be seen in Fig. 14 where the evolution of the deformation of the configuration over time is shown.
We also note that by staying on the fundamental branch (instead of following the bifurcation branch associated with the bifurcation point at F ¼ 264), we were able to detect another  bifurcation point for this problem at an approximate value of F ¼ 370. Fig. 15 shows the two bifurcation modes associated with this bifurcation point. High level of deformations were also reached by following these branches with our algorithm.
The proposed Moore-Penrose continuation algorithm has therefore proved to be a very powerful tool when dealing with these kinds of problem. Not only were we able to detect the different bifurcation points, but we were also able to attain much larger deformations than in the case of a standard Newton approach.

Conclusion
In this work, we have presented a new implementation technique based on a Schur complement approach for the Moore-Penrose continuation method. For this method, no additional matrices need to be constructed. It simply requires adding a postcondition to the Newton method solver. This postcondition modifies the Newton correction accordingly in order to obtain the desired Moore-Penrose correction. This postcondition can be activated when needed without modifying the implementation of the standard finite element method. Detection of bifurcation points and branch following are also simplified with this approach as most of the informations needed are already available. The numerical examples presented illustrated the good performance of this method.