A framework of the numerical solutions of the Landau-Lifshitz-Gilbert equation and a corresponding library are developed in this paper. The numerical framework is built on the tetrahedron meshes, consisting of the finite element method for the spatial discretization and the implicit midpoint scheme for the temporal discretization. The computational complexity on calculating the demagnetization field is effectively reduced by using a PDE approach, in which a gradient recovery technique is used for preserving the numerical accuracy. A C++ library is realized for the numerical simulations, and OpenMP is adopted for accelerating the simulations. The numerical convergence of the proposed method is studied in detail for the μMAG standard problem #3, in which a limit value for the side length of the cube is predicted. The ability of the proposed method on handling problems defined on the complex domain is successfully demonstrated by the quality numerical solutions from several examples, in which the computational domains are the thin films with irregular defects.