  Numerical integration of the stiffness matrix in higher order finite element
(FE) methods is recognized as one of the heaviest computational tasks in a FE
solver. The problem becomes even more relevant when computing the Gram matrix
in the algorithm of the Discontinuous Petrov Galerkin (DPG) FE methodology.
Making use of 3D tensor-product shape functions, and the concept of sum
factorization, known from standard high order FE and spectral methods, here we
take advantage of this idea for the entire exact sequence of FE spaces defined
on the hexahedron. The key piece to the presented algorithms is the exact
sequence for the one-dimensional element, and use of hierarchical shape
functions. Consistent with existing results, the presented algorithms for the
integration of $H^1$, $H(\text{curl})$, $H(\text{div})$, and $L^2$ inner
products, have the $O(p^7)$ computational complexity. Additionally, a modified
version of the algorithms is proposed when the element map can be simplified,
resulting in the reduced $O(p^6)$ complexity. Use of Legendre polynomials for
shape functions is critical in this implementation. Computational experiments
performed with $H^1$, $H(\text{div})$ and $H(\text{curl})$ test shape functions
show good correspondence with the expected rates.
