# Finite Element Method - 2D Mesh Generator - Metfem2D

## Finite Element Method. Mathematical Basis.

### Integration upon Elements.

We shall consider the linear expression derived in the Appendix B[1, 2, 3]

 δ φA T (KφA + f) = ∫ δφ ( -ε ∂2/∂x2 - ε∂2/∂y2) φ(x, y) dxdy + ∫ δφρ(x, y) dxdy = 0 Ω Ω

with the boundary condition φ = γ on Γ. In such a simply case of integral – differential problems with a differential operator L = -ε∂2/∂x2 - ε∂2/∂y2, the variable φ in the above – written equation only consists of one scalar function φ which is the sought solution, while the constant vector f is represented by the last term in that expression. To find the solution for such a problem means to determine the values of φ(x,y) in the whole domain Ω. The values of φ on its boundary Γ are already prescribed to γ. On the other hand, at the very beginning (see Eq. (2) in the previous page) we have postulated that a function φ could be approximated by an expansion φA given by means of some basis functions Nm(x,y), m = 1, …,n (for more details see Appendix A). Thus another possibility to deal with the Poisson problem is just to start from the functional Π and build a set of Euler equations ∂Π/∂φm = 0 where m = 1, …, n and φm approximates value of the solution φ calculated at the m-th mesh node.

 Π = ∫ [1/2ε ∑l (∂Nl/∂x φl)2 + 1/2ε ∑l (∂Nl/∂y φl)2 + ρ ∑l Nl φl ] dxdy + Ω
 + ∫ (γ - 1/2 ∑l Nlφl) ∑k Nk φk dΓ Γ

and after that we calculate the derivative ∂Π/∂φA,m. Moreover, let's simplify our problem by neglecting the last term in the above – presented equation and imposing φ = γ on the boundary Γ instead. In that manner, one obtains the expression

 ∂Π/∂φA,m = ∫ [ε ∑l (∂Nl/∂x φl) ∑k ∂ Nk/∂x δkm Ω
 + ε ∑l (∂ Nl/∂y φl) ∑k ∂Nk/∂y δkm + ρ ∑l Nl δlm] dxdy = 0

or in a simplified form

 ∂Π/∂φA, m = ∑l ( ∫ ε ( ∂Nl/∂x ∂Nm/∂x + ∂Nl/∂y ∂Nm/∂y ) dxdy ) φl + ∫ ρNm dxdy = 0. Ω Ω

It is worth mentioning that some requirements must be imposed on the shape functions N. Namely, if n-th order derivatives occur in any term ofL then the shape functions have to be such that their n-1 derivatives (pay an attention to the above – presented equation) are continuous and finite. Therefore, generally speaking Cn-1 continuity of shape functions must be preserve.
In turn, after substituting

 Kml = ε ∫ ( ∂Nl/∂x ∂Nm/∂x + ∂Nl/∂y ∂Nm/∂y )dxdy Ω fm = ∫ ρNm dxdy Ω

finally one obtains a set of equations

 ∂Π/∂ φm = ∑l Kml φl + fm = 0
for m, l = 1, …, n or in matrix description
∂Π/∂φ = K φ + f = 0.

It is worth noticing that the matrix K is a symmetric one because of the symmetry in exchange of subscripts l and m in the equation.
Now, we are obliged to employ the division of our domain Ω into a set of subdomains Ωi. It gives that

 Kml = ∑i Kiml = ∑i ∫ ε(x,y) ( ∂Nil(x,y)/∂x ∂Nim(x,y)/∂x + ∂Nil(x,y)/∂y ∂Nim(x,y)/∂y) dxdy Ωi fm = ∑i fim = ∑i ∫ ρ(x,y) Nim(x,y) dxdy. Ωi

Therefore, after the transformation to I subdomains the expression becomes

 ∑i Kiml φl + ∑i fim = 0

for i = 1, …, I and m = 1, …, n. Or in matrix notation

 φA = -K-1 f

In fact, the summation written above takes into account only these elements Ωi which contribute to m-th node, however, because of the consistency in notation all elements are included in the sum with the exception that those Nim functions for which node m does not occur in i-th element are put equal zero.
From now, the whole story is to calculate integrals

 Kiml = ∫ ε(x, y) (∂Nil (x, y)/∂x ∂Nim (x, y) /∂x + ∂Nil (x, y) /∂y ∂Nim (x, y)/∂y)dxdy Ωi
 1 1-L1 = ∫ dL1 ∫ dL2 ε(L1, L2, L3) |det Ji| (∇ Nil T TT ∇T Nim), 0 0
 1 1-L1 fim = ∫ ρ(x,y) Nim (x,y) dxdy = ∫ dL1 ∫ dL2 |det Ji| ρ(L1,L2,L3) Nim(L1, L2, L3) Ωi 0 0

where Nim = Lim, L3 = 1 - L1 - L2 (see Appendix A) whereas det(Ji) – the jacobian of i-th element, T matrix together with ∇ operator in new coordinates are evaluated in Appendix C.
An integration over the i-th subdomain Ωi, which is a triangular element with three nodes, enforces the transformation from n-dimensional global interpolation to the local interpolations given by means of Nik(x,y) functions where ik = 1, 2, 3. That is why in equations new indices il,im appear which further are allowed to take three possible values 1,2 and 3 for each element i (the local subspace).
As the next step, the Gauss quadrature is employed to compute above-written integrals numerically as it is described in the Appendix D.
And finally, after incorporating boundary conditions to equations by inserting appropriate boundary values of φ, the system of equations can be solved.

### References

1. ^ O. C. Zienkiewicz, R. L. Taylor and J. Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, Sixth edition., Elsevier 2005
2. ^ A. Kendall, H. Weimin, Theoretical Numerical analysis, A Functional Analysis Framework, Third Edition., Springer 2009
3. ^ R. Courant, D. Hilbert, Methods of Mathematical Physics, Volume 1, Interscience Publisher, New York, 1953