- PVSM.RU - https://www.pvsm.ru -

Программирование метода конечных элементов

Данная статья посвящена собственной реализации (солвер Joker FEM [1]) метода конечных элементов для систем уравнений диффузии-реакции.

Обычно предпочтительнее использовать готовые решения, однако если в задаче есть специфические особенности, то на основе простой библиотеки задачу решить легче.

Введение

Метод конечных элементов [1-4] — группа эффективных методов решения задач математической физики, механики сплошных сред. Вопросам программирования МКЭ посвящены источники [5-9].

Постановка задачи

Краевая задача Робина для системы $N$ уравнений диффузии-реакции в ограниченной области $Omega subset mathbb R^2$ с границей $Gamma$:

$ -a_iDelta u_i(x) + sum_{k=1}^N q_{ik}(x)u_k(x)=g_i(x), ;; x in Omega, $

$ a_ifrac{partial u_i(x)}{partial mathbf n} + b_i(x)u_i(x)=b_i(x)w_i(x), ;; x in Gamma. $

Здесь $i=1, 2, ldots, N$, $mathbf n$ — вектор внешней нормали к $Gamma$.

Слабая формулировка

Умножим каждое из $N$ уравнений на произвольную тестовую функцию $v$, проинтегрируем полученные тождества по $Omega$, применим формулу Грина и учтем граничные условия. Получим

$ a_iintlimits_Omega nabla u_i nabla v ,dx + intlimits_Gamma b_i u_i v ,dGamma + sum_{k=1}^N intlimits_Omega q_{ik} u_k v ,dx=intlimits_Omega g_i v ,dx + intlimits_Gamma b_i w_i v ,dGamma, ;;i=1, ldots, N. $

Введем функциональное пространство Соболева $V=H^1(Omega)$. Функции $u_1, ldots, u_N in V$ называются слабым решением исходной краевой задачи, если указанные тождества выполняются для произвольной тестовой функции $v in V$.

Слабую формулировку можно записать в следующем виде: найти функции $u_1, ldots, u_N in V$ такие, что

$ A_i(u_i, v) + sum_{k=1}^N B_{ik}(u_k, v)=F_i(v) ;; forall v in V, ; i=1, ldots, N, $

где $A_i, B_{ik}$ — билинейные формы, $F_i$ — линейная форма:

$ A_i(u, v)=a_iintlimits_Omega nabla u nabla v ,dx + intlimits_Gamma b_i u v ,dGamma, quad B_{ik}(u, v)=sum_{k=1}^N intlimits_Omega q_{ik} u v ,dx, $

$ F_i(v)=intlimits_Omega g_i v ,dx + intlimits_Gamma b_i w_i v ,dGamma. $

Алгоритм МКЭ

Введем в пространстве $V$ конечномерное подпространство $V_m$. Согласно методу Галёркина, приходим к дискретной задаче: найти функции $u_1, ldots, u_N in V_m$ такие, что

$ A_i(u_i, v) + sum_{k=1}^N B_{ik}(u_k, v)=F_i(v) ;; forall v in V_m, ; i=1, ldots, N. $

Пусть $phi_1, ldots, phi_m$ — базис $V_m$. Если положить $u_i(x)=displaystylesum_{s=1}^m c_{is}phi_s(x)$, $v(x)=phi_j(x)$, то задача сводится к СЛАУ относительно $Nm$ неизвестных $c_{ks}$:

$ sum_{s=1}^m A_i(phi_s, phi_j)c_{is} + sum_{k=1}^N sum_{s=1}^m B_{ik}(phi_s, phi_j)c_{ks}=F_i(phi_j), ;; i=1, ldots, N, ; j=1, ldots, m, $

или

$ sum_{s=1}^m left(a_iintlimits_{Omega} nablaphi_s cdot nablaphi_j ,dx + intlimits_Gamma b_iphi_sphi_j right)c_{is} + sum_{k=1}^Nsum_{s=1}^m left(int_Omega q_{ik}phi_sphi_j ,dx right)c_{ks}=$

$=intlimits_Omega g_iphi_j ,dx + intlimits_Gamma b_i w_i phi_j ,dGamma, ;; i=1, ldots, N, ; j=1, ldots, m. $

Простейший пример подпространства $V_m$ — пространство непрерывных кусочно-линейных функций. Область $Omega$ разбивается на конечные элементы (обычно треугольники), и на каждом треугольнике функции из $V_m$ линейны. В МКЭ базисные функции выбираются так, чтобы на большей части области они обращались в ноль, тогда в матрице СЛАУ будет много нулей. В качестве базисных функций $phi_j$ берут "пирамидки", которые равны 1 в одном узле сетки и 0 в остальных узлах. Каждому узлу сетки соответствует своя базисная функция, и размерность $m$ подпространства равна числу узлов сетки.

Отметим, что если $s$-й и $j$-й узлы не соединены ребром, то $varphi_svarphi_j equiv 0$, $nablavarphi_s cdot nablavarphi_j equiv 0$, поэтому при суммировании по $s$ достаточно перебрать только индексы $s$, для которых между $s$-м и $j$-м узлами есть ребро либо $s=j$.

