Programming project in TMA4220 – part 1

The programming project will be split into two parts. The first part is an introduction into the finite

element method and is designed to build up a solid code base for part 2, where you will actually

solve larger, real life problems. For the second part you will have to choose one of several tasks to

implement.

1 Gaussian quadrature

At the heart of every finite element code, lies the evaluation of an integral. This integral may be of

varying complexity depending on the problem at hand, and many of these integrals does not even

have a known analytical solution. Some integrals are possible to solve analytically, but of such

computational complexity that it is impractical to do so. As such, one often refers to numerical

integration schemes to do the core integration. One popular integration scheme is the Gaussian

quadrature.

In one dimension the gauss quadrature takes the form

Z 1

−1

g(z)dz ≈

X

Nq

q=1

ρqg(zq),

where Nq is the number of integration points, zq are the Gaussian quadrature points and ρq are the

associated Gaussian weights.

This extends to higher dimensions by

Z

Ωˆ

g(z)dz ≈

X

Nq

q=1

ρqg(zq),

and specifying the vector quadrature points zq as well as integrating over a suitable reference

domain Ωˆ (i.e. squares or triangles in 2D, tetrahedra or cubes in 3D).

a) 1D quadrature

Write a Python function quadrature1D(a, b, Nq, g) that returns a value I where the

variables are defined as:

I ∈ R value of the integral

a ∈ R integration start

b ∈ R integration end

Nq ∈ [1, 4] number of integration points

g : R → R function pointer∗

Verify that the function evaluates correctly by comparing with the analytical solution of the integral

Z 2

1

e

x

dx

Nq zq ρq

1-point-rule 0 2

2-point-rule −

p

p

1/3 1

1/3 1

−

p

3/5 5/9

3-point-rule 0 8/9

p

3/5 5/9

−

q

3+2√

6/5

7

18−

√

30

36

4-point-rule −

q

3−2

√

6/5

7

18+√

30

q

36

3−2

√

6/5

7

18+√

30

q

36

3+2√

6/5

7

18−

√

30

36

Table 1: 1D gauss quadrature

b) 2D quadrature

Using all numerical quadratures, it is important to first map the function to the reference domain. In

one dimension, this is the interval ζ ∈ [−1, 1]. In higher dimensions, we often map to barycentric

coordinates (or area coordinates as they are known in 2D). The gauss points are then given as

triplets in this coordinate system. The area coordinates are defined by

ζ1 =

A1

A

ζ2 =

A2

A

ζ3 =

A3

A

where A1, A2 and A3 are the area of the triangles depicted in figure 1 and A is the total area of the

triangle. Note that these do not form a linear independent basis as ζ1 + ζ2 + ζ3 = 1.

Nq (ζ1, ζ2, ζ3) ρ

1-point rule (1/3, 1/3, 1/3) 1

(1/2, 1/2, 0) 1/3

3-point rule (1/2, 0, 1/2) 1/3

(0, 1/2, 1/2) 1/3

(1/3, 1/3, 1/3) -9/16

4-point rule (3/5, 1/5, 1/5) 25/48

(1/5, 3/5, 1/5) 25/48

(1/5, 1/5, 3/5) 25/48

Table 2: 2D gauss quadrature

Write a Python function quadrature2D(p1, p2, p3, Nq, g) that returns a value I where

the variables are defined as:

Figure 1: Barycentric coordinates in two dimensions

I ∈ R value of the integral

p1 ∈ R

2 first corner point of the triangle

p2 ∈ R

2

second corner point of the triangle

p3 ∈ R

2

third corner point of the triangle

Nq ∈ {1, 3, 4} number of integration points

g : R

2 → R function pointer∗

Hint: An easy way of mapping barycentric coordinates ζ to physical coordinates x is by x =

ζ1p1 + ζ2p2 + ζ3p3

, where pi

, i = 1, 2, 3 is the corner points of the triangle.

Verify that the function evaluates correctly by comparing with the analytical solution of the integral

Z Z

Ω

log(x + y) dx dy,

where Ω is the triangle defined by the corner points (1, 0), (3, 1) and (3, 2).

(*) A function pointer in Python is a variable which represents a function instead of the usual

numerical values. In its simplest form it is declared as

f = lambda x: x**2 + 1

which would cause the variable f to contain a pointer to the function f(x) = x

2 + 1. The function

can then be evaluated by

y = f(4)

which should yield the result y = 17. A function may take in several arguments, i.e. f(x, y) =

x

2 + y

2 may be declared as

f = lambda x, y: x**2 + y**2

again the evaluation of the function is straightforward

y = f(2, 2)

Provided that the actual function body is capable of vector or matrix operations, then the input

