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

2017-12-23 в 2:31, admin, рубрики: c++, freefem++, вычислительная математика, дифференциальные уравнения, математика, математическая физика, математическое моделирование, метод конечных элементов, МКЭ, Программирование, численные методы

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

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

Введение

Метод конечных элементов [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++, триангуляция сохраняется в формате .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].

Источники

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

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

Источник

Информация

Обсуждаемое

Рекомендуем

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

Введение

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

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