Таким образом, решение краевой задачи методом конечных элементов сводится к построению и решению соответствующей СЛАУ.

Кусочно-линейные функции

Чтобы задать функцию $u in V_m$, нужно указать ее значения в узлах сетки. Обозначим координаты $i$-го узла через $(X_i, Y_i)$, а значение функции $u$ в $i$-м узле — через $U_i=u(X_i, Y_i)$.

Чтобы найти значение функции $u$ в произвольной точке $(x, y)$, принадлежащей треугольнику с $i$-й, $j$-й и $k$-й вершинами, введем барицентрические координаты $L_1, L_2, L_3$, которые принимают значения от 0 до 1 и связаны с декартовыми координатами $x, y$ следующими соотношениями (см. [3, разд. 4.7.1], [4, с. 35]):

$ begin{gathered} x=L_1 X_i + L_2 X_j + L_3 X_k, \ y=L_1 Y_i + L_2 Y_j + L_3 Y_k, \ 1=L_1 + L_2 + L_3. end{gathered} $

Отсюда

$ begin{gathered} L_1(x, y)=frac{1}{2A}(a_i + b_i x + c_i y), ; a_i=X_j Y_k - X_k Y_j, ; b_i=Y_j - Y_k, ; c_i=X_k - X_j, \ L_2(x, y)=frac{1}{2A}(a_j + b_j x + c_j y), ; a_j=X_k Y_i - X_i Y_k, ; b_j=Y_k - Y_i, ; c_j=X_i - X_k, \ L_3(x, y)=frac{1}{2A}(a_k + b_k x + c_k y), ; a_k=X_i Y_j - X_j Y_i, ; b_k=Y_i - Y_j, ; c_k=X_j - X_i, \ end{gathered} $

где $A=dfrac{1}{2}begin{vmatrix} 1 & X_i & Y_i \ 1 & X_j & Y_j \ 1 & X_k & Y_k end{vmatrix}$ — площадь треугольника со знаком.

Значение функции можно найти по формуле $u(x, y)=U_i L_1(x,y) + U_j L_2(x,y) + U_k L_3(x,y)$.

Численное интегрирование

Сведем интеграл по граничному ребру с $i$-м и $j$-м концами к интегралу по отрезку $[-1, 1]$:

$ intlimits_i^j f(x, y) ,dl=frac{|i, j|}{2} intlimits_{-1}^1 f(x(t), y(t)) ,dt. $

Здесь $|i, j|$ — длина ребра, $x(t)=dfrac{X_i + X_j}{2} + dfrac{X_j - X_i}{2}t$,
$y(t)=dfrac{Y_i + Y_j}{2} + dfrac{Y_j - Y_i}{2}t$.

Для вычисления последнего интеграла применим квадратурную формулу Гаусса:

$ intlimits_{-1}^1 varphi(s) ,ds approx sum_{j=1}^n w_j varphi(xi_j). $

Точки $xi_j in [-1, 1]$ и весовые множители $w_j$ берутся из таблицы [3, разд. 5.9.2].

Сведем интеграл по треугольнику с $i$-й, $j$-й и $k$-й вершинами к интегралу по треугольнику $D={(L_1, L_2) in [0,1] times [0,1] colon L_1 + L_2 leq 1}$:

$ intlimits_{i,j,k} f(x,y) ,dxdy=2 |i,j,k| intlimits_D f(x(L_1, L_2), y(L_1, L_2)) ,dL_1dL_2. $

Здесь $|i,j,k|=|A|$ — площадь треугольника, $x(L_1,L_2)=L_1X_i + L_2X_j + (1 - L_1 - L_2)X_k$,
$y(L_1, L_2)=L_1Y_i + L_2Y_j + (1 - L_1 - L_2)Y_k$.

Для вычисления последнего интеграла применим квадратурную формулу [3, разд. 5.11].

Программная реализация

Общие соображения

В данной предметной области тестирование затруднено [10], поэтому важно, чтобы код был написан предельно ясно. ООП позволяет повысить читабельность кода, для применения ООП нужно выделить набор понятий и сформулировать алгоритм в этих терминах. Это дает возможность в первую очередь обращать внимание на алгоритм, а только потом на детали реализации.

Триангуляция

Для построения сетки используется пакет FreeFem++ [2], триангуляция сохраняется в формате .mesh.

В библиотеке вводится класс Mesh, содержащий информацию о сетке, которая считывается из файла.

Построение СЛАУ

Для записи алгоритма вводятся классы базисных функций и подынтегральных функций (произведений функций), которые унаследованы от базовых классов Integrand и BoundaryIntegrand. В этих классах есть поля support (носитель) и boundary_support (граничный носитель) соответственно. Носитель — это список треугольников, на которых функция не равна 0. Граничный носитель — это список граничных ребер, на которых функция не равна 0. Также для классов подынтегральных функций определен метод Value, который вычисляет значение функции в заданной точке.

Функции Integrate и BoundaryIntegrate принимают на вход подынтегральную функцию, перебирают список-носитель и вызывают функцию интегрирования по заданному треугольнику (или граничному ребру), которая определена в базовом классе Integrand (BoundaryIntegrand) и вызывает виртуальную функцию Value.

