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

4.15.2 Case 1
4.15.3 $$f(x,y)=1$$
4.15.4 $$f(x,y)=xy$$
4.15.5 case 3

4.15.7 Case 4
4.15.8 Case 5
4.15.9 case 6

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 a bilinear trial function. Let width of square be 1.

Solution

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

The trial function is bilinear

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

Looking at one element, and using local coordinates systems with element having width $$2a$$ and height $$2b$$ gives

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

\begin {align*} \tilde {u}_{1} & =c_{1}-ac_{2}-bc_{3}+abc_{4}\\ \tilde {u}_{2} & =c_{1}+ac_{2}-bc_{3}-abc_{4}\\ \tilde {u}_{3} & =c_{1}+ac_{2}+by_{3}+abc_{4}\\ \tilde {u}_{4} & =c_{1}-ac_{2}+bc_{3}+abc_{4} \end {align*}

Or

$\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\\ \tilde {u}_{4}\end {Bmatrix} =\begin {pmatrix} 1 & -a & -b & ab\\ 1 & a & -b & -ab\\ 1 & a & b & ab\\ 1 & -a & b & ab \end {pmatrix}\begin {Bmatrix} c_{1}\\ c_{2}\\ c_{3}\\ c_{4}\end {Bmatrix}$

Hence

\begin {align} \begin {Bmatrix} c_{1}\\ c_{2}\\ c_{3}\\ c_{4}\end {Bmatrix} & =\begin {pmatrix} 1 & -a & -b & ab\\ 1 & a & -b & -ab\\ 1 & a & b & ab\\ 1 & -a & b & ab \end {pmatrix} ^{-1}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\\ \tilde {u}_{4}\end {Bmatrix} \nonumber \\ & =\frac {1}{4ab}\begin {pmatrix} ab & ab & ab & ab\\ -b & b & b & -b\\ -a & -a & a & a\\ 1 & -1 & 1 & -1 \end {pmatrix}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\\ \tilde {u}_{4}\end {Bmatrix} \tag {2} \end {align}

