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}/∂x^{2}  ε∂^{2}/∂y^{2})
φ(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}/∂x^{2}  ε∂^{2}/∂y^{2}, 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 N_{m}(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 mth mesh node.
Π

= 
∫

[1/2ε

∑_{l}

(∂N_{l}/∂x φ^{l})^{2} + 1/2ε

∑_{l}

(∂N_{l}/∂y φ^{l})^{2}
+ ρ

∑_{l}

N_{l} φ^{l}

] 
dxdy +


Ω

+ 
∫

(γ  1/2

∑_{l}

N_{l}φ^{l})

∑_{k}

N_{k} φ^{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}

(∂N_{l}/∂x φ^{l})

∑_{k}

∂ N_{k}/∂x δ^{k}_{m}


Ω


+ ε

∑_{l}

(∂ N_{l}/∂y φ^{l})

∑_{k}

∂N_{k}/∂y δ^{k}_{m}
+ ρ

∑_{l}

N_{l} δ^{l}_{m}] dxdy = 0

or in a simplified form
∂Π/∂φ_{A, m}
=

∑_{l}

(

∫

ε
(
∂N_{l}/∂x ∂N^{m}/∂x
+
∂N_{l}/∂y ∂N^{m}/∂y
)
dxdy
)
φ^{l} +

∫

ρN^{m} dxdy = 0.


Ω

Ω

It is worth mentioning that some requirements must be imposed on the shape functions N. Namely, if nth order derivatives occur in any term ofL then the shape functions have to be such that their n1 derivatives (pay an attention to the above – presented equation) are continuous and finite. Therefore, generally speaking C_{n1} continuity of shape functions must be preserve.
In turn, after substituting
K^{m}_{l} =
ε

∫

(
∂N_{l}/∂x ∂N^{m}/∂x
+
∂N_{l}/∂y ∂N^{m}/∂y
)dxdy


Ω

f^{m} =

∫

ρN^{m} dxdy


Ω

finally one obtains a set of equations
∂Π/∂ φ_{m}
=

∑_{l}

K^{m}_{l}
φ^{l} + f^{m}
= 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
K^{m}_{l}
=

∑_{i}

K^{im}_{l}
=

∑_{i}

∫

ε(x,y) (
∂N^{i}_{l}(x,y)/∂x
∂N^{im}(x,y)/∂x +
∂N^{i}_{l}(x,y)/∂y
∂N^{im}(x,y)/∂y) dxdy


Ω^{i}

f^{m} =

∑_{i}

f^{im} =

∑_{i}

∫

ρ(x,y) N^{im}(x,y) dxdy.


Ω^{i}

Therefore, after the transformation to I subdomains the expression becomes
∑_{i}

K^{im}_{l} φ^{l} +

∑_{i}

f^{im} = 0

for i = 1, …, I and m = 1, …, n. Or in matrix notation
In fact, the summation written above takes into account only these elements Ω^{i} which contribute to mth node, however, because of the consistency in notation all elements are included in the sum with the exception that those N^{i}_{m} functions for which node m does not occur in ith element are put equal zero.
From now, the whole story is to calculate integrals
K^{im}_{l}

=

∫

ε(x, y)

(∂N^{i}_{l}

(x, y)/∂x
∂N^{im}

(x, y)

/∂x +

∂N^{i}_{l}

(x, y)

/∂y ∂N^{im} (x, y)/∂y)dxdy


Ω^{i}


1

1L_{1}


=

∫

dL_{1}

∫

dL_{2} ε(L_{1}, L_{2},

L_{3})

det J^{i}
(∇ N_{il}
T T^{T} ∇^{T} N^{im}),


0

0


1

1L_{1}

f^{im}

=

∫

ρ(x,y)

N^{im}

(x,y) dxdy =

∫

dL_{1}

∫

dL_{2}

det J^{i} ρ(L_{1},L_{2},L_{3}) N^{im}(L_{1}, L_{2}, L_{3})


Ω^{i}

0

0

where N^{im} = L^{im}, L_{3} = 1  L_{1}  L_{2} (see Appendix A) whereas det(J^{i}) – the jacobian of ith element, T matrix together with ∇ operator in new coordinates are evaluated in Appendix C.
An integration over the ith subdomain Ω^{i}, which is a triangular element with three nodes, enforces the transformation from ndimensional global interpolation to the local interpolations given by means of N_{ik}(x,y) functions where i_{k} = 1, 2, 3. That is why in equations new indices i_{l},i_{m} 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 abovewritten 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