Discontinuous Galerkin Finite Element (DGFE) methods offer a mathematically
beautiful, computationally efficient, and efficiently parallelizable way to
solve hyperbolic PDEs. These properties make them highly desirable for
numerical calculations in relativistic astrophysics and many other fields. The
BSSN formulation of the Einstein equations has repeatedly demonstrated its
robustness. The formulation is not only stable but allows for puncture-type
evolutions of black hole systems. To-date no one has been able to solve the
full (3+1)-dimensional BSSN equations using DGFE methods. This is partly
because DGFE discretization often occurs at the level of the equations, not the
derivative operator, and partly because DGFE methods are traditionally
formulated for manifestly flux-conservative systems. By discretizing the
derivative operator, we generalize a particular flavor of DGFE methods, Local
DG methods, to solve arbitrary second-order hyperbolic equations. Because we
discretize at the level of the derivative operator, our method can be
interpreted as either a DGFE method or as a finite differences stencil with
non-constant coefficients.