FreeFEM Documentation on GitHub

stars - forks

# Mesh Generation

In this section, operators and tools on meshes are presented.

FreeFEM type for mesh variable:

• 1D mesh: meshL

• 2D mesh: mesh

• 3D volume mesh: mesh3

• 3D border meshes
• 3D surface meshS

• 3D curve meshL

Through this presentation, the principal commands for the mesh generation and links between mesh - mesh3 - meshS - meshL are described.

## The type mesh in 2 dimension

### Commands for 2d mesh Generation

The FreeFEM type to define a 2d mesh object is mesh.

#### The command square

The command square triangulates the unit square.

The following generates a $$4 \times 5$$ grid in the unit square $$[0,1]^2$$. The labels of the boundaries are shown in Fig. 57.

1 mesh Th = square(4, 5);


To construct a $$n\times m$$ grid in the rectangle $$[x_0,x_1]\times [y_0,y_1]$$, proceed as follows:

1 real x0 = 1.2;
2 real x1 = 1.8;
3 real y0 = 0;
4 real y1 = 1;
5 int n = 5;
6 real m = 20;
7 mesh Th = square(n, m, [x0+(x1-x0)*x, y0+(y1-y0)*y]);


Note

Adding the named parameter flags=icase with icase:

1. will produce a mesh where all quads are split with diagonal $$x-y=constant$$

2. will produce a Union Jack flag type of mesh

3. will produce a mesh where all quads are split with diagonal $$x+y=constant$$

4. same as in case 0, except two corners where the triangles are the same as case 2, to avoid having 3 vertices on the boundary

5. same as in case 2, except two corners where the triangles are the same as case 0, to avoid having 3 vertices on the boundary

1 mesh Th = square(n, m, [x0+(x1-x0)*x, y0+(y1-y0)*y], flags=icase);


Note

Adding the named parameter label=labs will change the 4 default label numbers to labs[i-1], for example int[int] labs=[11, 12, 13, 14], and adding the named parameter region=10 will change the region number to $$10$$, for instance (v 3.8).

To see all of these flags at work, check Square mesh example:

1 for (int i = 0; i < 5; ++i){
2    int[int] labs = [11, 12, 13, 14];
3    mesh Th = square(3, 3, flags=i, label=labs, region=10);
4    plot(Th, wait=1, cmm="square flags = "+i );
5 }


#### The command buildmesh

mesh building with border

Boundaries are defined piecewise by parametrized curves. The pieces can only intersect at their endpoints, but it is possible to join more than two endpoints. This can be used to structure the mesh if an area touches a border and create new regions by dividing larger ones:

 1 int upper = 1;
2 int others = 2;
3 int inner = 3;
4
5 border C01(t=0, 1){x=0; y=-1+t; label=upper;}
6 border C02(t=0, 1){x=1.5-1.5*t; y=-1; label=upper;}
7 border C03(t=0, 1){x=1.5; y=-t; label=upper;}
8 border C04(t=0, 1){x=1+0.5*t; y=0; label=others;}
9 border C05(t=0, 1){x=0.5+0.5*t; y=0; label=others;}
10 border C06(t=0, 1){x=0.5*t; y=0; label=others;}
11 border C11(t=0, 1){x=0.5; y=-0.5*t; label=inner;}
12 border C12(t=0, 1){x=0.5+0.5*t; y=-0.5; label=inner;}
13 border C13(t=0, 1){x=1; y=-0.5+0.5*t; label=inner;}
14
15 int n = 10;
16 plot(C01(-n) + C02(-n) + C03(-n) + C04(-n) + C05(-n)
17   + C06(-n) + C11(n) + C12(n) + C13(n), wait=true);
18
19 mesh Th = buildmesh(C01(-n) + C02(-n) + C03(-n) + C04(-n) + C05(-n)
20   + C06(-n) + C11(n) + C12(n) + C13(n));
21
22 plot(Th, wait=true);
23
24 cout << "Part 1 has region number " << Th(0.75, -0.25).region << endl;
25 cout << "Part 2 has redion number " << Th(0.25, -0.25).region << endl;


Borders and mesh are respectively shown in Fig. 58 and Fig. 59.

Border

Triangulation keywords assume that the domain is defined as being on the left (resp right) of its oriented parameterized boundary

$\Gamma_j = \{(x,y)\left|\; x=\varphi_x(t),\, y=\varphi_y(t),\, a_j\le t\le b_j\right.\}$

To check the orientation plot $$t\mapsto (\varphi_x(t),\varphi_y(t)),\, t_0\le t\le t_1$$. If it is as in Fig. 60, then the domain lies on the shaded area, otherwise it lies on the opposite side.

The general expression to define a triangulation with buildmesh is

1 mesh Mesh_Name = buildmesh(Gamma1(m1)+...+GammaJ(mj), OptionalParameter);


where $$m_j$$ are positive or negative numbers to indicate how many vertices should be on $$\Gamma_j,\, \Gamma=\cup_{j=1}^J \Gamma_J$$, and the optional parameter (see also References), separated with a comma, can be:

• nbvx= int, to set the maximum number of vertices in the mesh.

• fixedborder= bool, to say if the mesh generator can change the boundary mesh or not (by default the boundary mesh can change; beware that with periodic boundary conditions (see. Finite Element), it can be dangerous.

The orientation of boundaries can be changed by changing the sign of $$m_j$$.

The following example shows how to change the orientation. The example generates the unit disk with a small circular hole, and assigns “1” to the unit disk (“2” to the circle inside). The boundary label must be non-zero, but it can also be omitted.

1 border a(t=0, 2*pi){x=cos(t); y=sin(t); label=1;}
2 border b(t=0, 2*pi){x=0.3+0.3*cos(t); y=0.3*sin(t); label=2;}
3 plot(a(50) + b(30)); //to see a plot of the border mesh
4 mesh Thwithouthole = buildmesh(a(50) + b(30));
5 mesh Thwithhole = buildmesh(a(50) + b(-30));
6 plot(Thwithouthole, ps="Thwithouthole.eps");
7 plot(Thwithhole, ps="Thwithhole.eps");


Note

Notice that the orientation is changed by b(-30) in the 5th line. In the 7th line, ps="fileName" is used to generate a postscript file with identification shown on the figure.

Mesh with a hole

Note

Borders are evaluated only at the time plot or buildmesh is called so the global variables are defined at this time. In this case, since $$r$$ is changed between the two border calls, the following code will not work because the first border will be computed with r=0.3:

1 real r=1;
2 border a(t=0, 2*pi){x=r*cos(t); y=r*sin(t); label=1;}
3 r=0.3;
4 border b(t=0, 2*pi){x=r*cos(t); y=r*sin(t); label=1;}
5 mesh Thwithhole = buildmesh(a(50) + b(-30)); // bug (a trap) because
6     // the two circles have the same radius = :math:0.3


mesh building with array of border

Sometimes it can be useful to make an array of the border, but unfortunately it is incompatible with the FreeFEM syntax. To bypass this problem, if the number of segments of the discretization $$n$$ is an array, we make an implicit loop on all of the values of the array, and the index variable $$i$$ of the loop is defined after the parameter definition, like in border a(t=0, 2*pi; i)

A first very small example:

1 border a(t=0, 2*pi; i){x=(i+1)*cos(t); y=(i+1)*sin(t); label=1;}
2 int[int] nn = [10, 20, 30];
3 plot(a(nn)); //plot 3 circles with 10, 20, 30 points


And a more complex example to define a square with small circles:

 1 real[int] xx = [0, 1, 1, 0],
2           yy = [0, 0, 1, 1];
3 //radius, center of the 4 circles
4 real[int] RC = [0.1, 0.05, 0.05, 0.1],
5           XC = [0.2, 0.8, 0.2, 0.8],
6           YC = [0.2, 0.8, 0.8, 0.2];
7 int[int] NC = [-10,-11,-12,13]; //list number of :math:\pm segments of the 4 circles borders
8
9 border bb(t=0, 1; i)
10 {
11     // i is the index variable of the multi border loop
12     int ii = (i+1)%4;
13     real t1 = 1-t;
14     x = xx[i]*t1 + xx[ii]*t;
15     y = yy[i]*t1 + yy[ii]*t;
16     label = 0;
17 }
18
19 border cc(t=0, 2*pi; i)
20 {
21     x = RC[i]*cos(t) + XC[i];
22     y = RC[i]*sin(t) + YC[i];
23     label = i + 1;
24 }
25 int[int] nn = [4, 4, 5, 7]; //4 border, with 4, 4, 5, 7 segment respectively
26 plot(bb(nn), cc(NC), wait=1);
27 mesh th = buildmesh(bb(nn) + cc(NC));
28 plot(th, wait=1);


### Mesh Connectivity and data

The following example explains methods to obtain mesh information.

  1 // Mesh
2 mesh Th = square(2, 2);
3
4 cout << "// Get data of the mesh" << endl;
5 {
6     int NbTriangles = Th.nt;
7     real MeshArea = Th.measure;
8     real BorderLength = Th.bordermeasure;
9
10     cout << "Number of triangle(s) = " << NbTriangles << endl;
11     cout << "Mesh area = " << MeshArea << endl;
12     cout << "Border length = " << BorderLength << endl;
13
14     // Th(i) return the vextex i of Th
15     // Th[k] return the triangle k of Th
16     // Th[k][i] return the vertex i of the triangle k of Th
17     for (int i = 0; i < NbTriangles; i++)
18         for (int j = 0; j < 3; j++)
19             cout << i << " " << j << " - Th[i][j] = " << Th[i][j]
20                  << ", x = " << Th[i][j].x
21                  << ", y= " << Th[i][j].y
22                  << ", label=" << Th[i][j].label << endl;
23 }
24
25 cout << "// Hack to get vertex coordinates" << endl;
26 {
27     fespace femp1(Th, P1);
28     femp1 Thx=x,Thy=y;
29
30     int NbVertices = Th.nv;
31     cout << "Number of vertices = " << NbVertices << endl;
32
33     for (int i = 0; i < NbVertices; i++)
34         cout << "Th(" << i << ") : " << Th(i).x << " " << Th(i).y << " " << Th(i).label
35              << endl << "\told method: " << Thx[][i] << " " << Thy[][i] << endl;
36 }
37
38 cout << "// Method to find information of point (0.55,0.6)" << endl;
39 {
40     int TNumber = Th(0.55, 0.6).nuTriangle; //the triangle number
41     int RLabel = Th(0.55, 0.6).region; //the region label
42
43     cout << "Triangle number in point (0.55, 0.6): " << TNumber << endl;
44     cout << "Region label in point (0.55, 0.6): " << RLabel << endl;
45 }
46
47 cout << "// Information of triangle" << endl;
48 {
49     int TNumber = Th(0.55, 0.6).nuTriangle;
50     real TArea = Th[TNumber].area; //triangle area
51     real TRegion = Th[TNumber].region; //triangle region
52     real TLabel = Th[TNumber].label; //triangle label, same as region for triangles
53
54     cout << "Area of triangle " << TNumber << ": " << TArea << endl;
55     cout << "Region of triangle " << TNumber << ": " << TRegion << endl;
56     cout << "Label of triangle " << TNumber << ": " << TLabel << endl;
57 }
58
59 cout << "// Hack to get a triangle containing point x, y or region number (old method)" << endl;
60 {
61     fespace femp0(Th, P0);
62     femp0 TNumbers; //a P0 function to get triangle numbering
63     for (int i = 0; i < Th.nt; i++)
64         TNumbers[][i] = i;
65     femp0 RNumbers = region; //a P0 function to get the region number
66
67     int TNumber = TNumbers(0.55, 0.6); // Number of the triangle containing (0.55, 0,6)
68     int RNumber = RNumbers(0.55, 0.6); // Number of the region containing (0.55, 0,6)
69
70     cout << "Point (0.55,0,6) :" << endl;
71     cout << "\tTriangle number = " << TNumber << endl;
72     cout << "\tRegion number = " << RNumber << endl;
73 }
74
75 cout << "// New method to get boundary information and mesh adjacent" << endl;
76 {
77     int k = 0;
78     int l=1;
79     int e=1;
80
81     // Number of boundary elements
82     int NbBoundaryElements = Th.nbe;
83     cout << "Number of boundary element = " << NbBoundaryElements << endl;
84     // Boundary element k in {0, ..., Th.nbe}
85     int BoundaryElement = Th.be(k);
86     cout << "Boundary element " << k << " = " << BoundaryElement << endl;
87     // Vertice l in {0, 1} of boundary element k
88     int Vertex = Th.be(k)[l];
89     cout << "Vertex " << l << " of boundary element " << k << " = " << Vertex << endl;
90     // Triangle containg the boundary element k
91     int Triangle = Th.be(k).Element;
92     cout << "Triangle containing the boundary element " << k << " = " << Triangle << endl;
93     // Triangle egde nubmer containing the boundary element k
94     int Edge = Th.be(k).whoinElement;
95     cout << "Triangle edge number containing the boundary element " << k << " = " << Edge << endl;
96     // Adjacent triangle of the triangle k by edge e
97     int Adjacent = Th[k].adj(e); //The value of e is changed to the corresponding edge in the adjacent triangle
98     cout << "Adjacent triangle of the triangle " << k << " by edge " << e << " = " << Adjacent << endl;
99     cout << "\tCorresponding edge = " << e << endl;
100     // If there is no adjacent triangle by edge e, the same triangle is returned
102     // Else a different triangle is returned
104 }
105
106 cout << "// Print mesh connectivity " << endl;
107 {
108     int NbTriangles = Th.nt;
109     for (int k = 0; k < NbTriangles; k++)
110         cout << k << " : " << int(Th[k][0]) << " " << int(Th[k][1])
111              << " " << int(Th[k][2])
112              << ", label " << Th[k].label << endl;
113
114     for (int k = 0; k < NbTriangles; k++)
115         for (int e = 0, ee; e < 3; e++)
116             //set ee to e, and ee is change by method adj,
117             cout << k << " " << e << " <=> " << int(Th[k].adj((ee=e))) << " " << ee
119
120     int NbBoundaryElements = Th.nbe;
121     for (int k = 0; k < NbBoundaryElements; k++)
122         cout << k << " : " << Th.be(k)[0] << " " << Th.be(k)[1]
123              << " , label " << Th.be(k).label
124              << ", triangle " << int(Th.be(k).Element)
125              << " " << Th.be(k).whoinElement << endl;
126
127     real[int] bb(4);
128     boundingbox(Th, bb);
129     // bb[0] = xmin, bb[1] = xmax, bb[2] = ymin, bb[3] =ymax
130     cout << "boundingbox:" << endl;
131     cout << "xmin = " << bb[0]
132          << ", xmax = " << bb[1]
133          << ", ymin = " << bb[2]
134          << ", ymax = " << bb[3] << endl;
135 }


The output is:

  1 // Get data of the mesh
2 Number of triangle = 8
3 Mesh area = 1
4 Border length = 4
5 0 0 - Th[i][j] = 0, x = 0, y= 0, label=4
6 0 1 - Th[i][j] = 1, x = 0.5, y= 0, label=1
7 0 2 - Th[i][j] = 4, x = 0.5, y= 0.5, label=0
8 1 0 - Th[i][j] = 0, x = 0, y= 0, label=4
9 1 1 - Th[i][j] = 4, x = 0.5, y= 0.5, label=0
10 1 2 - Th[i][j] = 3, x = 0, y= 0.5, label=4
11 2 0 - Th[i][j] = 1, x = 0.5, y= 0, label=1
12 2 1 - Th[i][j] = 2, x = 1, y= 0, label=2
13 2 2 - Th[i][j] = 5, x = 1, y= 0.5, label=2
14 3 0 - Th[i][j] = 1, x = 0.5, y= 0, label=1
15 3 1 - Th[i][j] = 5, x = 1, y= 0.5, label=2
16 3 2 - Th[i][j] = 4, x = 0.5, y= 0.5, label=0
17 4 0 - Th[i][j] = 3, x = 0, y= 0.5, label=4
18 4 1 - Th[i][j] = 4, x = 0.5, y= 0.5, label=0
19 4 2 - Th[i][j] = 7, x = 0.5, y= 1, label=3
20 5 0 - Th[i][j] = 3, x = 0, y= 0.5, label=4
21 5 1 - Th[i][j] = 7, x = 0.5, y= 1, label=3
22 5 2 - Th[i][j] = 6, x = 0, y= 1, label=4
23 6 0 - Th[i][j] = 4, x = 0.5, y= 0.5, label=0
24 6 1 - Th[i][j] = 5, x = 1, y= 0.5, label=2
25 6 2 - Th[i][j] = 8, x = 1, y= 1, label=3
26 7 0 - Th[i][j] = 4, x = 0.5, y= 0.5, label=0
27 7 1 - Th[i][j] = 8, x = 1, y= 1, label=3
28 7 2 - Th[i][j] = 7, x = 0.5, y= 1, label=3
29 // Hack to get vertex coordinates
30 Number of vertices = 9
31 Th(0) : 0 0 4
32    old method: 0 0
33 Th(1) : 0.5 0 1
34    old method: 0.5 0
35 Th(2) : 1 0 2
36    old method: 1 0
37 Th(3) : 0 0.5 4
38    old method: 0 0.5
39 Th(4) : 0.5 0.5 0
40    old method: 0.5 0.5
41 Th(5) : 1 0.5 2
42    old method: 1 0.5
43 Th(6) : 0 1 4
44    old method: 0 1
45 Th(7) : 0.5 1 3
46    old method: 0.5 1
47 Th(8) : 1 1 3
48    old method: 1 1
49 // Method to find the information of point (0.55,0.6)
50 Triangle number in point (0.55, 0.6): 7
51 Region label in point (0.55, 0.6): 0
52 // Information of a triangle
53 Area of triangle 7: 0.125
54 Region of triangle 7: 0
55 Label of triangle 7: 0
56 // Hack to get a triangle containing point x, y or region number (old method)
57 Point (0.55,0,6) :
58    Triangle number = 7
59    Region number = 0
60 // New method to get boundary information and mesh adjacent
61 Number of boundary element = 8
62 Boundary element 0 = 0
63 Vertex 1 of boundary element 0 = 1
64 Triangle containing the boundary element 0 = 0
65 Triangle edge number containing the boundary element 0 = 2
66 Adjacent triangle of the triangle 0 by edge 1 = 1
67    Corresponding edge = 2
68 // Print mesh connectivity
69 0 : 0 1 4, label 0
70 1 : 0 4 3, label 0
71 2 : 1 2 5, label 0
72 3 : 1 5 4, label 0
73 4 : 3 4 7, label 0
74 5 : 3 7 6, label 0
75 6 : 4 5 8, label 0
76 7 : 4 8 7, label 0
77 0 0 <=> 3 1, adj: 1
78 0 1 <=> 1 2, adj: 1
79 0 2 <=> 0 2, adj: 0
80 1 0 <=> 4 2, adj: 1
81 1 1 <=> 1 1, adj: 0
82 1 2 <=> 0 1, adj: 1
83 2 0 <=> 2 0, adj: 0
84 2 1 <=> 3 2, adj: 1
85 2 2 <=> 2 2, adj: 0
86 3 0 <=> 6 2, adj: 1
87 3 1 <=> 0 0, adj: 1
88 3 2 <=> 2 1, adj: 1
89 4 0 <=> 7 1, adj: 1
90 4 1 <=> 5 2, adj: 1
91 4 2 <=> 1 0, adj: 1
92 5 0 <=> 5 0, adj: 0
93 5 1 <=> 5 1, adj: 0
94 5 2 <=> 4 1, adj: 1
95 6 0 <=> 6 0, adj: 0
96 6 1 <=> 7 2, adj: 1
97 6 2 <=> 3 0, adj: 1
98 7 0 <=> 7 0, adj: 0
99 7 1 <=> 4 0, adj: 1
100 7 2 <=> 6 1, adj: 1
101 0 : 0 1 , label 1, triangle 0 2
102 1 : 1 2 , label 1, triangle 2 2
103 2 : 2 5 , label 2, triangle 2 0
104 3 : 5 8 , label 2, triangle 6 0
105 4 : 6 7 , label 3, triangle 5 0
106 5 : 7 8 , label 3, triangle 7 0
107 6 : 0 3 , label 4, triangle 1 1
108 7 : 3 6 , label 4, triangle 5 1
109 boundingbox:
110 xmin = 0, xmax = 1, ymin = 0, ymax = 1


The real characteristic function of a mesh Th is chi(Th) in 2D and 3D where:

chi(Th)(P)=1 if $$P\in Th$$

chi(Th)(P)=0 if $$P\not\in Th$$

### The keyword “triangulate”

FreeFEM is able to build a triangulation from a set of points. This triangulation is a Delaunay mesh of the convex hull of the set of points. It can be useful to build a mesh from a table function.

The coordinates of the points and the value of the table function are defined separately with rows of the form: x y f(x,y) in a file such as:

1 0.51387 0.175741 0.636237
2 0.308652 0.534534 0.746765
3 0.947628 0.171736 0.899823
4 0.702231 0.226431 0.800819
5 0.494773 0.12472 0.580623
6 0.0838988 0.389647 0.456045
7 ...............


Triangulate

The third column of each line is left untouched by the triangulate command. But you can use this third value to define a table function with rows of the form: x y f(x,y).

The following example shows how to make a mesh from the file xyf with the format stated just above. The command triangulate only uses the 1st and 2nd columns.

 1 // Build the Delaunay mesh of the convex hull
2 mesh Thxy=triangulate("xyf"); //points are defined by the first 2 columns of file xyf
3
4 // Plot the created mesh
5 plot(Thxy);
6
7 // Fespace
8 fespace Vhxy(Thxy, P1);
9 Vhxy fxy;
10
11 // Reading the 3rd column to define the function fxy
12 {
13     ifstream file("xyf");
14     real xx, yy;
15     for(int i = 0; i < fxy.n; i++)
16         file >> xx >> yy >> fxy[][i]; //to read third row only.
17                                       //xx and yy are just skipped
18 }
19
20 // Plot
21 plot(fxy);


One new way to build a mesh is to have two arrays: one for the $$x$$ values and the other for the $$y$$ values.

1 //set two arrays for the x's and y's
2 Vhxy xx=x, yy=y;
3 //build the mesh
4 mesh Th = triangulate(xx[], yy[]);


### 2d Finite Element space on a boundary

To define a Finite Element space on a boundary, we came up with the idea of a mesh with no internal points (called empty mesh). It can be useful to handle Lagrange multipliers in mixed and mortar methods.

So the function emptymesh removes all the internal points of a mesh except points on internal boundaries.

1 {
2     border a(t=0, 2*pi){x=cos(t); y=sin(t); label=1;}
3     mesh Th = buildmesh(a(20));
4     Th = emptymesh(Th);
5     plot(Th);
6 }


It is also possible to build an empty mesh of a pseudo subregion with emptymesh(Th, ssd) using the set of edges from the mesh Th; an edge $$e$$ is in this set when, with the two adjacent triangles $$e =t1\cap t2$$ and $$ssd[T1] \neq ssd[T2]$$ where $$ssd$$ refers to the pseudo region numbering of triangles, they are stored in the int[int] array of size “the number of triangles”.

 1 {
2     mesh Th = square(10, 10);
3     int[int] ssd(Th.nt);
4     //build the pseudo region numbering
5     for(int i = 0; i < ssd.n; i++){
6         int iq = i/2; //because 2 triangles per quad
7         int ix = iq%10;
8         int iy = iq/10;
9         ssd[i] = 1 + (ix>=5) + (iy>=5)*2;
10     }
11     //build emtpy with all edges $e=T1 \cap T2$ and $ssd[T1] \neq ssd[T2]$
12     Th = emptymesh(Th, ssd);
13     //plot
14     plot(Th);
15     savemesh(Th, "emptymesh.msh");
16 }


Empty mesh

### Remeshing

#### The command movemesh

Meshes can be translated, rotated, and deformed by movemesh; this is useful for elasticity to watch the deformation due to the displacement $$\mathbf{\Phi}(x,y)=(\Phi_1(x,y),\Phi_2(x,y))$$ of shape.

It is also useful to handle free boundary problems or optimal shape problems.

If $$\Omega$$ is triangulated as $$T_h(\Omega)$$, and $$\mathbf{\Phi}$$ is a displacement vector then $$\mathbf{\Phi}(T_h)$$ is obtained by:

1 mesh Th = movemesh(Th,[Phi1, Phi2]);


Sometimes the transformed mesh is invalid because some triangles have flipped over (meaning it now has a negative area). To spot such problems, one may check the minimum triangle area in the transformed mesh with checkmovemesh before any real transformation.

For example:

$\begin{split}\begin{array}{rcl} \Phi_1(x,y) &=& x+k*\sin(y*\pi)/10)\\ \Phi_2(x,y) &=& y+k*\cos(y\pi)/10) \end{array}\end{split}$

