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 ^{[1]}.
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}  
Ω  Ω^{i}  
∫  L(u, ∂u/∂x, …) dΓ =  ∑_{i}  ∫  L(u, ∂u/∂x, …)dΓ^{i},  
Γ  Γ^{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}  φ_{A}^{j} N_{j}  =  φ_{A}^{j} N_{j}, 
where N_{j} (j = 1, …, n) are shape functions (basis functions or interpolation functions) ^{[1}^{, 2]} and all or the most of the parameters φ_{A}^{j} 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 ^{[1]} (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 L^{i}_{1}, L^{i}_{2} and L^{i}_{3}. There are general rules that govern the transformation from area coordinates to cartesian coordinates
x = L_{1}x_{1} + L_{2}x_{2} + L_{3}x_{3} y = L_{1}y_{1} + L_{2}y_{2} + L_{3}y_{3} 1 = L_{1} + L_{2} + L_{3} 
→  ( 
L_{1} L_{2} L_{3} 
)  =  ( 

)  1  ( 
x y 1 
) 
where set of pairs (x_{1}, y_{1}), (x_{2},y_{2}), (x_{3}, y_{3}) 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 N^{i}_{k}(x,y), N^{i}_{l}(x,y) for k,l=1,2,3 could be thought as a natural basis of the Ω^{i} triangular element.