coor = {{-a, -b}, {a, -b}, {a, b}, {-a, b}};
uVal = {u1, u2, u3, u4}; eq = c1 + c2 x + c3 y + c4 x y;
r = MapThread[eq /. {u -> #1, x -> #2[[1]], y -> #2[[2]]}&, {uVal, coor}];
mat = Last@Normal[CoefficientArrays[r, {c1, c2, c3, c4}]];
MatrixForm[Inverse[mat] // MatrixForm



Substituting (2) back into (1) and rearranging terms results in

\begin {align*} \tilde {u} & =c_{1}+c_{2}x+c_{3}y+c_{4}xy\\ & =\frac {1}{4ab}\left ( \left ( ab-bx-ay+xy\right ) \tilde {u}_{1}+\left ( ab+bx-ay-xy\right ) \tilde {u}_{2}+\left ( ab+bx+ay+xy\right ) +\left ( ab-bx+ay-xy\right ) \right ) \end {align*}

Since $$ab$$ is the $$\frac {1}{4}$$ of the area of element, then the above becomes

$\tilde {u}=\frac {1}{A}\left ( \left ( ab-bx-ay+xy\right ) \tilde {u}_{1}+\left ( ab+bx-ay-xy\right ) \tilde {u}_{2}+\left ( ab+bx+ay+xy\right ) +\left ( ab-bx+ay-xy\right ) \right )$

The above can now be written in term of what is called shape functions

$\tilde {u}(x,y)=N_{1}(x,y)\tilde {u_{1}}+N_{2}(x,y)\tilde {u_{2}}N_{3}(x,y)\tilde {u_{3}}+N_{4}(x,y)\tilde {u_{4}}$

Where

\begin {align*} N_{1} & =\frac {1}{A}\left ( ab-bx-ay+xy\right ) =\frac {1}{A}\left ( a-x\right ) \left ( b-y\right ) \\ N_{2} & =\frac {1}{A}\left ( ab+bx-ay-xy\right ) =\frac {1}{A}\left ( a+x\right ) \left ( b-y\right ) \\ N_{3} & =\frac {1}{A}\left ( ab+bx+ay+xy\right ) =\frac {1}{A}\left ( a+x\right ) \left ( b+y\right ) \\ N_{4} & =\frac {1}{A}\left ( ab-bx+ay-xy\right ) =\frac {1}{A}\left ( a-x\right ) \left ( b+y\right ) \end {align*}

Now that the shape functions are found, the next step is to determine the local element stiﬀness matrix. This can be found from the weak form integral over the area of the element.

\begin {equation} I_{i}=\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}w\left ( \nabla ^{2}u-f\left ( x,y\right ) \right ) \,dxdy \tag {3} \end {equation}

Where $$w\left ( x,y\right )$$ is the test function. For Galerkin method, the test function $$w_{i}=\frac {d\tilde {u}}{du_{i}}=N_{i}\left ( x,y\right )$$. Hence

$w\left ( x,y\right ) =\begin {Bmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \\ N_{4}\left ( x,y\right ) \end {Bmatrix} =\begin {Bmatrix} \frac {1}{A}\left ( a-x\right ) \left ( b-y\right ) \\ \frac {1}{A}\left ( a+x\right ) \left ( b-y\right ) \\ \frac {1}{A}\left ( a+x\right ) \left ( b+y\right ) \\ \frac {1}{A}\left ( a-x\right ) \left ( b+y\right ) \end {Bmatrix}$

Using weak form, integration by parts is applied to (2)

$I_{i}=\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\left ( -\frac {\partial w}{\partial x_{i}}\frac {\partial \tilde {u}}{\partial x_{i}}-wf\left ( x,y\right ) \right ) \,dxdy+\overbrace {\int \limits _{\Gamma }w\frac {\partial \tilde {u}}{\partial n}}^{\text {goes to zero}}$

The term $$\int \limits _{\Gamma }w\frac {\partial \tilde {u}}{\partial n}$$ is the integration over the boundary of the element. Since there is only an essential boundary condition over all the boundaries (this is the given Dirchilet boundary condition), $$w=0$$ on the boundary and this integral vanishes. There is no natural boundary conditions for this example.

For those elements not on the external edge of the overall grid (i.e. internal elements), each contribution to this integral will cancel from the adjacent internal element. What this means is that the above reduces to just the ﬁrst integral

\begin {align*} I_{i} & =\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\left ( -\frac {\partial w}{\partial x_{i}}\frac {\partial \tilde {u}}{\partial x_{i}}-wf\left ( x,y\right ) \right ) \,dxdy\\ & =\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}-\begin {Bmatrix} \frac {\partial N_{1}}{\partial x}\\ \frac {\partial N_{2}}{\partial x}\\ \frac {\partial N_{3}}{\partial x}\\ \frac {\partial N_{4}}{\partial x}\end {Bmatrix}\begin {Bmatrix} \frac {\partial N_{1}}{\partial x} & \frac {\partial N_{2}}{\partial x} & \frac {\partial N_{3}}{\partial x} & \frac {\partial N_{4}}{\partial x}\end {Bmatrix}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\\ \tilde {u}_{4}\end {Bmatrix} -\begin {Bmatrix} \frac {\partial N_{1}}{\partial y}\\ \frac {\partial N_{2}}{\partial y}\\ \frac {\partial N_{3}}{\partial y}\\ \frac {\partial N_{4}}{\partial y}\end {Bmatrix}\begin {Bmatrix} \frac {\partial N_{1}}{\partial y} & \frac {\partial N_{2}}{\partial y} & \frac {\partial N_{3}}{\partial y} & \frac {\partial N_{4}}{\partial y}\end {Bmatrix}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\\ \tilde {u}_{4}\end {Bmatrix} \\ & -\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\begin {Bmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \\ N_{4}\left ( x,y\right ) \end {Bmatrix} f\left ( x,y\right ) \,dxdy \end {align*}

Hence

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

In the above, we have the from $$I_{i}=\int \limits _{\Omega }-k_{i}\left \{ \tilde {u}\right \} d\Omega -\int \limits _{\Omega }\left \{ \tilde {u}\right \} f_{i}d\Omega$$, hence the element stiﬀness matrix is the ﬁrst integral, and the force vector comes from the second integral.

The integration is now carried out to obtain the element stiﬀness matrix. This was done using Mathematica. The local stiﬀness matrix for element $$i$$ is

\begin {align*} k_{i} & =\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\frac {\partial w}{\partial x_{i}}\frac {\partial \tilde {u}}{\partial x_{i}}dxdy\\ & =\begin {pmatrix} \frac {a^{2}+b^{2}}{3ab} & \frac {a}{6b}-\frac {b}{3a} & -\frac {a^{2}+b^{2}}{6ab} & -\frac {a}{3b}+\frac {b}{6a}\\ \frac {a}{6b}-\frac {b}{3a} & \frac {a^{2}+b^{2}}{3ab} & -\frac {a}{3b}+\frac {b}{6a} & -\frac {a^{2}+b^{2}}{6ab}\\ -\frac {a^{2}+b^{2}}{6ab} & -\frac {a}{3b}+\frac {b}{6a} & \frac {a^{2}+b^{2}}{3ab} & \frac {a}{6b}-\frac {b}{3a}\\ -\frac {a}{3b}+\frac {b}{6a} & -\frac {a^{2}+b^{2}}{6ab} & \frac {a}{6b}-\frac {b}{3a} & \frac {a^{2}+b^{2}}{3ab}\end {pmatrix} \end {align*}

For example, for element of width 1 and height 1, then $$a=\frac {1}{2},b=\frac {1}{2}$$ and the above becomes

$k_{i}= \begin {pmatrix} \frac {2}{3} & -\frac {1}{6} & -\frac {1}{3} & -\frac {1}{6}\\ -\frac {1}{6} & \frac {2}{3} & -\frac {1}{6} & -\frac {1}{3}\\ -\frac {1}{3} & -\frac {1}{6} & \frac {2}{3} & -\frac {1}{6}\\ -\frac {1}{6} & -\frac {1}{3} & -\frac {1}{6} & \frac {2}{3}\end {pmatrix}$

ClearAll[a, b];
area = 4 a b;
phi = (1/area) {(a - x) (b - y),
(a + x) (b - y),
(a + x) (b + y),
(a - x) (b + y)};
vx = {{D[phi[[1]], x]},
{D[phi[[2]], x]},
{D[phi[[3]], x]},
{D[phi[[4]], x]}};
vy = {{D[phi[[1]], y]},
{D[phi[[2]], y]},
{D[phi[[3]], y]},
{D[phi[[4]], y]}};
k = Integrate[vx.Transpose[vx]+vy.Transpose[vy],{x, -a, a}, {y, -b, b}]



Hence

$I_{i}=-k_{i}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\\ \tilde {u}_{4}\end {Bmatrix} -\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\begin {Bmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \\ N_{4}\left ( x,y\right ) \end {Bmatrix} f\left ( x,y\right ) \,dxdy$

Now the integration of the force vector is carried out.

\begin {align*} I_{i}^{f} & =\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\begin {Bmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \\ N_{4}\left ( x,y\right ) \end {Bmatrix} f\left ( x,y\right ) \,dxdy\\ & =\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\begin {Bmatrix} \frac {1}{A}\left ( a-x\right ) \left ( b-y\right ) \\ \frac {1}{A}\left ( a+x\right ) \left ( b-y\right ) \\ \frac {1}{A}\left ( a+x\right ) \left ( b+y\right ) \\ \frac {1}{A}\left ( a-x\right ) \left ( b+y\right ) \end {Bmatrix} f\left ( x,y\right ) \,dxdy \end {align*}

In the above, the integration depends on the physical location of the element. We used a local coordinate system initially to obtain the shape functions. This has now to be transformed to the global coordinates system. Taking the center of each element as $$\left ( x_{0},y_{0}\right )$$ then we need to replace $$a\rightarrow x_{0}+a$$ and $$b\rightarrow y_{0}+b$$ everywhere in the above integral. We did not have to do this for ﬁnding the local stiﬀness matrix, since that did not depend on the physical location of the element (it was constant). But for the integration of the force function, we need to do this mapping. In the code, the center of each element is found, and the replacement is done.

#### 4.15.1 Integrating the force vector

4.15.1.1 $$f\left ( x,y\right ) =xy$$

This is normally done using Gaussian quadrature method. In this example, the integral is found oﬀ-line using Mathematica.

Evaluating the force vector requires integration over the element. This was done oﬀ-line using Mathematica. The result is a function of the coordinate of the center of the element. Then when running the code, this was evaluated for each element in the main loop.

##### 4.15.1.1 $$f\left ( x,y\right ) =xy$$

For $$f\left ( x,y\right ) =xy$$ the above gives a result as function of the center of the physical element.

$I_{i}^{f}=\int \limits _{x=\left ( x_{0}-a\right ) }^{x=x_{0}+a}\int \limits _{y=\left ( y_{0}-b\right ) }^{y=y_{0}+b}\begin {Bmatrix} \frac {1}{A}\left ( \left ( x_{0}-a\right ) -x\right ) \left ( y_{0}+b-y\right ) \\ \frac {1}{A}\left ( \left ( x_{0}-a\right ) +x\right ) \left ( y_{0}+b-y\right ) \\ \frac {1}{A}\left ( \left ( x_{0}-a\right ) +x\right ) \left ( y_{0}+b+y\right ) \\ \frac {1}{A}\left ( \left ( x_{0}-a\right ) -x\right ) \left ( y_{0}+b+y\right ) \end {Bmatrix} xy\,dxdy$

This was done using Mathematica

$I_{i}^{f}=\frac {1}{9\left ( x_{0}+a\right ) \left ( y_{0}+b\right ) }\begin {Bmatrix} a^{2}b^{2}\left ( a-3x_{0}\right ) \left ( b-3y_{0}\right ) \\ -ab^{2}\left ( a^{2}+3ax_{0}+6x_{0}^{2}\right ) \left ( b-3y_{0}\right ) \\ ab\left ( a^{2}+3ax_{0}+6x_{0}^{2}\right ) \left ( b^{2}+3by_{0}+6y_{0}^{2}\right ) \\ -ab^{2}\left ( a-3x_{0}\right ) \left ( b^{2}+3by_{0}+y_{0}^{2}\right ) \end {Bmatrix}$

In[134]:= phi
Out[134]= {((a - x)*(b - y))/(4*a*b),
((a + x)*(b - y))/(4*a*b),
((a + x)*(b + y))/(4*a*b),
((a - x)*(b + y))/(4*a*b)}

In[135]:= phi /. {a :> x0 + a, b :> y0 + b}

Out[135]= {
((a - x + x0)*(b - y + y0))/(4*(a + x0)*(b + y0)),
((a + x + x0)*(b - y + y0))/(4*(a + x0)*(b + y0)),
((a + x + x0)*(b + y + y0))/(4*(a + x0)*(b + y0)),
((a - x + x0)*(b + y + y0))/(4*(a + x0)*(b + y0))}

Integrate[% x*y, {x, x0 - a, x0 + a}, {y, y0 - b, y0 + b}]

Out[136]= {
(a^2*b^2*(a - 3*x0)*(b - 3*y0))/(9*(a + x0)*(b + y0)),
-((a*b^2*(a^2 + 3*a*x0 + 6*x0^2)*(b - 3*y0))/(9*(a + x0)*(b + y0))),
(a*b*(a^2 + 3*a*x0 + 6*x0^2)*(b^2 + 3*b*y0 + 6*y0^2))/(9*(a + x0)*(b + y0)),
-((a^2*b*(a - 3*x0)*(b^2 + 3*b*y0 + 6*y0^2))/(9*(a + x0)*(b + y0)))}



#### 4.15.2 Case 1

4.15.2.1 case 2

Case $$f(x,y)=3\left ( \cos (4 \pi x)+\sin (3 \pi y)\right )$$

This was done using Mathematica

phi
phi /. {a :> (x0 + a), b :> (y0 + b)}
Integrate[%  3 (Cos[4 Pi x] + Sin[3 Pi y]),
{x, x0 - a, x0 + a}, {y, y0 - b, y0 + b}]

{(1/(96*Pi^2*(a + x0)*(b + y0)))*(9*b^2*Cos[4*Pi*(-a + x0)] -
9*b^2*Cos[4*Pi*(a + x0)] + 8*a*(12*a*b*Pi*Cos[3*Pi*(b - y0)] -
9*b^2*Pi*Sin[4*Pi*(-a + x0)] - 2*a*(Sin[3*Pi*(b - y0)] +
Sin[3*Pi*(b + y0)]))), (1/(96*Pi^2*(a + x0)*(b + y0)))*
(-9*b^2*Cos[4*Pi*(-a + x0)] + 9*b^2*Cos[4*Pi*(a + x0)] +
8*(12*a*b*Pi*(a + 2*x0)*Cos[3*Pi*(b - y0)] -
9*b^2*Pi*x0*Sin[4*Pi*(-a + x0)] + 9*a*b^2*Pi*Sin[4*Pi*(a + x0)] +
9*b^2*Pi*x0*Sin[4*Pi*(a + x0)] - 2*a^2*Sin[3*Pi*(b - y0)] -
4*a*x0*Sin[3*Pi*(b - y0)] - 2*a^2*Sin[3*Pi*(b + y0)] -
4*a*x0*Sin[3*Pi*(b + y0)])), (1/(96*Pi^2*(a + x0)*(b + y0)))*
(-9*b*(b + 2*y0)*Cos[4*Pi*(-a + x0)] + 9*b*(b + 2*y0)*Cos[4*Pi*(a + x0)] -
72*b*Pi*x0*(b + 2*y0)*Sin[4*Pi*(-a + x0)] + 72*b*Pi*(a + x0)*(b + 2*y0)*
Sin[4*Pi*(a + x0)] + 12*(a + x0)^2*(6*Pi*y0*Cos[3*Pi*(b - y0)] -
6*Pi*(b + y0)*Cos[3*Pi*(b + y0)] + Sin[3*Pi*(b - y0)] +
Sin[3*Pi*(b + y0)]) - 4*(-a + x0)*(a + 3*x0)*(6*Pi*y0*Cos[3*Pi*(b - y0)] -
6*Pi*(b + y0)*Cos[3*Pi*(b + y0)] + Sin[3*Pi*(b - y0)] +
Sin[3*Pi*(b + y0)])), (1/(96*Pi^2*(a + x0)*(b + y0)))*
(9*b*(b + 2*y0)*Cos[4*Pi*(-a + x0)] - 9*b*(b + 2*y0)*Cos[4*Pi*(a + x0)] +
8*a*(12*a*Pi*y0*Cos[3*Pi*(b - y0)] - 12*a*Pi*(b + y0)*Cos[3*Pi*(b + y0)] -
9*b^2*Pi*Sin[4*Pi*(-a + x0)] - 18*b*Pi*y0*Sin[4*Pi*(-a + x0)] +
2*a*Sin[3*Pi*(b - y0)] + 2*a*Sin[3*Pi*(b + y0)]))}



##### 4.15.2.1 case 2

Case $$f\left ( x,y\right ) =1$$

$I_{i}^{f}=\int \limits _{x=\left ( x_{0}-a\right ) }^{x=x_{0}+a}\int \limits _{y=\left ( y_{0}-b\right ) }^{y=y_{0}+b}\begin {Bmatrix} \frac {1}{A}\left ( \left ( x_{0}-a\right ) -x\right ) \left ( y_{0}+b-y\right ) \\ \frac {1}{A}\left ( \left ( x_{0}-a\right ) +x\right ) \left ( y_{0}+b-y\right ) \\ \frac {1}{A}\left ( \left ( x_{0}-a\right ) +x\right ) \left ( y_{0}+b+y\right ) \\ \frac {1}{A}\left ( \left ( x_{0}-a\right ) -x\right ) \left ( y_{0}+b+y\right ) \end {Bmatrix} \,dxdy$

This was done using Mathematica

$I_{i}^{f}= \begin {Bmatrix} \frac {a^2 b^2}{(a + x_0)(b + y_0)} \\ \frac {a b^2 (a + 2 x_0)}{(a + x_0) (b + y_0)}\\ \frac {a b (a + 2 x_0) (b + 2 y_0)}{(a + x_0)(b + y_0)}\\ \frac {a^2 b (b + 2 y_0)}{(a + x_0) (b + y_0)} \end {Bmatrix}$

phi /. {a :> (x0 + a), b :> (y0 + b)}
Integrate[% , {x, x0 - a, x0 + a}, {y, y0 - b, y0 + b}]
{(a^2*b^2)/((a + x0)*(b + y0)),
(a*b^2*(a + 2*x0))/((a + x0)*(b + y0)),
(a*b*(a + 2*x0)*(b + 2*y0))/((a + x0)*(b + y0)),
(a^2*b*(b + 2*y0))/((a + x0)*(b + y0))}



Now that the local $$k_{i}$$ and local $$f_{i}$$ are calculated, the next step is to assemble them to the global stiﬀness and global force vector. For 1D this process was easy and direct. For higher dimensions we need to make a table of mapping to help in the process. The following diagram shows this mapping table and how it is applied in this example. The mapping table will have as many rows as there are elements, and will have as many columns as the number of nodes per element. Hence in this example it will a matrix of size $$9\times 4$$. Each row gives the global node number for each local node number of that element. For example, for element $$1$$ the ﬁrst row will show the global node number for each of the local nodes for the ﬁrst element.

Local node numbers are always incremented counter clockwise, starting from $$1$$ to $$4$$ in this case.

Now that we have the mapping done, we can now assemble the global stiﬀness and vector matrix and solve for the unknowns.

To evaluate the force vector, we also need a table of physical coordinates of the center of each element relative to the origin of the global coordinates system, since the force vector integration depends on physical coordinates.

The integration was carried out oﬀ-line as shown above, but it was done in terms of physical coordinates of center of the element. The following diagram shows the mapping table of physical coordinates to elements. This uses $$\left ( 0,0\right )$$ as the origin and the coordinates of the lower left corner of the global coordinate system. Hence the coordinate the top right corner of the global grid will be $$\left ( 3h,3h\right )$$ where $$h$$ is the grid spacing.

Now the solution can start. We loop over each element and ﬁll in the global stiﬀness matrix using the mapping above, and also evaluate the global force vector using the physical coordinates.

When the global stiﬀness matrix and the global force vector is ﬁlled, we adjust the global stiﬀness matrix by putting zero in the row numbered $$i$$ that correspond to the edge node $$i$$, and by putting 1 in the diagonal entry $$(i,i)$$. We also put a zero in the force vector for each node on the boundary conditions since we are using Dirichlet boundary conditions which is zero. When this is done, the system $$Ax=f$$ is then solved for $$x$$ where $$x$$ now is the displacement or value of $$y$$ at the nodes. The following is an implementation of the above.

function matlab_97()
close all; clear all;

number_of_elements=[16 64 256 625];
function_string={'1','xy','3(\cos(4\pi x)+\sin(3 \pi y))'};
function_handles={@f_1,@f_xy,@f_2};

figure_image_file={'matlab_e97_1','matlab_e97_2','matlab_e97_3'};

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(a,b)

k11=(a^2+b^2)/(3*a*b);
k12=a/(6*b)-b/(3*a);
k13=-(a^2+b^2)/(6*a*b);
k14=-a/(3*b)+b/(6*a);
k_local = [k11  k12 k13 k14;
k12  k11 k14 k13;
k13  k14 k11 k12;
k14  k13 k12 k11];
end
%----------------------------------------------------
%use this for f(x,y)=x*y
function f_local=f_xy(x0,y0,a,b)

f_local=1/(9*(a + x0)*(b + y0))*...
[a^2*b^2*(a-3*x0)*(b -3*y0);
-a*b^2*(a^2 + 3*a*x0 + 6*x0^2)*(b - 3*y0);
a*b*(a^2 + 3*a*x0 + 6*x0^2)*(b^2 + 3*b*y0 + 6*y0^2);
-a^2*b*(a - 3*x0)*(b^2 + 3*b*y0 + 6*y0^2)...
];
end
%----------------------------------------------------
%use this for f(x,y)=1
function f_local=f_1(x0,y0,a,b)

f_local=[(a^2*b^2)/((a + x0)*(b + y0));
(a*b^2*(a + 2*x0))/((a + x0)*(b + y0));
(a*b*(a + 2*x0)*(b + 2*y0))/((a + x0)*(b + y0));
(a^2*b*(b + 2*y0))/((a + x0)*(b + y0))];

end
%----------------------------------------------------
%use this for f(x,y)= 3(cos(4 pi x)+sin(3 pi y))
function f_local=f_2(x0,y0,a,b)
f_local=[(1/(96*pi^2*(a + x0)*(b + y0)))*...
(9*b^2*cos(4*pi*(-a + x0)) -...
9*b^2*cos(4*pi*(a + x0)) + ...
8*a*(12*a*b*pi*cos(3*pi*(b - y0)) -...
9*b^2*pi*sin(4*pi*(-a + x0)) - ...
2*a*(sin(3*pi*(b - y0)) +...
sin(3*pi*(b + y0)))));
(1/(96*pi^2*(a + x0)*(b + y0)))*...
(-9*b^2*cos(4*pi*(-a + x0)) + 9*b^2*cos(4*pi*(a + x0)) +...
8*(12*a*b*pi*(a + 2*x0)*cos(3*pi*(b - y0)) -...
9*b^2*pi*x0*sin(4*pi*(-a + x0)) + ...
9*a*b^2*pi*sin(4*pi*(a + x0)) +...
9*b^2*pi*x0*sin(4*pi*(a + x0)) - 2*a^2*sin(3*pi*(b - y0)) -...
4*a*x0*sin(3*pi*(b - y0)) - 2*a^2*sin(3*pi*(b + y0)) -...
4*a*x0*sin(3*pi*(b + y0))));
(1/(96*pi^2*(a + x0)*(b + y0)))*...
(-9*b*(b + 2*y0)*cos(4*pi*(-a + x0)) + ...
9*b*(b + 2*y0)*cos(4*pi*(a + x0)) -...
72*b*pi*x0*(b + 2*y0)*sin(4*pi*(-a + x0)) + ...
72*b*pi*(a + x0)*(b + 2*y0)*...
sin(4*pi*(a + x0)) + 12*(a + x0)^2*...
(6*pi*y0*cos(3*pi*(b - y0)) -...
6*pi*(b + y0)*cos(3*pi*(b + y0)) + sin(3*pi*(b - y0)) +...
sin(3*pi*(b + y0))) - 4*(-a + x0)*(a + 3*x0)*...
(6*pi*y0*cos(3*pi*(b - y0)) -...
6*pi*(b + y0)*cos(3*pi*(b + y0)) + sin(3*pi*(b - y0)) +...
sin(3*pi*(b + y0))));
(1/(96*pi^2*(a + x0)*(b + y0)))*...
(9*b*(b + 2*y0)*cos(4*pi*(-a + x0)) - ...
9*b*(b + 2*y0)*cos(4*pi*(a + x0)) +...
8*a*(12*a*pi*y0*cos(3*pi*(b - y0)) - ...
12*a*pi*(b + y0)*cos(3*pi*(b + y0)) -...
9*b^2*pi*sin(4*pi*(-a + x0)) - ...
18*b*pi*y0*sin(4*pi*(-a + x0)) +...
2*a*sin(3*pi*(b - y0)) + 2*a*sin(3*pi*(b + y0))))];

end
%----------------------------------------------------
function solve_poisson_2D_square(fh,rhs,M,plot_file)
%2a is width of element
%2b is height of element
%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
%k=local stiffness matrix
L  = 1; %length of grid
N  = 4; %number of nodes per element
M1 = sqrt(M); %number of elements per row
a=(L/M1)/2;
b=(L/M1)/2;
k_local=getk(a,b);
N1 = M1+1; %number of nodes per edge
dof=4;
nodes = N1^2; %total number of nodes
kk    = zeros(nodes); %global stiffness vector
f     = zeros(nodes,1); %global force vector
mtb   = generate_coordinates_table(M1);
fmtb  = generate_physical_coordinates_table(M1,a,b);

for i = 1:M
x0 = fmtb(i,1); y0 = fmtb(i,2);
%call to load different force vector, pre-computed offline
f_local=fh(x0,y0,a,b);

%assemble force vecotr and global stiffness matrix
for ii = 1:N
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 (M1+2):(M1+1):nodes-N1 ...
(2*M1+2):(M1+1):(nodes-N1) nodes-N1:nodes];
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:2*a:1,0:2*b: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', '-r300', ['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 mtb=generate_coordinates_table(M)
%M=number_of_elements_per_row
mtb=zeros(M^2,4);
for i=1:M
for j=1:M
mtb(j+M*(i-1),:) =[j j+1 j+(M+2) j+(M+1)]+(M+1)*(i-1);
end
end
end
%----------------------------------------------------
function fmtb=generate_physical_coordinates_table(M,a,b)
%2a is width of element
%2b is height of element
%M=number_of_elements_per_row
fmtb=zeros(M^2,2);
for i=1:M
for j=1:M
fmtb(j+M*(i-1),:) = [(2*j-1)*a (2*i-1)*b];
end
end
end


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

#### 4.15.5 case 3

Case 3: $$f(x,y)=3\left ( \cos (4 \pi x)+\sin (3 \pi y)\right )$$

#### 4.15.6 Solving for reactangle grid

A small modiﬁcation is now made to allow one to specify diﬀerent width and height for the domain.

function matlab_97_2()
close all; clear all;

%number of elements per unit length, width, height
number_of_elements=[4 1 1;
8 1 2;
16 1 4;
25 1 6];
function_string={'1','xy','3(\cos(4\pi x)+\sin(3 \pi y))'};
function_handles={@f_1,@f_xy,@f_2};

figure_image_file={'matlab_e97_rectangle_1',...
'matlab_e97_rectangle_2','matlab_e97_rectangle_3'};

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

%----------------------------------
function k_local=getk(a,b)

k11=(a^2+b^2)/(3*a*b);
k12=a/(6*b)-b/(3*a);
k13=-(a^2+b^2)/(6*a*b);
k14=-a/(3*b)+b/(6*a);
k_local = [k11  k12 k13 k14;
k12  k11 k14 k13;
k13  k14 k11 k12;
k14  k13 k12 k11];
end
%----------------------------------------------------
%use this for f(x,y)=x*y
function f_local=f_xy(x0,y0,a,b)

f_local=1/(9*(a + x0)*(b + y0))*...
[a^2*b^2*(a-3*x0)*(b -3*y0);
-a*b^2*(a^2 + 3*a*x0 + 6*x0^2)*(b - 3*y0);
a*b*(a^2 + 3*a*x0 + 6*x0^2)*(b^2 + 3*b*y0 + 6*y0^2);
-a^2*b*(a - 3*x0)*(b^2 + 3*b*y0 + 6*y0^2)...
];
end
%----------------------------------------------------
%use this for f(x,y)=1
function f_local=f_1(x0,y0,a,b)

f_local=[(a^2*b^2)/((a + x0)*(b + y0));
(a*b^2*(a + 2*x0))/((a + x0)*(b + y0));
(a*b*(a + 2*x0)*(b + 2*y0))/((a + x0)*(b + y0));
(a^2*b*(b + 2*y0))/((a + x0)*(b + y0))];

end
%----------------------------------------------------
%use this for f(x,y)= 3(cos(4 pi x)+sin(3 pi y))
function f_local=f_2(x0,y0,a,b)
f_local=[(1/(96*pi^2*(a + x0)*(b + y0)))*...
(9*b^2*cos(4*pi*(-a + x0)) -...
9*b^2*cos(4*pi*(a + x0)) + 8*a*(12*a*b*pi*cos(3*pi*(b - y0)) -...
9*b^2*pi*sin(4*pi*(-a + x0)) - 2*a*(sin(3*pi*(b - y0)) +...
sin(3*pi*(b + y0)))));
(1/(96*pi^2*(a + x0)*(b + y0)))*...
(-9*b^2*cos(4*pi*(-a + x0)) + 9*b^2*cos(4*pi*(a + x0)) +...
8*(12*a*b*pi*(a + 2*x0)*cos(3*pi*(b - y0)) -...
9*b^2*pi*x0*sin(4*pi*(-a + x0)) + ...
9*a*b^2*pi*sin(4*pi*(a + x0)) +...
9*b^2*pi*x0*sin(4*pi*(a + x0)) - 2*a^2*sin(3*pi*(b - y0)) -...
4*a*x0*sin(3*pi*(b - y0)) - 2*a^2*sin(3*pi*(b + y0)) -...
4*a*x0*sin(3*pi*(b + y0))));
(1/(96*pi^2*(a + x0)*(b + y0)))*...
(-9*b*(b + 2*y0)*cos(4*pi*(-a + x0)) + 9*b*(b + 2*y0)*...
cos(4*pi*(a + x0)) -...
72*b*pi*x0*(b + 2*y0)*sin(4*pi*(-a + x0)) + ...
72*b*pi*(a + x0)*(b + 2*y0)*...
sin(4*pi*(a + x0)) + 12*(a + x0)^2*(6*pi*y0*...
cos(3*pi*(b - y0)) -...
6*pi*(b + y0)*cos(3*pi*(b + y0)) + ...
sin(3*pi*(b - y0)) +...
sin(3*pi*(b + y0))) - 4*(-a + x0)*(a + 3*x0)*...
(6*pi*y0*cos(3*pi*(b - y0)) -...
6*pi*(b + y0)*cos(3*pi*(b + y0)) + ...
sin(3*pi*(b - y0)) +...
sin(3*pi*(b + y0))));
(1/(96*pi^2*(a + x0)*(b + y0)))*...
(9*b*(b + 2*y0)*cos(4*pi*(-a + x0)) - 9*b*(b + 2*y0)*...
cos(4*pi*(a + x0)) +...
8*a*(12*a*pi*y0*cos(3*pi*(b - y0)) - 12*a*pi*(b + y0)*...
cos(3*pi*(b + y0)) -...
9*b^2*pi*sin(4*pi*(-a + x0)) - 18*b*pi*y0*sin(4*pi*(-a + x0)) +...
2*a*sin(3*pi*(b - y0)) + 2*a*sin(3*pi*(b + y0))))];

end
%----------------------------------------------------
function solve_poisson_2D_square(fh,rhs,M,Lx,Ly,plot_file)
%fh is function which evaluate the force vector
%rhs is string which is the f(x,y), for plotting only
%M=number of element per unit length
%Lx=length in x-direction
%Ly=length in y-direction
N  = 4; %number of nodes per element
Mx=M*Lx;
My=M*Ly;
a=(1/M)/2;
b=a;
k_local=getk(a,b);
dof=4;

nodes = (Mx+1)*(My+1); %total number of nodes
kk    = zeros(nodes); %global stiffness vector
f     = zeros(nodes,1); %global force vector
mtb   = generate_coordinates_table(Mx,My);
fmtb  = generate_physical_coordinates_table(Mx,My,a,b);

for i = 1:(Mx*My)
x0 = fmtb(i,1); y0 = fmtb(i,2);
%call to load different force vector, pre-computed offline
f_local=fh(x0,y0,a,b);

%assemble force vecotr and global stiffness matrix
for ii = 1:N
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:Mx+1 (Mx+2):(Mx+1):nodes-...
Mx (2*Mx+2):(Mx+1):(nodes-(Mx+1)) nodes-(Mx-1):nodes];
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,Mx+1,My+1);
[X,Y] = meshgrid(0:2*b:Ly,0:2*a:Lx);
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,Mx*My);
title(plot_title,'fontweight','bold','fontsize',14);

print(gcf, '-dpdf', '-r300', ['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 mtb=generate_coordinates_table(Mx,My)
%Mx= number of element x-wise
%My= number of element y-wise
mtb=zeros(Mx*My,4);
for i=1:My
for j=1:Mx
mtb(j+Mx*(i-1),:) =[j j+1 j+(Mx+2) j+(Mx+1)]+(Mx+1)*(i-1);
end
end
end
%----------------------------------------------------
function fmtb=generate_physical_coordinates_table(Mx,My,a,b)
%Mx= number of element x-wise
%My= number of element y-wise
%2a is width of element
%2b is height of element
fmtb=zeros(Mx*My,2);
for i=1:My
for j=1:Mx
fmtb(j+Mx*(i-1),:) = [(2*j-1)*a (2*i-1)*b];
end
end
end


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

#### 4.15.7 Case 4

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

#### 4.15.8 Case 5

Case : $$f(x,y)=xy$$

#### 4.15.9 case 6

Case $$f(x,y)=3\left ( \cos (4 \pi x)+\sin (3 \pi y)\right )$$.