for a big number $$k>1$$.

 1 verbosity = 4;
2
3 // Parameters
4 real coef = 1;
5
6 // Mesh
7 border a(t=0, 1){x=t; y=0; label=1;};
8 border b(t=0, 0.5){x=1; y=t; label=1;};
9 border c(t=0, 0.5){x=1-t; y=0.5; label=1;};
10 border d(t=0.5, 1){x=0.5; y=t; label=1;};
11 border e(t=0.5, 1){x=1-t; y=1; label=1;};
12 border f(t=0, 1){x=0; y=1-t; label=1;};
13 mesh Th = buildmesh(a(6) + b(4) + c(4) + d(4) + e(4) + f(6));
14 plot(Th, wait=true, fill=true, ps="Lshape.eps");
15
16 // Function
17 func uu = sin(y*pi)/10;
18 func vv = cos(x*pi)/10;
19
20 // Checkmovemesh
21 real minT0 = checkmovemesh(Th, [x, y]); //return the min triangle area
22 while(1){ // find a correct move mesh
23     real minT = checkmovemesh(Th, [x+coef*uu, y+coef*vv]);
24     if (minT > minT0/5) break; //if big enough
25     coef /= 1.5;
26 }
27
28 // Movemesh
29 Th = movemesh(Th, [x+coef*uu, y+coef*vv]);
30 plot(Th, wait=true, fill=true, ps="MovedMesh.eps");


Move mesh

Note

Consider a function $$u$$ defined on a mesh Th. A statement like Th=movemesh(Th...) does not change $$u$$ and so the old mesh still exists. It will be destroyed when no function uses it. A statement like $$u=u$$ redefines $$u$$ on the new mesh Th with interpolation and therefore destroys the old Th, if $$u$$ was the only function using it.

Now, we give an example of moving a mesh with a Lagrangian function $$u$$ defined on the moving mesh.

 1 // Parameters
2 int nn = 10;
3 real dt = 0.1;
4
5 // Mesh
6 mesh Th = square(nn, nn);
7
8 // Fespace
9 fespace Vh(Th, P1);
10 Vh u=y;
11
12 // Loop
13 real t=0;
14 for (int i = 0; i < 4; i++){
15     t = i*dt;
16     Vh f=x*t;
17     real minarea = checkmovemesh(Th, [x, y+f]);
18     if (minarea > 0) //movemesh will be ok
19     Th = movemesh(Th, [x, y+f]);
20
21     cout << " Min area = " << minarea << endl;
22
23     real[int] tmp(u[].n);
24     tmp = u[]; //save the value
25     u = 0;//to change the FEspace and mesh associated with u
26     u[] = tmp;//set the value of u without any mesh update
27     plot(Th, u, wait=true);
28 }
29 // In this program, since u is only defined on the last mesh, all the
30 // previous meshes are deleted from memory.


#### The command hTriangle

This section presents the way to obtain a regular triangulation with FreeFEM.

For a set $$S$$, we define the diameter of $$S$$ by

$\textrm{diam}(S)=\sup\{|\mathbf{x}-\mathbf{y}|; \; \mathbf{x},\, \mathbf{y}\in S\}$

The sequence $$\{\mathcal{T}_h\}_{h\rightarrow 0}$$ of $$\Omega$$ is called regular if they satisfy the following:

1. $$\lim_{h\rightarrow 0}\max\{\textrm{diam}(T_k)|\; T_k\in \mathcal{T}_h\}=0$$

2. There is a number $$\sigma>0$$ independent of $$h$$ such that $$\frac{\rho(T_k)}{\textrm{diam}(T_k)}\ge \sigma\quad \textrm{for all }T_k\in \mathcal{T}_h$$ where $$\rho(T_k)$$ are the diameter of the inscribed circle of $$T_k$$.

We put $$h(\mathcal{T}_h)=\max\{\textrm{diam}(T_k)|\; T_k\in \mathcal{T}_h\}$$, which is obtained by

1 mesh Th = ......;
2 fespace Ph(Th, P0);
3 Ph h = hTriangle;
4 cout << "size of mesh = " << h[].max << endl;


The function:

$f(x,y) = 10.0x^3+y^3+\tan^{-1}[\varepsilon/(\sin(5.0y)-2.0x)],\ \varepsilon = 0.0001$

sharply varies in value and the initial mesh given by one of the commands in the Mesh Generation part cannot reflect its sharp variations.

 1 // Parameters
2 real eps = 0.0001;
3 real h = 1;
4 real hmin = 0.05;
5 func f = 10.0*x^3 + y^3 + h*atan2(eps, sin(5.0*y)-2.0*x);
6
7 // Mesh
8 mesh Th = square(5, 5, [-1+2*x, -1+2*y]);
9
10 // Fespace
11 fespace Vh(Th,P1);
12 Vh fh = f;
13 plot(fh);
14
16 for (int i = 0; i < 2; i++){
18     fh = f; //old mesh is deleted
19     plot(Th, fh, wait=true);
20 }


FreeFEM uses a variable metric/Delaunay automatic meshing algorithm.

The command:

1 mesh ATh = adaptmesh(Th, f);


create the new mesh ATh adapted to the Hessian

$D^2f=(\partial^2 f/\partial x^2,\, \partial^2 f/\partial x\partial y, \partial^2 f/\partial y^2)$

of a function (formula or FE-function).

Mesh adaptation is a very powerful tool when the solution of a problem varies locally and sharply.

Here we solve the Poisson’s problem, when $$f=1$$ and $$\Omega$$ is an L-shape domain.

Tip

The solution has the singularity $$r^{3/2},\, r=|x-\gamma|$$ at the point $$\gamma$$ of the intersection of two lines $$bc$$ and $$bd$$ (see Fig. 70).

 1 // Parameters
2 real error = 0.1;
3
4 // Mesh
5 border ba(t=0, 1){x=t; y=0; label=1;}
6 border bb(t=0, 0.5){x=1; y=t; label=1;}
7 border bc(t=0, 0.5){x=1-t; y=0.5; label=1;}
8 border bd(t=0.5, 1){x=0.5; y=t; label=1;}
9 border be(t=0.5, 1){x=1-t; y=1; label=1;}
10 border bf(t=0, 1){x=0; y=1-t; label=1;}
11 mesh Th = buildmesh(ba(6) + bb(4) + bc(4) + bd(4) + be(4) + bf(6));
12
13 // Fespace
14 fespace Vh(Th, P1);
15 Vh u, v;
16
17 // Function
18 func f = 1;
19
20 // Problem
21 problem Poisson(u, v, solver=CG, eps=1.e-6)
22     = int2d(Th)(
23           dx(u)*dx(v)
24         + dy(u)*dy(v)
25     )
26     - int2d(Th)(
27           f*v
28     )
29     + on(1, u=0);
30
32 for (int i = 0; i < 4; i++){
33     Poisson;
34     Th = adaptmesh(Th, u, err=error);
35     error = error/2;
36 }
37
38 // Plot
39 plot(u);


To speed up the adaptation, the default parameter err of adaptmesh is changed by hand; it specifies the required precision, so as to make the new mesh finer or coarser.

The problem is coercive and symmetric, so the linear system can be solved with the conjugate gradient method (parameter solver=CG) with the stopping criteria on the residual, here eps=1.e-6).

By adaptmesh, the slope of the final solution is correctly computed near the point of intersection of $$bc$$ and $$bd$$ as in Fig. 71.

This method is described in detail in [HECHT1998]. It has a number of default parameters which can be modified.

If f1,f2 are functions and thold, Thnew are meshes:

1    Thnew = adaptmesh(Thold, f1 ... );
2    Thnew = adaptmesh(Thold, f1,f2 ... ]);
3    Thnew = adaptmesh(Thold, [f1,f2] ... );


The additional parameters of adaptmesh are:

• hmin= Minimum edge size.

Its default is related to the size of the domain to be meshed and the precision of the mesh generator.

• hmax= Maximum edge size.

It defaults to the diameter of the domain to be meshed.

• err= $$P_1$$ interpolation error level (0.01 is the default).

• errg= Relative geometrical error.

By default this error is 0.01, and in any case it must be lower than $$1/\sqrt{2}$$. Meshes created with this option may have some edges smaller than the -hmin due to geometrical constraints.

• nbvx= Maximum number of vertices generated by the mesh generator (9000 is the default).

• nbsmooth= number of iterations of the smoothing procedure (5 is the default).

• nbjacoby= number of iterations in a smoothing procedure during the metric construction, 0 means no smoothing, 6 is the default.

• ratio= ratio for a prescribed smoothing on the metric.

If the value is 0 or less than 1.1 no smoothing is done on the metric. 1.8 is the default. If ratio > 1.1, the speed of mesh size variations is bounded by $$log(\mathtt{ratio})$$.

Note

As ratio gets closer to 1, the number of generated vertices increases. This may be useful to control the thickness of refined regions near shocks or boundary layers.

• omega= relaxation parameter for the smoothing procedure. 1.0 is the default.

• iso= If true, forces the metric to be isotropic. false is the default.

• abserror= If false, the metric is evaluated using the criteria of equi-repartion of relative error.

false is the default. In this case the metric is defined by:

$\mathcal{M} = \left({1\over\mathtt{err}\,\, \mathtt{coef}^2} \quad { |\mathcal{H}| \over max(\mathtt{CutOff},|\eta|)}\right)^p$

Otherwise, the metric is evaluated using the criteria of equi-distribution of errors. In this case the metric is defined by:

$\mathcal{M} = \left({1\over \mathtt{err}\,\,\mathtt{coef}^2} \quad {|{\mathcal{H}|} \over {\sup(\eta)-\inf(\eta)}}\right)^p.\label{eq err abs}$
• cutoff= lower limit for the relative error evaluation. 1.0e-6 is the default.

• verbosity= informational messages level (can be chosen between 0 and $$\infty$$).

Also changes the value of the global variable verbosity (obsolete).

• inquire= To inquire graphically about the mesh. false is the default.

• splitpbedge= If true, splits all internal edges in half with two boundary vertices.

true is the default.

• maxsubdiv= Changes the metric such that the maximum subdivision of a background edge is bound by val.

Always limited by 10, and 10 is also the default.

• rescaling= if true, the function, with respect to which the mesh is adapted, is rescaled to be between 0 and 1.

true is the default.

• keepbackvertices= if true, tries to keep as many vertices from the original mesh as possible.

true is the default.

• IsMetric= if true, the metric is defined explicitly.

false is the default. If the 3 functions $$m_{11}, m_{12}, m_{22}$$ are given, they directly define a symmetric matrix field whose Hessian is computed to define a metric. If only one function is given, then it represents the isotropic mesh size at every point.

For example, if the partial derivatives fxx ($$=\partial^2 f/\partial x^2$$), fxy ($$=\partial^2 f/\partial x\partial y$$), fyy ($$=\partial^2 f/\partial y^2$$) are given, we can set Th = adaptmesh(Th, fxx, fxy, fyy, IsMetric=1, nbvx=10000, hmin=hmin);

• power= exponent power of the Hessian used to compute the metric.

1 is the default.

• thetamax= minimum corner angle in degrees.

Default is $$10^\circ$$ where the corner is $$ABC$$ and the angle is the angle of the two vectors $${AB}, {BC}$$, ($$0$$ imply no corner, $$90$$ imply perpendicular corner, …).

• splitin2= boolean value.

If true, splits all triangles of the final mesh into 4 sub-triangles.

• metric= an array of 3 real arrays to set or get metric data information.

The size of these three arrays must be the number of vertices. So if m11,m12,m22 are three P1 finite elements related to the mesh to adapt, you can write: metric=[m11[],m12[],m22[]] (see file convect-apt.edp for a full example)

• nomeshgeneration= If true, no adapted mesh is generated (useful to compute only a metric).

• periodic= Writing periodic=[[4,y],[2,y],[1,x],[3,x]]; builds an adapted periodic mesh.

The sample builds a biperiodic mesh of a square. (see periodic finite element spaces, and see the Sphere example for a full example)

We can use the command adaptmesh to build a uniform mesh with a constant mesh size. To build a mesh with a constant mesh size equal to $$\frac{1}{30}$$ try:

1 mesh Th=square(2, 2); //the initial mesh
2 plot(Th, wait=true, ps="square-0.eps");
3
4 Th = adaptmesh(Th, 1./30., IsMetric=1, nbvx=10000);
5 plot(Th, wait=true, ps="square-1.eps");
6
7 Th = adaptmesh(Th, 1./30., IsMetric=1, nbvx=10000); //More the one time du to
8 Th = adaptmesh(Th, 1./30., IsMetric=1, nbvx=10000); //Adaptation bound maxsubdiv=
9 plot(Th, wait=true, ps="square-2.eps");


#### The command trunc

Two operators have been introduced to remove triangles from a mesh or to divide them. Operator trunc has the following parameters:

• boolean function to keep or remove elements

• label= sets the label number of new boundary item, one by default.

• split= sets the level $$n$$ of triangle splitting.

Each triangle is split in $$n\times n$$, one by default.

To create the mesh Th3 where all triangles of a mesh Th are split in $$3{\times}3$$, just write:

1 mesh Th3 = trunc(Th, 1, split=3);


The following example construct all “trunced” meshes to the support of the basic function of the space Vh (cf. abs(u)>0), split all the triangles in $$5{\times} 5$$, and put a label number to $$2$$ on a new boundary.

 1 // Mesh
2 mesh Th = square(3, 3);
3
4 // Fespace
5 fespace Vh(Th, P1);
6 Vh u=0;
7
8 // Loop on all degrees of freedom
9 int n=u.n;
10 for (int i = 0; i < n; i++){
11     u[][i] = 1; // The basis function i
12     plot(u, wait=true);
13     mesh Sh1 = trunc(Th, abs(u)>1.e-10, split=5, label=2);
14     plot(Th, Sh1, wait=true, ps="trunc"+i+".eps");
15     u[][i] = 0; // reset
16 }


