четверг, 16 июня 2011 г.

Gmsh. Знакомство

Знакомство со Gmsh'ем я начинал с примеров. Поэтому сейчас на нескольких примерах мы и разберем самые основные и полезные особенности этой программы.
Для того, чтобы задать геометрию модели, триангуляцию которой необходимо получить, существует 3 способа:
  1. Интерактивное задание геометрии в окне Gmsh'а
  2. Написание .geo файла на встроенном языке Gmsh
  3. Импорт модели из CAD программ через различные форматы (.step, .brep, .stl и другие)
Первые 2 способа позволяют задать геометрию собственными средствами Gmsh. Для использования 3 способа необходимо иметь дополнительную программу, в которой вы нарисуете модель, а затем экспортируете ее в формате, который Gmsh сможет прочитать. Обзор таких программ вместе с форматами я непременно сделаю в будущем. Сейчас же остановимся на первых двух пунктах.

Интерактивное задание геометрии
Интерактивное задание геометрии в окне Gmsh'а подходит для тех, кто только-только начал работать с этой программой и хочет понять, что же она вообще из себя представляет. В этом случае полезно ознакомиться с этой методичкой. В остальных случаях я советую задавать геометрию сразу в .geo файле, используя встроенный язык Gmsh.

.geo файл
Как понятно из расширения, .geo (сокр. от geometry) файл - это файл, где задается геометрия для Gmsh'а. По формату это обычный текстовый файл. Поэтому те, кто создает файлы в Windows посредством проводника, могут создать файл с расширением .txt, а затем переименовать его. Загрузить .geo файл можно запустив Gmsh, а затем открыв через меню File->Open (или через горячую клавишу Ctrl+O). Другой способ - передать .geo файл в качестве параметра при запуске Gmsh'а в командной строке:

gmsh my.geo

При этом, если файла с именем my.geo не существует, он будет создан. Данный способ я считаю более предпочтительным, хотя, конечно, это на любителя.
Редактирование .geo файла возможно в любом редакторе текстов, например, в блокноте или редакторе Far'а. Можно отредактировать файл и из Gmsh'а. Для этого нужно нажать кнопку Edit в разделе Geometry окна меню.

Файл задания геометрии имеет следующую структуру:
  1. Точки
  2. Лини
  3. Замкнутые контуры из линий
  4. Поверхности
  5. Замкнутые контуры из поверхностей
  6. Объемы
  7. Физические сущности
Вообще говоря, это структура файла, задающего геометрию объемной фигуры. Если нужно задать геометрию только поверхности, то, очевидно, нужно опустить пункты 5 и 6.
В любом месте файла могут быть инициализированы переменные (объявлять их не нужно). Они используются для параметризации модели.

Пример двумерной модели
Зададим теперь геометрию модели, представленной на рисунке
Как говорится, рассмотрим область P, состоящую из двух непересекающихся подобластей: P1 и P2, соответствующих материалам с различными физическими свойствами. Пусть P1 - область внутри окружности, P2 - все остальное. При этом, чтобы воспользоваться классическим МКЭ, необходимо иметь согласованную сетку во всей области P. Помимо прочего, допустим, нам необходимо задать на нижней границе области краевое условие (неважно какое).
Итак, приступим:

0    // инициализация переменных
1    x0 = 0;
2    x1 = 1;
3    y0 = 0;
4    y1 = 1;
5    x_c = (x0 + x1) / 2;
6    y_c = (y0 + y1) / 2;
7    cl_1 = 0.2;
8    cl_2 = 0.1;
9    rad = 0.25;

10    If (2 * rad >= x1 - x0 || 2 * rad >= y1 - y0)
11      Printf("Error! Diameter is too large!");
12      //Exit;
13    EndIf

14    // точки прямоугольника
15    Point(1) = {x0, y0, 0, cl_1};
16    Point(2) = {x1, y0, 0, cl_1};
17    Point(3) = {x0, y1, 0, cl_1};
18    Point(4) = {x1, y1, 0, cl_1};

