В ранее опубликованном почта мы говорили об использовании библиотек с открытым исходным кодом Freefem++ и NetGen в программе аэродинамического моделирования.
В этой статье мы более подробно рассмотрим основные функции Freefem++ и представим небольшое введение в его язык ввода.
Это предоставит исходную информацию, которая часто нужна разработчикам при выборе сторонних компонентов для включения в разрабатываемое ими приложение.
Получение и запуск Freefem++
Freefem++ — программа, предназначенная для решения математических задач методом конечных элементов.Freefem++ можно использовать как в средах Windows, так и в Linux. Официальный сайт www.freefem.org/ff++ .
В нашем случае использовалась Freefem++ версии 3.18.-1. Архивный пакет, скачанный с официального сайта, содержит три загрузочных модуля, назначение которых следующее:
- FreeFem++.
exe — полноценная версия программы, которая при необходимости отображает графическое окно с результатами решения.
При необходимости при вызове можно указать два переключателя: -nowait, чтобы гарантировать завершение работы Freefem++ без подтверждения пользователя, и -nw, чтобы отказаться от отображения графического окна с результатами решения.
- FreeFem++-mpi.exe – для распараллеливания вычислений на основе MPI.
- FreeFem++-nw.exe – для вызова из приложений игнорирует графический вывод результатов.
Чтобы она не ждала подтверждения завершения работы с программой, нужно в командной строке указать ключ -nowait.
Используя операторы вывода или специальные функции сохранения в сценарии, результаты можно записать в файл для дальнейшего использования.
2D триангуляция
Ранее мы показали, что Freefem++ можно использовать для создания STL-описания 3D-модели, состоящей из элементов с плоскими гранями.Чтобы решить эту проблему, необходимо триангулировать грани.
Давайте посмотрим, как это делается с помощью Freefem++.
Пусть необходимо триангулировать поверхность, ограниченную границами (см.
рис.
1).
Рис.
1. Области поверхности Границы областей, для которых необходимо выполнить триангуляцию, задаются в Freefem++ в параметрическом виде.
Рассмотрим следующий код:
int n = 5;
border b0(t=0,1){x = 7 * t; y = 0; }
border b1(t=0,1){x = 7; y = 0 + 6 * t; }
border b2(t=0,1){x = 7 - 7 * t; y = 6; }
border b3(t=0,1){x = 0; y = 6 - 6 * t; }
border b4(t=0,1){x = 1 + 2 * t; y = 1.5; }
border b5(t=0,1){x = 5; y = 1.5 + 2 * t; }
border b6(t=0,1){x = 5 - 4 * t; y = 5.5; }
border b7(t=0,1){x = 1; y = 5.5 - 4 * t; }
border b8(t=0,1){x = 5; y = 3.5 + 2 * t; }
border b9(t=0,1){x = 3 + 2 * t; y = 1.5; }
border b10(t=0,1){x = 3 + 3 * t; y = 0.5; }
border b11(t=0,1){x = 6; y = 0.5 + 3 * t; }
border b12(t=0,1){x = 6 - 1 * t; y = 3.5; }
border b13(t=0,1){x = 3; y = 3.5 - 2 * t; }
border b14(t=0,1){x = 5 - 2 * t; y = 3.5; }
border b15(t=0,1){x = 3; y = 1.5 - 1 * t; }
plot(b0(n) + b1(n) + b2(n) + b3(n) + b4(n) + b5(n) + b6(n) + b7(n) + b8(n) + b9(n) + b10(n)
+ b11(n) + b12(n) + b13(n) + b14(n) + b15(n));
Freefem++ определяет типы скалярных переменных: int, real, string, bool. Объявление переменной выполняется так же, как в C++, но переменную можно инициализировать так же, как это делается для переменной n.
Границы определяются фрагментами с использованием параметризованных линий.
Для этого и нужен тип границы.
Фрагменты границ могут пересекаться только на своих концах.
Объявление фрагмента границы включает идентификатор границы, границу для изменения параметра (в данном случае параметра t) и выражения для изменения координат x и y в зависимости от параметра.
В примере все фрагменты являются сегментами.
Последний оператор в рассматриваемом коде —plot, он предназначен для графического отображения данных определенных типов.
В этом случае отображается граница, состоящая из фрагментов границы.
Фрагменты границы объединяются в общую границу с помощью оператора +.
Для каждого фрагмента указывается параметр, задающий количество частей, на которые следует разбить фрагмент. Этот вариант будет обсуждаться еще раз ниже в другом контексте.
Результат выполнения указанного оператора графика показан на рис.
2. 
Рис.
2. Границы участков поверхности, отображаемые оператором графика Теперь добавим код, выполняющий триангуляцию.
На рис.
3 показан результат его выполнения.
mesh Th = buildmesh(b0(n) + b1(n) + b2(n) + b3(n) + b4(n) + b5(n) + b6(n) + b7(n) + b8(n) + b9(n) + b10(n) + b11(n) + b12(n) + b13(n) + b14(n) + b15(n), fixeborder=1); plot(Th);
Рис.
3. Результат триангуляции Этот код вводит переменную Th типа mesh. Тип сетки представляет собой программную модель триангуляционной сетки.
Для получения объекта такого типа в данном случае используется функция buildmesh, которой на вход передается граница разделяемых областей.
В функцию buildmesh можно передать необязательный параметр fixeborder, который гарантирует, что точки на линии не будут изменены в процессе триангуляции.
Второй оператор –plot – в данном случае принимает на вход объект типа сетка и отображает сетку элементов в графическом окне (рис.
3).
Как уже было показано, граница представляется как сумма фрагментов границы типа border. Параметр, указанный здесь для каждого фрагмента, имеет двойное назначение.
Во-первых, он определяет, на сколько частей разбивается фрагмент границы для построения сетки.
Каждая часть фрагмента границы станет стороной треугольника сгенерированной сетки.
Во-вторых, знак параметра указывает область, которую следует триангулировать.
Если параметр положительный, то триангуляция производится слева от граничного фрагмента.
Если оно отрицательное, то триангуляция выполняется справа от границы.
Поясним это на примере.
Если бы в приведенном выше коде граница была указана следующим образом: b0(n) + b1(n) + b2(n) + b3(n) + b4(n) + b5(-2) + b6(n) + b7( n) + b8(n) + b9(-2) + b10(n) + b11(n) + b12(n) + b13(-2) + b14(-2) + b15(n), где для b5 b9, b13 и b14 разбиваются на два отрезка с отрицательным значением, тогда триангуляция будет выполняться, как показано на рис.
4. 
Рис.
4. Триангуляция с пустой подобластью.
Добавим код, демонстрирующий операции с переменной сетки, оператор цикла и вывод в файл.
ofstream fout("triangulation.out");
fout.precision(14);
fout << Th.nt << endl;
for(int i = 0; i < Th.nt; i++ )
{
fout
<< Th[i][0].
x << " "
<< Th[i][0].
y << " "
<< Th[i][1].
x << " "
<< Th[i][1].
y << " "
<< Th[i][2].
x << " "
<< Th[i][2].
y << " "
<< endl;
}
Для вывода данных в файл предусмотрен тип ofstream — аналог потока вывода файла из библиотеки fstream для C++.
При инициализации файловой переменной необходимо указать имя файла.
Если вам нужно только добавить данные в файл, то вам нужно указать второй параметр инициализации add. С файловой переменной можно выполнять ряд операций.
В примере функция точности используется для установки количества выводимых цифр после десятичной точки.
Есть еще несколько функций, управляющих форматом вывода, аналогичных манипуляторам из библиотеки iostream. << operator allows you to output variable values of scalar types to a file. In a character string, you can use the line end character '\n', or output the global variable endl for this purpose. Приведенный выше код содержит оператор цикла for, который полностью повторяет логику оператора for в C++.
Операторы увеличения и уменьшения также аналогичны операторам C++.
Freefem++ предоставляет доступ к различным свойствам сетки, представленным типом сетки.
Обратите внимание, что Freefem++ предоставляет одни и те же интерфейсы для работы как с двумерной, так и с трехмерной сеткой3. Поэтому здесь мы будем рассматривать этот интерфейс для двух указанных типов, чтобы не повторять далее его описания для mesh3. Как видно из примера, свойство nt возвращает количество элементов сетки (треугольников для сетки или тетраэдров для сетки3).
Также можно получить количество точек сетки — свойство nv, количество граничных элементов — nbe. Существует два типа индексации, которые можно применить к объекту типа сетка или сетка3: индексация с помощью квадратных скобок «[]» и индексация с круглыми скобками «()».
Первый из них дает элемент сетки (треугольник для сетки или тетраэдр для сетки3) с заданным индексом.
Таким образом, Th[i] — это элемент сетки (в данном примере треугольник), к которому, в свою очередь, применяются определенные функции, предусмотренные для элемента.
Индексация «()» предоставляет вершине заданный индекс.
К вершине можно применить операции по получению ее координат. К элементу сетки применяются следующие функции: индексация для получения вершин элемента, получение метки.
3D сетка
Freefem++ определяет тип mesh3 для 3D-сетки.Такую сетку можно получить с помощью функций генерации.
Например, функция buildlayers() позволяет создавать 3D-сетку путем вытягивания 2D-триангуляции вдоль оси Z. Однако у Freefem++ нет собственного генератора.
Для этих целей используется сторонний генератор TetGen (официальный сайт wias-berlin.de/software/tetgen ), а соответствующие функции Freefem++ являются программным интерфейсом этого генератора.
TetGen не может быть использован в коммерческих проектах без разрешения автора.
Но если такая лицензия создаёт проблемы, то можно использовать другой генератор сеток.
Freefem++ имеет возможность импортировать 3D-сетку из файла, соответствующего формату сетки.
Файл сетки имеет следующий формат:
- В начале файла записывается информация о версии формата (MeshVersionFormatted 1) и размерности пространства (Dimension 3).
MeshVersionFormatted 1 Dimension 3
- После ключевого слова Vertices указывается количество вершин сетки и перечисляются вершины, для которых указаны координаты X, Y, Z и метка вершины.
Vertices 65 7 5 1.5 0 7 3 1.5 0 7 3 0.5 0 7 5 3.5 0 .
- После ключевого слова Tetrahedra указывается количество элементов сетки и перечисляются элементы, для которых записаны четыре номера вершин и метка.
Tetrahedra 185 52 56 49 50 0 3 18 14 61 0 47 64 57 54 0 .
- После ключевого слова Triangles пишется количество граничных треугольников.
Для каждого треугольника указаны номера вершин узлов сетки, являющихся вершинами этих треугольников, и указана метка границы.
Triangles 96 2 3 8 2 2 8 1 2 4 1 7 2 .
- Файл заканчивается ключевым словом End.
mesh3 Th;
Th = readmesh3("example3d.mesh");
plot(Th);
Чтобы импортировать сетку из файла сетки, порядок номеров вершин тетраэдра имеет решающее значение.
В исходном коде Freefem++, в файле Mesh3dn.hpp, можно найти структуру DataTet, в которой есть метод вычисления объёма тетраэдра mesure().
Если вершины тетраэдра в файле сетки расположены в таком порядке, что рассчитанный объем отрицательный, то Freefem++ завершится сбоем.
Поэтому вершины в файле сетки должны быть перечислены, чтобы такой ситуации не возникло.
При необходимости сетку можно отобразить в графическом окне с помощью функции построения графика (см.
рис.
5).
Рис.
5. Отображение 3D-сетки с помощью функции построения графика.
Решение систем дифференциальных уравнений
Основная цель Freefem++ — решение систем уравнений в частных производных с использованием метода конечных элементов.Для этого функции, участвующие в системе уравнений, определяются на сетке конечных элементов.
Freefem++ позволяет вам определить тип пространства FE (пространство конечных элементов).
В зависимости от типа пространства вводятся функциональные объекты определенного пространства.
Например,
fespace Vh(Th, P23d);
Vh f = x^2 + y^2 + z^2;
где Vh — тип пространства на сетке Th, f — функция на этом пространстве.
При объявлении типа пространства конечных элементов указывается основа для элементов.
Она может принимать значения P03d — кусочно-постоянная, P13d — кусочно-линейная, P23d — кусочно-квадратичная и т. д. На рис.
6 изображена функция из примера с квадратичным базисом, а на рис.
7 — та же функция с линейной основе и с заметными искажениями.
Рис.
6. Функция конечных элементов с квадратичным базисом 
Рис.
7. Функция конечных элементов с линейным базисом.
Для решения системы уравнений в частных производных в Freefem++ задается вариационная формулировка задачи и при решении задачи формируются искомые функции.
На основе следующего скрипта мы рассмотрим, как пишется вариационная формулировка и другие элементы программирования.
mesh3 Th = readmesh3("example3d.mesh");
fespace VVh(Th, [P2,P2,P2,P1]);
fespace Vh(Th, P23d);
macro Grad(u) [dx(u),dy(u),dz(u)] //
macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) //
VVh [u1,u2,u3,p];
VVh [v1,v2,v3,q];
func fup = (1-x)*(x)*y*(1-y)*16;
solve Prob([u1,u2,u3,p],[v1,v2,v3,q]) =
int3d(Th,qforder=3)( Grad(u1)'*Grad(v1) + Grad(u2)'*Grad(v2) + Grad(u3)'*Grad(v3)
Теги: #Freefem++ #метод конечных элементов #сетка элементов #триангуляция #дифференциальные уравнения #разработка веб-сайтов #программирование #математика
-
5 Способов Привлечь Потенциальных Клиентов
19 Oct, 24 -
Видеоконференции: Будущее Уже Здесь
19 Oct, 24 -
Теорема Четырех Цветов
19 Oct, 24 -
Видимость Комментариев
19 Oct, 24 -
Вызов Globalsdb. 4-Я Гонка
19 Oct, 24