Trunc

#### The command change

This command changes the label of elements and border elements of a mesh.

Changing the label of elements and border elements will be done using the keyword change. The parameters for this command line are for two dimensional and three dimensional cases:

• refe= is an array of integers to change the references on edges

• reft= is an array of integers to change the references on triangles

• label= is an array of integers to change the 4 default label numbers

• region= is an array of integers to change the default region numbers

• renumv= is an array of integers, which explicitly gives the new numbering of vertices in the new mesh. By default, this numbering is that of the original mesh

• renumt= is an array of integers, which explicitly gives the new numbering of elements in the new mesh, according the new vertices numbering given by renumv=. By default, this numbering is that of the original mesh

• flabel= is an integer function given the new value of the label

• fregion= is an integer function given the new value of the region

• rmledges= is an integer to remove edges in the new mesh, following a label

• rmInternalEdges= is a boolean, if equal to true to remove the internal edges. By default, the internal edges are stored

These vectors are composed of $$n_{l}$$ successive pairs of numbers $$O,N$$ where $$n_{l}$$ is the number (label or region) that we want to change. For example, we have :

(8)$\begin{split}\mathtt{label} &= [ O_{1}, N_{1}, ..., O_{n_{l}},N_{n_{l}} ] \\ \mathtt{region} &= [ O_{1}, N_{1}, ..., O_{n_{l}},N_{n_{l}} ]\end{split}$

An application example is given here:

 1 // Mesh
2 mesh Th1 = square(10, 10);
3 mesh Th2 = square(20, 10, [x+1, y]);
4
5 int[int] r1=[2,0];
6 plot(Th1, wait=true);
7
8 Th1 = change(Th1, label=r1); //change the label of Edges 2 in 0.
9 plot(Th1, wait=true);
10
11 // boundary label: 1 -> 1 bottom, 2 -> 1 right, 3->1 top, 4->1 left boundary label is 1
12 int[int] re=[1,1, 2,1, 3,1, 4,1]
13 Th2=change(Th2,refe=re);
14 plot(Th2,wait=1) ;


#### The command splitmesh

Another way to split mesh triangles is to use splitmesh, for example:

1 // Mesh
2 border a(t=0, 2*pi){x=cos(t); y=sin(t); label=1;}
3 mesh Th = buildmesh(a(20));
4 plot(Th, wait=true, ps="NotSplittedMesh.eps");
5
6 // Splitmesh
7 Th = splitmesh(Th, 1 + 5*(square(x-0.5) + y*y));
8 plot(Th, wait=true, ps="SplittedMesh.eps");


Split mesh

### Meshing Examples

Tip

Two rectangles touching by a side

 1 border a(t=0, 1){x=t; y=0;};
2 border b(t=0, 1){x=1; y=t;};
3 border c(t=1, 0){x=t; y=1;};
4 border d(t=1, 0){x=0; y=t;};
5 border c1(t=0, 1){x=t; y=1;};
6 border e(t=0, 0.2){x=1; y=1+t;};
7 border f(t=1, 0){x=t; y=1.2;};
8 border g(t=0.2, 0){x=0; y=1+t;};
9 int n=1;
10 mesh th = buildmesh(a(10*n) + b(10*n) + c(10*n) + d(10*n));
11 mesh TH = buildmesh(c1(10*n) + e(5*n) + f(10*n) + g(5*n));
12 plot(th, TH, ps="TouchSide.esp");


Tip

NACA0012 Airfoil

1 border upper(t=0, 1){x=t; y=0.17735*sqrt(t) - 0.075597*t - 0.212836*(t^2) + 0.17363*(t^3) - 0.06254*(t^4);}
2 border lower(t=1, 0){x = t; y=-(0.17735*sqrt(t) -0.075597*t - 0.212836*(t^2) + 0.17363*(t^3) - 0.06254*(t^4));}
3 border c(t=0, 2*pi){x=0.8*cos(t) + 0.5; y=0.8*sin(t);}
4 mesh Th = buildmesh(c(30) + upper(35) + lower(35));
5 plot(Th, ps="NACA0012.eps", bw=true);


Tip

Cardioid

1 real b = 1, a = b;
2 border C(t=0, 2*pi){x=(a+b)*cos(t)-b*cos((a+b)*t/b); y=(a+b)*sin(t)-b*sin((a+b)*t/b);}
3 mesh Th = buildmesh(C(50));
4 plot(Th, ps="Cardioid.eps", bw=true);


Tip

Cassini Egg

1 border C(t=0, 2*pi) {x=(2*cos(2*t)+3)*cos(t); y=(2*cos(2*t)+3)*sin(t);}
2 mesh Th = buildmesh(C(50));
3 plot(Th, ps="Cassini.eps", bw=true);


Tip

By cubic Bezier curve

 1 // A cubic Bezier curve connecting two points with two control points
2 func real bzi(real p0, real p1, real q1, real q2, real t){
3     return p0*(1-t)^3 + q1*3*(1-t)^2*t + q2*3*(1-t)*t^2 + p1*t^3;
4 }
5
6 real[int] p00 = [0, 1], p01 = [0, -1], q00 = [-2, 0.1], q01 = [-2, -0.5];
7 real[int] p11 = [1,-0.9], q10 = [0.1, -0.95], q11=[0.5, -1];
8 real[int] p21 = [2, 0.7], q20 = [3, -0.4], q21 = [4, 0.5];
9 real[int] q30 = [0.5, 1.1], q31 = [1.5, 1.2];
10 border G1(t=0, 1){
11     x=bzi(p00[0], p01[0], q00[0], q01[0], t);
12     y=bzi(p00[1], p01[1], q00[1], q01[1], t);
13 }
14 border G2(t=0, 1){
15     x=bzi(p01[0], p11[0], q10[0], q11[0], t);
16     y=bzi(p01[1], p11[1], q10[1], q11[1], t);
17 }
18 border G3(t=0, 1){
19     x=bzi(p11[0], p21[0], q20[0], q21[0], t);
20     y=bzi(p11[1], p21[1], q20[1], q21[1], t);
21 }
22 border G4(t=0, 1){
23     x=bzi(p21[0], p00[0], q30[0], q31[0], t);
24     y=bzi(p21[1], p00[1], q30[1], q31[1], t);
25 }
26 int m = 5;
27 mesh Th = buildmesh(G1(2*m) + G2(m) + G3(3*m) + G4(m));
28 plot(Th, ps="Bezier.eps", bw=true);


Tip

Section of Engine

 1 real a = 6., b = 1., c = 0.5;
2
3 border L1(t=0, 1){x=-a; y=1+b-2*(1+b)*t;}
4 border L2(t=0, 1){x=-a+2*a*t; y=-1-b*(x/a)*(x/a)*(3-2*abs(x)/a );}
5 border L3(t=0, 1){x=a; y=-1-b+(1+b)*t; }
6 border L4(t=0, 1){x=a-a*t; y=0;}
7 border L5(t=0, pi){x=-c*sin(t)/2; y=c/2-c*cos(t)/2;}
8 border L6(t=0, 1){x=a*t; y=c;}
9 border L7(t=0, 1){x=a; y=c+(1+b-c)*t;}
10 border L8(t=0, 1){x=a-2*a*t; y=1+b*(x/a)*(x/a)*(3-2*abs(x)/a);}
11 mesh Th = buildmesh(L1(8) + L2(26) + L3(8) + L4(20) + L5(8) + L6(30) + L7(8) + L8(30));
12 plot(Th, ps="Engine.eps", bw=true);


Tip

Domain with U-shape channel

 1 real d = 0.1; //width of U-shape
2 border L1(t=0, 1-d){x=-1; y=-d-t;}
3 border L2(t=0, 1-d){x=-1; y=1-t;}
4 border B(t=0, 2){x=-1+t; y=-1;}
5 border C1(t=0, 1){x=t-1; y=d;}
6 border C2(t=0, 2*d){x=0; y=d-t;}
7 border C3(t=0, 1){x=-t; y=-d;}
8 border R(t=0, 2){x=1; y=-1+t;}
9 border T(t=0, 2){x=1-t; y=1;}
10 int n = 5;
11 mesh Th = buildmesh(L1(n/2) + L2(n/2) + B(n) + C1(n) + C2(3) + C3(n) + R(n) + T(n));
12 plot(Th, ps="U-shape.eps", bw=true);


Tip

Domain with V-shape cut

 1 real dAg = 0.02; //angle of V-shape
2 border C(t=dAg, 2*pi-dAg){x=cos(t); y=sin(t);};
3 real[int] pa(2), pb(2), pc(2);
4 pa[0] = cos(dAg);
5 pa[1] = sin(dAg);
6 pb[0] = cos(2*pi-dAg);
7 pb[1] = sin(2*pi-dAg);
8 pc[0] = 0;
9 pc[1] = 0;
10 border seg1(t=0, 1){x=(1-t)*pb[0]+t*pc[0]; y=(1-t)*pb[1]+t*pc[1];};
11 border seg2(t=0, 1){x=(1-t)*pc[0]+t*pa[0]; y=(1-t)*pc[1]+t*pa[1];};
12 mesh Th = buildmesh(seg1(20) + C(40) + seg2(20));
13 plot(Th, ps="V-shape.eps", bw=true);


Tip

Smiling face

 1 real d=0.1; int m = 5; real a = 1.5, b = 2, c = 0.7, e = 0.01;
2
3 border F(t=0, 2*pi){x=a*cos(t); y=b*sin(t);}
4 border E1(t=0, 2*pi){x=0.2*cos(t)-0.5; y=0.2*sin(t)+0.5;}
5 border E2(t=0, 2*pi){x=0.2*cos(t)+0.5; y=0.2*sin(t)+0.5;}
6 func real st(real t){
7     return sin(pi*t) - pi/2;
8 }
9 border C1(t=-0.5, 0.5){x=(1-d)*c*cos(st(t)); y=(1-d)*c*sin(st(t));}
10 border C2(t=0, 1){x=((1-d)+d*t)*c*cos(st(0.5)); y=((1-d)+d*t)*c*sin(st(0.5));}
11 border C3(t=0.5, -0.5){x=c*cos(st(t)); y=c*sin(st(t));}
12 border C4(t=0, 1){x=(1-d*t)*c*cos(st(-0.5)); y=(1-d*t)*c*sin(st(-0.5));}
13 border C0(t=0, 2*pi){x=0.1*cos(t); y=0.1*sin(t);}
14
15 mesh Th=buildmesh(F(10*m) + C1(2*m) + C2(3) + C3(2*m) + C4(3)
16     + C0(m) + E1(-2*m) + E2(-2*m));
17 plot(Th, ps="SmileFace.eps", bw=true);


Tip

3 points bending

 1 // Square for Three-Point Bend Specimens fixed on Fix1, Fix2
3 real a = 1, b = 5, c = 0.1;
4 int n = 5, m = b*n;
5 border Left(t=0, 2*a){x=-b; y=a-t;}
6 border Bot1(t=0, b/2-c){x=-b+t; y=-a;}
7 border Fix1(t=0, 2*c){x=-b/2-c+t; y=-a;}
8 border Bot2(t=0, b-2*c){x=-b/2+c+t; y=-a;}
9 border Fix2(t=0, 2*c){x=b/2-c+t; y=-a;}
10 border Bot3(t=0, b/2-c){x=b/2+c+t; y=-a;}
11 border Right(t=0, 2*a){x=b; y=-a+t;}
12 border Top1(t=0, b-c){x=b-t; y=a;}
14 border Top2(t=0, b-c){x=-c-t; y=a;}
15 mesh Th = buildmesh(Left(n) + Bot1(m/4) + Fix1(5) + Bot2(m/2)
16     + Fix2(5) + Bot3(m/4) + Right(n) + Top1(m/2) + Load(10) + Top2(m/2));
17 plot(Th, ps="ThreePoint.eps", bw=true);


## The type mesh3 in 3 dimension

Note

Up to the version 3, FreeFEM allowed to consider a surface problem such as the PDE is treated like boundary conditions on the boundary domain (on triangles describing the boundary domain). With the version 4, in particular 4.2.1, a completed model for surface problem is possible, with the definition of a surface mesh and a surface problem with a variational form on domain ( with triangle elements) and application of boundary conditions on border domain (describing by edges). The keywords to define a surface mesh is meshS.

### 3d mesh generation

Note

For 3D mesh tools, put load "msh3" at the top of the .edp script.

#### The command cube

The function cube like its 2d function square is a simple way to build cubic objects, it is contained in plugin msh3 (import with load "msh3").

The following code generates a $$3\times 4 \times 5$$ grid in the unit cube $$[0, 1]^3$$.

1 mesh3 Th = cube(3, 4, 5);


By default the labels are :

1. face $$y=0$$,

2. face $$x=1$$,

3. face $$y=1$$,

4. face $$x=0$$,

5. face $$z=0$$,

6. face $$z=1$$

and the region number is $$0$$.

A full example of this function to build a mesh of cube $$]-1,1[^3$$ with face label given by $$(ix + 4*(iy+1) + 16*(iz+1))$$ where $$(ix, iy, iz)$$ are the coordinates of the barycenter of the current face, is given below.

 1 load "msh3"
2
3 int[int] l6 = [37, 42, 45, 40, 25, 57];
4 int r11 = 11;
5 mesh3 Th = cube(4, 5, 6, [x*2-1, y*2-1, z*2-1], label=l6, flags =3, region=r11);
6
7 cout << "Volume = " << Th.measure << ", border area = " << Th.bordermeasure << endl;
8
9 int err = 0;
10 for(int i = 0; i < 100; ++i){
11     real s = int2d(Th,i)(1.);
12     real sx = int2d(Th,i)(x);
13     real sy = int2d(Th,i)(y);
14     real sz = int2d(Th,i)(z);
15
16     if(s){
17         int ix = (sx/s+1.5);
18         int iy = (sy/s+1.5);
19         int iz = (sz/s+1.5);
20         int ii = (ix + 4*(iy+1) + 16*(iz+1) );
21         //value of ix,iy,iz => face min 0, face max 2, no face 1
22         cout << "Label = " << i << ", s = " << s << " " << ix << iy << iz << " : " << ii << endl;
23         if( i != ii ) err++;
24     }
25 }
26 real volr11 = int3d(Th,r11)(1.);
27 cout << "Volume region = " << 11 << ": " << volr11 << endl;
28 if((volr11 - Th.measure )>1e-8) err++;
29 plot(Th, fill=false);
30 cout << "Nb err = " << err << endl;
31 assert(err==0);


The output of this script is:

 1 Enter: BuildCube: 3
2   kind = 3 n tet Cube = 6 / n slip 6 19
3 Cube  nv=210 nt=720 nbe=296
4 Out:  BuildCube
5 Volume = 8, border area = 24
6 Label = 25, s = 4 110 : 25
7 Label = 37, s = 4 101 : 37
8 Label = 40, s = 4 011 : 40
9 Label = 42, s = 4 211 : 42
10 Label = 45, s = 4 121 : 45
11 Label = 57, s = 4 112 : 57
12 Volume region = 11: 8
13 Nb err = 0


#### The command buildlayers

This mesh is obtained by extending a two dimensional mesh in the $$z$$-axis.

The domain $$\Omega_{3d}$$ defined by the layer mesh is equal to $$\Omega_{3d} = \Omega_{2d} \times [zmin, zmax]$$ where $$\Omega_{2d}$$ is the domain defined by the two dimensional meshes. $$zmin$$ and $$zmax$$ are functions of $$\Omega_{2d}$$ in $$\R$$ that defines respectively the lower surface and upper surface of $$\Omega_{3d}$$.

For a vertex of a two dimensional mesh $$V_{i}^{2d} = (x_{i},y_{i})$$, we introduce the number of associated vertices in the $$z-$$axis $$M_{i}+1$$.

We denote by $$M$$ the maximum of $$M_{i}$$ over the vertices of the two dimensional mesh. This value is called the number of layers (if $$\forall i, \; M_{i}=M$$ then there are $$M$$ layers in the mesh of $$\Omega_{3d}$$). $$V_{i}^{2d}$$ generated $$M+1$$ vertices which are defined by:

$\forall j=0, \ldots, M, \quad V_{i,j}^{3d} = ( x_{i}, y_{i}, \theta_{i}(z_{i,j}) ),$

where $$(z_{i,j})_{j=0,\ldots,M}$$ are the $$M+1$$ equidistant points on the interval $$[zmin( V_{i}^{2d} ), zmax( V_{i}^{2d})]$$:

$z_{i,j} = j \: \delta \alpha + zmin(V_{i}^{2d}), \quad \delta \alpha= \frac{ zmax( V_{i}^{2d} ) - zmin( V_{i}^{2d}) }{M}.$

The function $$\theta_{i}$$, defined on $$[zmin( V_{i}^{2d} ), zmax( V_{i}^{2d} )]$$, is given by:

$\begin{split}\theta_{i}(z) = \left \{ \begin{array}{cl} \theta_{i,0} & \mbox{if} \: z=zmin(V_{i}^{2d}), \\ \theta_{i,j} & \mbox{if} \: z \in ] \theta_{i,j-1}, \theta_{i,j}],\\ \end{array} \right.\end{split}$

with $$(\theta_{i,j})_{j=0,\ldots,M_{i}}$$ are the $$M_{i}+1$$ equidistant points on the interval $$[zmin( V_{i}^{2d} ), zmax( V_{i}^{2d} )]$$.

Set a triangle $$K=(V_{i1}^{2d}$$, $$V_{i2}^{2d}$$, $$V_{i3}^{2d})$$ of the two dimensional mesh. $$K$$ is associated with a triangle on the upper surface (resp. on the lower surface) of layer mesh:

$$( V_{i1,M}^{3d}, V_{i2,M}^{3d}, V_{i3,M}^{3d} )$$ (resp. $$( V_{i1,0}^{3d}, V_{i2,0}^{3d}, V_{i3,0}^{3d})$$).

Also $$K$$ is associated with $$M$$ volume prismatic elements which are defined by:

$\forall j=0,\ldots,M, \quad H_{j} = ( V_{i1,j}^{3d}, V_{i2,j}^{3d}, V_{i3,j}^{3d}, V_{i1,j+1}^{3d}, V_{i2,j+1}^{3d}, V_{i3,j+1}^{3d} ).$

Theses volume elements can have some merged point:

• 0 merged point : prism

• 1 merged points : pyramid

• 2 merged points : tetrahedra

• 3 merged points : no elements

The elements with merged points are called degenerate elements. To obtain a mesh with tetrahedra, we decompose the pyramid into two tetrahedra and the prism into three tetrahedra. These tetrahedra are obtained by cutting the quadrilateral face of pyramid and prism with the diagonal which have the vertex with the maximum index (see [HECHT1992] for the reason of this choice).

The triangles on the middle surface obtained with the decomposition of the volume prismatic elements are the triangles generated by the edges on the border of the two dimensional mesh. The label of triangles on the border elements and tetrahedra are defined with the label of these associated elements.

The arguments of buildlayers is a two dimensional mesh and the number of layers $$M$$.

The parameters of this command are:

• zbound= $$[zmin,zmax]$$ where $$zmin$$ and $$zmax$$ are functions expression.

Theses functions define the lower surface mesh and upper mesh of surface mesh.

• coef= A function expression between [0,1].

This parameter is used to introduce degenerate element in mesh.

The number of associated points or vertex $$V_{i}^{2d}$$ is the integer part of $$coef(V_{i}^{2d}) M$$.

• region= This vector is used to initialize the region of tetrahedra.

This vector contains successive pairs of the 2d region number at index $$2i$$ and the corresponding 3d region number at index $$2i+1$$, like change.

• labelmid= This vector is used to initialize the 3d labels number of the vertical face or mid face from the 2d label number.

This vector contains successive pairs of the 2d label number at index $$2i$$ and the corresponding 3d label number at index $$2i+1$$, like change.

• labelup= This vector is used to initialize the 3d label numbers of the upper/top face from the 2d region number.

This vector contains successive pairs of the 2d region number at index $$2i$$ and the corresponding 3d label number at index $$2i+1$$, like change.

• labeldown= Same as the previous case but for the lower/down face label.

Moreover, we also add post processing parameters that allow to moving the mesh. These parameters correspond to parameters transfo, facemerge and ptmerge of the command line movemesh.

The vector region, labelmid, labelup and labeldown These vectors are composed of $$n_{l}$$ successive pairs of number $$O_i,N_l$$ where $$n_{l}$$ is the number (label or region) that we want to get.

An example of this command is given in the Build layer mesh example.

Tip

Cube

 1 //Cube.idp
4
5 func mesh3 Cube (int[int] &NN, real[int, int] &BB, int[int, int] &L){
6     real x0 = BB(0,0), x1 = BB(0,1);
7     real y0 = BB(1,0), y1 = BB(1,1);
8     real z0 = BB(2,0), z1 = BB(2,1);
9
10     int nx = NN[0], ny = NN[1], nz = NN[2];
11
12     // 2D mesh
13     mesh Thx = square(nx, ny, [x0+(x1-x0)*x, y0+(y1-y0)*y]);
14
15     // 3D mesh
16     int[int] rup = [0, L(2,1)], rdown=[0, L(2,0)];
17     int[int] rmid=[1, L(1,0), 2, L(0,1), 3, L(1,1), 4, L(0,0)];
18     mesh3 Th = buildlayers(Thx, nz, zbound=[z0,z1],
19     labelmid=rmid, labelup = rup, labeldown = rdown);
20
21     return Th;
22 }


Tip

Unit cube

1 include "Cube.idp"
2
3 int[int] NN = [10,10,10]; //the number of step in each direction
4 real [int, int] BB = [[0,1],[0,1],[0,1]]; //the bounding box
5 int [int, int] L = [[1,2],[3,4],[5,6]]; //the label of the 6 face left,right, front, back, down, right
6 mesh3 Th = Cube(NN, BB, L);
7 medit("Th", Th);


Tip

Cone

An axisymtric mesh on a triangle with degenerateness

 1 load "msh3"
3
4 // Parameters
5 real RR = 1;
6 real HH = 1;
7
8 int nn=10;
9
10 // 2D mesh
11 border Taxe(t=0, HH){x=t; y=0; label=0;}
12 border Hypo(t=1, 0){x=HH*t; y=RR*t; label=1;}
13 border Vert(t=0, RR){x=HH; y=t; label=2;}
14 mesh Th2 = buildmesh(Taxe(HH*nn) + Hypo(sqrt(HH*HH+RR*RR)*nn) + Vert(RR*nn));
15 plot(Th2, wait=true);
16
17 // 3D mesh
18 real h = 1./nn;
19 int MaxLayersT = (int(2*pi*RR/h)/4)*4;//number of layers
20 real zminT = 0;
21 real zmaxT = 2*pi; //height 2*pi
22 func fx = y*cos(z);
23 func fy = y*sin(z);
24 func fz = x;
25 int[int] r1T = [0,0], r2T = [0,0,2,2], r4T = [0,2];
26 //trick function:
27 //The function defined the proportion
28 //of number layer close to axis with reference MaxLayersT
29 func deg = max(.01, y/max(x/HH, 0.4)/RR);
30 mesh3 Th3T = buildlayers(Th2, coef=deg, MaxLayersT,
31     zbound=[zminT, zmaxT], transfo=[fx, fy, fz],
32     facemerge=0, region=r1T, labelmid=r2T);
33 medit("cone", Th3T);


Tip

Buildlayer mesh

 1 load "msh3"
4
5 // Parameters
6 int C1 = 99;
7 int C2 = 98;
8
9 // 2D mesh
10 border C01(t=0, pi){x=t; y=0; label=1;}
11 border C02(t=0, 2*pi){ x=pi; y=t; label=1;}
12 border C03(t=0, pi){ x=pi-t; y=2*pi; label=1;}
13 border C04(t=0, 2*pi){ x=0; y=2*pi-t; label=1;}
14
15 border C11(t=0, 0.7){x=0.5+t; y=2.5; label=C1;}
16 border C12(t=0, 2){x=1.2; y=2.5+t; label=C1;}
17 border C13(t=0, 0.7){x=1.2-t; y=4.5; label=C1;}
18 border C14(t=0, 2){x=0.5; y=4.5-t; label=C1;}
19
20 border C21(t=0, 0.7){x=2.3+t; y=2.5; label=C2;}
21 border C22(t=0, 2){x=3; y=2.5+t; label=C2;}
22 border C23(t=0, 0.7){x=3-t; y=4.5; label=C2;}
23 border C24(t=0, 2){x=2.3; y=4.5-t; label=C2;}
24
25 mesh Th = buildmesh(C01(10) + C02(10) + C03(10) + C04(10)
26     + C11(5) + C12(5) + C13(5) + C14(5)
27     + C21(-5) + C22(-5) + C23(-5) + C24(-5));
28
29 mesh Ths = buildmesh(C01(10) + C02(10) + C03(10) + C04(10)
30     + C11(5) + C12(5) + C13(5) + C14(5));
31
32 // Construction of a box with one hole and two regions
33 func zmin = 0.;
34 func zmax = 1.;
35 int MaxLayer = 10;
36
37 func XX = x*cos(y);
38 func YY = x*sin(y);
39 func ZZ = z;
40
41 int[int] r1 = [0, 41], r2 = [98, 98, 99, 99, 1, 56];
42 int[int] r3 = [4, 12];  //the triangles of uppper surface mesh
43                         //generated by the triangle in the 2D region
44                         //of mesh Th of label 4 as label 12
45 int[int] r4 = [4, 45];  //the triangles of lower surface mesh
46                         //generated by the triangle in the 2D region
47                         //of mesh Th of label 4 as label 45.
48
49 mesh3 Th3 = buildlayers(Th, MaxLayer, zbound=[zmin, zmax], region=r1,
50     labelmid=r2, labelup=r3, labeldown=r4);
51     medit("box 2 regions 1 hole", Th3);
52
53 // Construction of a sphere with TetGen
54 func XX1 = cos(y)*sin(x);
55 func YY1 = sin(y)*sin(x);
56 func ZZ1 = cos(x);
57
58 real[int] domain = [0., 0., 0., 0, 0.001];
59 string test = "paACQ";
60 cout << "test = " << test << endl;
61 mesh3 Th3sph = tetgtransfo(Ths, transfo=[XX1, YY1, ZZ1],
62     switch=test, nbofregions=1, regionlist=domain);
63 medit("sphere 2 regions", Th3sph);


### Remeshing

Note

if an operation on a mesh3 is performed then the same operation is applyed on its surface part (its meshS associated)

#### The command change

This command changes the label of elements and border elements of a mesh. It’s the equivalent command in 2d mesh case.

Changing the label of elements and border elements will be done using the keyword change. The parameters for this command line are for two dimensional and three dimensional cases:

• reftet= is a vector of integer that contains successive pairs of the old label number to the new label number.

• refface= is a vector of integer that contains successive pairs of the old region number to new region number.

• flabel= is an integer function given the new value of the label.

• fregion= is an integer function given the new value of the region.

• rmInternalFaces= is a boolean, equal true to remove the internal faces.

• rmlfaces= is a vector of integer, where triangle’s label given are remove of the mesh

These vectors are composed of $$n_{l}$$ successive pairs of numbers $$O,N$$ where $$n_{l}$$ is the number (label or region) that we want to change. For example, we have:

$\begin{split}\mathtt{label} &= [ O_{1}, N_{1}, ..., O_{n_{l}},N_{n_{l}} ] \\ \mathtt{region} &= [ O_{1}, N_{1}, ..., O_{n_{l}},N_{n_{l}} ]\end{split}$

An example of use:

 1 // Mesh
2 mesh3 Th1 = cube(10, 10);
3 mesh3 Th2 = cube(20, 10, [x+1, y,z]);
4
5 int[int] r1=[2,0];
6 plot(Th1, wait=true);
7
8 Th1 = change(Th1, label=r1); //change the label of Edges 2 in 0.
9 plot(Th1, wait=true);
10
11 // boundary label: 1 -> 1 bottom, 2 -> 1 right, 3->1 top, 4->1 left boundary label is 1
12 int[int] re=[1,1, 2,1, 3,1, 4,1]
13 Th2=change(Th2,refe=re);
14 plot(Th2,wait=1) ;


#### The command trunc

This operator have been introduce to remove a piece of mesh or/and split all element or for a particular label element The three named parameter - boolean function to keep or remove elements - split= sets the level n of triangle splitting. each triangle is splitted in n × n ( one by default) - freefem:label= sets the label number of new boundary item (1 by default)

An example of use

1 load "msh3"
3 int nn=8;
4 mesh3 Th=cube(nn,nn,nn);
5 //  remove the small cube $]1/2,1[^2$
6 Th= trunc(Th,((x<0.5) |(y< 0.5)| (z<0.5)), split=3, label=3);
7 medit("cube",Th);


#### The command movemesh

3D meshes can be translated, rotated, and deformed using the command line movemesh as in the 2D case (see section movemesh). If $$\Omega$$ is tetrahedrized as $$T_{h}(\Omega)$$, and $$\Phi(x,y)=(\Phi1(x,y,z), \Phi2(x,y,z), \Phi3(x,y,z))$$ is the transformation vector then $$\Phi(T_{h})$$ is obtained by:

1 mesh3 Th = movemesh(Th, [Phi1, Phi2, Phi3], ...);
2 mesh3 Th = movemesh3(Th, transfo=[Phi1, Phi2, Phi3], ...);  (syntax with transfo=)


The parameters of movemesh in three dimensions are:

• transfo= sets the geometric transformation $$\Phi(x,y)=(\Phi1(x,y,z), \Phi2(x,y,z), \Phi3(x,y,z))$$

• region= sets the integer labels of the tetrahedra.

0 by default.

• label= sets the labels of the border faces.

This parameter is initialized as the label for the keyword change.

• facemerge= An integer expression.

When you transform a mesh, some faces can be merged. This parameter equals to one if the merges’ faces is considered. Otherwise it equals to zero. By default, this parameter is equal to 1.

• ptmerge = A real expression.

When you transform a mesh, some points can be merged. This parameter is the criteria to define two merging points. By default, we use

$ptmerge \: = \: 1e-7 \: \:Vol( B ),$

where $$B$$ is the smallest axis parallel boxes containing the discretion domain of $$\Omega$$ and $$Vol(B)$$ is the volume of this box.

• orientation = An integer expression equal 1, give the oientation of the triangulation, elements must be in the reference orientation (counter clock wise) equal -1 reverse the orientation of the tetrahedra

Note

The orientation of tetrahedra are checked by the positivity of its area and automatically corrected during the building of the adjacency.

An example of this command can be found in the Poisson’s equation 3D example.

 1 load "medit"