19    // точки окружности
20    Point(10) = {x_c, y_c, 0, cl_2};
21    Point(11) = {x_c-rad, y_c, 0, cl_2};
22    Point(12) = {x_c+rad, y_c, 0, cl_2};
21    Point(13) = {x_c, y_c-rad, 0, cl_2};
22    Point(14) = {x_c, y_c+rad, 0, cl_2};

Внизу страницы есть ссылки на рассмотренный здесь .geo файл.
Разберем, что же тут написано.
строка 0: В файлы задания геометрии можно вставлять комментарии - либо // (тогда будет комментироваться строка, начиная с // и до конца строки), либо /*...*/ (тогда комментируется часть текста). Т.е. использование комментариев здесь такое же как в С.
строка 1: Как я уже говорил, в процессе задания геометрии можно (и нужно!) использовать переменные для параметризации модели. Ведь, допустим, квадрат на приведенной выше картинке может оказать прямоугольником, или может измениться его расположение относительно начала координат. В любом случае использование переменных очень удобно. Кстати, их не нужно объявлять, как в С. И еще важный момент - все переменные в языке Gmsh имеют только 2 типа - вещественный и строковый (да-да, даже в операторе цикла For переменная-счетчик имеет вещественный тип. ужас, правда? :) ).
строки 1-4: В данном случае мы инициализируем координаты границ прямоугольника.
строки 5-6: Координаты центра окружности. В данном случае окружность находится в центре прямоугольника.
строки 7-8: Характеристическая длина (characteristic length = cl). Она определяет "мелкость" сетки. Чем больше cl, тем грубее сетка.
строка 9: Радиус окружности.
строки 10-13: Прежде чем строить геометрию, проверим основные параметры. Например, диаметр окружности должен быть меньше каждой из сторон прямоугольника. Аналогичным образом можно проверить, что x1 > x0 и т.д. Обратите внимание на ключевые слова условного оператора If EndIf. Регистр очень важен! Необходимо также отметить, что конструкция If Else EndIf не поддерживается! И еще. Обратите внимание на строку 12. Там закомментирован выход из Gmsh в случае неправильно заданных параметров. Если раскомментировать его и задать неверные исходные данные, Gmsh выдаст в консоль сообщений ошибку "Error! Diameter is too large!" и сразу же закроется, так что вы даже не успеете ничего прочитать. Если же оставить код таким как есть и задать неверные исходные данные, то Gmsh выдаст сообщение об ошибке, но не покажет его (!) (это актуально для версии Gmsh 2.4.2 - если вы пользуетесь другой версией, возможно, поведение будет другим). Чтобы увидеть сообщение, необходимо вручную открыть консоль сообщений. Делается это через меню Tools -> Message Console (или через горячую клавишу Ctrl+L).
строки 14-22: Собственно, задание точек геометрии. Формат задания точки: Point(номер точки) = {x, y, z, cl}; Номер должен быть уникальным (в рамках нумерации точек), но не обязательно последовательным. Так, например, при задании точек новой фигуры (в данном случае - окружности) удобно начать новую нумерацию. Основное правило задания любых элементов геометрии: () - круглые скобки сообщают новый номер, {} - фигурные скобки - известный номер или данные. Задавать точку можно и без указания характеристической длины: Point(n) = {x, y, z}; В этом случае Gmsh сам назначит ей характеристическую длину, однако управлять сгущением сетки вы уже не сможете.

Что вы получите в результате? Загрузите .geo файл с приведенным выше кодом в Gmsh. Вы получите следующую картинку:
Не слишком понятно. Чтобы стало чуть лучше, можно добавить координатные оси и вывести номера точек. Первое делается в меню Tools -> Options -> раздел General (слева) -> вкладка Axes -> в списке Axes mode выбирайте Simple axes, а в строке Axes tics в третьем столбце поставьте 0 (чтобы не выводил ось z - у нас же двумерная модель). Для ясности приведу скрин с обозначениями:
Чтобы вывести номера точек, необходимо в меню Tools -> Options -> раздел Geometry (слева) -> вкладка Visibility -> поставить галочку напротив Point numbers. Опять же скрин:
Итак, в результате мы получим:
Уже хоть что-то понятно. Иногда удобно выводить номера точек, линий, поверхностей или объемов (это все в том же меню, что и номера точек) для наглядности или чтобы не запутаться, когда геометрия сложная. Следующий этап - определение линий.

