Solving linear systems

Finite element computations make use of linear solvers to solve assembled linear systems of equations. The libraries Numpy [1], Scipy [2] and Eigen [3] offer capabilities to solve sparse linear systems by different methods. Muscat proposes a single class to wrap all this solvers into one unique simple interface.

Treatment of kinematic constraints

(7)\[\begin{split}\begin{equation} \begin{matrix} \text{minimize } \UFE \text{ for:} & \frac{1}{2} \UFE^T \KFE \UFE - \UFE^T \FFE \\ \end{matrix} \end{equation}\end{split}\]
(8)\[\begin{split}\begin{equation} \begin{matrix} \text{subject to} &\CFE \UFE = \GFE \\ \end{matrix} \end{equation}\end{split}\]

Equation (7) represents the assembled finite element system to be solved under the kinematic constraints defined the equation (8). This problem can be solved by four different approaches.

The constrains represented by the equation (8) are created from the Dirichlet like boundary conditions (like RBE2 and RBE3). In some cases, this system can have linearly dependent equations, for example when a user specified twice some boundary condition in the same entity (e.g. the shared edge in the case of \(u_x=0\) on 2 adjacent faces).

To be able to treat correctly the constraint defined by (8), we must clean the linearly dependent equations, and find (in some cases) a reduced number of dof to eliminate .

Let extract from \(\CFE\) only the non zero columns (\(\CFE^*\)) and define the augmented matrix of \(\CPFE\) as :

\[\begin{equation} \CPFE := \left[ \begin{array}{c|c} \CFE^* & \GFE \end{array} \right] \end{equation}\]

The \(\CPFE^T\) matrix can then be decomposed using a QR decomposition to extract the suitable (orthonormal) base. This is done using the routine of EIGEN SpQR. SpQR is only available if the C++ interface was compiled, a pure Python (using scipy.linalg.qr) backup algorithm is available. This backup approach uses dense linear and can be expensive (in terms of memory and cpu time). This algorithm is capable of finding the rank of the matrix \(\CPFE^T\) (using a prescribed tolerance).

\[\begin{equation} \CPFE \mathbf{P} = \mathbf{Q} \mathbf{R} \end{equation}\]

With, \(\mathbf{P}\) is the column permutation, \(\mathbf{Q}\) is the orthogonal matrix represented as Householder reflectors (in the case using Eigen). \(\mathbf{R}\) is the sparse triangular factor.

In many cases the constraints are on a single node or only on a reduced number of nodes. This make the matrix \(\CPFE^T\) very sparse, and we can compute the \(\mathbf{QR}\) decomposition block by block. For this the adjacency matrix of the matrix \(\CFE^*\) is computed (\(abs(\CFE^*)abs(\CFE^*)^T\)). A not zero \((i,j)\) term of the adjacency matrix means the rows \(i\) and \(j\) of \(\CPFE\) share at least one degree of freedom in common. Then we can calculate the connected components of the adjacency matrix and work on each block independently. As each block is treated (computation of the QR decomposition), the rank is identified to remove redundant equations. At the same time for each of the sub-matrices a subset of dofs equal to the rank is stored to created a list of slave dofs for later use.

For example:

2 Springs

Fig. 13 Two springs system with 4 degrees of freedom, A) initial state B) constrained solution.

For the system of 2 independent springs (4 degrees of freedom) represented in the previous figure, the tangent matrix and the right hand side member are,

\[\begin{split}\begin{equation} \KFE = \begin{bmatrix} 1000 & -1000 & 0 & 0 \\ -1000 & 1000 & 0 & 0 \\ 0 & 0 & 1000 &-1000 \\ 0 & 0 &-1000 & 1000 \end{bmatrix}, \FFE = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}. \end{equation}\end{split}\]

The constraints imposed to the problem are: 1) blockage of the first dof (\(u_0\) to the value 0), 2) kinematic relation between \(u_1\) and \(u_2\); \(u_2 - u_1 = 1\), 3) prescribed solution on \(u_3\) equal to 3. To demonstrate the treatment, each constraint was added twice to the matrix \(\CFE\) :

