FreeFEM Documentation on GitHub

stars - forks

Plotting in Matlab and Octave

Overview

In order to create a plot of a FreeFEM simulation in Matlab© or Octave two steps are necessary:

  • The mesh, the finite element space connectivity and the simulation data must be exported into files

  • The files must be imported into the Matlab / Octave workspace. Then the data can be visualized with the ffmatlib library

The steps are explained in more detail below using the example of a stripline capacitor.

Note

Finite element variables must be in P1 or P2. The simulation data can be 2D or 3D.

2D Problem

Consider a stripline capacitor problem which is also shown in Fig. 52. On the two boundaries (the electrodes) \(C_{A}\), \(C_{K}\) a Dirichlet condition and on the enclosure \(C_{B}\) a Neumann condition is set. The electrostatic potential \(u\) between the two electrodes is given by the Laplace equation

\[\Delta u(x,y) = 0\]

and the electrostatic field \(\mathbf{E}\) is calculated by

\[\mathbf{E} = -\nabla u\]
 1int CA=3, CK=4, CB=5;
 2real w2=1.0, h=0.4, d2=0.5;
 3
 4border bottomA(t=-w2,w2){ x=t; y=d2; label=CA;};
 5border rightA(t=d2,d2+h){ x=w2; y=t; label=CA;};
 6border topA(t=w2,-w2){ x=t; y=d2+h; label=CA;};
 7border leftA(t=d2+h,d2){ x=-w2; y=t; label=CA;};
 8
 9border bottomK(t=-w2,w2){ x=t; y=-d2-h; label=CK;};
10border rightK(t=-d2-h,-d2){ x=w2; y=t; label=CK;};
11border topK(t=w2,-w2){ x=t; y=-d2; label=CK;};
12border leftK(t=-d2,-d2-h){ x=-w2; y=t; label=CK;};
13
14border enclosure(t=0,2*pi){x=5*cos(t); y=5*sin(t); label=CB;}
15
16int n=15;
17mesh Th = buildmesh(enclosure(3*n)+
18             bottomA(-w2*n)+topA(-w2*n)+rightA(-h*n)+leftA(-h*n)+
19             bottomK(-w2*n)+topK(-w2*n)+rightK(-h*n)+leftK(-h*n));
20
21fespace Vh(Th,P1);
22
23Vh u,v;
24real u0=2.0;
25
26problem Laplace(u,v,solver=LU) =
27          int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v))
28        + on(CA,u=u0)+on(CK,u=0);
29
30real error=0.01;
31for (int i=0;i<1;i++){
32   Laplace;
33   Th=adaptmesh(Th,u,err=error);
34   error=error/2.0;
35}
36Laplace;
37
38Vh Ex, Ey;
39Ex = -dx(u);
40Ey = -dy(u);
41
42plot(u,[Ex,Ey],wait=true);

Exporting Data

The mesh is stored with the FreeFEM command savemesh(), while the connectivity of the finite element space and the simulation data are stored with the macro commands ffSaveVh() and ffSaveData(). These two commands are located in the ffmatlib.idp file which is included in the ffmatlib. Therefore, to export the stripline capacitor data the following statement sequence must be added to the FreeFEM code:

 1include "ffmatlib.idp"
 2
 3//Save mesh
 4savemesh(Th,"capacitor.msh");
 5//Save finite element space connectivity
 6ffSaveVh(Th,Vh,"capacitor_vh.txt");
 7//Save some scalar data
 8ffSaveData(u,"capacitor_potential.txt");
 9//Save a 2D vector field
10ffSaveData2(Ex,Ey,"capacitor_field.txt");

Importing Data

The mesh file can be loaded into the Matlab / Octave workspace using the ffreadmesh() command. A mesh file consists of three main sections:

  1. The mesh points as nodal coordinates

  2. A list of boundary edges including boundary labels

  3. List of triangles defining the mesh in terms of connectivity

The three data sections mentioned are returned in the variables p, b and t. The finite element space connectivity and the simulation data can be loaded using the ffreaddata() command. Therefore, to load the example data the following statement sequence must be executed in Matlab / Octave:

 1%Add ffmatlib to the search path
 2addpath('add here the link to the ffmatlib');
 3%Load the mesh
 4[p,b,t,nv,nbe,nt,labels]=ffreadmesh('capacitor.msh');
 5%Load the finite element space connectivity
 6vh=ffreaddata('capacitor_vh.txt');
 7%Load scalar data
 8u=ffreaddata('capacitor_potential.txt');
 9%Load 2D vector field data
10[Ex,Ey]=ffreaddata('capacitor_field.txt');

2D Plot Examples

ffpdeplot() is a plot solution for creating patch, contour, quiver, mesh, border, and region plots of 2D geometries. The basic syntax is:

1[handles,varargout] = ffpdeplot(p,b,t,varargin)

varargin specifies parameter name / value pairs to control the plot behaviour. A table showing all options can be found in the ffmatlib documentation. A small selection of possible plot commands is given as follows:

  • Plot of the boundary and the mesh:

1ffpdeplot(p,b,t,'Mesh','on','Boundary','on');
../_images/capacitor_boundary_mesh_500x400.png

Fig. 52 Boundary and Mesh

  • Patch plot (2D map or density plot) including mesh and boundary:

1ffpdeplot(p,b,t,'VhSeq',vh,'XYData',u,'Mesh','on','Boundary','on', ...
2          'XLim',[-2 2],'YLim',[-2 2]);
../_images/capacitor_patch_500x400.png

Fig. 53 Patch Plot with Mesh

  • 3D surf plot:

1ffpdeplot(p,b,t,'VhSeq',vh,'XYData',u,'ZStyle','continuous', ...
2          'Mesh','off');
3lighting gouraud;
4view([-47,24]);
5camlight('headlight');
../_images/capacitor_surf_500x400.png

Fig. 54 3D Surf Plot

  • Contour (isovalue) and quiver (vector field) plot:

1ffpdeplot(p,b,t,'VhSeq',vh,'XYData',u,'Mesh','off','Boundary','on', ...
2          'XLim',[-2 2],'YLim',[-2 2],'Contour','on','CColor','b', ...
3          'XYStyle','off', 'CGridParam',[150, 150],'ColorBar','off', ...
4          'FlowData',[Ex,Ey],'FGridParam',[24, 24]);
../_images/capacitor_contour_quiver_500x400.png

Fig. 55 Contour and Quiver Plot

Download run through example:

Matlab / Octave file

FreeFEM script

3D Plot Examples

3D problems are handled by the ffpdeplot3D() command, which works similarly to the ffpdeplot() command. In particular in three-dimensions cross sections of the solution can be created. The following example shows a cross-sectional problem of a three-dimensional parallel plate capacitor.

../_images/capacitor3d_slice_500x400.png

Fig. 56 Slice on a 3D Parallel Plate Capacitor

Download run through example:

Matlab / Octave file

FreeFEM script

References

Table of content