Como implementar de forma eficiente as condições de contorno de Dirichlet em matrizes globais de rigidez esparsa de elementos finitos

Estou me perguntando como as condições de contorno de Dirichlet em matrizes esparsas de elementos finitos globais são realmente implementadas com eficiência. Por exemplo, digamos que nossa matriz global de elementos finitos fosse:

$$K = \begin{bmatrix} 5 & 2 & 0 & -1 & 0 \\ 2 & 4 & 1 & 0 & 0 \\ 0 & 1 & 6 & 3 & 2 \\ -1 & 0 & 3 & 7 & 0 \\ 0 & 0 & 2 & 0 & 3 \end{bmatrix} \hspace{5mm}\text{and right-hand side vector}\hspace{5mm} b = \begin{bmatrix} b1 \\ b2 \\ b3 \\ b4 \\ b5 \\ \end{bmatrix}$$

Em seguida, para aplicar uma condição Dirichlet no primeiro nó ( ), a primeira linha, colocamos 1 em e a primeira coluna do lado direito. Por exemplo, nosso sistema se tornaria: K $x_{1}=c$$K_{11}$

$$K = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 4 & 1 & 0 & 0 \\ 0 & 1 & 6 & 3 & 2 \\ 0 & 0 & 3 & 7 & 0 \\ 0 & 0 & 2 & 0 & 3 \end{bmatrix} \hspace{5mm}\text{and right-hand side vector}\hspace{5mm} b = \begin{bmatrix} c \\ b2-2\times{c} \\ b3-0\times{c} \\ b4+1\times{c} \\ b5-0\times{c} \\ \end{bmatrix}$$

Isso é bom em teoria, mas se nossa matriz K for armazenada em formato de linha compactada (CRS), mover as colunas para o lado direito se tornará caro para sistemas grandes (com muitos nós sendo dirichlet). Uma alternativa seria não mover as colunas correspondentes a uma condição de Dirichlet para o lado direito, ou seja, nosso sistema se tornaria:

$$K = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 2 & 4 & 1 & 0 & 0 \\ 0 & 1 & 6 & 3 & 2 \\ -1 & 0 & 3 & 7 & 0 \\ 0 & 0 & 2 & 0 & 3 \end{bmatrix} \hspace{5mm}\text{and right-hand side vector}\hspace{5mm} b = \begin{bmatrix} c \\ b2 \\ b3 \\ b4 \\ b5 \\ \end{bmatrix}$$

No entanto, isso tem uma grande desvantagem, pois o sistema não é mais simétrico e, portanto, não podemos mais usar gradiente conjugado pré-condicionado (ou outros solucionadores simétricos). Uma solução interessante que me deparei é o "Método de grandes números", que encontrei no livro "Programação de elementos finitos em Java", de Gennadiy Nikishkov. Este método usa o fato de que a precisão dupla contém apenas cerca de 16 dígitos de precisão. Em vez de colocar um 1 na posição , colocamos um número grande. Por exemplo, nosso sistema se torna: K = [ 1.0 e 64 2 0 - 1 0 2 4 1 0 0 0 1 6 3 2 - 1 0 3 7 0 0 0 2 0 3$K_{11}$

$$K = \begin{bmatrix} 1.0e64 & 2 & 0 & -1 & 0 \\ 2 & 4 & 1 & 0 & 0 \\ 0 & 1 & 6 & 3 & 2 \\ -1 & 0 & 3 & 7 & 0 \\ 0 & 0 & 2 & 0 & 3 \end{bmatrix} \hspace{5mm}\text{and right-hand side vector}\hspace{5mm} b = \begin{bmatrix} c\times{1.0e64} \\ b2 \\ b3 \\ b4 \\ b5 \\ \end{bmatrix}$$

As vantagens deste método são que ele mantém a simetria da matriz e também é muito eficiente para formatos de armazenamento esparsos. Minhas perguntas são as seguintes:

Como as condições de contorno de Dirichlet são normalmente implementadas em códigos de elementos finitos para calor / fluidos? As pessoas costumam usar o método de grandes números ou fazem outra coisa? Existe alguma desvantagem no método de grandes números que alguém pode ver? Suponho que provavelmente exista algum método eficiente padrão usado na maioria dos códigos comerciais e não comerciais que resolva esse problema (obviamente, não espero que as pessoas conheçam todo o funcionamento interno de todos os solucionadores de elementos finitos comerciais, mas esse problema parece básico / fundamental o suficiente para que alguém tenha trabalhado em tais projetos e possa fornecer orientação).

Em deal.II ( http://www.dealii.org - exoneração de responsabilidade: sou um dos principais autores dessa biblioteca), eliminamos linhas e colunas inteiras e, em geral, não é muito caro. O truque é usar o fato de que o padrão de dispersão é tipicamente simétrico, para que você saiba em quais linhas você deve procurar ao eliminar uma coluna inteira.

A melhor abordagem, na minha opinião, é eliminar essas linhas e colunas nas matrizes de células, antes de serem adicionadas à matriz global. Lá você trabalha com matrizes completas, para que tudo seja eficiente.

Nunca ouvi falar da abordagem de grandes números e não a usaria porque certamente levaria a problemas terrivelmente mal condicionados.

Para referência, os algoritmos que usamos em deal.II são descritos conceitualmente nas aulas 21.6 e 21.65 em http://www.math.colostate.edu/~bangerth/videos.html . Eles correspondem à sua descrição.

Zero BCs são fáceis. Para BCs diferentes de zero, você também pode usar multiplicadores Lagrange. Por exemplo, veja aqui . Uma vantagem dos LMs é que você pode usar qualquer equação de restrição, embora o sistema se torne indefinido, portanto você precisa de um solucionador apropriado.