  The static offsets caused by earthquakes are well described by elastostatic
models with a discontinuity in the displacement along the fault. A traditional
approach to model this discontinuity is to align the numerical mesh with the
fault and solve the equations using finite elements. However, this distorted
mesh can be difficult to generate and update. We present a new numerical
method, inspired by the Immersed Interface Method, for solving the elastostatic
equations with embedded discontinuities. This method has been carefully
designed so that it can be used on parallel machines on an adapted finite
difference grid. We have implemented this method in Gamra, a new code for earth
modelling. We demonstrate the correctness of the method with analytic tests,
and we demonstrate its practical performance by solving a realistic earthquake
model to extremely high precision.
