The ability to effectively adapt a mesh is a very important feature of high fidelity finite element modeling. In a finite element analysis, a relatively high node density is desired in areas of the model where there are high error estimates from an initial analysis. Providing a higher node density in such areas improves the accuracy of the model and reduces the computational time compared to having a high node density over the entire model. Node densities can be determined for any model using the sizing functions based on the geometry of the model or the error estimates from the finite element analysis. Robust methods for mesh adaptation using sizing functions are available for refining triangular, tetrahedral, and quadrilateral elements. However, little work has been published for adaptively refining all hexahedral meshes using sizing functions. This thesis describes a new approach to drive hexahedral refinement based upon an error sizing function and a mechanism to compare the sizes of the node after refinement.