SFEMaNS  version 5.3
Reference documentation for SFEMaNS
Mesh Generator

Download

First download the mesh generator of SFEMaNS by typing the command:

wget http://math.tamu.edu/~guermond/DOWNLOADS/MESH_GENERATOR_SFEMaNS.tar.gz

Extract this file in the directory where you want to install the mesh generator:

tar -xvf MESH_GENERATOR_SFEMaNS.tar.gz

It creates a directory called MESH_GENERATOR_SFEMaNS. To compile the executables that allow to generate P1-P2 finite element meshes, type the following commands:

cd MESH_GENERATOR_SFEMaNS
make clean
make all

It generates three executables called maill.exe, colle_p1p1.exe and symmetrize.exe. Note that you may have to modify the file make.inc if you don't have access to mpi.f90. For pratical use, we advise the user to create aliases by typing:

alias SFEMANS_MESH_GEN_DIR="path to directory of the mesh generator"
alias maill.exe = "($SFEMANS_MESH_GEN_DIR)/Code/maill.exe"
alias colle_p1p2.exe = "($SFEMANS_MESH_GEN_DIR)/Collage/colle_p1p2.exe"
alias symmetrize.exe = "($SFEMANS_MESH_GEN_DIR)/Collage/symmetrize.exe"

Such aliases should be written in your .bashrc (or equivalent) to avoid redefining them everytime you log in.

Visualization

When generating a mesh, as detailled below, many files are created so one can check if the mesh is generated correctly. These files, that have the extension ".plt", can be visualized with the executable plotmtv with the command:

plotmtv file_to_visualize

This executable and its documentation are available in ($SFEMANS_MESH_GEN_DIR)/PLOTMTV. For practical use, you can create an alias by typing:

alias plotmtv="($SFEMANS_MESH_GEN_DIR)/PLOTMTV/plotmtv"

Plotmtv offers various mode of visualization (zooming, rotation, 2D and 3D rendering, etc.). You can also create postscript file by clicking on print to file. Such a postscript file is always named dataplot.ps so don't forget to modify its name after creation. Note that you can also edit files with the extension ".plt".

Generate a mesh

We now describe the procedure to be followed to create mesh of a two-dimensional domain of your choice. To do that you need three files named:

my_project.in to set the boundary of the mesh.
topology.name_of_mesh to set the topology of the mesh at the boundary.
backgrid.name_of_mesh to set specific topology at points inside the domain.

These files allow to create a grid with P1 elements with the following command:

maill.exe < my_project.in

In order to create a grid with P1-P2 elements, you need a file named data_colle and type the following command:

colle_p1p2.exe

Eventually the P1-P2 mesh can be symmetrized with respect to a given line. It is done by editing the file data_symmetrize and typing the command:

symmetrize.exe

The following sections describe the above procedure on a template example. Before going on, you need to go to your work directory and type the following command:

cp ($SFEMANS_MESH_GEN_DIR)/EXAMPLES/TEMPLATE/* .

It copies the following files my_example.in, topology.example, backgrid.example, data_colle and data_symmetrize in your directory. The following explains how to edit these files and how to use the above executables correctly.

Generation of a mesh P1

We describe how to modify the files my_example.in, topology.example and backgrid.example that allows to generate a mesh P1. We note that such a mesh can't be used by the code SFEMaNS which requires a P1-P2 mesh. Such a P1-P2 mesh is generated later using this P1 mesh.

The my_project.in file

The first step consists of defining the boundary of your domain. Draw your domain on a piece of paper. Divide the boundary into elementary paths. Starting from \(1\), number the paths composing the outer boundary going clockwise. Do the same for the path extremities. The path extremities are called vertices. If any is present, continue to number the inner boundaries and corresponding inner vertices by going anti-clockwise. The rule is that the domain must always be on your right-hand side as you go along boundaries.

The information on paths and vertices is stored in my_project.in. For instance copy my_example.in to my_project.in as follows

cp my_example.in my_project.in
  1. The first line of my_project.in is composed of the name of the project and an integer called integer_domain_index that, for the time being, you set to \(1\). It will be used later when assigning domain where each equations need to be solved in the finite element code.

    name_of_mesh integer_domain_index

  2. The second line contains an integer that is equal to the number of paths that composes the boundary of the domain, say \(n_p\).

  3. The following \(2 n_p\) lines describe the paths. Four types of paths are possible:

    1. Line segments: line.
    2. Circle pieces: circle.
    3. Ellipse pieces: ellipse.
    4. Paths tabulated by the user in a data file.

  4. If you want to specify a line segment, proceed as follows by typing:

    line integer_bdy_index

    \(x_1\) \(y_1\) \(x_2\) \(y_2\)

    Here "integer_bdy_index" is an integer that you assign to the path. It will be used latter when assigning boundary conditions in the finite element code. For the time being you can always set this number to \(1\) or to whatever value you want.

    \(x_1\) is the \(x\)-coordinate of the first vertex of the path. \(y_1\) is the \(y\)-coordinate of the first vertex of the path. \(x_2\) and \(y_2\) are the \(x\)- and \(y\)-coordinates of the second vertex of the path.

  5. If you want to specify a portion of circle, proceed as follows by typing:

    circle integer_bdy_index

    \(x_c\) \(y_c\) \(r\) \(\theta_1\) \(\theta_2\)

    where \(x_c\) and \(y_c\) are the \(x\)- and \(y\)-coordinates of the center of the circle; \(r\) is the radius of the circle; \(\theta_1\) is the angle, measured is degrees, of the first vertex and \(\theta_2\) is the angle of the second vertex. The origin of the angles is the horizontal half line starting at the circle center and pointing to the right.

  6. If you want to specify a piece of ellipse, proceed as follows by typing:

    ellipse integer_bdy_index

    \(x_c\) \(y_c\) \(a_x\) \(a_y\) \(\theta_1\) \(\theta_2\)

    where \(x_c\) and \(y_c\) are the \(x\)- and \(y\)-coordinates of the center of the ellipse; \(a_x\) is the first semi-axis and \(a_y\) is the second semi-axis; \(\theta_1\) is the angle, measured is degrees, of the first vertex and \(\theta_2\) is the angle of the second vertex. The origin of the angles is the horizontal half line starting at the center of the ellipse and pointing to the right.

  7. If you want to specify a path tabulated in a file, proceed as follows by typing:

    data integer_bdy_index

    file_containing_my_data

    where file_containing_my_data is a file written with the following format

    \[ \begin{matrix} \text{number_of_points} & &\\ x_1 & & y_1 \\ x_2 & & y_2 \\ \vdots & & \vdots \\ x_{\text{number_of_points}} & & y_{\text{number_of_points}} \end{matrix} \]

    The first line contains an integer that is equal to the number of points that are tabulated. The other lines contain the coordinates of the points composing the path. The first point is the first vertex of the path and the last point is the second vertex. The points must be ordered so that as one progresses along the path the domain is on the right-hand side.

The topology.name_of_mesh file

The second step of the work consists of defining the mesh size at the boundary of the domain. This is done by specifying the corresponding data in the file topology.name_of_mesh. The first part of the name, i.e. topology is fixed, whereas the second part, i.e. name_of_mesh in the present case, is specified by the user in the first line in the my_project.in file.

To prepare your own file copy topology.example in topology.name_of_mesh as follows

cp topology.example topology.name_of_mesh

and edit topology.name_of_mesh as follows.

  1. The first line of the file is a comment that explains the meaning of the two integers that are on the second line.

  2. The second line contains two integers. The first one is the number of paths composing the boundary and the second one is the number of vertices. These two numbers must be equal to \(n_p\).
  3. The third line is a blank line.
  4. Then there are packets of lines limited by the words BEGIN and END and these packets are separated by one blank line. There are as many packets as paths, i.e. \(n_p\) in all. The BEGIN and END lines are comments.
    1. The first line after BEGIN is a comment recalling the meaning of the integers in the second line.
    2. The second line contains four integers. The first integer is the number of the paths. This number must be identical to the order of appearance of the path in my_project.in. The second integer says at how many curvilinear abscissas one will give the value of the mesh size; let \(n_h\) be this integer. The third integer is the number of the first vertex of the path and the fourth integer is the number of the second vertex of the path.
    3. The third line is a comment to remind the meaning of the following lines.
    4. Then there are \(n_h\) lines before reaching the END line. Each line is composed of one integer and two real numbers. Always set the first integer to \(1\). The second number is the curvilinear abscissa where one wants to enforce the value of the mesh size. The curvilinear abscissa \(s\) is comprised between \(0\) and \(1\). One must always start with \(s=0\) and finish with \(s=1\). The third number is the value of the mesh size that one wants to set.

The backgrid.name_of_mesh file

It happens that one needs to to refine the grid inside the domain to increase the accuracy of the computation. This can be done by using the file backgrid.name_of_mesh. The first part of the name, i.e. backgrid is fixed, whereas the second part, i.e. name_of_mesh in the present case, is specified by the user in the first line of the my_project.in file. This files must exists even if the user does not want to refine the grid inside the domain.

Copy the file backgrid.example in backgrid.name_of_mesh as follows:

cp backgrid.example backgrid.name_of_mesh

Then edit backgrid.name_of_mesh for your own purposes.

  1. The first line of the file is a comment recalling the meaning of the integer in second line.
  2. The second line contains an integer that is equal to the number of points inside the domain where the user wants to enforce the mesh size, say \(n_r\). The number \(n_r\) must be set to zero if the user does not want to refine the mesh.
  3. The third line has to be the following:

    X Y H

    It is a comment recalling the structure of the following lines.

  4. The \(n_r\) lines thereafter contain three real numbers each, say

    \[ x \qquad y \qquad h \]

    \(x\) and \(y\) are the coordinate of the point where one wants to refine the mesh. \(h\) is the mesh size that one wants to impose.

We note that if \(n_r=0\), the mesh generator only reads the first three lines of the file backgrid.name_of_mesh.

The executable maill.exe

Once you have modified the three files described above, you can type the following command:

maill.exe < my_project.in

It will create a mesh with P1 finite element called FEM.name_of_mesh. Moreover a file called gridplot.name_of_mesh is also created. It can be visualized with plotmtv to check the correct construction of the mesh.

Generation of a P1-P2 mesh

To create a P1-P2 finite element mesh that the code SFEMaNS can used, you need to use the file named data_colle and the executable colle_p1p1.exe. After this step, you can symmetrize your P1-P2 mesh with respect to a given line with a file called data_symmetrize and the executable symmetrize.exe. We draw the attention of the reader that the symmetrization step is optional.

The data_colle file

This file allows to define the number of subdomains, meaning P1 meshes, that you are you going to glue together so you can generate a global P1-P2 mesh that will be used by the code SFEMaNS.

  1. The first line is a logical (.t. to generate formatted mesh and .f. to generate unformatted mesh)
  2. The second line is the number of P1-mesh to glue.
  3. The third line contains the path and the name of the first P1 mesh.
  4. If you want to create P1-P2 mesh that is the union of two or more P1 mesh, you need to add the three following lines for each P1 mesh to add:
    1. The first line is an integer \(\text{nb}_\text{inter}\) that represents the number of interfaces to keep between the two meshes.
    2. The second line contains \(\text{nb}_\text{inter}\) integer that are the index of the interfaces to keep (defined in the files *.in).
    3. The third line contains the path and the name of the P1 mesh to add.

Here is an example of a file data_colle that generates a formatted P1-P2 mesh by gluing two P1-meshes. We assume that theses meshes are in your working directory and are called FEM.name_of_mesh1 and FEM.name_of_mesh2. Moreover we decide to keep one interface of index 3.

.t. ! .t.= formatted .f.=unformatted
2 ! nb of subdomain
. FEM.name_of_mesh1
1 ! nb of interfaces to keep
3 ! index of interfaces
. FEM.name_of_mesh2

The executable colle_p1p2.exe

Once you modified your data_colle, you can type the following command:

colle_p1p2.exe

To the question " Attention to cavities: number of cavities = ", type the number of cavities present in your mesh. A cavity is a hole in the domain; it is not part of your mesh.

The executable generates a P1-P2 mesh called COLLE.FEM that can be used by the code SFEMaNS. Moreover it also generates files that can be visualized with plotmtv such as colle.plt (final P1 mesh), colle_dom1.plt (final P2 mesh). We note the file sides.plt, that represent the final P1 mesh, also displays the boundary and the interface you conserved. When visualizing sides.plt with plotmtv, you can click on "3D plot" to see the index of all of the boundaries and interfaces (previously defined in the file my_project.in).

Optional symmetrization step

The last feature of SFEMaNS's mesh generator is to symmetrize a P1-P2 mesh with respect to a given line. First you need to modify the file data_symmetrize that you previously copied in your work directory.

  1. The first line contains a logical (.t. if the unsymmetrized mesh is formatted and .f. if it is unformatted).
  2. The second line contains the directory of the mesh to symmetrize and the name of the mesh to symmetrize.
  3. The third line contains the three coefficient a, b and c where ax+by+c=0 is the equation of the symmetry line.
  4. The fourth (and last) line contains a logical (.t. to keep the interface between sub-meshes, .f. else).

Once you modified the file data_symmetrize, you can type the following command:

symmetrize.exe

It will generate a symmetrized P1-P2 mesh called symm_COLLE.FEM and "*.plt" files (such as sides_p2.plt, i_d_p2.plt, etc.) that can be visualized with plotmtv.

Example

The goal of this example is to use all the executables described above to generate a P1-P2 symmetrized mesh. We plan to construct a mesh that correspond to a cylinder with radius R=1 and height H=2 that is included in a sphere of center (0,0) and of radius 11. The symmetrization will be done with respect to the line y=0. The coordinates (x,y) correspond to the cylindrical coordinates (r,z) of the cylinder, meaning r for the radius and the z for the height in the direction of the revolution axis of the cylinder. All the files required to do such a mesh can be found in the directory:

($SFEMaNS_MESH_GEN_DIR)/EXAMPLES/EXAMPLE_DOC .

Generation of the P1 meshes

First we set the boundaries of the cylinder in the file rect_in.in as follows:

rect_in 1
4
line 1
0. 0. 0. 1.
line 2
0. 1. 1. 1.
line 3
1. 1. 1. 0.
line 4
1. 0. 0. 0.

Since we plan to symmetrize the mesh with respect to the line \(y=0\), this files represent a meridian plane of the cylinder's upper half.

Then we choose a mesh size of \(0.05\) in P1 by creating the file topology.rect_in as follows:

NV NE
4 4
BEGIN
E N BV EV
1 2 1 2
IDATA SDATA HDATA
1 0.0 0.05
1 1.0 0.05
END
BEGIN
E N BV EV
2 2 2 3
IDATA SDATA HDATA
1 0.0 0.05
1 1.0 0.05
END
BEGIN
E N BV EV
3 2 3 4
IDATA SDATA HDATA
1 0.0 0.05
1 1.0 0.05
END
BEGIN
E N BV EV
4 2 4 1
IDATA SDATA HDATA
1 0.0 0.05
1 1.0 0.05
END

Moreover we impose a mesh size of \( 0.025\) inside the domain at the point \((x,y)=(0.5,0.5)\) by setting the backgrid.rect_in as follows:

N
1
X Y H
0.5 0.5 0.025

We can now generate the P1 mesh with the command:

maill.exe < rect_in.in

It creates a file called FEM.rect_in that contains the information of the upper cylinder P1 mesh. In addition, a file called gridplot.rect_in is generated and can be visualized with plotmtv.

We now repeat the process for the spherical domain. While doing so, we need to ensure that the mesh size (topology) of the sphere and the cylinder coincide on their interfaces and also that these interfaces have the same index. So we create a the file circle_ext.in as follows:

circle_ext 2
5
circle 5
0 0 11 90.d0 0.d0
line 4
11 0. 1. 0
line 3
1.0 0. 1.0 1.0
line 2
1.0 1.0 0.0 1.0
line 1
0.0 1.0 0 11.0

where we keep the index 2 and 3 for the interfaces between the cylinder and the sphere. Note that the lines 1 and 4 have the same index that in rect.in but they are different boundaries.

Regarding the topology of the sphere, we use a mesh size of \(1\) on the outter boundary, meaning on \( x^2+y^2=11^2\). Moreover, we enforce a mesh size of \(0.05\) for \( x\leq 2 \) and \( y \leq 2\), which correspond to \( 10 \%\) of the arc lenght. The resulting file topology.circle_ext is written as follows:

NV NE
5 5
BEGIN
E N BV EV
1 2 1 2
IDATA SDATA HDATA
1 0.0 1.0
1 1.0 1.0
END
BEGIN
E N BV EV
2 3 2 3
IDATA SDATA HDATA
1 0.00 1.0
1 0.90 0.05
1 1.00 0.05
END
BEGIN
E N BV EV
3 2 3 4
IDATA SDATA HDATA
1 0.0 0.05
1 1.0 0.05
END
BEGIN
E N BV EV
4 2 4 5
IDATA SDATA HDATA
1 0.00 0.05
1 1.00 0.05
END
BEGIN
E N BV EV
5 3 5 1
IDATA SDATA HDATA
1 0.00 0.05
1 0.10 0.05
1 1.00 1.00
END

We do not impose specific refinement inside the domain so we define the file backgrid.circle_ext as follows:

N
0
X Y H

We create the P1 mesh by executing the following commands:

maill.exe < circle_ext.in

It generates a P1 mesh called FEM.circle_ext and a plotmtv file called backgrid.circle_ext.

Generation of the P1-P2 global mesh

To create a global P1-P2 mesh of the domain, we define the following data_colle:

.t. ! .t.= formatted .f.=unformatted
2 ! nb of subdomain
. FEM.rect_in
2 ! nb of interfaces to keep
2 3 ! index of interfaces
. FEM.circle_ext

where we specified that we want to keep the interfaces of index 2 and 3 between the cylinder and the circle (so one can later impose boundary or interface conditions).

We can now execute the following commands to generate the global P1-P2 mesh:

colle_p1p2.exe

We answer 0 to the question: "Attention to cavities: number of cavities = ". It generates a P1-P2 mesh called COLLE.FEM that can be used by the code SFEMaNS. The correct generation of the mesh can be checked by visualizing various outputs (such as side_p1.plt, sides.plt, etc.). with plotmtv.

Eventually we symmetrize the mesh with respect to the line \( y=0 \) with the following data_symmetrize:

.t. ! .t.= formatted .f.=unformatted
'.' 'COLLE.FEM' ! directory of the mesh to symmetrize and name of the mesh to symmetrize
0.d0 1.d0 .0d0 ! a,b,c, where ax+by+c=0 is the equation of symmetry line
.t. ! keep sides set to true

and by typing the command:

symmetrize.exe

It generates a symmetrized P1-P2 mesh of COLLE.FEM called symm_COLLE.FEM . The symmetrized mesh can be visualized with the files sides_p1.plt, i_d_p1.plt, etc.

It concludes this example. We note other examples of mesh generation are present in ($SFEMaNS_MESH_GEN_DIR)/EXAMPLES. We refer to the sections Examples with manufactured solutions and Examples on physical problems for more details.