This work is concerned with the numerical solution of Friedrichs systems by discontinuous Galerkin finite element methods (DGFEMs). Friedrichs systems are boundary value problems with symmetric, positive, linear first-order partial differential operators and allow the unified treatment of a wide range of elliptic, parabolic, hyperbolic and mixed-type equations. We do not assume that the exact solution of a Friedrichs system belongs to a Sobolev space, but only require that it is contained in the associated graph space, which amounts to differentiability in the characteristic direction. We show that the numerical approximations to the solution of a Friedrichs system by the DGFEM converge in the energy norm under hierarchical h- and p- refinement. We introduce a new compatibility condition for the boundary data, from which we can deduce, for instance, the validity of the integration-by-parts formula. Consequently, we can admit domains with corners and allow changes of the inertial type of the boundary, which corresponds in special cases to the componentwise transition from in- to outflow boundaries. To establish the convergence result we consider in equal parts the theory of graph spaces, Friedrichs systems and DGFEMs. Based on the density of smooth functions in graph spaces over Lipschitz domains, we study trace and extension operators and also investigate the eigensystem associated with the differential operator. We pay particular attention to regularity properties of the traces, that limit the applicability of energy integral methods, which are the theoretical underpinning of Friedrichs systems. We provide a general framework for Friedrichs systems which incorporates a wide range of singular boundary conditions. Assuming the aforementioned compatibility condition we deduce well-posedness of admissible Friedrichs systems and the stability of the DGFEM. In a separate study we prove hp-optimality of least-squares stabilised DGFEMs.