Генерация 3D-сетки с предопределенными регионами поверхности при помощи NetGen

Введение

В предыдущем посте мы рассмотрели особенности использования сторонних библиотек с открытым кодом Freefem++ и NetGen в программе моделирования аэродинамических процессов. Речь шла о возможности включения этих библиотек в коммерческий проект с позиции лицензирования, об особенностях выполнения функций и включения в программную архитектуру. Данная статья является дополнением к предыдущей, в ней более детально рассмотрим библиотеку NetGen. Интерес представляют функции генерации 3D-сетки конечных элементов с заданными регионами на поверхности модели.

Подключение NetGen в MS Visual Studio

Покажем, как можно подключить NetGen к программе на C++. С официального сайта проекта NetGen скачиваем архив. В нашем случае был доступен NetGen версии 4.9.13. Для подключения библиотеки нужны три папки из архива: libsrc, nglib, windows. В Visual Studio создаем проект и решение для демонстрационной программы и подключаем существующий файл проекта nglib.vcxproj из папки windows. Так как в нашем эксперименте проект nglib был создан в более ранней версии Visual Studio, то попутно выполняем конвертацию проекта под новую версию. В настройках проекта nglib добавляем include-папку libsrc\include, добавляем символ NO_PARALLEL_THREADS, удаляем pthreadVC2.lib из дополнительных зависимостей компоновщика, удаляем пост-событие компоновщика и настраиваем целевые папки так, чтобы проект собирался без ошибок.

Включение заголовочного файла nglib.h надо записывать следующим не очень привычным способом:

namespace nglib
{
#include "nglib.h"
}
using namespace nglib;

Тестовая модель

Пространство, которое должно быть заполнено конечными элементами, может задаваться для NetGen двумя способами:

1. В виде операторов конструктивной блочной геометрии (CSG).
2. В виде описания границы пространства в файле формата STL.

Здесь рассмотрим только второй вариант описания поверхности – описание трехмерной модели в текстовом формате STL. Графическое изображение тестовой модели показано на рис. 1. Модель представляет собой параллелепипед размером 10 Х 10 Х 5. На одной из граней располагается прямоугольный регион размером 5 Х 2 со смещением (4;2), границы которого должны оставаться неизменными

Рис. 1. Триангуляция поверхности для формирования STL файла. Красным цветом выделен регион, границы которого должны оставаться без изменений

Ниже приводится содержимое STL-файла, описывающего тестовую модель:

solid
facet normal 1 0 0 outer loop vertex 10 4 2 vertex 10 6.5 3 vertex 10 4 4 endloop endfacet
facet normal 1 0 0 outer loop vertex 10 6.5 3 vertex 10 9 4 vertex 10 4 4 endloop endfacet
facet normal 1 0 0 outer loop vertex 10 9 2 vertex 10 6.5 3 vertex 10 4 2 endloop endfacet
facet normal 1 0 0 outer loop vertex 10 9 2 vertex 10 9 4 vertex 10 6.5 3 endloop endfacet
facet normal 1 0 0 outer loop vertex 10 0 5 vertex 10 0 0 vertex 10 4 2 endloop endfacet
facet normal 1 0 0 outer loop vertex 10 4 2 vertex 10 4 4 vertex 10 0 5 endloop endfacet
facet normal 1 0 0 outer loop vertex 10 10 5 vertex 10 0 5 vertex 10 4 4 endloop endfacet
facet normal 1 0 0 outer loop vertex 10 9 4 vertex 10 10 5 vertex 10 4 4 endloop endfacet
facet normal 1 0 0 outer loop vertex 10 9 2 vertex 10 10 5 vertex 10 9 4 endloop endfacet
facet normal 1 0 0 outer loop vertex 10 10 0 vertex 10 10 5 vertex 10 9 2 endloop endfacet
facet normal 1 0 0 outer loop vertex 10 10 0 vertex 10 9 2 vertex 10 4 2 endloop endfacet
facet normal 1 0 0 outer loop vertex 10 10 0 vertex 10 4 2 vertex 10 0 0 endloop endfacet
facet normal -1 0 0 outer loop vertex 0 10 5 vertex 0 10 0 vertex 0 0 5 endloop endfacet
facet normal -1 0 0 outer loop vertex 0 0 0 vertex 0 0 5 vertex 0 10 0 endloop endfacet
facet normal 0 -1 0 outer loop vertex 0 0 5 vertex 0 0 0 vertex 10 0 5 endloop endfacet
facet normal 0 -1 0 outer loop vertex 10 0 0 vertex 10 0 5 vertex 0 0 0 endloop endfacet
facet normal 0 1 0 outer loop vertex 10 10 5 vertex 10 10 0 vertex 0 10 5 endloop endfacet
facet normal 0 1 0 outer loop vertex 0 10 0 vertex 0 10 5 vertex 10 10 0 endloop endfacet
facet normal 0 0 1 outer loop vertex 0 0 5 vertex 10 0 5 vertex 0 10 5 endloop endfacet
facet normal 0 0 1 outer loop vertex 10 10 5 vertex 0 10 5 vertex 10 0 5 endloop endfacet
facet normal 0 0 -1 outer loop vertex 10 0 0 vertex 0 0 0 vertex 10 10 0 endloop endfacet
facet normal 0 0 -1 outer loop vertex 0 10 0 vertex 10 10 0 vertex 0 0 0 endloop endfacet
endsolid  

