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.
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 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
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.
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)} }
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.
Case \(f(x,y)=1\)
|
|
|
|
|
|
|
|
|
|
|
|