arguments may be of vector or matrix form. The syntax remains unchanged by this.

2 Poisson in 2 dimensions

We are going to solve the two-dimensional Poisson problem, given by

∇2u(x, y) = −f(x, y) (1)

u(x, y)|r=1 = 0,

with f given in polar coordinates as

f(r, θ) = −8π cos

2πr2

+ 16π

2

r

2

sin

2πr2

on the domain Ω given by the unit disk, i.e. Ω =

(x, y) : x

2 + y

2 ≤ 1

.

a) Analytical solution

Verify that the following expression is in fact a solution to the problem (3)

u(x, y) = sin

2π(x

2 + y

2

)

. (2)

b) Weak formulation

Show that the problem can be rewritten as

a(u, v) = l(v), ∀v ∈ X.

with the bilinear functional a and the linear functional l given by

a(u, v) = Z Z

Ω

∇u · ∇v dx dy,

l(v) = Z Z

Ω

fv dx dy.

What is the definition of the space X?

c) Galerkin projection

Instead of searching for a solution u in the entire space X we are going to be looking for a

solution in a much smaller space Xh ⊂ X. Let Ω be discretized into M triangles such that our

computational domain is the union of all of these Ω = ∪M

k=1Kk. Each triangle Kk is then defined

by its three corner nodes xi

. For each of these nodes there corresponds one basis function. The

space Xh is then defined by

Xh = {v ∈ X : v|Kk ∈ P1(Kk), 1 ≤ k ≤ M}

for which the basis functions {ϕi}

n

i=1 satisfy

Xh = span{ϕi}

n

i=1 ϕj (xi) =

and δij is the Kronecker delta. By searching for a solution uh ∈ Xh, it is then possible to write

this as a weighted sum of the basis functions, i.e. uh =

Pn

i=1 u

i

hϕi(x, y).

Show that the problem ”Find uh ∈ Xh such that a(uh, v) = l(v) ∀v ∈ Xh” is equivalent to the

following problem

Find u such that

Au = f (3)

with

A = [Aij ] = [a(ϕi

, ϕj )]

u = [u

i

h

]

f = [fi

] = [l(ϕi)].

d) Implementation

We are now going to actually solve the system (3). First we are going to take a look at the triangulation {Kk}. From the webpage https://wiki.math.ntnu.no/tma4220/2019h/

project you may download the mesh generators.

The function getDisk is generating the unit disk Ω. Plot at least three meshes of different sizes

using the mesh generated from this function.

e) Stiffness matrix

Build the stiffness matrix A. Use the Gaussian quadrature from exercise 1 to do this.

The matrix A should now be singular. Verify this in your code and explain why this is the case.

f) Right hand side

Build the right hand side vector f in the same manner as A.

g) Boundary conditions

Implement the homogeneous dirichlet boundary conditions. Describe what method you used for

this and how you did it.

h) Verification

Solve the system (3) and verify that you are getting (approximately) the same result as the analytical solution (2).

3 Neumann boundary conditions

We are going to change to boundary conditions of our problem to

∇2u(x, y) = −f(x, y) (4)

u(x, y)|∂ΩD

= 0,

∂u(x, y)

∂n

∂ΩN

= g(x, y),

with the source term f and exact solution u given as above, and g as

g(r, θ) = 4πr cos

2πr2

. (5)

The dirichlet boundary condition is defined on ∂ΩD =

x

2 + y

2 = 1, y < 0

, and the neumann

boundary condition as ∂ΩN =

x

2 + y

2 = 1, y > 0

shown in figure 2.

Figure 2: Dirichlet and Neumann boundary conditions

a) Boundary condition

Verify that (2) satisfies (5) at the boundary.

b) Variational formulation

How does a(·, ·) and l(·) change with the introduction of Neumann boundary conditions

c) Gauss quadrature

The neumann boundary condition is given as an integral and should be evaluated using Gaussian

quadrature. Modify your quadrature methods from task 1 to solve line integrals in two dimensions,

i.e. the method signature in quadrature1D(a, b, Nq, g) should change to a ∈ R

2

and b

∈ R

2

.

d) Implementation

Change your code from task 2 to solve this new boundary value problem. How does your solution

in the interior compare to the one you got in task 2? How does your solution at the boundary

compare?

EasyDue™ 支持PayPal, AliPay, WechatPay, Taobao等各种付款方式!

**E-mail:** easydue@outlook.com **微信:**easydue

**EasyDue™是一个服务全球中国留学生的专业代写公司
专注提供稳定可靠的北美、澳洲、英国代写服务
专注提供CS、统计、金融、经济、数学等覆盖100+专业的作业代写服务**