Информация в STL-файле имеет следующий формат. В начале и конце файла указываются ключевые слова solid и endsolid, между которыми перечисляются треугольники, задающие требуемую поверхность в формате:

facet normal nx ny nz 
outer loop 
vertex v1x v1y v1z 
vertex v2x v2y v2z 
vertex v3x v3y v3z 
endloop 
endfacet

где nx, ny, nz – составляющие вектора нормали по осям X, Y, Z, направленного наружу из моделируемого объекта; (v1x, v1y, v1z), (v2x, v2y, v2z), (v3x, v3y, v3z) – координаты трех вершин треугольника.

Пример программы, генерирующей сетку конечных элементов

Ниже приведен текст программы, выполняющей все действия, необходимые для получения сетки конечных элементов. Данная программа читает STL-файл OneRegionBar.stl и записывает сетку конечных элементов в файл OneRegionBar.vol.

#include "stdafx.h"
#include <iostream>
#include <fstream>

namespace nglib
{
#include "..\NetGen\nglib\nglib.h"
}

using namespace std;
using namespace nglib;

void main()
{
	const char *stlFileName = "..\\Data\\OneRegionBar.stl";
	const char *volFileName = "..\\Data\\OneRegionBar.vol";

	// Ребра на поверхности
	const int EDGES_NUM = 4;
	const int POINTS_NUM = 4;
	typedef double Point[3];
	Point points[POINTS_NUM] = {{10, 4, 2}, {10, 9, 2}, {10, 9, 4}, {10, 4, 4}};
	typedef pair<int, int> Edge;
	Edge edges[EDGES_NUM] = {Edge(0,1), Edge(1, 2), Edge(2,3), Edge(3, 0)};

	// Параметры генерации сетки
	Ng_Meshing_Parameters ngMeshParameters;
	
	ngMeshParameters.maxh = 10; 
	ngMeshParameters.fineness = 0.4;
	ngMeshParameters.secondorder = 0;
	ngMeshParameters.quad_dominated = 0;

	ngMeshParameters.grading = 0.0;
	ngMeshParameters.optsurfmeshenable = false;
	ngMeshParameters.optsteps_2d = 0;
	ngMeshParameters.closeedgeenable = false;

	// Результат выполнения операций
	Ng_Result ng_res;

	//
	// Начало построения сетки
	//

	// Чтение STL-файла, создание объекта "STL-геометрии"
	Ng_STL_Geometry *stlGeometry = Ng_STL_LoadGeometry(stlFileName);

	// Добавление в "геометрию" неизменяемых ребер
	for(int i = 0; i < EDGES_NUM; i++)
	{
		double* p1 = points[edges[i].first];
		double* p2  = points[edges[i].second];
		Ng_STL_AddEdge(stlGeometry, p1, p2);
	}

	// Инициализация объекта с описанием STL
	ng_res = Ng_STL_InitSTLGeometry(stlGeometry);
	if(ng_res != NG_OK)
	{
		cout << "Error in Ng_STL_InitSTLGeometry" << endl;
	}

	// Пока что пустой объект сетки элементов
   	Ng_Mesh *mesh = Ng_NewMesh();

	// Формировать ребра на поверхности модели
	ng_res = Ng_STL_MakeEdges(stlGeometry, mesh, &ngMeshParameters);
	if(ng_res != NG_OK)
	{
		cout << "Error in Ng_STL_MakeEdges" << endl;
	}

	// Генерировать поверхность пространства
	ng_res = Ng_STL_GenerateSurfaceMesh(stlGeometry, mesh, &ngMeshParameters);
	if(ng_res != NG_OK)
	{
		cout << "Error in Ng_STL_GenerateSurfaceMesh" << endl;
	}

	// Генерировать сеть конечных элементов
	ng_res = Ng_GenerateVolumeMesh(mesh, &ngMeshParameters);
	if(ng_res != NG_OK)
	{
		cout << "Error in Ng_GenerateVolumeMesh" << endl;
	}

	// Сохранить сетку конечных элементов в файл
	Ng_SaveMesh(mesh,volFileName);
} 

