One important characteristic of mixed finite element methods is their ability to provide accurate and locally conservative fluxes, an advantage over standard H1-finite element discretizations. However, the development of p or hp adaptive strategies for mixed formulations presents a challenge in terms of code complexity for the construction of Hdiv basis functions of high order on non-conforming meshes, and compatibility verification of the approximation spaces for primal and dual variables (inf − sup condition). In the present paper, a methodology is presented for the assembly of such approximation spaces based on quadrilateral and triangular meshes. In order to validate the computational implementations, and to show their consistent applications to mixed formulations, elliptic model problems are simulated to show optimal convergence rates for h and p refinements, using uniform and non-uniform (non-conformal) settings for a problem with smooth solution, and using adaptive hp-meshes for the approximation of a solution with strong gradients. Results for similar simulations using H1 are also presented, and both methods are compared in terms of accuracy and required number of degrees of freedom using static condensation.