  We have extended Cosmos++, a multi-dimensional unstructured adaptive mesh
code for solving the covariant Newtonian and general relativistic radiation
magnetohydrodynamic (MHD) equations, to accommodate both discrete finite volume
and arbitrarily high order finite element structures. The new finite element
implementation, called CosmosDG, is based on a discontinuous Galerkin (DG)
formulation, using both entropy-based artificial viscosity and slope limiting
procedures for regularization of shocks. High order multi-stage forward Euler
and strong stability preserving Runge-Kutta time integration options complement
high order spatial discretization. We have also added flexibility in the code
infrastructure allowing for both adaptive mesh and adaptive basis order
refinement to be performed separately or simultaneously in a local
(cell-by-cell) manner. We discuss in this report the DG formulation and present
tests demonstrating the robustness, accuracy, and convergence of our numerical
methods applied to special and general relativistic MHD, though we note an
equivalent capability currently also exists in CosmosDG for Newtonian systems.
