Numerical simulations play a key role in the study of systems where gravity is strong enough that it must be treated relativistically. In the first portion of this work, I apply a recent numerical method, the discontinuous Galerkin (DG) method, to improve the accuracy of simulations in which matter couples to strong gravity. I use the DG method to simultaneously evolve both the spacetime geometry and the matter on the same computational grid, and I deform this grid to conform to the problem symmetries. I show results for 3D evolutions of a Kerr black hole, a neutron star in which the spacetime metric is held fixed, and, finally, a neutron star where the spacetime and matter are both dynamical. These results mark the first application of the DG method to simultaneous evolution of the spacetime geometry and matter. The evolutions show long-term stability, good accuracy, and an improved rate of convergence as compared to evolutions with a comparable-resolution finite volume method. In the second portion of this work, I contribute to the development of a new visualization technique for systems with strong gravitational fields. We present a method of calculating the strong-field gravitational lensing caused by many analytic and numerical spacetimes. We then use this method to simulate the visual distortions near isolated black holes and black hole binary systems; we produce both demonstrative images that illustrate details of the spatial distortion, and realistic images of collections of stars, taking both lensing amplification and redshift into account.