Далее дадим подробные комментарии к тексту программы:

1. Координаты точек, являющихся границами ребер, находятся в массиве points[POINTS_NUM]. В массиве edges[EDGES_NUM] заданы ребра региона поверхности в виде соответствующих пар номеров точек. Эти ребра надо будет указать NetGen, чтобы они сохранились на поверхности модели.
   
2. Управление процессом генерации сетки осуществляется при помощи параметров, собранных в объекте ngMeshParameters класса Ng_Meshing_Parameters. Один из параметров – maxh ограничивает максимальный размер элементов сетки. С помощью других параметров можно делать тонкую настройку процесса генерации: ограничить минимальный размер элементов, инициировать оптимизацию сетки и т.д.
   
3. STL-файл считывается функцией Ng_STL_LoadGeometry(), а возвращается объект STL-геометрии:
   

   
         Ng_STL_Geometry *stlGeometry = Ng_STL_LoadGeometry(stlFileName); 
   
   

   где Ng_STL_Geometry следует понимать как тип указателя на объект с STL-геометрией. На самом деле Ng_STL_Geometry, как и некоторые другие типы, является синонимом void*. Авторы NetGen таким способом полностью скрыли особенности реализации сложных объектов.
   
4. Далее в объект STL-геометрии надо добавить ребра. Они добавляются в цикле при помощи функции Ng_STL_AddEdge():
   

   
         Ng_STL_AddEdge(stlGeometry, p1, p2);
   
   

   где stlGeometry – объект STL-геометрии; p1 и p2 – вершины ребра. Каждая вершина задается массивом из компонент по осям X, Y, Z. Эти вершины выбираются из массива ребер edges, который был сформирован ранее.
   
5. Выполнить вызов функции Ng_STL_InitSTLGeometry() для инициализации объекта STL-геометрии:
   

   
         ng_res = Ng_STL_InitSTLGeometry(stlGeometry);
   

   
   
6. Создать пустой объект сетки элементов, на который указывает mesh. Для этого выполняется функция Ng_NewMesh():
   

   
          Ng_Mesh *mesh = Ng_NewMesh(); 
   

   
   
7. Формировать ребра на поверхности модели в соответствии с STL-геометрией и параметрами генерации:
   

   
          ng_res = Ng_STL_MakeEdges(stlGeometry, mesh, &ngMeshParameters); 
   

   
   
8. Выполнить триангуляцию поверхности модели функцией Ng_STL_GenerateSurfaceMesh(). Здесь уже учитываются ранее заданные ребра, которые останутся неискаженными, и будут учтены параметры генерации:
   

   	
          ng_res = Ng_STL_GenerateSurfaceMesh(stlGeometry, mesh, &ngMeshParameters);
   

   
   
