Personal Website:
Scott Congreve

<< Back

Finite Element Methods 1 (NMNV405) — Practicals (Winter Semester 2024/2025)

Practicals

Tuesday 9:00 – 10:30, K7 Sokolovská 83 Karlín

Homework

There will be four homeworks during the course of the year. Obtaining credit for this course will involve obtaining at least 50% of the marks available from these homeworks.

1. Homework 1; Deadline: 19.11.2024; Solutions
2. Homework 2; Deadline: 03.12.2024; Solutions
3. Homework 3; Deadline: 17.12.2024; Solutions
4. Homework 4; Deadline: 07.01.2025; Solutions

Basis Functions

The following documents contains details of various finite elements: examples_of_finite_elements.pdf.

3D visualisations of various basis functions on triangles and rectangles can be found here.

The displayed basis function can be changed via the selection boxes in the top left corner. The bottom left corner allows switching the view between 3D (orthogonal and perspective) and 2D views. In 3D view the basis function can be rotated by dragging. In the bottom left corner the principal lattice (white lines) can be toggled on/off.

The following website also lists details about all the basis functions we have discussed (plus a lot more): defelement.com

Notes on Gaussian Quadrature

The following journal article has information on the choice of points and weights for Gaussian quadrature on a triangle:

Dunavant, D.A. (1985), High degree efficient symmetrical Gaussian quadrature rules for the triangle. Int. J. Numer. Meth. Engng., 21:1129-1148. https://doi.org/10.1002/nme.1620210612

Example FEM Code

An example MATLAB finite element code which solves the problem

$$-\Delta u = f,$$

on the unit square $[0,1]^2$ with forcing function

$$f=2\pi^2\sin\pi x\sin\pi y$$

and homogeneous Dirichlet boundary conditions is available in the following zip file:

• femcode.zip.

The main important files are:

fem.m

The main finite element code

basis_*.m and grad_basis_*.m

Evaluates the basis functions and gradients, repectively, for the various element types

quadrature_rect.m

Computes numerical quadrature points and weights for a rectangle from 1D Guass-Legendre quadrature points. The 1D points are dynamically computed by using Newton's method to find the roots of the relevant Legendre polynomial, where an initial guess for the ith point is given by $\cos \left(\frac{i-.25}{n+.5}\pi\right)$ for $n$ total points (this guess is good enough to ensure all points are found uniquely)

quadrature_tri.m

Looks up the numerical quadrature points and weights for a triangle from the tables in Dunavant, D.A. (1985)

To run the finite solver for the above problem on a $16\times 16$ mesh of squares with Lagrange polynomials of polynomial order 3 ($Q_3$) enter

fem(16,'rect',3,1)

and to run on a $16\times 16$ mesh of squares, where each square is divided into two triangles, with Lagrange polynomials of polynomial order 3 ($P_3$) enter

fem(16,'tri',3,1)

In general, the solver can be run on an $n\times n$ mesh of squares (or triangles) using different finite element types at a specific polynomial degree $p$ by running

fem(n,ele_type,p,1)

where ele_type is one of the following:

• 'tri' - For Lagrange polynomials of order $p$ on triangles; i.e. $P_p$
• 'rect' - For Lagrange polynomials of order $p$ on squares; i.e. $Q_p$
• 'reduced_tri' - For reduced 2-simplex polynomials of order 3 on triangles; i.e. $P_p'$ (p must be 3)
• 'reduced_quad' - For reduced rectangle polynomials of order $p$ on square; i.e. $Q_p'$ (p must be 2 or 3)
• 'cr' - For non-conforming Crouzeix-Raviart on triangles (p must be 1)
• 'rotated_bilinear' - For non-conforming rotated bilinear on squares (p must be 1)

Note the last argument to fem is optional. If it is specified as 1 then the resulting solution is displayed as a surface plot, otherwise the solution is not shown and just the values of the degrees of freedom are returned.

The above problem has a known analytical solution

$$u=\sin\pi x\sin\pi y$$ and, hence, we can compute the actual errors in the $L^2$-norm and $H^1$-seminorm; i.e., $$\| u-u_h \|_{0,\Omega} \qquad \text{and} \qquad |u-u|_{1,\Omega}$$

We know that

$$\|u-u_h\|_{1,\Omega} \leq C h^k|u|_{k,\Omega}=\overline{C}h^k,$$ and letting $\mathcal{E}_N = \|u-u_h\|_{1,\Omega}$ we have, asymptotically, that $$\begin{aligned} \mathcal{E}_N &\approx \overline{C} h^k \\ \log_{10} \mathcal{E}_N &\approx \log_{10}\overline{C} + k \log_{10} h \\ \log_{10} \mathcal{E}_N &\approx \widetilde{C} + k \log_{10} h \end{aligned}$$

Therefore, if we compute the solution for various values of $h=n^{-1}$, plot $\log_{10} \mathcal{E}_N$ against $\log_{10} \mathcal{h}_N$, and perform linear regression we will get an estimate for $k$ (Note, a similar result exists for the L2-error where we expect $k+1$ order of convergence).

To run a test on a sequence of meshes with sizes $h=2^{-1},...,(2^{6})^{-1}$, and perform the above linear regression, for a specific element type and polynomial degree simply run the following command in MATLAB:

fem_convergence(ele_type,p)

where ele_type and p are as for fem.

You should see approximately order $k=p$ convergence for all element types except reduced 2-simplex where you should see approximately order $k=p-1$ convergence.