In my opinion the basics stay the same but there are some differences worth commenting upon.
The basic principles of defining a local coordinate system, integrating local dense contribution matrices using quadrature rules, and assembling them into global sparse matrices are identical. Often, hexahedral basis functions arise from taking tensor products of some underlying one-dimensional basis, and are integrated using tensor products of one-dimensional quadrature rules, which exposes some opportunities for optimizing code.
One notable difference is polynomial completeness. High order tetrahedral basis tend to arise from some sort of pascal's triangle / pascal's tetrahedron arrangement. For instance, a 2D gradient-conforming basis needs to span {$1$, $x$, $y$, $x^2$, $xy$, $y^2$} to be complete to second order, and these 6 functions can be arranged very tidily on a 6-node triangle (3 vertices, 3 edge midpoints). In contrast, the most natural way to write a second order basis for a 2D quadrilateral is via a tensor product of two 1D lines, leading to a span of {$1$, $x$, $y$, $x^2$, $xy$, $y^2$, $x^2y$, $xy^2$, $x^2y^2$}. This contains more functions / requires more work, but is not any more accurate in the asymptotic sense. This over-completeness becomes more exaggerated as you raise the order or migrate from 2D to 3D.
Despite this overcomplete-ness at the element level, hexahedral elements are generally more efficient in total cost, simply because they are "larger" and fill space more efficiently. Model accuracy is typically a function of mesh size $h$. A hexahedron with sides of length $h$ takes up $h^3$ volume, while a tetrahedron with sides of length $h$ only takes up a volume of $h^3/6$. So given some computational domain with volume $V$ and an desired accuracy/mesh size $h$, I expect to need 6x as many tetrahedra than hexahedra. 976
Even though this difference is just a constant factor, we're talking about the size of the input here, so it will be further magnified by any downstream algorithmic complexities (for instance, if you're using solver/preconditioning techniques that scale superlinearly).
I do want to comment briefly on accuracy. A surprising aspect of some basis functions is that element shape / distortion can effect their order of accuracy, because of how they depend on the pullback/jacobian. As a concrete example, consider 2D Nedelec (curl-conforming) functions. The curl of these functions (ie what's going into the "stiffness" matrix, to borrow a term) is proportional to the jacobian. For a triangle, the jacobian is constant, so you get curl complete to a constant (this is what you want, a lowest order "$\star$-conforming" element should have "$\star$" complete to a constant). In contrast, the jacobian for a quadrilateral is only constant if the shape mapping is affine (ie the quadrilateral is a parallelogram), which means this element is non-convergent on general meshes (ie you can refine $h$ but the error will not improve). Raviart-Thomas elements are similarly affected, as they have a similar dependence on jacobian. Note that a similar thing happens (non-constant jacobian / non-convergence) if you use curvilinear shape mappings (for either triangles or quadrilaterials), I only call out this case specifically because it's surprising that a linear/straight-line quadrilateral can exhibit this problem too (when it's not a parallelogram).
All this said, I find the biggest obstacle to using hexahedra for FEM are not formulation or convergence issues, but mesh generation and computational geometry issues. There are robust general-purpose algorithms for automatic tetrahedral mesh generation (delaunay refinement in particular), but hexahedral meshing is more brittle and often requires user interaction to decompose complicated geometry into simple parts that are eligible for structured strategies or pave/sweep strategies. Note that you can fuse tetrahedral grids with hexahedral grids using pyramidal elements, but the basis functions on pyramids can be pretty technical (pyramids are often implemented in terms of degenerate/distorted hexadra, and the distortion function ends up contributing non-polynomial terms to the jacobian/basis functions, which can make it difficult to obtain the correct order of accuracy). 441