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 для работы с сеткой.