There are different possibilities of choosing balanced pairs of approximation spaces for dual (flux) and primal (pressure) variables to be used in discrete versions of the mixed finite element method for elliptic problems arising in fluid simulations. Three cases shall be studied and compared for discretized three dimensional formulations based on tetrahedral, hexahedral and prismatic meshes. The principle guiding the constructions is the property that the divergence of the dual space and the primal approximation space should coincide, while keeping the same order of accuracy for the flux variable and varying the accuracy order of the primal variable. There is the classic case of BDMk spaces based on tetrahedral meshes and polynomials of total degree k for the dual variable, and k-1 for the primal variable, showing stable simulations with optimal convergence rates of orders k+1 and k, respectively. Another case is related to RTk and BDFMk+1 spaces for hexahedral and tetrahedral meshes, respectively, but holding for prismatic elements as well. It gives identical approximation order k+1 for both primal and dual variables, an improvement in accuracy obtained by increasing the degree of primal functions to k, and by enriching the dual space with some properly chosen internal shape functions of degree k+1, while keeping degree k for the border fluxes. A new type of approximation is proposed by further incrementing the order of some internal flux functions to k+2, and matching primal functions to k+1 (higher than the border fluxes of degree k). Thus, higher convergence rate of order k+2 is obtained for the primal variable. Using static condensation, the global condensed system to be solved in all the cases has same dimension (and structure), which is proportional to the space dimension of the border fluxes for each element geometry. Illustrating results comparing the three different space configurations are presented for simulations based on hierarchical high order shape functions for H(div)-conforming spaces, which are specially constructed for affine tetrahedral, hexahedral and prismatic meshes. Expected convergence rates are obtained for the flux, pressure and divergence variables. © 2016 Elsevier B.V.306 479502