### 4.16 How to solve Poisson PDE in 2D using ﬁnite elements methods using triangle element?

4.16.1 case 1
4.16.2 Program output
4.16.3 case 2

Solve $$\nabla ^{2}u=f\left ( x,y\right )$$ on square using FEM. Assume $$u=0$$ on boundaries and solve using $$f\left ( x,y\right ) =xy$$ and also $$f\left ( x,y\right ) =1$$ and $$f\left ( x,y\right ) =3\left ( \cos (4 \pi x)+\sin (3 \pi y)\right )$$

Use Galerkin method and weak form, Using triangle element.

Solution

Using as an example with 9 elements to illustrate the method. The program below can be called with diﬀerent number of elements.

Looking at one element

The trial function is

\begin {equation} \tilde {u}=c_{1}+c_{2}x+c_{3}y \tag {1} \end {equation}

Evaluating the trial function (1) at each corner node of the above element gives

\begin {align*} \tilde {u}_{1} & =c_{1}+c_{2}x_{1}+c_{3}y_{1}\\ \tilde {u}_{2} & =c_{1}+c_{2}x_{2}+c_{3}y_{2}\\ \tilde {u}_{3} & =c_{1}+c_{2}x_{3}+c_{3}y_{3} \end {align*}

Hence

$\begin {pmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {pmatrix} =\begin {pmatrix} 1 & x_{1} & y_{1}\\ 1 & x_{2} & y_{2}\\ 1 & x_{3} & y_{3}\end {pmatrix}\begin {pmatrix} c_{1}\\ c_{2}\\ c_{3}\end {pmatrix}$

Therefore

\begin {align} \begin {pmatrix} c_{1}\\ c_{2}\\ c_{3}\end {pmatrix} & =\begin {pmatrix} 1 & x_{1} & y_{1}\\ 1 & x_{2} & y_{2}\\ 1 & x_{3} & y_{3}\end {pmatrix} ^{-1}\begin {pmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {pmatrix} \nonumber \\ & =\frac {1}{\left \vert D\right \vert }\begin {pmatrix} x_{2}y_{3}-x_{3}y_{2} & x_{3}y_{1}-x_{1}y_{3} & x_{1}y_{2}-y_{1}x_{2}\\ y_{2}-y_{3} & y_{3}-y_{1} & y_{1}-y_{2}\\ x_{3}-x_{2} & x_{1}-x_{3} & x_{2}-x_{1}\end {pmatrix}\begin {pmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {pmatrix} \tag {2} \end {align}

Where $$\left \vert D\right \vert$$ is the determinant of $$\begin {pmatrix} 1 & x_{1} & y_{1}\\ 1 & x_{2} & y_{2}\\ 1 & x_{3} & y_{3}\end {pmatrix}$$ which is $$\left \vert D\right \vert =-y_{1}x_{2}+x_{3}y_{1}+x_{1}y_{2}-x_{3}y_{2}-x_{1}y_{3}+x_{2}y_{3}$$.

This quantity is twice the area of the triangle $$A=\frac {1}{2}\left ( x_{1}\left ( y_{2}-y_{3}\right ) +x_{2}\left ( y_{3}-y_{1}\right ) +x_{3}\left ( y_{1}-y_{2}\right ) \right )$$.

Substituting (2) into (1) in order to ﬁnd the shape functions gives\begin {align*} \tilde {u} & =c_{1}+c_{2}x+c_{3}y\\ & =\frac {1}{2A}\left [ \left ( x_{2}y_{3}-x_{3}y_{2}\right ) \tilde {u}_{1}+\left ( x_{3}y_{1}-x_{1}y_{3}\right ) \tilde {u}_{2}+\left ( x_{1}y_{2}-y_{1}x_{2}\right ) \tilde {u}_{3}\right ] \\ & +\frac {1}{2A}\left [ \left ( y_{2}-y_{3}\right ) \tilde {u}_{1}+\left ( y_{3}-y_{1}\right ) \tilde {u}_{2}+\left ( y_{1}-y_{2}\right ) \tilde {u}_{3}\right ] x\\ & +\frac {1}{2A}\left [ \left ( x_{3}-x_{2}\right ) \tilde {u}_{1}+\left ( x_{1}-x_{3}\right ) \tilde {u}_{2}+\left ( x_{2}-x_{1}\right ) \tilde {u}_{3}\right ] y \end {align*}

Collecting terms in $$\tilde {u}_{i}$$

\begin {align*} \tilde {u} & =\tilde {u}_{1}\frac {1}{2A}\left ( \left ( x_{2}y_{3}-x_{3}y_{2}\right ) +\left ( y_{2}-y_{3}\right ) x+\left ( x_{3}-x_{2}\right ) y\right ) \\ & +\tilde {u}_{2}\frac {1}{2A}\left ( \left ( x_{3}y_{1}-x_{1}y_{3}\right ) +\left ( y_{3}-y_{1}\right ) x+\left ( x_{1}-x_{3}\right ) y\right ) \\ & +\tilde {u}_{3}\frac {1}{2A}\left ( \left ( x_{1}y_{2}-y_{1}x_{2}\right ) +\left ( y_{1}-y_{2}\right ) x+\left ( x_{2}-x_{1}\right ) y\right ) \end {align*}

Hence

\begin {align*} N_{1}\left ( x,y\right ) & =\frac {1}{2A}\left ( \left ( x_{2}y_{3}-x_{3}y_{2}\right ) +\left ( y_{2}-y_{3}\right ) x+\left ( x_{3}-x_{2}\right ) y\right ) \\ N_{2}\left ( x,y\right ) & =\frac {1}{2A}\left ( \left ( x_{3}y_{1}-x_{1}y_{3}\right ) +\left ( y_{3}-y_{1}\right ) x+\left ( x_{1}-x_{3}\right ) y\right ) \\ N_{3}\left ( x,y\right ) & =\frac {1}{2A}\left ( \left ( x_{1}y_{2}-y_{1}x_{2}\right ) +\left ( y_{1}-y_{2}\right ) x+\left ( x_{2}-x_{1}\right ) y\right ) \end {align*}

Therefore (1) can be written as

$\tilde {u}=N_{1}\left ( x,y\right ) \tilde {u}_{1}+N_{2}\left ( x,y\right ) \tilde {u}_{2}+N_{3}\left ( x,y\right ) \tilde {u}_{3}$

Now that the shape functions are found, the test or weight function $$w\left ( x,y\right )$$ can be found. Since $$w_{i}\left ( x,y\right ) =\frac {\partial \tilde {u}}{\partial \tilde {u}_{i}}$$ for $$i=1,2,3$$, therefore this results in

$w\left ( x,y\right ) =\begin {pmatrix} w_{1}\\ w_{2}\\ w_{3}\end {pmatrix} =\begin {pmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \end {pmatrix} =\frac {1}{2A}\begin {pmatrix} \left ( x_{2}y_{3}-x_{3}y_{2}\right ) +\left ( y_{2}-y_{3}\right ) x+\left ( x_{3}-x_{2}\right ) y\\ \left ( x_{3}y_{1}-x_{1}y_{3}\right ) +\left ( y_{3}-y_{1}\right ) x+\left ( x_{1}-x_{3}\right ) y\\ \left ( x_{1}y_{2}-y_{1}x_{2}\right ) +\left ( y_{1}-y_{2}\right ) x+\left ( x_{2}-x_{1}\right ) y \end {pmatrix}$

The weak form formulation now gives

$I_{i}=-\int \limits _{\Omega }\frac {\partial N_{i}}{\partial x}\frac {\partial N_{j}}{\partial x}d\Omega -\int \limits _{\Omega }\frac {\partial N_{i}}{\partial y}\frac {\partial N_{j}}{\partial y}d\Omega +\int \limits _{\Gamma }w\frac {\partial \tilde {u}}{\partial n}d\Gamma -\int \limits _{\Omega }\begin {pmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \end {pmatrix} f\left ( x,y\right ) d\Omega$

Hence \begin {align*} I_{i} & =-\int \limits _{\Omega }\begin {pmatrix} \frac {\partial N_{1}\left ( x,y\right ) }{\partial x}\\ \frac {\partial N_{2}\left ( x,y\right ) }{\partial x}\\ \frac {\partial N_{3}\left ( x,y\right ) }{\partial x}\end {pmatrix}\begin {Bmatrix} \frac {\partial N_{1}\left ( x,y\right ) }{\partial x} & \frac {\partial N_{2}\left ( x,y\right ) }{\partial x} & \frac {\partial N_{3}\left ( x,y\right ) }{\partial x}\end {Bmatrix}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} d\Omega \\ & -\int \limits _{\Omega }\begin {pmatrix} \frac {\partial N_{1}\left ( x,y\right ) }{\partial y}\\ \frac {\partial N_{2}\left ( x,y\right ) }{\partial y}\\ \frac {\partial N_{3}\left ( x,y\right ) }{\partial y}\end {pmatrix}\begin {Bmatrix} \frac {\partial N_{1}\left ( x,y\right ) }{\partial y} & \frac {\partial N_{2}\left ( x,y\right ) }{\partial y} & \frac {\partial N_{3}\left ( x,y\right ) }{\partial y}\end {Bmatrix}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} d\Omega \\ & +\int \limits _{\Gamma }\begin {pmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \end {pmatrix} \frac {\partial \tilde {u}}{\partial n}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} d\Gamma \\ & -\int \limits _{\Omega }\begin {pmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \end {pmatrix} f\left ( x,y\right ) d\Omega \end {align*}

Since $$w=0$$ on the boundary with Dirichlet conditions, and since there are no natural boundary conditions, the above simpliﬁes to

\begin {align*} I_{i} & =-\int \limits _{\Omega }\begin {pmatrix} \frac {\partial N_{1}\left ( x,y\right ) }{\partial x}\\ \frac {\partial N_{2}\left ( x,y\right ) }{\partial x}\\ \frac {\partial N_{3}\left ( x,y\right ) }{\partial x}\end {pmatrix}\begin {Bmatrix} \frac {\partial N_{1}\left ( x,y\right ) }{\partial x} & \frac {\partial N_{2}\left ( x,y\right ) }{\partial x} & \frac {\partial N_{3}\left ( x,y\right ) }{\partial x}\end {Bmatrix}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} d\Omega \\ & -\int \limits _{\Omega }\begin {pmatrix} \frac {\partial N_{1}\left ( x,y\right ) }{\partial y}\\ \frac {\partial N_{2}\left ( x,y\right ) }{\partial y}\\ \frac {\partial N_{3}\left ( x,y\right ) }{\partial y}\end {pmatrix}\begin {Bmatrix} \frac {\partial N_{1}\left ( x,y\right ) }{\partial y} & \frac {\partial N_{2}\left ( x,y\right ) }{\partial y} & \frac {\partial N_{3}\left ( x,y\right ) }{\partial y}\end {Bmatrix} d\Omega \begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} -\int \limits _{\Omega }\begin {pmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \end {pmatrix} f\left ( x,y\right ) d\Omega \end {align*}

Now each integral is calculated for the triangle element. First the derivatives are found

$\begin {pmatrix} \frac {\partial N_{1}\left ( x,y\right ) }{\partial x}\\ \frac {\partial N_{2}\left ( x,y\right ) }{\partial x}\\ \frac {\partial N_{3}\left ( x,y\right ) }{\partial x}\end {pmatrix} =\frac {1}{2A}\begin {pmatrix} \left ( y_{2}-y_{3}\right ) \\ \left ( y_{3}-y_{1}\right ) \\ \left ( y_{1}-y_{2}\right ) \end {pmatrix}$

and

$\begin {pmatrix} \frac {\partial N_{1}\left ( x,y\right ) }{\partial y}\\ \frac {\partial N_{2}\left ( x,y\right ) }{\partial y}\\ \frac {\partial N_{3}\left ( x,y\right ) }{\partial y}\end {pmatrix} =\frac {1}{2A}\begin {pmatrix} \left ( x_{3}-x_{2}\right ) \\ \left ( x_{1}-x_{3}\right ) \\ \left ( x_{2}-x_{1}\right ) \end {pmatrix}$

Hence

\begin {align*} \begin {pmatrix} \frac {\partial N_{1}\left ( x,y\right ) }{\partial x}\\ \frac {\partial N_{2}\left ( x,y\right ) }{\partial x}\\ \frac {\partial N_{3}\left ( x,y\right ) }{\partial x}\end {pmatrix}\begin {Bmatrix} \frac {\partial N_{1}\left ( x,y\right ) }{\partial x} & \frac {\partial N_{2}\left ( x,y\right ) }{\partial x} & \frac {\partial N_{3}\left ( x,y\right ) }{\partial x}\end {Bmatrix} & =\frac {\partial N_{i}}{\partial x}\frac {\partial N_{j}}{\partial x}\\ & =\frac {1}{\left ( 2A\right ) ^{2}}\begin {pmatrix} \left ( y_{2}-y_{3}\right ) ^{2} & \left ( y_{2}-y_{3}\right ) \left ( y_{3}-y_{1}\right ) & \left ( y_{2}-y_{3}\right ) \left ( y_{1}-y_{2}\right ) \\ \left ( y_{3}-y_{1}\right ) \left ( y_{2}-y_{3}\right ) & \left ( y_{3}-y_{1}\right ) ^{2} & \left ( y_{3}-y_{1}\right ) \left ( y_{1}-y_{2}\right ) \\ \left ( y_{1}-y_{2}\right ) \left ( y_{2}-y_{3}\right ) & \left ( y_{1}-y_{2}\right ) \left ( y_{3}-y_{1}\right ) & \left ( y_{1}-y_{2}\right ) ^{2}\end {pmatrix} \end {align*}

And

\begin {align*} \begin {pmatrix} \frac {\partial N_{1}\left ( x,y\right ) }{\partial y}\\ \frac {\partial N_{2}\left ( x,y\right ) }{\partial y}\\ \frac {\partial N_{3}\left ( x,y\right ) }{\partial y}\end {pmatrix}\begin {Bmatrix} \frac {\partial N_{1}\left ( x,y\right ) }{\partial y} & \frac {\partial N_{2}\left ( x,y\right ) }{\partial y} & \frac {\partial N_{3}\left ( x,y\right ) }{\partial y}\end {Bmatrix} & =\frac {\partial N_{i}}{\partial y}\frac {\partial N_{j}}{\partial y}\\ & =\frac {1}{\left ( 2A\right ) ^{2}}\begin {pmatrix} \left ( x_{3}-x_{2}\right ) ^{2} & \left ( x_{3}-x_{2}\right ) \left ( x_{1}-x_{3}\right ) & \left ( x_{3}-x_{2}\right ) \left ( x_{2}-x_{1}\right ) \\ \left ( x_{1}-x_{3}\right ) \left ( x_{3}-x_{2}\right ) & \left ( x_{1}-x_{3}\right ) ^{2} & \left ( x_{1}-x_{3}\right ) \left ( x_{2}-x_{1}\right ) \\ \left ( x_{2}-x_{1}\right ) \left ( x_{3}-x_{2}\right ) & \left ( x_{2}-x_{1}\right ) \left ( x_{1}-x_{3}\right ) & \left ( x_{2}-x_{1}\right ) ^{2}\end {pmatrix} \end {align*}

We see the above terms do not depend on $$x$$ nor on $$y$$. Hence the next integration is done over the area of the triangle as follow

$\int \limits _{\Omega }\frac {\partial N_{i}}{\partial x}\frac {\partial N_{j}}{\partial x}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} d\Omega =\frac {\partial N_{i}}{\partial x}\frac {\partial N_{j}}{\partial x}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} \int \limits _{\Omega }d\Omega$

But $$\int \limits _{\Omega }d\Omega$$ is the area of the triangle. Hence

$\int \limits _{\Omega }\frac {\partial N_{i}}{\partial x}\frac {\partial N_{j}}{\partial x}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} d\Omega =\frac {\partial N_{i}}{\partial x}\frac {\partial N_{j}}{\partial x}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} A$

