3. Principle of resolution by elimination

3.1. Definition of the problem to be solved

The system to be solved can take the following general form:

(3.1)\[ \ left [K\ right]\ {u\} =\ {f\}\]

under duress

(3.2)\[ \ left [B\ right]\ {u\} =\ {{u} _ {0}\}\]

With the introduction of Lagrange multipliers, the problem becomes

(3.3)\[\begin{split} \ {\ begin {array} {c}\ left [K\ right]\ {k\ right]\ {\ right]\ {u\} + {\ left [B\ right]} ^ {T}\ {\ mathrm {\ lambda}\}\\ left [K\ right]\\\\ left [K\ right]\ {u\ right]} ^ {T}\ {\ mathrm {\ lambda}\}\} =\ {f\}}\ {f\}}\\ f\}\}\\\\\ left\\\ left [B\\ right]\\\ left [B\ right]\ {u\ right]} =\ {{u} _ {0}\}\ end {array}\end{split}\]

and is put in the general form

\[\]

: label: eq-4

left [begin {array} {cc} K& {B} ^ {B} ^ {T}\ B& 0end {array}right]left{begin {array} {c} u\ mathrm {lambda}\ lambda}\lambda}\{lambda}\\ lambda}\\lambda}\}\ lambda}\\{lambda}\\\lambda}right}

In the case of double dualization (AFFE_CHAR_MECA), this matrix is increased, and the system to be solved becomes

\[\]

: label: eq-5

left [begin {array} {ccc} {ccc} K& {B}} ^ {B} ^ {T}\ B& -mathrm {alpha} {I} _ {d} &mathrm {alpha} {I} {I} {I} _ {d} _ {d} & -mathrm {alpha} {I} _ {d}end {array}right]left{begin {array} {c} u\ {mathrm {lambda}} _ {1}\ {mathrm {lambda}}} _ {2}end {array}right}} =left{begin {array} {lambda}}} _ {lambda}} _ {lambda}} _ {lambda}} _ {lambda}} _ {lambda}} _ {2}\ end {array}right} =left{begin {array} {lambda}}} _ {lambda}}} _ {lambda}} _ {lambda}} _ {0}\ {u} _ {0}end {array}right}

after naturally posing \({\lambda }_{1}={\lambda }_{2}=\lambda /2\).

To make its semi-definite character positive to the generic problem \(\mathit{Ku}=f\), in the case where the matrix \(K\) presents a topology resulting from the layout in the form () or (), we will try to decompose the solution \(u\) of the problem in the form \(u={u}_{p}+{u}_{g}\), where

• \({u}_{p}\) is the particular solution to the problem,

• \({u}_{g}\) is the general solution to the problem, associated with the \([B]\{{u}_{g}\}=\{0\}\) boundary conditions.

3.2. Search for the particular solution

The search for the particular solution is simply finding \({u}_{p}\) such as \([B]\{{u}_{p}\}=\{{u}_{0}\}\). To calculate \({u}_{p}\), we consider the decomposition on the image and the core of \(B\), i.e. \({u}_{p}={u}_{p-i}+{u}_{p-n}\), with \({u}_{p-i}\) belonging to the image of \(B\) and \({u}_{p-n}\) belonging to the core of \(B\).

This problem can be solved in a least-squares sense, and in this case, we consider the solution of () with the minimum \({L}_{2}\) norm, that is \({u}_{p}={u}_{p-i}\). \(\{{u}_{p-i}\}\) is the solution to the minimization problem:

(3.4)\[\begin{split} \ {{u} _ {p-i}\} =\ underset {\ {u} _ {{u} _ {p}\}} {\ text {argMin}}\ left ({\ Vert [B]\ {{u} _ {u} _ {p}\\} _ {{u}\}}\ green} ^ {2}\ right)\end{split}\]

To efficiently calculate \({u}_{p-i}\), provided that \(B\) is of full rank (rank equal to \(\mathit{Nb}\), so no redundant constraints), the following property is used:

(3.5)\[ \ left ([B] {[B]} ^ {T}\ right) {\ left ([B] {[B]} ^ {T}\ right)} ^ {-1} = {\ left [{I} _ {I} _ {d}\ right]}} _ {\ mathit {nbxNB}} = {\ left [{I}} _ {d}\ right]}\]

And so

(3.6)\[ \ left ([B] {[B]} ^ {T}\ right) {\ left ([B] {[B]} ^ {T}\ right)} ^ {-1}\ {{u} _ {0}\} _ {0}\} _ {0}\}\]

The solution sought is naturally given by:

(3.7)\[\begin{split} \ {{u} _ {p}\} =\ {{u} _ {{u} _ {p-i}\} = {[B]} ^ {T} {\ left ([B]}} ^ {T}\ right)} ^ {T}\\ right)} ^ {T}\ right)} ^ {-1}\ -1}\ {-1}\}\end{split}\]

In practice, we perform a Cholesky factorization of the symmetric matrix \({\mathit{BB}}^{T}\), using the MUMPS library, then we solve

(3.8)\[ [B] {[B]} ^ {T}\ {y\} =\ {{u} _ {0}\}\]

Finally, we get \(\{{u}_{p}\}\):

