We discuss the creation and implementation of a generalized library, named ELEMENTS, of mathematical functions for supporting a very broad range of element types including: linear, quadratic, and cubic serendipity elements in 2D and 3D; high-order spectral elements; and a linear 4D element. The ELEMENTS library can be used for research and development of both continuous and discontinuous finite element methods for solving a diverse range of partial differential equations. The library has functions for calculating quantities that are commonly used in finite element methods such as the gradient of a basis function, the Jacobi matrix, the inverse Jacobi matrix, the determinant of the Jacobi matrix, and a physical position inside the element, to name a few examples. The library also supports both Gauss–Legendre and Gauss–Lobatto quadrature rules up to 8 quadrature points in each coordinate direction. The examples and discussions in this paper will focus on Lagrangian solid mechanics and dynamics, but ELEMENTS can be used for many other applications.