The conventional approach to construct quadratic elements for an n-sided polygon will yield n(n+1)/2 shape functions, which increases the computational effort. It is well known that the serendipity elements based on isoparametric formulation suffers from mesh distortion. Floater and Lai proposed a systematic way to construct higher-order approximations over arbitrary polygons using the generalized barycentric and triangular coordinates. This approach ensures 2n shape functions with nodes only on the boundary of the polygon. In this paper, we extend the polygonal splines approach to 3 dimensions and construct serendipity shape functions over hexahedra and convex polyhedra. This is done by expressing the shape functions using the barycentric coordinates and the local tetrahedral coordinates. The quadratic shape functions possess Kronecker delta property and satisfy constant, linear, and quadratic precision. The accuracy and the convergence properties of the quadratic serendipity shape elements are demonstrated with a series of standard patch tests. The numerical results show that the quadratic serendipity elements pass the patch test, yield optimal convergence rate, and can tolerate extreme mesh distortion. Copyright © 2017 John Wiley & Sons, Ltd.