Геометрия объекта со сложной структурой микровключений

Моделирование процессов в объектах со сложной структурой микровключений (сюда относятся композиты, геофизические среды, различные смеси, бетон и др), которое в последнее время становится весьма популярным (или оно давно стало популярным, а я заметил это только сейчас), побудило меня сконструировать какой-нибудь код, которым могли бы пользоваться для определения геометрии таких объектов. Конечно, у каждой задачи свои особенности и создать универсальный код вряд ли получится, но, может, какие-нибудь его части или идеи могут быть кому-нибудь полезны (а кому-то, кто моделирует бетон с кусочками щебня - точно пригодится ;-) ).
Начнем с простого. Зададим геометрию объекта, имеющего большое количество сферических непересекающихся микровключений. Слово "непересекающихся" здесь ключевое, потому что мы все будем делать в Gmsh'е, а он, как известно, не умеет осуществлять логические операции над объектами (объединение, пересечение и т.д.). Для этих целей нужно воспользоваться дополнительным инструментом, некоторые из них я перечислял в посте Триангуляция сложных моделей.
Многие приемы задания геометрии уже были рассмотрены, например, здесь. Некоторые мы рассмотрим впервые.

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

brick.geo

// начальная точка параллелепипеда
x0 = 0;
y0 = 0;
z0 = 0;
// размеры параллелепипеда
lenX = 1;
lenY = 1;
lenZ = 1;
// характеристическая длина элементов разбиения параллелепипеда
clBrick = lenX / 2;

// задаем точку
Point(1) = { x0, y0, z0, clBrick };

// получаем линию
lineBrick[] = Extrude { lenX, 0, 0 } { Point{1}; };

// получаем поверхность
surfBrick[] = Extrude { 0, lenY, 0 } { Line{lineBrick[1]}; };

// получаем объем
volBrick[]  = Extrude { 0, 0, lenZ } { Surface{surfBrick[1]}; };

// сохраняем поверхность параллелепипеда в виде массива
// для удобства дальнейшего использования
brickSurface[] = Boundary { Volume{volBrick[1]}; };

// сохраняем также номер замкнутого контура поверхностей
sl = newreg;
Surface Loop(sl) = { brickSurface[] };
brickSurfaceLoop = sl;

Теперь зададим шар, символизирующий микровключения. Оформим его в виде функции. Как и для параллелепипеда, создадим новый файл.

sphere.geo

// центр шара - (x_c, y_c, z_c)
// радиус шара - rad
// характеристическая длина элементов разбиения шара - clSphere
// порядковый номер шара - number

Function SphereTemplate
  // задаем 3 точки сектора сферы
  p1 = newp; // команда newp определяет новый уникальный номер точки
  Point(p1) = { x_c, y_c, z_c, clSphere };
  p2 = newp;
  Point(p2) = { x_c - rad, y_c, z_c, clSphere };
  p3 = newp;
  Point(p3) = { x_c, y_c - rad, z_c, clSphere };

  // определяем одну линию - полуокружность
  line = newl; // команда newl определяет новый уникальный номер линии
  Circle(line) = { p2, p1, p3 };

  // получаем поверхность - 1/8 часть всей поверхности шара
  surfSphere1[] = Extrude { { rad, 0, 0 }, { x_c, y_c, z_c }, Pi / 2 }
                          { Line{line}; };
  // получаем остальные части поверхности, копируя исходную поверхность
  // с зеркальным отображением по трем плоскостям
  surfSphere2[] = Symmetry { 1, 0, 0, -x_c }
                           { Duplicata { Surface{surfSphere1[1]}; } };
  surfSphere3[] = Symmetry { 0, 1, 0, -y_c }
                           { Duplicata { Surface{surfSphere1[1],
                                                 surfSphere2[0]}; } };
  surfSphere4[] = Symmetry { 0, 0, 1, -z_c }
                           { Duplicata { Surface{surfSphere1[1],
                                                 surfSphere2[0],
                                                 surfSphere3[0],
                                                 surfSphere3[1]}; } };

  // замкнутый контур поверхностей шара
  sl1 = newsl; // команда newsl определяет новый уникальный номер
               // контура поверхностей
  Surface Loop(sl1) = { surfSphere1[1], surfSphere2[0],
                        surfSphere3[0], surfSphere3[1],
                        surfSphere4[0], surfSphere4[1],
                        surfSphere4[2], surfSphere4[3] };
  surfaceLoops[number] = sl1; // массив номеров контуров поверхностей

  // задаем шар
  v1 = newv; // команда newv определяет новый уникальный номер объема
  Volume(v1) = { sl1 };
  sphereVolumes[number] = v1; // массив номеров объемов

Return

