The paper discusses the problems of a finite element method programming. Modern C++ functional programming and template metaprogramming approach to finite element analysis is presented. This approach simplifies the implementation of an effective assembly of the stiffness matrix for a problem defined by a weak form. This method is tested by a solution of the Poisson equation on an unstructured 3D tetrahedral mesh using FEM C++ library FEMEngine developed by authors. The function which calculates the matrix is generated by higher order functions during the compilation stage. The performance of the computation is analyzed by studying of a disassembled code and by comparison with the popular open source FEM software.