Variational integrators are geometric structure-preserving numerical methods that preserve the symplectic structure, satisfy a discrete Noether's theorem, and exhibit exhibit excellent long-time energy stability properties. An exact discrete Lagrangian arises from Jacobi's solution of the Hamilton-Jacobi equation, and it generates the exact flow of a Lagrangian system. By approximating the exact discrete Lagrangian using an appropriate choice of interpolation space and quadrature rule, we obtain a systematic approach for constructing variational integrators. The convergence rates of such variational integrators are related to the best approximation properties of the interpolation space. Many gauge field theories can be formulated variationally using a multisymplectic Lagrangian formulation, and we will present a characterization of the exact generating functionals that generate the multisymplectic relation. By discretizing these using group-equivariant spacetime finite element spaces, we obtain methods that exhibit a discrete multimomentum conservation law. We will then briefly describe an approach for constructing group-equivariant interpolation spaces that take values in the space of Lorentzian metrics that can be efficiently computed using a generalized polar decomposition. The goal is to eventually apply this to the construction of variational discretizations of general relativity, which is a second-order gauge field theory whose configuration manifold is the space of Lorentzian metrics.