Lesson 2¶

We consider the finite element discretization of the following problem which can be written in a weak formulation, for any test function that vanishes on , It is assumed that the domain has a smooth boundary such that We also assume the existence of a mapping such that where is the unit line/square. is called a patch, logical domain or reference element. The weak formulation can be written as where and Let us denote the following bilinear differential operator Galerkin-Ritz approximation¶

Let us consider the following sets: Let be finite dimensional subspaces of , where is a parameter intended to tend to zero. Let be a finite dimensional approximation of . Remark that the last set is not necessarily a vector space.

The discrete Galerkin formulation writes:

Find such that: where here are given functions, with .

We suppose a rearrangement of the basis functions such that . This helps us to construct a basis of s.t :  is such that and so, Therefore, writes Note

From now on, we will drop the subscript .

therefore For , we denote and therefore we get: Finally, let us introduce the matrix: and the vectors :  this leads to the linear system Assembling process¶

We have seen that  because form a partition of the physical domain , where is the set of all elements

The next step is to compute each of the element contribution Element contribution¶

Element contribution relies on computing integrals over a given element. As our basis functions (and all others functions are smooth) we can use a quadrature formulae.

Note

The choice of a quadrature formulae depends on the basis functions we want to integrate. Generally, when their derivatives presents a discontinuity, a natural choice is the use of the Gauss-Legendre points.

Quadrature formulae in 1D¶

You can find, here the complete script that generates quadratures points on a given line.

The next plot shows the qudrature points formulae for a mesh for a cubic polynomial degree. The following example shows how to use the module quadratures

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 # -*- coding: UTF-8 -*- #! /usr/bin/python import numpy as np from quadratures import * import matplotlib.pyplot as plt import matplotlib matplotlib.rcParams.update({'font.size': 9}) qd = quadratures() # ... create the 1D mesh n = 3 x = np.linspace(0.,1., n+2) # ... # ... polynomial degree p = 2 # ... # ... generate the quadrature points # and their corresponding weights # on the whole mesh [x_leg,w_leg] = qd.generate(x, p, "legendre") [x_lob,w_lob] = qd.generate(x, p, "lobatto") [x_rdl,w_rdl] = qd.generate(x, p, "radau_left") [x_rdr,w_rdr] = qd.generate(x, p, "radau_right") # ... # ... x_leg = np.ravel(x_leg) x_lob = np.ravel(x_lob) x_rdl = np.ravel(x_rdl) x_rdr = np.ravel(x_rdr) y_lob = 0.*np.ones_like(x_lob) y_leg = 0.5*np.ones_like(x_leg) y_rdl = 1.*np.ones_like(x_rdl) y_rdr = 1.5*np.ones_like(x_rdr) plt.plot(x_leg, y_leg, 'ob', label="Gauss-Legendre") plt.plot(x_lob, y_lob, 'or', label="Gauss-Lobatto") plt.plot(x_rdl, y_rdl, 'og', label="Gauss-Radau-left") plt.plot(x_rdr, y_rdr, 'oc', label="Gauss-Radau-right") for xi in x: plt.axvline(x=xi, ymin=0., ymax = 1., linewidth=0.5, linestyle='--', color='k') plt.legend() plt.xlim(-0.1, 1.1) plt.ylim(-0.1, 2.1) plt.title("Gauss quadrature points on " + str(x) + " for polynomial degree of " + str(p) ) plt.savefig("ex4_quadrature_1d.png") # ...

Note the shape of the arrays x_leg and w_leg

>>> x_leg.shape
(4,3)
>>> w_leg.shape
(4,3)

In fact the first dimension of x_leg and w_leg refeers to the element, while the second one refeers to the quadrature point.

Quadrature formulae in 2D¶

The following script shows how to construct a 2D grid by calling two time the generate method and using the numpy-meshgrid function.

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 # -*- coding: UTF-8 -*- #! /usr/bin/python import numpy as np from quadratures import * import matplotlib.pyplot as plt import matplotlib matplotlib.rcParams.update({'font.size': 9}) qd = quadratures() # ... create the 1D meshes nx = 3 ; ny = 4 x = np.linspace(0.,1., nx+2) y = np.linspace(0.,1., ny+2) # ... # ... polynomial degree p = 2 # ... # ... generate the quadrature points # and their corresponding weights # on the whole mesh [x_leg,wx_leg] = qd.generate(x, p, "legendre") [y_leg,wy_leg] = qd.generate(y, p, "legendre") # ... # ... x_leg = np.ravel(x_leg) y_leg = np.ravel(y_leg) X,Y = np.meshgrid(x_leg, y_leg) plt.plot(X, Y, 'ob') for xi in x: plt.axvline(x=xi, ymin=0., ymax = 1., linewidth=0.5, linestyle='--', color='k') for yi in y: plt.axhline(y=yi, xmin=0., xmax = 1., linewidth=0.5, linestyle='--', color='k') plt.xlim(-0.1, 1.1) plt.ylim(-0.1, 1.1) plt.title("Gauss-Legendre quadrature points for polynomial degree of " + str(p) ) plt.savefig("ex4_quadrature_2d.png") # ...

The resulting grid is shown in the next figure Implementing the Finite Element Method¶

You can find the compressed project directory here. This project consists of the following files

1. quadratures.py includes different quadrature formulae, as described in the previous section.
2. grid.py defines 1D and 2D grids.
3. connectivity.py defines 1D and 2D connectivities depending on the number of elements, polynomial degrees and regularity.

Todo

Initialize the arrays ID, LM, IEN depending on the number of elements, polynomial degree and regularity. You can start by considering a periodic or homogeneous Dirichlet boundary condition

1. basis.py defines 1D and 2D basis functions. Bernstein polynomials are provided.

Todo

add B-Splines using Bezier extraction matrices. This will be adressed in the next lesson.

1. space.py defines the vectorial space .

Todo

to finish

1. matrix.py defines the discretization matrix .

Todo

to finish

1. field.py defines both the right-hand-side and the unknown.

Todo

to finish

1. pde.py this is the main class of our project. A PDE object contains the assembly procedure. You must provide the element_contribution as an argument.

Todo

to finish