# A beginning example for the finite element method


We consider the following partial differential equation

$$\left[ \begin{array}{ll} -\Delta u + u = 1,& \text{on } \Omega\\ u=0,& \text{on } \partial \Omega \end{array} \right.$$

on the domain

$$\Omega = [0, 1]^2$$.

We use a regular triangular mesh with 32 congruent right triangles.
We follow our normal procedure and arrive at the weak form

$$\intb{\nabla u \cdot \nabla v} + \intb{uv} = \intb{v}$$

We now have the partial differential equation written in terms of its bilinear form

$$\mathcal{A}(u,v) = \mathcal{F}(v)$$ where
$$\displaystyle u_h(x,y) = \sum_{j = 1}^{Nh} \eta_j \phi_j(x,y) \; \forall u_h \in V_h$$

and $$v_h = \phi_i(x,y)$$

Substituting we have

$$\displaystyle \mathcal{A}( \sum_{j = 1}^{Nh} \eta_j \phi_j(x,y), \phi_i(x,y))$$
$$\displaystyle \sum_{j = 1}^{Nh} \eta_j\mathcal{A}(\phi_j(x,y), \phi_i(x,y))$$
and
$$A_{ij} = \mathcal{A}(\phi_j(x,y), \phi_i(x,y))$$

This can be written discretely as $$A \eta = F$$

Where $$A$$ is made up of the entries $$A(\phi_i, \phi_j)$$. These are only non-zero where the basis functions overlap. The connectivity given as well as the defintion of our linear basis functions dictate that only a first order stencil will contribute to a row in $$A$$. That is, only nearest neighbors contribute.

The symmetry in our regular grid allows us to consider only a single triangle as all triangles are simply a rotation of this base. This triangle has coordinates (0,0), (0.25, 0), (0, 0.25). We know from the rules of integration that no matter the orientation, all integrals against other basis functions in this triangle will be identical in other locations.

Consider our bilinear form here

$$\mathcal{A}(\phi_1, \phi_2) = \intb{\nabla(\phi_1) \cdot \nabla(\phi_2)} + \intb{(\phi_1 \phi_2)}$$

where each $$\phi$$ is a linear basis function with a value of one at its respective corner of the triangle and zero at the other two corners with linear variation between. These can be written as

$$\phi_a(x, y) = a_0 + a_1x + a_2y$$
$$\phi_b(x, y) = b_0 + b_1x + b_2y$$
$$\phi_c(x, y) = c_0 + c_1x + c_2y$$

where $$\phi_a(x,y)$$ is defined for the point (0,0),  $$\phi_b$$ for (0.25, 0), and $$\phi_c$$ for (0, 0.25).The coefficients can be solved for each basis using the restriction set about how the value varies internally. For example, $$\phi_a$$‘s coefficients are given by

$$\left[ \begin{array}{ccc} 1 & 0.25 & 0\\ 1 & 0 & 0.25\\ 1 & 0 & 0\\ \end{array} \right] \left[ \begin{array}{c} a_0\\ a_1\\ a_2\\ \end{array} \right] = \left[ \begin{array}{c} 1 \\ 0 \\ 0 \\ \end{array} \right]$$

where the matrix is the values of the coordinates and the right hand side is the value of the basis function as we defined. In this case we get that

$$\phi_(x,y) = 1 – 4x – 4y$$
$$\phi_b(x,y) = 4x$$
$$\phi_c(x,y) = 4y$$
In general we get

$$\phi_a(x,y) = 1 – \phi_b(x,y) – \phi_c(x,y)$$

Over our given triangle we have a possibility of nine total evaluations of $$A(\phi_i, \phi_j)$$. They are $$A(\phi_a, \phi_b), A(\phi_b, \phi_a), A(\phi_a, \phi_c), A(\phi_c, \phi_a), A(\phi_b, \phi_c), A(\phi_c, \phi_b), A(\phi_a, \phi_a), A(\phi_b, \phi_b), A(\phi_c, \phi_c)$$. However, because of symmetry we only have six integrals to evaluate. These are shown below.

$$A(\phi_a, \phi_a) = \intb{\nabla \phi_a \cdot \nabla \phi_a} + \intb{\phi_a \phi_a}$$
$$= \itinta{a}{a} + \itintb{a}{a}$$
$$= \dint{[-4,-4]\cdot[-4,-4]} + \dint{(1-4x-4y)^2}$$
$$= 1 + \frac{1}{192}$$

$$A(\phi_b, \phi_b) = \intb{\nabla \phi_b \cdot \nabla \phi_b} + \intb{\phi_b \phi_b}$$
$$= \itinta{b}{b} + \itintb{b}{b}$$
$$= \dint{[4,0]\cdot[4,0]} + \dint{(4x)^2}$$
$$= \frac{1}{2} + \frac{1}{192}$$

$$A(\phi_c, \phi_c) = \intb{\nabla \phi_c \cdot \nabla \phi_c} + \intb{\phi_c \phi_c}$$
$$= \itinta{c}{c} + \itintb{c}{c}$$
$$= \dint{[0,4]\cdot[0,4]} + \dint{(4y)^2}$$
$$= \frac{1}{2} + \frac{1}{192}$$

$$A(\phi_a, \phi_b) = A(\phi_b, \phi_a) \intb{\nabla \phi_a \cdot \nabla \phi_b} + \intb{\phi_a \phi_b}$$
$$= \itinta{a}{b} + \itintb{a}{b}$$
$$= \dint{[-4,4]\cdot[4,0]} + \dint{(1-4x-4y)(4x)}$$
$$= -\frac{1}{2} + \frac{1}{384}$$

$$A(\phi_a, \phi_c) = A(\phi_c, \phi_a) \intb{\nabla \phi_a \cdot \nabla \phi_c} + \intb{\phi_a \phi_c}$$
$$= \itinta{a}{c} + \itintb{a}{c}$$
$$= \dint{[-4,4]\cdot[0,4]} + \dint{(1-4x-4y)(4y)}$$
$$= -\frac{1}{2} + \frac{1}{384}$$

$$A(\phi_b, \phi_c) = A(\phi_c, \phi_b) \intb{\nabla \phi_b \cdot \nabla \phi_c} + \intb{\phi_b \phi_c}$$
$$= \itinta{b}{c} + \itintb{b}{c}$$
$$= \dint{[4,0]\cdot[0,4]} + \dint{(4x)(4y)}$$
$$= \frac{1}{384}$$

From these computations it should be noted that any edge connected to the right angle results in the same integral. Also, each non-right-angle point against itself gives the same integral. This implies additional symmetry in the problem. It should not be surprising given the grid.

Now to determine each matrix entry we must sum the contributions from each triangle in the correct way. For example, looking at $$\phi_1$$
we notice that there are six triangles surrounding it. In two triangles, point 1 is a $$\phi_a$$ contribution in the other four it is either a $$\phi_b$$ or $$\phi_c$$  contribution. Therefore, the diagonal entry on row one looks like

$$A_{11} = 2A(\phi_a, \phi_a) + 2A(\phi_b, \phi_b) + 2A(\phi_c, \phi_c)$$

$$= 2\frac{193}{192} + 4\frac{97}{192} = \frac{129}{32}$$

and likewise for the other two basis functions in the local neighborhood of point 1.

$$A_{12} = A(\phi_a, \phi_b) + A(\phi_b, \phi_a)$$
$$= – 2 \frac{191}{384} = -\frac{191}{192}$$

and lastly

$$A_{14} = A(\phi_a, \phi_c) + A(\phi_c, \phi_a)$$
$$= -2 \frac{191}{384} = -\frac{191}{192}$$

The rest of the nine nodes follow this pattern due to grid symmetry. However, we did not have a diagonal edge in node 1. Let’s look at node 2.

$$A_{22} = A_{11}$$
$$A_{21} = A_{23} = A_{12}$$
$$A_{24} = A(\phi_b, \phi_c) + A(\phi_c, \phi_b)$$
$$= 2\frac{1}{384} = \frac{1}{192}$$
$$A_{25} = A_{14}$$

These patterns continue and greatly reduce the number of computations required to fill out the left hand side matrix. The last piece required is the right hand side of the system. The full expansion is shown below.

Due to the choice of the source term being unity, each of the right hand side pieces is simply the summation of the basis functions for each triangle connected to a given node. That is,

$$\mathcal{F}(v) = \intb{v}$$

We need to evaluate this for our three cases $$\phi_a$$, $$\phi_b$$, and $$\phi_c$$. These are

$$\intb{\phi_a} = \dint{(1-4x-4y)} = \frac{1}{96}$$
$$\intb{\phi_b} = \dint{(4x)} = \frac{1}{96}$$
$$\intb{\phi_c} = \dint{(4y)} = \frac{1}{96}$$

Since all of the contributions are the same and each node is connected to six triangles in our mesh,
we can right each right hand side vector entry as

$$F_i = 6\frac{1}{96} = \frac{1}{16}$$

These pieces give the following system. The full system as solved can be written as

$$A = \left[ \begin{array}{ccccccccc} \frac{129}{32}& -\frac{191}{192}& 0& -\frac{191}{192}& & & & & \\ -\frac{191}{192}& \frac{129}{32}& -\frac{191}{192}& \frac{1}{192}& -\frac{191}{192}& & & & \\ 0& -\frac{191}{192}& \frac{129}{32}& 0& \frac{1}{192}& -\frac{191}{192}& & & \\ -\frac{191}{192}& \frac{1}{192}& 0& \frac{129}{32}& -\frac{191}{192}& 0& -\frac{191}{192}& & \\ & -\frac{191}{192}& \frac{1}{192}& -\frac{191}{192}& \frac{129}{32}& -\frac{191}{192}& \frac{1}{192}& -\frac{191}{192}& \\ & & -\frac{191}{192}& 0& -\frac{191}{192}& \frac{129}{32}& 0& \frac{1}{192}& -\frac{191}{192}\\ & & & -\frac{191}{192}& \frac{1}{192}& 0& \frac{129}{32}& -\frac{191}{192}& 0\\ & & & & -\frac{191}{192}& \frac{1}{192}& -\frac{191}{192}& \frac{129}{32}& -\frac{191}{192}\\ & & & & & -\frac{191}{192}& 0& -\frac{191}{192}& \frac{129}{32}\\ \end{array} \right]$$
$$\eta = \left[ \begin{array}{c} 0.0414\\ 0.0524\\ 0.0413\\ 0.0524\\ 0.0671\\ 0.0524\\ 0.0413\\ 0.0524\\ 0.0414\\ \end{array} \right]$$
$$F = \left[ \begin{array}{c} \frac{1}{16}\\ \frac{1}{16}\\ \frac{1}{16}\\ \frac{1}{16}\\ \frac{1}{16}\\ \frac{1}{16}\\ \frac{1}{16}\\ \frac{1}{16}\\ \frac{1}{16} \end{array} \right]$$

The plotted solution is shown in figure 1. As one might expect from the system of PDEs given. The $$\eta$$ values show symmetry about the line $$y=x$$.

Figure 1: Surface plot $$\eta$$ values

If the linear basis functions are swapped for quadratic ones we arrive at a much more complicated formulation. These basis functions are then

$$\phi_a(x,y) = a_0 + a_1x + a_2y + a_3xy + a_4x^2 + a_5y^2$$
$$\phi_b(x,y) = b_0 + b_1x + b_2y + b_3xy + b_4x^2 + b_5y^2$$
$$\phi_c(x,y) = c_0 + c_1x + c_2y + c_3xy + c_4x^2 + c_5y^2$$
$$\phi_d(x,y) = d_0 + d_1x + d_2y + d_3xy + d_4x^2 + d_5y^2$$
$$\phi_e(x,y) = e_0 + e_1x + e_2y + e_3xy + e_4x^2 + e_5y^2$$
$$\phi_f(x,y) = f_0 + f_1x + f_2y + f_3xy + f_4x^2 + f_5y^2$$

Also, the number of control point (quadrature points) in each triangle goes up from three to six. This is consistent with the six unknowns in our above basis functions. Thus, the problem becomes much more substantial computationally.

2011-11-17