среда, 22 июня 2011 г.

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

Продолжаем знакомиться с одним из самых популярных триангуляторов современности - Gmsh'ем. В прошлый раз мы рассматривали способ задания геометрии простой двумерной модели в .geo файле с помощью встроенного языка Gmsh. Сегодня продолжим исследовать возможности этого языка на примере более сложной трехмерной геометрии. Двумерная модель, как вы помните, имела следующий вид:
Пусть двумерная модель является сечением некоторой трехмерной модели. Конечно, таких моделей можно придумать несколько. Например:
1. Сфера в параллелепипеде
2. Цилиндр в параллелепипеде (при этом цилиндр может проходить через параллелепипед насквозь)
Однако, для того, чтобы рассмотреть некоторые особенности определения трехмерной геометрии, остановимся лишь на первом случае - сфера в параллелепипеде.
Кстати, обратите внимание на направление координатных осей, чтобы потом не запутаться (хотя на сферу с какой стороны ни смотри - все равно выглядит она одинаково):
Если вы не хотите каждый раз подолгу возить мышью, чтобы привести координатные оси в похожий вид, можете нажать Alt+Shift+Y и оси повернутся как нужно.

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

Способ 1
Вообще говоря, для того, чтобы задать эту модель "в лоб", почти ничего нового относительно описания двумерной геометрии делать не нужно. Просто больше писанины. Поэтому сразу приведу текст .geo файла (его можно скачать отсюда).

1    // номера физических объектов
2    LEFT_BOUNDARY = 3001;
3    RIGHT_BOUNDARY = 3002;
4    FRONT_BOUNDARY = 3003;
5    BACK_BOUNDARY = 3004;
6    BOTTOM_BOUNDARY = 3005;
7    TOP_BOUNDARY = 3006;
8    CUBE = 4001;
9    SPHERE = 4002;

10    // опции
11    Mesh.Algorithm = 6; // двумерный алгоритм - Фронтальный
12    Mesh.Algorithm3D = 4; // трехмерный алгоритм - Фронтальный
13    Mesh.OptimizeNetgen = 1; // включить Netgen-оптимизацию сетки

14    // параметры
15    clCube = 0.4; // характеристические длины
16    clSphere = 0.2;
17    rad = 0.25; // радиус сферы
18    X_BEG = 0; X_END = 1; // координаты границ параллелепипеда
19    Y_BEG = 0; Y_END = 1;
20    Z_BEG = 0; Z_END = 1;
21    lenX = X_END - X_BEG; // длины сторон параллелепипеда
22    lenY = Y_END - Y_BEG;
23    lenZ = Z_END - Z_BEG;
24    x_cen = (X_BEG + X_END) / 2.0; // координаты центра сферы
25    y_cen = (Y_BEG + Y_END) / 2.0;
26    z_cen = (Z_BEG + Z_END) / 2.0;
27    // наименьшая сторона параллелепипеда
28    min = (lenX < lenY ? lenX : lenY);
29    minSide = (lenZ < min ? lenZ : min);

30    // ---------- ТОЧКИ -----------

31    // параллелепипед
32    Point(1) = { X_BEG, Y_BEG, Z_BEG, clCube * minSide };
33    Point(2) = { X_END, Y_BEG, Z_BEG, clCube * minSide };
34    Point(3) = { X_BEG, Y_BEG, Z_END, clCube * minSide };
35    Point(4) = { X_END, Y_BEG, Z_END, clCube * minSide };
36    Point(5) = { X_BEG, Y_END, Z_BEG, clCube * minSide };
37    Point(6) = { X_END, Y_END, Z_BEG, clCube * minSide };
38    Point(7) = { X_BEG, Y_END, Z_END, clCube * minSide };
39    Point(8) = { X_END, Y_END, Z_END, clCube * minSide };
40    // сфера
41    Point(11) = { x_cen, y_cen, z_cen, clSphere * rad };
42    Point(12) = { x_cen, y_cen, z_cen + rad, clSphere * rad };
43    Point(13) = { x_cen - rad, y_cen, z_cen, clSphere * rad };
44    Point(14) = { x_cen, y_cen, z_cen - rad, clSphere * rad };
45    Point(15) = { x_cen + rad, y_cen, z_cen, clSphere * rad };
46    Point(16) = { x_cen, y_cen - rad, z_cen, clSphere * rad };
47    Point(17) = { x_cen, y_cen + rad, z_cen, clSphere * rad };

48    // ---------- ЛИНИИ -----------

49    // параллелепипед
50    Line(1) = { 1, 2 };
51    Line(2) = { 1, 3 };
52    Line(3) = { 2, 4 };
53    Line(4) = { 3, 4 };
54    Line(5) = { 5, 6 };
55    Line(6) = { 5, 7 };
56    Line(7) = { 6, 8 };
57    Line(8) = { 7, 8 };
58    Line(9) = { 1, 5 };
59    Line(10) = { 2, 6 };
60    Line(11) = { 3, 7 };
61    Line(12) = { 4, 8 };
62    // сфера
63    Circle(21) = { 12, 11, 13 };
64    Circle(22) = { 13, 11, 14 };
65    Circle(23) = { 14, 11, 15 };
66    Circle(24) = { 15, 11, 12 };
67    Circle(25) = { 13, 11, 16 };
68    Circle(26) = { 16, 11, 15 };
69    Circle(27) = { 15, 11, 17 };
70    Circle(28) = { 17, 11, 13 };
71    Circle(29) = { 12, 11, 16 };
72    Circle(30) = { 16, 11, 14 };
73    Circle(31) = { 14, 11, 17 };
74    Circle(32) = { 17, 11, 12 };

