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

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


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


Введение


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


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


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


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


$ a_i\frac{\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_i\int\limits_\Omega \nabla u_i \nabla v \,dx + \int\limits_\Gamma b_i u_i v \,d\Gamma + \sum_{k=1}^N \int\limits_\Omega q_{ik} u_k v \,dx = \int\limits_\Omega g_i v \,dx + \int\limits_\Gamma b_i w_i v \,d\Gamma, \;\;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_i\int\limits_\Omega \nabla u \nabla v \,dx + \int\limits_\Gamma b_i u v \,d\Gamma, \quad B_{ik}(u, v) = \sum_{k=1}^N \int\limits_\Omega q_{ik} u v \,dx, $


$ F_i(v) = \int\limits_\Omega g_i v \,dx + \int\limits_\Gamma b_i w_i v \,d\Gamma. $


Алгоритм МКЭ


Введем в пространстве $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) = \displaystyle\sum_{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_i\int\limits_{\Omega} \nabla\phi_s \cdot \nabla\phi_j \,dx + \int\limits_\Gamma b_i\phi_s\phi_j \right)c_{is} + \sum_{k=1}^N\sum_{s=1}^m \left(\int_\Omega q_{ik}\phi_s\phi_j \,dx \right)c_{ks} = $


$ = \int\limits_\Omega g_i\phi_j \,dx + \int\limits_\Gamma b_i w_i \phi_j \,d\Gamma, \;\; i = 1, \ldots, N, \; j = 1, \ldots, m. $


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


Отметим, что если $s$-й и $j$-й узлы не соединены ребром, то $\varphi_s\varphi_j \equiv 0$, $\nabla\varphi_s \cdot \nabla\varphi_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]$:


$ \int\limits_i^j f(x, y) \,dl = \frac{|i\, j|}{2} \int\limits_{-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$.


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


$ \int\limits_{-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\}$:


$ \int\limits_{i\,j\,k} f(x,y) \,dxdy = 2 |i\,j\,k| \int\limits_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++, триангуляция сохраняется в формате .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.


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


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


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


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


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


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


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


Заключение


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


Источники


  1. Алексеев Г.В. Численные методы математической физики. Введение в метод конечных элементов. Владивосток: Дальнаука, 1999.
  2. Даутов Р.З., Карчевский М.М. Введение в теорию метода конечных элементов. Казань: КГУ, 2004.
  3. Zienkiewicz O.C., Taylor R.L., Zhu J.Z. The finite element method: its basis and fundamentals. Elsevier, 2005.
  4. Сегерлинд Л. Применение метода конечных элементов. М.: Мир, 1979.
  5. Метод конечных элементов на примере уравнения Пуассона // Хабрахабр, 27.07.2015.
  6. Написание МКЭ расчетчика в менее чем 180 строк кода //
    Хабрахабр, 01.12.2015.
  7. Практика использования Freefem++ // Хабрахабр, 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. Почему юнит-тесты не работают в научных приложениях // Хабрахабр, 26.04.2010.
  11. Гренкин Г.В. Алгоритм решения задачи граничного оптимального управления в модели сложного теплообмена // Дальневост. матем. журн. 2016. Т. 16. № 1. С. 24–38.

© Habrahabr.ru