Principle of resolution by elimination ========================================= Definition of the problem to be solved --------------------------------- The system to be solved can take the following general form: .. math:: :label: eq-1 \ left [K\ right]\ {u\} =\ {f\} under duress .. math:: :label: eq-2 \ left [B\ right]\ {u\} =\ {{u} _ {0}\} With the introduction of Lagrange multipliers, the problem becomes .. math:: :label: eq-3 \ {\ 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} and is put in the general form .. math:: : label: eq-4 \ left [\ begin {array} {cc} K& {B} ^ {B} ^ {T}\\ B& 0\ end {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 .. math:: : 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 :math:`{\lambda }_{1}={\lambda }_{2}=\lambda /2`. To make its semi-definite character positive to the generic problem :math:`\mathit{Ku}=f`, in the case where the matrix :math:`K` presents a topology resulting from the layout in the form () or (), we will try to decompose the solution :math:`u` of the problem in the form :math:`u={u}_{p}+{u}_{g}`, where * :math:`{u}_{p}` is the particular solution to the problem, * :math:`{u}_{g}` is the general solution to the problem, associated with the :math:`[B]\{{u}_{g}\}=\{0\}` boundary conditions. Search for the particular solution ------------------------------------- The search for the particular solution is simply finding :math:`{u}_{p}` such as :math:`[B]\{{u}_{p}\}=\{{u}_{0}\}`. To calculate :math:`{u}_{p}`, we consider the decomposition on the image and the core of :math:`B`, i.e. :math:`{u}_{p}={u}_{p-i}+{u}_{p-n}`, with :math:`{u}_{p-i}` belonging to the image of :math:`B` and :math:`{u}_{p-n}` belonging to the core of :math:`B`. This problem can be solved in a least-squares sense, and in this case, we consider the solution of () with the minimum :math:`{L}_{2}` norm, that is :math:`{u}_{p}={u}_{p-i}`. :math:`\{{u}_{p-i}\}` is the solution to the minimization problem: .. math:: :label: eq-6 \ {{u} _ {p-i}\} =\ underset {\ {u} _ {{u} _ {p}\}} {\ text {argMin}}\ left ({\ Vert [B]\ {{u} _ {u} _ {p}\\} _ {{u}\}}\ green} ^ {2}\ right) To efficiently calculate :math:`{u}_{p-i}`, provided that :math:`B` is of full rank (rank equal to :math:`\mathit{Nb}`, so no redundant constraints), the following property is used: .. math:: :label: eq-7 \ left ([B] {[B]} ^ {T}\ right) {\ left ([B] {[B]} ^ {T}\ right)} ^ {-1} = {\ left [{I} _ {I} _ {d}\ right]}} _ {\ mathit {nbxNB}} = {\ left [{I}} _ {d}\ right]} And so .. math:: :label: eq-8 \ left ([B] {[B]} ^ {T}\ right) {\ left ([B] {[B]} ^ {T}\ right)} ^ {-1}\ {{u} _ {0}\} _ {0}\} _ {0}\} The solution sought is naturally given by: .. math:: :label: eq-9 \ {{u} _ {p}\} =\ {{u} _ {{u} _ {p-i}\} = {[B]} ^ {T} {\ left ([B]}} ^ {T}\ right)} ^ {T}\\ right)} ^ {T}\ right)} ^ {-1}\ -1}\ {-1}\} In practice, we perform a Cholesky factorization of the symmetric matrix :math:`{\mathit{BB}}^{T}`, using the MUMPS library, then we solve .. math:: :label: eq-10 [B] {[B]} ^ {T}\ {y\} =\ {{u} _ {0}\} Finally, we get :math:`\{{u}_{p}\}`: .. csv-table:: ":math:`\{{u}_{p}\}=\{{u}_{p-i}\}={[B]}^{T}\{y\}`" Finding the general solution --------------------------------- Since the general solution :math:`{u}_{g}` checks for :math:`[B]\{{u}_{g}\}=\{0\}`, it belongs to the core of :math:`B`. The search for the general solution therefore involves a step of building a :math:`Z` base of the :math:`B` core. We can then project the problem onto the core, and solve this problem, which is better posed and smaller in size. Obtaining a core base ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Formally, building the :math:`B` core can be done in several ways: * Decomposition into single values of :math:`B`, which is :math:`B={\mathit{USV}}^{T}`. :math:`\mathit{Ker}(B)` is the space generated by the columns of :math:`{V}^{T}` associated with a zero singular value. * Decomposition :math:`\mathit{QR}` of :math:`{B}^{T}`, which is :math:`B={R}^{T}{Q}^{T}`. :math:`\mathit{Ker}(B)` is the space generated by the columns of :math:`Q` associated with a zero value on the diagonal of :math:`{R}^{T}`. * Decomposition :math:`\mathit{LU}` of :math:`{B}^{T}`. The kernel calculation is not a direct step, and is obtained from sub-blocks of :math:`L`. It is this third method that is implemented in code_aster. We use the SuperLu library to build factorization :math:`\mathit{LU}` of :math:`{B}^{T}`: .. math:: :label: eq-11 {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`: .. math:: :label: eq-12 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 :math:`\mathit{Ker}(B)` is given by: .. math:: :label: eq-13 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 :math:`B` is of full rank. If this is not the case, then simply build :math:`{L}_{1}` from the block of :math:`L` associated with the non-zero terms of the diagonal of :math:`U`. Projection and resolution ~~~~~~~~~~~~~~~~~~~~~~~~~~ Once the core base is calculated, the resolution is direct. We ask :math:`{u}_{g}=Z\overline{u}`, and the solution of the problem is therefore directly given by .. math:: :label: eq-14 {[Z]} ^ {T} [K] [Z] [Z]\ {\ overline {u}\} = {Z} ^ {T}\ {f\} since, by construction, .. math:: :label: eq-15 {[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 :math:`{u}_{p}` and :math:`{u}_{g}`. To do this, we need to solve the first line of the system (), or .. math:: :label: eq-16 {[B]} ^ {T}\ {\ mathrm {\ lambda}\}} =\ {f\} - [K]\ left ([Z]\ {\ overline {u}\}} +\ {{u}\} +\ {{u} _ {p}\}\ right) If :math:`B` is of full rank, the solution is formally written: .. math:: :label: eq-17 \ {\ 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 :math:`B{B}^{T}` calculated for the search for the particular solution. Similarly, if :math:`B` is not of a full rank, it is possible to remove redundant constraints.