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 different 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 stiffness 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 first 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 stiffness matrix is the first integral, and the force vector comes from the second integral.
The integration is now carried out to obtain the element stiffness matrix. This was done using Mathematica. The local stiffness 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 finding the local stiffness 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.
This is normally done using Gaussian quadrature method. In this example, the integral is found off-line using Mathematica.
Evaluating the force vector requires integration over the element. This was done off-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.
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)))}
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)]))}
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 stiffness 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 first row will show the global node number for each of the local nodes for the first 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 stiffness 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 off-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 fill in the global stiffness matrix using the mapping above, and also evaluate the global force vector using the physical coordinates.
When the global stiffness matrix and the global force vector is filled, we adjust the global stiffness 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); 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', '-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 different number of elements.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Case 3: \(f(x,y)=3\left ( \cos (4 \pi x)+\sin (3 \pi y)\right ) \)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
A small modification is now made to allow one to specify different 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); shading interp 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 different number of elements.
Case \(f(x,y)=1\).
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Case : \(f(x,y)=xy\)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Case \(f(x,y)=3\left ( \cos (4 \pi x)+\sin (3 \pi y)\right ) \).
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
