Practice using Freefem ++
In a previously published post, we talked about using the open source libraries Freefem ++ and NetGen in a program for modeling aerodynamic processes. This article will examine in more detail the basic features of Freefem ++ as a small introduction to its input language. This will provide initial information that developers often need when choosing third-party components for inclusion in the designed application.
Freefem ++ is a program designed to solve mathematical problems based on the finite element method. Freefem ++ can be used in both Windows and Linux environments. The official website is www.freefem.org/ff++ . In our case, Freefem ++ version 3.18.-1 was used. The archive bundle, downloaded from the official site, has three boot modules, the purpose of which is as follows:
To run FreeFem ++ in the background as part of another program from the archive, in addition to the boot module, the following files are required: libff.dll, libgcc_s_dw2-1.dll, libgfortran-3.dll, libgoto2.dll, libstdc ++ - 6.dll.
When the Freefem ++ boot module is launched, the command line will be given the name of the script file written in a special C-like programming language. Using the output operators or special save functions in the script, the results can be written to a file for further use.
Earlier, we showed that Freefem ++ can be used to form an STL description of a three-dimensional model consisting of elements with flat faces. To solve this problem, it is necessary to perform triangulation of faces. Let's see how this is done using Freefem ++. Let it be required to perform triangulation of the surface defined by the boundaries (see Fig. 1).
Fig. 1. Surface Regions
The boundaries of the areas for which triangulation is necessary are specified in Freefem ++ in a parametric form. Consider the following code:
Freefem ++ defines scalar types of variables: int, real, string, bool. The variable declaration is performed as in C ++, and the variable can be initialized, as is done for the variable n.
Borders are defined by fragments using parameterized lines. The border type is for this purpose intended. Border fragments can only intersect at the ends. The declaration of a boundary fragment includes the identifier of the boundary, the boundary of the parameter change (in this case, parameter t) and expressions for changing the x and y coordinates depending on the parameter. In the example, all fragments are segments.
The last operator in the considered code is plot, it is intended for graphical display of data of certain types. In this case, it displays a border consisting of fragments of borders. Border fragments are combined into a common border using the + operator. For each fragment, a parameter indicating the number of parts into which the fragment should be divided is indicated. This parameter will be considered again in another context below. The result of the specified plot operator is shown in Fig. 2.
Fig. 2. The boundaries of the regions of the surface displayed by the plot operator.
Now we add the code that performs the triangulation. In fig. 3 shows the result of its implementation.
Fig. 3. The result of triangulation
In this code, a variable Th of type mesh is introduced. The mesh type is a software model of a triangulation mesh. To obtain an object of this type, in this case, the buildmesh function is used, to which the boundary of the regions to be divided is transmitted to the input. The buildmesh function can be passed an optional fixeborder parameter, the use of which ensures that points on the line will not be changed during the triangulation process. The second operator - plot - in this case, takes an object of type mesh at the input and shows the grid of elements in the graphics window (Fig. 3).
As already shown, a border is represented as the sum of fragments of borders that are of type border. The parameter specified for each fragment has a dual purpose here. Firstly, it indicates how many parts the fragment of the border is divided into to build the grid. Each part of the border fragment will become a side of the triangle of the formed mesh. Secondly, the sign of the parameter indicates the area that should undergo triangulation. If the parameter is positive, then triangulation is performed to the left of the border fragment. If it is negative, then triangulation is performed to the right of the border. Let us illustrate this with an example. If the boundary were given in the above code as follows: b0 (n) + b1 (n) + b2 (n) + b3 (n) + b4 (n) + b5 (-2) + b6 (n) + b7 (n) + b8 (n) + b9 (-2) + b10 (n) + b11 (n) + b12 (n) + b13 (-2) + b14 (-2) + b15 (n), where for b5, b9, b13 and b14 are divided into two segments with a negative value, then the triangulation would be performed as shown in Fig. 4.
Fig. 4. Triangulation with an empty subdomain.
Add code that demonstrates operations with a variable of type mesh, a loop operator, and output to a file.
To output data to a file, the ofstream type is provided, an analog of the file output stream from the fstream library for C ++. When initializing a file variable, you must specify the file name. If you only need to add data to the file, you need to specify the second append initialization parameter. You can perform a number of operations with a file variable. In the example, the precision function sets the number of displayed digits after the decimal separator. There are several more functions that control the output format, similar to the manipulators from the iostream library. The << operator allows you to output values of variables of scalar types to a file. In the character string, you can use the end of line character '\ n', or for this purpose output the global variable endl.
The above code has a for loop statement that completely repeats the logic of the for statement in C ++. The increment and decrement operators also repeat the similar C ++ operators.
Freefem ++ provides access to various properties of the mesh of elements represented by the mesh type. Here we note that Freefem ++ provides the same interfaces for working with two-dimensional mesh type mesh, and with three-dimensional mesh type mesh3. Therefore, here we consider this interface for the two types indicated, so as not to repeat its descriptions for mesh3 later. As you can see from the example, the nt property returns the number of grid elements (triangles for mesh or tetrahedrons for mesh3). You can also get the number of grid points - the nv property, the number of boundary elements - nbe.
Two types of indexing can be applied to an object of type mesh or mesh3: indexing with square brackets "[]" and indexing with parentheses "()". The first of them gives a grid element (a triangle for mesh or a tetrahedron for mesh3) with a given index. Thus, Th [i] is the grid element (triangle in this example), to which, in turn, certain functions provided for the element are applied. "()" - indexing provides a vertex with a given index. Operations to obtain its coordinates can be applied to a vertex. The following functions are applied to the grid element: indexing to get the vertices of the element, getting a label.
Freefem ++ defines the type of mesh3 for a 3D mesh. Such a grid can be obtained using the generation functions. For example, the function buildlayers () allows you to create a 3D mesh by pulling 2D triangulation along the Z axis. However, Freefem ++ does not have its own generator. For these purposes, a third-party TetGen generator is used (the official website wias-berlin.de/software/tetgen ), and the corresponding Freefem ++ functions use the program interface to this generator. TetGen may not be used in commercial projects without the permission of the author. But if such a license creates a problem, then you can use another grid generator. Freefem ++ has the ability to import a 3D mesh from a file that conforms to the mesh format. The mesh file has the following format:
The mesh file is imported using the readmesh3 function, which specifies the file name:
For importing a mesh from a mesh file, the sequence order of the vertex numbers of the tetrahedrons is critical. In the Freefem ++ source code, in the Mesh3dn.hpp file, you can find the DataTet structure, which has a method for calculating the volume of the mesure () tetrahedron. If the vertices of the tetrahedron in the mesh file follow in such a way that the calculated volume is negative, then Freefem ++ fails. Therefore, the vertices in the mesh file must be listed so that this situation does not occur.
If necessary, the grid can be displayed in the graphic window with the plot function (see Fig. 5).
Fig. 5. 3D mesh display with plot function
The main purpose of Freefem ++ is to solve systems of partial differential equations using the finite element method. For this, the functions involved in the system of equations are determined on a finite element grid. Freefem ++ allows you to determine the type of FE space (finite element space). In accordance with the type of space, objects of functions of a certain space are introduced. For instance,
where Vh is the type of space on the grid Th, f is a function on this space.
When declaring the type of a finite element space, a basis for the elements is indicated. It can take the values P03d - piecewise constant, P13d - piecewise linear, P23d - piecewise quadratic, etc. In fig. Figure 6 shows the function from the example with a quadratic basis, and Fig. 7 - the same function with a linear basis and with noticeable distortion.
Fig. 6. Finite-element function with a quadratic basis
. 7. Finite element function with linear basis
To solve the system of partial differential equations in Freefem ++, a variational statement of the problem is specified, and when the problem is solved, the desired functions are formed. Based on the following script, we will consider how the variational formulation and other programming elements are written.
The following constructions are presented in this example. A finite element space can be specified for a vector field. In the example, the space VVh is specified for a vector field with four components, and Vh for a scalar field.
Freefem ++ provides the ability to declare macros that begin with the macro keyword. During macro substitution, the macro name will be replaced with the macro body with the replacement of formal parameters with actual ones. In the example, there are two macros: Grad (u) and div (u1, u2, u3). These macros use the special notation dx, du, dz, which denote differentials by their respective coordinates.
A function is declared using the func keyword. Moreover, if the arguments of the function are denoted by x, y, z, which are reserved words, then you do not need to specify formal parameters.
The solve construction is used to record the variational formulation and solution of the problem. After the task name, unknown and test functions are listed. In this case, the Prob problem has two vector functions defined on the space VVh: [u1, u2, u3, p] - the desired one and [v1, v2, v3, q] - the trial one. After the equal sign, the expression for the variational formulation is written. In the example at this point, there are two specific Freefem ++ constructs. To record a three-dimensional integral, a special int3d construction is used in which the integration space is specified, in this case, the grid Th. In addition, the matrix transposition operator, denoted by the apostrophe, is used here. At the end of the variational formulation, boundary conditions are written in the solve construction. The boundary condition is indicated by the special word on, after which the label corresponding to the area from the object of class mesh3 is indicated in brackets. Writing these labels to the mesh file was discussed above. For the region, the components of the desired field are listed and a constraint is written for each of them. A constraint may be an arbitrary expression. In the example, numbers and the fup function are specified as a constraint.
When the field component is output using the plot function, the finished element space with isosurfaces is shown in the graphic window, as shown in Fig. 8. The number of isosurfaces can be set with the nbiso parameter in the script, and then manually changed in the graphics window.
Fig. 8. Result of the decision
Freefem ++ is a powerful tool for solving mathematical problems, which can be used in the development of applications, including commercial ones. In the latter case, a useful feature of Freefem ++ is the ability to import a 3D mesh of finite elements, which allows the use of third-party mesh generators.
Getting and executing Freefem ++
Freefem ++ is a program designed to solve mathematical problems based on the finite element method. Freefem ++ can be used in both Windows and Linux environments. The official website is www.freefem.org/ff++ . In our case, Freefem ++ version 3.18.-1 was used. The archive bundle, downloaded from the official site, has three boot modules, the purpose of which is as follows:
- FreeFem ++. Exe is a full-fledged version of the program, which, if necessary, displays a graphical window with the results of the solution. If necessary, when calling, you can specify two keys: -nowait so that Freefem ++ exits without confirming completion by the user, and -nw - to refuse to display a graphical window with the solution results.
- FreeFem ++ - mpi.exe - for parallelizing MPI-based computing.
- FreeFem ++ - nw.exe - for calling from applications, ignores the graphical output of the results. In order for him not to wait for confirmation of completion of work with the program, you need to specify the -nowait switch on the command line.
To run FreeFem ++ in the background as part of another program from the archive, in addition to the boot module, the following files are required: libff.dll, libgcc_s_dw2-1.dll, libgfortran-3.dll, libgoto2.dll, libstdc ++ - 6.dll.
When the Freefem ++ boot module is launched, the command line will be given the name of the script file written in a special C-like programming language. Using the output operators or special save functions in the script, the results can be written to a file for further use.
2D triangulation
Earlier, we showed that Freefem ++ can be used to form an STL description of a three-dimensional model consisting of elements with flat faces. To solve this problem, it is necessary to perform triangulation of faces. Let's see how this is done using Freefem ++. Let it be required to perform triangulation of the surface defined by the boundaries (see Fig. 1).
Fig. 1. Surface Regions
The boundaries of the areas for which triangulation is necessary are specified in Freefem ++ in a parametric form. Consider the following code:
int n = 5;
border b0 (t = 0,1) {x = 7 * t; y = 0; }
border b1 (t = 0,1) {x = 7; y = 0 + 6 * t; }
border b2 (t = 0,1) {x = 7 - 7 * t; y = 6; }
border b3 (t = 0,1) {x = 0; y = 6 - 6 * t; }
border b4 (t = 0,1) {x = 1 + 2 * t; y = 1.5; }
border b5 (t = 0,1) {x = 5; y = 1.5 + 2 * t; }
border b6 (t = 0,1) {x = 5 - 4 * t; y = 5.5; }
border b7 (t = 0,1) {x = 1; y = 5.5 - 4 * t; }
border b8 (t = 0,1) {x = 5; y = 3.5 + 2 * t; }
border b9 (t = 0,1) {x = 3 + 2 * t; y = 1.5; }
border b10 (t = 0,1) {x = 3 + 3 * t; y = 0.5; }
border b11 (t = 0,1) {x = 6; y = 0.5 + 3 * t; }
border b12 (t = 0,1) {x = 6 - 1 * t; y = 3.5; }
border b13 (t = 0,1) {x = 3; y = 3.5 - 2 * t; }
border b14 (t = 0,1) {x = 5 - 2 * t; y = 3.5; }
border b15 (t = 0,1) {x = 3; y = 1.5 - 1 * t; }
plot (b0 (n) + b1 (n) + b2 (n) + b3 (n) + b4 (n) + b5 (n) + b6 (n) + b7 (n) + b8 (n) + b9 (n ) + b10 (n)
+ b11 (n) + b12 (n) + b13 (n) + b14 (n) + b15 (n));Freefem ++ defines scalar types of variables: int, real, string, bool. The variable declaration is performed as in C ++, and the variable can be initialized, as is done for the variable n.
Borders are defined by fragments using parameterized lines. The border type is for this purpose intended. Border fragments can only intersect at the ends. The declaration of a boundary fragment includes the identifier of the boundary, the boundary of the parameter change (in this case, parameter t) and expressions for changing the x and y coordinates depending on the parameter. In the example, all fragments are segments.
The last operator in the considered code is plot, it is intended for graphical display of data of certain types. In this case, it displays a border consisting of fragments of borders. Border fragments are combined into a common border using the + operator. For each fragment, a parameter indicating the number of parts into which the fragment should be divided is indicated. This parameter will be considered again in another context below. The result of the specified plot operator is shown in Fig. 2.
Fig. 2. The boundaries of the regions of the surface displayed by the plot operator.
Now we add the code that performs the triangulation. In fig. 3 shows the result of its implementation.
mesh Th = buildmesh (b0 (n) + b1 (n) + b2 (n) + b3 (n) + b4 (n) + b5 (n) + b6 (n) + b7 (n) + b8 (n) + b9 (n) + b10 (n) + b11 (n) + b12 (n) + b13 (n) + b14 (n) + b15 (n), fixeborder = 1); plot (Th);
Fig. 3. The result of triangulation
In this code, a variable Th of type mesh is introduced. The mesh type is a software model of a triangulation mesh. To obtain an object of this type, in this case, the buildmesh function is used, to which the boundary of the regions to be divided is transmitted to the input. The buildmesh function can be passed an optional fixeborder parameter, the use of which ensures that points on the line will not be changed during the triangulation process. The second operator - plot - in this case, takes an object of type mesh at the input and shows the grid of elements in the graphics window (Fig. 3).
As already shown, a border is represented as the sum of fragments of borders that are of type border. The parameter specified for each fragment has a dual purpose here. Firstly, it indicates how many parts the fragment of the border is divided into to build the grid. Each part of the border fragment will become a side of the triangle of the formed mesh. Secondly, the sign of the parameter indicates the area that should undergo triangulation. If the parameter is positive, then triangulation is performed to the left of the border fragment. If it is negative, then triangulation is performed to the right of the border. Let us illustrate this with an example. If the boundary were given in the above code as follows: b0 (n) + b1 (n) + b2 (n) + b3 (n) + b4 (n) + b5 (-2) + b6 (n) + b7 (n) + b8 (n) + b9 (-2) + b10 (n) + b11 (n) + b12 (n) + b13 (-2) + b14 (-2) + b15 (n), where for b5, b9, b13 and b14 are divided into two segments with a negative value, then the triangulation would be performed as shown in Fig. 4.
Fig. 4. Triangulation with an empty subdomain.
Add code that demonstrates operations with a variable of type mesh, a loop operator, and output to a file.
ofstream fout ("triangulation.out");
fout.precision (14);
fout << Th.nt << endl;
for (int i = 0; i <Th.nt; i ++)
{
fout
<< Th [i] [0] .x << ""
<< Th [i] [0] .y << ""
<< Th [i] [1] .x << ""
<< Th [i] [1] .y << ""
<< Th [i] [2] .x << ""
<< Th [i] [2] .y << ""
<< endl;
}
To output data to a file, the ofstream type is provided, an analog of the file output stream from the fstream library for C ++. When initializing a file variable, you must specify the file name. If you only need to add data to the file, you need to specify the second append initialization parameter. You can perform a number of operations with a file variable. In the example, the precision function sets the number of displayed digits after the decimal separator. There are several more functions that control the output format, similar to the manipulators from the iostream library. The << operator allows you to output values of variables of scalar types to a file. In the character string, you can use the end of line character '\ n', or for this purpose output the global variable endl.
The above code has a for loop statement that completely repeats the logic of the for statement in C ++. The increment and decrement operators also repeat the similar C ++ operators.
Freefem ++ provides access to various properties of the mesh of elements represented by the mesh type. Here we note that Freefem ++ provides the same interfaces for working with two-dimensional mesh type mesh, and with three-dimensional mesh type mesh3. Therefore, here we consider this interface for the two types indicated, so as not to repeat its descriptions for mesh3 later. As you can see from the example, the nt property returns the number of grid elements (triangles for mesh or tetrahedrons for mesh3). You can also get the number of grid points - the nv property, the number of boundary elements - nbe.
Two types of indexing can be applied to an object of type mesh or mesh3: indexing with square brackets "[]" and indexing with parentheses "()". The first of them gives a grid element (a triangle for mesh or a tetrahedron for mesh3) with a given index. Thus, Th [i] is the grid element (triangle in this example), to which, in turn, certain functions provided for the element are applied. "()" - indexing provides a vertex with a given index. Operations to obtain its coordinates can be applied to a vertex. The following functions are applied to the grid element: indexing to get the vertices of the element, getting a label.
3D mesh
Freefem ++ defines the type of mesh3 for a 3D mesh. Such a grid can be obtained using the generation functions. For example, the function buildlayers () allows you to create a 3D mesh by pulling 2D triangulation along the Z axis. However, Freefem ++ does not have its own generator. For these purposes, a third-party TetGen generator is used (the official website wias-berlin.de/software/tetgen ), and the corresponding Freefem ++ functions use the program interface to this generator. TetGen may not be used in commercial projects without the permission of the author. But if such a license creates a problem, then you can use another grid generator. Freefem ++ has the ability to import a 3D mesh from a file that conforms to the mesh format. The mesh file has the following format:
- Information about the format version (MeshVersionFormatted 1), space dimension (Dimension 3) is recorded at the beginning of the file.
MeshVersionFormatted 1 Dimension 3
- After the Vertices keyword, the number of grid vertices is indicated and the vertices are listed for which X, Y, Z coordinates and the vertex label are indicated.
Vertices 65 7 5 1.5 0 7 3 1.5 0 7 3 0.5 0 7 5 3.5 0 ...
- After the Tetrahedra keyword, the number of grid elements is indicated and the elements for which four vertex numbers and a label are written are listed.
Tetrahedra 185 52 56 49 50 0 3 18 14 61 0 47 64 57 54 0 ...
- After the Triangles keyword, the number of boundary triangles is recorded. For each triangle, the numbers of the vertices of the grid nodes that are the vertices of these triangles are listed, and the border label is indicated.
Triangles 96 2 3 8 2 2 8 1 2 4 1 7 2 ...
- The file ends with the End keyword.
The mesh file is imported using the readmesh3 function, which specifies the file name:
mesh3 Th;
Th = readmesh3 ("example3d.mesh");
plot (Th);
For importing a mesh from a mesh file, the sequence order of the vertex numbers of the tetrahedrons is critical. In the Freefem ++ source code, in the Mesh3dn.hpp file, you can find the DataTet structure, which has a method for calculating the volume of the mesure () tetrahedron. If the vertices of the tetrahedron in the mesh file follow in such a way that the calculated volume is negative, then Freefem ++ fails. Therefore, the vertices in the mesh file must be listed so that this situation does not occur.
If necessary, the grid can be displayed in the graphic window with the plot function (see Fig. 5).
Fig. 5. 3D mesh display with plot function
Solution of systems of differential equations
The main purpose of Freefem ++ is to solve systems of partial differential equations using the finite element method. For this, the functions involved in the system of equations are determined on a finite element grid. Freefem ++ allows you to determine the type of FE space (finite element space). In accordance with the type of space, objects of functions of a certain space are introduced. For instance,
fespace Vh (Th, P23d);
Vh f = x ^ 2 + y ^ 2 + z ^ 2; where Vh is the type of space on the grid Th, f is a function on this space.
When declaring the type of a finite element space, a basis for the elements is indicated. It can take the values P03d - piecewise constant, P13d - piecewise linear, P23d - piecewise quadratic, etc. In fig. Figure 6 shows the function from the example with a quadratic basis, and Fig. 7 - the same function with a linear basis and with noticeable distortion.
Fig. 6. Finite-element function with a quadratic basis
. 7. Finite element function with linear basis
To solve the system of partial differential equations in Freefem ++, a variational statement of the problem is specified, and when the problem is solved, the desired functions are formed. Based on the following script, we will consider how the variational formulation and other programming elements are written.
mesh3 Th = readmesh3 ("example3d.mesh");
fespace VVh (Th, [P2, P2, P2, P1]);
fespace Vh (Th, P23d);
macro Grad (u) [dx (u), dy (u), dz (u)] //
macro div (u1, u2, u3) (dx (u1) + dy (u2) + dz (u3)) //
VVh [u1, u2, u3, p];
VVh [v1, v2, v3, q];
func fup = (1-x) * (x) * y * (1-y) * 16;
solve Prob ([u1, u2, u3, p], [v1, v2, v3, q]) =
int3d (Th, qforder = 3) (Grad (u1) '* Grad (v1) + Grad (u2)' * Grad (v2) + Grad (u3) '* Grad (v3)
- div (u1, u2, u3) * q - div (v1, v2, v3) * p + 1e-10 * q * p)
+ on (5, u1 = 0, u2 = 0, u3 = 0)
+ on (1, u1 = fup, u2 = 0, u3 = 0);
plot (p, nbiso = 10);
The following constructions are presented in this example. A finite element space can be specified for a vector field. In the example, the space VVh is specified for a vector field with four components, and Vh for a scalar field.
Freefem ++ provides the ability to declare macros that begin with the macro keyword. During macro substitution, the macro name will be replaced with the macro body with the replacement of formal parameters with actual ones. In the example, there are two macros: Grad (u) and div (u1, u2, u3). These macros use the special notation dx, du, dz, which denote differentials by their respective coordinates.
A function is declared using the func keyword. Moreover, if the arguments of the function are denoted by x, y, z, which are reserved words, then you do not need to specify formal parameters.
The solve construction is used to record the variational formulation and solution of the problem. After the task name, unknown and test functions are listed. In this case, the Prob problem has two vector functions defined on the space VVh: [u1, u2, u3, p] - the desired one and [v1, v2, v3, q] - the trial one. After the equal sign, the expression for the variational formulation is written. In the example at this point, there are two specific Freefem ++ constructs. To record a three-dimensional integral, a special int3d construction is used in which the integration space is specified, in this case, the grid Th. In addition, the matrix transposition operator, denoted by the apostrophe, is used here. At the end of the variational formulation, boundary conditions are written in the solve construction. The boundary condition is indicated by the special word on, after which the label corresponding to the area from the object of class mesh3 is indicated in brackets. Writing these labels to the mesh file was discussed above. For the region, the components of the desired field are listed and a constraint is written for each of them. A constraint may be an arbitrary expression. In the example, numbers and the fup function are specified as a constraint.
When the field component is output using the plot function, the finished element space with isosurfaces is shown in the graphic window, as shown in Fig. 8. The number of isosurfaces can be set with the nbiso parameter in the script, and then manually changed in the graphics window.
Fig. 8. Result of the decision
conclusions
Freefem ++ is a powerful tool for solving mathematical problems, which can be used in the development of applications, including commercial ones. In the latter case, a useful feature of Freefem ++ is the ability to import a 3D mesh of finite elements, which allows the use of third-party mesh generators.