3D mesh generation with predefined surface regions using NetGen

    Introduction


    In a previous post, we examined the features of using third-party open source libraries Freefem ++ and NetGen in a program for modeling aerodynamic processes. It was about the possibility of including these libraries in a commercial project from the standpoint of licensing, about the features of performing functions and inclusion in the software architecture. This article is a supplement to the previous one, in it we will examine the NetGen library in more detail. Of interest are the functions of generating a 3D mesh of finite elements with specified regions on the surface of the model.

    NetGen connection in MS Visual Studio


    We show how you can connect NetGen to a C ++ program. Download the archive from the official site of the NetGen project . In our case, NetGen version 4.9.13 was available. To connect the library, you need three folders from the archive: libsrc, nglib, windows. In Visual Studio, create a project and a solution for the demo program and attach the existing project file nglib.vcxproj from the windows folder. Since in our experiment, the nglib project was created in an earlier version of Visual Studio, we simultaneously convert the project to a new version. In the nglib project settings, add the libsrc \ include include folder, add the NO_PARALLEL_THREADS symbol, remove pthreadVC2.lib from the additional linker dependencies, delete the linker post-event, and configure the target folders so that the project builds without errors.

    The inclusion of the nglib.h header file should be written in the following not very familiar way:

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

    Test model


    The space that should be filled with finite elements can be specified for NetGen in two ways:
    1. In the form of constructive block geometry (CSG) operators.
    2. In the form of a description of the space boundary in the STL file.

    Here we consider only the second variant of the surface description — the description of the three-dimensional model in the STL text format. A graphic image of the test model is shown in Fig. 1. The model is a parallelepiped of size 10 X 10 X 5. On one of the faces is a rectangular region of size 5 X 2 with an offset (4; 2), the boundaries of which should remain unchanged
    image
    . 1. Surface triangulation to form an STL file. The region whose borders should remain unchanged is highlighted in red.

    The following is the contents of the STL file describing the test model:
    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  

    The information in the STL file has the following format. At the beginning and end of the file, the keywords solid and endsolid are indicated, between which triangles are listed that specify the required surface in the format:
    facet normal nx ny nz 
    outer loop 
    vertex v1x v1y v1z 
    vertex v2x v2y v2z 
    vertex v3x v3y v3z 
    endloop 
    endfacet
    where nx, ny, nz are the components of the normal vector along the axes X, Y, Z, directed outward from the simulated object; (v1x, v1y, v1z), (v2x, v2y, v2z), (v3x, v3y, v3z) - the coordinates of the three vertices of the triangle.

    Example of a finite element mesh generating program


    Below is the text of a program that performs all the steps necessary to obtain a finite element mesh. This program reads the OneLegionBar.stl STL file and writes a finite element grid to the OneRegionBar.vol file.

    #include "stdafx.h"
    #include 
    #include 
    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";
    	// Ribs on the surface
    	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 Edge
    	Edge edges [EDGES_NUM] = {Edge (0,1), Edge (1, 2), Edge (2,3), Edge (3, 0)};
    	// Parameters of mesh generation
    	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;
    	// The result of operations
    	Ng_Result ng_res;
    	//
    	// Start meshing
    	//
    	// Read the STL file, create the "STL geometry" object
    	Ng_STL_Geometry * stlGeometry = Ng_STL_LoadGeometry (stlFileName);
    	// Adding immutable edges to the "geometry"
    	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);
    	}
    	// Initialize an object with an STL description
    	ng_res = Ng_STL_InitSTLGeometry (stlGeometry);
    	if (ng_res! = NG_OK)
    	{
    		cout << "Error in Ng_STL_InitSTLGeometry" << endl;
    	}
    	// So far, the empty element mesh object
       	Ng_Mesh * mesh = Ng_NewMesh ();
    	// Form edges on the model surface
    	ng_res = Ng_STL_MakeEdges (stlGeometry, mesh, & ngMeshParameters);
    	if (ng_res! = NG_OK)
    	{
    		cout << "Error in Ng_STL_MakeEdges" << endl;
    	}
    	// Generate a space surface
    	ng_res = Ng_STL_GenerateSurfaceMesh (stlGeometry, mesh, & ngMeshParameters);
    	if (ng_res! = NG_OK)
    	{
    		cout << "Error in Ng_STL_GenerateSurfaceMesh" << endl;
    	}
    	// Generate a finite element network
    	ng_res = Ng_GenerateVolumeMesh (mesh, & ngMeshParameters);
    	if (ng_res! = NG_OK)
    	{
    		cout << "Error in Ng_GenerateVolumeMesh" << endl;
    	}
    	// Save the mesh of finite elements to a file
    	Ng_SaveMesh (mesh, volFileName);
    } 

    Next, we give detailed comments on the program text:
    1. The coordinates of the points that are the edges of the edges are in the points [POINTS_NUM] array. The edges [EDGES_NUM] array contains edges of the surface region in the form of corresponding pairs of point numbers. These edges will need to be specified by NetGen so that they are preserved on the surface of the model.
    2. The process of generating the grid is controlled by the parameters collected in the ngMeshParameters object of the Ng_Meshing_Parameters class. One of the parameters - maxh limits the maximum size of the grid elements. Using other parameters, you can fine-tune the generation process: limit the minimum size of elements, initiate grid optimization, etc.
    3. The STL file is read by the Ng_STL_LoadGeometry () function, and the STL geometry object is returned:
            Ng_STL_Geometry * stlGeometry = Ng_STL_LoadGeometry (stlFileName); 
      
      where Ng_STL_Geometry should be understood as a type of pointer to an object with STL geometry. In fact, Ng_STL_Geometry, like some other types, is synonymous with void *. NetGen authors in this way completely hid the implementation features of complex objects.
    4. Next, add edges to the STL geometry object. They are added in a loop using the Ng_STL_AddEdge () function:
            Ng_STL_AddEdge (stlGeometry, p1, p2);
      
      where stlGeometry is an object of STL geometry; p1 and p2 are the vertices of the edge. Each vertex is defined by an array of components along the X, Y, Z axes. These vertices are selected from the array of edges edges, which was formed earlier.
    5. Call the Ng_STL_InitSTLGeometry () function to initialize the STL geometry object:
            ng_res = Ng_STL_InitSTLGeometry (stlGeometry);
      

    6. Create an empty element mesh object that mesh points to. To do this, the Ng_NewMesh () function is executed:
             Ng_Mesh * mesh = Ng_NewMesh (); 
      

    7. Form edges on the model surface in accordance with the STL geometry and generation parameters:
             ng_res = Ng_STL_MakeEdges (stlGeometry, mesh, & ngMeshParameters); 
      

    8. Triangulate the model surface with the Ng_STL_GenerateSurfaceMesh () function. Here, the previously defined edges that remain undistorted are already taken into account, and the generation parameters will be taken into account:
      	
             ng_res = Ng_STL_GenerateSurfaceMesh (stlGeometry, mesh, & ngMeshParameters);
      

    9. And, finally, the element grid is directly generated by the Ng_GenerateVolumeMesh () function taking into account the generation parameters:
             ng_res = Ng_GenerateVolumeMesh (mesh, & ngMeshParameters);
      

    10. In this program, the grid generation result is output to the VOL-file by means of the Ng_SaveMesh () function, which transfers the pointer to the grid and the file name:
             Ng_SaveMesh (mesh, volFileName);
      


    Grid generation result


    The result of mesh generation is a program object by a pointer of type Ng_Mesh, in a demo program it is a mesh pointer. There is a set of functions that provide programmatic access to the properties of the grid of elements: get the number of points in the grid Ng_GetNP (), get the number of triangles on the surface Ng_GetNSE (), get the number of finite elements Ng_GetNE (), get the coordinates of the point Ng_GetPoint (), get the surface element Ng_GetSurfaceElement ( ), get the final element Ng_GetVolumeElement (), save the grid in the Ng_SaveMesh file.

    The VOL file saved in the program can be opened in the netgen.exe visualizer. The result of viewing the file is shown in Fig. 2. As can be seen in the figure, the boundaries of the surface region remained unchanged.


    Fig. 2. The surface of the model after generating a mesh of elements

    The VOL file generated by NetGen has the following format:
    1. At the beginning of the file, the format code (mesh3d), the dimension of the model (dimension 3), and the geometry code (geomtype 0) are indicated:
      mesh3d
      dimension
      3
      geomtype
      0 
      

    2. Next, after the keyword surfaceelements, the quantity is indicated and the surface elements are listed. For each element, the surface number, the number of “material” of the surface, the reserved integer field, the number of points that describe the surface element (for triangles - 3), and the numbers of the three points of the vertices of the surface element are recorded. Hereinafter, points are numbered from 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. The file lists the finite elements, for each of which the number of “material” is recorded, the number of vertices of one element (there are four for tetrahedrons) and the numbers of points that are vertices:
      # matnr np p1 p2 p3 p4
      volumeelements
      1009
             1 4 213 85 153 214
             1 4 298 307 301 305
      # ...
      

    4. The information about the edges on the surface is indicated:
      # 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. The points that are the nodes of the grid are listed:
      # XYZ
      points
      333
         10.0000000000000000 4.0000000000000000 4.0000000000000000
         10.0000000000000000 4.0000000000000000 2.0000000000000000
      # ...
      

    6. At the end of the file are written color components for surfaces identified by numbers. These colors are used by the visualizer to paint the surface. Here, for the surface with number 2, the color is set so as to highlight the rectangular region in Fig. 2.
      # Surfnr Red Green Blue
      face_colours
      7
             2 1.00000000 1.00000000 0.00000000
             3 0.00000000 1.00000000 0.00000000
      # ...
      


    2D meshing with NetGen


    NetGen has a feature for building a 2D mesh of elements. In the particular case, this function can be useful, for example, for generating an STL file for a 3D model with flat faces. The geometry of the model is described in the input file as the boundaries of the regions of space. The boundaries of the subdomains of the partition are specified either using line segments or using quadratic splines. When describing the border, indicate the number of the area to the left and to the right of the border. The area outside the grid is encoded by number 0.


    Fig. 3. Coding of points and areas for generating a 2D mesh of elements.

    We give an example of an input file corresponding to Fig. 3.
    splinecurves2dv2
    2
    points
    100
    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 
    

    The splinecurves2dv2 keyword is indicated at the beginning of this file. Following it is the grid rebuild factor (here - 2). After the keyword points, the model description points are listed: the point number and coordinates along the X, Y axes.

    Next, after the segments keyword, a list of boundary segments follows in the following format:
    1. Area number to the left of the border.
    2. Area number to the right of the border.
    3. The number of points to describe the boundary segment (in example 2).
    4. The number of the starting point of the border.
    5. The endpoint number of the boundary.
    6. Generation control flags. The bc flag (boundary condition number) must be indicated even when not in use.

    After the keyword materials, the areas (“materials”) are listed in the following format:
    1. Area number.
    2. The name of the material.
    3. Generation control flags. The maxh flag can be specified here, limiting the maximum size of the element grid.

    The following is a program that generates a 2D mesh:
    #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);
    }
    

    To generate a 2D mesh, an object for describing 2D geometry using the geom pointer is used. The result of mesh generation is in the object by the mesh pointer. In the example, the same function Ng_SaveMesh () is used to display the result of the generation, as in the example of generating the 3D mesh. Therefore, in this case, the output file has a structure for a 3D mesh, but the value of the Z coordinate in it is always zero. The program interface for working with a 2D grid contains the functions Ng_GetNP_2D () - get the number of grid nodes, Ng_GetNE_2D () - get the number of elements in the grid, Ng_GetPoint_2D () - get the grid node coordinates by the node number, Ng_GetElement_2D () - get the coordinates of the element vertices by item number. A visualization of the result of generating a 2D mesh is shown in Fig. 4. As can be seen from the figure, the borders of the region are not distorted.


    Fig. 4. The result of generating a 2D mesh

    Conclusion


    This article gives a brief overview of NetGen, from which it follows that this library can be conveniently used in applications that require the generation of a 3D mesh of elements. Library functions can be learned and used without writing code. NetGen comes with a netgen.exe application that implements a visual interface to all library functions. The drawings in this article, depicting a grid of elements, are obtained using this application.

    Also popular now: