428 The finite element method for partial differential equations
This equation is defined by a matrixkwhich is called thelocal stiffness matrixfor
the triangle under consideration. Introducing the vectors
b=
ybc
yca
yab
(13.19)
and
c=
xcb
xac
xbc
(13.20)
we see, after some calculation, that the stiffness matrixkcan be evaluated as
k=bbT+ccT, (13.21)which leads to the result
[k]=1
4 A
y^2 bc+xcb^2 ybcyca+xcbxac ybcyab+xcbxba
ybcyca+xcbxac y^2 ca+x^2 ac ycayab+xacxba
ybcyab+xcbxba ycayab+xacxba y^2 ab+xba^2
. (13.22)
We not only need the matrix representing the Laplacian operator, but we must also
evaluate the integral containing the source termf(r)inEq. (13.4). The continuous
functionf(r)is approximated by a piecewise linear function ̃fon the triangles –
just like the solutionφ(r). For a particular triangle with verticesa,b, andc,wehave
̃f(r)=faξa+fbξb+fcξc. (13.23)We must multiply this by the linear approximation forφ,Eq. (13.12), and then
integrate over the element, using(13.13). The result must then be differentiated
with respect toφa,φb,φc, which results in a vector element
ra= 2 A(
fa
12+
fb
24+
fc
24)
, (13.24)
and similar forrbandrc.
The matrix–vector multiplication can be carried out as a loop over all triangular
elements where for each element the stiffness matrix is applied to the three ver-
tices of that triangle. Note that the stiffness matrix should always act on theold
vector containing the field valuesφaand that the result should be added to the
new vector (which initially is set to zero; see the next section). If we have a matrix–
vector multiplication and a right hand side of the form(13.15), we can apply the
conjugate gradients method to solve the matrix equation.
We have overlooked one aspect of the problem: if a triangle contains vertices on
the boundary where the value of the solution is given (Dirichlet boundary condition),
the corresponding values ofφ(r)are known and therefore not included in the vector.