I have not implemented these elements before, but I like their reduced cardinality compared to (e.g.) a tensor product of lagrange interpolants, which is very "overcomplete" (especially for orders>2)
I think they are not complete. My test program takes a quad and builds the 14x14 mass (projection?) matrix consisting the 8 serendipity functions + the 6 functions that comprise the second order polynomial basis (1 x y xx xy yy). Matlab reports that matrix has a rank of 9 (not 8), indicating there is a function inside the polynomial basis that the serendipity functions fail to capture. The thing is, the 9'th eigenvalue is like 10^5 down from the others, the serendipity basis "almost" captures all the polynomials, just barely doesn't. When I deliberately make the quad rectangular, the rank is reported as 8, the 9th eigenvalue vanishes to machine precision. That 9th eigenvalue grows as the quad is skewed more and more.
My suspicion has grown further upon reexamination of my reference material more closely (Huebner's "The Finite Element Method for Engineers", pp 167-171) - the serendipity functions are in the section entitled "Rectangular Elements" (e.g. orthogonal edges), and there's no mention of them in the "Arbitrary Quadrilateral Elements".
I have been googling and regoogling but thought somebody here might know off the top of their head, being more familiar with the space (seems like its more common in structures - I'm more an electromagnetics guy).
Thanks in advance for any help.