Gmsh. Compound

Появились интересные вопросы по командам Compound Line | Surface | Volume. Все это совпало с моей собственной попыткой (сразу скажу - неудачной) использовать эти команды для одного хитрого задания геометрии. Таким образом, накопился небольшой, но вполне достаточный для поста объем информации, который решено было выложить в блог.

Команды Compound довольно новые. Если верить истории версий Gmsh, они появились в версии 2.5.0. Поэтому любители версии 2.4.2, к которым не так давно относился и я, ими воспользоваться не смогут.

Начнем с описания этих комад, приведенного в документации:

Compound Line ( expression ) = { expression-list };
Creates a compound line from several elementary lines. When meshed, a compound line will be reparametrized as a single line, whose mesh can thus cross internal boundaries. The expression inside the parentheses is the compound line's identification number; the expression-list on the right hand side contains the identification number of the elementary lines that should be reparametrized as a single line.

Compound Surface ( expression ) = { expression-list }
      < Boundary { { expression-list }, { expression-list }, { expression-list }, { expression-list } } > ;
Creates a compound surface from several elementary surfaces. When meshed, a compound surface will be reparametrized as a single surface, whose mesh can thus cross internal boundaries. The expression inside the parentheses is the compound surface's identification number; the mandatory expression-list on the right hand side contains the identification number of the elementary surfaces that should be reparametrized as a single surface.

Compound Volume ( expression ) = { expression-list };
Creates a compound volume from several elementary volumes. When meshed, a compound volume will be reparametrized as a single volume, whose mesh can thus cross internal boundaries. The expression inside the parentheses is the compound volume's identification number; the expression-list on the right hand side contains the identification number of the elementary volumes that should be reparametrized as a single volume.

Как видно, в описании очень много общего, поэтому, абстрагируясь от слов Line | Surface | Volume, и заменив их одним словом "объект", сделаем перевод:

Составной Объект ( номер ) = { список номеров };
Создает составной объект из нескольких элементарных объектов. Во время построения сетки составной объект будет перепараметризован как один объект, сетка на котором может пересекать внутренние границы. Номер в скобках - это идентификационный номер составного объекта; список номеров в правой части содержит идентификационные номера элементарных объектов, которые должны быть перепараметризованы как один объект.

Я умышленно опустил часть команды Compound Surface, содержащую опциональную часть < Boundary... >, т.к. о ней даже в документации ничего не сказано.

Итак, начитавшись подобных прелестей, я бросился помогать коллеге решать мучившую его проблему. Проблема такая: необходимо задать Ruled поверхность, натянутую на более чем 4 линии. Как известно, для Plane поверхности не важно на скольки линиях строиться. Для Ruled - только 3 или 4. И никак иначе. У нас же геометрия имела следующий вид:

model1.geo

cl = 0.1; // характеристическая длина сетки

// p1[] - точки одной дуги
p1[0] = newp; Point(p1[0]) = { 0, 0, 0, cl };
// p2[] - точки другой дуги
p2[0] = newp; Point(p2[0]) = { 0, 3, 0, cl };

N = 10; // количество линий на дуге
For i In {0 : N - 1}
  ex[] = Extrude{ { 0, 1, 0 }, { 1, 0, 0 }, Pi/N } { Point{p1[i]}; };
  p1[i + 1] = ex[0];
  l1[i] = ex[1];
EndFor

ex[] = Extrude{ { 0, 1, 0 }, { 1, 0, 0 }, Pi } { Point{p2[0]}; };
p2[1] = ex[0];
l2[0] = ex[1];

// линии между дугами
l3[0] = newl; Line(l3[0]) = { p1[0], p2[0] };
l3[1] = newl; Line(l3[1]) = { p1[N], p2[1] };

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

ll1 = newll; // контур линий
Line Loop(ll1) = { l1[], l3[1], -l2[], -l3[0] };
s1 = news; // поверхность
Plane Surface(s1) = { ll1 };

мы получим, нечто такое

А написать так:

s1 = news; // поверхность
Ruled Surface(s1) = { ll1 };

нам не даст Gmsh:

Error : Wrong definition of Ruled Surface 15: 13 borders instead of 3 or 4

Проблема решается просто: необходимо задать дополнительные линии на нужной поверхности таким образом, чтобы получилось несколько Ruled поверхностей.

Например, соединив все заданные точки следующим образом:

Понятно, что ни к чему хорошему это не приведет:

Другая идея - разбить вторую дугу на такое же количество линий, что и первую:

model2.geo

cl = 0.1; // характеристическая длина сетки

