A three-dimensional, incompressible, Rayleigh damped magnetic resonance elastography (MRE) material property reconstruction algorithm capable of reconstructing the spatial distribution of both the real and imaginary parts of the shear modulus, density and bulk modulus from full-field MR-detected harmonic motion data was developed. The algorithm uses a subzone-based implementation of motion error minimization techniques, using 27 hexahedral finite elements, and is written in FORTRAN to run on high performance distributed computing systems. The theory behind the methods used is presented in a form that is directly applicable to the code's structure, to serve as a reference for future research building on this algorithm. Globally defined Rayleigh damping parameter reconstructions using simulated data showed that it is possible to reconstruct the correct combination of Rayleigh parameters under noise levels comparable to MR measurements. The elastic wave equation is used to demonstrate that use of a one parameter damping model to fit a Rayleigh damped material can lead to artefacts in the reconstructed damping parameter images, a prediction that is verified using simulated reconstructions. Initial results using MR-detected motion data from both gelatine phantoms and in-vivo cases produced good reconstructions of real shear modulus, as well as showing promise for successful imaging of damping properties. An initial investigation into an alternative elemental basis function approach to supporting the material property distribution produced some promising results, as well as highlighting some significant issues with large variations across the elements.