  We generalize the two dimensional mixed finite elements of Arbogast and
Correa [T. Arbogast and M. R. Correa, SIAM J. Numer. Anal., 54 (2016), pp.
3332--3356] defined on quadrilaterals to three dimensional cuboidal hexahedra.
The construction is similar in that polynomials are used directly on the
element and supplemented with functions defined on a reference element and
mapped to the hexahedron using the Piola transform. The main contribution is
providing a systematic procedure for defining supplemental functions that are
divergence-free and have any prescribed polynomial normal flux. General
procedures are also presented for determining which supplemental normal fluxes
are required to define the finite element space. Both full and reduced
$H(\rm{div})$-approximation spaces may be defined, so the scalar variable,
vector variable, and vector divergence are approximated optimally. The spaces
can be constructed to be of minimal local dimension, if desired.