Replacing $$\frac {\partial N_{i}}{\partial x}\frac {\partial N_{j}}{\partial x}$$ by the expression we found for it above, gives

\begin {align} \int \limits _{\Omega }\frac {\partial N_{i}}{\partial x}\frac {\partial N_{j}}{\partial x}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} d\Omega & =\frac {A}{\left ( 2A\right ) ^{2}}\begin {pmatrix} \left ( y_{2}-y_{3}\right ) ^{2} & \left ( y_{2}-y_{3}\right ) \left ( y_{3}-y_{1}\right ) & \left ( y_{2}-y_{3}\right ) \left ( y_{1}-y_{2}\right ) \\ \left ( y_{3}-y_{1}\right ) \left ( y_{2}-y_{3}\right ) & \left ( y_{3}-y_{1}\right ) ^{2} & \left ( y_{3}-y_{1}\right ) \left ( y_{1}-y_{2}\right ) \\ \left ( y_{1}-y_{2}\right ) \left ( y_{2}-y_{3}\right ) & \left ( y_{1}-y_{2}\right ) \left ( y_{3}-y_{1}\right ) & \left ( y_{1}-y_{2}\right ) ^{2}\end {pmatrix}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} \nonumber \\ & =\frac {1}{4A}\begin {pmatrix} \left ( y_{2}-y_{3}\right ) ^{2} & \left ( y_{2}-y_{3}\right ) \left ( y_{3}-y_{1}\right ) & \left ( y_{2}-y_{3}\right ) \left ( y_{1}-y_{2}\right ) \\ \left ( y_{3}-y_{1}\right ) \left ( y_{2}-y_{3}\right ) & \left ( y_{3}-y_{1}\right ) ^{2} & \left ( y_{3}-y_{1}\right ) \left ( y_{1}-y_{2}\right ) \\ \left ( y_{1}-y_{2}\right ) \left ( y_{2}-y_{3}\right ) & \left ( y_{1}-y_{2}\right ) \left ( y_{3}-y_{1}\right ) & \left ( y_{1}-y_{2}\right ) ^{2}\end {pmatrix}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} \tag {3} \end {align}

