In this paper, a finite element reconstruction algorithm for three-dimensional photoacoustic tomography is described. The algorithm is based on rigorous iterative solution to the Helmholtz photoacoustic wave equation coupled with regularization techniques and is able to recover both the images of absorbed optical energy density and acoustic speed simultaneously. The algorithm is tested using various numerical examples that mimic cancer detection and joint imaging. The results show that the algorithm is able to reconstruct photoacoustic images quantitatively in terms of the location, size, optical and acoustic properties of the target, and background media for various examples examined.