23    // линии прямоугольника
24    Line(1) = {1, 2};
25    Line(2) = {1, 3};
26    Line(3) = {2, 4};
27    Line(4) = {3, 4};

28    // линии окружности
29    Circle(11) = {11, 10, 13};
30    Circle(12) = {13, 10, 12};
31    Circle(13) = {12, 10, 14};
32    Circle(14) = {14, 10, 11};

строки 23-27: Прямые линии задаются командой Line(номер линии) = {начальная точка, конечная точка}; Номер линии должен быть уникальным (в рамках нумерации линий), но не обязательно последовательным. Все так же, как и с нумерацией точек. Обратите внимание, что Line(1) = {1, 2}; не та же самая линия, что и Line(1) = {2, 1}; То есть у всех линий есть направление и оно важно! В связи с этим я для себя определил, что всегда (ну или почти всегда) задаю линии от точки с меньшими координатами к точке с большими координатами (в том смысле, что координаты тем меньше, чем ближе они к точке (0, 0, 0)). И тогда я всегда знаю, куда направлена моя линия.
строки 28-32: Окружность задается дугами с помощью команды Circle(n) = {начальная точка дуги, центр окружности, конечная точка дуги}; Здесь также важно направление, поэтому я задаю линии в едином обходе - по(против) часовой стрелке(и) - все зависит от того, как смотреть. И еще одна важная особенность, на которой спотыкаются почти все, кто начинает работать со Gmsh'ем: радиус дуги должен быть меньше 180 градусов!!! В самом деле, вы не задумались, зачем для задания окружности мы определили 4 точки на границе? А для того, чтобы задать 4 дуги по 90 градусов каждая! Конечно, можно было определить 3 точки и задать 3 дуги по 120 градусов, но, я считаю, что первый вариант проще.
Итак, что же мы увидим, загрузив .geo файл с приведенным кодом в Gmsh:
Уже похоже на правду. Однако мы еще далеки от конечного результата.
Для того, чтобы построить треугольную сетку, необходимо определить поверхности. А для того, чтобы задать поверхности, нужно определить замкнутые контуры линий. Этим мы сейчас и займемся:

33    // замкнутый контур прямоугольника
34    Line Loop(21) = {1, 3, -4, -2};

35    // замкнутый контур окружности
36    Line Loop(22) = {11, 12, 13, 14};

Код довольно короткий, однако требующий подробного изучения.
строки 33-34: Замкнутый контур из любых линий (прямых, дуг, сплайнов - каких угодно) задается командой (обратите внимание на пробел и регистр!) Line Loop(номер контура) = {последовательность линий}; Номер контура должен быть уникален в рамках нумерации линий (!). По крайней мере, так сказано в документации. Однако замечено, что если поставить номер контура, совпадающий с номером какой-нибудь линии, ничего не изменится. А теперь самое важное - последовательность линий. Она не может быть произвольной. Линии должны следовать одна за другой так, чтобы конец одной линии совпадал с началом другой. Вот где нужно вспомнить про ориентацию линий! Так, например, начав с линии 1 (оканчивающейся точкой 2) продолжить можно только линией 3 (начинающейся с точки 2). Однако линия 3 оканчивается точкой 4, а с 4-ой точки ни одна из линий не начинается. В этом случае нужно "повернуть" линию, которая также оканчивается на 4-ую точку. Это линия 4. Для "поворота" линии служит знак "-". Линию 2, кстати, тоже придется повернуть. Этот момент задания геометрии, на самом деле, является самым сложным.
строки 35-36: Если вы, как и я, для задания линий окружности выберите единое направление, то мучиться с выстраиванием последовательности вам не придется. Здесь все довольно просто.

В результате добавления строк 33-36 в .geo файл модель в сущности не изменится. Мы также еще не можем построить сетку, а вид геометрии останется прежним.
Поэтому самое время определить поверхности:

37    // поверхности
38    Plane Surface(1) = {21, 22};
39    Plane Surface(2) = {22};