Now we do the same for the second integral

$\int \limits _{\Omega }\frac {\partial N_{i}}{\partial y}\frac {\partial N_{j}}{\partial y}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} d\Omega =\frac {\partial N_{i}}{\partial y}\frac {\partial N_{j}}{\partial y}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} \int \limits _{\Omega }d\Omega$

But $$\int \limits _{\Omega }d\Omega$$ is the area of the triangle. Hence

$\int \limits _{\Omega }\frac {\partial N_{i}}{\partial y}\frac {\partial N_{j}}{\partial y}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} d\Omega =\frac {\partial N_{i}}{\partial y}\frac {\partial N_{j}}{\partial y}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} A$

Replacing $$\frac {\partial N_{i}}{\partial y}\frac {\partial N_{j}}{\partial y}$$ by the expression we found for it above, gives

\begin {align} \int \limits _{\Omega }\frac {\partial N_{i}}{\partial x}\frac {\partial N_{j}}{\partial x}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} d\Omega & =\frac {A}{\left ( 2A\right ) ^{2}}\begin {pmatrix} \left ( x_{3}-x_{2}\right ) ^{2} & \left ( x_{3}-x_{2}\right ) \left ( x_{1}-x_{3}\right ) & \left ( x_{3}-x_{2}\right ) \left ( x_{2}-x_{1}\right ) \\ \left ( x_{1}-x_{3}\right ) \left ( x_{3}-x_{2}\right ) & \left ( x_{1}-x_{3}\right ) ^{2} & \left ( x_{1}-x_{3}\right ) \left ( x_{2}-x_{1}\right ) \\ \left ( x_{2}-x_{1}\right ) \left ( x_{3}-x_{2}\right ) & \left ( x_{2}-x_{1}\right ) \left ( x_{1}-x_{3}\right ) & \left ( x_{2}-x_{1}\right ) ^{2}\end {pmatrix}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} \nonumber \\ & =\frac {1}{4A}\begin {pmatrix} \left ( x_{3}-x_{2}\right ) ^{2} & \left ( x_{3}-x_{2}\right ) \left ( x_{1}-x_{3}\right ) & \left ( x_{3}-x_{2}\right ) \left ( x_{2}-x_{1}\right ) \\ \left ( x_{1}-x_{3}\right ) \left ( x_{3}-x_{2}\right ) & \left ( x_{1}-x_{3}\right ) ^{2} & \left ( x_{1}-x_{3}\right ) \left ( x_{2}-x_{1}\right ) \\ \left ( x_{2}-x_{1}\right ) \left ( x_{3}-x_{2}\right ) & \left ( x_{2}-x_{1}\right ) \left ( x_{1}-x_{3}\right ) & \left ( x_{2}-x_{1}\right ) ^{2}\end {pmatrix}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} \tag {4} \end {align}