75    // ---------- КОНТУРЫ ЛИНИЙ -----------

76    // параллелепипед
77    Line Loop(41) = { 2, 11, -6, -9 };
78    Line Loop(42) = { 3, 12, -7, -10 };
79    Line Loop(43) = { 1, 3, -4, -2 };
80    Line Loop(44) = { 5, 7, -8, -6 };
81    Line Loop(45) = { 1, 10, -5, -9 };
82    Line Loop(46) = { 4, 12, -8, -11 };
83    // сфера
84    Line Loop(51) = { 21, 25, -29 };
85    Line Loop(52) = { 24, 29, 26 };
86    Line Loop(53) = { 22, -30, -25 };
87    Line Loop(54) = { 23, -26, 30 };
88    Line Loop(55) = { 21, -28, 32 };
89    Line Loop(56) = { 24, -32, -27 };
90    Line Loop(57) = { 22, 31, 28 };
91    Line Loop(58) = { 23, 27, -31 };

92    // ---------- ПОВЕРХНОСТИ -----------

93    // параллелепипед
94    Plane Surface(1) = { 41 };
95    Plane Surface(2) = { 42 };
96    Plane Surface(3) = { 43 };
97    Plane Surface(4) = { 44 };
98    Plane Surface(5) = { 45 };
99    Plane Surface(6) = { 46 };
100    // сфера
101    Ruled Surface(7) = { 51 };
102    Ruled Surface(8) = { 52 };
103    Ruled Surface(9) = { 53 };
104    Ruled Surface(10) = { 54 };
105    Ruled Surface(11) = { 55 };
106    Ruled Surface(12) = { 56 };
107    Ruled Surface(13) = { 57 };
108    Ruled Surface(14) = { 58 };

109    // ---------- КОНТУРЫ ПОВЕРХНОСТЕЙ -----------

110    // параллелепипед
111    Surface Loop(21) = { 1, 2, 3, 4, 5, 6 };
112    // сфера
113    Surface Loop(22) = { 7, 8, 9, 10, 11, 12, 13, 14 };

114    // ---------- ОБЪЕМЫ -----------

115    // параллелепипед без сферы
116    Volume(1) = { 21, 22 };
117    // сфера
118    Volume(2) = { 22 };

119    // ---------- ФИЗИЧЕСКИЕ СУЩНОСТИ -----------

120    Physical Surface(LEFT_BOUNDARY) = { 1 };
121    Physical Surface(RIGHT_BOUNDARY) = { 2 };
122    Physical Surface(FRONT_BOUNDARY) = { 3 };
123    Physical Surface(BACK_BOUNDARY) = { 4 };
124    Physical Surface(BOTTOM_BOUNDARY) = { 5 };
125    Physical Surface(TOP_BOUNDARY) = { 6 };
126    Physical Volume(CUBE) = { 1 };
127    Physical Volume(SPHERE) = { 2 };

Разберем некоторые интересные моменты.
строки 1-9: В самом начале .geo файла мы определяем номера физических объектов, которые мы хотим получить в результате в .msh файле. Кстати говоря, если убрать эти строки, а в объявлении физических сущностей названия объектов поставить в кавычки:

Physical Surface("LEFT_BOUNDARY") = { 1 };
Physical Surface("RIGHT_BOUNDARY") = { 2 };
Physical Surface("FRONT_BOUNDARY") = { 3 };
Physical Surface("BACK_BOUNDARY") = { 4 };
Physical Surface("BOTTOM_BOUNDARY") = { 5 };
Physical Surface("TOP_BOUNDARY") = { 6 };
Physical Volume("CUBE") = { 1 };
Physical Volume("SPHERE") = { 2 };

то в .msh файле вы получите дополнительную секцию вот такого вида:

$PhysicalNames
8
2 1 "LEFT_BOUNDARY"
2 2 "RIGHT_BOUNDARY"
2 3 "FRONT_BOUNDARY"
2 4 "BACK_BOUNDARY"
2 5 "BOTTOM_BOUNDARY"
2 6 "TOP_BOUNDARY"
3 7 "CUBE"
3 8 "SPHERE"
$EndPhysicalNames

