Lesson 2

We consider the finite element discretization of the following problem

- \nabla \cdot (A \nabla u) + ( \mathbf{v} - \mathbf{w} )  \cdot \nabla u + b u &= f, ~~~~ \Omega
  \\
  u &= g, ~~~~ \Gamma_D
  \\
  A \nabla u \cdot \mathbf{n} &= k, ~~~~ \Gamma_N

which can be written in a weak formulation, for any test function \phi that vanishes on \Gamma_D,

\int_{\Omega} A \nabla u \cdot \nabla \phi ~d\Omega + \int_{\Omega} \left( \mathbf{v} \cdot \nabla u \right) \phi ~d\Omega + \int_{\Omega}  u \mathbf{w} \cdot \nabla \phi~d\Omega + \int_{\Omega} ( b + \nabla \cdot \mathbf{w} ) u \phi  ~d\Omega =
\\
\int_{\Omega} f \phi ~d\Omega + \int_{\Gamma_N}  A \nabla u \cdot \mathbf{n}  \phi ~d\Gamma_N + \int_{\Gamma_D}  u \mathbf{w} \cdot \mathbf{n}  \phi ~d\Gamma_D

It is assumed that the domain \Omega has a smooth boundary \partial \Omega such that

\Gamma_D \cup \Gamma_N &=\partial \Omega
\\
\Gamma_D \cap \Gamma_N &= \emptyset

We also assume the existence of a mapping \mathbf{F} such that \Omega=\mathbf{F}(\mathcal{P}) where \mathcal{P} is the unit line/square. \mathcal{P} is called a patch, logical domain or reference element.

for the moment, we will focus on the simple case where :math:`\mathbf{F}` is the identity mapping.

The weak formulation can be written as

a(u,\phi) = l(\phi)

where

a(u,\phi) = \int_{\Omega} \{ A \nabla u \cdot \nabla \phi + \left( \mathbf{v} \cdot \nabla u \right) \phi + u \mathbf{w} \cdot \nabla \phi + ( b + \nabla \cdot \mathbf{w} ) u \phi \} ~d\Omega

and

l(\phi) = \int_{\Omega} f \phi ~d\Omega + \int_{\Gamma_N}  k \phi ~d\Gamma_N + \int_{\Gamma_D}  g \mathbf{w} \cdot \mathbf{n} \phi ~d\Gamma_D

Let us denote \mathcal{D} the following bilinear differential operator

\mathcal{D}(u,\phi) = A \nabla u \cdot \nabla \phi + \left( \mathbf{v} \cdot \nabla u \right) \phi + u \mathbf{w} \cdot \nabla \phi + ( b + \nabla \cdot \mathbf{w} ) u \phi

Galerkin-Ritz approximation

Let us consider the following sets:

\mathcal{X} &= H^1(\Omega)
\\
\mathcal{S} &= \{ v/ v \in \mathcal{X} , v|_{\Gamma_D} = g \}
\\
\mathcal{V} &= \{ v/ v \in \mathcal{X} , v|_{\Gamma_D} = 0 \}

Let \mathcal{X}_h, \mathcal{V}_h be finite dimensional subspaces of \mathcal{X}, \mathcal{V}, where h is a parameter intended to tend to zero. Let \mathcal{S}_h be a finite dimensional approximation of \mathcal{S}. Remark that the last set is not necessarily a vector space.

The discrete Galerkin formulation writes:

Find u_h = u_h^0 + g_h  \in \mathcal{S}_h such that:

a(u_h,\phi_h) = l(\phi_h), \hspace*{1cm} \forall \phi_h \in \mathcal{V}_h

where here g_h, k are given functions, with g_h \in \mathcal{S}_h.

We suppose a rearrangement of the basis functions such that \varphi_i|_{\Gamma_D} = 0, \forall i \in \{ 1, .., n - n_{D} \}. This helps us to construct a basis of \mathcal{V}_h s.t :

\forall \phi_h \in \mathcal{V}_h , \phi_h = \sum_{i=1}^{n - n_D} \left[ \phi_h \right]^{i} \varphi_i

g^h \in \mathcal{S}_h is such that \left[ g_h \right]^{i} = 0,  \forall i \in \{ 1, .., n - n_{D} \} and so, g^h  =  \sum_{i=n - n_D + 1}^{n} \left[ g_h \right]^{i} \varphi_i Therefore, u_h \in  \mathcal{S}_h writes

u_h =  \sum_{i=1}^{n - n_D} \left[ u_h \right]^{i} \varphi_i + \sum_{i=n - n_D + 1}^{n} \left[ g_h \right]^{i} \varphi_i~.

Note

From now on, we will drop the subscript h.

therefore

\sum_{j=1}^{n - n_D} \left[ u \right]^{j} a( \varphi_i , \varphi_j ) =
l(\varphi_i) -   \sum_{m=n - n_D + 1}^{n} \left[ g \right]^{m} a( \varphi_i , \varphi_m ),~~~~ \forall i \in \{ 1, .., n - n_{D} \}

For i,j \in \{ 1, .., n - n_{D} \}, we denote \Sigma_{i,j} = a(\varphi_i, \varphi_j) and L_{i} = l(\varphi_i) - a(\varphi_i,g_h)

therefore we get:

\sum_{j=1}^{n - n_D}  \Sigma_{i,j} \left[ u \right]^{j} = L_i, ~~~~\forall i \in \{ 1, .., n - n_{D} \}

Finally, let us introduce the matrix:

\Sigma = (\Sigma_{i,j})_{1\leq i,j \leq n - n_{D}}

and the vectors :

L = (L_{i})^T_{1\leq i \leq n - n_{D}}

\left[ u \right] = ( \left[ u \right]^{i})^T_{1\leq i \leq n - n_{D}}

this leads to the linear system

\Sigma \left[ u \right] = L~.

Assembling process

We have seen that

\sum_{j=1}^{n - n_D}  \Sigma_{i,j} \left[ u \right]^{j} = L_i,~~~~ \forall i \in \{ 1, .., n - n_{D} \}

\Sigma_{i,j} &= a(\varphi_i, \varphi_j) =  \int_{\Omega} \mathcal{D}(\varphi_i, \varphi_j) ~ d\Omega
\\
&= \sum_{e \in \mathcal{E}} \int_{Q_e} \mathcal{D}(\varphi_i, \varphi_j) ~ d\Omega
\\
&=  \sum_{e \in \mathcal{E}} a_e(\varphi_i, \varphi_j)

because \{Q_e, e \in \mathcal{E} \} form a partition of the physical domain \Omega, where \mathcal{E} is the set of all elements

The next step is to compute each of the element contribution a_e

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 \{x_0=0, x_1=\frac{1}{4},x_2=\frac{1}{2},x_3=\frac{3}{4}, x_4=1 \} for a cubic polynomial degree.

_images/quadrature_1d.png

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

_images/quadrature_2d.png

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 \mathcal{V}_h.

Todo

to finish

  1. matrix.py defines the discretization matrix \Sigma.

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