строки 37-39: Поверхность можно задать 2-мя командами (обратите внимание на пробел и регистр!): Plane Surface либо Ruled Surface. Plane Surface используется для задания плоских поверхностей, а Ruled Surface - для изогнутых (в документации сказано: "Ruled Surface - это поверхность, которая может быть проинтерполирована с помощью конечноразностной интерполяции"). Поскольку при работе с двумерными моделями последнее не может встретиться в принципе, для задания поверхностей используем команду Plane Surface(номер поверхности) = {основной контур, вложенные контуры...}; Номер поверхности должен быть уникальным в рамках нумерации поверхностей, не обязательно последовательным. А теперь самое интересное. Со строкой 39 все понятно - плоская поверхность задается как замкнутый контур линий, который определяет окружность. Но если мы зададим плоскую поверхность 1 просто как контур 21, то получим несогласованную сетку! То есть в результате мы будем иметь две не связанных друг с другом триангуляции - для прямоугольника и для окружности. Для того, чтобы сетка была согласованной, необходимо из прямоугольного контура "вырезать" контур окружности! Таким образом, мы определим поверхность области P2 (еще не забыли постановку?). Для этого все вложенные контуры, которые мы хотим "вырезать", нужно перечислить через запятую в объявлении поверхности после основного контура. Надеюсь, понятно объяснил :)

Загрузим полученный .geo файл в Gmsh и отключим показ номеров точек (чтобы не мешались). Затем в разделе Mesh, расположенном там же, где и раздел Geometry, найдем кнопку 2D и нажмем ее. Кстати, в раздел Mesh можно перейти, нажав клавишу m. Вот что получим в результате:
Симпатично. Но я уверен, что вы недостаточно удовлетворены качеством сетки (особенно в части аппроксимации границ окружности). Для того, чтобы измельчить сетку, измените значения cl_1 и cl_2. Например, уменьшите их в 2 раза. Вы убедитесь, что сетка станет гораздо лучше.
Для того, чтобы в дальнейшем использовать сетку, сохраним ее. Для этого существует 2 способа.
1. В меню File выбрать Save As (горячая клавиша Ctrl+S), ввести имя файла с расширением .msh (вообще форматов очень много, но это - родной Gmsh'евский формат) и выбрать тип файла (текстовый версии 1.0 или 2.0, или бинарный версии 2.0). Советую не сохранять в формате ASCII 1.0 - это сильно устаревший вариант. Если сохраните в формате Binary 2.0 - то не сможете его просмотреть глазами, а хотелось бы. Хотя бы пару раз. Поэтому выбирайте формат - ASCII 2.0.
2. В меню File выбрать Save Mesh (горячая клавиша Shift+Ctrl+S). В этом случае сетка сохранится в формате ASCII 2.0 с именем my.msh (при условии, что вы использовали файл геометрии my.geo), т.е. никакой свободы выбора :)
Можно, кстати, построить сетку и сохранить ее с минимальными затратами времени и щелканий мышью. Для этого в командной строке наберите:

gmsh my.geo -2

Все просто. Двумерная сетка построится и сохранится с именем my.msh в формате ASCII 2.0. Если вы хотите еще и имя файла сетки задать, то наберите

gmsh my.geo -2 -o anothername.msh

Добавьте в командную строку -bin и получите файл сетки в формате Binary 2.0. Конечно, этот способ лишит вас наслаждения наблюдать сетку на экране монитора... Но иногда этого и не нужно :)

Что же, думаете мы все сделали? Как бы не так!
Ведь мы хотим различать материалы из областей P1 и P2. А как определить какой треугольник из какой области? Давайте заглянем в .msh файл.

Структура .msh файла
Будем рассматривать структуру файла в формате ASCII 2.0. Как уже говорилось, формат ASCII 1.0 сильно устарел, а Binary 2.0 не даст возможности убедиться в значениях интересующих нас параметров просто заглянув в файл. Кстати, подробное описание всех этих форматов есть в документации к Gmsh'у. Итак,

$MeshFormat // в этой секции для нет ничего интересного (по крайней мере - пока)
2.1 0 8
$EndMeshFormat
$Nodes
n // количество узлов
i x y z
// узлы в формате: номер (i =1,..,n) координаты
...
$EndNodes
$Elements
m // количество элементов
i type ntag t1 t2 t3 nodes
// элементы в формате: номер (i = 1,..,m) тип элемента (линия, треугольник и т.д.) число тэгов (в версии Gmsh 2.4.2 их 3, в версии Gmsh 2.5.0 - не обязательно) тэг1 (физическая область) тэг2 (геометрическая область) тэг3 (область из разделения) узлы
...
$EndElements