\(\{{u}_{p}\}=\{{u}_{p-i}\}={[B]}^{T}\{y\}\)

3.3. Finding the general solution

Since the general solution \({u}_{g}\) checks for \([B]\{{u}_{g}\}=\{0\}\), it belongs to the core of \(B\). The search for the general solution therefore involves a step of building a \(Z\) base of the \(B\) core. We can then project the problem onto the core, and solve this problem, which is better posed and smaller in size.

3.3.1. Obtaining a core base

Formally, building the \(B\) core can be done in several ways:

• Decomposition into single values of \(B\), which is \(B={\mathit{USV}}^{T}\). \(\mathit{Ker}(B)\) is the space generated by the columns of \({V}^{T}\) associated with a zero singular value.

• Decomposition \(\mathit{QR}\) of \({B}^{T}\), which is \(B={R}^{T}{Q}^{T}\). \(\mathit{Ker}(B)\) is the space generated by the columns of \(Q\) associated with a zero value on the diagonal of \({R}^{T}\).

• Decomposition \(\mathit{LU}\) of \({B}^{T}\). The kernel calculation is not a direct step, and is obtained from sub-blocks of \(L\).

It is this third method that is implemented in code_aster. We use the SuperLu library to build factorization \(\mathit{LU}\) of \({B}^{T}\):

(3.9)\[ \begin{align}\begin{aligned} {P} _ {r} {B} ^ {T} {P} {P} _ {c}} ^ {T} =\ left [L\ right]\ left [U\ right] =\ left [\ begin {array} {cc} {cc} {cc} {cc}} {[C}} {{L} _ {1}]}} {[L} _ {1}]}}} _ {{N} _ {1}]}} _ {{N} _ {b}}\ times {N} _ {b}} & [:ref:\\The matrices :math:`{P}_{r}` and :math:`{P}_{c}` respectively designate the permutation matrices imposed on the rows and columns of :math:`{B}^{T}` for the construction of the factorization :math:`\mathit{LU}`. Blocks :math:`{L}_{1}` and :math:`{L}_{3}` are lower triangular and block :math:`{U}_{1}` is upper triangular.\\Note:\\In practice, :math:`{P}_{c}=\mathit{Id}`, there is no permutation of the columns in :math:`{B}^{T}`.\\We derive the following expression for :math:`B`:\end{aligned}\end{align} \]

(3.10)\[ B= {U} ^ {T} {L} ^ {T} {T} {P} {P} {P} _ {r} =\ left [\ begin {array} {cc} {[U} _ {1}]}} _ {{N} _ {1}]}}} _ {{N} _ {1}]}} _ {{N} _ {1}]}} _ {{N} _ {1}]}} _ {{N} _ {1}]} _ {N} _ {b} _ {b}}\ times {N} _ {b} _ {b}}\ times {N} _ {b}\]

It is then easy to see that a base of \(\mathit{Ker}(B)\) is given by:

(3.11)\[ Z= {P} _ {r} ^ {-1}}\ left [\ begin {array} {c}} - {[{L} _ {1}]} ^ {-T} {[{L} _ {2}]}} ^ {2}]} ^ {T}}} ^ {T}}} ^ {T}}} ^ {T}}}} ^ {L} _ {2}]}} ^ {2}]} ^ {2}]} ^ {2}]} ^ {2}]} ^ {2}]} ^ {2}]} ^ {2}]} ^ {T}]} ^ {2}]} ^ {T}]} ^ {2}]} ^ {T}]} ^ {T}]} ^ {2}]} ^ {T}]} ^ {T}}\ end {array}\ right]\]

Here we are making the assumption that \(B\) is of full rank. If this is not the case, then simply build \({L}_{1}\) from the block of \(L\) associated with the non-zero terms of the diagonal of \(U\).

3.3.2. Projection and resolution

Once the core base is calculated, the resolution is direct. We ask \({u}_{g}=Z\overline{u}\), and the solution of the problem is therefore directly given by

(3.12)\[ {[Z]} ^ {T} [K] [Z] [Z]\ {\ overline {u}\} = {Z} ^ {T}\ {f\}\]

since, by construction,

(3.13)\[ {[Z]} ^ {T} {[B]} {[B]} ^ {T}\ {\ mathrm {\ lambda}\} =\ {0\}\]

In some cases, it is necessary to access Lagrange multipliers. They must therefore be rebuilt, based on solutions \({u}_{p}\) and \({u}_{g}\). To do this, we need to solve the first line of the system (), or

(3.14)\[ {[B]} ^ {T}\ {\ mathrm {\ lambda}\}} =\ {f\} - [K]\ left ([Z]\ {\ overline {u}\}} +\ {{u}\} +\ {{u} _ {p}\}\ right)\]

If \(B\) is of full rank, the solution is formally written:

(3.15)\[ \ {\ mathrm {\ lambda}\}} = {\ left ([B] {[B]} ^ {T}\ right)} ^ {-1} [B]\ {f\} - [K]\ left ([T] [T]\ {\ left ([T]\ {\ overline {u}\}\ right)\]

and is calculated in practice by reusing the factorization of \(B{B}^{T}\) calculated for the search for the particular solution. Similarly, if \(B\) is not of a full rank, it is possible to remove redundant constraints.