Adding (3) and (4) to obtain the local stiﬀness matrix

\begin {align*} k_{i}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} _{i} & =\frac {1}{4A}\begin {pmatrix} \left ( y_{2}-y_{3}\right ) ^{2} & \left ( y_{2}-y_{3}\right ) \left ( y_{3}-y_{1}\right ) & \left ( y_{2}-y_{3}\right ) \left ( y_{1}-y_{2}\right ) \\ \left ( y_{3}-y_{1}\right ) \left ( y_{2}-y_{3}\right ) & \left ( y_{3}-y_{1}\right ) ^{2} & \left ( y_{3}-y_{1}\right ) \left ( y_{1}-y_{2}\right ) \\ \left ( y_{1}-y_{2}\right ) \left ( y_{2}-y_{3}\right ) & \left ( y_{1}-y_{2}\right ) \left ( y_{3}-y_{1}\right ) & \left ( y_{1}-y_{2}\right ) ^{2}\end {pmatrix}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} \\ & +\frac {1}{4A}\begin {pmatrix} \left ( x_{3}-x_{2}\right ) ^{2} & \left ( x_{3}-x_{2}\right ) \left ( x_{1}-x_{3}\right ) & \left ( x_{3}-x_{2}\right ) \left ( x_{2}-x_{1}\right ) \\ \left ( x_{1}-x_{3}\right ) \left ( x_{3}-x_{2}\right ) & \left ( x_{1}-x_{3}\right ) ^{2} & \left ( x_{1}-x_{3}\right ) \left ( x_{2}-x_{1}\right ) \\ \left ( x_{2}-x_{1}\right ) \left ( x_{3}-x_{2}\right ) & \left ( x_{2}-x_{1}\right ) \left ( x_{1}-x_{3}\right ) & \left ( x_{2}-x_{1}\right ) ^{2}\end {pmatrix}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\end {Bmatrix} \end {align*}

