The finite elements method (FEM) is based on the idea of division the whole domain Ω into a number of finite sized elements or subdomains Ωi in order to approximate a continuum problem by a behavior of an equivalent assembly of discrete finite elements .
In the presence of a set of elements Ωi the total integral over the domain Ω is represented by the sum of integrals over individual subdomains Ωi
|∫||L (u, ∂u/∂x, …)dΩ =||∑i||∫||L(u, ∂u/∂x, …) dΩi|
|∫||L(u, ∂u/∂x, …) dΓ =||∑i||∫||L(u, ∂u/∂x, …)dΓi,|
where L(u, ∂u/∂x, …) denotes a differential operator.
The continuum problem is posed by appropriate differential equations (e. g. Laplace or Poisson equation) and boundary conditions that are imposed on the unknown solution φ. The general procedure of FEM is aimed at finding an approximate solution φA given by the expansion:
|φ ∼ φA =||∑ j||φAj Nj||=||φAj Nj,|
where Nj (j = 1, …, n) are shape functions (basis functions or interpolation functions) [1, 2] and all or the most of the parameters φAj remain unknown. After dividing the domain Ω, the shape functions are defined locally for elements Ωi. A typical finite element is triangular in shape and thus has three main nodes. It is easy to demonstrate that triangular subdomains fit better to the boundary Γ than others e. g. rectangular ones. Among the triangular elements family one can find linear, quadratic and cubic elements  (see also Appendix A). A choice of an appropriate type of subdomains depends on a desired order of approximation and thus arises directly from the continuum problem. The higher order of element, the better approximation. Each triangular element can be described in terms of its area coordinates Li1, Li2 and Li3. There are general rules that govern the transformation from area coordinates to cartesian coordinates
x = L1x1 + L2x2 + L3x3
y = L1y1 + L2y2 + L3y3
1 = L1 + L2 + L3
where set of pairs (x1, y1), (x2,y2), (x3, y3) represents cartesian nodal coordinates. In turn the area coordinates are related to shape functions in a manner that depends on the element order. In further analysis only the linear triangular elements will be used. For them, the shape functions are simply the area coordinates (see Appendix A). Therefore, each pair of shape functions Nik(x,y), Nil(x,y) for k,l=1,2,3 could be thought as a natural basis of the Ωi triangular element.