A framework for the numerical solution of the Landau-Lifshitz-Gilbert equation is developed in this paper. The numerical framework is based on the finite element method on tetrahedral meshes for the spatial discretization and the implicit midpoint scheme for the temporal discretization. The computational complexity for 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. The numerical convergence of the proposed method is studied in detail for the μMAG standard problem #3, from which a limit is predicted for the desired side length. The capability of the proposed method on handling problems defined on complex domains is successfully demonstrated by several examples, in which the computational domains are thin films with irregular defects.