Таким образом, код построения СЛАУ имеет простой вид:

for (int i = 0; i < data.N; ++i) {
  for (int j = 0; j < data.mesh.nodes_num; ++j) {
    for (auto s : data.mesh.adjacent_nodes[j]) {
      sys.AddCoeff(i, j, i, s,
        data.a[i] * Integrate(grad(phi[s]) * grad(phi[j]))
          + BoundaryIntegrate(data.b[i] * phi[s] * phi[j])
      );
    }
    for (int k = 0; k < data.N; ++k) {
      for (auto s : data.mesh.adjacent_nodes[j]) {
        sys.AddCoeff(i, j, k, s,
          Integrate(data.q[i][k] * phi[s] * phi[j])
        );
      }
    }
    sys.AddRhs(i, j,
      Integrate(data.g[i] * phi[j])
        + BoundaryIntegrate(data.b[i] * data.w[i] * phi[j])
    );
  }
}

Решение СЛАУ

Для решения СЛАУ применяется библиотека MTL4 [3].

Доступ к результату

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

Разобьем объемлющий прямоугольник области $Omega$ на $nm$ блоков, где $n=[W/a]$, $m=[H/b]$, $W, H$ — длина и ширина объемлющего прямоугольника, $a, b$ — средние размеры объемлющих прямоугольников треугольников триангуляции.

Для каждого из $nm$ блоков заполним список треугольников, которые его пересекают. Для этого переберем все треугольники и найдем их объемлющие прямоугольники,
для которых легко найти пересекающие их блоки.

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

Тестирование

Тестирование программы проводилось на задачах с известным точным решением. При удвоенном измельчении сетки погрешность уменьшается примерно в 4 раза или в 2 раза, что свидетельствует о том, что метод имеет 2-й или 1-й порядок точности. Уменьшение порядка сходимости до 1-го может быть связано с несогласованностью граничных условий на стыках частей границы, на которых заданы разные значения функции $w_i$.

Заключение

Преимущество разработанной библиотеки перед крупными солверами — ее простота, а также возможность учесть специфические особенности задачи (кусочно-постоянный коэффициент в граничном условии). В дальнейшем планируется реализовать аналогичную библиотеку для 3-мерных задач, которая будет использована для решения задач оптимального управления граничным коэффициентом в модели сложного теплообмена, где данный коэффициент является кусочно-постоянным [11].

Источники

  1. Алексеев Г.В. Численные методы математической физики. Введение в метод конечных элементов. [4] Владивосток: Дальнаука, 1999.
  2. Даутов Р.З., Карчевский М.М. Введение в теорию метода конечных элементов. [5] Казань: КГУ, 2004.
  3. Zienkiewicz O.C., Taylor R.L., Zhu J.Z. The finite element method: its basis and fundamentals. Elsevier, 2005.
  4. Сегерлинд Л. Применение метода конечных элементов. М.: Мир, 1979.
  5. Метод конечных элементов на примере уравнения Пуассона [6] // Хабрахабр, 27.07.2015.
  6. Написание МКЭ расчетчика в менее чем 180 строк кода [7] //
    Хабрахабр, 01.12.2015.
  7. Практика использования Freefem++ [8] // Хабрахабр, 06.06.2013.
  8. Gockenbach M.S. Understanding and implementing the finite element method. SIAM, 2006.
  9. Smith I.M., Griffiths D.V., Margetts L. Programming the finite element method. Wiley, 2014.
  10. Почему юнит-тесты не работают в научных приложениях [9] // Хабрахабр, 26.04.2010.
  11. Гренкин Г.В. Алгоритм решения задачи граничного оптимального управления в модели сложного теплообмена [10] // Дальневост. матем. журн. 2016. Т. 16. № 1. С. 24-38.

Автор: Глеб Гренкин

Источник [11]

Сайт-источник PVSM.RU: https://www.pvsm.ru

Путь до страницы источника: https://www.pvsm.ru/programmirovanie/271460

Ссылки в тексте:

[1] Joker FEM: https://github.com/grenkin/joker-fem

[2] FreeFem++: http://www.freefem.org/

[3] MTL4: http://www.simunova.com/mtl4

[4] Численные методы математической физики. Введение в метод конечных элементов.: https://vk.com/doc195326829_201397317

[5] Введение в теорию метода конечных элементов.: http://kpfu.ru/staff_files/F1229619272/MKEbookDRZ.pdf

[6] Метод конечных элементов на примере уравнения Пуассона: https://habrahabr.ru/post/263597/

[7] Написание МКЭ расчетчика в менее чем 180 строк кода: https://habrahabr.ru/post/271723/

[8] Практика использования Freefem++: https://habrahabr.ru/post/182396/

[9] Почему юнит-тесты не работают в научных приложениях: https://habrahabr.ru/post/92038/

[10] Алгоритм решения задачи граничного оптимального управления в модели сложного теплообмена: http://mi.mathnet.ru/dvmg319

[11] Источник: https://habrahabr.ru/post/344564/?utm_source=habrahabr&utm_medium=rss&utm_campaign=344564