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

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

Введение

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

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

Краевая задача Робина для системы $$ уравнений диффузии-реакции в ограниченной области $$ с границей $$:

Здесь $$, $$ — вектор внешней нормали к $$.

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

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

Введем функциональное пространство Соболева $$. Функции $$ называются слабым решением исходной краевой задачи, если указанные тождества выполняются для произвольной тестовой функции $$.

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

где $$ — билинейные формы, $$ — линейная форма:

Алгоритм МКЭ

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

Пусть $$ — базис $$. Если положить $$, $$, то задача сводится к СЛАУ относительно $$ неизвестных $$:

или

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

Отметим, что если $$-й и $$-й узлы не соединены ребром, то $$, $$, поэтому при суммировании по $$ достаточно перебрать только индексы $$, для которых между $$-м и $$-м узлами есть ребро либо $$.

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

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

Чтобы задать функцию $$, нужно указать ее значения в узлах сетки. Обозначим координаты $$-го узла через $$, а значение функции $$ в $$-м узле — через $$.

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

Отсюда

где $$ — площадь треугольника со знаком.

Значение функции можно найти по формуле $$.

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

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

Здесь $$ — длина ребра, $$,
$$.

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

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

Сведем интеграл по треугольнику с $$-й, $$-й и $$-й вершинами к интегралу по треугольнику $$:

Здесь $$ — площадь треугольника, $$,
$$.

Для вычисления последнего интеграла применим квадратурную формулу [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.

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

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

Разобьем объемлющий прямоугольник области $$ на $$ блоков, где $$, $$, $$ — длина и ширина объемлющего прямоугольника, $$ — средние размеры объемлющих прямоугольников треугольников триангуляции.

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

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

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

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

Заключение

Преимущество разработанной библиотеки перед крупными солверами — ее простота, а также возможность учесть специфические особенности задачи (кусочно-постоянный коэффициент в граничном условии). В дальнейшем планируется реализовать аналогичную библиотеку для 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.