Finite Element Method(FEM) for Two Dimensional Laplace
Equation with Dirichlet Boundary Conditions

April 9, 2007

1

Variational Formulation of the Laplace Equation

The problem is to solve the Laplace equation

rPu = 0

(1)

in domain subject to Dirichlet boundary conditions on @
. We know from our study of the uniqueness of the solution of the Laplace equation that nding the solution is equivalent to nding u that minimizes
1 Z jjrujjPd
W=
(2)
2
subject to the same boundary conditions. Here the dierential d denotes the volume dierential and stands for dxdy for a plane region. W has interpretations such as stored energy or dissipated power in various problems.

2

Meshing

First we approximate the boundary of by polygons. Then can be divided into small triangles called triangular elements. There is a great deal of exibility in this division process. The term meshing is used for this division. For the resulting FEM matrices to be well-conditioned it is important that the triangles produced by meshing should not have angles which are too small.

At the end of the meshing process the following quantities are created.
 Nv : number of vertices or nodes.
 Nv ¢ 2 array of real numbers holding the x and y coordinates of the vertices.
 Ne: number of triangular elements.
1








Ne ¢ 3 array of integers holding the vertices of the triangular elements.

Nvf : Number of vertices on which the u values are not specied or free.
Nvf ¢ 1 array of integers holding free vertex indices.

Nvp : Number of vertices on which the u values are specied or prescribed.
Nvp ¢ 1 array of integers holding prescribed vertex indices.
Nvp ¢ 1 array of real numbers holding prescribed u values.

Data structures holding adjacency information for the vertices, edges, and the triangular elements are also generated by sophisticated FEM meshing subroutines.

3

Planar Approximation over a Triangle

£ ££ u = a + bx + cy d‚ ££ U d $$£ $$$ U
$
U
2
: (x ; y ) : (x ; y ) r 2 2
2r 2 ¢ ¢ rrr ¢ r¢ : (x ; y )
Let us consider a triangular element whose node numbers are ,  , and . The node coordinates are (x; y), (x ; y ), (x ; y ). On this triangle u(x; y) is assumed to have a planar variation: h u(x; y) = a + bx + cy = 1 x y

2 i6 4

a b c

3
7
5

(3)

First a, b, c are to be determined in terms of the values of u(x; y) at the nodes U, U , U , using the three equations u(xj ; yj ) = Uj ; j P f; ; g
(4)
In matrix form the equations are
2
32 3 2
3
1 x y a U
6
76 7 6
7
(5)
4 1 x y  5 4 b 5 = 4 U  5
1 x y c U
So
23
2
3 I 2
3
h i6 a 7 h i 6 1 x y 7 6 U 7 u(x; y) = 1 x y 4 b 5 = 1 x y 4 1 x y 5 4 U 5 c 1 x y
U
2
32
3 h i 1 6 x y x y x y x y x y x y 7 6 U 7 y y y y 5 4 U 5
= 1 x y 2A 4 y  y x x x x x x
U
(6)
= U (x; y) + U  (x; y) + U (x; y)

where and 2A = x y x y + x y xy + xy x y

= 21 [x y x y + (y y )x + (x x )y]
A
etc. 2A is twice the area of the triangle  . The functions j are interpolatory in nature. j; k P f; ; g j (xk ; yk ) = jk ;

 : (x ; y )

 (x; y )

(7)
(8)
(9)

2
22
: (x ; y )
¨¢
22222 ¨¨¢
22 ¨¨ ¢
2
¨
2r222
r
¢
rrr Pee
¢
rrr ee ¢¢ rrr ¢ ee r¢  : (x ; y )

The functions are also known as natural coordinates, simplex coordinates or areal coordinates.
Let P be the point (x; y). Then
Area of triangle P 
(10)
 (x; y ) =
Area of triangle 
Exercises:
 Show that (x; y) +  (x; y) + (x; y) = 1.
P
P
P
 At which point is (x; y) +  (x; y) + (x; y) minimum? What is the minimum value?
Since the original domain is now approximated as a union of small triangular elements, the total
W corresponding to stored energy or power dissipation can be expressed as a sum of element W 's.
For the  element
Z
(e) = 1
P
(11)
W
2 @
A jjrujj d
Now
ru = Ur  + U r  + U r
(12)
It should be noted that due to the planar variation of u(x; y) r (x; y) = 21 [(y y )^ + (x x )^ ] x y
(13)
A etc. are constant over the entire element. So
Z
1 X X U S (e)U
(e) = 1
W
jjrujjPd = 2
(14)
j jk k
2 @
A
j a;; ka;; where (e)
Sjk =

Z

@ A

(r j ) ¡ (r k )d; j; k P f; ; g

(15)

Let

2

U (e) = 6
4

Then

U
U
U

3
7
5

(16)

