Literature survey of acupuncture study
2.2. FREEFEM++ AND ITS INTERPRETED LANGUAGE
2.2.2 Meshing tools and mesh exemples
The initial step, as in all finite element method numerical simulations, is to discretize a geometry Ω into the associate mesh, that is in FreeFem++ a triangulation Th. Two-dimensional and three-Two-dimensional meshes can be generated directly using the FreeFem++
built-in tools. FreeFem++ can also import and export mesh files that are in the FreeFem++
required format. The tools savemesh and readmesh are used to save or load, respec-tively, a mesh file ".msh". Third-party software such as Gmsh [73] or TetGen [74] can be used together with FreeFem++.
2.2.2.1 Uniform square mesh
A square mesh can be defined with the built-in function square. It allows the user to triangulate a square domain and generate the mesh from a set of boundary points. The code in script 2.2 can be written to build a mesh of square domain. The keyword plot is used to plot the mesh (see figure 2.1). In the same way, the code in script 2.3 can be written to build a mesh of rectangle domain plotted in figure 2.1.
1 i n t nn = 1 0 ; / / number of points on each border
2 mesh Th = s q u a r e ( nn , nn ) ;
2.2. FREEFEM++ AND ITS INTERPRETED LANGUAGE
3 p l o t( Th ) ;
Script 2.2: Script to build a mesh of square domain plotted in figure 2.1.
1 i n t nn = 1 0 ;
2 r e a l x0 = 0 , x1 = 1 0 ;
3 r e a l y0 = 0 , y1 = 5 ;
4 mesh Th=s q u a r e( nn , 2 ∗ nn , [ x0 + ( x1−x0 ) ∗x, y0 + ( y1−y0 ) ∗y] ) ;
5 p l o t( Th ) ;
Script 2.3: Script to build the mesh of a rectangle domain plotted in figure 2.1.
Figure 2.1: Uniform square mesh and rectangle mesh generated by the scripts 2.2 and 2.3, respectively
2.2.2.2 Analytic description of a domain
The geometry can also be defined by an analytic description of boundaries by pieces from which the mesh is automatically generated. Each piece of boundary is described by its parametric equation and is declared with the border type. The keyword label labels one curve or a group of curves of the boundary. The label type can be either an integer or a name. The mesh is then generated by calling the tool buildmesh that will automatically generate the Delaunay triangulation from a set of boundary points (see script 2.4). An example of generated mesh is shown in figure 2.2.
Note that the orientation of the curve is important in order for the domain to be defined on the correct side of the curve. In the buildmesh command, a number of vertices must be provided each piece of boundary. These numbers are positive or negative. The change of sign changes the orientation of the curve. It should also be noted that the boundaries can only intersect at their end points.
2.2. FREEFEM++ AND ITS INTERPRETED LANGUAGE
1 b o r d e r b1 ( t = 0 , 2 ∗p i / 5 . )
2 {x= ( 6 ∗c o s( t )−c o s( 6 ∗ t ) ) ;y= ( 6 ∗s i n( t )−s i n( 6 ∗ t ) ) ; l a b e l= 1 ; }
3 b o r d e r b2 ( t =2∗p i/ 5 . , 4 ∗p i / 5 . )
4 {x= ( 6 ∗c o s( t )−c o s( 6 ∗ t ) ) ;y= ( 6 ∗s i n( t )−s i n( 6 ∗ t ) ) ; l a b e l= 1 ; }
5 b o r d e r b3 ( t =4∗p i/ 5 . , 6 ∗p i / 5 . )
6 {x= ( 6 ∗c o s( t )−c o s( 6 ∗ t ) ) ;y= ( 6 ∗s i n( t )−s i n( 6 ∗ t ) ) ; l a b e l= 1 ; }
7 b o r d e r b4 ( t =6∗p i/ 5 . , 8 ∗p i / 5 . )
8 {x= ( 6 ∗c o s( t )−c o s( 6 ∗ t ) ) ;y= ( 6 ∗s i n( t )−s i n( 6 ∗ t ) ) ; l a b e l= 1 ; }
9 b o r d e r b5 ( t =8∗p i/ 5 . , 1 0 ∗p i / 5 . )
10 {x= ( 6 ∗c o s( t )−c o s( 6 ∗ t ) ) ;y= ( 6 ∗s i n( t )−s i n( 6 ∗ t ) ) ; l a b e l= 1 ; }
11 i n t np = 2 5 ; / / mesh resolution
12 mesh Th=b u i l d m e s h( b1 ( np ) +b1 ( np ) +b1 ( np ) +b1 ( np ) +b1 ( np ) ) ;
13 p l o t( Th ) ;
14 s a v e m e s h( Th , " m e s h f i l e . msh " ) ;
Script 2.4: Script to build the mesh plotted in figure 2.2.
Ranunculoid mesh
Figure 2.2: Mesh generated by the script in 2.4 representing a plum blossom of the Taiwan national flower
2.2.2.3 Mesh construction from an image
FreeFem++ embeds tools that allow the generation of meshes from a given image.
The code in script 2.5 illustrates the capability to generate a two-dimensional mesh of Taiwan and a three-dimensional mesh of the elevation of Taiwan (see figures 2.3 and 2.4).
2.2. FREEFEM++ AND ITS INTERPRETED LANGUAGE
2.2. FREEFEM++ AND ITS INTERPRETED LANGUAGE
42 d e e p = d e e p ∗ 2 0 /a b s( d e e p [ ] .max) ;
43 p l o t( deep ,w a i t= 1 ) ;
44 i n t nn = 3 0 ;
45 r e a l maxdeep = d e e p [ ] .max;
46 i n t[ i n t] r u p = [ 0 , 2 0 0 ] , rdown = [ 0 , 1 0 0 ] , r m i d ( 1 7 ∗ 2 ) ;
47 f o r( i n t i = 0 ; i < r m i d . n ; + + i ) r m i d [ i ]=1+ i / 2 ;
48 mesh3 Th3=b u i l d l a y e r s( Sh , nn , / / 3d mesh with layers
49 c o e f= d e e p / maxdeep , z b o u n d= [ 0 , d e e p ] , l a b e l m i d= rmid ,
50 r e f f a c e u p = r u p , r e f f a c e l o w = rdown ) ;
51 m e d i t( " Th3 " , Th3 ) ; / / plot in medit
52 s a v e m e s h( Th3 , " Taiwan3D . mesh " ) ;
Script 2.5: Script to build 2d and 3d meshes from an image
The keyword load is used to include an external library such as msh3 to handle three-dimensional meshes and medit to visualize 3d meshes and results in the software MEDIT [75]. FreeFem++ natively includes three-dimensional visualization but MEDIT has advanced visualization tools.
Figure 2.3: Meshes of Taiwan generated from the image on the left
2.2.2.4 Mesh manipulation
Truncating and adding meshes Another way to build meshes is to use the built-in truncation tool or the addition operator. To truncate a mesh, i.e. to remove triangles from the mesh, one can use the function trunc. The code in script 2.6 gives an exemple of the use of the function trunc. The resulting mesh is plotted in figure 2.5.
1 mesh Th=s q u a r e( nn , l y / l x ∗ nn , [ x0 + ( x1−x0 ) ∗x, y0 + ( y1−y0 ) ∗y] ) ;
2.2. FREEFEM++ AND ITS INTERPRETED LANGUAGE
Figure 2.4: 3d mesh of Taiwan generated from the elevation map on the left
2 mesh Th1=t r u n c( Th , (max(a b s(x) ,a b s(y) ) >=1) ,l a b e l= 5 ) ; / / new mesh and label number of new boundary
Script 2.6: Script to build a truncated mesh. An exemple is plotted in figure 2.5.
To add or collate two meshes together, one can simply use the + operator. Note that to collate two meshes, the vertices on the adjoining boundary must correspond (see figure 2.10).
1 mesh Th=Th1+Th2 ;
2 p l o t( Th ) ;
Script 2.7: Script to add two meshes to form a single mesh
Mesh adaptation A basic tool that can be used for mesh adaptation is the function splitmesh. This function divides each triangle by the value at the center of the triangle of a spatial function. The splitmesh function can be used to divide some or all triangles in a mesh. In fact, one can use a characteristic function to divide only the triangles in the support of the characteristic function. An example of the use of the previous function is given in script 2.8. The resulting is shown in figure 2.7.
2.2. FREEFEM++ AND ITS INTERPRETED LANGUAGE
Figure 2.5: Mesh generated with the truncation tool
Figure 2.6: Mesh of two collated rectangles
2.2. FREEFEM++ AND ITS INTERPRETED LANGUAGE
1 b o r d e r b ( t = 0 , 4 . ∗p i) {x= −.5∗c o s( t ) +c o s( − 0 . 5 ∗ t ) + 1 . 5 ;y=−s i n( − 0 . 5 ∗ t ) ; l a b e l= 1 ; }
2 mesh Th=b u i l d m e s h( b ( 5 0 ) ) ;
3 p l o t( Th , w a i t= 1 ,p s= " n o s p l i t m e s h . e p s " ) ;
4
5 mesh Sh= s p l i t m e s h ( Th , 1 +x∗x) ; / / split in int(1+x*x) triangles
6 p l o t( Sh , w a i t= 1 ,p s= " s p l i t m e s h . e p s " ) ; / / see figure 2.7 (middle)
7
8 mesh Qh= s p l i t m e s h ( Th , 5 ) ; / / all triangles are split in 5
9 p l o t( Qh , w a i t= 1 ,p s= " s p l i t a l l m e s h . e p s " ) ; / / see figure 2.7 (right) Script 2.8: Script to split mesh triangles with splitmesh.
Figure 2.7: Mesh adaptation with the splitmesh function. Initial mesh (left), mesh with triangles split in 1 + x ∗ x triangles (middle), and mesh with all triangle split in 5 triangles (right).
The best tool for mesh adaptation is the adaptmesh function. This tool is based on a variable metric/Delaunay automatic meshing algorithm. Mesh adaptation is very useful to generate a finer mesh where the solution of a problem varies sharply. It is also very useful to make the mesh coarser in given regions where the solution varies very slowly in order to reduce computational cost. The generated mesh can be adapted to a function or finite element function. The function adaptmesh takes various arguments, including the required precision err, the minimum hmin or maximum hmax edge size of the triangles, or the maximum number of vertices nbvx. See the full documentation for the other arguments. Mesh adaptation can be performed following the example given in script 2.9 and illustrated in figure 2.8.
1 b o r d e r b ( t = 0 , 4 . ∗p i) {x= −.5∗c o s( t ) +c o s( − 0 . 5 ∗ t ) + 1 . 5 ;y=−s i n( − 0 . 5 ∗ t ) ; l a b e l= 1 ; }
2 mesh Th=b u i l d m e s h( b ( 7 0 ) ) ;
2.2. FREEFEM++ AND ITS INTERPRETED LANGUAGE
3 f e s p a c e Vh ( Th ,P1) ;
4 f u n c f = 10∗exp( −50∗s q u a r e(x−1) / −50∗s q u a r e(y) ) ; / / f = 10 ∗ e(−50∗(x−1)2+y2)
5 Vh u= f ; / / u is the projection of f to Vh
6 Vh u0 ;
7 p l o t( Th , u ,w a i t= 1 ,p s= " a d a p t m e s h 0 . e p s " ) ;
8 f o r (i n t i t = 0 ; i t < 5 ; ++ i t ) {
9 Th = a d a p t m e s h( Th , u , e r r = 1 . e −3) ; / / mesh adaptation
10 u= f ;
11 p l o t( Th , u ,w a i t= 1 ,p s= " a d a p t m e s h . e p s " ) ; / / see figure 2.8.
12 }
Script 2.9: Script to adapt a mesh with the function adaptmesh.
Note that the command adaptmesh does not destroy the old mesh and all finite element functions using the old mesh remain unchanged. In the script 2.9, the command u=f redefines the variable u on the new mesh and, as a result, destroys the old mesh because uwas the only variable defined on it.
Figure 2.8: Mesh adaptation with the adaptmesh function. Initial mesh (left), mesh adapted after one iteration (middle), and mesh adapted after five iterations (right).
Movemesh Meshes can be translated, rotated and deformed by the function movemesh.
This tool is very useful to follow the deformation due to displacement or to deal with mov-ing boundary problems. If Th is the triangulation and d the displacement vector, then the displaced mesh can be obtained with the script 2.10 and illustrated in figure 2.9.
2.2. FREEFEM++ AND ITS INTERPRETED LANGUAGE
1 b o r d e r b ( t = 0 , 4 . ∗p i) {x= −.5∗c o s( t ) +c o s( − 0 . 5 ∗ t ) + 1 . 5 ;y=−s i n( − 0 . 5 ∗ t ) ; l a b e l= 1 ; }
2 mesh Th=b u i l d m e s h( b ( 6 0 ) ) ;
3 p l o t( Th , w a i t= 1 ,p s= " movemesh0 . e p s " ) ;
4
5 Th = movemesh( Th , [y,−x] ) ; / / rotation
6 p l o t( Th , w a i t= 1 ,p s= " m o v e m e s h r o t a t i o n . e p s " ) ;
7
8 Th = movemesh( Th , [x+c o s(x) ,y] ) ; / / deformation
9 p l o t( Th , w a i t= 1 ,p s= " m o v e m e s h d e f o r m a t i o n . e p s " ) ;
Script 2.10: Mesh manipulation with the movemesh function
Figure 2.9: Mesh manipulation with the movemesh function. Initial mesh (left), mesh rotation of π rad (middle), and mesh deformation with a displacement of cos(x) in the x-direction (right).