4.16 How to solve Poisson PDE in 2D using finite elements methods using triangle element?

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 different number of elements.

Figure 4.2: node numbering

Looking at one element

Figure 4.3: element in physical coordinates

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 find 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 simplifies 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 stiffness 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 final 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 off-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 difference. 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

Figure 4.4: Integrating over elements

Note that this is not how integration is normally done in actual finite 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 first carried out in order to map each elements local node numbers to the global node numbering so we know where to add each element stiffness matrix to the global stiffness 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); 
shading interp 
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 different number of elements.

4.16.3 case 2

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