9. И, наконец, выполняется непосредственно генерация сетки элементов функцией Ng_GenerateVolumeMesh() с учетом параметров генерации:
   

   
          ng_res = Ng_GenerateVolumeMesh(mesh, &ngMeshParameters);
   

   
   
10. В данной программе результат генерации сетки выводится в VOL-файл посредством функции Ng_SaveMesh(), которой передается указатель на сетку и имя файла:
    

    
           Ng_SaveMesh(mesh,volFileName);
    

    
    

Результат генерации сетки

Результатом генерации сетки является программный объект по указателю типа Ng_Mesh, в демонстрационной программе это указатель mesh. Имеется набор функций, обеспечивающих программный доступ к свойствам сетки элементов: получить количество точек в сетке Ng_GetNP(), получить количество треугольников на поверхности Ng_GetNSE(), получить количество конечных элементов Ng_GetNE(), получить координаты точки Ng_GetPoint(), получить элемент поверхности Ng_GetSurfaceElement(), получить конечный элемент Ng_GetVolumeElement(), сохранить сетку в файл Ng_SaveMesh.

VOL-файл, сохраненный в программе, может быть открыт в визуализаторе netgen.exe. Результат просмотра файла показан на рис. 2. Как видно на рисунке, границы региона поверхности остались неизменными.


Рис. 2. Поверхность модели после генерации сетки элементов

VOL-файл, формируемый NetGen, имеет следующий формат:

1. В начале файла указывается код формата (mesh3d), размерность модели (dimension 3), код геометрии (geomtype 0):
   

   mesh3d
   dimension
   3
   geomtype
   0 
   

   
   
2. Далее после ключевого слова surfaceelements указывается количество и перечисляются элементы поверхности. Для каждого элемента записаны номер поверхности, номер «материала» поверхности, зарезервированное целочисленное поле, количество точек, которыми описывается элемент поверхности (для треугольников — 3), номера трех точек вершин элемента поверхности. Здесь и далее точки нумеруются от 1.
   

   # surfnr    bcnr   domin  domout      np      p1      p2      p3
   surfaceelements
   546
          2       1       1       0       3       3       4      13
          2       1       1       0       3       4       5     107 
   #...
   

   
   
3. В файле перечисляются конечные элементы, для каждого из которых записаны номер «материала», количество вершин одного элемента (для тетраэдров их четыре) и номера точек, являющихся вершинами:
   

   #  matnr      np      p1      p2      p3      p4
   volumeelements
   1009
          1       4     213      85     153     214
          1       4     298     307     301     305
   #...
   

   
   
4. Указывается информация о ребрах на поверхности:
   

   # surfid  0   p1   p2   trignum1    trignum2   domin/surfnr1    domout/surfnr2   ednr1   dist1   ednr2   dist2 
   edgesegmentsgi2
   220
          1       0       1       2        1        1        0        0        1            0        1            2
          2       0       2       1        6        6        0        0        1            2        1            0
   #...
   

   
   
5. Перечисляются точки, являющиеся узлами сетки:
   

   #          X             Y             Z
   points
   333
      10.0000000000000000      4.0000000000000000      4.0000000000000000
      10.0000000000000000      4.0000000000000000      2.0000000000000000
   #...
   

   
   
6. В конце файла записаны компоненты цвета для поверхностей, идентифицируемых номерами. Эти цвета используются визуализатором для раскрашивания поверхности. Здесь для поверхности с номером 2 задан цвет так, чтобы выделить прямоугольный регион на рис. 2.
   

   #   Surfnr     Red     Green     Blue
   face_colours
   7
          2   1.00000000   1.00000000   0.00000000
          3   0.00000000   1.00000000   0.00000000
   #...
   

   
   

Формирование 2D-сетки при помощи NetGen

