Как эффективно реализовать граничные условия Дирихле в глобальных разреженных матрицах жесткости конечных элементов

Мне интересно, как на самом деле эффективно реализуются граничные условия Дирихле в глобальных разреженных матрицах конечных элементов. Например, допустим, что наша глобальная конечно-элементная матрица была:

$$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}$$

Затем применить условие Дирихле к первому узлу ($x_{1}=c$) мы обнуляем первый ряд, ставим 1 на $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}$$

Это все хорошо в теории, но если наша K-матрица хранится в формате сжатых строк (CRS), то перемещение столбцов в правую часть становится дорогим для больших систем (многие узлы являются дирихлетами). Альтернативой было бы не перемещать столбцы, соответствующие условию Дирихле, в правую часть, то есть наша система стала бы:

$$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}$$

Это, однако, имеет существенный недостаток в том, что система больше не является симметричной, и поэтому мы больше не можем использовать предобусловленный сопряженный градиент (или другие симметричные решатели). Одним интересным решением, с которым я столкнулся, является «Метод больших чисел», который я нашел в книге Геннадий Никишкова «Программирование конечных элементов на Java». Этот метод использует тот факт, что двойная точность содержит только около 16 цифр точности. Вместо того, чтобы поставить 1 в$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}$$

Преимущества этого метода в том, что он поддерживает симметрию матрицы, а также очень эффективен для разреженных форматов хранения. Мои вопросы тогда следующие:

Как граничные условия Дирихле обычно реализуются в кодах конечных элементов для тепла / жидкости? Люди обычно используют метод больших чисел или делают что-то еще? Есть ли недостаток в методе больших чисел, который кто-то может увидеть? Я предполагаю, что, вероятно, существует какой-то стандартный эффективный метод, используемый в большинстве коммерческих и некоммерческих кодов, который решает эту проблему (очевидно, я не ожидаю, что люди будут знать все внутренние механизмы каждого коммерческого решателя конечных элементов, но эта проблема кажется базовой / фундаментальной Достаточно того, что кто-то, вероятно, работал над такими проектами и мог обеспечить руководство).

В сделке II ( http://www.dealii.org - : я один из основных авторов этой библиотеки) мы исключаем целые строки и столбцы, и это не слишком дорого в целом. Хитрость заключается в том, чтобы использовать тот факт, что шаблон разреженности, как правило, симметричен, поэтому вы знаете, какие строки нужно просматривать при удалении целого столбца.

На мой взгляд, лучший подход - исключить эти строки и столбцы в матрицах ячеек, прежде чем они будут добавлены в глобальную матрицу. Там вы работаете с полными матрицами, так что все эффективно.

Я никогда не слышал о методе «больших чисел» и не буду его использовать, потому что, несомненно, он приведет к ужасно плохо обусловленным проблемам.

Для справки, алгоритмы, которые мы используем в deal.II, концептуально описаны в лекциях 21.6 и 21.65 на http://www.math.colostate.edu/~bangerth/videos.html . Они точно соответствуют вашему описанию.

Ноль БК легки. Для ненулевых БК вы также можете использовать множители Лагранжа. Например, смотрите здесь . Одним из преимуществ LM является то, что вы можете использовать любое уравнение ограничения, хотя система становится неопределенной, поэтому вам нужен соответствующий решатель.