• Author / Uploaded
• Harshan Arumugam

Citation preview

Computational PDE’s II Tutorial Prof. Carsten Carstensen Humboldt-Universit¨at zu Berlin, Germany Hella Rabus

Tutorial 1

Basics to finite element methods

Goals +

Download and run the programs FEM15.m and FEM50.m. Understand the implementation and data structures of FEM15.m. + Learn to use Octave to implement a FE method. + Learn how to solve and display examples for different domains. + Learn how to compute and display basis functions.

1.1

Continuous formulations

The weak formulation of a PDE reads: Given a bounded linear functional b ∈ V ∗ in some Hilbert space V with scalar product a, seek its Riesz representation u, namely u ∈ V with a(u, v) = b(v) for all v ∈ V. The Laplace Problem In this tutorial, we will demonstrate the procedure of solving a PDE with FEM for the following model problem. Strong formulation Find u ∈ H 2 (Ω) such that −4u = f

in Ω,

u = 0 on ∂Ω.

Weak formulation Find u ∈ V := H01 (Ω) such that a(u, v) = b(v) for all v ∈ H01 (Ω) R R with a(u, v) = Ω ∇u · ∇v dx and b(v) = Ω f v dx.

1.2

Approximation in finite-dimensional subspace

Let V` = span {ϕ1 , ϕ2 , . . . , ϕm } be a finite-dimensional subspace of the Hilbert space V . Thus V` is also a Hilbert space. Let {ϕk }k∈K , with K = {1, . . . , m}, denote a set of basis functions of V` . Thus a|V` ×V` : V` × V` −→ is bilinear and V -elliptic, and b|V` ∈ V`∗ . The discrete solution u` ∈ V` satisfies

R

a(u` , v` ) = b(v` ) for all v` ∈ V` .

(1)