В NetGen имеется функция построения 2D-сетки элементов. В частном случае эта функция может быть полезна, например, для формирования STL файла для 3D-модели, имеющей плоские грани. Геометрия модели описывается во входном файле в виде границ областей пространства. Границы подобластей разбиения задаются либо с помощью отрезков прямых, либо с помощью квадратических сплайнов. При описании границы надо указывать номер области слева и справа от границы. Область вне сетки кодируется номером 0.


Рис. 3. Кодирование точек и областей для генерации 2D-сетки элементов

Приведем пример входного файла, соответствующего рис. 3.

splinecurves2dv2
2
points
1 0 0
2 3 0
3 3 2
4 0 2
5 1 0
6 2 0
7 2 1
8 1 1
segments
1 0 2 1 5 -bc=1
2 0 2 5 6 -bc=1
1 0 2 6 2 -bc=1
1 0 2 2 3 -bc=1
1 0 2 3 4 -bc=1
1 0 2 4 1 -bc=1
2 1 2 6 7 -bc=1
2 1 2 7 8 -bc=1
2 1 2 8 5 -bc=1
 
materials
1 domain1 -maxh=1
2 domain2 -maxh=1 

В начале этого файла указывается ключевое слово splinecurves2dv2. Следом за ним коэффициент перестроения сетки (здесь – 2). После ключевого слова points перечисляются точки описания модели: номер точки и координаты по осям X, Y.

Далее после ключевого слова segments следует перечень сегментов границ в следующем формате:

1. Номер области слева от границы.
2. Номер области справа от границы.
3. Количество точек для описания сегмента границы (в примере 2).
4. Номер начальной точки границы.
5. Номер конечной точки границы.
6. Флаги управления генерацией. Флаг bc (номер граничного условия) приходится указывать даже тогда, когда он не используется.

После ключевого слова materials перечисляются области («материалы») в следующем формате:

1. Номер области.
2. Название материала.
3. Флаги управления генерацией. Здесь может быть указан флаг maxh, ограничивающий максимальный размер сетки элементов.

Ниже приведена программа, генерирующая 2D-сетку:

#include "stdafx.h"

namespace nglib
{
#include "..\NetGen\nglib\nglib.h"
}

using namespace nglib;
void main(){
	const char *in2DFileName = "..\\Data\\triangulation.in2d";
	const char *volFileName = "..\\Data\\triangulation.vol";

	Ng_Geometry_2D * geom;

	Ng_Init();

	geom = Ng_LoadGeometry_2D(in2DFileName);

	Ng_Meshing_Parameters mp;
	mp.maxh = 10000;
	mp.fineness = 1;
	mp.secondorder = 0;

	Ng_Mesh * mesh;
	Ng_GenerateMesh_2D (geom, &mesh, &mp);

	Ng_SaveMesh(mesh, volFileName);
}

Для генерации 2D-сетки используется объект описания 2D-геометрии по указателю geom. Результат генерации сетки находится в объекте по указателю mesh. В примере для вывода результата генерации используется та же самая функция Ng_SaveMesh(), что и в примере генерации 3D-сетки. Поэтому в данном случае выходной файл имеет структуру для 3D-сетки, но значение координаты Z в нем везде равно нулю. Программный интерфейс для работы с 2D-сеткой содержит функции Ng_GetNP_2D() – получить количество узлов сетки, Ng_GetNE_2D() – получить количество элементов в сетке, Ng_GetPoint_2D() – получить координаты узла сетки по номеру узла, Ng_GetElement_2D() – получить координаты вершин элемента по номеру элемента. Визуализация результата генерации 2D-сетки показана на рис. 4. Как видно из рисунка, границы региона не искажены.


Рис. 4. Результат генерации 2D-сетки

Заключение

В данной статье дано краткое представление NetGen, из которого следует, что эта библиотека удобным образом может быть использована в приложениях, требующих генерацию 3D-сетки элементов. Функции библиотеки могут быть изучены и использованы и без написания кода. В поставку NetGen входит приложение netgen.exe, которое реализует визуальный интерфейс ко всем функциям библиотеки. Рисунки в данной статье, изображающие сетки элементов, получены с помощью этого приложения.