A novel method for calculating the Nusselt number, Nu, in a steady-state Rayleigh-Benard convection problem is presented, in which calculations are done for a square box with constant temperature, free-slip boundary conditions at the top and bottom, and a reflection symmetry along the side walls. The element heat flux is obtained by averaging over the entire element; element heat fluxes are then projected to the adjacent nodes. Compared with previous methods, the approach reduces the calculated depth variation in horizontally averaged flux by more than a factor of 10 and shows more rapid convergence of Nu as a function of grid size.