1T
(17)
2
S (e) is called the element Dirichlet matrix. Show that,
(e) 1
(18)
S = [(y y )(y y ) + (x x )(x x )]
2A
etc. The Dirichlet matrix is symmetric, positive semi-denite. The expression for element energy is a quadratic form in the nodal values U, U , U . Since W is expressed as a sum of all the element energies, it follows that W is given by a positive semi-denite quadratic form in the nodal values
UI , UP , . . . UN .
1
W = U T SU
(19)
2 where 2
3
W (e) = U (e) S (e) U (e)

U

6
6
=6
6
4

UI
UP

...

UN

7
7
7
7
5

(20)

and the Dirichlet matrix S has contributions from all element S (e)'s. In the program developed
(e)
here, S is initialized to zeros and as one loops over all elements any S contribution is added to
S . This is in contrast to other FEM programs which actually generate 3 ¢ 3 element matrices and then combine all of them during element assembly.

4

Minimization of

W

P of the N U values are prescribed. N P are free to vary. For simplicity of presentation we assume

that the prescribed values come after the free values. (MATLAB and modern matrix libraries allow sub-matrices to be selected based on arbitrary index sets. So there is no need to actually number the prescribed nodes after the free nodes.) Then the matrices U and S may be partitioned as follows:
"

U = Uf
Up
"

#

(21)
#

S = Sff Sfp
(22)
Spf Spp
T
Note that due to the symmetry of S , Sfp = Spf . In terms of the partitioned matrices W may be

expressed as

1 U T SU
2
"
#"
#
1 h U T U T i Sff Sfp Uf
=2 f p
Spf Spp
Up
1
= 1 UfT Sff Uf + 2 UpT Spf Uf + 1 UfT SfpUp + 1 UpT SppUp
2
2
2

W=

(23)
(24)
(25)

But since UpT Spf Uf is a scalar it equals its own transpose.
T
T
T
Up Spf Uf = (Up Spf Uf )T = UfT Spf Up = UfT Sfp Up

(26)
T
The last equality is due to the fact that Sfp = Spf . So by combining the two equal middle terms of
(25) we get
1T
1
(27)
W = UfT Sff Uf + UfT Sfp Up + Up Spp Up
2
2
Now this W is to be minimized with respect to Uf . At the minimum the gradient of W with respect to Uf equals zero. rUf W = Sff Uf + SfpUp = 0
(28)
So the solution is Uf = SffI Sfp Up
(29)
I
Once Uf is known, U is known, and one computes W = P U T SU . In physical problems, quantities related to W are usually important.
In eigenvalue problems like the Helmholtz equation, and loading problems like the Poisson
R
R
IR
equation, not only P  jjrujjPd , but also I  uPd is of importance. Just as I  jjrujjPd is
P
RPP expressed as I U T SU , the quantity I  u d can be expressed as I U T T U , where the T matrix is
P
P
P
called the metric matrix. Like the Dirichlet matrix S , T can be constructed from individual element contributions. 1 X X U T (e)U
1Z
uP d =
(30)
2 ¡@
A
2 ja;; ka;; j jk k where (e) = Z
Tjk
(31) j (x; y ) k (x; y )d; j; k P f; ; g
¡@
A
(e) R
First we evaluate Tjj = ¡@
A jP(x; y)d . We consider a dierential strip parallel to the side

opposite node j in the triangle  If j changes by d j on the strip, its width is hd j and the length of the strip is (1 j )b, where b is the length of the side opposite to node j and h is the height of the perpendicular from node j to the opposite side. So
Z
(e) Z
P d = I P (1 j )bhd j = bh = A
(32)
Tjj =
12 6
¡@
A j
Hj
(e) (e) where A is the area of triangle  . Thus the diagonal terms of the element T matrix, T , T ,
(e)
and T

, are equal to A=6. What about the o-diagonal terms?
(e) (e) (e)
T + T + T = since + 

+

Z

¡@ A

P+

= 1. So
(e) + T (e) + T = Z
(e)
T

+

d =

Z

¡@ A

 (  +  + )d =

Z

¡@ A

 d

(33)

(e)
But T = A=6. So
Similarly we also have

¡@ A

 d

(e) (e)
T + T =

=

Z

I
H

 (1

A

)bhd  = bh = A
63

A=A
366

(e) (e)
T + T =

A

6

(34)
(35)
(36)

(e) (e)
T
 + T =

A

(37)
6
(e)
(e)
But T = T etc. by the symmetry of the element T matrix. So adding the above three equations and dividing by 2 we see that
(e) (e) (e) A
(38)
T + T + T =
4
It it then seen that each o-diagonal element of the element T matrix is equal to A=12.

5

Example Code

The example code available on the web site solves two problems.
 Solution of Laplace equation on a rectangle: On three sides of the rectangle u = 0, while on the other side u = 1.
 Capacitance of a cable with elliptic cross sections for both the conductors.
Please download the code and run these examples.