[CFD Workshop] Model mesh (triangular mesh)

Preface
mesh geometry
grid number
Programming implementation
- Data reading
- Grid data construction

This series of blog posts is about my process of learning the theory of two-dimensional shallow water equations and compiling a practical two-dimensional shallow water flow model. (Review of the previous blog)

Foreword

This 2D shallow water model will use an unstructured grid. To keep things simple, I’ll build the model and write the solver starting from a triangular mesh.
Many code libraries and software can generate triangle meshes. For example, DHI MIKE software, Triangle library, etc. In this regard, I personally recommend SMS software; SMS software can easily generate triangular meshes by importing boundary grid points, and can check the quality of the generated meshes.
Of course, we need to export these grid data, edit them into the data we need, and finally input them into this shallow water model. This blog post will clarify what data we export.

Grid Geometry

First, I defined the parts of the triangle mesh as follows:

Grid node (point) refers to the three vertices that form a triangle, as shown in the hollow circle in the figure below;
The mesh surface (cell) refers to the entire triangle;
Cell center refers to the geometric center of the triangle, as shown in the red circle in the figure below;
Grid edges (edges), an edge is determined by two grid nodes (points).

Generally speaking, the two necessary elements that constitute a mesh are mesh nodes (points) and mesh edges (edges); when the vertices of the triangle are drawn and the mesh edges are connected, a set of triangle meshes can be Sure.
But in general, our commonly used input data are mesh nodes and the vertex numbers that make up each triangle. For grid nodes, the input data is the serial number, x and y coordinate values of each vertex; after that, the serial number sequence that constitutes each triangle is given. Taking the above figure as an example, the point data is as follows:

Node number	x coordinate	y coordinate
1	x₁	y₁
2	x₂	y₂
3	x₃	y₃
4	x₄	y₄
5	x ₅	y₅

The data that constitutes each cell is as follows (this type of data is referred to as “cell” below, that is, the default input elements are point and cell data):

Grid number	The three grid points (grid point number) that constitute the grid
1	1, 2, 5
2	2, 3, 5
3	3, 4, 5

Grid number

Based on the above example, we know that point and cell data do not share the same numbering system. For example, there is no necessary relationship between grid No. 1 and node No. 1, but the relationship between grid triangles and nodes is described by cell data.
Here, we emphasize that the points, edges and cells in the triangular grid each use a numbering system; only the cell center and cell data can share the same numbering system. system. Therefore, in this model, the number of the grid triangle cell is equal to the number of its center point cell center.

Programming implementation

Data reading

First, we read the grid coordinates from the txt file through functions such as fopen and fgets, recorded as xp and yp. Both xp and yp are one-dimensional values of size Np (xp[0 ~ (Np-1)]), and Np represents the total number of grids. It should be noted that in C language, the starting point of array index is 0.
After that, read in the node number cell that constitutes each triangle and record it as a node array. The node data is a two-dimensional array with a size of node[0 ~ 2][0 ~ (Nc-1)]; where Nc represents the total number of triangles.

Here we give a sample code:

