The finite element solution of the hyperbolic heat conduction equations is addressed. The governing system of equations is solved for the temperature and heat fluxes as independent variables. A standard Galerkin method is used for the spatial discretization and a Crank-Nicolson method is adopted for marching in the time domain. It is shown that the proposed method can easily evaluate the entropy production within the domain and assess the thermodynamic equilibrium of the system. The performance of the proposed algorithm is verified by solving a 1D test case. A 2D test case is also studied and some interesting features of the hyperbolic heat conduction are demonstrated. Copyright (C) 1999 John Wiley & Sons, Ltd.