\[\begin{split}\begin{equation} \CFE = \begin{bmatrix} 1 & 0 &0 &0 \\ 1 & 0 &0 &0 \\ 0 &-1 &1 &0 \\ 0 &-2 &2 &0 \\ 0 & 0 &0 &1\\ 0 & 0 &0 &1 \end{bmatrix}, \GFE = \begin{bmatrix} 0 \\ 0 \\ 1 \\ 2 \\ 3 \\ 3 \end{bmatrix} \end{equation}\end{split}\]

In this case all degrees of freedom are present in matrix \(\CFE\), this means \(\CFE^* = \CFE\). We start by building the adjacency matrix \(\CFE\CFE^T\), we are interested in only the non zero values so we really calculate \(abs(\CFE)abs(\CFE^T) != 0\).

\[\begin{split}\begin{equation} \mathbf{AM} := abs(\CFE)abs(\CFE^T)!=0 \to \begin{bmatrix} 1 & 1 &0 &0 &0 &0 \\ 1 & 1 &0 &0 &0 &0 \\ 0 & 0 &1 &1 &0 &0\\ 0 & 0 &1 &1 &0 &0\\ 0 &0 &0 & 0 &1 &1\\ 0 &0 &0 & 0 &1 &1\\ \end{bmatrix} \end{equation}\end{split}\]

The computation of the connected components can be expressed by the solution vector:

\[\begin{equation} Connected Component Of(\mathbf{AM}) = \begin{bmatrix} 0 & 0& 1 & 1 &2 &2 \end{bmatrix} \end{equation}\]

We can see 3 connected components (rows of A). Each of these components can be treated independent (they are orthogonal). For the first component (line 1 and 2), it is evident that the lines are linearly dependent and we have only one emph{real} constraint. This can be detected by the SpQr routine and only the relevant base vectors are kept (we also normalize the vectors). We also store the rank first dofs of each component to define the slave indices. At the end we obtain 1) a vector of slave dofs of the system \([0,1,3]\), and a clean constants system defined by :

\[\begin{split}\begin{equation} \MFE = \begin{bmatrix} 1 &0 &0 &0 \\ 0 &-\sqrt{\frac{1}{2}} &\sqrt{\frac{1}{2}} &0 \\ 0 &0 &0 &1 \\ \end{bmatrix}, \VFE = \begin{bmatrix} 0 \\ \sqrt{\frac{1}{2}} \\ -3 \\ \end{bmatrix} \end{equation}\end{split}\]

Penalization

The penalization approach is the simplest way to solve the constrained system (in terms of programming). The idea of this method is to add to the original system to be solve a penalization term scaled by a really large value (\(\alpha_p = 10^8\)). The perturbed system becomes :

\[\begin{split}\begin{equation} \KFE_p = \KFE + maxdiag(\KFE)*\alpha_p*(\MFE^T \MFE) \\ \end{equation}\end{split}\]
\[\begin{split}\begin{equation} \FFE_p = \FFE + maxdiag(\KFE)*\alpha_p*(\MFE^T \VFE) \\ \end{equation}\end{split}\]

Because the modified system has the same number of dofs and no base transformation was done the solution of this system gives directly an original system (equations (7) and (8)) solution approximation

Lagrange Multipliers

This method imposes the constraints by adding one extra dof for each row of equation (8). This yields to the following modified system:

\[\begin{split}\begin{equation} \KFE_{\lambda} = \begin{bmatrix} \KFE & \MFE \\ \MFE & 0 \\ \end{bmatrix}, \FFE_{\lambda} = \begin{bmatrix} \FFE \\ \VFE \end{bmatrix}, \end{equation}\end{split}\]

The resolution of this linear system gives the solution vector:

\[\begin{split}\begin{equation} \UFE_{\lambda} = \begin{bmatrix} \UFE \\ \lambda \end{bmatrix}, \end{equation}\end{split}\]

The solution to the original system is then contained in the first part of the vector \(\UFE_{\lambda}\).

In the case that we have an initial a solution \(\UFE^0\), and we want to solve the system \(\KFE_{\lambda}\UFE_{\lambda} = \FFE_{\lambda}\) using this initial guess, we cand define the initial guess \(\UFE_{\lambda}\):