// p1[] - точки одной дуги
p1[0] = newp; Point(p1[0]) = { 0, 0, 0, cl };
// p2[] - точки другой дуги
p2[0] = newp; Point(p2[0]) = { 0, 3, 0, cl };

N = 10; // количество линий на дуге
For i In {0 : N - 1}
  // одна дуга
  ex[] = Extrude{ { 0, 1, 0 }, { 1, 0, 0 }, Pi/N } { Point{p1[i]}; };
  p1[i + 1] = ex[0];
  l1[i] = ex[1];
  // другая дуга
  ex[] = Extrude{ { 0, 1, 0 }, { 1, 0, 0 }, Pi/N } { Point{p2[i]}; };
  p2[i + 1] = ex[0];
  l2[i] = ex[1];
EndFor

// линии между дугами
For i In {0 : N}
  l3[i] = newl;
  Line(l3[i]) = { p1[i], p2[i] };
EndFor

В этом случае ясно, что, определив Ruled поверхности:

For i In {0 : N - 1}
  ll = newll;
  Line Loop(ll) = { l1[i], l3[i + 1], -l2[i], -l3[i] };
  s[i] = news;
  Ruled Surface(s[i]) = { ll };
EndFor

мы получим, то, что хотели:

Но вот было бы хорошо не создавать новые точки и такое количество поверхностей, а обойтись исходным набором геометрических примитивов. Вот здесь и появилась идея использовать Compound Line для определения первой полуокружности как единой линии для того, чтобы затем с легкостью использовать ее для создания одной Ruled поверхности, натянутой на 4 линии.

Пробуем:

model3.geo

cl = 0.1; // характеристическая длина сетки

// p1[] - точки одной дуги
p1[0] = newp; Point(p1[0]) = { 0, 0, 0, cl };
// p2[] - точки другой дуги
p2[0] = newp; Point(p2[0]) = { 0, 3, 0, cl };

N = 10; // количество линий на дуге
For i In {0 : N - 1}
  ex[] = Extrude{ { 0, 1, 0 }, { 1, 0, 0 }, Pi/N } { Point{p1[i]}; };
  p1[i + 1] = ex[0];
  l1[i] = ex[1];
EndFor

ex[] = Extrude{ { 0, 1, 0 }, { 1, 0, 0 }, Pi } { Point{p2[0]}; };
p2[1] = ex[0];
l2[0] = ex[1];

// линии между дугами
l3[0] = newl; Line(l3[0]) = { p1[0], p2[0] };
l3[1] = newl; Line(l3[1]) = { p1[N], p2[1] };

cline = newl; // составная линия
Compound Line(cline) = { l1[] };

bound1[] = Boundary { Line{l1[0]}; };
bound2[] = Boundary { Line{l1[N - 1]}; };
bound3[] = Boundary { Line{l3[1]}; };
bound4[] = Boundary { Line{-l2[0]}; };
bound5[] = Boundary { Line{-l3[0]}; };

Printf("line loop:");
Printf("%f %f", bound1[0], bound2[1]);
Printf("%f %f", bound3[0], bound3[1]);
Printf("%f %f", bound4[0], bound4[1]);
Printf("%f %f", bound5[0], bound5[1]);

ll = newl; // контур линий
Line Loop(ll) = { cline, l3[1], -l2[0], -l3[0] };
s = news; // поверхность
Ruled Surface(s) = { ll };

Блок кода Boundary+Printf на геометрию не влияет и приведен только для того, чтобы показать, что контур линий, который мы хотим определить, действительно замкнутый. Смотрим в консоль сообщений Gmsh:

line loop:
1.000000 13.000000
13.000000 14.000000
14.000000 2.000000
2.000000 1.000000

Все нормально. К сожалению, команда Boundary неприменима к составным геометрическим элементам, т.е. написать Boundary{ Line{cline}; } не получилось, ибо Gmsh ругался так: "Error: Impossible to take boundary of entity 14 (of type 214)"

Вообще говоря, уже это должно насторожить, потому как, если Gmsh не может определить границы составной линии, то сможет ли он ее использовать далее для построения модели, т.е. в более сложных командах нежели Boundary? Видимо нет, потому что, несмотря на правильность контура, в результате обработки скрипта, приведенного ранее, Gmsh выдает:

Error : Line Loop 15 is wrong
Error : Wrong definition of Ruled Surface 16: 1 borders instead of 3 or 4

Я, на всякий случай, попробовал построить 2D сетку, Gmsh упал.

В общем, задумка не удалась. Таким образом, единственным способом задания сетки на поверхности, определенной изначально, оказался способ с добавлением новых точек (скрипт model2.geo).