2 include "cube.idp"
3 int[int]  Nxyz=[20,5,5];
4 real [int,int]  Bxyz=[[0.,5.],[0.,1.],[0.,1.]];
5 int [int,int]  Lxyz=[[1,2],[2,2],[2,2]];
6 real E = 21.5e4;
7 real sigma = 0.29;
8 real mu = E/(2*(1+sigma));
9 real lambda = E*sigma/((1+sigma)*(1-2*sigma));
10 real gravity = -0.05;
11 real sqrt2=sqrt(2.);
12
13 mesh3 Th=Cube(Nxyz,Bxyz,Lxyz);
14 fespace Vh(Th,[P1,P1,P1]);
15 Vh [u1,u2,u3], [v1,v2,v3];
16
17 macro epsilon(u1,u2,u3)  [dx(u1),dy(u2),dz(u3),(dz(u2)+dy(u3))/sqrt2,(dz(u1)+dx(u3))/sqrt2,(dy(u1)+dx(u2))/sqrt2] // EOM
18 macro div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM
19
20 solve Lame([u1,u2,u3],[v1,v2,v3])=
21   int3d(Th)(
22             lambda*div(u1,u2,u3)*div(v1,v2,v3)
23             +2.*mu*( epsilon(u1,u2,u3)'*epsilon(v1,v2,v3) )
24               )
25   - int3d(Th) (gravity*v3)
26   + on(1,u1=0,u2=0,u3=0);
27
28 real dmax= u1[].max;
29 real coef= 0.1/dmax;
30
31 int[int] ref2=[1,0,2,0]; // array
32 mesh3 Thm=movemesh(Th,[x+u1*coef,y+u2*coef,z+u3*coef],label=ref2);
33 // mesh3 Thm=movemesh3(Th,transfo=[x+u1*coef,y+u2*coef,z+u3*coef],label=ref2); older syntax
34 Thm=change(Thm,label=ref2);
35 plot(Th,Thm, wait=1,cmm="coef  amplification = "+coef );


movemesh doesn’t use the prefix tranfo= [.,.,.], the geometric transformation is directly given by [.,.,.] in the arguments list

#### The command extract

This command offers the possibility to extract a boundary part of a mesh3

• refface , is a vector of integer that contains a list of triangle face references, where the extract function must be apply.

• label , is a vector of integer that contains a list of tetrahedra label

1 load"msh3"
2 int nn = 30;
3 int[int] labs = [1, 2, 2, 1, 1, 2]; // Label numbering
4 mesh3 Th = cube(nn, nn, nn, label=labs);
5 // extract the surface (boundary) of the cube
6 int[int] llabs = [1, 2];
7 meshS ThS = extract(Th,label=llabs);


#### The command buildSurface

This new function allows to build the surface mesh of a volume mesh, under the condition the surface is the boundary of the volume. By definition, a mesh3 is defined by a list of vertices, tetrahedron elements and triangle border elements. buildSurface function create the meshS corresponding, given the list vertices which are on the border domain, the triangle elements and build the list of edges. Remark, for a closed surface mesh, the edges list is empty.

#### The command movemesh23

A simple method to tranform a 2D mesh in 3D Surface mesh. The principe is to project a two dimensional domain in a three dimensional space, 2d surface in the (x,y,z)-space to create a surface mesh 3D, meshS.

Warning

Since the release 4.2.1, the FreeFEM function movemesh23 returns a meshS type.

This corresponds to translate, rotate or deforme the domain by a displacement vector of this form $$\mathbf{\Phi(x,y)} = (\Phi1(x,y), \Phi2(x,y), \Phi3(x,y))$$.

The result of moving a two dimensional mesh Th2 by this three dimensional displacement is obtained using:

1 **meshS** Th3 = movemesh23(Th2, transfo=[Phi(1), Phi(2), Phi(3)]);


The parameters of this command line are:

• transfo= [$$\Phi 1$$, $$\Phi 2$$, $$\Phi 3$$] sets the displacement vector of transformation $$\mathbf{\Phi(x,y)} = [\Phi1(x,y), \Phi2(x,y), \Phi3(x,y)]$$.

• label= sets an integer label of triangles.

• orientation= sets an integer orientation to give the global orientation of the surface of mesh. Equal 1, give a triangulation in the reference orientation (counter clock wise) equal -1 reverse the orientation of the triangles

• ptmerge= A real expression.

When you transform a mesh, some points can be merged. This parameter is the criteria to define two merging points. By default, we use

$ptmerge \: = \: 1e-7 \: \:Vol( B ),$

where $$B$$ is the smallest axis, parallel boxes containing the discretized domain of $$\Omega$$ and $$Vol(B)$$ is the volume of this box.

We can do a “gluing” of surface meshes using the process given in Change section. An example to obtain a three dimensional mesh using the command line tetg and movemesh23 is given below.

 1 load "msh3"
3
4 // Parameters
5 real x10 = 1.;
6 real x11 = 2.;
7 real y10 = 0.;
8 real y11 = 2.*pi;
9
10 func ZZ1min = 0;
11 func ZZ1max = 1.5;
12 func XX1 = x;
13 func YY1 = y;
14
15 real x20 = 1.;
16 real x21 = 2.;
17 real y20=0.;
18 real y21=1.5;
19
20 func ZZ2 = y;
21 func XX2 = x;
22 func YY2min = 0.;
23 func YY2max = 2*pi;
24
25 real x30=0.;
26 real x31=2*pi;
27 real y30=0.;
28 real y31=1.5;
29
30 func XX3min = 1.;
31 func XX3max = 2.;
32 func YY3 = x;
33 func ZZ3 = y;
34
35 // Mesh
36 mesh Thsq1 = square(5, 35, [x10+(x11-x10)*x, y10+(y11-y10)*y]);
37 mesh Thsq2 = square(5, 8, [x20+(x21-x20)*x, y20+(y21-y20)*y]);
38 mesh Thsq3 = square(35, 8, [x30+(x31-x30)*x, y30+(y31-y30)*y]);
39
40 // Mesh 2D to 3D surface
41 meshS Th31h = movemesh23(Thsq1, transfo=[XX1, YY1, ZZ1max], orientation=1);
42 meshS Th31b = movemesh23(Thsq1, transfo=[XX1, YY1, ZZ1min], orientation=-1);
43
44 meshS Th32h = movemesh23(Thsq2, transfo=[XX2, YY2max, ZZ2], orientation=-1);
45 meshS Th32b = movemesh23(Thsq2, transfo=[XX2, YY2min, ZZ2], orientation=1);
46
47 meshS Th33h = movemesh23(Thsq3, transfo=[XX3max, YY3, ZZ3], orientation=1);
48 meshS Th33b = movemesh23(Thsq3, transfo=[XX3min, YY3, ZZ3], orientation=-1);
49
50 // Gluing surfaces
51 meshS Th33 = Th31h + Th31b + Th32h + Th32b + Th33h + Th33b;
52 plot(Th33, cmm="Th33");
53
54 // Tetrahelize the interior of the cube with TetGen
55 real[int] domain =[1.5, pi, 0.75, 145, 0.0025];
56 meshS Thfinal = tetg(Th33, switch="paAAQY", regionlist=domain);
57 plot(Thfinal, cmm="Thfinal");
58
59 // Build a mesh of a half cylindrical shell of interior radius 1, and exterior radius 2 and a height of 1.5
60 func mv2x = x*cos(y);
61 func mv2y = x*sin(y);
62 func mv2z = z;
63 meshS Thmv2 = movemesh(Thfinal, transfo=[mv2x, mv2y, mv2z], facemerge=0);
64 plot(Thmv2, cmm="Thmv2");


### 3d Meshing examples

Tip

Lake

 1 load "msh3"
3
4 // Parameters
5 int nn = 5;
6
7 // 2D mesh
8 border cc(t=0, 2*pi){x=cos(t); y=sin(t); label=1;}
9 mesh Th2 = buildmesh(cc(100));
10
11 // 3D mesh
12 int[int] rup = [0, 2], rlow = [0, 1];
13 int[int] rmid = [1, 1, 2, 1, 3, 1, 4, 1];
14 func zmin = 2-sqrt(4-(x*x+y*y));
15 func zmax = 2-sqrt(3.);
16
17 mesh3 Th = buildlayers(Th2, nn,
18     coef=max((zmax-zmin)/zmax, 1./nn),
19     zbound=[zmin,zmax],
20     labelmid=rmid,
21     labelup=rup,
22     labeldown=rlow);
23
24 medit("Th", Th);


Tip

Hole region

 1 load "msh3"
4
5 // 2D mesh
6 mesh Th = square(10, 20, [x*pi-pi/2, 2*y*pi]); // ]-pi/2, pi/2[X]0,2pi[
7
8 // 3D mesh
9 //parametrization of a sphere
10 func f1 = cos(x)*cos(y);
11 func f2 = cos(x)*sin(y);
12 func f3 = sin(x);
13 //partial derivative of the parametrization
14 func f1x = sin(x)*cos(y);
15 func f1y = -cos(x)*sin(y);
16 func f2x = -sin(x)*sin(y);
17 func f2y = cos(x)*cos(y);
18 func f3x = cos(x);
19 func f3y = 0;
20 //M = DF^t DF
21 func m11 = f1x^2 + f2x^2 + f3x^2;
22 func m21 = f1x*f1y + f2x*f2y + f3x*f3y;
23 func m22 = f1y^2 + f2y^2 + f3y^2;
24
25 func perio = [[4, y], [2, y], [1, x], [3, x]];
26 real hh = 0.1;
27 real vv = 1/square(hh);
28 verbosity = 2;
29 Th = adaptmesh(Th, m11*vv, m21*vv, m22*vv, IsMetric=1, periodic=perio);
30 Th = adaptmesh(Th, m11*vv, m21*vv, m22*vv, IsMetric=1, periodic=perio);
31 plot(Th, wait=true);
32
33 //construction of the surface of spheres
34 real Rmin = 1.;
35 func f1min = Rmin*f1;
36 func f2min = Rmin*f2;
37 func f3min = Rmin*f3;
38
39 meshS ThSsph = movemesh23(Th, transfo=[f1min, f2min, f3min]);
40
41 real Rmax = 2.;
42 func f1max = Rmax*f1;
43 func f2max = Rmax*f2;
44 func f3max = Rmax*f3;
45
46 meshS ThSsph2 = movemesh23(Th, transfo=[f1max, f2max, f3max]);
47
48 //gluing meshes
49 meshS ThS = ThSsph + ThSsph2;
50
51 cout << " TetGen call without hole " << endl;
52 real[int] domain2 = [1.5, 0., 0., 145, 0.001, 0.5, 0., 0., 18, 0.001];
53 mesh3 Th3fin = tetg(ThS, switch="paAAQYY", nbofregions=2, regionlist=domain2);
54 medit("Sphere with two regions", Th3fin);
55
56 cout << " TetGen call with hole " << endl;
57 real[int] hole = [0.,0.,0.];
58 real[int] domain = [1.5, 0., 0., 53, 0.001];
59 mesh3 Th3finhole = tetg(ThS, switch="paAAQYY",
60     nbofholes=1, holelist=hole, nbofregions=1, regionlist=domain);
61 medit("Sphere with a hole", Th3finhole);


Tip

Build a 3d mesh of a cube with a balloon

 1 load "msh3"
4 include "MeshSurface.idp"
5
6 // Parameters
7 real hs = 0.1; //mesh size on sphere
8 int[int] N = [20, 20, 20];
9 real [int,int] B = [[-1, 1], [-1, 1], [-1, 1]];
10 int [int,int] L = [[1, 2], [3, 4], [5, 6]];
11
12 // Meshes
13 meshS ThH = SurfaceHex(N, B, L, 1);
14 meshS ThS = Sphere(0.5, hs, 7, 1);
15
16 meshS ThHS = ThH + ThS;
17 medit("Hex-Sphere", ThHS);
18
19 real voltet = (hs^3)/6.;
20 cout << "voltet = " << voltet << endl;
21 real[int] domain = [0, 0, 0, 1, voltet, 0, 0, 0.7, 2, voltet];
22 mesh3 Th = tetg(ThHS, switch="pqaAAYYQ", nbofregions=2, regionlist=domain);
23 medit("Cube with ball", Th);


Cube sphere

## The type meshS in 3 dimension

Warning

Since the release 4.2.1, the surface mesh3 object (list of vertices and border elements, without tetahedra elements) is remplaced by meshS type.

### Commands for 3d surface mesh generation

#### The command square3

The function square3 like the function square in 2d is the simple way to a build the unit square plan in the space $$\mathbb{R^3}$$. To use this command, it is necessary to load the pluging msh3 (need load "msh3"). A square in 3d consists in building a 2d square which is projected from $$\mathbb{R^2}$$ to $$\mathbb{R^3}$$. The parameters of this command line are:

• n,m generates a n×m grid in the unit square

• [.,.,.] is [ $$\Phi 1$$, $$\Phi 2$$, $$\Phi 3$$ ] is the geometric transformation from $$\mathbb{R^2}$$ to $$\mathbb{R^3}$$. By default, [ $$\Phi 1$$, $$\Phi 2$$, $$\Phi 3$$ ] = [x,y,0]

• orientation= equal 1, gives the orientation of the triangulation, elements are in the reference orientation (counter clock wise) equal -1 reverse the orientation of the triangles it’s the global orientation of the surface 1 extern (-1 intern)

 1      real R = 3, r=1;
2      real h = 0.2; //
3      int nx = R*2*pi/h;
4      int ny = r*2*pi/h;
5      func torex= (R+r*cos(y*pi*2))*cos(x*pi*2);
6      func torey= (R+r*cos(y*pi*2))*sin(x*pi*2);
7      func torez= r*sin(y*pi*2);
8
9
10      meshS ThS=square3(nx,ny,[torex,torey,torez],orientation=-1) ;


The following code generates a $$3\times 4 \times 5$$ grid in the unit cube $$[0, 1]^3$$ with a clock wise triangulation.

#### surface mesh builders

Adding at the top of a FreeFEM script include "MeshSurface.idp", constructors of sphere, ellipsoid, surface mesh of a 3d box are available.

• SurfaceHex(N, B, L, orient)

• this operator allows to build the surface mesh of a 3d box

• int[int] N=[nx,ny,nz]; // the number of seg in the 3 direction

• real [int,int] B=[[xmin,xmax],[ymin,ymax],[zmin,zmax]]; // bounding bax

• int [int,int] L=[[1,2],[3,4],[5,6]]; // the label of the 6 face left,right, front, back, down, right

• orient the global orientation of the surface 1 extern (-1 intern),

• returns a meshS type

• Ellipsoide (RX, RY, RZ, h, L, OX, OY, OZ, orient)

• h is the mesh size

• L is the label

• orient the global orientation of the surface 1 extern (-1 intern)

• OX, OY, OZ are real numbers to give the Ellipsoide center ( optinal, by default is (0,0,0) )

• where RX, RY, RZ are real numbers such as the parametric equations of the ellipsoid is:

• returns a meshS type

$\forall u \in [- \frac{\pi}{2},\frac{\pi}{2} [ \text{ and } v \in [0, 2 \pi], \vectthree{x=\text{Rx } cos(u)cos(v) + \text{Ox }}{y=\text{Ry } cos(u)sin(v) + \text{Oy }}{z = \text{Rz } sin(v) + \text{Oz } }$
• Sphere(R, h, L, OX, OY, OZ, orient)

• where R is the raduis of the sphere,

• OX, OY, OZ are real numbers to give the Ellipsoide center ( optinal, by default is (0,0,0) )

• h is the mesh size of the shpere

• L is the label the the sphere

• orient the global orientation of the surface 1 extern (-1 intern)

• returns a meshS type

 1 func meshS SurfaceHex(int[int] & N,real[int,int] &B ,int[int,int] & L,int orientation){
2     real x0=B(0,0),x1=B(0,1);
3     real y0=B(1,0),y1=B(1,1);
4     real z0=B(2,0),z1=B(2,1);
5
6     int nx=N[0],ny=N[1],nz=N[2];
7
8     mesh Thx = square(ny,nz,[y0+(y1-y0)*x,z0+(z1-z0)*y]);
9     mesh Thy = square(nx,nz,[x0+(x1-x0)*x,z0+(z1-z0)*y]);
10     mesh Thz = square(nx,ny,[x0+(x1-x0)*x,y0+(y1-y0)*y]);
11
12     int[int] refx=[0,L(0,0)],refX=[0,L(0,1)];   //  Xmin, Ymax faces labels renumbering
13     int[int] refy=[0,L(1,0)],refY=[0,L(1,1)];   //  Ymin, Ymax faces labesl renumbering
14     int[int] refz=[0,L(2,0)],refZ=[0,L(2,1)];   //  Zmin, Zmax faces labels renumbering
15
16     meshS Thx0 = movemesh23(Thx,transfo=[x0,x,y],orientation=-orientation,label=refx);
17     meshS Thx1 = movemesh23(Thx,transfo=[x1,x,y],orientation=+orientation,label=refX);
18     meshS Thy0 = movemesh23(Thy,transfo=[x,y0,y],orientation=+orientation,label=refy);
19     meshS Thy1 = movemesh23(Thy,transfo=[x,y1,y],orientation=-orientation,label=refY);
20     meshS Thz0 = movemesh23(Thz,transfo=[x,y,z0],orientation=-orientation,label=refz);
21     meshS Thz1 = movemesh23(Thz,transfo=[x,y,z1],orientation=+orientation,label=refZ);
22     meshS Th= Thx0+Thx1+Thy0+Thy1+Thz0+Thz1;
23
24  return Th;
25  }
26
27  func meshS Ellipsoide (real RX,real RY,real RZ,real h,int L,real Ox,real Oy,real Oz,int orientation) {
28     mesh  Th=square(10,20,[x*pi-pi/2,2*y*pi]);  //  $]\frac{-pi}{2},frac{-pi}{2}[\times]0,2\pi[$
29     //  a parametrization of a sphere
30     func f1 =RX*cos(x)*cos(y);
31     func f2 =RY*cos(x)*sin(y);
32     func f3 =RZ*sin(x);
33         //    partiel derivative
34     func f1x= -RX*sin(x)*cos(y);
35     func f1y= -RX*cos(x)*sin(y);
36     func f2x= -RY*sin(x)*sin(y);
37     func f2y= +RY*cos(x)*cos(y);
38     func f3x=-RZ*cos(x);
39     func f3y=0;
40     // the metric on the sphere  $M = DF^t DF$
41     func m11=f1x^2+f2x^2+f3x^2;
42     func m21=f1x*f1y+f2x*f2y+f3x*f3y;
43     func m22=f1y^2+f2y^2+f3y^2;
44     func perio=[[4,y],[2,y],[1,x],[3,x]];  // to store the periodic condition
45     real hh=h;// hh  mesh size on unite sphere
46     real vv= 1/square(hh);
51     int[int] ref=[0,L];
52     meshS ThS=movemesh23(Th,transfo=[f1,f2,f3],orientation=orientation,refface=ref);
54   return ThS;
55  }
56
57  func meshS Ellipsoide (real RX,real RY,real RZ,real h,int L,int orientation) {
58   return Ellipsoide (RX,RY,RZ,h,L,0.,0.,0.,orientation);
59  }
60  func meshS Sphere(real R,real h,int L,int orientation) {
61   return Ellipsoide(R,R,R,h,L,orientation);
62  }
63  func meshS Sphere(real R,real h,int L,real Ox,real Oy,real Oz,int orientation) {
64   return Ellipsoide(R,R,R,h,L,Ox,Oy,Oz,orientation);
65  }


#### 2D mesh generators combined with movemesh23

FreeFEM ‘s meshes can be built by the composition of the movemesh23 command from a 2d mesh generation. The operation is a projection of a 2d plane in $$\mathbb{R^3}$$ following the geometric transformation [ $$\Phi 1$$, $$\Phi 2$$, $$\Phi 3$$ ].

 1 load "msh3"
2 real l = 3;
3 border a(t=-l,l){x=t; y=-l;label=1;};
4 border b(t=-l,l){x=l; y=t;label=1;};
5 border c(t=l,-l){x=t; y=l;label=1;};
6 border d(t=l,-l){x=-l; y=t;label=1;};
7 int n = 100;
8 border i(t=0,2*pi){x=1.1*cos(t);y=1.1*sin(t);label=5;};
9 mesh th= buildmesh(a(n)+b(n)+c(n)+d(n)+i(-n));
10 meshS Th= movemesh23(th,transfo=[x,y,cos(x)^2+sin(y)^2]);


### Remeshing

#### The command trunc

This operator allows to define a meshS by truncating another one, i.e. by removing triangles, and/or by splitting each triangle by a given positive integer s. In a FreeFEM script, this function must be called as follows:

meshS TS2= trunc (TS1, boolean function to keep or remove elements, split = s, label = …)

The command has the following arguments:

• boolean function to keep or remove elements

• split= sets the level n of triangle splitting. each triangle is splitted in n × n ( one by default)

• label= sets the label number of new boundary item (1 by default)

An example of how to call the function

 1 real R = 3, r=1;
2 real h = 0.2; //
3 int nx = R*2*pi/h;
4 int ny = r*2*pi/h;
5 func torex= (R+r*cos(y*pi*2))*cos(x*pi*2);
6 func torey= (R+r*cos(y*pi*2))*sin(x*pi*2);
7 func torez= r*sin(y*pi*2);
8 // build a tore
9 meshS ThS=square3(nx,ny,[torex,torey,torez]) ;
10 ThS=trunc(ThS, (x < 0.5) | (y < 0.5) | (z > 1.), split=4);


#### The command movemesh

Like 2d and 3d type meshes in FreeFEM, meshS can be translated, rotated or deformated by an application [$$\Phi 1$$, $$\Phi 2$$, $$\Phi 3$$]. The image $$T_{h}(\Omega)$$ is obtained by the command movemeshS.

The parameters of movemeshS are:

• transfo= sets the geometric transformation $$\Phi(x,y)=(\Phi1(x,y,z), \Phi2(x,y,z), \Phi3(x,y,z))$$

• region= sets the integer labels of the triangles.

0 by default.

• label= sets the labels of the border edges.

This parameter is initialized as the label for the keyword change.

• edgemerge= An integer expression.

When you transform a mesh, some triangles can be merged and fix the parameter to 1, else 0 By default, this parameter is equal to 1.

• ptmerge = A real expression.

When you transform a mesh, some points can be merged. This parameter is the criteria to define two merging points. By default, we use

$ptmerge \: = \: 1e-7 \: \:Vol( B ),$

where $$B$$ is the smallest axis parallel boxes containing the discretion domain of $$\Omega$$ and $$Vol(B)$$ is the volume of this box.

• orientation = An integer expression

equal 1, give the oientation of the triangulation, elements must be in the reference orientation (counter clock wise) equal -1 reverse the orientation of the triangles. It’s the global orientation of the normals at the surface 1 extern (-1 intern)

Example of using

1 meshS Th1 = square3(n,n,[2*x,y,1],orientation=-1);
2 meshS Th2=movemeshS(Th1, transfo=[x,y,z]);
3 meshS Th3=movemesh(Th1, [x,y,z]);


#### The command change

Equivalent for a 2d or 3d mesh, the command change changes the label of elements and border elements of a meshS.

The parameters for this command line are:

• reftri= is a vector of integer that contains successive pairs of the old label number to the new label number for elements.

• refedge= is a vector of integer that contains successive pairs of the old region number to new region number for boundary elements.

• flabel= is an integer function given the new value of the label.

• fregion= is an integer function given the new value of the region.

• rmledges= is a vector of integer, where edge’s label given are remove of the mesh

These vectors are composed of $$n_{l}$$ successive pairs of numbers $$O,N$$ where $$n_{l}$$ is the number (label or region) that we want to change. For example, we have:

$\begin{split}\mathtt{label} &= [ O_{1}, N_{1}, ..., O_{n_{l}},N_{n_{l}} ] \\ \mathtt{region} &= [ O_{1}, N_{1}, ..., O_{n_{l}},N_{n_{l}} ]\end{split}$

## The type meshL in 3 dimension

### Commands for 3d curve mesh generation

#### The command segment

The function segment is a basic command to define a curve in 3D space.

The parameters of this command line are:

• n generates a n subsegments from the unit line

• [.,.,.] is [ $$\Phi 1$$, $$\Phi 2$$, $$\Phi 3$$ ] is the geometric transformation from $$\mathbb{R^1}$$ to $$\mathbb{R^3}$$. By default, [ $$\Phi 1$$, $$\Phi 2$$, $$\Phi 3$$ ] = [x,0,0]

• orientation= equal 1, gives the orientation of the triangulation, elements are in the reference orientation (counter clock wise) equal -1 reverse the orientation of the triangles it’s the global orientation of the surface 1 extern (-1 intern)

• cleanmesh= is a boolean, allowing remove the duplicated nodes

• removeduplicate= is a boolean, allowing remove the duplicated elements and border elements

• precismesh this parameter is the criteria to define two merging points.

By default, it value is 1e-7 and define the smallest axis parallel boxes containing the discretion domain of $$\Omega$$

By defaut, the border points are marked by label 1 and 2.

1      real R = 3, r=1;
2      real h = 0.1; //
3      int nx = R*2*pi/h;
4      func torex= (R+r*cos(y*pi*2))*cos(x*pi*2);
5      func torey= (R+r*cos(y*pi*2))*sin(x*pi*2);
6      func torez= r*sin(y*pi*2);
7      meshL Th=segment(nx,[torex,torey,torez],removeduplicate=true) ;


The following code generates a 10 subsegments from the unit line with a clock wise triangulation, according to the geometric transformation [torex,torey,torez] and removing the duplicated points/elements

#### The command buildmesh

This operator allows to define a curve mesh from multi-borders. The domain can be defined by a parametrized curve (keyword border), such as Th1 in the following example or piecewise by parametrized curves, such as the construction of the mesh Th2.

The pieces can only intersect at their endpoints, but it is possible to join more than two endpoints.

 1 load "msh3"
2
3 // conical helix
4 border E1(t=0, 10.*pi){x=(1.)*t*cos(t); y=-(1.)*t*sin(t); z=t;}
5 meshL Th1=buildmeshL(E1(1000));
6
7 int upper = 1, others = 2, inner = 3, n = 10;
8 border D01(t=0, 1) {x=0; y=-1+t; }
9 border D02(t=0, 1){x=1.5-1.5*t; y=-1; z=3;label=upper;}
10 border D03(t=0, 1){x=1.5; y=-t; z=3;label=upper;}
11 border D04(t=0, 1){x=1+0.5*t; y=0; z=3;label=others;}
12 border D05(t=0, 1){x=0.5+0.5*t; y=0; z=3;label=others;}
13 border D06(t=0, 1){x=0.5*t; y=0; z=3;label=others;}
14 border D11(t=0, 1){x=0.5; y=-0.5*t; z=3;label=inner;}
15 border D12(t=0, 1){x=0.5+0.5*t; y=-0.5; z=3;label=inner;}
16 border D13(t=0, 1){x=1; y=-0.5+0.5*t; z=3;label=inner;}
17
18 meshL Th2=buildmeshL(D01(-n) + D02(-n) + D03(-n) + D04(-n) + D05(-n)
19    + D06(-n) + D11(n) + D12(n) + D13(n));


### Remeshing

#### The command trunc

This operator allows to define a meshL by truncating another one, i.e. by removing segments, and/or by splitting each element by a given positive integer s. Here, an example to use this function:

meshL ThL2= trunc (ThL1, boolean function to keep or remove elements, split = s, label = …)

The command has the following arguments:

• boolean function to keep or remove elements

• split= sets the level n of edge splitting, each edge is splitted in n subpart( one by default)

• label= sets the label number of new boundary item (1 by default)

• new2old

• old2new

• renum

• orientation=

equal 1, gives the orientation of the triangulation, elements are in the reference orientation (counter clock wise) equal -1 reverse the orientation of the triangles it’s the global orientation of the surface 1 extern (-1 intern)

• cleanmesh= is a boolean, allowing remove the duplicated nodes

• removeduplicate= is a boolean, allowing remove the duplicated elements and border elements

• precismesh this parameter is the criteria to define two merging points.

By default, it value is 1e-7 and define the smallest axis parallel boxes containing the discretion domain of $$\Omega$$

An example of how to call this function

1 int nx=10;
2 meshL Th=segment(nx,[5.*x,cos(pi*x),sin(pi*x)]);
3 Th=trunc(Th, (x < 0.5) | (y < 0.5) | (z > 1.), split=4);


#### The command movemesh

This is the classical mesh transformation FreeFEM function, meshL can be deformed by an application [ $$\Phi 1$$, $$\Phi 2$$, $$\Phi 3$$ ]. The image $$T_{h}(\Omega)$$ is obtained by the command movemeshL.

The parameters of movemesh are:

• transfo= sets the geometric transformation $$\Phi(x,y)=(\Phi1(x,y,z), \Phi2(x,y,z), \Phi3(x,y,z))$$

• refedge= sets the integer labels of the triangles.

0 by default.

• refpoint= sets the labels of the border points.

This parameter is initialized as the label for the keyword change.

• precismesh this parameter is the criteria to define two merging points.

By default, it value is 1e-7 and define the smallest axis parallel boxes containing the discretion domain of $$\Omega$$

• orientation = An integer expression

equal 1, give the oientation of the triangulation, elements must be in the reference orientation (counter clock wise) equal -1 reverse the orientation of the triangles. It’s the global orientation of the normals at the surface 1 extern (-1 intern)

• cleanmesh= is a boolean, allowing remove the duplicated nodes

• removeduplicate= is a boolean, allowing remove the duplicated elements and border elements

Note

The definition of the geometric transformation depends on the space dimension of the studied problem. It means that, with curve FEM, it’s possible to treat a real 1D problem (space coordinate is x) then the transformation is given by x: ->F(x), that means [F_x] and F_y=F_z=0 in FreeFEM function.

Example of using

1 int nx=100;
2 meshL Th=Sline(nx);
3 meshL Th31=movemesh(Th, [x]);
4 meshL Th32=movemesh(Th, [x,-x*(x-1)]);
5 meshL Th3=Th31+Th32;


#### The command change

The command change changes the label of elements and border elements of a meshL.

The parameters for this command line are:

• refedge= is a vector of integer that contains successive pairs of the old label number to the new label number for elements.

• refpoint= is a vector of integer that contains successive pairs of the old region number to new region number for boundary elements.

• flabel= is an integer function given the new value of the label.

• fregion= is an integer function given the new value of the region.

• rmlpoint= is a vector of integer, where edge’s label given are remove of the mesh

These vectors are composed of $$n_{l}$$ successive pairs of numbers $$O,N$$ where $$n_{l}$$ is the number (label or region) that we want to change. For example, we have:

$\begin{split}\mathtt{label} &= [ O_{1}, N_{1}, ..., O_{n_{l}},N_{n_{l}} ] \\ \mathtt{region} &= [ O_{1}, N_{1}, ..., O_{n_{l}},N_{n_{l}} ]\end{split}$

#### The commands buildBdMesh and Gamma

The command Gamma allows to extract the border mesh independly of a surface mesh. With this function, the constructed border mesh contains the full geometric description of th eboundary surface. In case where the border mesh doesn’t exist, before calling Gamma, must build it by calling the buildBdMesh function (see the next function description).

1 load "msh3"
2 int n= 10;
3 meshS Th=square3(n,n);
4 Th = buildBdMesh(Th); // build the border mesh
5 // build Th1, the border of Th, defined by edges elements and point border elements
6 meshL Th1 = Th.Gamma;


#### Glue of meshL meshes

An assembling of meshL is possible thanks to the operator +. The result returns a meshL, with caution of the right orientation at the merged interfaces. Here, the function checkmesh can be called.

1  int n=10;
2  meshL Th1 = segment(n);
3  meshL Th2 = segment(n,[0,x,0],orientation=1);
4  meshL Th3 = segment(n,[x,0,1],orientation=1);
5  meshL Th4 = segment(n,[0,0,x],orientation=-1);
6
7  meshL Th = Th1+Th2+Th3+Th4;
8  Th=rebuildBorder(Th, ridgeangledetection=pi/2.+0.0001);


Warning

For the moment, the case of no manifold mesh are not considered in FreeFEM. To check if the meshL contains no manifold elements, the command nbnomanifold.

#### The command extract

This operator allows to extract a labeled piece or the entire border of a 2D mesh and project it in 3D. Optionally, a geometic transformation can be applied.

1      mesh Th=square(10,10);
2      int[int] ll=[4];
3      meshL ThL = extract(Th,[x+2,y*5],refedge=ll);


#### The commands rebuildBorder

This operator, used in the last example, allows to reconstruted the border elements following a special criteria ridgeangledetection. By default, it value is $$\frac{8}{9}*arctan(1)\approx40°$$, the diedral angle for a decahedron.

### The commands checkmesh

This function is avalaible for all 3D meshes. It checkes and validates the a given mesh, allows to remove duplicate vertices and/or elements and border elements. The possible arguments are

• precismesh= this parameter is the criteria to define two merging points.

By default, it value is 1e-7 and define the smallest axis parallel boxes containing the discretion domain of $$\Omega$$

• removeduplicate= is a boolean, allowing remove the duplicated elements and border elements

• rebuildboundary= is a boolean, allowing rebuild the border elements (in case of incomplete list given by the mesh)

Example:

1  mesh3 Th = checkmesh(Th);


### TetGen: A tetrahedral mesh generator

TetGen is a software developed by Dr. Hang Si of Weierstrass Institute for Applied Analysis and Stochastics in Berlin, Germany [HANG2006]. TetGen is free for research and non-commercial use. For any commercial license utilization, a commercial license is available upon request to Hang Si.

This software is a tetrahedral mesh generator of a three dimensional domain defined by its boundary (a surface). The input domain takes into account a polyhedral or a piecewise linear complex. This tetrahedralization is a constrained Delaunay tetrahedralization.

The method used in TetGen to control the quality of the mesh is a Delaunay refinement due to Shewchuk [SHEWCHUK1998]. The quality measure of this algorithm is the Radius-Edge Ratio (see Section 1.3.1 [HANG2006] for more details). A theoretical bound of this ratio of the Shewchuk algorithm is obtained for a given complex of vertices, constrained segments and facets of surface mesh, with no input angle less than 90 degrees. This theoretical bound is 2.0.

The launch of TetGen is done with the keyword tetg. The parameters of this command line is:

• reftet= sets the label of tetrahedra.

• label= is a vector of integers that contains the old labels number at index $$2i$$ and the new labels number at index $$2i+1$$ of Triangles.

This parameter is initialized as a label for the keyword change.

• switch= A string expression.

This string corresponds to the command line switch of TetGen see Section 3.2 of [HANG2006].

• nbofholes= Number of holes (default value: “size of holelist / 3”).

• holelist= This array corresponds to holelist of TetGenio data structure [HANG2006].

A real vector of size 3 * nbofholes. In TetGen, each hole is associated with a point inside this domain. This vector is $$x_{1}^{h}, y_{1}^{h}, z_{1}^{h}, x_{2}^{h}, y_{2}^{h}, z_{2}^{h}, \cdots,$$ where $$x_{i}^{h},y_{i}^{h},z_{i}^{h}$$ is the associated point with the $$i^{\mathrm{th}}$$ hole.

• nbofregions= Number of regions (default value: “size of regionlist / 5”).

• regionlist= This array corresponds to regionlist of TetGenio data structure [HANG2006].

The attribute and the volume constraint of region are given in this real vector of size 5 * nbofregions. The $$i^{\mathrm{th}}$$ region is described by five elements: $$x-$$coordinate, $$y-$$coordinate and $$z-$$coordinate of a point inside this domain ($$x_{i},y_{i},z_{i}$$); the attribute ($$at_{i}$$) and the maximum volume for tetrahedra ($$mvol_{i}$$) for this region.

The regionlist vector is: $$x_{1}, y_{1}, z_{1}, at_{1}, mvol_{1}, x_{2}, y_{2}, z_{2}, at_{2}, mvol_{2}, \cdots$$.

• nboffacetcl= Number of facets constraints “size of facetcl / 2”).

• facetcl= This array corresponds to facetconstraintlist of TetGenio data structure [HANG2006].

The $$i^{th}$$ facet constraint is defined by the facet marker $$Ref_{i}^{fc}$$ and the maximum area for faces $$marea_{i}^{fc}$$. The facetcl array is: $$Ref_{1}^{fc}, marea_{1}^{fc}, Ref_{2}^{fc}, marea_{2}^{fc}, \cdots$$.

This parameters has no effect if switch q is not selected.

Principal switch parameters in TetGen:

• p Tetrahedralization of boundary.

• q Quality mesh generation.

The bound of Radius-Edge Ratio will be given after the option q. By default, this value is 2.0.

• a Constructs with the volume constraints on tetrahedra.

These volumes constraints are defined with the bound of the previous switch q or in the parameter regionlist.

• A Attributes reference to region given in the regionlist.

The other regions have label 0.

The option AA gives a different label at each region. This switch works with the option p. If option r is used, this switch has no effect.

• r Reconstructs and Refines a previously generated mesh.

This character is only used with the command line tetgreconstruction.

• Y This switch preserves the mesh on the exterior boundary.

This switch must be used to ensure a conformal mesh between two adjacent meshes.

• YY This switch preserves the mesh on the exterior and interior boundary.

• C The consistency of the result’s mesh is testing by TetGen.

• CC The consistency of the result’s mesh is testing by TetGen and also constrained checks of Delaunay mesh (if p switch is selected) or the consistency of Conformal Delaunay (if q switch is selected).

• V Give information of the work of TetGen.

More information can be obtained in specified VV or VVV.

• Q Quiet: No terminal output except errors

• M The coplanar facets are not merging.

• T Sets a tolerance for coplanar test.

The default value is $$1e-8$$.

• d Intersections of facets are detected.

To obtain a tetrahedral mesh with TetGen, we need the surface mesh of a three dimensional domain. We now give the command line in FreeFEM to construct these meshes.

The keyword tetgtransfo

This keyword corresponds to a composition of command line tetg and movemesh23.

1 tetgtransfo(Th2, transfo=[Phi(1), Phi(2), Phi(3)]), ...) = tetg(Th3surf, ...),


where Th3surf = movemesh23(Th2, transfo=[Phi(1), Phi(2), Phi(3)]) and Th2 is the input two dimensional mesh of tetgtransfo.

The parameters of this command line are, on one hand, the parameters label, switch, regionlist, nboffacetcl, facetcl of keyword tetg and on the other hand, the parameter ptmerge of keyword movemesh23.

Note

To use tetgtransfo, the result’s mesh of movemesh23 must be a closed surface and define one region only. Therefore, the parameter regionlist is defined for one region.

An example of this keyword can be found in line 61 of the Build layer mesh example.

The keyword tetgconvexhull

FreeFEM, using TetGen, is able to build a tetrahedralization from a set of points. This tetrahedralization is a Delaunay mesh of the convex hull of the set of points.

The coordinates of the points can be initialized in two ways. The first is a file that contains the coordinate of points $$X_{i}=(x_{i}, y_{i}, z_{i})$$. This file is organized as follows:

$\begin{split}\begin{array}{ccc} n_{v} & & \\ x_{1} & y_{1} & z_{1} \\ x_{2} & y_{2} & z_{2} \\ \vdots &\vdots & \vdots \\ x_{n_v} & y_{n_v} & z_{n_v} \end{array}\end{split}$

The second way is to give three arrays that correspond respectively to the $$x-$$coordinates, $$y-$$coordinates and $$z-$$coordinates.

The parameters of this command line are :

• switch= A string expression.

This string corresponds to the command line switch of TetGen see Section 3.2 of [HANG2006].

• reftet= An integer expression.

Set the label of tetrahedra.

• label= An integer expression.

Set the label of triangles.

In the string switch, we can’t used the option p and q of TetGen.

### Reconstruct/Refine a 3d mesh with TetGen

Meshes in three dimension can be refined using TetGen with the command line tetgreconstruction.

The parameter of this keyword are

• region= an integer array that changes the region number of tetrahedra.

This array is defined as the parameter reftet in the keyword change.

• label= an integer array that changes the label of boundary triangles.

This array is defined as the parameter label in the keyword change.

• sizeofvolume= a reel function.

This function constraints the volume size of the tetrahedra in the domain (see Isotrope mesh adaption section to build a 3d adapted mesh).

The parameters switch, nbofregions, regionlist, nboffacetcl and facetcl of the command line which call TetGen (tetg) is used for tetgrefine.

In the parameter switch=, the character r should be used without the character p.

For instance, see the manual of TetGen [HANG2006] for effect of r to other character.

The parameter regionlist defines a new volume constraint in the region. The label in the regionlist will be the previous label of region.

This parameter and nbofregions can’t be used with the parameter sizeofvolume.

**Example refinesphere.edp**

 1 load "msh3"
4
5 mesh Th=square(10,20,[x*pi-pi/2,2*y*pi]);  //  $]\frac{-pi}{2},frac{-pi}{2}[\times]0,2\pi[$
6 //  a parametrization of a sphere
7 func f1 =cos(x)*cos(y);
8 func f2 =cos(x)*sin(y);
9 func f3 = sin(x);
10 //  partiel derivative of the parametrization DF
11 func f1x=sin(x)*cos(y);
12 func f1y=-cos(x)*sin(y);
13 func f2x=-sin(x)*sin(y);
14 func f2y=cos(x)*cos(y);
15 func f3x=cos(x);
16 func f3y=0;
17 // $M = DF^t DF$
18 func m11=f1x^2+f2x^2+f3x^2;
19 func m21=f1x*f1y+f2x*f2y+f3x*f3y;
20 func m22=f1y^2+f2y^2+f3y^2;
21
22 func perio=[[4,y],[2,y],[1,x],[3,x]];
23 real hh=0.1;
24 real vv= 1/square(hh);
25 verbosity=2;
28 plot(Th,wait=1);
29
30 verbosity=2;
31
32 // construction of the surface of spheres
33 real Rmin  = 1.;
34 func f1min = Rmin*f1;
35 func f2min = Rmin*f2;
36 func f3min = Rmin*f3;
37
38 meshS ThS=movemesh23(Th,transfo=[f1min,f2min,f3min]);
39
40 real[int] domain = [0.,0.,0.,145,0.01];
41 mesh3 Th3sph=tetg(ThS,switch="paAAQYY",nbofregions=1,regionlist=domain);
42
43 int[int] newlabel = [145,18];
44 real[int] domainrefine = [0.,0.,0.,145,0.0001];
45 mesh3 Th3sphrefine=tetgreconstruction(Th3sph,switch="raAQ",region=newlabel,nbofregions=1,regionlist=domainrefine,sizeofvolume=0.0001);
46
47 int[int] newlabel2 = [145,53];
48 func fsize = 0.01/(( 1 + 5*sqrt( (x-0.5)^2+(y-0.5)^2+(z-0.5)^2) )^3);
49 mesh3 Th3sphrefine2=tetgreconstruction(Th3sph,switch="raAQ",region=newlabel2,sizeofvolume=fsize);
50
51  medit("sphere",Th3sph,wait=1);
52  medit("sphererefinedomain",wait=1,Th3sphrefine);
53  medit("sphererefinelocal",wait=1,Th3sphrefine2);
54
55 // FFCS: testing 3d plots
56 plot(Th3sph);
57 plot(Th3sphrefine);
58 plot(Th3sphrefine2);


### 2d case

#### format of mesh data

Users who want to read a triangulation made elsewhere should see the structure of the file generated below:

1 border C(t=0, 2*pi){x=cos(t); y=sin(t);}
2 mesh Th1 = buildmesh(C(10));
3 savemesh(Th1, "mesh.msh");


The mesh is shown on Fig. 95.

The information about Th are saved in the file mesh.msh whose structure is shown on Table 1. An external file contains a mesh at format .mesh can be read by the ommand readmesh(file_name).

There, $$n_v$$ denotes the number of vertices, $$n_t$$ the number of triangles and $$n_s$$ the number of edges on boundary.

For each vertex $$q^i,\, i=1,\cdots,n_v$$, denoted by $$(q^i_x,q^i_y)$$ the $$x$$-coordinate and $$y$$-coordinate.

Each triangle $$T_k, k=1,\cdots,n_t$$ has three vertices $$q^{k_1},\, q^{k_2},\,q^{k_3}$$ that are oriented counter-clockwise.

The boundary consists of 10 lines $$L_i,\, i=1,\cdots,10$$ whose end points are $$q^{i_1},\, q^{i_2}$$.

In the Fig. 95, we have the following.

$$n_v=14, n_t=16, n_s=10$$

$$q^1=(-0.309016994375, 0.951056516295)$$

$$\dots$$

$$q^{14}=(-0.309016994375, -0.951056516295)$$

The vertices of $$T_1$$ are $$q^9, q^{12},\, q^{10}$$.

$$\dots$$

The vertices of $$T_{16}$$ are $$q^9, q^{10}, q^{6}$$.

The edge of the 1st side $$L_1$$ are $$q^6, q^5$$.

$$\dots$$

The edge of the 10th side $$L_{10}$$ are $$q^{10}, q^6$$.

Table 1 The structure of mesh_sample.msh

Content of the file

Explanation

14 16 10

$$n_v\quad n_t\quad n_e$$

-0.309016994375 0.951056516295 1

0.309016994375 0.951056516295 1

-0.309016994375 -0.951056516295 1

$$q^1_x\quad q^1_y\quad$$ boundary label $$=1$$

$$q^2_x\quad q^2_y\quad$$ boundary label $$=1$$

$$q^{14}_x\quad q^{14}_y\quad$$ boundary label $$=1$$

9 12 10 0

5 9 6 0

9 10 6 0

$$1_1\quad 1_2\quad 1_3\quad$$ region label $$=0$$

$$2_1\quad 2_2\quad 2_3\quad$$ region label $$=0$$

$$16_1\quad 16_2\quad 16_3\quad$$ region label $$=0$$

6 5 1

5 2 1

10 6 1

$$1_1\quad 1_2\quad$$ boundary label $$=1$$

$$2_1\quad 2_2\quad$$ boundary label $$=1$$

$$10_1\quad 10_2\quad$$ boundary label $$=1$$

In FreeFEM there are many mesh file formats available for communication with other tools such as emc2, modulef, … (see Mesh format chapter ).

The extension of a file implies its format. More details can be found on the file format .msh in the article by F. Hecht “bamg : a bidimensional anisotropic mesh generator” [HECHT1998_2].

A mesh file can be read into FreeFEM except that the names of the borders are lost and only their reference numbers are kept. So these borders have to be referenced by the number which corresponds to their order of appearance in the program, unless this number is overwritten by the keyword label. Here are some examples:

 1 // Parameters
2 int n = 10;
3
4 // Mesh
5 border floor(t=0, 1){x=t; y=0; label=1;};
6 border right(t=0, 1){x=1; y=t; label=5;};
7 border ceiling(t=1, 0){x=t; y=1; label=5;};
8 border left(t=1, 0){x=0; y=t; label=5;};
9
10 mesh th = buildmesh(floor(n) + right(n) + ceiling(n) + left(n));
11
12 //save mesh in different formats
13 savemesh(th, "toto.am_fmt"); // format "formated Marrocco"
14 savemesh(th, "toto.Th"); // format database db mesh "bamg"
15 savemesh(th, "toto.msh"); // format freefem
16 savemesh(th, "toto.nopo"); // modulef format
17
18 // Fespace
19 fespace femp1(th, P1);
20 femp1 f = sin(x)*cos(y);
21 femp1 g;
22
23 //save the fespace function in a file
24 {
25   ofstream file("f.txt");
26   file << f[] << endl;
27 } //the file is automatically closed at the end of the block
28 //read a file and put it in a fespace function
29 {
30   ifstream file("f.txt");
31   file >> g[] ;
32 }//the file is equally automatically closed
33
34 // Plot
35 plot(g);
36
37 // Mesh 2
40
41 // Fespace 2
42 fespace Vh2(th2, P1);
43 Vh2 u, v;
44
45 // Problem
46 //solve:
47 //  $u + \Delta u = g$ in $\Omega$
48 //  $u=0$ on $\Gamma_1$
49 //  $\frac{\partial u }{\partial n} = g$ on $\Gamma_2$
50 solve Problem(u, v)
51   = int2d(th2)(
52       u*v
53     - dx(u)*dx(v)
54     - dy(u)*dy(v)
55   )
56   + int2d(th2)(
57     - g*v
58   )
59   + int1d(th2, 5)(
60       g*v
61   )
62   + on(1, u=0)
63   ;
64
65 // Plot
66 plot(th2, u);


#### Input/output for a mesh

The function readmesh allows to build a mesh from a data file
1 mesh Th=readmeshS("Th.mesh");
2 mesh Thff = readmesh("Thff.msh"); // FreeFEM format

• the command savemesh

The function savemesh allows to export a mesh
1 savemesh(Th,"Th.mesh")
2 savemesh(Thff,"Thff.msh") // FreeFEM format
3
4 savemesh(th, "toto.msh"); // format freefem  savemesh(th, "toto.am_fmt"); // format "formated Marrocco"
5 savemesh(th, "toto.Th"); // format database db mesh "bamg"
6 savemesh(th, "toto.nopo"); // modulef format
7 // allows to save the 2d mesh with the 3rd space coordinate as a scalar solution for visualise
8 savemesh(Th,"mm",[x,y,u]); //  save surface mesh for medit, see for example minimal-surf.edp
9 exec("medit mm;rm mm.bb mm.faces mm.points");


The function vtkload allows to build a mesh from a data mesh at vtk format mesh

1 load "iovtk"

• the command savevtk

The function savevtk allows to export a mesh to a data mesh at vtk format mesh
1  load "iovtk"
2  savevtk("Th.vtk",Th);


The function gmshloadS allows to build a mesh from a data mesh file at formatmsh (GMSH)
1 load "gmsh"

• the command savegmsh

The function savegmsh allows to export a mesh to a data mesh msh (GMSH)
1     load "gmsh"
2 savegmsh(Th, "Th");


### 3d case

#### format of mesh data

In three dimensions, the file mesh format supported for input and output files by FreeFEM are the extension .msh and .mesh. These formats are described in the Mesh Format section.

Extension file .msh The structure of the files with extension .msh in 3D is given by:

$\begin{split}\begin{array}{cccccc} n_v & n_{tet} & n_{tri} & & \\ q^1_x & q^1_y & q^1_z & Vertex label & \\ q^2_x & q^2_y & q^2_z & Vertex label & \\ \vdots & \vdots & \vdots & \vdots & \\ q^{n_v}_x & q^{n_v}_y & q^{n_v}_z & Vertex label & \\ 1_1 & 1_2 & 1_3 & 1_4 & region label \\ 2_1 & 2_2 & 2_3 & 2_4 & region label \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ (n_{tet})_1 & (n_{tet})_2 & (n_{tet})_3 & (n_{tet})_4 & region label \\ 1_1 & 1_2 & 1_3 & boundary label & \\ 2_1 & 2_2 & 2_3 & boundary label & \\ \vdots & \vdots & \vdots & \vdots & \\ (n_tri)_{1} & (n_{tri})_2 & (n_{tri})_3 & boundary label & \\ \end{array}\end{split}$

In this structure, $$n_v$$ denotes the number of vertices, $$n_{tet}$$ the number of tetrahedra and $$n_{tri}$$ the number of triangles.

For each vertex $$q^i,\, i=1,\cdots,n_v$$, we denote by $$(q^i_x,q^i_y,q^i_z)$$ the $$x$$-coordinate, the $$y$$-coordinate and the $$z$$-coordinate.

Each tetrahedra $$T_k, k=1,\cdots,n_{tet}$$ has four vertices $$q^{k_1},\, q^{k_2},\,q^{k_3}, \,q^{k_4}$$.

The boundary consists of a union of triangles. Each triangle $$tri_j, j=1,\cdots,n_{tri}$$ has three vertices $$q^{j_1},\, q^{j_2},\,q^{j_3}$$.

extension file .mesh The data structure for a three dimensional mesh is composed of the data structure presented in Mesh Format section and a data structure for the tetrahedra. The tetrahedra of a three dimensional mesh are referred using the following field:

 1 Tetrahedra
2 NbTetrahedra
3 Vertex1 Vertex2 Vertex3 Vertex4 Label
4 ...
5 Vertex1 Vertex2 Vertex3 Vertex4 Label
6 Triangles
7 NbTriangles
8 Vertex1 Vertex2 Vertex3 Label
9 ...
10 Vertex1 Vertex2 Vertex3 Label


This field is express with the notation of Mesh Format section.

#### Input/output for a mesh3

The function readmesh3 allows to build a mesh3 from a data file
1 mesh3 Th3=readmesh3("Th3.mesh");
2 mesh3 Th3ff = readmesh3("Th3ff.msh"); // FreeFEM format

• the command savemesh

The function savemesh allows to export a mesh3
1 savemesh(Th3,"Th3.mesh")
2 savemesh(Th3ff,"Th3ff.msh") // FreeFEM format


The function vtkload3 allows to build a mesh3 from a data mesh at vtk format mesh

1 load "iovtk"

• the command savevtk

The function savevtk allows to export a mesh3 to a data mesh at vtk format mesh
1  load "iovtk"
2  savevtk("Th3.vtk",Th3);


The function gmshload3 allows to build a mesh3 from a data mesh file at formatmsh (GMSH)
1 load "gmsh"

• the command savegmsh

The function savegmsh allows to export a mesh3 to a data mesh msh (GMSH)
1     load "gmsh"
2 savegmsh(Th3, "Th3");


### Surface 3d case

#### format of mesh data

Like 2d and 3d, the input and output format files supported by FreeFEM are the extension .msh and .mesh. These formats are described in the Mesh Format section.

Extension file .msh The structure of the files with extension .msh in surface 3D is given by:

$\begin{split}\begin{array}{cccccc} n_v & n_{tri} & n_{edges} & & \\ q^1_x & q^1_y & q^1_z & Vertex label & \\ q^2_x & q^2_y & q^2_z & Vertex label & \\ \vdots & \vdots & \vdots & \vdots & \\ q^{n_v}_x & q^{n_v}_y & q^{n_v}_z & Vertex label & \\ 1_1 & 1_2 & 1_3 & region label \\ 2_1 & 2_2 & 2_3 & region label \\ \vdots & \vdots & \vdots & \vdots \\ (n_{tri})_1 & (n_{tri})_2 & (n_{tri})_3 & region label \\ 1_1 & 1_2 & boundary label & \\ 2_1 & 2_2 & boundary label & \\ \vdots & \vdots & \vdots & \\ (n_edge)_{1} & (n_{edge})_2 & boundary label & \\ \end{array}\end{split}$

In this structure, $$n_v$$ denotes the number of vertices, $$n_{tet}$$ the number of tetrahedra and $$n_{tri}$$ the number of triangles.

For each vertex $$q^i,\, i=1,\cdots,n_v$$, we denote by $$(q^i_x,q^i_y,q^i_z)$$ the $$x$$-coordinate, the $$y$$-coordinate and the $$z$$-coordinate.

Each tetrahedra $$T_k, k=1,\cdots,n_{tet}$$ has four vertices $$q^{k_1},\, q^{k_2},\,q^{k_3}, \,q^{k_4}$$.

The boundary consists of a union of triangles. Each triangle $$be_j, j=1,\cdots,n_{tri}$$ has three vertices $$q^{j_1},\, q^{j_2},\,q^{j_3}$$.

extension file .mesh The data structure for a three dimensional mesh is composed of the data structure presented in Mesh Format section and a data structure for the tetrahedra. The tetrahedra of a three dimensional mesh are referred using the following field:

 1 MeshVersionFormatted 2
2 Dimension 3
3
4 Vertices
5 NbVertices
6 (v0)x (v0)y (v0)z
7 ...
8 (vn)x (vn)y (vn)z
9
10 Triangles
11 NbTriangles
12 Vertex1 Vertex2 Vertex3 Label
13 ...
14 Vertex1 Vertex2 Vertex3 Label
15
16 Edges
17 NbEdges
18 Vertex1 Vertex2 Label
19 ...
20 Vertex1 Vertex2 Label
21
22 End


This field is express with the notation of Mesh Format section.

#### Input/output for a meshS

The function readmeshS allows to build a meshS from a data file
1 meshS ThS=readmeshS("ThS.mesh");
2 meshS Th3ff = readmeshS("ThSff.msh"); // FreeFEM format

• the command savemesh

The function savemesh allows to export a meshS
1 savemesh(ThS,"ThS.mesh")
2 savemesh(ThSff,"ThSff.msh") // FreeFEM format


The function vtkloadS allows to build a meshS from a data mesh at vtk format mesh

1 load "iovtk"

• the command savevtk

The function savevtk allows to export a meshS to a data mesh at vtk format mesh
1  load "iovtk"
2  savevtk("ThS.vtk",ThS);


The function gmshloadS allows to build a meshS from a data mesh file at formatmsh (GMSH)
1 load "gmsh"

• the command savegmsh

The function savegmsh allows to export a meshS to a data mesh msh (GMSH)
1     load "gmsh"
2 savegmsh(ThS, "ThS");


## Medit

The keyword medit allows to display a mesh alone or a mesh and one or several functions defined on the mesh using the Pascal Frey’s freeware medit. medit opens its own window and uses OpenGL extensively. Naturally to use this command medit must be installed.

A vizualisation with medit of scalar solutions $$f1$$ and $$f2$$ continuous, piecewise linear and known at the vertices of the mesh Th is obtained using:

1 medit("sol1 sol2", Th, f1, f2, order=1);


The first plot named sol1 display f1. The second plot names sol2 display f2.

The arguments of the function medit are the name of the differents scenes (separated by a space) of medit, a mesh and solutions.

Each solution is associated with one scene. The scalar, vector and symmetric tensor solutions are specified in the format described in the section dealing with the keyword savesol.

The parameters of this command line are :

• order= 0 if the solution is given at the center of gravity of elements.

1 is the solution is given at the vertices of elements.

• meditff= set the name of execute command of medit.

By default, this string is medit.

• save= set the name of a file .sol or .solb to save solutions.

This command line allows also to represent two differents meshes and solutions on them in the same windows. The nature of solutions must be the same. Hence, we can vizualize in the same window the different domains in a domain decomposition method for instance. A vizualisation with medit of scalar solutions $$h1$$ and $$h2$$ at vertices of the mesh Th1 and Th2 respectively are obtained using:

1 medit("sol2domain", Th1, h1, Th2, h2, order=1);


Tip

Medit

 1 load "medit"
2
3 // Initial Problem:
4 // Resolution of the following EDP:
5 // -Delta u_s = f on \Omega = { (x,y) | 1 <= sqrt(x^2+y^2) <= 2 }
6 // -Delta u_1 = f1 on \Omega_1 = { (x,y) | 0.5 <= sqrt(x^2+y^2) <= 1. }
7 // u = 1 on Gamma
8 // Null Neumman condition on Gamma_1 and on Gamma_2
9 // We find the solution u by solving two EDP defined on domain Omega and Omega_1
10 // This solution is visualize with medit
11
12 verbosity=3;
13
14 // Mesh
15 border Gamma(t=0, 2*pi){x=cos(t); y=sin(t); label=1;};
16 border Gamma1(t=0, 2*pi){x=2*cos(t); y=2*sin(t); label=2;};
17 border Gamma2(t=0, 2*pi){x=0.5*cos(t); y=0.5*sin(t); label=3;};
18
19 mesh Th = buildmesh(Gamma1(40) + Gamma(-40)); //Omega
20 mesh Th1 = buildmesh(Gamma(40) + Gamma2(-40)); //Omega_1
21
22 // Fespace
23 fespace Vh(Th, P2);
24 func f = sqrt(x*x + y*y);
25 Vh us, v;
26
27 fespace Vh1(Th1, P2);
28 func f1 = 10*sqrt(x*x+y*y);
29 Vh1 u1, v1;
30
31 // Macro
32 macro Grad2(us) [dx(us), dy(us)] // EOM
33
34 // Problem
35 problem Lap2dOmega (us, v, init=false)
36     = int2d(Th)(
38     )
39     - int2d(Th)(
40         f*v
41     )
42     +on(1, us=1)
43     ;
44
45 problem Lap2dOmega1 (u1, v1, init=false)
46     = int2d(Th1)(
48     )
49     - int2d(Th1)(
50         f1*v1
51     )
52     + on(1, u1=1)
53     ;
54
55 // Solve
56 Lap2dOmega;
57 Lap2dOmega1;
58
59 // Plot with medit
60 medit("solution", Th, us, Th1, u1, order=1, save="testsavemedit.solb");


## Mshmet

Mshmet is a software developed by P. Frey that allows to compute an anisotropic metric based on solutions (i.e. Hessian-based). This software can return also an isotropic metric. Moreover, mshmet can also construct a metric suitable for levelset interface capturing. The solution can be defined on 2D or 3D structured/unstructured meshes. For example, the solution can be an error estimate of a FE solution.

Solutions for mshmet are given as an argument. The solution can be a func, a vector func, a symmetric tensor, a fespace function, a fespace vector function and a fespace symmetric tensor. The symmetric tensor argument is defined as this type of data for datasol argument. This software accepts more than one solution.

For example, the metric $$M$$ computed with mshmet for the solution $$u$$ defined on the mesh $$Th$$ is obtained by writing:

1 fespace Vh(Th, P1);
2 Vh u; //a scalar fespace function
3 real[int] M = mshmet(Th, u);


The parameters of the keyword mshmet are :

• normalization = (b) do a normalization of all solution in $$[0,1]$$.

• aniso = (b) build anisotropic metric if 1 (default 0: isotropic)

• levelset = (b) build metric for levelset method (default: false)

• verbosity = (l) level of verbosity

• nbregul = (l) number of regularization’s iteration of solutions given (default 0).

• hmin = (d)

• hmax = (d)

• err = (d) level of error.

• width = (d) the width

• metric = a vector of double.

This vector contains an initial metric given to mshmet. The structure of the metric vector is described in the next paragraph.

• loptions = a vector of integer of size 7.

This vector contains the integer parameters of mshmet (for expert only).

• loptions(0): normalization (default 1).

• loptions(1): isotropic parameters (default 0).

1 for isotropic metric results otherwise 0.

• loptions(2): level set parameters (default 0).

1 for building level set metric otherwise 0.

• loptions(3): debug parameters (default 0).

1 for turning on debug mode otherwise 0.

• loptions(4): level of verbosity (default 10).

• loptions(5): number of regularization’s iteration of solutions given (default 0).

• loptions(6): previously metric parameter (default 0).

1 for using previous metric otherwise 0.

• doptions= a vector of double of size 4.

This vector contains the real parameters of mshmet (for expert only).

• doptions(0): hmin : min size parameters (default 0.01).

• doptions(1): hmax : max size parameters (default 1.0).

• doptions(2): eps : tolerance parameters (default 0.01).

• doptions(2): width : relative width for Level Set ($$0<w<1$$) (default 0.05).

The result of the keyword mshmet is a real[int] which contains the metric computed by mshmet at the different vertices $$V_{i}$$ of the mesh.

With $$nv$$ is the number of vertices, the structure of this vector is:

$M_{iso}= (m(V_0), m(V_1), \ldots, m(V_{nv}))^t$

for a isotropic metric $$m$$. For a symmetric tensor metric $$h=\left(\begin{array}{ccc} m_{1 1} & m_{1 2} & m_{1 3}\\ m_{2 1} & m_{2 2} & m_{2 3} \\ m_{3 1} & m_{3 2} & m_{3 3} \end{array}\right)$$ , the parameters metric is:

$M_{aniso}= (H(V_{0}), \ldots, H(V_{nv}) )^t$

where $$H(V_{i})$$ is the vector of size 6 defined by $$[m11,m21,m22,m31,m32,m33]$$

Tip

mshmet

 1 load "mshmet"
4
5 // Parameters
6 real error = 0.01;
7 func zmin = 0;
8 func zmax = 1;
9 int MaxLayer = 10;
10
11 // Mesh
12 border a(t=0, 1.0){x=t; y=0; label=1;};
13 border b(t=0, 0.5){x=1; y=t; label=2;};
14 border c(t=0, 0.5){x=1-t; y=0.5; label=3;};
15 border d(t=0.5, 1){x=0.5; y=t; label=4;};
16 border e(t=0.5, 1){x=1-t; y=1; label=5;};
17 border f(t=0.0, 1){x=0; y=1-t; label=6;};
18 mesh Th = buildmesh(a(6) + b(4) + c(4) + d(4) + e(4) + f(6));
19 mesh3 Th3 = buildlayers(Th, MaxLayer, zbound=[zmin, zmax]);
20
21 // Fespace
22 fespace Vh3(Th3, P2);
23 Vh3 u3, v3;
24
25 fespace Vh3P1(Th3, P1);
26 Vh3P1 usol;
27
28 // Problem
29 problem Problem2(u3, v3, solver=sparsesolver)
30     = int3d(Th3)(
31         u3*v3*1.0e-10
32         + dx(u3)*dx(v3)
33         + dy(u3)*dy(v3)
34         + dz(u3)*dz(v3)
35     )
36     - int3d(Th3)(
37         v3
38     )
39     +on(0, 1, 2, 3, 4, 5, 6, u3=0)
40     ;
41
42 // Solve
43 Problem2;
44 cout << u3[].min << " " << u3[].max << endl;
45
46 medit("Sol", Th3, u3);
47
48 real[int] bb = mshmet(Th3,u3);
49 cout << "Metric:" << bb << endl;
50 for (int ii = 0; ii < Th3.nv; ii++)
51     usol[][ii] = bb[ii];
52
53 medit("Metric", Th3, usol);


## FreeYams

FreeYams is a surface mesh adaptation software which is developed by P. Frey. This software is a new version of yams. The adapted surface mesh is constructed with a geometric metric tensor field. This field is based on the intrinsic properties of the discrete surface.

Also, this software allows to construct a simplification of a mesh. This decimation is based on the Hausdorff distance between the initial and the current triangulation. Compared to the software yams, FreeYams can be used also to produce anisotropic triangulations adapted to levelset simulations. A technical report on freeYams documentation is available here.

To call FreeYams in FreeFEM, we used the keyword freeyams. The arguments of this function are the initial mesh and/or metric. The metric with freeyams are a func, a fespace function, a symmetric tensor function, a symmetric tensor fespace function or a vector of double (real[int]). If the metric is a vector of double, this data must be given in metric parameter. Otherwise, the metric is given in the argument.

For example, the adapted mesh of Thinit defined by the metric $$u$$ defined as fespace function is obtained by writing:

1 fespace Vh(Thinit, P1);
2 Vh u;
3 mesh3 Th = freeyams(Thinit, u);


The symmetric tensor argument for freeyams keyword is defined as this type of data for datasol argument.

• aniso= (b) aniso or iso metric (default 0, iso)

• mem= (l) memory of for freeyams in Mb (default -1, freeyams choose)

• hmin= (d)

• hmax= (d)

• gradation= (d)

• option= (l)

• 0 : mesh optimization (smoothing+swapping)

• 1 : decimation+enrichment adaptated to a metric map. (default)

• -1 : decimation adaptated to a metric map.

• 2 : decimation+enrichment with a Hausdorff-like method

• -2 : decimation with a Hausdorff-like method

• 4 : split triangles recursively.

• 9 : No-Shrinkage Vertex Smoothing

• ridgeangle= (d)

• absolute= (b)

• verbosity= (i)

• metric= vector expression.

This parameters contains the metric at the different vertices on the initial mesh. With $$nv$$ is the number of vertices, this vector is:

$M_{iso}= ( m(V_0), m(V_1), \ldots, m(V_{nv}) )^t$

for a scalar metric $$m$$. For a symmetric tensor metric $$h=\left(\begin{array}{ccc} m_{1 1} & m_{1 2} & m_{1 3}\\ m_{2 1} & m_{2 2} & m_{2 3} \\ m_{3 1} & m_{3 2} & m_{3 3} \end{array}\right)$$, the parameters metric is:

$M_{aniso}= ( H(V_{0}), \ldots, H(V_{nv}) )^t$

where $$H(V_{i})$$ is the vector of size 6 defined by $$[m11,m21,m22,m31,m32,m33]$$

• loptions= a vector of integer of size 13.

This vectors contains the integer options of FreeYams. (just for the expert)

• loptions(0): anisotropic parameter (default 0).

If you give an anisotropic metric 1 otherwise 0.

• loptions(1): Finite Element correction parameter (default 0).

1 for no Finite Element correction otherwise 0.

• loptions(2): Split multiple connected points parameter (default 1).

1 for splitting multiple connected points otherwise 0.

• loptions(3): maximum value of memory size in Mbytes (default -1: the size is given by freeyams).

• loptions(4): set the value of the connected component which we want to obtain.

(Remark: freeyams give an automatic value at each connected component).

• loptions(5): level of verbosity

• loptions(6): Create point on straight edge (no mapping) parameter (default 0).

1 for creating point on straight edge otherwise 0.

• loptions(7): validity check during smoothing parameter.

This parameter is only used with No-Shrinkage Vertex Smoothing optimization (optimization option parameter 9). 1 for No validity checking during smoothing otherwise 0.

• loptions(8): number of desired’s vertices (default -1).

• loptions(9): number of iteration of optimizations (default 30).

• loptions(10): no detection parameter (default 0).

1 for detecting the ridge on the mesh otherwise 0. The ridge definition is given in the parameter doptions(12).

• loptions(11): no vertex smoothing parameter (default 0).

1 for smoothing the vertices otherwise 0.

• loptions(12): Optimization level parameter (default 0).

• 0 : mesh optimization (smoothing+swapping)

• 1 : decimation+enrichment adaptated to a metric map.

• -1: decimation adaptated to a metric map.

• 2 : decimation+enrichment with a Hausdorff-like method

• -2: decimation with a Hausdorff-like method

• 4 : split triangles recursively.

• 9 : No-Shrinkage Vertex Smoothing

• doptions= a vector of double of size 11.

This vectors contains the real options of freeyams.

• doptions(0): Set the geometric approximation (Tangent plane deviation) (default 0.01).

• doptions(1): Set the lamda parameter (default -1).

• doptions(2): Set the mu parmeter (default -1).

• doptions(3): Set the gradation value (Mesh density control) (default 1.3).

• doptions(4): Set the minimal size(hmin) (default -2.0: the size is automatically computed).

• doptions(5): Set the maximal size(hmax) (default -2.0: the size is automatically computed).

• doptions(6): Set the tolerance of the control of Chordal deviation (default -2.0).

• doptions(7): Set the quality of degradation (default 0.599).

• doptions(8): Set the declic parameter (default 2.0).

• doptions(9): Set the angular walton limitation parameter (default 45 degree).

• doptions(10): Set the angular ridge detection (default 45 degree).

Tip

freeyams

 1 load "msh3"
4
5 // Parameters
6 int nn = 20;
7 real zmin = 0;
8 real zmax = 1;
9
10 // Mesh
11 mesh Th2 = square(nn, nn);
12 int[int] rup = [0, 2], rdown = [0, 1];
13 int[int] rmid = [1, 1, 2, 1, 3, 1, 4, 1];
14 mesh3 Th = buildlayers(Th2, nn, zbound=[zmin, zmax], reffacemid=rmid, reffaceup=rup, reffacelow=rdown);
15 mesh3 Th3 = freeyams(Th);
16
17 medit("SurfaceMesh", Th3);


## mmg3d

Todo

mmg3d-v4.0

Mmg3d is a 3D remeshing software developed by C. Dobrzynski and P. Frey.

This software allows to remesh an initial mesh made of tetrahedra. This initial mesh is adapted to a geometric metric tensor field or to a displacement vector (moving rigid body). The metric can be obtained with mshmet.

Note

• If no metric is given, an isotropic metric is computed by analyzing the size of the edges in the initial mesh.

• If a displacement is given, the vertices of the surface triangles are moved without verifying the geometrical structure of the new surface mesh.

The parameters of mmg3d are :

• options= vector expression.

This vector contains the option parameters of mmg3d. It is a vector of 6 values, with the following meaning:

• Optimization parameters : (default 1)

0 : mesh optimization.

1 : adaptation with metric (deletion and insertion vertices) and optimization.

-1 : adaptation with metric (deletion and insertion vertices) without optimization.

4 : split tetrahedra (be careful modify the surface).

9 : moving mesh with optimization.

-9 : moving mesh without optimization.

• Debug mode : (default 0)

1 : turn on debug mode.

0 : otherwise.

• Specify the size of bucket per dimension (default 64)

• Swapping mode : (default 0)

1 : no edge or face flipping.

0 : otherwise.

• Insert points mode : (default 0)

1 : no edge splitting or collapsing and no insert points.

0 : otherwise.

1. Verbosity level (default 3)

• memory= integer expression.

Set the maximum memory size of new mesh in Mbytes. By default the number of maximum vertices, tetrahedra and triangles are respectively 500 000, 3000 000, 100000 which represent approximately a memory of 100 Mo.

• metric= vector expression.

This vector contains the metric given at mmg3d. It is a vector of size $$nv$$ or 6 $$nv$$ respectively for an isotropic and anisotropic metric where $$nv$$ is the number of vertices in the initial mesh. The structure of metric vector is described in the mshmet.

• displacement= $$[\Phi1, \Phi2, \Phi3]$$ set the displacement vector of the initial mesh $$\mathbf{\Phi(x,y)} = [\Phi1(x,y), \Phi2(x,y), \Phi3(x,y)]$$.

• displVect= sets the vector displacement in a vector expression.

This vector contains the displacement at each point of the initial mesh. It is a vector of size 3 $$nv$$.

Tip

mmg3d

 1 load "msh3"
4 include "Cube.idp"
5
6 // Parameters
7 int n = 6;
8 int[int] Nxyz = [12, 12, 12];
9 real [int, int] Bxyz = [[0., 1.], [0., 1.], [0., 1.]];
10 int [int, int] Lxyz = [[1, 1], [2, 2], [2, 2]];
11
12 // Mesh
13 mesh3 Th = Cube(Nxyz, Bxyz, Lxyz);
14
15 real[int] isometric(Th.nv);
16 for (int ii = 0; ii < Th.nv; ii++)
17     isometric[ii] = 0.17;
18
19 mesh3 Th3 = mmg3d(Th, memory=100, metric=isometric);
20
21 // Plot
22 medit("Initial", Th);
23 medit("Isometric", Th3);


Tip

Falling spheres

 1 load "msh3"
5 include "MeshSurface.idp"
6
7 // Parameters
8 real hs = 0.8;
9 int[int] N = [4/hs, 8/hs, 11.5/hs];
10 real [int, int] B = [[-2, 2], [-2, 6], [-10, 1.5]];
11 int [int, int] L = [[311, 311], [311, 311], [311, 311]];
12
13 int[int] opt = [9, 0, 64, 0, 0, 3];
14 real[int] vit=[0, 0, -0.3];
15 func zero = 0.;
16 func dep = vit[2];
17
18 // Meshes
19 meshS ThH = SurfaceHex(N, B, L, 1);
20 meshS ThSg = Sphere(1, hs, 300, -1);
21 meshS ThSd = Sphere(1, hs, 310, -1);
22 ThSd = movemesh(ThSd, [x, 4+y, z]);
23 meshS ThHS = ThH + ThSg + ThSd;
24 medit("ThHS", ThHS);
25
26 real voltet = (hs^3)/6.;
27 real[int] domain = [0, 0, -4, 1, voltet];
28 real [int] holes = [0, 0, 0, 0, 4, 0];
29 mesh3 Th = tetg(ThHS, switch="pqaAAYYQ", nbofregions=1, regionlist=domaine, nbofholes=2, holelist=holes);
30 medit("Box-With-two-Ball", Th);
31
32 // Fespace
33 fespace Vh(Th, P1);
34 Vh uh,vh;
35
36 // Macro
38
39 // Problem
40 problem Lap (uh, vh, solver=CG)
41     = int3d(Th)(
43     )
44     + on(310, 300, uh=dep)
45     + on(311, uh=0.)
46     ;
47
48 // Falling loop
49 for(int it = 0; it < 29; it++){
50     cout << " ITERATION " << it << endl;
51
52     // Solve
53     Lap;
54
55     // Plot
56     plot(Th, uh);
57
58     // Sphere falling
59     Th = mmg3d(Th, options=opt, displacement=[zero, zero, uh], memory=1000);
60 }


## A first 3d isotrope mesh adaptation process

Tip

 1 load "msh3"
5
6 // Parameters
7 int nn = 6;
8 int[int] l1111 = [1, 1, 1, 1]; //labels
9 int[int] l01 = [0, 1];
10 int[int] l11 = [1, 1];
11
12 real errm = 1e-2; //level of error
13
14 // Mesh
15 mesh3 Th3 = buildlayers(square(nn, nn, region=0, label=l1111),
16 nn, zbound=[0, 1], labelmid=l11, labelup=l01, labeldown=l01);
17
18 Th3 = trunc(Th3, (x<0.5) | (y < 0.5) | (z < 0.5), label=1); //remove the ]0.5,1[^3 cube
19
20 // Fespace
21 fespace Vh(Th3, P1);
22 Vh u, v, usol, h;
23
24 // Macro
25 macro Grad(u) [dx(u), dy(u), dz(u)] // EOM
26
27 // Problem
28 problem Poisson (u, v, solver=CG)
29     = int3d(Th3)(
31     )
32     - int3d(Th3)(
33         1*v
34     )
35     + on(1, u=0)
36     ;
37
38 // Loop
39 for (int ii = 0; ii < 5; ii++){
40     // Solve
41     Poisson;
42     cout << "u min, max = " << u[].min << " "<< u[].max << endl;
43
44     h=0.; //for resizing h[] because the mesh change
45     h[] = mshmet(Th3, u, normalization=1, aniso=0, nbregul=1, hmin=1e-3, hmax=0.3, err=errm);
46     cout << "h min, max = " << h[].min << " "<< h[].max << " " << h[].n << " " << Th3.nv << endl;
47     plot(u, wait=true);
48
49     errm *= 0.8; //change the level of error
50     cout << "Th3 " << Th3.nv < " " << Th3.nt << endl;
51     Th3 = tetgreconstruction(Th3, switch="raAQ", sizeofvolume=h*h*h/6.); //rebuild mesh
53 }


## Build a 2d mesh from an isoline

The idea is to get the discretization of an isoline of fluid meshes, this tool can be useful to construct meshes from image. First, we give an example of the isovalue meshes $$0.2$$ of analytical function $$\sqrt{(x-1/2)^2 +(y-1/2)^2}$$, on unit square.

 1 load "isoline"
2
3 real[int,int] xy(3, 1); //to store the isoline points
4 int[int] be(1); //to store the begin, end couple of lines
5 {
6     mesh Th = square(10, 10);
7     fespace Vh(Th, P1);
8     Vh u = sqrt(square(x-0.5) + square(y-0.5));
9     real iso = 0.2 ;
10     real[int] viso = [iso];
11     plot(u, viso=viso,Th);//to see the iso line
12
13     int nbc = isoline(Th, u, xy, close=1, iso=iso, beginend=be, smoothing=0.1);


The isoline parameters are Th the mesh, the expression $$u$$, the bidimentionnal array xy to store the list coordinate of the points. The list of named parameter are :

• iso= value of the isoline to compute (0 is the default value)

• close= close the isoline with the border (default true), we add the part of the mesh border such the value is less than the isovalue

• smoothing= number of smoothing process is the $${l} ^{r} {s}$$ where $$l$$ is the length of the current line component, $$r$$ the ratio, $$s$$ is smoothing value. The smoothing default value is 0.

• ratio= the ratio (1 by default).

• eps= relative $$\varepsilon$$ (default 1e-10)

• beginend= array to get begin, end couple of each of sub line (resize automatically)

• file= to save the data curve in data file for gnuplot

In the array xy you get the list of vertices of the isoline, each connex line go from $$i= i_0^c ,\dots, i_1^c-1$$ with $$i_0^c =be(2*c)$$ $$i_1^c =be(2*c+1)$$, and where $$x_i= xy(0,i), y_i=yx( 1,i), l_i=xy(2,i)$$.

Here $$l_i$$ is the length of the line (the origin of the line is point $$i_0^c$$).

The sense of the isoline is such that the upper part is at the left size of the isoline. So here : the minimum is a point $$0.5,05$$ so the curve 1 turn in the clockwise sense, the order of each component are sort such that the number of point by component is decreasing.

 1       cout << "Number of the line component = " << nbc << endl;
2       cout << "Number of points = " << xy.m << endl;
3       cout << "be = " << be << endl;
4
5       // shows the lines component
6       for (int c = 0; c < nbc; ++c){
7          int i0 = be[2*c], i1 = be[2*c+1]-1;
8          cout << "Curve " << c << endl;
9          for(int i = i0; i <= i1; ++i)
10             cout << "x= " << xy(0,i) << " y= " << xy(1,i) << " s= " << xy(2, i) << endl;
11          plot([xy(0, i0:i1), xy(1, i0:i1)], wait=true, viso=viso, cmm=" curve "+c);
12    }
13 }
14
15 cout << "length of last curve = " << xy(2, xy.m-1) << endl;


We also have a new function to easily parametrize a discrete curve defined by the couple $$be, xy$$.

 1 border Curve0(t=0, 1){
2     int c=0; //component 0
3     int i0=be[2*c], i1=be[2*c+1]-1;
4     P=Curve(xy, i0, i1, t); //Curve 0
5     label=1;
6 }
7
8 border Curve1(t=0, 1){
9     int c=1; //component 1
10     int i0=be[2*c], i1=be[2*c+1]-1;
11     P=Curve(xy, i0, i1, t); //Curve 1
12     label=1;
13 }
14
15 plot(Curve1(100)); //show curve
16 mesh Th = buildmesh(Curve1(-100));
17 plot(Th, wait=true);


Secondly, we use this idea to build meshes from an image, we use the plugins ppm2rnm to read pgm a gray scale image and then we extract the gray contour at level 0.25.

Tip

Leman lake

 1 load "ppm2rnm"
3
4 // Parameters
5 string leman = "LemanLake.pgm";
6 real AreaLac = 580.03; //in km^2
7 real hsize = 5;
8 real[int, int] Curves(3, 1);
9 int[int] be(1);
10 int nc; //nb of curve
11 {
12     real[int, int] ff1(leman); //read image
13     //and set it in a rect. array
14     int nx = ff1.n, ny = ff1.m;
15     //build a Cartesian mesh such that the origin is in the right place.
16     mesh Th = square(nx-1, ny-1, [(nx-1)*(x), (ny-1)*(1-y)]);
17     //warning the numbering of the vertices (x,y) is
18     //given by $i = x/nx + nx* y/ny$
19     fespace Vh(Th, P1);
20     Vh f1;
21     f1[] = ff1; //transform array in finite element functions.
22     nc = isoline(Th, f1, iso=0.25, close=1, Curves, beginend=be, smoothing=.1, ratio=0.5);
23 }
24
25 //The longest isoline: the lake
26 int ic0 = be(0), ic1 = be(1)-1;
27 plot([Curves(0, ic0:ic1), Curves(1, ic0:ic1)], wait=true);
28
29 int NC = Curves(2, ic1)/hsize;
30 real xl = Curves(0, ic0:ic1).max - 5;
31 real yl = Curves(1, ic0:ic1).min + 5;
32 border G(t=0, 1){P=Curve(Curves, ic0, ic1, t); label=1+(x>xl)*2+(y<yl);}
33 plot(G(-NC), wait=true);
34
35 mesh Th = buildmesh(G(-NC));
36 plot(Th, wait=true);
37
38 real scale = sqrt(AreaLac/Th.area);
39 Th = movemesh(Th, [x*scale, y*scale]);
40 cout << "Th.area = " << Th.area << " Km^2 " << " == " << AreaLac << " Km^2 " << endl;
41 plot(Th, wait=true, ps="leman.eps");


Isoline