Anomalous X-ray pulsars (AXPs) and soft gamma-ray repeaters (SGRs) are recognized as the most promising magnetar candidates, as indicated by their energetic bursts and rapid spin-downs. It is expected that the strong magnetic field leaves distinctive imprints on the emergent radiation both by affecting the radiative processes in atmospheres of magnetars and by scattering in the upper magnetospheres. We construct a self-consistent physical model that incorporates emission from the magnetar surface and its reprocessing in the three-dimensional (3D) twisted magnetosphere using a Monte Carlo technique. The synthetic spectra are characterized by four parameters: surface temperature kT, surface magnetic field strength B, magnetospheric twist angle.., and the normalized electron velocity beta. We also create a tabular model (STEMS3D) and apply it to a large sample of XMM-Newton spectra of magnetars. The model successfully fits nearly all spectra, and the obtained magnetic field for 7 out of the 11 sources are consistent with the values inferred from the spin-down rates. We conclude that the continuum-fitting by our model is a robust method to measure the magnetic field strength and magnetospheric configuration of AXPs and SGRs. Investigating the multiple observations of variable sources, we also study the mechanism of their spectral evolution. Our results suggest that the magnetospheres in these sources are highly twisted (Delta phi > 1), and the behavior of magnetospheric twisting and untwisting is revealed in the 2002 outburst of 1E 2259+586.