То есть Gmsh сам присвоит уникальный номер каждому физическому объекту (как вы догадались, этот номер - второе число в строке. Первое число - размерность геометрического элемента: 2 - для треугольников, 3 - для тетраэдров).
строки 10-13: В .geo файле очень удобно определять опции Gmsh'а. Список настроек, которые можно изменить, очень внушительный. Я даже планирую как-нибудь остановиться на нем подробнее. А пока рассмотрим такие опции как алгоритмы построения сетки. Для построения двумерной и трехмерной сеток можно использовать разные алгоритмы. Вот их список с номерами, которые можно задавать в .geo файле:
2D: Mesh.Algorithm = 1 (Адаптивный - по умолчанию), 5 (Делоне), 6 (Фронтальный)
3D: Mesh.Algorithm3D = 1 (Делоне - по умолчанию), 4 (Фронтальный)
Можно также оптимизировать сетку для повышения качества элементов:
Mesh.Optimize = 1 (по умолчанию - 0)
Mesh.OptimizeNetgen = 1 (по умолчанию - 0)
строки 14-108: Здесь все аналогично заданию двумерной модели. Единственное отличие - использование изогнутой поверхности Ruled Surface для определения частей сферы.
строки 109-113: Здесь определяются замкнутые контуры поверхностей. С этим мы не сталкивались в двумерном случае. Однако этот этап определения модели является несложным. Для определения контура поверхностей, независимо от того плоские ли это поверхности или изогнутые, используется команда Surface Loop(номер) = {последовательность поверхностей}; Номер контура должен быть уникальным в рамках нумерации поверхностей. Gmsh "не смотрит" на ориентацию поверхностей, поэтому их можно задавать в любой последовательности (в отличие от контура линий). Единственное условие - чтобы получившийся контур был замкнутым.
строки 114-118: Определение объемов осуществляется командой Volume(номер объема) = {основной контур поверхностей, вложенные контуры, ...}; Как и в случае с объявлением поверхностей, вложенные объемы можно "вырезать" из основного, используя их контуры. Номер объема должен быть уникальным в рамках нумерации объемов.
строки 119-127: Определение физических сущностей, я думаю, вопросов вызвать уже не должно.

Способ 2
"Продвинутый" способ задания модели подразумевает использование, помимо основных элементов геометрии (точка, линия, поверхность и т.д.), также функций, которые во многом упрощают код .geo файла. Начнем, как и прежде, с текста Gmsh-скрипта (скачать его можно здесь):

1    // номера физических объектов
2    LEFT_BOUNDARY = 3001;
3    RIGHT_BOUNDARY = 3002;
4    FRONT_BOUNDARY = 3003;
5    BACK_BOUNDARY = 3004;
6    BOTTOM_BOUNDARY = 3005;
7    TOP_BOUNDARY = 3006;
8    CUBE = 4001;
9    SPHERE = 4002;

10    // опции
11    Mesh.Algorithm = 6; // двумерный алгоритм - Фронтальный
12    Mesh.Algorithm3D = 4; // трехмерный алгоритм - Фронтальный
13    Mesh.OptimizeNetgen = 1; // включить Netgen-оптимизацию сетки
14    Geometry.ExtrudeReturnLateralEntities = 0;

15    // параметры
16    clCube = 0.4; // характеристические длины
17    clSphere = 0.2;
18    rad = 0.25; // радиус сферы
19    X_BEG = 0; X_END = 1; // координаты границ параллелепипеда
20    Y_BEG = 0; Y_END = 1;
21    Z_BEG = 0; Z_END = 1;
22    lenX = X_END - X_BEG; // длины сторон параллелепипеда
23    lenY = Y_END - Y_BEG;
24    lenZ = Z_END - Z_BEG;
25    x_cen = (X_BEG + X_END) / 2.0; // координаты центра сферы
26    y_cen = (Y_BEG + Y_END) / 2.0;
27    z_cen = (Z_BEG + Z_END) / 2.0;
28    // наименьшая сторона параллелепипеда
29    min = (lenX < lenY ? lenX : lenY);
30    minSide = (lenZ < min ? lenZ : min);

31    // ---------- ТОЧКИ -----------

32    // прямоугольник
33    Point(1) = { X_BEG, Y_BEG, Z_BEG, clCube * minSide };
34    Point(2) = { X_END, Y_BEG, Z_BEG, clCube * minSide };
35    Point(3) = { X_BEG, Y_BEG, Z_END, clCube * minSide };
36    Point(4) = { X_END, Y_BEG, Z_END, clCube * minSide };
37    // полуокружность
38    Point(11) = { x_cen, y_cen, z_cen, clSphere * rad };
39    Point(12) = { x_cen - rad, y_cen, z_cen, clSphere * rad };
40    Point(13) = { x_cen, y_cen, z_cen - rad, clSphere * rad };
41    Point(14) = { x_cen + rad, y_cen, z_cen, clSphere * rad };

42    // ---------- ЛИНИИ -----------

43    // прямоугольник
44    Line(1) = { 1, 2 };
45    Line(2) = { 1, 3 };
46    Line(3) = { 2, 4 };
47    Line(4) = { 3, 4 };
48    // полуокружность
49    Circle(21) = { 12, 11, 13 };
50    Circle(22) = { 13, 11, 14 };

51    // ---------- КОНТУРЫ ЛИНИЙ -----------

52    // прямоугольник
53    Line Loop(41) = { 1, 3, -4, -2 };

54    // ---------- ПОВЕРХНОСТИ -----------

