Mesh Generation

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.

1mesh 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:

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

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

1for (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:

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

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

1mesh 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.

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

Mesh with a hole

1real r=1;
2border a(t=0, 2*pi){x=r*cos(t); y=r*sin(t); label=1;}
3r=0.3;
4border b(t=0, 2*pi){x=r*cos(t); y=r*sin(t); label=1;}
5mesh 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:

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

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

 1real[int] xx = [0, 1, 1, 0],
 2          yy = [0, 0, 1, 1];
 3//radius, center of the 4 circles
 4real[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];
 7int[int] NC = [-10,-11,-12,13]; //list number of :math:`\pm` segments of the 4 circles borders
 8
 9border 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
19border 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}
25int[int] nn = [4, 4, 5, 7]; //4 border, with 4, 4, 5, 7 segment respectively
26plot(bb(nn), cc(NC), wait=1);
27mesh th = buildmesh(bb(nn) + cc(NC));
28plot(th, wait=1);

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:

1mesh 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:

for a big number \(k>1\).

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

Move mesh

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

 1// Parameters
 2int nn = 10;
 3real dt = 0.1;
 4
 5// Mesh
 6mesh Th = square(nn, nn);
 7
 8// Fespace
 9fespace Vh(Th, P1);
10Vh u=y;
11
12// Loop
13real t=0;
14for (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

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

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

The command adaptmesh

The function:

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
 2real eps = 0.0001;
 3real h = 1;
 4real hmin = 0.05;
 5func f = 10.0*x^3 + y^3 + h*atan2(eps, sin(5.0*y)-2.0*x);
 6
 7// Mesh
 8mesh Th = square(5, 5, [-1+2*x, -1+2*y]);
 9
10// Fespace
11fespace Vh(Th,P1);
12Vh fh = f;
13plot(fh);
14
15// Adaptmesh
16for (int i = 0; i < 2; i++){
17    Th = adaptmesh(Th, fh);
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:

1mesh ATh = adaptmesh(Th, f);

create the new mesh ATh adapted to the Hessian

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.

Mesh adaptation

 1// Parameters
 2real error = 0.1;
 3
 4// Mesh
 5border ba(t=0, 1){x=t; y=0; label=1;}
 6border bb(t=0, 0.5){x=1; y=t; label=1;}
 7border bc(t=0, 0.5){x=1-t; y=0.5; label=1;}
 8border bd(t=0.5, 1){x=0.5; y=t; label=1;}
 9border be(t=0.5, 1){x=1-t; y=1; label=1;}
10border bf(t=0, 1){x=0; y=1-t; label=1;}
11mesh Th = buildmesh(ba(6) + bb(4) + bc(4) + bd(4) + be(4) + bf(6));
12
13// Fespace
14fespace Vh(Th, P1);
15Vh u, v;
16
17// Function
18func f = 1;
19
20// Problem
21problem 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
31// Adaptmesh loop
32for (int i = 0; i < 4; i++){
33    Poisson;
34    Th = adaptmesh(Th, u, err=error);
35    error = error/2;
36}
37
38// Plot
39plot(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:

See Reference part for more inforamtions

• 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:

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

Mesh adaptation

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:

1mesh 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
 2mesh Th = square(3, 3);
 3
 4// Fespace
 5fespace Vh(Th, P1);
 6Vh u=0;
 7
 8// Loop on all degrees of freedom
 9int n=u.n;
10for (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}