Несколько замечаний по поводу функций в Gmsh'е. Начинается функция со строки Function funcName Заканчивается словом Return (без ';' в конце). Никаких входных и выходных параметров нет. Просто-напросто те переменные, которые вы используете в функции, должны быть инициализированы к моменту вызова функции, который осуществляется командой Call funcName; (здесь уже ';' присутствует). В нашей функции определения шара это переменные x_c, y_c, z_c, rad, clSphere, number.

Приведем также список всех команд типа newp, newl, использующихся в нашем коде:

• newp определяет новый уникальный номер точки (point)
  
• newl определяет новый уникальный номер линии (line)
  
• newll определяет новый уникальный номер контура линий (line loop)
  
• news определяет новый уникальный номер поверхности (surface)
  
• newsl определяет новый уникальный номер контура поверхностей (surface loop)
  
• newv определяет новый уникальный номер объема (volume)
  
• newreg определяет новый уникальный номер региона (region) - может использоваться вместо любой из вышеперечисленных команд
  

Теперь несколько комментариев по поводу выбора такой реализации шара. Есть, конечно, и другие варианты. Например, можно определить шаблонный шар с центром в начале координат и единичным радиусом. Тогда наиболее простым подходом к созданию большого количества сферических включений в различных местах и желательно с различным радиусом было бы просто копировать исходный шар в заданные координаты (с помощью связки команд Translate и Duplicata). А поскольку команда Translate в качестве аргумента содержит вектор перемещения от одного центра к другому, очевидно, при таком расположении исходного шара вектор перемещения совпадал бы с координатами центра нового шара. А радиус равный 1 позволил бы в дальнейшем с легкостью создавать разномасштабные шары просто масштабируя на величину нового радиуса (командой Dilate). Этот вариант тоже имеет право на жизнь, однако сходу реализовать его у меня не получилось. А создавать что-то менее универсальное не хотелось. Поэтому в качестве шаблона определения шара выбрана функция.
Соберем теперь наши готовые геометрические элементы в одном файле, который и определит целиком геометрию исследуемого объекта.

model.geo

Include "brick.geo";
Include "sphere.geo";

// пусть центры включений расположены не случайно,
// а в узлах некоторой трехмерной решетки.
// подход, основанный на отдельном задании координат центров включений,
// позволит без большого изменения кода перейти
// к случаю, когда включения разбросаны по параллелепипеду случайным образом

// количество включений по каждому измерению
nInclusionsX = 2;
nInclusionsY = 2;
nInclusionsZ = 2;

// радиус сферического включения
radIncl = lenX / 10; // величина lenX взята из "brick.geo"

// общее количество включений
nIncl = nInclusionsX * nInclusionsY * nInclusionsZ;

// координаты центров включений
n = 0;
For i In { 1:nInclusionsZ }
  For j In { 1:nInclusionsY }
    For k In { 1:nInclusionsX }
      // x0, y0, z0, lenX, lenY, lenZ взяты из "brick.geo"
      x[n] = x0 + k * lenX / (nInclusionsX + 1);
      y[n] = y0 + j * lenY / (nInclusionsY + 1);
      z[n] = z0 + i * lenZ / (nInclusionsZ + 1);
      n = n + 1;
    EndFor
  EndFor
EndFor

// создание включений
For i In { 1:nIncl }
  x_c = x[i-1];
  y_c = y[i-1];
  z_c = z[i-1];
  rad = radIncl; // здесь можно варьировать радиус шаров,
                 // создавая разномасштабные включения
  clSphere = rad; // "мелкость" сетки в каждом шаре также можно варьировать
  number = i; // номер включения
  Call SphereTemplate; // вызов функции
EndFor

// окончательно оформляем массив номеров контуров поверхностей
surfaceLoops[0] = brickSurfaceLoop;

// создаем объем на основе массива контуров.
// этот объем - то, что не занято включениями
v = newreg;
Volume(v) = { surfaceLoops[] };

Physical Volume(1001) = { v }; // скелет (матрица)
Physical Volume(1002) = { sphereVolumes[] }; // поры (включения)

Для построения сетки остается только открыть в Gmsh'е файл "model.geo".
Отмечу, что код во многом основан на примере t5.geo из туториала Gmsh'а. Также можно заметить, что код получился вполне универсальным в том смысле, что добавив функции-шаблоны, задающие включения цилиндрического, эллипсоидного и других типов, можно создать вполне интересный инструмент определения геометрии объектов со сложной структурой микровключений. Единственное условие - включения не должны пересекаться. Приведенный код это не отслеживает.
Ну и на десерт - что получилось у меня с использованием данного кода (привожу только включения: слева 2х2х2, справа - 3х3х3)