or

$k_{i}=\frac {1}{4A}\begin {pmatrix} \left ( y_{2}-y_{3}\right ) ^{2}+\left ( x_{3}-x_{2}\right ) ^{2} & \left ( y_{2}-y_{3}\right ) \left ( y_{3}-y_{1}\right ) +\left ( x_{3}-x_{2}\right ) \left ( x_{1}-x_{3}\right ) & \left ( y_{2}-y_{3}\right ) \left ( y_{1}-y_{2}\right ) +\left ( x_{3}-x_{2}\right ) \left ( x_{2}-x_{1}\right ) \\ \left ( y_{3}-y_{1}\right ) \left ( y_{2}-y_{3}\right ) +\left ( x_{1}-x_{3}\right ) \left ( x_{3}-x_{2}\right ) & \left ( y_{3}-y_{1}\right ) ^{2}+\left ( x_{1}-x_{3}\right ) ^{2} & \left ( y_{3}-y_{1}\right ) \left ( y_{1}-y_{2}\right ) +\left ( x_{1}-x_{3}\right ) \left ( x_{2}-x_{1}\right ) \\ \left ( y_{1}-y_{2}\right ) \left ( y_{2}-y_{3}\right ) +\left ( x_{2}-x_{1}\right ) \left ( x_{3}-x_{2}\right ) & \left ( y_{1}-y_{2}\right ) \left ( y_{3}-y_{1}\right ) +\left ( x_{2}-x_{1}\right ) \left ( x_{1}-x_{3}\right ) & \left ( y_{1}-y_{2}\right ) ^{2}+\left ( x_{2}-x_{1}\right ) ^{2}\end {pmatrix}$

