Virtual Quake is a boundary element code that performs simulations of fault systems based on stress interactions between fault elements to understand long term statistical behavior. It uses field observations to define fault topology, long-term slip rates and frictional parameters. The faults are meshed into interacting elements and quasi-static elastic interactions are calculated between these elements. Stress is then applied to each element at geologically-observed rates until frictional parameters are exceeded. At this point an element will fail and transfer stress to the rest of the system. Under the correct conditions, transferred stress results in propagating ruptures throughout the system, i.e. a simulated earthquake. The design of Virtual Quake allows for fast execution so many thousands of events can be generated over very long simulated time periods. The result is a rich dataset from which to study the statistical properties of the rupturing fault system. Additionally, many data visualization and analysis tools are provided in the PyVQ python script. v3.1.1 A mesher bug has been fixed that was causing incorrect interpolation of fault properties between trace points. v3.1.0 Improved mesher to better handle faults with a curved trace and shallow dip. v3.0.0 This update adds a great deal of new functionality and stability: The rupture model has been overhauled for stability. Event slip matrix solutions have been removed in favor of a purely cellular automaton method. Faults are now separate objects from sections. This allows ruptures to spread between sections belonging to the same larger fault, while allowing each section to have its properties defined independently. Many new PyVQ plotting and filtering options have been added. Several major and minor bugs have been fixed.