55    // прямоугольник
56    Plane Surface(1) = { 41 };

57    // ---------- ФУНКЦИИ -----------

58    // поверхности параллелепипеда
59    s1[] = Extrude { 0, Y_END, 0 } { Line{1, 2, 3, 4}; };
60    s2[] = Translate { 0, Y_END, 0 } { Duplicata{ Surface{1}; } };

61    // поверхности сферы
62    sph1[] = Extrude { { 1, 0, 0 }, { x_cen, y_cen, z_cen }, 2 * Pi / 3 } { Line{21, 22}; };
63    sph2[] = Extrude { { 1, 0, 0 }, { x_cen, y_cen, z_cen }, 2 * Pi / 3 } { Line{59, 62}; };
64    sph3[] = Extrude { { 1, 0, 0 }, { x_cen, y_cen, z_cen }, 2 * Pi / 3 } { Line{65, 68}; };

65    // ---------- КОНТУРЫ ПОВЕРХНОСТЕЙ -----------

66    // параллелепипед
67    Surface Loop(21) = { 1, 45, 49, 53, 57, 58 };
68    // сфера
69    Surface Loop(22) = { 61, 64, 67, 70, 73, 76 };

70    // ---------- ОБЪЕМЫ -----------

71    // параллелепипед без сферы
72    Volume(1) = { 21, 22 };
73    // сфера
74    Volume(2) = { 22 };

75    // ---------- ФИЗИЧЕСКИЕ СУЩНОСТИ -----------

76    Physical Surface(LEFT_BOUNDARY) = { 49 };
77    Physical Surface(RIGHT_BOUNDARY) = { 53 };
78    Physical Surface(FRONT_BOUNDARY) = { 1 };
79    Physical Surface(BACK_BOUNDARY) = { 58 };
80    Physical Surface(BOTTOM_BOUNDARY) = { 45 };
81    Physical Surface(TOP_BOUNDARY) = { 57 };
82    Physical Volume(CUBE) = { 1 };
83    Physical Volume(SPHERE) = { 2 };