And the ﬁnal integral (the force vector) depends on the nature of $$f\left ( x,y\right )$$. For $$f\left ( x,y\right ) =1$$ we obtain

$\int \limits _{\Omega }\begin {pmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \end {pmatrix} f\left ( x,y\right ) d\Omega =\int \limits _{\Omega }\frac {1}{2A}\begin {pmatrix} \left ( x_{2}y_{3}-x_{3}y_{2}\right ) +\left ( y_{2}-y_{3}\right ) x+\left ( x_{3}-x_{2}\right ) y\\ \left ( x_{3}y_{1}-x_{1}y_{3}\right ) +\left ( y_{3}-y_{1}\right ) x+\left ( x_{1}-x_{3}\right ) y\\ \left ( x_{1}y_{2}-y_{1}x_{2}\right ) +\left ( y_{1}-y_{2}\right ) x+\left ( x_{2}-x_{1}\right ) y \end {pmatrix} dydx$

And for $$f\left ( x,y\right ) =xy\,\$$similarly

$\int \limits _{\Omega }\begin {pmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \end {pmatrix} f\left ( x,y\right ) d\Omega =\int \limits _{\Omega }\frac {1}{2A}\begin {pmatrix} \left ( x_{2}y_{3}-x_{3}y_{2}\right ) +\left ( y_{2}-y_{3}\right ) x+\left ( x_{3}-x_{2}\right ) y\\ \left ( x_{3}y_{1}-x_{1}y_{3}\right ) +\left ( y_{3}-y_{1}\right ) x+\left ( x_{1}-x_{3}\right ) y\\ \left ( x_{1}y_{2}-y_{1}x_{2}\right ) +\left ( y_{1}-y_{2}\right ) x+\left ( x_{2}-x_{1}\right ) y \end {pmatrix} xydydx$

And similarly for other $$f\left ( x,y\right )$$. These were all integrated oﬀ-line and the resulting expression evaluated during the run of the program. These expression are all functions of the node coordinates $$\left ( x_{1},y_{1},x_{2},y_{2},x_{3},y_{3}\right )$$. No actual integration is done in the code, but only evaluation of the integration result obtained. The integration was done using Mathematica. These are the results for the three functions of interest. The integration was done on two triangles. The odd numbered ones and the even numbered ones due to the node numbering diﬀerence. For odd numbered triangles, the integration limits are $$\int \limits _{x_1}^{x_2}\int \limits _{y_1}^{y_{3}-x}dydx$$ while on the even numbered elements the limits are $$\int \limits _{x_{2}}^{x_{1}}\int \limits _{y_1}^{y_{2}-x}dydx$$ as shown in this diagram

Note that this is not how integration is normally done in actual ﬁnite element methods. Using Gaussian quadrature method is the recommended way. This method was done here for learning and to compare.

#### 4.16.1 case 1

Case $$f\left ( x,y\right )=1$$. For odd numbered elements

n1 = (1/(2*area))*((x2*y3 - x3*y2) + (y2 - y3)*x + (x3 - x2)*y);
n2 = (1/(2*area))*((x2*y1 - x1*y3) + (y3 - y1)*x + (x1 - x3)*y);
n3 = (1/(2*area))*((x1*y2 - y1*x2) + (y1 - y2)*x + (x2 - x1)*y);
w = {{n1}, {n2}, {n3}};
Integrate[w, {x, x1, x2}, {y, y1, y3 - x}]

Out[11]= {{(1/(12*area))*((x1 - x2)*(x2^3 + x1^2*(x2 - x3 + 2*y2 - 2*y3) +
3*x3*(y1 - y3)*(y1 - 2*y2 + y3) - x2^2*(x3 - 2*y2 + 2*y3) -
3*x2*(y1^2 + x3*y2 - x3*y3 + y2*y3 - y1*(y2 + y3)) + x1*
(x2^2 - 3*(y2 - y3)*(x3 - y1 + y3) - x2*(x3 - 2*y2 + 2*y3))))},
{-(((x1 - x2)*(x1^3 + x1^2*(x2 - x3 + 2*y1 - 2*y3) - x2^2*(x3 + y1 + 2*y3) +
3*x3*(y1^2 - y3^2) + 3*x2*(-y1^2 + y3*(x3 + y3)) +
x1*(x2^2 + 3*x3*y3 - x2*(x3 + y1 + 2*y3))))/(12*area))},
{((x1 - x2)^2*(x1^2 + x2^2 + x1*(x2 + 2*y1 + y2 - 3*y3) + x2*(y1 + 2*y2 - 3*y3) +
3*(y1 - y3)*(y2 - y3)))/(12*area)}}



For even numbered elements

Integrate[w, {x, x1, x2}, {y, y1, y2 - x}]
{{((x1 - x2)*(x2^3 + 3*x3*(y1 - y2)^2 + x1*(x2^2 + 3*(y1 - y2)*(y2 - y3) -
x2*(x3 + y2 - y3)) + x1^2*(x2 - x3 + 2*y2 - 2*y3) -
3*x2*(y1 - y2)*(y1 - y3) - x2^2*(x3 + y2 - y3)))/(12*area)},
{-((1/(12*area))*((x1 - x2)*(x1^3 - 3*x3*(y1 - y2)^2 - 3*x2*(y1 - y2)*
(x3 - y1 + y3) + x1^2*(x2 - x3 + 2*y1 - 3*y2 + y3) -
x2^2*(x3 - 2*y1 + 2*y3) + x1*(x2^2 - 3*(y1 - y2)*(x3 + y2 - y3) +
x2*(-x3 + 2*y1 - 3*y2 + y3)))))},
{((x1 - x2)^2*(x1^2 + x1*(x2 + 2*y1 - 2*y2) + x2*(x2 + y1 - y2)))/(12*area)}
}



#### 4.16.2 Program output

This is an implementation of the above for $$f(x,y)=1$$. The coordinates mapping was ﬁrst carried out in order to map each elements local node numbers to the global node numbering so we know where to add each element stiﬀness matrix to the global stiﬀness matrix. Also we need to build the mapping of each element to the physical coordinates of its local node number 1 to use for calculating the force vector since that depends on physical location of each element.

The mapping of physical coordinate location relative to the global frame of reference is shown below

function matlab_98()
close all; clear all;
number_of_elements=[16*2 64*2 256*2 625*2];
function_string={'1'};
%function_handles={@f_1,@f_xy,@f_2};%something wrong with the other
%cases, need to work more on it.
function_handles={@f_1};
figure_image_file={'matlab_e98_1'};

for i=1:length(function_handles)
for j=1:length(number_of_elements)
plot_file=sprintf('%s_%d',figure_image_file{i},...
number_of_elements(j));
close all;
solve_poisson_2D_square(function_handles{i},...
function_string{i},number_of_elements(j),plot_file);
end
end
end

%----------------------------------
function k_local=getk(area,x1,y1,x2,y2,x3,y3)

k11=(y2-y3)^2+(x3-x2)^2;
k12=(y2-y3)*(y3-y1)+(x3-x2)*(x1-x3);
k13=(y2-y3)*(y1-y2)+(x3-x2)*(x2-x1);
k22=(y3-y1)^2+(x1-x3)^2;
k23=(y3-y1)*(y1-y2)+(x1-x3)*(x2-x1);
k33=(y1-y2)^2+(x2-x1)^2;
k_local = [k11  k12 k13;
k12  k22 k23;
k13  k23 k33];

k_local = k_local/(4*area);
end
%----------------------------------------------------
%use this for f(x,y)=1
function f_local=f_1(area,x1,y1,x2,y2,x3,y3,ele)

if mod(ele,2)
f_local=[(1/(12*area))*((x1 - x2)*(x2^3 + x1^2*...
(x2 - x3 + 2*y2 - 2*y3) +...
3*x3*(y1 - y3)*(y1 - 2*y2 + y3) - x2^2*(x3 - 2*y2 + 2*y3) -...
3*x2*(y1^2 + x3*y2 - x3*y3 + y2*y3 - y1*(y2 + y3)) + x1*(x2^2 -...
3*(y2 - y3)*(x3 - y1 + y3) - x2*(x3 - 2*y2 + 2*y3))));
-((1/(12*area))*((x1 - x2)*(x1^3 + x1^2*(x2 - x3 + 2*y1 - 2*y3) -...
3*x3*(y1 - y3)^2 - 3*x2*(y1 - y3)*(x3 - y1 + y3) -...
x2^2*(x3 - 2*y1 + 2*y3) + x1*(x2^2 + 3*x3*(-y1 + y3) -...
x2*(x3 - 2*y1 + 2*y3)))));
((x1 - x2)^2*(x1^2 + x2^2 + x1*(x2 + 2*y1 + y2 - 3*y3) +...
x2*(y1 + 2*y2 - 3*y3) + 3*(y1 - y3)*(y2 - y3)))/(12*area)];
else
f_local=[((x1 - x2)*(x2^3 + 3*x3*(y1 - y2)^2 + x1*(x2^2 +...
3*(y1 - y2)*(y2 - y3) - x2*(x3 + y2 - y3)) + x1^2*...
(x2 - x3 + 2*y2 - 2*y3) -...
3*x2*(y1 - y2)*(y1 - y3) - x2^2*(x3 + y2 - y3)))/(12*area);
-((1/(12*area))*((x1 - x2)*(x1^3 - 3*x3*(y1 - y2)^2 -...
3*x2*(y1 - y2)*(x3 - y1 + y3) + x1^2*...
(x2 - x3 + 2*y1 - 3*y2 + y3) -...
x2^2*(x3 - 2*y1 + 2*y3) + x1*(x2^2 - 3*(y1 - y2)*(x3 + y2 - y3)...
+ x2*(-x3 + 2*y1 - 3*y2 + y3)))));
((x1 - x2)^2*(x1^2 + x1*(x2 + 2*y1 - 2*y2) + ...
x2*(x2 + y1 - y2)))/(12*area)];

end

end
%----------------------------------------------------
function solve_poisson_2D_square(fh,rhs,M,plot_file)
%fh is function which evaluate the force vector
%rhs is string which is the f(x,y), for plotting only
%M=total number of elements
L  = 1; %length of grid
N  = 3; %number of nodes per element
h=(L/sqrt(M/2));
%area is same for each element, otherwise use
%area = 0.5*(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2));
area= (1/2)*h^2;

N1 = sqrt(M/2)+1; %number of nodes per edge
dof=3;
nodes = N1^2; %total number of nodes
kk    = zeros(nodes); %global stiffness vector
f     = zeros(nodes,1); %global force vector
mtb   = generate_coordinates_table(M);
fmtb  = generate_physical_coordinates_table(M,h);

for i = 1:M
x1 = fmtb(i,1); y1 = fmtb(i,2);
x2 = fmtb(i,3); y2 = fmtb(i,4);
x3 = fmtb(i,5); y3 = fmtb(i,6);
k_local=getk(area,x1,y1,x2,y2,x3,y3);
%call to load different force vector, pre-computed offline
f_local=fh(area,x1,y1,x2,y2,x3,y3,i);

for ii = 1:N %assemble force vector and global stiffness matrix
f(mtb(i,ii)) = f(mtb(i,ii)) + f_local(ii);
for jj = 1:N
kk(mtb(i,ii),mtb(i,jj)) = kk(mtb(i,ii),mtb(i,jj))+...
k_local(ii,jj);
end
end
end

%fix the global stiffness matrix and force vector for B.C.
edge_nodes=[1:N1 (N1+1):N1:nodes-N1+1 ...
(2*N1):N1:nodes nodes-N1+2:nodes-1];
for i=1:length(edge_nodes)
z=edge_nodes(i);
kk(z,:)=0;
f(z)=0;
end

for i=1:length(edge_nodes)
z=edge_nodes(i);
kk(z,z)=1;
end

y=kk\f; %solve
Z=reshape(y,N1,N1);
[X,Y] = meshgrid(0:h:1,0:h:1);
h=surf(X,Y,Z);
set(h','edgecolor','k');
set(gcf,'defaulttextinterpreter','latex');
plot_title=sprintf('$\\nabla^2 u=%s$, number of elements $%d$',rhs,M);
title(plot_title,'fontweight','bold','fontsize',14);

print(gcf, '-dpdf', '-r600', ['images/',plot_file]);
print(gcf, '-dpng', '-r300', ['images/',plot_file]);

figure;
[C,h] = contourf(X,Y,Z);
clabel(C,h)
set(gcf,'defaulttextinterpreter','latex');
title(plot_title,'fontweight','bold','fontsize',14);
colormap cool

print(gcf, '-dpdf', '-r300', ['images/',plot_file,'_c']);
print(gcf, '-dpng', '-r300', ['images/',plot_file,'_c']);

end
%----------------------------------------------------
function A=generate_coordinates_table(M)
%M=number_of_elements
A=zeros(M,3);
n=2*sqrt(M/2); %number of elements per row

for i=1:n/2
for j=1:n
ele=j+n*(i-1);
if j==1
A(ele,:)=[(n/2+1)*(i-1)+1 (n/2+1)*(i-1)+2 (n/2+1)*i+1];
else
if mod(j,2)
A(ele,:) =[A(ele-1,3) A(ele-1,3)+1 A(ele-1,1)];
else
A(ele,:) =[A(ele-1,3)+1 A(ele-1,3) A(ele-1,2)];
end
end
end
end
end
%----------------------------------------------------
function A=generate_physical_coordinates_table(M,h)
%M=number_of_elements
A=zeros(M,6);
n=2*sqrt(M/2); %number of elements per row

for i=1:n/2
for j=1:n
ele=j+n*(i-1);
if j==1
A(ele,:)=[0 (i-1)*h  h (i-1)*h  0 i*h];
else
if mod(j,2)
A(ele,:) =[A(ele-1,5) A(ele-1,6) ...
A(ele-1,5)+h A(ele-1,6)  A(ele-1,1) A(ele-1,2)];
else
A(ele,:) =[A(ele-1,5)+h A(ele-1,6) ...
A(ele-1,5) A(ele-1,6)  A(ele-1,3) A(ele-1,4)];
end
end
end
end
end


Here is the output for the three force functions, for diﬀerent number of elements.

#### 4.16.3 case 2

Case $$f(x,y)=1$$