/*
* Function: ReadGrid
* Usage: Read the grid data from the input files
* Called by: main
*/
void ReadGrid() {
int num=0, i;

// 1. to read grid points
FILE* fp1 = fopen("points.txt", "r");
if (fp1 == NULL) {
fprintf(stderr, "Fail to read grid files (points.txt)! \\
");
exit(EXIT_FAILURE);
}
char* row1 = (char*)malloc(sizeof(char) * N_str);
fgets(row1, N_str, fp1);
sscanf(row1, "%d", & amp;num);

xp = (double *)malloc(sizeof(double) * num);
yp = (double *)malloc(sizeof(double) * num);
zbp = (double *)malloc(sizeof(double) * num);
for (i = 0; i < num; i + + ) {
fgets(row1, N_str, fp1);
sscanf(row1, "%lf %lf \\
", & amp;xp[i], & amp;yp[i]);
}
Np = num;
free(row1);
fclose(fp1);

// 2. to read geometric data of all cells
FILE* fp2 = fopen("cells.txt", "r");
if (fp2 == NULL) {
fprintf(stderr, "Fail to read grid files (elements.txt)! \\
");
exit(EXIT_FAILURE);
}
char* row2 = (char*)malloc(sizeof(char) * N_str);
fgets(row2, N_str, fp2);
sscanf(row2, "%d", & amp;num);

node = (int **)malloc(sizeof(int) * 3);
for (int i = 0; i < 3; i + + )
{
node[i] = (int *)malloc(sizeof(int) * num);
}
for (i = 0; i < num; i + + ) {
fgets(row2, N_str, fp1);
sscanf(row2, "%d %d %d \\
", & amp;node[0][i], & amp;node[1][i], & amp;node[2][i] );
\t\t
// redefine the index from zero instead of one
node[0][i] = node[0][i] - 1;
node[1][i] = node[1][i] - 1;
node[2][i] = node[2][i] - 1;

}
Nc = num;
free(row2);
fclose(fp2);

fprintf(stderr, "The grid data has been read... \\
");
}

An example of the code required to read into the file points.txt is as follows (the first line represents Np):

5
0.0 0.0
1.0 1.0
4.0 1.0
5.0 0.0
3.0 0.0

An example of the code required to read into the file cells.txt is as follows (the first line represents Nc):

Grid data construction

In the model, we also need to know the coordinates of the center point, the grid area, the length of each side, the adjacent relationship of the grid, etc. (this will be introduced later). Here, we first take these quantities as examples to introduce the processing of grid data.

/*
* Function: ConstructGrid
* Usage: Construct the grids for reading other data
* Called by: main
*/
void ConstructGrid() {
int i, j, k, p1, p2, p3, count, count2;
double x1, x2, x3, y1, y2, y3;

xc = (double *)malloc(sizeof(double) * Nc);
yc = (double *)malloc(sizeof(double) * Nc);
area = (double *)malloc(sizeof(double) * Nc);
neigh = (int **)malloc(sizeof(int) * 3);
for (int i = 0; i < 3; i + + )
{
neigh[i] = (int *)malloc(sizeof(int) * Nc);
}
edgeL = (double **)malloc(sizeof(int) * 3);
for (int i = 0; i < 3; i + + )
{
edgeL[i] = (double *)malloc(sizeof(int) * Nc);
}

for (i = 0; i < Nc; i + + ) {
p1 = node[0][i];
p2 = node[1][i];
p3 = node[2][i];

x1 = xp[p1];
x2 = xp[p2];
x3 = xp[p3];
y1 = yp[p1];
y2 = yp[p2];
y3 = yp[p3];

// 1. to define the cell centers and calculate the area
xc[i] = (x1 + x2 + x3) / 3.0;
yc[i] = (y1 + y2 + y3) / 3.0;
area[i] = fabs(x1*y2 - x2*y1 + x2*y3 - x3*y2 + x3*y1 - x1*y3);

// 2. to find the neiborn cells
count = 0;
for (j = 0; j < Nc; j + + ) {
if (count >= 3) break; // "count=3" means that all neighboring cells have been found
if (j == i) continue; // skip the itself
count2 = 0;
for (k = 0; k < 3; k + + ) {
if (node[k][j] == p1 || node[k][j] == p2 || node[k][j] == p3) count2 = count2 + 1;
}
if (count2 == 2) {
neigh[count][i] = j;
count = count + 1;
}
}
\t\t
if (count == 1) {
neigh[1][i] = -1;
neigh[2][i] = -1;
}
else if(count == 2){
neigh[2][i] = -1;
}

// 3. to calculate the lengths of all edges
edgeL[0][i] = sqrt((x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1));
edgeL[1][i] = sqrt((x3 - x2)*(x3 - x2) + (y3 - y2)*(y3 - y2));
edgeL[2][i] = sqrt((x1 - x3)*(x1 - x3) + (y1 - y3)*(y1 - y3));

}
fprintf(stderr, "The computational grid has been constructed... \\
");
}

First of all, in the first part, we defined the arrays xc and yc to record the coordinates of the grid center point. These two variables can open up storage space through the malloc function (lines 10 to 11 of the above code). The same applies to other variables, and the malloc function can be used to define the space.
After that, we calculate xc and yc at the center of the grid (lines 37 ~ 38 of the above code); then calculate the area area of each grid (line 39 of the above code). The area formula is as follows:

∣

Area = |x_1 y_2 – x_2 y_1 + x_2 y_3 – x_3 y_2 + x_3 y_1 – x_1 y_3|

Area=∣x1?y2x2?y1? + x2?y3x3?y2? + x3?y1x1?y3?∣
After that, we search the adjacent grids of each grid, which is very helpful for us to solve the flux between grids. Because within each time step, only adjacent grids have mass and momentum flux transfer. We store the adjacency of the grid through the array neigh[0 ~ 2][1 ~ (Nc-1)]. In lines 42 ~ 62 of the above code, we search adjacent grids. The idea is as follows:

For grid i and grid j, determine whether they have two identical vertices; this can be done through their node[ ][i] and node[ ][j]
to judge. When it is found that node[ ][i] and node[ ][j] have two identical elements (lines 46 ~ 53 of the above code), record i and j at this time (lines 50 ~ 53 of the above code).
A grid may have 1, 2 or 3 adjacent grids. For the case where the number of adjacent grids is 3, the neighbor array does not require additional processing; when the number of adjacent grids is 1 or 2, the excess elements of the neighbor array will be filled with -1

Taking the grid in the above figure as an example, after the search is completed, the neigh array results obtained are as follows (note the number here, the starting point of the array in C language is 0):

1 -1 -1
0 2 -1
1 -1 -1

In addition, please note that the processing method of the second point above is temporary. The edges of grids with adjacent grid numbers of 1 or 2 will correspond to the boundaries of the calculation area. Therefore, the neigh array needs to be processed later.
Finally, we calculate the length of each edge and store it in the edgeL array (lines 65 to 67 of the above code).