\[\begin{split}\begin{equation} \UFE_{\lambda} = \begin{bmatrix} \UFE^0 \\ \MFE(\FFE-\KFE) \end{bmatrix} \end{equation}\end{split}\]

Note

This technique has the disadvantage of losing the positive definite properties of the original system (zeros on the diagonal of the operator \(\KFE_{\lambda}\))

Substitution

Let decompose the \(\MFE\) (using the indices collected during the QR decomposition) into slave (\(s\)) and master (\(m\)) dofs :

\[\begin{equation} \begin{matrix} \MFE = \left[ \begin{array}{c|c} \MFE_s & \MFE_m \end{array} \right] \end{matrix} \end{equation}\]

The original system becomes:

\[\begin{split}\begin{equation} \begin{bmatrix}%[l] \KFE_s & \KFE_{s,m} \\ \KFE_{m,s} & \KFE_{m} \end{bmatrix} \begin{bmatrix} \UFE_{s} \\ \UFE_{m} \end{bmatrix} = \begin{bmatrix} \FFE_{s} \\ \FFE_{m} \end{bmatrix} \end{equation}\end{split}\]

Then the constraint equation (8) can be rewritten to solve \(\UFE_s\) :

\[\begin{split}\begin{equation} \begin{array}{l} %\MFE_{s}\UFE_s + \MFE_m \UFE_m = \VFE \\ \UFE_s = \MFE_{s}^{-1} (\VFE - \MFE_m \UFE_m ) \end{array} \end{equation}\end{split}\]

Now we can rewrite the \(\UFE\) in the form of :

(9)\[\begin{split}\begin{equation} \UFE = \begin{bmatrix} \UFE_s \\ \UFE_m \\ \end{bmatrix}= \begin{bmatrix} \MFE_{s}^{-1} (\VFE - \MFE_m \UFE_m ) \\ \UFE_m \\ \end{bmatrix}= \begin{bmatrix} \MFE_{s}^{-1} (\VFE ) \\ 0 \\ \end{bmatrix}+ \begin{bmatrix} -\MFE_{s}^{-1} ( \MFE_m ) \\ I \\ \end{bmatrix}\UFE_m \end{equation}\end{split}\]

By defying:

\[\begin{split}\begin{equation} \XFE = \begin{bmatrix} -\MFE_{s}^{-1} ( \MFE_m ) \\ I \\ \end{bmatrix}, \DFE= \begin{bmatrix} \MFE_{s}^{-1} (\VFE ) \\ 0 \\ \end{bmatrix} \end{equation}\end{split}\]

And injecting this expression on the original system we obtain:

\[\begin{equation} \XFE^T \KFE \XFE \UFE_m = \XFE^T \left( \FFE - \KFE \DFE \right) \end{equation}\]
\[\begin{split}\begin{equation} \begin{matrix} \underbrace{ \begin{bmatrix} -\MFE_{s}^{-1} ( \MFE_m ) \\ I \\ \end{bmatrix}^T }_{X^T} \begin{bmatrix}%[l] \KFE_s & \KFE_{s,m} \\ \KFE_{m,s} & \KFE_{m} \end{bmatrix} \underbrace{ \begin{bmatrix} -\MFE_{s}^{-1} ( \MFE_m ) \\ I \\ \end{bmatrix} }_{X} \UFE_m = \dotsb \\ \dotsb \underbrace{ \begin{bmatrix} -\MFE_{s}^{-1} ( \MFE_m ) \\ I \\ \end{bmatrix}^T }_{X^T} \left(\begin{bmatrix} \FFE_{s}\\ \FFE_{m} \\ \end{bmatrix} - \begin{bmatrix}%[l] \KFE_s & \KFE_{s,m} \\ \KFE_{m,s} & \KFE_{m} \end{bmatrix} \underbrace{\begin{bmatrix} \MFE_{s}^{-1} (\VFE ) \\ 0 \\ \end{bmatrix} }_D \right) \\ \end{matrix} \end{equation}\end{split}\]

The modified system involves only \(\UFE_m\) dofs, and can be solved using a standard solver. The solution to the original system is then calculated using equation (9).

Ainsworth Method

The Ainsworth method [4] proposes a general approach to solve the original system of equations without changing the number of degrees of freedom of the system.

Footnotes