Заглянем в наш файл my.msh.
В первую очередь нас интересует секция $Elements.
Итак, какие типы элементов нам встречаются:  15, 1 и 2, т.е. точки, линии и треугольники соответственно. Но, какие еще точки, и что за линии, спросите вы. Ведь нам нужны только треугольники. При этом посмотрев на тэги треугольников, мы увидим либо 0 1 0, либо 0 2 0. Здесь 1 и 2 - это как раз номера поверхностей, которые мы задавали в .geo файле. Таким образом, чтобы различить какой треугольник относится к какой области можно воспользоваться вторым тэгом. Но что делать, если у вас не 2 области, а несколько десятков, имеющих разные геометрические номера, но характеризующие материалы с одинаковыми свойствами? Что делать с ненужными элементами сетки (точки, линии)? И наконец, как нам "выловить" нижнюю грань области для задания краевого условия (не забыли?)? На все эти вопросы ответ один - использовать физические сущности (physical entities).

Физические сущности
Для того, чтобы сохранить в файл только то, что нам требуется, и объединить несколько разных геометрических элементов, используются, так называемые, физические сущности (или объекты).
Определить их можно с помощью команды Physical Point/Line/Surface/Volume(номер) = {геометрические объекты};
Рассмотрим использование физических объектов на нашем примере.
Добавим в конец нашего .geo файла следующие строки:

40    // физические сущности
41   Physical Line(2001) = {1};
42   Physical Surface(1001) = {1};
43   Physical Surface(1002) = {2};

строка 41: Определяется физическая линия 2001, состоящая из геометрической линии 1 - нижней границы области. Теперь эта линия попадет в .msh файл и мы легко сможем определить точки, в которых нужно задавать краевое условие, по номеру 2001.
строки 42-43: Определяются физические поверхности 1001 и 1002. По этим номерам мы всегда отыщем треугольники, принадлежащие разным областям и соответствующие материалам с различными физическими свойствами. Кроме того,  никакой другой элемент сетки, кроме указанных как физические, в файл сетки .msh не попадет.

Вы можете скачать готовые .geo файлы с приведенным выше кодом в двух версиях:
1. с полными комментариями
2. почти без комментариев

Ну вот теперь, кажется, все. Рассмотрение трехмерной модели сделаем в следующей раз.

2 комментария:

  1. Вот незадача, у меня почему-то не выходит поставаить физическую сущность (точку)
    Вот такой скрипт
    Merge "untitled.0-Windows\untitled";
    //Size
    a=1;
    //characteristic line
    c11=a/10;
    с12=Sqrt(2a^2)/100;
    c13=a/2;
    c14=a/5;
    c15=a/10;
    c16=a/100;
    Point(1) = {0, a, 0, c16};
    //+
    Point(2) = {a, 0, 0, c15};
    //+
    Point(3) = {0, 0, 0, c15};
    //+
    //Line_inter
    Line(1) = {1, 2};
    //+
    Line(2) = {2, 3};
    //+
    Line(3) = {3, 1};
    //+
    //Borders lines
    Line Loop(4) = {3, 1, 2};
    //+
    //surfase of triangle
    Plane Surface(5) = {4};


    Если приписать к нему
    Physical Poin(6) = {1};

    то файл сохранаяется в формате msh, но в самой структуре (если открыть через блокнот) становится только 1н узел и 1 элемент.
    Что я делаю не так?

    ОтветитьУдалить
    Ответы
    1. > Что я делаю не так?

      Если Вы хотите сохранить двумерную сетку. определите поверхность как физическую. См. также следующее:

      "строки 42-43: Определяются физические поверхности 1001 и 1002. По этим номерам мы всегда отыщем треугольники, принадлежащие разным областям и соответствующие материалам с различными физическими свойствами. *Кроме того, никакой другой элемент сетки, кроме указанных как физические, в файл сетки .msh не попадет.*"

      Удалить