We introduce a new family of discontinuous Galerkin (DG) finite element schemes for the discretization of first order systems of hyperbolic partial differential equations (PDE) on unstructured simplex meshes in two and three space dimensions that respect the two basic vector calculus identities exactly also at the discrete level, namely that the curl of the gradient is zero and that the divergence of the curl is zero. The key ingredient here is the construction of two compatible discrete nabla operators, a primary one and a dual one, both defined on general unstructured simplex meshes in multiple space dimensions. Our new schemes extend existing cell-centered finite volume methods based on corner fluxes to arbitrary high order of accuracy in space. An important feature of our new method is the fact that only two different discrete function spaces are needed to represent the numerical solution, and the choice of the appropriate function space for each variable is related to the origin and nature of the underlying PDE. The first class of variables is discretized at the aid of a discontinuous Galerkin approach, where the numerical solution is represented via piecewise polynomials of degree N and which are allowed to jump across element interfaces. This set of variables is related to those PDE which are mere consequences of the definitions, derived from some abstract scalar and vector potentials, and for which involutions like the divergence-free or the curl-free property must hold if satisfied by the initial data. The second class of variables is discretized via classical continuous Lagrange finite elements of approximation degree M=N+1 and is related to those PDE which can be derived as the Euler-Lagrange equations of an underlying variational principle. The primary nabla operator takes as input the data from the FEM space and returns data in the DG space, while the dual nabla operator takes as input the data from the DG space and produces output in the FEM space. The two discrete nabla operators satisfy a discrete Schwarz theorem on the symmetry of discrete second derivatives. From there, both discrete vector calculus identities follow automatically. We apply our new family of schemes to three hyperbolic systems with involutions: the system of linear acoustics, in which the velocity field must remain curl-free and the vacuum Maxwell equations, in which the divergence of the magnetic field and of the electric field must remain zero. In our approach, only the magnetic field will remain exactly divergence free. As a third model we study the Maxwell-GLM system of Munz et al., which contains a unique mixture of curl-curl and div-grad operators and in which the magnetic field may be either curl-free or divergence-free, depending on the choice of the initial data. In all cases we prove that the proposed schemes are exactly total energy conservative and thus nonlinearly stable in the L2 norm.