Presentation + project: link, Matlab files to download: link

If you have no access to Matlab, try to install Octave instead (link here).

Download the package 'Matlab_FEM_course.zip' (link above) and unzip it. Enter

cd basics

to get to the folder, where intoductionary demos are located. Type

start_basics

to run the first example. The code generates the sequence of uniformly refined meshes up to a given maximal level (defined by levels). For illustration, levels 0 and 1 are displayed below.

It also compares two implementations (serial and vectorized) of areas computations. To keep your CPU busy, we can run computations on meshes with several million of elements (you can try it yourself). Can we compare the mentioned serial and vectorized implementation. Which is faster? If you are running Matlab and have Parallel Computing Toolbox enabled (together NVIDIA graphics card), you can even shift some computations from CPU to GPU.

Exercise 1: Can you modify the code to compute the total area of a (convex) polyhedral set given by point (0,0), (1,0), (1,1), (0.5,2), (-1,1)?

Triangular meshes shown above can be easily used eg. for visualizations of functions. Let us consider for instance the function y=sin(\pi x) defined over the L-shape domain above (here \pi denotes Ludolph's number). The, the command

show_nodal_scalar(sin(pi*coordinates(:,1)),elements,coordinates)

displays the shape of the function. Note that if variables 'elements' and 'coordinates' correspond to very coarse meshes, the resulting shape of the function does not look very smooth. In fact, it is never smooth whit our approximation, but for a reasonable fine resolution you will not notice it.

Exercise 2: How would you modify the command above to visualize the function y=sin(\pi x) sin(\pi y) over the rectangle [0,2] x [0,1] ?

Parts of this code were based on our package

'Fast assembly of stiffness and matrices in finite element method using nodal elements' (link here).

The package allows for fast assembly of matrices needed to solve certain classes of partial differential equations (PDEs) by using the finite elements method (FEM).

Recent versions of Matlab contain a PDE toolbox that computes several classes of PDEs of interest. A user does not need to go to details of FEM, but should be able to setup a model geometry, boundary conditions, various coefficients, create a triangular mesh and run the solver. We demonstrate the PDE package of a simple example. All commands showing the functionality of the PDE toolbox are collected in the file 'test_PDE_toolbox' below.

Type two commands

[p,e,t]=initmesh('lshapeg');

pdemesh(p,e,t)

to create and display an initial triangular mesh of an L-shape model (which is already define for us in MATLAB). Note that the mesh is defined by three objects: points p (similar to coordinates above), edges e (they were used but hidden up to now) and triangles t (similar to elemens above). Their exact structure can be found in Matlab documentation. If the obtain refinement is not sufficiently fine, we can refine it again by the command

[p,e,t] = refinemesh('lshapeg',p,e,t)

Both meshes are visualized below:

Note that the mesh is defined by three objects: points p (similar to coordinates above), edges e (they were used but hidden up to now) and triangles t (similar to elements above). Their exact structure can be found in Matlab documentation. Image having an elastic membrane attached to the frame with the shape given by a boundary of the L-shape geometry. You can compute the deformation of the membrane from a Poisson equation (containing famous Laplace operator, you learnt about in the vector analysis course). The computation can be run by the command

u=assempde('lshapeb',p,e,t,1,0,1);

and a displacement can be displayed by commands

figure;

pdeplot(p,e,t,'xydata',u,'zdata',u,'mesh','on');

axis equal

Displacements are displayed below: