In this talk we present the results of a comparison numerical study [1] of the efficiency of different finite element formulations for elliptic problems. Namely, they are: a) the classical H¹-conforming Galerkin, and b) the Mixed Formulations, c) the Symmetric Interior Penalty DG method, the hybrid formulations d) Stabilized Hybrid Discontinuous Galerkin method, and  e)  Locally Discontinuous but Globally Continuous method.  Excepting the Symmetric Interior Penalty DG method, all the other cases allow static condensation, reducing the size of the global system of equations to be solved, and consequently the execution time. The different formulations are implemented by the NeoPZ software and compared in terms of the L²-norm of the approximation errors, number of degrees-of-freedom, with and without static condensation, and CPU times required for the simulations. A brief description of the general structure of NeoPZ will be presented by emphasizing its essential features.  For a 3D test problem with smooth solution, the simulations are performed with h refinement, and constant polynomial degree. For a singular 2D problem, the results are for approximation spaces based on a given set of hp-refined meshes adapted to the internal layer.  Finally, we will present how triangular tensor spaces have been incorporated to NeoPZ for approximations of two dimensional elasticity problems.