Для чего же тогда нужны составные элементы. Обратимся еще раз к описанию:

"Во время построения сетки составной объект будет перепараметризован как один объект, сетка на котором может пересекать внутренние границы."

Попробуем применить это к нашему случаю. Дополним скрипт model2.geo следующими строками:

cline1 = newl;
Compound Line(cline1) = { l1[] };
csurf = news;
Compound Surface(csurf) = { s[] };

Мы определили составную линию вместо всех линий первой дуги и составную поверхность вместо всех Ruled поверхностей.

В этом случае геометрия изменится следующим образом (слева - без Compound):

Хотя, на самом деле, изменение внешнего вида модели определяется параметром Geometry.HideCompounds, который равен 1 по умолчанию. В этом случае те геометрические примитивы, которые образуют составной объект, не показываются. Если принудительно поставить Geometry.HideCompounds = 0; в тексте скрипта, то все линии и точки вернутся на место.

А вот что произойдет с поверхностной сеткой (слева - без Compound):

На правом рисунке присутствуют несколько сеток - одна для Compound Surface, другие для всех Ruled Surface, которые образуют составную. Чтобы этого избежать, можно полазить в Tools->Visibility и выбрать показ только составной поверхности, а можно переключить параметр Geometry.HideCompounds = 1;

Получаем следующее:

Увеличим изображение в районе первой полуокружности (верхний рисунок - без Compound):

Видно, что без Compound, сетка согласована с узлами, определенными на первой дуге. В свою очередь составная линия строит 1D сетку без учета геометрических узлов образующих ее линий. То же касается и составной поверхности, которая не учитывает линий Ruled поверхностей. Таким образом, видно "как сетка пересекает внутренние границы".

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

Рассмотрим следующий пример. Допустим, в приведенной выше модели имеется сгущение сетки к одному из узлов первой дуги, заданное уменьшенным значением характеристической длины. Например, так:

// изменяем характеристическую длину одной точки
Characteristic Length{ p1[N - 2] } = cl / 10;

После этого посмотрим как изменились поверхностные сетки (слева - без Compound):

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

Однако преимущества составных объектов все же есть, но они немного из другой области. Более подробно ознакомиться с ними позволяет статья "J.-F. Remacle, C. Geuzaine, G. Compere and E. Marchandise, High Quality Surface Remeshing Using Harmonic Maps, International Journal for Numerical Methods in Engineering, 2009", которая имеется в свободном доступе. Попробую кратко изложить основную идею. Понятно, что, если мы задаем геометрию в Gmsh'е, трудно придумать ситуацию, в которой Compound объекты имели бы большое преимущество над обычными. Например, в приведенном выше примере сгущение сетки может быть вызвано только какой-нибудь необходимостью. Составной объект это сгущение обойдет стороной. Однако Gmsh работает не только с моделями, заданными в .geo файле, но и с более сложными объектами, импортированными из CAD программ, а также полученными в ходе, например, сканирования. Многие объекты, для которых требуется получить качественную триангуляцию, представляются в виде STL (stereolithography) разбиения (поверхностное симплициальное разбиение). На основе поверхностной STL сетки Gmsh может построить трехмерную сетку, однако во многих случаях, STL разбиение имеет плохое качество. Так вот одна из целей внедрения Compound объектов - улучшить качество поверхностных STL сеток за счет перепараметризации геометрических объектов. Примером может служить скрипт t13.geo, приведенный в документации. Другой целью является улучшение сетки для моделей, созданных в CAD системах. Большинство таких систем представляет поверхности с виде NURBS (Non Uniform Rational B-Splines). Такое представление, вообще говоря, позволяет строить качественные сетки. Но бывают модели, части которых заданы отдельными геометрическими поверхностями, в то время как физически это один объект. Обычно сетку в этом случае строят для каждой геометрической поверхности отдельно. И это понятно. Однако если одна геометрическая поверхность имеет малые размеры, получается довольно сильное сгущение сетки, которое, вообще-то, никому не нужно. Это аналогично случаю, рассмотренному нами ранее, но только наоборот - когда от локального сгущения нужно как-то избавиться. В этом случае полезно несколько геометрических поверхностей объединить с помощью составного объекта, и сетка будет более равномерная и, соответственно, качественная.

В заключение можно отметить, что приемы с составными объектами, описанные в документации, работают не для каждой модели. По крайней мере, для трех STL файлов, методом тыка выбранных мной в интернете, ни одна из сеток не построилась. Если что-то получится, выложу отдельным постом.