This reads a(u` , ϕ` ) = b(ϕ` ) for all k ∈ K and is equivalent to X X u= xj ϕj and b(ϕk ) = xj a(ϕj , ϕk ) for all k ∈ K, j∈K

AIMS

j∈K

2007/12/10-2007/12/15

1

Computational PDE’s II Tutorial Prof. Carsten Carstensen Humboldt-Universit¨at zu Berlin, Germany Hella Rabus    i.e. A  

x1 x2 .. . xm





    =  

b(ϕ1 ) b(ϕ2 ) .. .

    

with A := (Ajk ) = (a(ϕj , ϕk ))j,k=1,··· ,m .

b(ϕm )

Since a(·, ·) is bilinear, bounded, and V -elliptic, there exists a unique discrete solution u` ∈ V` satisfying (1), i.e. the stiffness matrix A is positive definite and regular. For Laplace problem V` = S01 (T ) ⊂ H01 (Ω) is the set of piecewise linear functions on the triangulation T` of Ω, being zero at the boundary. Let NΩ denote the set of all inner nodes of the triangulation T` . The nodal basis functions {ϕz }z∈NΩ are given by ( 1 for p = z, ϕz (p) = for all p ∈ N , 0 else, such that ϕz is an affine function on each element T ∈ T` .

1.3

Stiffness matrix for the Laplace problem

The entries of the stiffness matrix can be determined by Z Z X ∇ϕj · ∇ϕk dx, Ajk = a(ϕj , ϕk ) = ∇ϕj · ∇ϕk dx = Ω

T ⊆ ωzj ∪ ωzk

T

(2)

where ωzj denotes the patch of the node zj , more precisely ωzj :=

[

T.

T ∈Tj , ϕ zj | T ≡ /0

Choosing V` = S01 (T ) as above, there are three nodal basis functions on each triangle, which are not identically zero and their gradients are constant here. For fixed T ∈ T , let {z1 , z2 , z3 } := T ∩ N denote the nodes of T . The gradients can be computed by       0 0 1 1 1 ∇ϕz1 Q :=  ∇ϕz2  = P −1  1 0  where P :=  xz1 xz2 xz3  . y z1 y z2 y z3 ∇ϕz3 0 1

AIMS

2007/12/10-2007/12/15

2

Computational PDE’s II Tutorial Prof. Carsten Carstensen Humboldt-Universit¨at zu Berlin, Germany Hella Rabus Thus, the local stiffness matrix AT for the triangle T = conv {z1 , z2 , z3 } is given by   R R R ∇ϕ · ∇ϕ dx ∇ϕ · ∇ϕ dx ∇ϕ · ∇ϕ dx z z z z z z 1 1 1 2 1 3 T T T   .. .. ... AT =   . . R R ∇ϕz3 · ∇ϕz1 dx ··· ∇ϕz3 · ∇ϕz3 dx T T 1 =|T | QQT = det(P ) QQT . 2 According to (2), the stiffness matrix can be computed by summing the corresponding entries of the local stiffness matrices. This can be done in a loop over the elements of the triangulation T as in FEM15.

1.4

Computation of right-hand side (RHS)

We need to compute the values of b(ϕj ), j = 1, . . . , m Z X Z b(ϕj ) = f ϕj dx = f ϕj dx. Ω

T ∈T` ϕ j |T ≡ /0

T

In the Laplace problem the basis functions are defined above and f ≡ 1, thus Z Z f ϕj dx = ϕj dx = |T |/3 = det(P )/6. T

T

The RHS can be computed analogously by summing up the corresponding local values. See FEM15 line 12.

1.5

Boundary conditions

Boundary conditions are constraints to u or its derivative Du on the boundary ∂Ω of Ω. Conditions for u on ∂Ω and for Du on ∂Ω are called Dirichlet boundary conditions and Neumann boundary conditions, respectively.

1.6

Example

The implementation of boundary conditions will be explained with respect to the following example. −4u = f u = uD ∂u =g ∂ν AIMS

in Ω, on ΓD ,

(3) (4)

on ΓN = ∂Ω\ΓD ,

(5)

2007/12/10-2007/12/15

3

Computational PDE’s II Tutorial Prof. Carsten Carstensen Humboldt-Universit¨at zu Berlin, Germany Hella Rabus where (4) and (5) represent the boundary conditions. To guarantee the existence and uniqueness of a solution (in a weak sense), ΓD has to be a nonempty, relatively open set compared to ∂Ω.

1.7

Dirichlet boundary conditions

The realisation of Dirichlet boundary conditions is done with a partition of the unknown solution u. Assuming there exists a u˜D ∈ H 1 (Ω) which realises the Dirichlet boundary condition γ(˜ uD )|ΓD = uD , where γ : H 1 (Ω) −→ L2 (∂Ω) is the trace operator, then  1 u = w + u˜D for w ∈ HD (Ω) = v ∈ H 1 (Ω) with γ(v)|ΓD = 0 . (6) The choice of the function space H 1 (Ω) stems from the theory of weak solutions.

1.8

Weak formulation with Neumann boundary conditions

Multiplying (3) by an arbitrary test function v ∈ C ∞ (Ω), integrating over the domain Ω, and then integrating by parts, leads to Z Z Z ∂u ∇u · ∇v dx = f v dx + v dx. Ω Ω ∂Ω ∂ν Hence, Z

Z

Z

∂u ∇w · ∇v dx = f v dx + v dx − Ω Ω ∂Ω ∂ν

Z ∇˜ uD · ∇v dx

(7)

Ω

1 for all v ∈ C ∞ (Ω). Since the unknown is w ∈ HD (Ω), such a solution w is uniquely defined if (7) is satisfied for all test functions with zero boundary values on ΓD . Hence Z Z Z Z ∂u ∞ dx− ∇˜ uD ·∇v dx for all v ∈ CD (Ω). ∇w·∇v dx = f v dx+ v ∂ν Ω ∂Ω Ω Ω

Since the test function v vanishes on the Dirichlet part of the boundary, and the gradient of u is known in direction of the outer unit normal through (5), 1 we obtain the weak formulation of (3): Find w ∈ HD (Ω) such that Z Z Z Z 1 ∇w · ∇v dx = f v dx + vg dx − ∇˜ uD · ∇v dx for all v ∈ HD (Ω), Ω

Ω

ΓN

Ω

(8) AIMS

2007/12/10-2007/12/15

4

Computational PDE’s II Tutorial Prof. Carsten Carstensen Humboldt-Universit¨at zu Berlin, Germany Hella Rabus 1 i.e. find w ∈ V = HD (Ω) such that a(w, v) = b(v) for all v ∈ V , where Z a(w, v) := ∇w · ∇v dx Ω Z Z Z b(v) := f v dx + vg dx − ∇˜ uD · ∇v dx. Ω

ΓN

Ω

Note that a(·, ·) is a bounded, V -elliptic bilinear form on V and b ∈ V ∗ .

1.9

Boundary conditions in FEM

As before, the corresponding formulation for an approximation of w in a 1 finite-dimensional subspace V` reads: Find w` ∈ V` = SD (T ) such that a(w` , v` ) = b(v` ) for all v` ∈ V` , where Z a(w` , v` ) = ∇w` · ∇v` dx, ZΩ Z Z b(v` ) = f v` dx + v` g dx, − ∇ˆ uD · ∇v` dx. Ω

ΓN

Ω

For simplicity we assume the Dirichlet boundary values to be piecewise affine on the outer edges. In this case uˆD = u˜D ∈ S 1 (T ). In general, the quality of uˆD ∈ S 1 (T ) realising the Dirichlet boundary condition uˆD |ΓD ≈ uD has an effect on the quality of u` = w` + uˆD approximating u.

1.10

Linear system of equations

As before we obtain a linear system of equations X X w= xj ϕj , u˜D = xj ϕj , zj ∈ND

j∈NΩ

   A  

x1 x2 .. . xm





b(ϕ1 )   b(ϕ2 )    = ..   . b(ϕm )

    

with A = (Ajk ) = (a (ϕj , ϕk ))1≤j,k≤m .

Here the nodes 1, 2, . . . , m are the inner nodes, ϕj are nodal basis functions and ND the set of nodes at the Dirichlet boundary. The coefficients {xj }zj ∈ND have to be determined with respect to the boundary values.

AIMS

2007/12/10-2007/12/15

5

Computational PDE’s II Tutorial Prof. Carsten Carstensen Humboldt-Universit¨at zu Berlin, Germany Hella Rabus

Exercises E 1.1

Getting to know Octave

If you are short in experience with the programming languages Octave and Matlab, go through the short introduction linked on the website, start Octave on your computer and use the short introduction about Octave for getting familiar with that language.

E 1.2

Start first example

Download archive FEM15.tar.gz to a new directory called CPDE2 and extract the archive by using the command tar -xvzf FEM15.tar.gz on your terminal. Change to the new directory FEM15, start octave on your computer by typing octave e.g. in the embedded terminal of gedit. Execute the script FEM15example.m by typing FEM15example on the Octave prompt.

E 1.3

Understand the program FEM15.m

Open the file FEM15.m with your favourite editor and have a look at the code. Make sure, that you understand each line of FEM15.m.

P3 Note the special structure of the data: Each vertex corresponds to one row of coordinates in c4n. Each triangle corresponds to one row of three vertices in n4e, numbered counterclockwise, starting with the vertices of the reference edge.

E3

@



P1

@ E2 @ @ @

P2

E1

To understand line 10 you may execute the following lines on the command line of Octave. A=[0,1,2;3,4,5] A(2,1)=10 A=zeros(4,4) A([4,2],[3,1])=[5,6;7,8] Observe how the last example acts on selected entries of A. AIMS

2007/12/10-2007/12/15

6

Computational PDE’s II Tutorial Prof. Carsten Carstensen Humboldt-Universit¨at zu Berlin, Germany Hella Rabus

E 1.4

Create triangulation

Modify the program FEM15example.m such that the solution of the problem −4u = 1 in Ω, u = 0 on ∂Ω with Ω as shown in Figure 1 will be computed.

1

0.5

E 1.5

P1 basis function

Write an m-file program that plots a P1 -basis function on a given grid and for a given node using triplot3 (see help triplot3 for usage).

E 1.6

50 lines of FEM code

0

−0.5

−1

−1

−0.5

0

0.5

1

Figure 1: L-shaped domain

1. Download the archive FEM50.tar.gz to your CPDE2 directory and extract it. 2. Run fem_50 and examine the output. 3. Create data for an initial grid on an L-shaped domain (cf. Figure 1), which can be used by fem_50. Then, solve the model problem 2 on that grid using fem_50. For better approximations use red_refine34.m to refine your mesh. Be aware of the correct structure of the geometric data, explained in exercise E 1.3.

AIMS

2007/12/10-2007/12/15

7