Обратите внимание, что текст .geo файла сократился на 44 строки! Неплохо для такой незамысловатой модели. Теперь разберем ключевые моменты.
строка 14: В блоке опций появилась новая строка Geometry.ExtrudeReturnLateralEntities = 0; Ее смысл поясним чуть позднее.
строки 31-41: Теперь мы не задаем все точки модели. Для параллелепипеда нам достаточно определить точки только одной грани. А для сферы - одной полуокружности.
строки 42-50: Аналогично дело обстоит и с линиями.
строки 51-56: Что касается контуров линий и поверхностей, то в явном виде нам достаточно задать лишь одну грань параллелепипеда (на самом деле - и этого можно не делать, однако слишком перегружать .geo файл функциями Gmsh'а я не советую. Там, где несложно задать элементы геометрии явно - лучше сделать это).
строка 59: Ну а теперь поговорим непосредственно о функциях Gmsh'а для задания геометрии. В данном случае использовалась команда Extrude, что означает "вытягивание" или "выдавливание". То есть, вытягивая точку, вы получите линию, вытягивая линию - поверхность, вытягивая поверхность - объем. У этой команды 3 формата:
1. Extrude{X, Y, Z} {список фигур}
Это вытягивание вдоль вектора, заданного координатами X, Y, Z. Список фигур для вытягивания можно задать так: Point{1}; Line{15, 16}; Surface{34, 35, 36}; Различные типы элементов разделяются точкой с запятой. Обратите внимание, точка с запятой после команды Extrude не ставится. Однако большинство функций Gmsh'а (если не все) над элементами геометрии возвращают номера новых элементов, полученных в результате действия функции. И если вы хотите сохранить возвращаемый массив, тогда точку с запятой после команды придется поставить. Например, out[] = Extrude{1, 0, 0} {Line{1};};
2. Extrude{{X1, Y1, Z1}, {X2, Y2, Z2}, Angle} { список фигур }
Это вытягивание с поворотом. X1, Y1, Z1 - координаты линии, вокруг которой осуществляется вращение. X2, Y2, Z2 - координаты точки на этой линии. Angle - угол поворота в радианах (делайте угол поворота не больше Pi - кстати, в Gmsh уже есть такая встроенная константа). Список фигур задается аналогичным первому случаю образом.
3. Extrude{{X1, Y1, Z1}, {X2, Y2, Z2}, {X3, Y3, Z3}, Angle} { список фигур }
Вытягивание вдоль линии с поворотом. X1, Y1, Z1 - координаты вектора, вдоль которого осуществляется перемещение. Остальные переменные имеют такой же смысл что и во втором случае.
Посмотрим теперь, что мы задали в .geo файле:
s1[] = Extrude { 0, Y_END, 0 } { Line{1, 2, 3, 4}; };
Таким образом, мы использовали первый тип вытягивания - вдоль линии. Поскольку для выбранной ориентации координатных осей (см. выше) мы задали переднюю грань параллелепипеда, вытягиваем 4 линии вдоль оси Y на заданную длину. В результате мы получим 4 боковых грани. О том, какие значения получил массив s1, смотри ниже.
строка 60: Команда Translate используется для перемещения элементов геометрии на новое место. Формат ее таков: Translate{X, Y, Z} {список фигур}. X, Y, Z - координаты вектора перемещения. Ситуация с точкой с запятой после команды обстоит здесь абсолютно так же, как и с командой Extrude. В список фигур для Translate помимо Point, Line и Surface также можно поставить и объемы Volume. Если в списке фигур перед перемещаемым элементом поставить команду Duplicata{элементы}, тогда переместится не сам элемент, а его копия (дубликат). Таким образом, команда Translate{X, Y, Z} {Duplacata{фигура1}; Duplicata{фигура2};} попросту копирует нужные элементы и ставит их на нужное место. В нашем случае скопирована передняя грань параллелепипеда и получена задняя.
строки 62-64: Определяем поверхности сферы, "вытягивая" четверти окружности. В данном случае используется команда Extrude в своем втором формате - вытягивание вращением. При этом вращаются 2 линии, каждая из которых представляет собой четверть окружности. Для получения полной поверхности сферы используется 3 вращения на угол 120 градусов (как я уже говорил, угол поворота лучше брать меньше Pi). Однако, я уверен, больше всего вас в данный момент тревожат номера линий, которые я использовал в строках 63 и 64. Это номера линий, ограничивающих поверхности, полученные в результате предыдущего вытягивания. Чтобы их узнать, можно, например, вывести значения массива sph1[] в консоль сообщений Gmsh'а (Tools->Message Console или Ctrl+L). Делается это так: Printf("%g %g %g %g", sph1[0], sph1[1], sph1[2], sph1[3]); Для данного примера вы получите: 59 61 62 64 Однако сами по себе номера скажут вам немного. Далее, нужно будет зайти в Tools -> Visibility (или Shift+Ctrl+V) и на вкладке Elementary entities посмотреть, что эти номера обозначают:
Таким образом, для следующего вращения будут использоваться линии 59 и 62. Аналогично поступаем со следующим вытягиванием.
строки 66-83: В оставшихся строках происходит объявление контуров поверхностей, объемов и физических объектов. В принципе - ничего нового. Однако здесь также встречаются номера поверхностей, полученных вытягиванием линий. Эти номера, как и ранее, можно посмотреть в Gmsh'е. Однако, как вы понимаете, это не совсем правильный подход. Все-таки, мы боремся за универсальность задания нашей модели. Поэтому, если вдруг, Gmsh сгенерирует для вытянутой поверхности новый порядковый номер, его нужно будет отловить и переписать соответствующие места .geo файла. Конечно, Gmsh выбирает номера не случайным образом, и чтобы номер изменился, код файла геометрии должен претерпеть значительные изменения. Однако всегда хочется сделать свой код независимым от дальнейших изменений. Поэтому покопаемся в структуре массивов, возвращаемых функциями Gmsh'а. По умолчанию команда
out[] = Extrude{0, 0, 1} {Line{15};}; // или другие 2 ее вариации
возвращает 4 числа:
1. номер линии, на которой заканчивается перемещение заданной линии 15
2. номер полученной поверхности
3 и 4. номера линий, которые ограничивают полученную поверхность "с боков"
Приведем простой наглядный пример.

Point(1) = {0, 0, 0};
Point(2) = {1, 0, 0};
Line(1) = {1, 2};
out[] = Extrude{0, 1, 0}{Line{1};};
Printf("%g %g %g %g", out[0], out[1], out[2], out[3]);

В консоли сообщений мы видим: 2 5 4 -3, а в окне Gmsh'а:
Таким образом, для того, чтобы в дальнейшем использовать номер вытянутой поверхности, вместо ее явного номера 5 можно пользоваться out[1]. Если вам мешают номера боковых линий, их можно отключить командой Geometry.ExtrudeReturnLateralEntities = 0; В этом случае массив out[] будет состоять из 2-х чисел: 2 5, т.е. номера линии и номера поверхности.
Теперь, зная это, перепишем часть кода .geo файла для "продвинутого" способа (код целиком можно скачать здесь):

57    // ---------- ФУНКЦИИ -----------

58    // поверхности параллелепипеда
59    s1[] = Extrude { 0, Y_END, 0 } { Line{1, 2, 3, 4}; };
60    s2[] = Translate { 0, Y_END, 0 } { Duplicata{ Surface{1}; } };
61    // поверхности сферы
62    sph1[] = Extrude { { 1, 0, 0 }, { x_cen, y_cen, z_cen }, 2 * Pi / 3 } { Line{21, 22}; };
63    sph2[] = Extrude { { 1, 0, 0 }, { x_cen, y_cen, z_cen }, 2 * Pi / 3 } { Line{sph1[0], sph1[2]}; };
64    sph3[] = Extrude { { 1, 0, 0 }, { x_cen, y_cen, z_cen }, 2 * Pi / 3 } { Line{sph2[0], sph2[2]}; };

65    // ---------- КОНТУРЫ ПОВЕРХНОСТЕЙ -----------

66    // параллелепипед
67    Surface Loop(21) = { 1, s1[1], s1[3], s1[5], s1[7], s2[0] };
68    // сфера
69    Surface Loop(22) = { sph1[1], sph1[3], sph2[1], sph2[3], sph3[1], sph3[3] };

70    // ---------- ОБЪЕМЫ -----------

71    // параллелепипед без сферы
72    Volume(1) = { 21, 22 };
73    // сфера
74    Volume(2) = { 22 };

75    // ---------- ФИЗИЧЕСКИЕ СУЩНОСТИ -----------

76    Physical Surface(LEFT_BOUNDARY) = { s1[3] };
77    Physical Surface(RIGHT_BOUNDARY) = { s1[5] };
78    Physical Surface(FRONT_BOUNDARY) = { 1 };
79    Physical Surface(BACK_BOUNDARY) = { s2[0] };
80    Physical Surface(BOTTOM_BOUNDARY) = { s1[1] };
81    Physical Surface(TOP_BOUNDARY) = { s1[7] };
82    Physical Volume(CUBE) = { 1 };
83    Physical Volume(SPHERE) = { 2 };

строка 67 и 79: Команда Translate возвращает массив, в котором перемещенный элемент стоит на нулевой позиции.

Я надеюсь, теперь основные моменты задания трехмерной геометрии стали понятны.
Для того, чтобы построить саму сетку, теперь вам достаточно нажать кнопку 3D в разделе Mesh, как показано на рисунке:
Кстати, напомню, что все рисунки, в том числе и этот вид меню, актуальны для версии Gmsh 2.4.2. В последней версии (на данный момент это 2.5.0) вид немного другой.

В следующий раз мы рассмотрим некоторые элементы управления, которые предоставляет Gmsh для работы с сеткой.

16 комментариев:

  1. Можно ли какой-либо командой сделать копии объекта множество раз в заданные координаты?

    ОтветитьУдалить
  2. Допустим, в приведенном примере (сфера в параллелепипеде) нужно разместить несколько сфер. Как это сделать, скопировав исходную таким образом, чтобы новые также были отдельными телами в параллелепипеде?

    ОтветитьУдалить
  3. Сделать это можно довольно просто с помощью связки команд Duplicata и Translate.

    Формат команды такой:

    ret[] = Translate { X, Y, Z } { Duplicata { Point | Line | Surface | Volume { v1, v2, ... }; } };

    где X, Y, Z - координаты вектора перемещения, v1, v2 - номера объектов, которые надо скопировать и переместить.

    Пример со сферой, описанной в первом способе:
    Уменьшим радиус, скажем, до 0.05. Центр сферы в точке (0.5, 0.5, 0.5). Создадим такую же сферу с центром в точке (0.6, 0.6, 0.6) и "вырежем" ее из куба, чтобы получить согласованную сетку (привожу здесь часть кода, которая подверглась изменениям):

    // ...
    rad = 0.05;
    // ...

    // ---------- ПОВЕРХНОСТИ -----------

    // параллелепипед
    // ...

    // сфера
    // ...

    newSph[] = Translate { 0.1, 0.1, 0.1 } { Duplicata { Surface { 7, 8, 9, 10, 11, 12, 13, 14 }; } };

    // ---------- КОНТУРЫ ПОВЕРХНОСТЕЙ -----------

    // параллелепипед
    // ...
    // сфера
    // ...
    Surface Loop(23) = { newSph[] };

    // ---------- ОБЪЕМЫ -----------

    // параллелепипед без сфер
    Volume(1) = { 21, 22, 23 };
    // сферы
    Volume(2) = { 22 };
    Volume(3) = { 23 };

    // ---------- ФИЗИЧЕСКИЕ СУЩНОСТИ -----------

    // ...
    Physical Volume(SPHERE) = { 2, 3 };

    ОтветитьУдалить
  4. Заданы контуры поверхностей шара (44) в кубе (25). Затем из поверхностей (44) сделан соответствующий объем. После этого скопирован этот объем. Далее задан объем куба за исключением двух шаров (первоначального и скопированного).
    Вопрос: так как отображаются только объемы шаров, но не куба, что сделано не так?
    И можно ли сократить как-то код?
    ---------------------------------
    // Создание объема образца
    // Размеры образца
    a=15;
    b=15;
    c=15;
    // Координаты образца
    x0 = 0;
    y0 = 0;
    z0 = 0;
    x = 1*a;
    y = 1*b;
    z = 1*c;
    cl1 = c/20;
    // Построение образца
    Point(1) = {x0 , y0 , z0, cl1};
    Point(2) = {x, y0, z0, cl1};
    Point(3) = {x, y, z0, cl1};
    Point(4) = {x0, y, z0, cl1};
    Point(5) = {x0, y, z, cl1};
    Point(6) = {x, y, z, cl1};
    Point(7) = {x, y0, z, cl1};
    Point(8) = {x0, y0, z, cl1};
    Line(1) = {1, 2};
    Line(2) = {2, 3};
    Line(3) = {3, 4};
    Line(4) = {4, 1};
    Line(5) = {2, 7};
    Line(6) = {7, 6};
    Line(7) = {6, 3};
    Line(8) = {7, 8};
    Line(9) = {8, 5};
    Line(10) = {5, 6};
    Line(11) = {4, 5};
    Line(12) = {1, 8};
    Line Loop(13) = {1, 2, 3, 4};
    Line Loop(14) = {5, 6, 7, -2};
    Plane Surface(15) = {14};
    Line Loop(16) = {8, 9, 10, -6};
    Plane Surface(17) = {16};
    Line Loop(18) = {8, -12, 1, 5};
    Plane Surface(19) = {18};
    Line Loop(20) = {12, 9, -11, 4};
    Plane Surface(21) = {20};
    Line Loop(22) = {11, 10, 7, 3};
    Plane Surface(23) = {22};
    Plane Surface(24) = {13};

    // Создание пор
    rad = 2;
    cl2 = c/2;
    // Координаты ц.т. поры
    x_c1=5; y_c1=5; z_c1=5;
    Point(9) = {x_c1, y_c1, z_c1, cl2}; // точка ц.т. поры
    // Точки окружности поры
    Point(10) = {x_c1-rad, y_c1, z_c1, cl2};
    Point(11) = {x_c1+rad, y_c1, z_c1, cl2};
    Point(12) = {x_c1, y_c1-rad, z_c1, cl2};
    Point(13) = {x_c1, y_c1+rad, z_c1, cl2};
    Point(14) = {x_c1, y_c1, z_c1-rad, cl2};
    Point(15) = {x_c1, y_c1, z_c1+rad, cl2};
    // Линии окружности поры
    Circle(16) = {10, 9, 12};
    Circle(17) = {12, 9, 11};
    Circle(18) = {11, 9, 13};
    Circle(19) = {13, 9, 10};
    Circle(20) = {12, 9, 14};
    Circle(21) = {14, 9, 13};
    Circle(22) = {13, 9, 15};
    Circle(23) = {15, 9, 12};
    Circle(24) = {11, 9, 14};
    Circle(25) = {14, 9, 10};
    Circle(26) = {10, 9, 15};
    Circle(27) = {15, 9, 11};
    // замкнутый контур окружности
    Line Loop(28) = {22, -26, -19};
    Plane Surface(29) = {28};
    Line Loop(30) = {26, 23, -16};
    Plane Surface(31) = {30};
    Line Loop(32) = {22, 27, 18};
    Plane Surface(33) = {32};
    Line Loop(34) = {23, 17, -27};
    Plane Surface(35) = {34};
    Line Loop(36) = {25, -19, -21};
    Plane Surface(37) = {36};
    Line Loop(38) = {18, -21, -24};
    Plane Surface(39) = {38};
    Line Loop(40) = {25, 16, 20};
    Plane Surface(41) = {40};
    Line Loop(42) = {17, 24, -20};
    Plane Surface(43) = {42};

    Surface Loop(25) = {24, 19, 17, 21, 23, 15};
    Surface Loop(44) = {29, 31, 33, 35, 37, 39, 41, 43};

    Volume(45) = {44};
    Translate {5, 5, 5} {Duplicata{Volume{45};}}
    Volume(47) = {25, 44};

    ОтветитьУдалить
  5. >>отображаются только объемы шаров, но не куба

    Не понял, что вы имеете ввиду.
    Запустил ваш код. Вывел номера объемов (Tools -> Options -> Geometry -> Visibility -> Volume numbers). Отображаются номера 45, 46, 47.

    По поводу короткого кода: короткий - не значит, простой. А в случае с Gmsh'ем зачастую - наоборот. Но вы, конечно, можете укоротить код, используя команду Extrude для создания линий, поверхностей и объема куба (задав лишь одну точку), и команду Symmetry для создания шара (задав лишь одну восьмую его часть) или Extrude с поворотом (задав тоже только одну точку, ну или, может, две). Вариантов много.

    Теперь конструктивная критика относительно вашего кода. Посмотрите как выглядит дискретизация ваших шаров. У меня (gmsh 2.4.2) она просто чудовищна. А все потому, что вы используете Plane Surface вместо Ruled Surface для задания поверхностей шара. Смотрите внимательнее готовые примеры.
    Теперь такой вопрос. Вы создали новый шар, скопировав объем. Все здорово, но как вы собираетесь строить согласованную сетку? Или ваш метод решения не требует этого (разрывный Галеркин и др)? Чтобы сетка была согласована, нужно "вырезать" объем шара из объема куба. А вырезать можно только контуры поверхностей.
    И последний вопрос - чем вам не угодил мой пример из комментов? :) По-моему, ваш код делает то же самое.

    ОтветитьУдалить
  6. >>отображаются только объемы шаров, но не куба

    Ага, теперь понял. В кубе сетка не строится. Чтобы построилась, нужно "вырезать" вложенные объемы, используя контуры поверхностей. Это относится к скопированному шару. Тогда все будет нормально.

    ОтветитьУдалить
  7. >> А все потому, что вы используете Plane Surface вместо Ruled Surface для
    задания поверхностей шара. Смотрите внимательнее готовые примеры.

    Спасибо! Я попробую.

    >> Все здорово, но как вы собираетесь строить согласованную сетку? Или ваш метод решения не требует этого (разрывный Галеркин и др)? Чтобы сетка была согласована, нужно "вырезать" объем шара из объема куба. А вырезать можно только контуры поверхностей.

    Требует. Сетка должна быть согласована. А что сделано не так?

    >> И последний вопрос - чем вам не угодил мой пример из комментов? :) По-моему, ваш код делает то же самое.

    Он угодил. Просто я раньше потихоньку строить уже начал. Не стал переделывать.

    >> ...нужно "вырезать" вложенные объемы, используя контуры поверхностей. Это относится к скопированному шару. Тогда все будет нормально...

    То есть, если я правильно понял, последняя строка должна иметь вид:
    Volume(47) = {25, 44, 46};?
    Но программа выдает ошибку: Unknown surfaceloop 46.
    Я думаю, что как-то надо скопированный объем обозначить? Правильно? Подскажите, пожалуйста!


    И такой вопрос в тему:
    Я это пробую пока на примере одного-двух шаров. В окончательном виде у меня в кубе должно быть несколько тысяч шаров. Я готовлю подпрограммку, которая выдаст мне все их координаты. Как Вы думаете, приведенный код для двух шаров подойдет? Т.е., сначала задается один шар, а потом множество раз копируется в заданные координаты.
    Если такое возможно, то какой самый короткий способ это сделать?

    Заранее спасибо!

    ОтветитьУдалить
  8. >> ...нужно "вырезать" вложенные объемы, используя контуры поверхностей. Это относится к скопированному шару. Тогда все будет нормально...

    А как узнать поверхности скопированного шара? И если их будет несколько тысяч, каким образом можно их будет вырезать?

    ОтветитьУдалить
  9. >> То есть, если я правильно понял,
    >> последняя строка должна иметь вид:
    >> Volume(47) = {25, 44, 46};

    Нет. Вырезать можно только контуры поверхностей. Т.е. 25 - это контур поверхностей куба, 44 - контур поверхности первого шара, а 46 - это объем. Поэтому так делать нельзя.
    Посмотрите мой пример из комментов - там копируются не объемы, а поверхности исходного шара. Затем из скопированных поверхностей создается контур, который вырезается из контура куба. Таким образом, вы получаете согласованную сетку.
    Можно придумать и другой способ.
    Вы скопировали объем. Но если посмотреть в массив, который возвращается командой Translate, можно заметить, что номера скопированных поверхностей он не содержит. Чтобы их узнать можно воспользоваться командой Boundary.
    Формат такой:

    surf[] = Boundary { Volume {nv }; };

    nv - номер скопированного объема ( в вашем случае это 46 )
    массив surf[] содержит номера скопированных поверхностей

    После этого организуйте контур поверхностей нового шара:

    Surface Loop(123) = { surf[] };

    Теперь его можно вырезать из контура поверхностей куба:

    Volume(47) = {25, 44, 123};

    В итоге сетка также согласована.

    Чтобы задать много шаров, конечно, проще всего их все копировать из одного. Для простоты можно задать шар с центром в начале координат (в этом случае вектор перемещения в команде Translate совпадает с координатами центра нового шара), а потом исходный шар удалить (см. команду Delete).
    Я, например, все этого себе представляю в виде такого цикла:

    For i In {1:N} // N - количество шаров
    // копирование шара в заданные координаты
    EndIf

    ОтветитьУдалить
  10. Спасибо! Попробую.
    Как Вы считаете, по силам GMSH разбить сетку такого кубика с количеством шариков под 3000? Они будут размерами в 2 раза меньше.

    ОтветитьУдалить
  11. Сложно сказать. Все зависит от технических возможностей. Пока не запустите, не узнаете.

    ОтветитьУдалить
  12. >> For i In {1:N} // N - количество шаров
    // копирование шара в заданные координаты
    EndIf

    А можно, пожалуйста, подробнее? Как скопировать шар множество раз? Где и на каком этапе я должен задать координаты? Просто с циклами не сталкивался...

    Т.е., первая строчка:
    For i In {1:N} // N - количество шаров
    понятна. Задается некоторый массив с определенным количеством шаров

    Дальше:
    // копирование шара в заданные координаты
    Каким образом можно скопировать шар множество раз?
    Сколько шаров, столько и строчек с командой Translate задавать?
    Или можно в одной команде Translate как-то координаты местоположения всех шаров указать?
    Или как-то координаты изначально указываются, а командой Translate в цикле как-то вызываются потом?

    В конце:
    EndIf
    Это, как я понимаю, конец цикла?

    ОтветитьУдалить
  13. >> В конце:
    >> EndIf
    >> Это, как я понимаю, конец цикла?

    Опечатка. Должно быть EndFor.
    По поводу всего остального - я, наверное, попозже пример какой-нибудь сооружу.

    ОтветитьУдалить
  14. Спасибо! Был бы очень признателен!

    ОтветитьУдалить
  15. Пример опубликован отдельным постом.

    ОтветитьУдалить